王文文 李 諾 韓建強(qiáng) 羅 濤 肖 濤

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

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

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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

圖1 利用MGGP 進(jìn)化的一種典型模型結(jié)構(gòu)Fig.1 A typical structure of model evolved by MGGP

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

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

(2)計(jì)算種群中每個(gè)個(gè)體的適應(yīng)度；

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

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

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

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

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

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

3 實(shí)驗(yàn)測量與結(jié)果

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

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

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

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

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

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

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

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

4 MGGP預(yù)測

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

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

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

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

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

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

5 結(jié)論

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