李雨生 吳國忱

（山東青島266580中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院）

?

一種三維正交方位各向異性介質(zhì)巖石物理建模及彈性波正演模擬方法

（山東青島266580中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院）

通過線性滑動理論和巖石物理等效理論， 將兩組正交直立裂隙介質(zhì)等效為一種正交方位各向異性介質(zhì)進行三維巖石物理建模， 通過高階交錯網(wǎng)格有限差分求解彈性波動方程模擬地震波在該種介質(zhì)中的傳播過程． 在建模過程中改變物性參數(shù)， 分析不同裂隙密度條件下的炮集和波場特征， 以及正交各向異性的方位特征． 研究結(jié)果表明， 各向異性強度隨裂隙密度等物性增大而增強， 而且這些特征在共炮點道集和波場中均有所體現(xiàn)．

正交方位各向異性 裂隙介質(zhì) Hudson理論 線性滑動理論 交錯網(wǎng)格有限差分

引言

地球介質(zhì)廣泛存在波動各向異性， 地震各向異性主要表現(xiàn)在地震波傳播速度是傳播方向的函數(shù)、 體波間相互耦合、 橫波發(fā)生分裂等（吳國忱， 2006）． 近幾十年地震各向異性研究取得了很大的進展， 特別是Backus（1962）研究表明周期性薄互層可引起視各向異性， Crampin證實了裂隙誘導(dǎo)各向異性和橫波分裂的存在， 并于1987年提出了廣泛擴容各向異性（extensive dilatancy anisotropy， 簡寫為EDA）模型（Crampin， 1978， 1981， 1987）． 裂隙介質(zhì)不僅是油氣運移的重要通道， 還是重要的油氣存儲空間， 因此在油氣勘探領(lǐng)域越來越受到重視（金抒辛， 何樵登， 2005； 吳國忱， 秦海旭， 2014）． 裂隙中彈性波在傳播過程中呈現(xiàn)出方位各向異性特征． 具有方位各向異性特征的介質(zhì)包括HTI （horizontal transverse isotropic）介質(zhì)、 正交各向異性介質(zhì)和單斜各向異性介質(zhì)等， 本文著重將兩組正交直立裂隙等效為一種正交各向異性介質(zhì)． 雖然早在20世紀(jì)80年代初， Crampin等（1980）就提出雙平面裂隙系統(tǒng)， 但后續(xù)的相關(guān)研究只側(cè)重于在薄互層背景上發(fā)育有垂直定向排列的裂隙（Schoenberg， Hilbig, 1997； 張文生， 宋海濱， 2001； 何燕， 2008）．

線性滑動理論（Schoenberg， 1980； Coates， Schoenberg， 1995； Schoenberg， Muir， 1989； Schoenberg， Sayers， 1995; Bakulinetal， 2000a， b， c）將裂隙看作無限的、 很薄的且具有很強屈服性的層面， 或?qū)⑵淇醋骶哂芯€性滑動邊界條件的柔性平面． 在此基礎(chǔ)上， 可將裂隙模型等效為正交各向異性來進行巖石物理建模， 為地震波數(shù)值模擬提供速度模型．

地震波正演模擬是模擬地震波在地球介質(zhì)中的傳播過程， 并研究地震波的傳播特性與地球介質(zhì)參數(shù)的關(guān)系（孫成禹， 2007）， 波動方程數(shù)值模擬是正演模擬的一種非常重要的方法． 數(shù)值模擬是通過數(shù)值計算達(dá)到對實際地震剖面的最優(yōu)逼近， 是應(yīng)用數(shù)值方法研究地下儲層地震反射特征的基礎(chǔ)工具． 常見的波動方程數(shù)值解法包括有限差分法和有限元法等． 有限差分法是將波動方程中的介質(zhì)參數(shù)和波場函數(shù)離散化， 以差分算子代替微分算子， 在有限精度內(nèi)對地震波傳播問題進行模擬． Virieux（1987）給出了波動方程一階速度-應(yīng)力方程并用交錯網(wǎng)格模擬地震波在非均勻介質(zhì)中的傳播過程． 交錯網(wǎng)格差分法相對規(guī)則網(wǎng)格差分法, 其網(wǎng)格頻散顯著減小， 精度明顯提高， 而且可以取較大的空間步長， 提高計算效率（董良國等， 2000）． 利用有限差分法模擬地震波的傳播過程主要有3個需要注意的問題， 即邊界條件處理、 頻散壓制和差分穩(wěn)定性條件． 其中邊界條件的處理本文采用最佳匹配層（perfectly matched layer， 簡寫為PML）法． 該方法由Berenger（1994）針對電磁波傳播情況提出， Collino和Tsogka（2001）將其應(yīng)用于各向異性介質(zhì)地震波數(shù)值模擬． Boris和Book（1973）將求解流體動力學(xué)問題時的通量校正傳輸方法應(yīng)用于聲波方程求解， 實現(xiàn)了聲波方程差分計算中數(shù)值頻散的有效壓制．

本文擬基于線性滑動理論結(jié)合Hudson理論， 將兩組垂直正交裂隙介質(zhì)等效為正交方位各向異性介質(zhì)， 利用多組物性參數(shù)， 計算得到不同裂隙密度下的彈性參數(shù)， 并以此來建立三維模型． 通過有限差分正演模擬彈性波在這些介質(zhì)中的傳播過程， 分析不同物性條件下共炮點道集和波場的特征及其各向異性的方位特征．

1 基于線性滑動理論結(jié)合Hudson理論的巖石物理建模

Schoenberg（1980）指出， 可以忽略裂隙的形狀而將其看作滿足線性滑動邊界條件的柔性平面． 所謂線性滑動， 是指兩個彈性固體接觸面上的應(yīng)變是不連續(xù)的， 但引起應(yīng)變的應(yīng)力在接觸面上是連續(xù)的， 且應(yīng)力與應(yīng)變成線性關(guān)系．

正交裂隙是實際存在的一種重要正交介質(zhì)， 是各向同性介質(zhì)背景或VTI（vertical transverse isotropic）介質(zhì)背景中兩組正交的直立裂隙． 圖1給出了上述正交介質(zhì)的實際發(fā)育情況及相應(yīng)模型. 可以看出， 該正交介質(zhì)由兩組正交的直立裂隙組成．

圖1 正交介質(zhì)模型， 圖中紅線和虛線均表示裂隙

Bakulin等（2000a, b, c）結(jié)合線性滑動理論和Hudson（1981）理論給出了不同裂隙介質(zhì)彈性矩陣的計算方法. 各向同性介質(zhì)背景下的兩組正交直立裂隙的彈性系數(shù)矩陣可表示為（Bakulinetal, 2000b）：

（1）

其中，

（2）

（3）

l1=1-ΔN1，l2=1-rΔN1，l3=1-r2ΔN1，l4=4r2g2ΔN1ΔN2，

m1=1-ΔN2，m2=1-rΔN2，m3=1-r2ΔN2，

（4）

式中：λ和μ分別為裂隙所在背景介質(zhì)的拉梅參數(shù)，ΔN和ΔT分別為法向柔度和切向柔度， 其大小與裂隙充填物有關(guān)， 具體表達(dá)式詳見Hudson（1981）.

2 三維正交方位各向異性介質(zhì)彈性波方程及高階交錯網(wǎng)格有限差分正演模擬

2.1 三維正交介質(zhì)彈性波一階速度-應(yīng)力方程

根據(jù)本構(gòu)方程， 幾何方程和運動微分方程推導(dǎo)出三維正交介質(zhì)一階速度-應(yīng)力方程為

（5）

（6）

式中,v為速度，τ為應(yīng)力，c為剛度矩陣元素，ρ為密度，F(xiàn)為震源項.

2.2 三維正交介質(zhì)彈性波交錯網(wǎng)格有限差分正演模擬

利用泰勒展開得到式（5）和式（6）的時間和空間差分形式：

時間1階差分：

圖2 變量定義的網(wǎng)格

（7）

空間2N階差分：

）.

（8）

交錯網(wǎng)格變量在網(wǎng)格上的定義形式與常規(guī)網(wǎng)格不同． 圖2給出了變量定義的網(wǎng)格， 7種網(wǎng)格點上定義不同的變量， 其具體定義形式列于表1.

2.3 正演模擬過程中的問題

表1 標(biāo)準(zhǔn)交錯網(wǎng)格中彈性波場分量和彈性參數(shù)的空間位置 Table 1 Position of elastic wavefield components and elastic parameters in standard staggered grid

3個問題， 即邊界條件處理， 頻散壓制和差分穩(wěn)定性條件．

鑒于只能在有限區(qū)域內(nèi)進行波動方程求解， 因此在區(qū)域邊界處會產(chǎn)生強烈的人為邊界反射. 該邊界反射對研究地震波的傳播沒有意義且會干擾對正常波場的認(rèn)識， 故相關(guān)人員近些年一直致力于消除邊界反射的研究， 產(chǎn)生了許多邊界條件， 本文采用PML邊界條件．

為避免直接進行傅里葉反變換需要在時間域作褶積， 首先對波場進行分解， 將不同波場分量分別沿x、y、z軸方向分裂， 并假定體力為0， 得到

（9）

這里以

（10）

為例， 推導(dǎo)其PML邊界方程. 根據(jù)復(fù)數(shù)坐標(biāo)變換規(guī)律

（11）

可得空間方向微分變換關(guān)系

（12）

式中dn（n）為不同方向的衰減系數(shù). 對式（10）進行傅里葉變換得

英格曼神甫的懇求得到了少佐的批準(zhǔn)。他的部隊在寒冷中靜默地多候了二十分鐘。英格曼給的理由是說得過去的：唱詩禮服很久沒被穿過，有的需要釘鈕扣，有的需要縫補、熨燙。士兵們站在圍墻外，一個挨一個，刺刀直指前方。多二十分鐘就多二十分鐘吧，好東西是值得等待的。日本人是最講究儀式的。一盤河豚上桌，都裝點成藝術(shù)品，何況美味的處女。

（13）

將上式代入式（12）并在時間和空間方向進行泰勒展開， 化簡可得

（14）

此則為式（10）的PML邊界表達(dá)式. 類似可推導(dǎo)出三維正交介質(zhì)PML表達(dá)式為

（15a）

（15b）

（15c）

（16a）

（16b）

（16c）

（16d）

（16e）

（16f）

圖3為震源位于一均勻立方正交介質(zhì)模型中心激發(fā)所得到的正演波場. 圖3a中縱波經(jīng)過邊界反射回來的能量強于橫波， 對波場分析造成嚴(yán)重干擾； 圖3b中縱波傳到邊界處被消除干凈， 相比圖3a中橫波能量較強． 可以看出PML邊界條件消除人工邊界的效果明顯． 這里注意， 由于u，v，w等3個波場分量上均有一個界面垂直于縱波傳播方向而導(dǎo)致截面上縱波很弱， 故這里每個波場分量只取兩個截面．

圖3 加最佳匹配層條件前（a）、 后（b）的三維正交各向異性介質(zhì)波場

3 模型試算

3.1 均勻模型正交介質(zhì)波場分析

利用基于線性滑動理論結(jié)合Hudson理論求取的對應(yīng)不同裂隙密度的彈性系數(shù)矩陣建立兩組正交直立裂隙的正交介質(zhì)均勻模型． 模型大小均為3.0 km×3.0 km×3.0 km； 混合密度ρ=2.625×103kg/m3； 拉梅參數(shù)分別為λ1=16.3 GPa,μ1=7.0 GPa，λ2=7.7 GPa,μ2=44.0 GPa； 背景介質(zhì)為砂巖和泥巖各占一半， 干燥裂隙， 裂隙縱橫比為0.01， 裂隙密度為0.04． 在模型中心采用等能量震源激發(fā)， 時間網(wǎng)格dt=1 ms， 空間網(wǎng)格dx=dy=dz=10 m， 分別在過震源的3個垂面截取三分量0.2 s時的波場快照進行對比．

圖4給出了兩組正交直立裂隙正交介質(zhì)的三維波場． 結(jié)合圖1可以看出， 由于裂隙的存在， 平行于裂隙的方向為各向同性面， 其波速大于垂直裂隙方向．xoy面內(nèi)波場傳播受到兩套相互正交的裂隙影響， 裂隙密度相等， 縱波波場近似成正方形；xoz面和yoz面分別受單套裂隙影響， 縱波波場近似成矩形， 垂直裂隙方向邊長小于平行裂隙方向． 所有截面內(nèi)橫波分裂嚴(yán)重， 說明當(dāng)裂隙密度達(dá)到0.1時介質(zhì)各向異性較強．

圖4 裂隙密度為0.1時的三維正交各向異性介質(zhì)波場

3.2 裂隙密度對波場的影響

通過改變裂隙密度建立4種裂隙模型， 裂隙密度分別給定為0.02， 0.04， 0.06和0.1． 圖5為正交介質(zhì)波場隨裂隙密度的變化情況． 可以看到， 在裂隙密度從左到右依次增大的過程中， 波場傳播速度減小， 縱波波場形狀由圓形逐步轉(zhuǎn)化為類似方形． 其中u分量的xoy面波場受到兩組裂隙的影響， 當(dāng)裂隙密度增大為0.1時其近似成正方形；u分量xoz面波場受到單組裂隙的影響， 當(dāng)裂隙密度增大為0.1時其近似成矩形． 從橫波波場變化可以看出， 隨著裂隙密度增大， 介質(zhì)各向異性增強， 其橫波分裂現(xiàn)象加重， 三叉區(qū)現(xiàn)象變得尤為明顯．

圖5 正交各向異性介質(zhì)波場隨裂隙密度ρF變化

圖6 層狀正交介質(zhì)模型

3.3 層狀模型分析

為分析不同裂隙密度對道集的影響， 我們建立多組三維正交介質(zhì)模型， 如圖6所示． 圖中黃色為各向同性層, 棕色為裂隙層， 大小為3.0 km×3.0 km×1.5 km， 時間網(wǎng)格dt=1 ms， 空間網(wǎng)格dx=dy=dz=10 m. 這些模型均為3層： 第一層和第三層為各向同性介質(zhì)， 給定縱、 橫波速度分別為4.0 km/s和2.5 km/s； 中間層為裂隙密度分別為0.02， 0.04， 0.06和0.1的正交介質(zhì).

首先分析裂隙密度為0.02時的共炮點道集和波場． 圖7為3個波場分量0°（沿x軸方向）測線的共炮點道集． 切除直達(dá)波后， 可以明顯看出兩個分界面的反射縱波和反射橫波．

炮點放置于模型頂面中心處， 過炮點截取u分量3個相互垂直面上0.2 s時的波場, 如圖8所示. 由于截取傳播時間的原因， 在波場上我們只能區(qū)分一個界面的透射縱橫波和反射縱橫波．

分別抽取與x軸方向成不同夾角的u分量道集, 繪制成圖9， 方位角分別為0°， 30°， 45°， 60°， 90°， 120°， 135°， 150°和180°等9個方位道集． 可以看出， 以90°為中心左右對稱（兩組裂隙密度相同）， 紅色橢圓標(biāo)出了不同方位角反射旅行時的變化， 其中45°和135°時道集上反射波旅行時最大， 這與其同時受到兩組裂隙的影響有關(guān)．

圖7 裂隙密度為0.02時層狀模型u分量正交介質(zhì)共炮點道集

圖8 裂隙密度為0.02時層狀模型正交介質(zhì)0.5 s時的波場

圖9 裂隙密度為0.02時層狀模型正交介質(zhì)u分量不同方位共炮點道集

不同裂隙密度下層狀模型u分量波場如圖10所示. 可以看出, 隨著裂隙密度的增大, 紅色橢圓中標(biāo)出的裂隙層底界面反射波旅行時的增大， 橫波分裂現(xiàn)象更突出， 說明裂隙對波場傳播和介質(zhì)各向異性性質(zhì)有較大影響．

圖10 不同裂隙密度ρF下層狀模型正交介質(zhì)u分量共炮點道集

4 討論與結(jié)論

本文利用線性滑動理論結(jié)合Hudson理論將兩組正交直立裂隙介質(zhì)等效為一種新的正交方位各向異性， 并用該方法進行巖石物理建模， 通過改變物性參數(shù)得到多組彈性參數(shù)． 在此基礎(chǔ)上， 推導(dǎo)了該正交介質(zhì)的彈性波動方程， 并采用有限差分法進行求解． 正演過程中使用PML邊界條件處理人為邊界反射． 最后通過分析正演波場和共炮點道集來模擬地震波在該介質(zhì)中的傳播過程， 分析裂隙密度等物性條件對波場和道集的影響以及方位各向異性特征．

本文結(jié)果表明， 裂隙密度的增大會使正交裂隙介質(zhì)的各向異性強度明顯增強， 從而導(dǎo)致垂直和平行裂隙方向上的縱波速度差別增大， 橫波變得更加復(fù)雜， 分裂嚴(yán)重． 通過對層狀模型正演模擬， 在共炮點道集和正演波場中可以清晰地分辨出直達(dá)波和反射波中的縱波和分裂后的橫波． 抽取不用方位測線記錄的共炮點道集， 通過比較反射波旅行時來描述介質(zhì)的方位特征， 結(jié)果表明不同方位上地震波傳播過程不盡相同， 這與介質(zhì)中發(fā)育的兩組直立正交裂隙不無關(guān)聯(lián)．

董良國， 馬在田， 曹景忠， 王華忠， 耿建華， 雷冰， 許世勇. 2000. 一階彈性波方程交錯網(wǎng)格高階差分解法[J]. 地球物理學(xué)報， 43（3）: 411--419.

Dong L G， Ma Z T， Cao J Z， Wang H Z， Geng J H， Lei B， Xu S Y. 2000. A staggered-grid high-order difference method of one-order elastic wave equation[J].ChineseJournalofGeophysics， 43（3）: 411--419 （in Chinese）.

何燕. 2008. 正交各向異性彈性波高階有限差分正演模擬研究[D]. 東營: 中國石油大學(xué)（華東）地球資源與信息學(xué)院: 4.

He Y. 2008.High-OrderFinite-DifferenceForwardModelingofElastic-WaveinOrthorhombicAnisotropicMedia[D]. Dongying: College of Georesources and Information， China University of Petroleum: 4 （in Chinese）.

金抒辛， 何樵登. 2005. 泥巖裂隙的正演模擬方法研究[J]. 石油物探， 44（2）: 119--123.

Jin S X， He Q D. 2005. The forward modeling of mudstone fractures[J].GeophysicalProspectingofPetroleum， 44（2）: 119--123 （in Chinese）.

孫成禹. 2007. 地震波理論與方法[M]． 東營： 中國石油大學(xué)出版社: 211.

Sun C Y.2007.TheoryandMethodsofSeismicWaves[M]. Dongying: China University of Petroleum Press: 211 （in Chinese）.

吳國忱. 2006. 各向異性介質(zhì)地震波傳播與成像[M]． 東營： 中國石油大學(xué)出版社: 1.

Wu G C. 2006.SeismicWavesPropagationandImaginginAnisotropicMedium[M]. Dongying: China University of Petroleum Press: 1 （in Chinese）.

吳國忱， 秦海旭. 2014. 裂縫介質(zhì)旋轉(zhuǎn)交錯網(wǎng)格正演模擬[J]. 地震學(xué)報， 36（6）: 1075--1088.

Wu G C， Qin H X. 2014. Rotated staggering grid forward modeling in fractured medium[J].ActaSeismologicaSinica， 36（6）: 1075--1088 （in Chinese）.

張文生， 宋海濱. 2001. 三維正交各向異性介質(zhì)三分量高精度有限差分正演模擬[J]. 石油地球物理勘探， 36（4）: 422--432

Zhang W S， Song H B. 2001. Finite-difference forward modeling of three component records with high precision in 3-D orthorhombic anisotropic media[J].OilGeophysicalProspecting， 36（4）: 422--432 （in Chinese）.

Backus G E. 1962. Long-wave elastic anisotropy produced by horizontal layering[J].JGeophysRes， 67（6）: 4427--4440.

Bakulin A， Grechka V， Tsvankin L. 2000a. Estimation of fracture parameters from reflection seismic data， part 1: HTI model due to a single fracture set[J].Geophysics， 65（6）: 1788--1802.

Bakulin A， Grechka V， Tsvankin L. 2000b. Estimation of fracture parameters from reflection seismic data， part 2: Fracture models with orthorhombic symmetry[J].Geophysics， 65（6）: 1803--1817.

Bakulin A， Grechka V， Tsvankin L. 2000c. Estimation of fracture parameters from reflection seismic data， part 3: Fracture models with monoclinic symmetry[J].Geophysics， 65（6）: 1818--1830.

Berenger J P. 1994. A perfectly matched layer for absorption of electromagnetic waves[J].JComputPhys， 114（2）: 185--200.

Boris J P， Book D L. 1973. Flux-corrected transport. I: SHASTA， a fluid transport algorithm that works[J].JComputPhys， 11（1）: 38--69.

Coates R T， Schoenberg M. 1995. Finite-difference modeling of faults and fractures[J].Geophysics， 60（5）: 1514--1526.

Collino F， Tsogka C. 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics， 66（1）: 294--307.

Crampin S. 1978. Seismic-wave propagation through a cracked solid: Polarization as a possible dilatancy diagnostic[J].GeophysJInt， 53（3）: 467--496.

Crampin S， McGonigle R， Bamford D. 1980. Estimating crack parameters from observations of P-wave velocity anisotropy[J].Geophysics， 45（3）: 345--360.

Crampin S. 1981. A review of wave motion in anisotropic and cracked elastic-media[J].WaveMotion， 3（4）: 343--391.

Crampin S. 1987. Geological and industrial implications of extensive-dilatancy anisotropy[J].Nature， 328（6130）: 491--496.

Hudson J A. 1981. Wave speeds and attenuation of elastic waves in material containing cracks[J].GeophysJInt， 64（1）: 133--150.

Schoenberg M. 1980. Elastic wave behavior across linear slip interfaces[J].JAcoustSocAm， 68（5）: 1516--1521.

Schoenberg M， Muir F. 1989. A calculus for finely layered anisotropic media[J].Geophysics， 54（5）: 581--589.

Schoenberg M， Sayers C. 1995. Seismic anisotropy of fractured rock[J].Geophysics， 60（1）: 204--221.

Schoenberg M， Hilbig K. 1997. Orthorhombic media: Modeling elastic wave behavior in a vertically fractured earth[J].Geophysics， 62（6）: 1954--1974.

Virieux J. 1987. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[J].Geophy-sics， 51（4）: 889--901.

Rock physics modeling and elastic wave forward modeling methods in a three-dimensional azimuthally anisotropic medium of orthorhombic symmetry

（SchoolofGeosciences,ChinaUniversityofPetroleum，ShandongQingdao266580,China）

In this paper， by rock physical equivalent theory and linear slip theory， two sets of orthorhombic fractures are equivalent to a new three-dimensional azimuthally anisotropic medium of orthorhombic symmetry for rock physical modeling. And then elastic wave equations are solved to simulate the propagation process in this medium by high-order staggered-grid finite-difference method. During the process of modeling， with the variation of physical parameters， we analyze CSP （common shot points） sets and the wave field characteristics on the conditions of different fracture densities， as well as azimuthal characteristics of orthorhombic anisotropy. The research results show that the anisotropic strengths enhance with fracture density increasing， which is reflected in the CSP sets and forwarding wavefield.

azimuthal anisotropy of orthorhombic symmetry; fracture medium; Hudson theory; linear slip theory; staggered-grid finite-difference

10.11939/jass.2015.04.013.

國家973項目（2013CB228604）資助.

2014-11-19收到初稿， 2015-04-27決定采用修改稿.

e-mail: 15610031211@163.com

4.013

P315.3+1

A

李雨生, 吳國忱. 2015. 一種三維正交方位各向異性介質(zhì)巖石物理建模及彈性波正演模擬方法. 地震學(xué)報, 37（4）: 678--689.

Li Y S, Wu G C. 2015. Rock physics modeling and elastic wave forward modeling methods in a three-dimensional azimuthally anisotropic medium of orthorhombic symmetry.ActaSeismologicaSinica, 37（4）: 678--689.

doi:10.11939/jass.2015.04.013.