基于曲率屬性的構(gòu)造應(yīng)力預(yù)測方法

馬 妮 印興耀 宗兆云 孫成禹 王世星

(①中國石化石油物探技術(shù)研究院，江蘇南京 211103； ②中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ③海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266071)

0 引言

構(gòu)造應(yīng)力導(dǎo)致巖層發(fā)生形變，形成不同形態(tài)的地質(zhì)構(gòu)造，控制著裂縫的形成和分布，影響油氣的運移。因此，構(gòu)造應(yīng)力的預(yù)測不僅有助于人們認(rèn)識和理解地質(zhì)構(gòu)造的形成、分布和演化規(guī)律，還可以為油氣田注水開發(fā)過程中井網(wǎng)的布置和調(diào)整提供依據(jù)[1]。

早期學(xué)者采用物理模擬方法對構(gòu)造應(yīng)力場分布特征進(jìn)行研究，以了解構(gòu)造變形的過程，但該方法存在時空尺度的局限性和實驗條件的限制。其后數(shù)值模擬技術(shù)被引入到構(gòu)造應(yīng)力場的研究中，建立的地質(zhì)模型實現(xiàn)了二維和三維空間的構(gòu)造應(yīng)力場數(shù)值模擬[2]，克服了時空尺度的局限性。譚成軒等[3]利用Super SAP有限元程序探索了含油氣盆地三維構(gòu)造應(yīng)力場數(shù)值模擬方法。陳書平等[4]通過設(shè)計的4個地質(zhì)模型，采用有限元法模擬分析了盆地演化過程中的四期應(yīng)力場。張帆等[5]針對特定的地震反射層位，利用連續(xù)介質(zhì)彈性有限元模型計算主應(yīng)力和剪切應(yīng)力，并預(yù)測了研究區(qū)內(nèi)裂縫發(fā)育情況。王紅才等[6]基于有限元分析方法利用計算機模擬三維構(gòu)造應(yīng)力場，并將該方法應(yīng)用于遼河油田灘海地區(qū)。雷景生[7]研究了基于結(jié)構(gòu)優(yōu)化的有限元網(wǎng)格自動劃分方法，并對某油藏的應(yīng)力場進(jìn)行數(shù)值模擬。王連捷等[8]利用有限元法模擬的三維構(gòu)造應(yīng)力場預(yù)測油氣聚集有利區(qū)。付玉華等[9]提出了構(gòu)造應(yīng)力場邊界載荷的有限元反演方法，用于反演構(gòu)造應(yīng)力場的邊界力。董培育等[10]利用蒙特卡洛方法，根據(jù)庫倫─摩爾破裂準(zhǔn)則反演巴顏喀拉塊體的初始構(gòu)造應(yīng)力場。構(gòu)造應(yīng)力的數(shù)值模擬方法雖克服了時空尺度的局限性，但在模型構(gòu)建、邊界條件和所受載荷類型等方面存在主觀性，且模擬的結(jié)果存在多解性和可靠性等問題[2]。

曲率屬性和構(gòu)造應(yīng)力都與裂縫相關(guān)，國內(nèi)外學(xué)者研究了曲率屬性與構(gòu)造應(yīng)力的關(guān)系。Price等[11]利用曲率屬性估算構(gòu)造應(yīng)力，首次建立了應(yīng)力與曲率、彈性力學(xué)參數(shù)之間的定量關(guān)系。Sheorey[12]在考慮地層曲率和巖層彈性參數(shù)、密度與熱膨脹系數(shù)變化的前提下，利用靜彈性熱應(yīng)力模型估算地層應(yīng)力。Roberts[13]推導(dǎo)了層面曲率的計算公式，并利用彎曲的巖層模型描述曲率與應(yīng)力的關(guān)系。Sigismondi等[14]研究了曲率與構(gòu)造之間的相關(guān)性。劉建偉等[15]基于測井巖石物理參數(shù)分析，明確了地應(yīng)力、破裂壓力與巖石彈性參數(shù)的關(guān)系，實現(xiàn)了非常規(guī)油藏應(yīng)力場及應(yīng)力甜點的地震預(yù)測。孟陽等[16]應(yīng)用GeoEast系統(tǒng)提取識別斷層的敏感分方位地震屬性，精細(xì)識別了ZJB區(qū)塊五峰—龍馬溪組頁巖儲層的小斷層。王玲玲等[17]基于多裂縫參數(shù)地震物理模型，分析疊后屬性，優(yōu)選敏感地震屬性進(jìn)行融合以識別裂縫。何英[18]利用曲率屬性估算地層的構(gòu)造應(yīng)力，克服了復(fù)雜構(gòu)造應(yīng)力場模型構(gòu)建的難點，但是沒有給出最大正曲率和最小負(fù)曲率與構(gòu)造應(yīng)力的直接關(guān)系。

利用曲率屬性估算構(gòu)造應(yīng)力的傳統(tǒng)方法是通過提取x、y方向的曲率和xOy面的扭率計算相應(yīng)的應(yīng)力分量而求取構(gòu)造應(yīng)力，但沒有構(gòu)建曲率屬性與構(gòu)造應(yīng)力之間的定量關(guān)系，且目前的商業(yè)軟件通常不能直接提取x、y方向的曲率和xOy面的扭率，只能提取最大正曲率、最小負(fù)曲率、高斯曲率、平均曲率、極大曲率和極小曲率等屬性。為此，本文在薄板彎曲理論的假設(shè)條件下，通過研究曲率與構(gòu)造應(yīng)力之間的定量關(guān)系，推導(dǎo)由最大正曲率、最小負(fù)曲率直接表征構(gòu)造應(yīng)力的計算公式，即通過提取體曲率屬性得到最大正曲率和最小負(fù)曲率，進(jìn)而計算地層構(gòu)造應(yīng)力。

1 方法原理

1.1 彈性薄板彎曲理論

薄板是指板的厚度遠(yuǎn)小于其他方向尺度的兩個平行平面所圍成的板狀構(gòu)件。設(shè)薄板的中平面為坐標(biāo)平面，即z=0。坐標(biāo)原點O為中平面內(nèi)一點，x和y軸在中平面內(nèi)。在x軸正方向上產(chǎn)生的位移為u，y軸正方向上產(chǎn)生的位移為v，z軸正方向上產(chǎn)生的位移為擾度w(x,y)。z軸向上，薄板的時間厚度為z(單位為s)。薄板模型如圖1所示。

一般來說，σx、σy和σz為正應(yīng)力分量(下標(biāo)表示受力面的法線方向)，對應(yīng)的應(yīng)變分量分別為εx、εy和εz；τxy、τxz、τyx、τyz、τzy、τzx為切應(yīng)力分量(下標(biāo)中前者表示受力面的法線方向，后者表示應(yīng)力分量的方向)，對應(yīng)的應(yīng)變分量分別為γxy、γxz、γyx、γyz、γzy、γzx[19]。

薄板小擾度彎曲的近似理論[20]是以下述三個基本假設(shè)條件為基礎(chǔ)。

(1)在平行于中平面的截面上，應(yīng)力分量τzx和τzy、σz相比于其余應(yīng)力分量要小得多，因此它們所引起的應(yīng)變γzx、γzy可忽略不計，即

(1)

由式(1)可得

(2)

(2)中平面垂直方向上的正應(yīng)變εz可以不計，即

(3)

式(3)表明薄板中平面上法線的各點都具有相同的位移w，即為擾度。

(3)薄板中平面內(nèi)的各點都沒有平行于中平面的位移，表示為

(4)

由式(4)可得

(5)

圖1 薄板模型

1.2 薄板彎曲理論假設(shè)下的廣義胡克定律

薄板的小擾度彎曲問題是對位移進(jìn)行求解，因此選取擾度w(x,y)作為未知函數(shù)，對式(2)中的z變量積分可得擾度w表示的縱向位移u和v[19]，即

(6)

式中f1(x,y)和f2(x,y)為任意函數(shù)，進(jìn)一步由式(4)可得

(7)

彈性力學(xué)的基本方程包括平衡方程、物理方程、幾何方程和協(xié)調(diào)方程等。本文選用其中的物理方程和幾何方程進(jìn)行相關(guān)推導(dǎo)。根據(jù)薄板小擾度彎曲近似理論的假設(shè)條件(1)和(2)可知，應(yīng)變εz、γzx和γzy忽略不計，則由彈性力學(xué)的幾何方程可得到擾度w表示的主要應(yīng)變分量，即

(8)

根據(jù)薄板彎曲理論的假設(shè)條件，彈性力學(xué)的廣義胡克定律可進(jìn)一步表示為

(9)

式中：E表示楊氏模量；μ表示泊松比。將式(9)變形得到

(10)

1.3 地層曲率屬性的提取

Roberts[13]提出層面曲率的計算方法，即采用最小二乘法或應(yīng)用其他的近似法，擬合一個二次曲面，得到地層面的趨勢函數(shù)，進(jìn)而估算地層面上某點的曲率分量。其二元二次方程為

w(x,y)=ax2+by2+cxy+dx+ey+f

(11)

計算曲面某一點曲率的方法有多種，其中分波數(shù)的傅氏分析法，利用傅里葉變換的微分性質(zhì)，將空間域的某變量變換到波數(shù)域，乘以-ik因子后再經(jīng)過傅氏反變換得到該變量在空間域的一階導(dǎo)數(shù)[21-26]。根據(jù)地層傾角的定義可知，曲面的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)分別表示地層的視傾角及其一階導(dǎo)數(shù)，利用地層的視傾角及其一階導(dǎo)數(shù)可以計算擬合曲面方程中系數(shù)a～e。計算公式[13]為

(12)

式中p和q為x和y方向地層的視傾角屬性。由于各種曲率(如高斯曲率、平均曲率、極大曲率、極小曲率、最大正曲率和最小負(fù)曲率等)的計算公式中沒有包含系數(shù)f，因此無需對f進(jìn)行計算。在求取擬合曲面方程的系數(shù)后，可以計算各種曲率。本文只需要計算最大正曲率和最小負(fù)曲率

(13)

式中Kpos、Kneg分別表示最大正曲率、最小負(fù)曲率。在薄板小擾度彎曲理論的假設(shè)條件下，薄板中平面內(nèi)各點在不同方向上的曲率與擾度w的關(guān)系可以近似表示為[27]

(14)

式中Kx、Ky和Kxy分別為x、y方向的曲率和xOy面的扭率，由式(8)、式(14)可得

(15)

由式(10)可得

(16)

由曲率的性質(zhì)可知，當(dāng)?shù)貙酉蛏贤蛊饡r，曲率為正，地層所受的應(yīng)力為張應(yīng)力；當(dāng)?shù)貙酉蛳掳枷輹r，曲率為負(fù)，地層所受的應(yīng)力為壓應(yīng)力。為了與地質(zhì)力學(xué)體系一致，本文規(guī)定張應(yīng)力為負(fù)，壓應(yīng)力為正[28]。

1.4 基于曲率屬性的構(gòu)造應(yīng)力計算公式

地層中的主應(yīng)力及其方向計算公式為[18,28]

(17)

式中：σmax、σmin分別為最大、最小主應(yīng)力；α為σmax與x軸的夾角；β為σmin與y軸的夾角。

最大正曲率和最小負(fù)曲率是描述裂縫、撓曲、褶皺和斷層較有效的曲率屬性[29]。為了更好地利用曲率屬性表征地層中的主應(yīng)力，將式(14)代入式(16)，得

(18)

由式(13)、式(17)和式(18)經(jīng)過一系列推導(dǎo)可以得到地層主應(yīng)力與最大正曲率、最小負(fù)曲率之間關(guān)系的表達(dá)式

(19)

由于式(17)、式(19)是在線彈性模型(薄板彎曲理論的假設(shè)條件)下推導(dǎo)得到的，而實際地層在構(gòu)造形變過程中是彈塑性形變，因此在薄板彎曲理論假設(shè)條件下估算的構(gòu)造應(yīng)力是一種簡單的近似簡化過程，忽略了地層塑性形變的影響。由式(19)可以得到最大和最小水平主應(yīng)力變化差與最大水平主應(yīng)力的比值

(20)

最大、最小主應(yīng)力之間的差異越大，巖石越易于發(fā)生破裂而形成裂縫。

1.5 巖石力學(xué)參數(shù)反演

根據(jù)式(19)估算構(gòu)造應(yīng)力，除了要提取曲率屬性，還要利用疊前地震資料反演巖石力學(xué)參數(shù)。本文采用基于彈性阻抗反演的巖石力學(xué)參數(shù)疊前地震反演方法獲得楊氏模量和泊松比。

宗兆云等[30-31]、印興耀等[32-33]在Aki-Richards反射系數(shù)近似方程的基礎(chǔ)上，推導(dǎo)了YPD近似方程，該方程是縱波反射系數(shù)與楊氏模量、泊松比及密度等反射系數(shù)的線性關(guān)系

(21)

(22)

式中BEI表示脆性彈性阻抗。在界面兩側(cè)介質(zhì)的參數(shù)差異較小時，可用微分轉(zhuǎn)換關(guān)系表示，即

(23)

(24)

其中

式中BEI0、E0、μ0和ρ0分別是入射角為0°時BEI、E、μ和ρ的平均值。

對式(24)兩邊取對數(shù)，得到易于求解的線性化形式為[30-31]

(25)

基于式(25)構(gòu)建的方程組，利用疊前彈性阻抗反演得到的彈性阻抗數(shù)據(jù)體及測井?dāng)?shù)據(jù)等資料，可以得到任意采樣點的楊氏模量和泊松比。將得到的楊氏模量和泊松比代入到構(gòu)造應(yīng)力計算公式(式(19))中，結(jié)合基于局部τ-p變換的曲率屬性提取方法[22]得到的曲率及地層厚度z，就可以估算地層的構(gòu)造應(yīng)力。

與傳統(tǒng)的構(gòu)造應(yīng)力求取公式(式(16)和式(17))相比較，本文方法不需要計算Kx、Ky和Kxy參數(shù)。目前Kx、Ky和Kxy參數(shù)通常不能由商業(yè)軟件直接提取，需編寫相關(guān)程序計算擬合曲面方程(式(11))的系數(shù)進(jìn)而得到。而最大正曲率和最小負(fù)曲率屬性可通過主流商業(yè)軟件直接獲得。因此，本文方法實用、簡單易行。

2 應(yīng)用實例

將本文方法應(yīng)用于中國東部某頁巖段(S層)裂縫型油氣藏區(qū)。研究區(qū)內(nèi)發(fā)育一系列北東向、北西向及近東西向的斷裂，斷裂附近的裂縫較為發(fā)育。裂縫呈“X”型或樹狀，裂縫密度西部較東部大。

為了驗證本文方法，對比兩種方法求取構(gòu)造應(yīng)力的效果。

按照前人傳統(tǒng)方法，提取Kx、Ky和Kxy(圖2)；結(jié)合疊前地震反演得到的楊氏模量和泊松比，通過式(16)求得σx、σy和σxy(圖3)；最后基于式(17)求取S層最大主應(yīng)力(圖4a)。

圖2 S層Kx(a)、Ky(b)和Kxy(c)沿層切片

利用本文地層曲率屬性的提取方法求得最大正曲率、最小負(fù)曲率屬性(圖5)。應(yīng)用疊前地震反演的巖石力學(xué)參數(shù)獲取儲層的楊氏模量和泊松比，根據(jù)式(19)求取S層最大主應(yīng)力(圖4b)。

從圖4可以看出，兩種方法計算得到的構(gòu)造應(yīng)力沿層切片結(jié)果基本一致。

由于本文方法所用的式(19)是基于前人方法的式(17)推導(dǎo)得到的，在推導(dǎo)過程中沒有采用近似簡化的過程，直接表達(dá)了曲率屬性與構(gòu)造應(yīng)力之間的關(guān)系，因此本文方法所得到的構(gòu)造應(yīng)力沿層切片與前人方法得到的構(gòu)造應(yīng)力切片基本一致，從而證明了本文方法的正確性。

圖3 S層σx(a)、σy(b)和σxy(c)沿層切片

圖4 不同方法得到的S層最大主應(yīng)力沿層切片對比(a)傳統(tǒng)方法； (b)本文方法

3 結(jié)束語

在薄板彎曲理論的假設(shè)條件下，研究了曲率屬性與構(gòu)造應(yīng)力之間的定量關(guān)系，推導(dǎo)了利用最大正曲率和最小負(fù)曲率屬性表征構(gòu)造應(yīng)力的計算公式，提出了直接利用提取的最大正曲率和最小負(fù)曲率屬性估算地層構(gòu)造主應(yīng)力的方法。

該方法利用相關(guān)商業(yè)軟件就可以得到所需的最大正曲率和最小負(fù)曲率屬性，不需要通過提取x、y方向的曲率和xOy面的扭率計算相應(yīng)的應(yīng)力分量進(jìn)而求取構(gòu)造應(yīng)力，更加實用、簡單易行。