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

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

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

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

1 方法原理

1.1 基于地層傾角的地應力計算方法

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

(1)

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

(2)

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

(3)

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

(4)

(5)

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

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

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

(6)

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

(7)

由廣義胡克定律，進一步得到應力與應變關系式

(8)

式(8)中:σx和σy分別表示x和y方向的應力;E代表楊氏模量;τxy表示剪應力；v代表泊松比。將式(7)代入式(8)可得曲率和應力的關系式為

(9)

至此在求得各方向的應力大小之后，可進一步求得目的層段的最大主應力大小及方向為

(10)

(11)

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

因此，當我們計算得到目的層段的地層傾角p和q之后，可進一步基于曲率-應變-應力的關系得到地層傾角-應變-應力之間的關系，從而可以求得目的層段地層各個方向構造應力的大小，并進一步計算地層主應力的大小及方向。

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

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

Cij(t,p,q)=

(12)

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

1.3 基于傾角屬性的地應力場計算流程

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

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

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

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

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

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

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

2 實例應用

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

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

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

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

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

3 結論

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

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

3) 渤海某工區(qū)花崗巖儲層實際資料的計算結果表明，本文方法的預測結果與測井解釋的一致性較好，能夠較為準確地建立地層構造應力場，可以有效地應用于研究區(qū)的裂縫預測，并且具有一定的推廣意義。

[1] 黃保綱,趙春明,楊慶紅,等.等效應力場模擬與疊前彈性波阻抗反演綜合預測錦州25-1南潛山裂縫儲層[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] 賈立宏,李造鼎.地應力場灰色模擬方法研究[J].東北大學學報(自然科學版),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] 張帆，賀振華，黃德濟，等.預測裂隙發(fā)育帶的構造應力場數(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] 張國強，王桂萱.基于神經網絡結構分界與FLAC的初始地應力場反演[J].大連大學學報，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] 何英.高精度曲率分析及其在構造識別中的應用[D].成都：成都理工大學,2011.HE Ying.High precision curvature analysis and its application of structural identification [D].Chengdu:Chengdu University of Technology,2011.

[6] 印興耀，高京華，宗兆云.基于離心窗傾角掃描的曲率屬性提取[J].地球物理學報,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] 楊威，賀振華，陳學華.三維體曲率屬性在斷層識別中的應用[J].地球物理學進展,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].北京:國防科技大學出版社,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] 譚成軒,王連捷,孫寶珊,等.含油氣盆地三維構造應力場數(shù)值模擬方法[J].地質力學學報，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] 王雷，陳海清，陳國文,等.應用曲率屬性預測裂縫發(fā)育帶及其產狀[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.