熊曉軍 崔澤飛 巫芙蓉 龔思宇 李 翔 劉 陽

(1.成都理工大學(xué)油氣藏地質(zhì)與開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室 四川成都 610059；2.中國石油川慶鉆探工程公司地球物理勘探公司 四川成都 610213)

地層斷裂帶的發(fā)育與構(gòu)造應(yīng)力場(chǎng)關(guān)系密切，建立地層構(gòu)造應(yīng)力場(chǎng)可用于裂縫預(yù)測(cè)[1]。常規(guī)地應(yīng)力場(chǎng)構(gòu)建方法可以分為4類：第1類是基于經(jīng)驗(yàn)關(guān)系式的方法，如賈立宏 等[2]提出的地應(yīng)力場(chǎng)灰色模擬方法；第2類是進(jìn)行應(yīng)力場(chǎng)的數(shù)值模擬，如張帆 等[3]提出的構(gòu)造應(yīng)力場(chǎng)數(shù)值模擬方法；第3類是進(jìn)行應(yīng)力場(chǎng)的反演，如張國強(qiáng) 等[4]提出的基于神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)分界的應(yīng)力場(chǎng)反演方法；第4類是基于三維地震資料進(jìn)行應(yīng)力場(chǎng)的計(jì)算，如何英[5]提出的基于二維面曲率的應(yīng)力場(chǎng)計(jì)算方法等。其中，第1類和第3類基于測(cè)井?dāng)?shù)據(jù)通過數(shù)學(xué)方法推算整個(gè)區(qū)域的地應(yīng)力場(chǎng)，其結(jié)果缺乏準(zhǔn)確性；第2類和第4類通過地質(zhì)解釋的構(gòu)造信息計(jì)算地應(yīng)力場(chǎng)，其結(jié)果分辨率較低。

為了克服上述4類方法的計(jì)算精度受目的層層位解釋精度影響的缺陷，本文引入三維傾角屬性對(duì)其進(jìn)行改進(jìn)。該方法首先采用基于離心窗掃描的傾角計(jì)算方法[6]，通過自動(dòng)掃描地震數(shù)據(jù)獲得地層傾角，進(jìn)而計(jì)算地層曲率(相對(duì)于一般的中心窗和多窗掃描方法，離心窗掃描方法對(duì)地層細(xì)微構(gòu)造的刻畫更加精確[7])；然后再采用有限差分方法和Price(1990年)提出的“曲率-應(yīng)變-應(yīng)力”方程[8]計(jì)算地層構(gòu)造應(yīng)力場(chǎng)。渤海某工區(qū)花崗巖儲(chǔ)層實(shí)際資料應(yīng)用表明，本文方法的預(yù)測(cè)結(jié)果與測(cè)井解釋的一致性較好，可以有效應(yīng)用于研究區(qū)裂縫預(yù)測(cè)，并且具有一定的推廣意義。

1 方法原理

1.1 基于地層傾角的地應(yīng)力計(jì)算方法

地層曲率K是地層傾角Ψ關(guān)于弧長S的偏導(dǎo)數(shù)，即K=?Ψ/?s。當(dāng)?shù)貙觾A角Ψ較小時(shí)，Ψ≈tanΨ=?z/?x,z(x,y)表示地層曲面,弧長S≈x，則曲率可近似表達(dá)為K=?2z/?x2。為了方便后續(xù)的公式推導(dǎo)，在此定義變量Kxy，它的物理意義是地層的扭曲度。則地層在x方向(inline方向)和y方向(crossline方向)的曲率以及地層的扭曲度為

(1)

定義地層沿x和y軸方向的視傾角分別為p和q[9]。實(shí)際上p和q是地層在地震剖面上x和y方向的斜率[10],即曲面z(x,y)沿x和y方向的一階偏導(dǎo)數(shù)

(2)

根據(jù)上述關(guān)系式可以推導(dǎo)出傾角表示的地層曲率方程為

(3)

因?yàn)閜和q在剖面上是離散的，所以p和q對(duì)x和y的偏導(dǎo)數(shù)可以通過簡(jiǎn)單的顯式有限差分方法[11]求出，即

(4)

(5)

其中，p1～p9在網(wǎng)格點(diǎn)[12]的位置如圖1所示,q1～q9位置可參考圖1。

圖1 離散數(shù)據(jù)點(diǎn)位置關(guān)系Fig .1 Location of discrete data points

在彈性力學(xué)中，根據(jù)變形幾何方程和薄板假設(shè)理論，令曲面代表的地層時(shí)間厚度為0，則該地層在構(gòu)造應(yīng)力作用下的應(yīng)變可以表示為

(6)

式(6)中：εx和εy分別表示x和y方向的應(yīng)變;γxy表示剪應(yīng)變;h表示地層厚度，在這里可以設(shè)為1。將式(3)代入式(6)中得應(yīng)變與曲率的關(guān)系式為

(7)

由廣義胡克定律，進(jìn)一步得到應(yīng)力與應(yīng)變關(guān)系式

(8)

式(8)中:σx和σy分別表示x和y方向的應(yīng)力;E代表?xiàng)钍夏Ａ?τxy表示剪應(yīng)力；v代表泊松比。將式(7)代入式(8)可得曲率和應(yīng)力的關(guān)系式為

(9)

至此在求得各方向的應(yīng)力大小之后，可進(jìn)一步求得目的層段的最大主應(yīng)力大小及方向?yàn)?/p>

(10)

(11)

式(10)、(11)中：σmax為主應(yīng)力;α為主應(yīng)力與x軸的夾角。

因此，當(dāng)我們計(jì)算得到目的層段的地層傾角p和q之后，可進(jìn)一步基于曲率-應(yīng)變-應(yīng)力的關(guān)系得到地層傾角-應(yīng)變-應(yīng)力之間的關(guān)系，從而可以求得目的層段地層各個(gè)方向構(gòu)造應(yīng)力的大小，并進(jìn)一步計(jì)算地層主應(yīng)力的大小及方向。

1.2 基于離心窗傾角掃描的地層傾角計(jì)算方法

印興耀 等[6]提出的離心窗掃描技術(shù)可以直接從地震資料中提取地層傾角，該方法將分析點(diǎn)周圍的4個(gè)點(diǎn)依次作為中心點(diǎn)在x、y、t方向上開一個(gè)窗口掃描地震波形的相似性,掃描窗口在x、y、t方向上的延伸長度需要根據(jù)實(shí)際情況確定。離心窗掃描的計(jì)算方法如下：

Cij(t,p,q)=

(12)

在確定了4個(gè)窗口在x、y、t方向上的延伸長度后，通過不斷調(diào)整p和q的大小使相似性達(dá)到最大，這時(shí)的p和q即為該分析點(diǎn)的地層傾角。

1.3 基于傾角屬性的地應(yīng)力場(chǎng)計(jì)算流程

在通過彈性力學(xué)的變形幾何方程和薄板假設(shè)理論推導(dǎo)出的地層地應(yīng)力場(chǎng)方程即式(9)中，楊氏模量E和泊松比v可以利用地層的縱波速度(Vp)、橫波速度(Vs)和密度(ρ)求解[13]，即

類似于疊后波阻抗反演中的速度建模方法[14]，選用多井層約束方法進(jìn)行楊氏模量和泊松比的計(jì)算。

通過離心窗掃描獲得地層傾角屬性以及上述方法獲得地層楊氏模量和泊松比，就可以計(jì)算地層的應(yīng)力場(chǎng)，具體步驟如下：

1) 通過離心窗掃描目的層段的地震數(shù)據(jù)，計(jì)算地層傾角p和q;

2) 基于測(cè)井統(tǒng)計(jì)的縱、橫波速度及密度參數(shù)，計(jì)算各井的楊氏模量和泊松比大小，并基于多井楊氏模量和泊松比參數(shù)進(jìn)行目的層段的多井楊氏模量和泊松比建模;

3) 將研究區(qū)各點(diǎn)的傾角屬性和楊氏模量代入傾角-應(yīng)力方程即式(9)，計(jì)算目的層段各方向的構(gòu)造應(yīng)力大小，并進(jìn)一步計(jì)算主應(yīng)力的大小及方向；

4) 分析統(tǒng)計(jì)主應(yīng)力大小及方向。

2 實(shí)例應(yīng)用

將本文提出的基于傾角屬性的三維地應(yīng)力場(chǎng)分析方法應(yīng)用于渤海某工區(qū)。該工區(qū)為大型花崗巖潛山油氣藏區(qū)，風(fēng)化殼儲(chǔ)層在形成過程中受到區(qū)域構(gòu)造應(yīng)力的作用形成大量內(nèi)幕小斷層[15]，再由風(fēng)化淋濾作用與風(fēng)化剝蝕作用形成溶蝕縫和風(fēng)化破碎帶。裂縫系統(tǒng)是花崗巖潛山油氣藏重要的油氣運(yùn)移通道[16]，所以研究地應(yīng)力對(duì)該地區(qū)的油氣開發(fā)具有重要意義。

圖2為本文得到的研究區(qū)應(yīng)力場(chǎng)方向示意圖，每個(gè)小玫瑰圖為剖分網(wǎng)格內(nèi)的主應(yīng)力發(fā)育方向，短線顏色代表方位角，藍(lán)色代表低值，紅色代表高值。圖3為研究區(qū)測(cè)井方法得到的井點(diǎn)位置地應(yīng)力玫瑰圖。對(duì)比圖2、3中井點(diǎn)位置處應(yīng)力的大小與方向可以發(fā)現(xiàn)，除P2井外，2種方法結(jié)果一致。分析認(rèn)為P2井位于工區(qū)邊緣，可能受資料處理邊緣效應(yīng)的影響。 研究區(qū)裂縫發(fā)育與應(yīng)力大小有密切關(guān)系，按照應(yīng)力值的大小可以將工區(qū)分為3個(gè)區(qū)域：應(yīng)力值在300～1 600 N/m2范圍內(nèi)為第1級(jí)強(qiáng)應(yīng)力區(qū)，在0～300 N/m2范圍內(nèi)為第2級(jí)中應(yīng)力區(qū)，在-300～0 N/m2范圍內(nèi)為第3級(jí)弱應(yīng)力區(qū)。圖4為本文方法得到的研究區(qū)應(yīng)力場(chǎng)值分段統(tǒng)計(jì)分析結(jié)果，可以看出：強(qiáng)應(yīng)力主要分布在工區(qū)西南部的2個(gè)條帶狀區(qū)域內(nèi)；中應(yīng)力為該研究區(qū)主要應(yīng)力，呈條帶狀分布于整個(gè)區(qū)域；與強(qiáng)應(yīng)力和中應(yīng)力不同，弱應(yīng)力連片分布于研究區(qū)若干個(gè)區(qū)域內(nèi)，這也符合潛山溝谷的構(gòu)造特征。研究區(qū)中有一些區(qū)域位于多方向應(yīng)力交會(huì)的區(qū)域，這些地方應(yīng)力場(chǎng)比較復(fù)雜，更容易發(fā)育裂縫。

圖2 本文方法得到的研究區(qū)應(yīng)力場(chǎng)方向示意圖Fig .2 Orientation diagram of stress field in the study area by the present approach

圖3 研究區(qū)井點(diǎn)位置地應(yīng)力玫瑰圖Fig .3 Stress rose diagram of each well in the study area

圖4 本文方法得到的研究區(qū)應(yīng)力場(chǎng)值分段統(tǒng)計(jì)分析結(jié)果Fig .4 Statistical analysis of different distribution of stress field in the study area by the present approach

3 結(jié)論

1) 本文方法根據(jù)地層傾角和曲率的關(guān)系，將曲率-應(yīng)變-應(yīng)力方程變形為傾角-應(yīng)變-應(yīng)力方程，實(shí)現(xiàn)了通過地層傾角計(jì)算地應(yīng)力。

2) 本文方法借用了三維體曲率分析的思路，直接針對(duì)三維地震數(shù)據(jù)體進(jìn)行計(jì)算，有效克服了常規(guī)方法受層位解釋精度影響的缺陷，可以提高地應(yīng)力場(chǎng)計(jì)算結(jié)果的準(zhǔn)確性。

3) 渤海某工區(qū)花崗巖儲(chǔ)層實(shí)際資料的計(jì)算結(jié)果表明，本文方法的預(yù)測(cè)結(jié)果與測(cè)井解釋的一致性較好，能夠較為準(zhǔn)確地建立地層構(gòu)造應(yīng)力場(chǎng)，可以有效地應(yīng)用于研究區(qū)的裂縫預(yù)測(cè)，并且具有一定的推廣意義。

[1] 黃保綱,趙春明,楊慶紅,等.等效應(yīng)力場(chǎng)模擬與疊前彈性波阻抗反演綜合預(yù)測(cè)錦州25-1南潛山裂縫儲(chǔ)層[J].中國海上油氣,2012,24(1):17-20.HUANG Baogang,ZHAO Chunming,YANG Qinghong,et al.Predicting fractured reservoirs in Jinzhou25-1S buried hill by comprehensively using equivalent stress-field simulation and pre-stack elastic impedance inversion [J].China Offshore Oil and Gas,2012,24(1):17-20.

[2] 賈立宏,李造鼎.地應(yīng)力場(chǎng)灰色模擬方法研究[J].東北大學(xué)學(xué)報(bào)(自然科學(xué)版),1995,16(6):559-563.JIA Lihong,LI Zaoding.Methodical studies on grey simulation of In-Situ stress field [J].Journal of Northeastern University (Natural Science),1995,16(6):559-563.

[3] 張帆，賀振華，黃德濟(jì)，等.預(yù)測(cè)裂隙發(fā)育帶的構(gòu)造應(yīng)力場(chǎng)數(shù)值模擬技術(shù)[J].石油地球物理勘探,2000,35(2):154-163.ZHANG Fan,HE Zhenhua,HUANG Deji,et al.Structural stress field numerical simulation technique for fracture zone prediction [J].OGP,2000,35(2):154-163.

[4] 張國強(qiáng)，王桂萱.基于神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)分界與FLAC的初始地應(yīng)力場(chǎng)反演[J].大連大學(xué)學(xué)報(bào)，2007,28(6):43-47.ZHANG Guoqiang,WANG Guixuan.Back-analysis of initial ground stress field based on neural network ensemble and generation of the initial ground stress field by using FLAC [J].Journal of Dalian University,2007,28(6):43-47.

[5] 何英.高精度曲率分析及其在構(gòu)造識(shí)別中的應(yīng)用[D].成都：成都理工大學(xué),2011.HE Ying.High precision curvature analysis and its application of structural identification [D].Chengdu:Chengdu University of Technology,2011.

[6] 印興耀，高京華，宗兆云.基于離心窗傾角掃描的曲率屬性提取[J].地球物理學(xué)報(bào),2014,57(10)：3411-3421.YIN Xingyao,GAO Jinghua,ZONG Zhaoyun.Curvature attribute based on dip scan with eccentric window[J].Chinese Journal of Geophysics,2014,57(10):3411-3421.

[7] 楊威，賀振華，陳學(xué)華.三維體曲率屬性在斷層識(shí)別中的應(yīng)用[J].地球物理學(xué)進(jìn)展,2011,26(1):110-115.YANG Wei,HE Zhenhua,CHEN Xuehua.Application of three-dimensional volumetric curvature attributes to fault identification [J].Progress in Geophysics,2011,26(1):110-115.

[8] 黃炎.彈性薄板理論[M].北京:國防科技大學(xué)出版社,1992:10-50.

[9] MARFURT K J.Robust estimates of 3D reflector dip and azimuth[J].Geophysics,2006,71(4):29-40.

[10] Al-DOSSSARY S,MARFURT K J.3D volumetric multispectral estimates of reflector curvature and rotation [J].Geophysics,2006,71(5):41-51.

[11] ROBERT A.Curvature attributes and their application to 3D interpreted horizons [J].First Break,2001,19(2):85-100.

[12] COOKE D A,SCHNEIDER W A.Generalized linear inversion of reflection seismic data [J].Geophysics,1983,48(6):665-676.

[13] STEWART S A,WYNN T J.Mapping spatial variation in rock properties in relationship to scale-dependent structure using spectral curvature [J].Geology,2000,28(3):691-694.

[14] BERGBAUER S,MUKERJI T,HENNINGS P.Improving curvature analyses of deformed horizons using scale-dependent filtering techniques [J].AAPG Bulletin,2003,87(2):1255-1272.

[15] 譚成軒,王連捷,孫寶珊,等.含油氣盆地三維構(gòu)造應(yīng)力場(chǎng)數(shù)值模擬方法[J].地質(zhì)力學(xué)學(xué)報(bào)，1997，3(1)：71-79.TAN Chengxuan,WANG Lianjie,SUN Baoshan,et al.An approach to numerical simulation of 3-D tectonic stress field of the oil-gas-bearing basin [J].Journal of Geomechanics,1997,3(1):71-79.

[16] 王雷，陳海清，陳國文,等.應(yīng)用曲率屬性預(yù)測(cè)裂縫發(fā)育帶及其產(chǎn)狀[J].石油地球物理勘探，2010,45(6):885-889.WANG Lei,CHEN Haiqing,CHEN Guowen,et al.Application of curvature attributes in predicting fracture-developed zone and its orientation [J].OGP,2010,45(6):885-889.