周進(jìn)舉,王德利,靳 松,馮 飛

(1.吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林長春130026;2.中國石油大學(xué)(北京)物探重點(diǎn)實(shí)驗(yàn)室，北京102249)

頁巖作為一種典型的橫向各向同性介質(zhì)，不同于傳統(tǒng)的地震勘探中各向同性的地震模型，在地震勘探中必須考慮其各向異性的特點(diǎn)，建立相應(yīng)的地質(zhì)模型。另外，在實(shí)際地質(zhì)條件下，頁巖中常常含有垂直或傾斜裂隙，而裂隙受多種因素控制，有很明顯的各向異性特征，同時(shí)裂隙也是頁巖氣儲(chǔ)存的重要空間和運(yùn)移的重要通道[1]。因此，若對(duì)頁巖氣儲(chǔ)層的地質(zhì)條件進(jìn)行模擬，我們應(yīng)當(dāng)考慮含有裂隙的各向異性模型，而研究的基礎(chǔ)就是獲得各種裂隙各向異性等效模型的彈性參數(shù)。

對(duì)于含裂隙介質(zhì)彈性系數(shù)的計(jì)算，Hudson[2-4]提出了微裂隙理論，用這個(gè)理論可以計(jì)算出廣泛擴(kuò)容性各向異性(Extensive Dilatancy Anisotropy，簡稱EDA)介質(zhì)的彈性常數(shù);Schoenberg等[5-6]提出的線性滑動(dòng)理論解決了方位各向異性介質(zhì)和含垂直裂隙橫向各向同性介質(zhì)的等效彈性系數(shù)的計(jì)算問題;2002年，傅旦丹等[7]根據(jù)Backus方法并結(jié)合Hudson微裂隙理論，推出了周期性薄互層(Periodic Thin Layer，簡稱PTL)與EDA組合等效正交各向異性介質(zhì)彈性常數(shù)的計(jì)算公式;同年，王德利等[8]根據(jù)Hudson等關(guān)于裂隙介質(zhì)彈性系數(shù)計(jì)算的擾動(dòng)理論和Bond變換矩陣原理，計(jì)算了含兩組斜交的垂直裂隙形成的單斜各向異性介質(zhì)的等效彈性系數(shù)。

對(duì)于含垂直裂隙的頁巖介質(zhì)，可以建立PTL與EDA的組合模型，這種組合模型在長波長情況下已經(jīng)證明具有正交各向異性的性質(zhì);同樣，當(dāng)裂隙傾斜時(shí)，可以建立一個(gè)PTL+傾斜EDA模型，在長波長情況下，即波長遠(yuǎn)大于薄互層層厚時(shí)，這種介質(zhì)是否具有傾斜橫向各向同性(TTI)介質(zhì)的性質(zhì)就成了值得研究的問題。目前，國內(nèi)外針對(duì)TTI介質(zhì)的正演模擬做了很多研究工作。吳國忱等[9]作了TTI介質(zhì)彈性波頻率-空間域有限差分模擬，并分析了TTI介質(zhì)彈性系數(shù)變化時(shí)的波形變化規(guī)律。祝賀君等[10]用高階同位網(wǎng)格有限差分法對(duì)TTI介質(zhì)地震波場進(jìn)行了模擬;Saenger等[11]用旋轉(zhuǎn)網(wǎng)格有限差分法對(duì)各向異性介質(zhì)進(jìn)行了正演模擬;Han等[12]對(duì)傾斜橫向各向同性介質(zhì)進(jìn)行了彈性波模擬。以上的研究都是直接旋轉(zhuǎn)VTI介質(zhì)彈性系數(shù)得到TTI介質(zhì)彈性系數(shù)，然后進(jìn)行正演模擬，他們的模擬結(jié)果可以作為對(duì)照，來證明PTL+傾斜EDA模型是否具有TTI各向異性特征。

本文在前人對(duì)含裂隙介質(zhì)彈性系數(shù)計(jì)算的研究基礎(chǔ)上，針對(duì)PTL+傾斜EDA組合模型，利用Hudson理論計(jì)算出每層EDA介質(zhì)的彈性系數(shù)，然后通過Bond變換得到傾斜EDA介質(zhì)的彈性系數(shù)，在長波長的情況下，用Backus加權(quán)平均法計(jì)算了整體的等效彈性系數(shù)。進(jìn)而用高階旋轉(zhuǎn)網(wǎng)格有限差分法對(duì)該模型的地震波場響應(yīng)進(jìn)行了正演模擬，證明了所給出的含傾斜裂隙頁巖介質(zhì)彈性系數(shù)計(jì)算方法的正確性。模擬結(jié)果與TTI介質(zhì)的模擬結(jié)果相對(duì)照，證明了PTL+傾斜EDA組合模型，即含傾斜裂隙頁巖介質(zhì)，具有和TTI介質(zhì)類似的各向異性特征。另外，這種方法還可以用于研究含傾斜裂隙薄互層介質(zhì)中裂隙傾角、裂隙密度和裂隙填充物性質(zhì)對(duì)地震波傳播的影響。

1 彈性系數(shù)的計(jì)算

頁巖介質(zhì)是典型的PTL介質(zhì)，在PTL模型中加入傾斜裂隙，即可很好地模擬含傾斜裂隙頁巖介質(zhì)。我們把這種模型稱為PTL+傾斜EDA介質(zhì)，如圖1所示。

圖1 PTL+傾斜EDA介質(zhì)

對(duì)于PTL介質(zhì)的每一層都可以看作是傾斜的EDA介質(zhì)，則每一層的彈性系數(shù)都可以用計(jì)算EDA介質(zhì)的方法求得，然后再進(jìn)行Bond變換，傾斜任意角度。這樣PTL+傾斜EDA介質(zhì)就可以看作周期性的傾斜EDA介質(zhì)。然后采用加權(quán)平均法求得整體的等效彈性系數(shù)。

1.1 EDA介質(zhì)彈性系數(shù)的計(jì)算

EDA介質(zhì)如圖2所示，計(jì)算其彈性系數(shù)可以使用Hudson理論。

圖2 EDA介質(zhì)

Hudson給出的二階攝動(dòng)形式為

(1)

經(jīng)過推導(dǎo)，含定向排列裂隙介質(zhì)的彈性系數(shù)矩陣為

(2)

其中，C23=C33-2C44，因此描述EDA介質(zhì)的彈性系數(shù)有5個(gè)[7]。

1.2 Bond變換

模型每一層的EDA的裂隙面都不垂直于x軸，其方向相當(dāng)于繞y軸在xoz面旋轉(zhuǎn)了θ角。所以在計(jì)算各層的彈性系數(shù)加權(quán)之前，需要對(duì)每一層進(jìn)行Bond變換，變換矩陣為Mθ[13]。原坐標(biāo)系下的各層彈性系數(shù)矩陣(2)經(jīng)Bond變換后，其形式為

(3)

1.3 加權(quán)平均

根據(jù)Backus加權(quán)平均法[7]和公式(3)以及附錄A中求得的應(yīng)力與應(yīng)變的本構(gòu)關(guān)系，可以得到對(duì)應(yīng)的等效彈性系數(shù)矩陣:

(4)

其中，

(4)式給出了PTL+傾斜EDA模型的彈性參數(shù)計(jì)算表達(dá)式，而符號(hào)“〈〉”表示在垂向(z軸方向)上的加權(quán)平均。根據(jù)White的理論[14]可以知道，對(duì)各薄層彈性參數(shù)來說，其對(duì)應(yīng)權(quán)函數(shù)的值正比于此薄層在整個(gè)地質(zhì)體模型中的厚度的比例。因此，當(dāng)知道PTL介質(zhì)各層的縱、橫波速度，厚度，密度以及裂隙密度、裂隙縱橫比、裂隙填充物后，即可用(4)式求其對(duì)應(yīng)的13個(gè)彈性參數(shù)。

以上的推導(dǎo)過程中，裂隙的方位角為0，當(dāng)方位角不為0時(shí)，即觀測坐標(biāo)系與本征坐標(biāo)系之間的夾角為φ，上述彈性系數(shù)矩陣還要再旋轉(zhuǎn)一次，旋轉(zhuǎn)矩陣為Mφ，旋轉(zhuǎn)后的彈性系數(shù)矩陣為

(5)

這種情況下，共有21個(gè)獨(dú)立的彈性系數(shù)，與三斜介質(zhì)的彈性系數(shù)個(gè)數(shù)相同，此時(shí)，準(zhǔn)P波、準(zhǔn)SH波和準(zhǔn)SV波耦合[15]。

2 模擬算例

本文在驗(yàn)證上述含傾斜裂隙頁巖介質(zhì)彈性系數(shù)計(jì)算方法的正確性時(shí)，采用Saenger和Bohlen提出的旋轉(zhuǎn)網(wǎng)格有限差分法，把速度定義在網(wǎng)格單元頂點(diǎn)，應(yīng)力和彈性系數(shù)定義在網(wǎng)格單元中心，這樣在計(jì)算應(yīng)力時(shí)的所有空間偏導(dǎo)數(shù)可以直接求解[16]；采用的吸收邊界條件是Cerjan等[17]提出的簡單吸收邊界。網(wǎng)格點(diǎn)數(shù)600×600，網(wǎng)格間距5m，時(shí)間步長1ms，爆炸震源，主頻30Hz。

如圖1所示，方位角(φ)為0，此時(shí)觀測方向與裂隙走向一致，在這種情況下接收到的y分量始終為0，因此，只給出了裂隙傾角(θ)分別為30°,45°和60°時(shí)波場快照中的x分量和z分量。圖3到圖5分別為裂隙傾角為30°,45°和60°時(shí)的波場快照，可以看出,縱波波前面的長軸方向與介質(zhì)的裂隙方向一致,隨著裂隙傾角的旋轉(zhuǎn)而旋轉(zhuǎn);另外，波場快照中都可以觀察到橫波分裂現(xiàn)象，產(chǎn)生振幅奇異性，符合吳國忱等[9]和牟永光等[18]文獻(xiàn)的論述。由此證明了本文所給出的彈性系數(shù)計(jì)算方法的正確性和PTL+傾斜EDA模型具有傾角各向異性特征。

圖3 含傾斜裂隙(傾角為30°)PTL介質(zhì)波場快照a x分量; b z分量

圖4 含傾斜裂隙(傾角為45°)PTL介質(zhì)波場快照a x分量; b z分量

圖5 含傾斜裂隙(傾角為60°)PTL介質(zhì)波場快照a x分量; b z分量

當(dāng)裂隙方位角為30°和45°，傾角為30°時(shí)，正演模擬結(jié)果分別如圖6和圖7所示。圖6和圖7中P波波前面為橢圓形，橫波波前面出現(xiàn)了明顯的橫波分裂現(xiàn)象。y分量不再為0，出現(xiàn)了SH波;在x分量和z分量上出現(xiàn)SV波，而且SV波和SH波的波形和相位都有所不同。另外，圖6和圖7的對(duì)比結(jié)果顯示方位角的變化對(duì)P波波前面的影響不大，但是會(huì)明顯影響S波波前面，當(dāng)方位角為45°時(shí)，在垂直于對(duì)稱軸方向橫波不發(fā)生分裂，在對(duì)稱軸方向出現(xiàn)明顯快橫波。圖3到圖5中裂隙方位角為0時(shí)，縱波波前面呈橢圓狀，兩段比較尖銳;而圖6和圖7中，縱波波前面越來越圓。以上這些特點(diǎn)證明了PTL+傾斜EDA介質(zhì)模型具有方位各向異性特征。

當(dāng)裂隙含氣，裂隙方位角為45°，裂隙傾角為30°時(shí)的波場快照如圖8所示。對(duì)比裂隙飽含水和干燥時(shí)的波場快照(圖7和圖8)可以看出:干燥裂隙，即裂隙含氣時(shí)，縱波波前面呈明顯的橢圓形，橫波波前面具有更強(qiáng)的能量，說明裂隙含氣時(shí)對(duì)地震波的傳播有顯著影響，介質(zhì)具有更強(qiáng)的各向異性特征。這也說明本文所給出的彈性系數(shù)計(jì)算方法還可用于研究含裂隙頁巖介質(zhì)中裂隙傾角、裂隙密度和裂隙填充物性質(zhì)對(duì)地震波傳播的影響。

圖6 含傾斜裂隙(方位角為30°，傾角為30°)PTL介質(zhì)波場快照a x分量; b y分量; c z分量

圖7 含傾斜裂隙(方位角為45°，傾角為30°)PTL介質(zhì)波場快照a x分量; b y分量; c z分量

圖8 含氣傾斜裂隙(方位角為45°，傾角為30°)PTL介質(zhì)波場快照a x分量; b y分量; c z分量

3 結(jié)束語

1) 正演模擬結(jié)果證明了本文根據(jù)Hudson微裂隙理論、Bond變換和Backus加權(quán)平均法推導(dǎo)出的含傾斜裂隙頁巖介質(zhì)模型等效彈性系數(shù)計(jì)算方法的正確性，同時(shí)也說明了Backus加權(quán)平均法在計(jì)算薄互層介質(zhì)等效彈性系數(shù)時(shí)的有效性。

2) 正演模擬波場快照中橢圓形波形對(duì)稱軸隨裂隙傾角旋轉(zhuǎn)的情況和隨裂隙方位角不同而不同的橫波分裂現(xiàn)象，證明了含傾斜裂隙頁巖介質(zhì)模型具有和TTI介質(zhì)類似的傾角各向異性和方位各向異性特征。

3) 因?yàn)楸疚慕o出的含傾斜裂隙頁巖介質(zhì)彈性系數(shù)由介質(zhì)各層的縱、橫波速度，厚度，密度以及裂隙密度、裂隙縱橫比、裂隙填充物等參數(shù)直接求出，因此本文方法還可應(yīng)用于研究薄互層頁巖介質(zhì)中裂隙傾角、裂隙密度和裂隙填充物性質(zhì)對(duì)地震波傳播的影響。

附錄A

各層的應(yīng)力和應(yīng)變的本構(gòu)關(guān)系有

(1)

把σxx，σyy，σxy，uxz，uyz，uzz用σxz，σyz，σzz，uξx，uξy(ξ=x,y,z)表示，由(1)式可得

(2)

其中，

(3)

根據(jù)加權(quán)平均函數(shù)的特殊性質(zhì)[7]，對(duì)(2)式進(jìn)行加權(quán)平均運(yùn)算，并再次變換本構(gòu)方程，使應(yīng)力用應(yīng)變表示:

(4)

其中，

(5)

根據(jù)以上求得的應(yīng)力與應(yīng)變的本構(gòu)關(guān)系，即可得到對(duì)應(yīng)的等效彈性系數(shù)矩陣。

參 考 文 獻(xiàn)

[1] 李新景,胡素云,程克明.北美裂縫性頁巖氣勘探開發(fā)的啟示[J].石油勘探與開發(fā),2007,34(4):392-400

Li J X,Hu S Y,Cheng K M.Suggestions from the development of fractured shale gas in North America[J].Petroleum Exploration and Development,2007,34(4):392-400

[2] Hudson J A.Wave speeds and attenuation of elastic waves in material containing cracks[J].Geophysical Journal of the Royal Astronomical Society,1981,64(1):133-150

[3] Hudson J A.A higher order approximation to the wave propagation constants for a cracked solid [J].Geophysical Journal of the Royal Astronomical Society,1986,87(1):265-274

[4] Hudson J A.Crack distributions which account for a given seismic anisotropy[J].Geophysical Journal International,1991,104(3):517-521

[5] Schoenberg M,Douma J.Elastic wave propagation in media with parallel fractures and aligned cracks[J].Geophysical Prospecting,1988,36(6):571-590

[6] Schoenberg M,Helbig K.Orthorhombic media:modeling elastic wave behavior in a vertically fractured earth[J].Geophysics,1997,62(6):1954-1974

[7] 傅旦丹,何樵登.含垂直裂解隙PTL介質(zhì)長波長等效正交各向異性彈性常數(shù)[J].地球物理學(xué)報(bào),2002,45(增刊):307-320

Fu D D,He Q D.Elastic constants of long-wavelength effective orthorhombic anisotropy for PTL media containing vertical cracks[J].Chinese Journal of Geophysics(in Chinese),2002,45(S1):307-320

[8] 王德利,何樵登.裂隙型單斜介質(zhì)中彈性系數(shù)的計(jì)算及波的傳播特性研究[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),2002,32(2):91-96

Wang D L,He Q D.Elastic constants calculation method and propagation properties for cracked monoclinic media[J].Journal of Jilin University(Earth Science Edition)，2002,32(2):91-96

[9] 吳國忱,羅彩明,梁楷.TTI介質(zhì)彈性波頻率-空間域有限差分?jǐn)?shù)值模擬[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),2007,37(5):1023-1033

Wu G C,Luo C M,Lang K.Frequency-space domain finite difference numerical simulation of elastic wave in TTI media[J].Journal of Jilin University (Earth Science Edition),2007,37(5):1023-1033

[10] 祝賀君,張偉,陳曉非.二維各向異性介質(zhì)中地震波

場的高階同位網(wǎng)格有限差分模擬[J].地球物理學(xué)報(bào),2009,52(6):1536-1546

Zhu H J,Zhang W,Chen X F.Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method[J].Chinese Journal of Geophysics(in Chinese),2009,52(6):1536-1546

[11] Saenger E H,Bohlen T.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J].Geophysics,2004,69(2):583-591

[12] Han B,Seol S J,Byun J.Elastic modelling in tilted transversely isotropic media with convolutional perfectly matched layer boundary conditions[J].Exploration Geophysics,2012,43(2):77-86

[13] 孫瑞艷.TTI介質(zhì)旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分及其組合邊界條件[D].青島:中國石油大學(xué)(華東),2010

Sun R Y.Rotated staggered grid finite-difference and the combined boundary conditions in TTI media[D].Qingdao:China University of Petroleum(EastChina),2010

[14] White J E.地下的聲音-地震波的應(yīng)用[M]. 陸其鵠，齊國英,譯.北京:地震出版社，1986:1-167

White J E.Underground sound:application of seismic wave[M].Lu Q H,Qi G Y,Translator.Beijing: Seismological Press,1986:1-167

[15] 郝奇.三維TTI介質(zhì)地震波場正演模擬[D].長春:吉林大學(xué),2007

Hao Q.Seismic wavefield modeling in 3D TTI media[D].Changchun:Jilin University,2007

[16] 嚴(yán)紅勇,劉洋.黏彈TTI介質(zhì)中旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)報(bào),2012,55(4):1354-1365

Yan H Y,Liu Y.Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media[J].Chinese Journal of Geophysics(in Chinese),2012,55(4):1354-1365

[17] Cerjan C,Kosloff D,Kosloff R,et al.A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics,1985,50(4):705-708

[18] 牟永光，裴正林.三維復(fù)雜介質(zhì)地震數(shù)值模擬[M].北京:石油工業(yè)出版社，2005:75-86

Mu Y G,Pei Z L.Seismic numerical modeling for 3-D complex media[M].Beijing:Petroleum Industry Press,2005:75-86