王文文 李 諾 韓建強 羅 濤 肖 濤

(1 中海油田服務股份有限公司油田技術事業(yè)部 北京 101149)

(2 中國科學院大學 北京 100049)

(3 北京市海洋深部鉆探測量工程技術中心 北京 100190)

(4 中國電子科技集團公司第三研究所 北京 100015)

(5 中國石油渤海鉆探工程有限公司工程技術研究院 任丘 062550)

0 引言

掌握巖石的力學性質有助于解決眾多地球科學、巖土工程等領域的問題，特別是在石油和天然氣的勘探和開發(fā)領域。楊氏模量是最常見的力學參數(shù)，這些參數(shù)可以從巖石壓縮測試(也稱靜態(tài)測量)中測量的應力-應變關系中獲得。但是靜態(tài)測量實驗裝置復雜，對人員的操作要求高，其過程通常十分耗時而且測試費用昂貴，因此在井中的整個深度區(qū)間實施連續(xù)取心以獲得巖石力學參數(shù)的連續(xù)剖面是不現(xiàn)實的，但該剖面在分析井壁穩(wěn)定性、出砂、壓裂等方面至關重要[1]。實驗室彈性波速度的動態(tài)測試相對來說十分高效，而且成本很低，一般利用超聲波脈沖傳輸法測量縱波和橫波速度來計算動態(tài)彈性性質[2]。在現(xiàn)場，也可以利用地震勘探或聲波測井得到的地震波或聲波速度來獲得深度連續(xù)的動態(tài)彈性剖面[3]。

由于儲層巖石復雜的孔隙結構，并且孔隙空間通常含有流體，其動態(tài)和靜態(tài)性質之間往往存在很大差異[4]。通常對實驗室測量的少量巖心的動靜態(tài)數(shù)據(jù)進行統(tǒng)計分析，利用擬合方法(一般是線性擬合)得到二者的經(jīng)驗關系，以此對現(xiàn)場的動態(tài)數(shù)據(jù)進行校正來得到靜態(tài)性質[5]。然而，由于信號衰減較大，計算巖石動態(tài)性質的必要參數(shù)——橫波速度在實驗室中往往難以測量；而在測井中，由于聲波測井方法和地層性質的限制，橫波速度也難以獲得[6-7]。這些限制使得動態(tài)、靜態(tài)性質的轉換問題更加復雜。此外，在回歸建模中，需要用有限的參數(shù)和方程預定義經(jīng)驗模型的結構，現(xiàn)有的經(jīng)驗關系存在精度較低、適用性不高的問題。為了克服這些局限性，提出一種高精度的方法，從有限的測量數(shù)據(jù)中預測巖石靜態(tài)性質。

近年來，人工神經(jīng)網(wǎng)絡、模糊推理系統(tǒng)等人工智能方法已經(jīng)用于解決此類問題[8]，但是它們僅適用大數(shù)據(jù)集，同時不能提供靜態(tài)性質預測的解析表達式[9-10]。本文提出一種基于多基因遺傳規(guī)劃(Multi-gene genetic programming，MGGP)的儲層巖石靜態(tài)性質預測方法，可以不需要先驗信息建立預測方程。首先，介紹了動態(tài)和靜態(tài)模量的定義及其區(qū)別的內在原因。其次，描述了MGGP算法的特點和流程。最后，通過一組實驗數(shù)據(jù)分別研究了線性擬合和MGGP 在靜態(tài)模量計算中的應用。

1 靜態(tài)和動態(tài)楊氏模量

靜態(tài)楊氏模量是從單軸壓縮測試中獲得的應力-應變曲線直線段的斜率，所謂單軸壓測試，指的是圓柱形樣品在無側向圍壓下承受端部平面上的力。值得注意的是，由于巖石具有非均質、非線性的特點，其應力-應變曲線并不是一條直線，由此可按照斜率的不同定義3 種模量：切線模量、割線模量和平均模量。本文討論的是平均模量及應力-應變曲線近似直線段的斜率，在此段微裂隙完全閉合，可以認為巖石發(fā)生的是線彈性變形。此時靜態(tài)楊氏模量Esta定義為軸向應力增量Δσz與軸向應變Δεz的比值：

實驗室通常采用超聲脈沖透射法測量樣品的縱橫波速度，根據(jù)彈性波傳播理論，動態(tài)楊氏模量Edyn定義如下：

其中，ρb是體密度，VP和VS分別為縱橫波速度。

儲層巖石由于具有復雜的孔隙結構，靜態(tài)加載過程中應力加載速率較慢，可以近似為0 頻，變形幅度較大(＞10-4)，變形包括彈性和非彈性成分[11]。而在動態(tài)測量中，由于頻率較高，巖石只承受瞬時應力，而且變形幅度很小(＜10-6)，測量過程中的變形可以看作是彈性的。如果巖石飽和流體，由于作用頻率的差異導致靜態(tài)測量過程中孔隙壓力處處平衡，動態(tài)測量過程中巖石處于不排水狀態(tài)，孔壓不平衡，情況更加復雜[12-13]。因此，無論是實驗室超聲透射法還是現(xiàn)場聲波測井方法，得到的動態(tài)模量都與巖石的靜態(tài)模量都會有差異，甚至差異較大。

2 多基因遺傳規(guī)劃

遺傳規(guī)劃(Genetic programming，GP)是由Koza[14]提出的一種進化算法。與遺傳算法(Genetic algorithm，GA)類似，GP 以達爾文進化論為基礎，借鑒了生物界的自然選擇和遺傳機制。在一個由多個個體組成的群體中，每個個體都代表一個可行的解，通過選擇、復制、交換、變異等遺傳操作，根據(jù)最優(yōu)適應度逐步迭代，得到最優(yōu)解。然而不同于只能對字符串進行操作的GA，GP是一種符號回歸方法，可以對樹狀結構的計算機程序進行操作，將輸入與輸出以解析方程的形式聯(lián)系起來。多基因遺傳規(guī)劃(MGGP)是GP 的變體，結合了GP 的模型構建能力和統(tǒng)計回歸方法的參數(shù)估計能力[15]。在傳統(tǒng)的遺傳規(guī)劃中，進化模型由一棵樹組成，而在多基因遺傳規(guī)劃中，每一個回歸模型都是多棵樹的加權線性組合。利用MGGP 進化得到的模型結構示例如圖1所示，函數(shù)集和符號集分別構成樹的節(jié)點和葉。在圖1中，函數(shù)集包括“l(fā)og，+，×，sin”，終止符集包括“A，x1，x2”，d0是偏置項，d1和d2是權重系數(shù)。與GP 相比，MGGP 能夠使預測模型具備非線性特征，從而具有更高的精度和效率[15-16]。

圖1 利用MGGP 進化的一種典型模型結構Fig.1 A typical structure of model evolved by MGGP

圖2給出了MGGP的工作流程，可以概括如下：

(1)隨機產(chǎn)生初代種群，包含一組由函數(shù)和變量組成的個體；

(2)計算種群中每個個體的適應度；

(3)選擇適應度好的優(yōu)良個體為母體；

(4)通過突變、交叉、復制等方式產(chǎn)生新的個體，創(chuàng)造新的一代(子代)；

(5)重復步驟(2)和步驟(3)，直到滿足終止準則。

圖2 MGGP 算法的一般流程Fig.2 The general workflow of MGGP algorithm

用于評估模型的適應度函數(shù)是實際值ydata和預測值ypre之間的均方根誤差(Root mean square error，RMSE)：

其中，N為參與評估的樣本數(shù)目。用相關系數(shù)R來表示測量值與預測值之間的相關關系密切程度，即兩個變量之間的協(xié)方差和兩者標準差乘積的比值。

3 實驗測量與結果

實驗中，共對32 塊井下取心的砂巖樣品進行實驗測量與分析，相關實驗在中國科學院聲學研究所的高溫高壓巖石動靜態(tài)三軸實驗系統(tǒng)GCTS RTR-2000 進行，如圖3所示。該裝置包括主鋼架、壓力控制器、超聲波控制器等，能滿足樣品干燥和飽和態(tài)下的聲學和力學參數(shù)測量。該系統(tǒng)最大作用力2000 kN，內置精度為0.15%的載荷傳感器，能夠進行軸向應力加載；儀器配置有圍壓孔壓加載裝置，最大加載量140 MPa，精度為0.25%，可以進行圍壓加載；實驗中軸向應變加載速率為每秒0.001%，一分鐘測量一次；采用量程±3 mm的線性可變差動變壓器(Linear variable differential transformer，LVDT)變形傳感器測量應變，精度為0.25%。靜態(tài)楊氏模量Esta是應力-應變曲線在30 MPa～35 MPa的應力區(qū)間段內的擬合直線斜率，因為在該應力區(qū)間內曲線近似為一條直線，可認為巖石發(fā)生線彈性變形，見圖4。采用超聲脈沖透射法測量了相同應力范圍內的縱波速度VP和橫波速度VS，結果取為整數(shù)，單位為m/s。

圖3 GCTS RTR-2000 巖石力學三軸測試系統(tǒng)Fig.3 The GCTS RTR-2000 rock mechanic triaxial test system

圖4 一塊巖心樣品的單軸應力-應變曲線Fig.4 Stressstrain curve of a core sample

用作為完全彈性體的標準鋁樣進行測試，得到系統(tǒng)測量誤差如下：(1)應力應變測量誤差為0.25%；(2)縱波速度測量誤差為0.52%，橫波速度測量誤差為0.77%；(3)動態(tài)模量測試誤差為1%，靜態(tài)模量測試誤差為0.5%。

結合測量的體密度ρb，用式(2)計算動態(tài)楊氏模量Edyn，結果保留小數(shù)點后兩位，單位為GPa。通過SCMS-C3 型覆壓孔滲測量儀測得樣品孔隙度φ，結果保留小數(shù)點后四位。所有測量結果在表1中給出，前30塊樣品為飽和測量結果(-S表示)，后30 塊樣品為干燥測量結果(-D表示)。由于測量過程中信號的強衰減，其中4 塊飽和樣品和1 塊干燥樣品的橫波信號由于信噪比較低，沒有獲得可靠的測量值。而且樣品31 和32 在干燥測量后，樣品出現(xiàn)明顯裂隙，沒有進行飽和測量。

表1 樣品參數(shù)測量結果Table1 Parameters measurement results of samples

繪制干燥和飽和條件下的動靜態(tài)模量交會圖如圖5所示。動態(tài)楊氏模量Edyn和靜態(tài)楊氏模量Esta之間存在較好的線性相關性，通過線性擬合得到二者之間的線性關系：

圖5 干燥和飽和樣品的動靜態(tài)模量Fig.5 Dynamic Young’s modulus compared to static Young’s modulus of dry and saturated rocks

4 MGGP預測

在MGGP 方法中，靜態(tài)楊氏模量為待預測參數(shù)，而終止符集中的體密度、孔隙度和縱波速度是預測參數(shù)，函數(shù)集包含運算符“+、-、*、/、sqrt、ln(x)、log(x)、exp(x)、abs(x)、square、cubic、power、sinh(x)、cosh(x)、sin(x)、cos(x)”。終止符集和函數(shù)集用來構建MGGP 模型。隨機選取70%的數(shù)據(jù)進行模型的訓練，其余30%的數(shù)據(jù)進行測試。表2列出了構建楊氏模預測模型所用的MGGP 模型參數(shù)。

表2 預測靜態(tài)楊氏模量的MGGP 模型參數(shù)Table2 MGGP model parameters used to predict static Young’s modulus

在進化300 代之后適應度函數(shù)均方根誤差RMSE 值為3.44，得到了一個關于孔隙度φ、體密度ρb和縱波速度VP的靜態(tài)楊氏模量方程：

為了驗證式(5)的效果和適用性，將MGGP 預測的結果和動態(tài)數(shù)據(jù)標定值(式(4))分別與實際測量值進行了比較，結果如圖6和圖7所示。在訓練集和測試集中，MGGP 模型都能有效地預測楊氏模量，訓練集的相關系數(shù)為0.865，測試集的相關系數(shù)為0.726。對于線性擬合方程預測結果，5 塊樣品由于橫波速度缺失而缺少經(jīng)驗預測值(在圖中顯示為0 值)。除去缺少數(shù)據(jù)的5 塊樣品，線性擬合結果與測量值的相關系數(shù)為0.829。由此可見，與利用動態(tài)數(shù)據(jù)得到的線性模型相比，MGGP 訓練集中預測的楊氏模量更接近實測靜態(tài)值，并且對于測試集也有很高的精度，同時避免了橫波數(shù)據(jù)缺失無法直接計算的影響。

圖6 訓練集線性擬合、MGGP 模型靜態(tài)楊氏模量預測值與實驗測量值對比Fig.6 Static Young’s moduli from experiment，linear fit and MGGP in training set

圖7 測試集線性擬合、MGGP 模型靜態(tài)楊氏模量預測值與實驗測量值對比Fig.7 Static Young’s moduli from experiment，linear fit and MGGP in testing set

5 結論

本文采用MGGP 方法對儲層砂巖樣品的靜態(tài)楊氏模量進行了預測。該模型比用動態(tài)數(shù)據(jù)線性擬合標定的經(jīng)驗關系有更好的適應性。此外，在缺失橫波速度的情況下，它也能夠利用較少的輸入?yún)?shù)，即孔隙度、體積密度和縱波速度來進行預測，從而減小了橫波波速度不準確或丟失所引起的計算誤差。結果表明，MGGP 在實驗中用來預測靜態(tài)楊氏模量是可行的，這為現(xiàn)場數(shù)據(jù)(測井或地震)在解決巖石力學問題提供了更好的數(shù)據(jù)處理手段。