嚴(yán)紅勇，劉 洋＊

1中國(guó)石油大學(xué)（北京）油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249 2中國(guó)石油大學(xué)（北京）CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249

黏彈TTI介質(zhì)中旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬

嚴(yán)紅勇1，2，劉 洋1，2＊

1中國(guó)石油大學(xué)（北京）油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249 2中國(guó)石油大學(xué)（北京）CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249

以Carcione黏彈各向異性理論為基礎(chǔ)，給出了適用于黏彈性具有任意傾斜對(duì)稱軸橫向各向同性介質(zhì)（黏彈TTI介質(zhì)）的二維三分量一階速度－應(yīng)力方程，采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格任意偶數(shù)階精度有限差分格式求解該方程，并推導(dǎo)出了二維黏彈TTI介質(zhì)完全匹配層（PML）吸收邊界條件公式和相應(yīng)的旋轉(zhuǎn)交錯(cuò)網(wǎng)格任意偶數(shù)階精度有限差分格式，實(shí)現(xiàn)了該類介質(zhì)的地震波場(chǎng)數(shù)值模擬．?dāng)?shù)值模擬結(jié)果表明：該方法模擬精度高，邊界吸收效果好，可以得到高精度的波場(chǎng)快照和合成記錄；并且波場(chǎng)快照和合成記錄能較好地反映地下介質(zhì)的各向異性特征和黏彈性特征．

黏彈性，各向異性，旋轉(zhuǎn)交錯(cuò)網(wǎng)格，有限差分，完全匹配層

1 引 言

地震波在地下介質(zhì)中傳播時(shí)，不僅受介質(zhì)的各向異性影響，還受地層介質(zhì)內(nèi)在的黏滯性影響［1］．因此，為了準(zhǔn)確地描述地震波在地下介質(zhì)中的傳播特征和更精確地指導(dǎo)地震數(shù)據(jù)采集、處理和解釋，應(yīng)選取能夠同時(shí)模擬地層各向異性特征和黏滯性特征的黏彈各向異性介質(zhì)模型來(lái)進(jìn)行波場(chǎng)數(shù)值模擬與分析．

目前，黏彈各向異性介質(zhì)的地震波場(chǎng)數(shù)值模擬研究已經(jīng)取得了一些進(jìn)展．Carcione［1－3］基于廣義線性固體模型引入新的黏彈性各向異性本構(gòu)關(guān)系，建立了黏彈各向異性介質(zhì)中隨時(shí)間而變化的松弛矩陣，而且推導(dǎo)出了黏彈VTI介質(zhì)的波動(dòng)方程，并模擬了縱波、橫波的各向異性和衰減特征．張忠杰等［4－5］用近似解析分析法模擬了黏彈EDA介質(zhì)中的波場(chǎng)．杜啟振等［6－8］在Carcione的基礎(chǔ)上給出了以各向異性主軸方位為參數(shù)的各向異性黏彈性介質(zhì)波動(dòng)方程，并先后采用有限元法、偽譜法和有限差分法對(duì)方位各向異性黏彈性介質(zhì)中的波場(chǎng)進(jìn)行了模擬與分析．郭智奇等［9］利用偽譜法模擬了黏彈VTI介質(zhì)中的波場(chǎng)．譚未一等［10］采用規(guī)則網(wǎng)格有限差分法模擬了黏彈VTI介質(zhì)中的波場(chǎng)．李軍等［11］利用交錯(cuò)網(wǎng)格高階有限差分對(duì)黏彈性VTI介質(zhì)中的波場(chǎng)進(jìn)行了數(shù)值模擬．然而，這些研究對(duì)黏彈各向異性介質(zhì)的波場(chǎng)模擬主要針對(duì)兩種特殊情況：黏彈VTI介質(zhì)或黏彈HTI介質(zhì)．而具有任意傾斜對(duì)稱軸橫向各向同性（TTI）介質(zhì)在實(shí)際地層各向異性的表現(xiàn)中更具有廣泛性［12－14］．因此，如果能給出黏彈TTI介質(zhì)波動(dòng)方程，并針對(duì)其進(jìn)行波場(chǎng)數(shù)值模擬與分析，則更具有一般性．

黏彈TTI介質(zhì)可以看成是黏彈VTI介質(zhì)經(jīng)過(guò)旋轉(zhuǎn)一定的角度后得到的，其對(duì)稱軸與在地面建立的觀測(cè)坐標(biāo)系存有一定的夾角，如果利用普通交錯(cuò)網(wǎng)格法求解二維三分量的黏彈TTI介質(zhì)波動(dòng)方程，為計(jì)算應(yīng)力而求質(zhì)點(diǎn)振動(dòng)速度的空間偏導(dǎo)數(shù)時(shí)，有些偏導(dǎo)不能直接求解，需要進(jìn)行漂移或借助輔助網(wǎng)格．而且在普通交錯(cuò)網(wǎng)格中，要求彈性模量定義在所有的半網(wǎng)格點(diǎn)和整網(wǎng)格點(diǎn)上，而實(shí)際中通常只定義在半網(wǎng)格點(diǎn)上，這就需要對(duì)彈性模量進(jìn)行插值［15］．因此，當(dāng)彈性模量變化很大時(shí)，普通交錯(cuò)網(wǎng)格法就容易出現(xiàn)計(jì)算不穩(wěn)定的問(wèn)題［16］．Saengere等［17－18］首先提出了旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分算法，并指出其可以克服普通交錯(cuò)網(wǎng)格差分的上述缺陷．而后，一些學(xué)者利用旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分法分別對(duì)孔隙介質(zhì)、各向異性介質(zhì)等的波場(chǎng)進(jìn)行了數(shù)值模擬［19－21］．

旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分的所有物理量只定義在兩個(gè)不同的位置．其中，應(yīng)力（或質(zhì)點(diǎn)振動(dòng)速度）定義在離散網(wǎng)格單元的中心，質(zhì)點(diǎn)振動(dòng)速度（或應(yīng)力）定義在離散網(wǎng)格單元的頂點(diǎn)，故可將所有的彈性模量定義在相同的位置，且可與應(yīng)力所在的位置相同［17，19］．因此，采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分法，在模擬非均勻介質(zhì)時(shí)，彈性模量不需要進(jìn)行平均；在模擬TTI介質(zhì)時(shí)，計(jì)算應(yīng)力的所有物理量的空間偏導(dǎo)數(shù)可以直接求解，且彈性模量也不再需要進(jìn)行插值處理［19］．所以旋轉(zhuǎn)交錯(cuò)網(wǎng)格差分法更適合黏彈TTI介質(zhì)的波場(chǎng)模擬，并且相對(duì)于普通交錯(cuò)網(wǎng)格差分法，旋轉(zhuǎn)交錯(cuò)網(wǎng)格差分法可以得到更加穩(wěn)定、更可信的解．

本文在Carcione［1－2］提出的黏彈各向異性本構(gòu)關(guān)系基礎(chǔ)上，給出黏彈TTI介質(zhì)的二維三分量一階速度－應(yīng)力方程，然后采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格任意偶數(shù)階精度有限差分格式求解該方程，并將PML引入到黏彈TTI介質(zhì)波動(dòng)方程的旋轉(zhuǎn)交錯(cuò)網(wǎng)格差分?jǐn)?shù)值模擬中，實(shí)現(xiàn)了該類介質(zhì)波場(chǎng)的旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬，并且根據(jù)模擬結(jié)果具體分析了地震波在傳播過(guò)程中所表現(xiàn)出的各向異性和黏彈性特征．

2 黏彈TTI介質(zhì)波動(dòng)方程

2．1 觀測(cè)坐標(biāo)系內(nèi)的本構(gòu)關(guān)系

Carcione的黏彈各向異性理論認(rèn)為介質(zhì)對(duì)地震波的吸收衰減作用是由大量耗散機(jī)制引起的，且各種耗散機(jī)制可以用黏彈性的本構(gòu)關(guān)系進(jìn)行描述，采用其建立的應(yīng)力應(yīng)變關(guān)系［1］，即

其中，xv（0）（v＝1，2）表示0時(shí)刻的弛豫模量，當(dāng)v＝1時(shí)，對(duì)應(yīng)于縱波，當(dāng)v＝2時(shí)，對(duì)應(yīng)于橫波；TT＝（T1，T2，T3，T4，T5，T6）＝（σxx，σyy，σzz，σyz，σzx，σxy）表示應(yīng)力向量；ST＝（S1，S2，S3，S4，S5，S6）＝（εxx，εyy，εzz，εyz，εzx，εxy）表示應(yīng)變向量；K＝1，2，…，6；v＝1，2）分別表示非弛豫空間函數(shù)和弛豫空間函數(shù)，并且其與彈性參數(shù)有關(guān)；Lv（v＝1，2）表示縱波或橫波弛豫機(jī)制的總數(shù)是記憶變量的時(shí)間導(dǎo)數(shù)；φvl（0）（l＝1，2，…，6；v＝1，2）表示第l個(gè)弛豫機(jī)制的響應(yīng)函數(shù)；表示第l個(gè)機(jī)制的應(yīng)力弛豫時(shí)間．上述變量及函數(shù)的具體表達(dá)式可以參見(jiàn)參考文獻(xiàn)［1］．

根據(jù)黏彈各向異性介質(zhì)的本構(gòu)關(guān)系，坐標(biāo)旋轉(zhuǎn)以后，可以得到觀測(cè)坐標(biāo)系下以任意方位角和傾角為參數(shù)的應(yīng)力應(yīng)變關(guān)系方程［1］，即

式中，M為旋轉(zhuǎn)變換矩陣［12］，且

其中，θ和φ分別為在觀測(cè)坐標(biāo)系下的方位角和傾角．將式（3）變形可得

式中矩陣定義如下

2．2 一階速度－應(yīng)力方程

由黏彈性各向異性介質(zhì)的本構(gòu)關(guān)系，給出黏彈TTI介質(zhì)的二維三分量一階速度－應(yīng)力方程［1，2］．

運(yùn)動(dòng)方程：

其中，ρ為介質(zhì)的密度，vx、vy、vz分別是質(zhì)點(diǎn)速度x、y、z的分量，σxx、σzz為正應(yīng)力，σxz、σxy、σyz為切應(yīng)力．

應(yīng)力－速度關(guān)系：

記憶變量方程：

由式（7）、（8）、（9）構(gòu)成黏彈TTI介質(zhì)中波動(dòng)方程組，以此為基礎(chǔ)對(duì)黏彈TTI介質(zhì)做二維三分量波場(chǎng)數(shù)值模擬．

3 旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分格式

圖1解釋了旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分算法的二維空間交錯(cuò)方法．在旋轉(zhuǎn)交錯(cuò)網(wǎng)格中，所有的速度分量都被定義在同一位置，而所有的應(yīng)力分量位于另外的同一位置（黏彈TTI介質(zhì)波動(dòng)方程中的記憶變量的導(dǎo)數(shù)都定義在與應(yīng)力分量相同的位置），因此，也可以將介質(zhì)密度和彈性模量定義在相同的位置［19］，本文中將它們定義在與應(yīng)力分量相同的位置上．

圖1 旋轉(zhuǎn)交錯(cuò)網(wǎng)格的二維網(wǎng)格單元示意圖［17］Fig．1 Stencil of discretization scheme in rotated staggered grid scheme in 2D

旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分法是先沿網(wǎng)格對(duì)角線方向計(jì)算相關(guān)物理量的空間導(dǎo)數(shù)，然后再將這些對(duì)角線方向的計(jì)算結(jié)果進(jìn)行線性疊加來(lái)得到沿坐標(biāo)軸方向的空間導(dǎo)數(shù)［19］．在二維空間的情形下，旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分沿坐標(biāo)軸方向的差分算子的計(jì)算表達(dá)式為［17］

式中，N為差分階數(shù)的一半；m＝1，2，…，N；am為差分系數(shù)，其值詳細(xì)推導(dǎo)與普通交錯(cuò)網(wǎng)格的差分系數(shù)的推導(dǎo)是一致的，詳見(jiàn)參考文獻(xiàn)［22］．

由式（12）—（15）可以推導(dǎo)出在x和z方向的偏導(dǎo)數(shù)，即

利用式（14）—（17）直接求解黏彈TTI介質(zhì)的波動(dòng)方程組，可以得到方程的旋轉(zhuǎn)交錯(cuò)網(wǎng)格任意偶數(shù)階精度差分格式．

結(jié)合各向異性介質(zhì)中彈性波方程差分算法的穩(wěn)定性條件和旋轉(zhuǎn)交錯(cuò)網(wǎng)格的特性，給出黏彈TTI介質(zhì)中二維一階速度－應(yīng)力方程的時(shí)間二階精度、空間2 N階精度的旋轉(zhuǎn)交錯(cuò)網(wǎng)格差分穩(wěn)定性條件近似表達(dá)式［23］，即

4 PML吸收邊界條件

Berenger于1994年首次提出了PML吸收邊界條件，并將其運(yùn)用來(lái)消除電磁波場(chǎng)數(shù)值模擬中的人為截?cái)噙吔缢a(chǎn)生的反射波［24］．隨后許多學(xué)者將此方法引入到了地震波場(chǎng)的數(shù)值模擬中，并取得了比較好的效果［25］．本文將PML引入到用旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分求解的黏彈TTI介質(zhì)波動(dòng)方程中．在二維空間情形下，旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分法中實(shí)現(xiàn)PML的思路，與普通交錯(cuò)網(wǎng)格有限差分中PML算法的基本思路是相同的．根據(jù)PML吸收邊界條件的基本原理，采用時(shí)間域變量分裂的方法，對(duì)一階速度－應(yīng)力黏彈TTI介質(zhì)波動(dòng)方程進(jìn)行波場(chǎng)變量分離，可以得到應(yīng)用于黏彈TTI介質(zhì)波動(dòng)方程的PML吸收邊界條件．在二維空間的條件下，對(duì)于式（7）、（8）、（9），每一個(gè)波場(chǎng)變量可分為兩部分，即

式中上標(biāo)x和z代表該項(xiàng)只與相應(yīng)的空間導(dǎo)數(shù)有關(guān)．同時(shí)，可以得到分裂后的PML吸收邊界中帶有衰減因子的波動(dòng)方程（下面只以σxx為例，逐步推導(dǎo)給出其吸收層邊界方程的具體形式），即

其中，函數(shù)d（x）和d（z）是與x和z方向有關(guān)的衰減因子［25］．

在數(shù)值模擬過(guò)程中PML吸收層被分成了不同的區(qū)域，如圖2所示，不同區(qū)域的衰減因子的取值是不同的．在b1和b3代表的吸收層區(qū)域，d（x）＝0，d（z）＞0；在b2和b4代表的吸收層區(qū)域，d（x）＞0，d（z）＝0；在a1、a2、a3和a4代表的PML吸收層區(qū)域，d（x）＞0，d（z）＞0．

圖2 完全匹配層吸收邊界示意圖［25］Fig．2 Stencil of perfectly matched layer absorbing boundary conditions

限于篇幅，僅以式（20）中σxx分量和式（22）為例進(jìn)行差分離散，給出其二階時(shí)間差分精度、高階空間差分的旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分格式，即D和D分別表示?／?和?／?在和方向的任意偶數(shù)階偏導(dǎo)數(shù)式，故將式（14）和式（15）代入式（25）中，其余分裂變量采用同樣的方式推導(dǎo)，即可求得PML邊界條件的任意偶數(shù)階旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分格式．

5 數(shù)值模擬與分析

5．1 PML吸收邊界條件效果驗(yàn)證

為了驗(yàn)證PML吸收邊界的效果，設(shè)計(jì)了一個(gè)大小為2000m×2000m、方位角為45°、對(duì)稱軸傾角為60°的均勻黏彈TTI介質(zhì)模型，空間采樣間隔Δx＝Δz＝5m，時(shí)間步長(zhǎng)為1ms．模型介質(zhì)在自然坐標(biāo)系下的彈性參數(shù)如表1，模型的縱波Q值為80，橫波Q值為60．震源子波采用Ricker子波，震源位置（1000m，1000m），主頻為20Hz，震源為x方向的橫波源．

在數(shù)值模擬過(guò)程中采用二十階空間差分精度、二階時(shí)間差分精度的旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分格式對(duì)上述模型做二維三分量波場(chǎng)數(shù)值模擬．圖3（a，b，c）分別是考慮吸收邊界條件時(shí)模擬得到的x、y、z三個(gè)分量的波場(chǎng)快照，在每個(gè)分量，四個(gè)波場(chǎng)快照從左到右對(duì)應(yīng)的時(shí)間依次為600ms、840ms、1200ms和1800ms．從圖3可以看出，x、y、z三個(gè)分量入射到截?cái)噙吔缣幍母黝愅庑胁ǘ急煌耆?，這表明本文在用旋轉(zhuǎn)交錯(cuò)網(wǎng)格差分法模擬黏彈TTI介質(zhì)中波場(chǎng)時(shí)所采用的PML邊界條件對(duì)各類外行波都具有良好的吸收性能．

表1 均勻黏彈性TI介質(zhì)在自然坐標(biāo)系下的彈性參數(shù)Table 1 The elastic parameters of the homogeneous viscoelastic TTI media in the normal axis

5．2 各向異性特征分析

設(shè)計(jì)了一個(gè)大小為2000m×2000m均勻黏彈TI介質(zhì)模型，空間采樣間隔Δx＝Δz＝5m，時(shí)間步長(zhǎng)為1ms．模型的縱波Q值為80，橫波Q值為60，模型介質(zhì)在自然坐標(biāo)系下的彈性參數(shù)如表1．震源子波采用Ricker子波，震源位置（1000m，1000m），主頻為20Hz，震源為x方向的橫波源．?dāng)?shù)值模擬差分精度為：十八階空間差分精度、二階時(shí)間差分精度．圖4是黏彈TTI介質(zhì)在方位角和對(duì)稱軸傾角都為0時(shí)模擬得到的波場(chǎng)快照和單炮記錄．圖5是黏彈TTI介質(zhì)在方位角為45°、對(duì)稱軸傾角為60°時(shí)模擬得到的波場(chǎng)快照和單炮記錄．

從圖4中的各個(gè)分量可以看出：當(dāng)黏彈TTI介質(zhì)的對(duì)稱軸傾角和方位角都等于0時(shí)，均勻黏彈TI介中的地震波記錄包括擬縱波qP波波形和擬橫波qS波波形；波場(chǎng)快照中qS波波前面的形狀與單炮傳播記錄中qS波波形初至所形成的三叉區(qū)很相似，qS波波形及同相軸三叉區(qū)的出現(xiàn)導(dǎo)致波場(chǎng)變得更加復(fù)雜化；且單炮記錄和波場(chǎng)快照中y分量均為0．

圖3 黏彈TTI介質(zhì)中加PML吸收邊界的波場(chǎng)快照Fig．3 Snapshots of wavefield in viscoelastic TTI media with PML absorbing boundary conditions

圖4 黏彈各向異性介質(zhì)中地震波波場(chǎng)快照（左圖）和傳播記錄（右圖）介質(zhì)方位角和傾角都為0，（a）x分量；（b）z分量．Fig．4 Snapshots（left）and records（right）of seismic wave caused in viscoelastic and anisotropic media Media azimuth and dip angle are both zero degree（a）x－component，（b）z－component．

圖5 黏彈TTI介質(zhì)中地震波波場(chǎng)快照（左圖）和單炮傳播記錄（右圖）介質(zhì)方位角為45°、對(duì)稱軸傾角為60°；（a）x分量；（b）y分量；（c）z分量．Fig．5 Snapshots（left）and records（right）of seismic wave caused in viscoelastic and anisotropic media Media azimuth is 45degree and dip angle of symmetry axis is 60degree（a）x－component；（b）y－component；（c）z－component．

從圖5中的各個(gè)分量可以看到：當(dāng)黏彈TTI介質(zhì)的對(duì)稱軸傾角不為0時(shí)，均勻黏彈TTI介質(zhì)中地震波記錄的三個(gè)分量中都有qP波波形、快橫波qS1波形和慢橫波qS2波形，即存在橫波分裂現(xiàn)象；從單炮傳播記錄中也可以看到，各類波的同相軸變得不對(duì)稱，且不一定是雙曲線形態(tài)，并且隨著偏移距的增大，在其對(duì)稱軸上傾的方向上，qP波形的初至?xí)r間要逐漸小于下傾向上qP波形的初至?xí)r間，而兩類qS波的情況則相反［27］．

為了簡(jiǎn)明地分析黏彈TTI介質(zhì)中地震波的反射特征，設(shè)計(jì)了一個(gè)大小為5000m×3000m兩層介質(zhì)模型．模型上層為黏彈TTI介質(zhì)，其方位角為45°、對(duì)稱軸傾角為60°；下層為黏彈各向同性介質(zhì)，空間采樣間隔Δx＝Δz＝5m，時(shí)間步長(zhǎng)為1ms．模型介質(zhì)參數(shù)如表2．震源子波采用Ricker子波，震源位置（2500m，1200m），主頻為25Hz，震源為脹縮源．?dāng)?shù)值模擬差分精度為二十階空間差分精度、二階時(shí)間差分精度．

表2 兩層模型介質(zhì)在自然坐標(biāo)系下的介質(zhì)參數(shù)Table 2 The media′s parameters of the double layers in the normal axis

圖6是采用表2中介質(zhì)參數(shù)模擬得到的波場(chǎng)快照和單炮記錄．從圖6中可以看到，當(dāng)上層黏彈TI介質(zhì)的方位角和對(duì)稱軸傾角都不為0時(shí)：（1）經(jīng)過(guò)黏彈TTI介質(zhì)與黏彈各向同性介質(zhì)分界面的反射和透射后，其中地震反射波包括：擬縱波qP波、轉(zhuǎn)換的快橫波qPS1波、轉(zhuǎn)換的慢橫波qPS2波、擬快橫波qS1波、擬慢橫波qS2波、轉(zhuǎn)換的擬縱波qSP波，且可看到明顯的橫波分裂現(xiàn)象；而在下層黏彈各向同性介質(zhì)中的透射波只存在P波、S波、轉(zhuǎn)換波SP波和轉(zhuǎn)換波PS波；（2）在黏彈TTI介質(zhì)中，qS波的波形初至?xí)纬扇鎱^(qū)［26］，這種特殊現(xiàn)象使地震波場(chǎng)變得更為復(fù)雜化，并特別容易引起誤解；（3）黏彈TTI介質(zhì)中，單炮記錄上各類反射波的同相軸不再是雙曲線型，且qS反射波的同相軸偏離最為明顯．

圖6 利用表2中介質(zhì)參數(shù)模擬得到的地震波波場(chǎng)快照（左圖）和單炮記錄（右圖）介質(zhì)方位角為45°、對(duì)稱軸傾角為60°；（a）為x分量，（b）為y分量，（c）為z分量；圖中各種波前的字母：D、R和T分別表示直達(dá)波、反射波和透射波．Fig．6 Snapshots（left）and records（right）of seismic wave simulated by media parameters according to table 2 Media azimuth is 45degree and dip angle of symmetry axis is 60degree（a）x－component；（b）y－component；（c）z－component；D，Rand Tstand for directed wave，reflected wave and transmitted wave，respectively．

5．3 黏彈性特征分析

保持介質(zhì)的彈性參數(shù)不變，但改變其縱波和橫波的Q值，設(shè)計(jì)四種黏彈TTI介質(zhì)模型以及一個(gè)TTI彈性介質(zhì)模型（Q值無(wú)窮大），具體的介質(zhì)參數(shù)如表3．設(shè)計(jì)模型的大小為2000m×2000m、方位角為45°、對(duì)稱軸傾角為60°，空間采樣間隔Δx＝Δz＝5m，時(shí)間步長(zhǎng)為1ms．震源子波采用Ricker子波，震源位置（1000m，1000m），主頻為20Hz，震源為x方向的橫波源．利用二十階空間差分精度、二階時(shí)間差分精度對(duì)四個(gè)模型分別做正演模擬，?。?500m，500m）處質(zhì)點(diǎn)的x、y、z三個(gè)分量振動(dòng)記錄進(jìn)行分析．圖7、8、9分別是對(duì)表3中介質(zhì)1、2、3、4模擬得到的質(zhì)點(diǎn)（1500m，500m）處的x、y、z分量記錄．圖10是對(duì)表3中介質(zhì)1、5、4分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的x分量記錄對(duì)比圖．而圖11是對(duì)表3中介質(zhì)1、5、3分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的x分量記錄對(duì)比圖．

表3 黏彈TTI介質(zhì)在自然坐標(biāo)系下的介質(zhì)參數(shù)Table 3 The media′s parameters of the viscoelastic TTI media in the normal axis

從圖7、8、9中可以看到，無(wú)論是x、y、z任一分量，由于介質(zhì)的黏滯性而使質(zhì)點(diǎn)的振動(dòng)振幅都明顯衰減，并且qS波比qP波的衰減幅度要大，而且品質(zhì)因子越小，振幅衰減的越多，接收到的子波的峰值時(shí)間越往后延遲，并且波形變寬、畸變也越嚴(yán)重．從圖10中可以看出，縱波品質(zhì)因子QP對(duì)應(yīng)于膨脹滯彈性的形變，且QP越小qP波振幅衰減越嚴(yán)重．從圖11中可以看出，橫波品質(zhì)因子QS對(duì)應(yīng)于剪切滯彈性的形變，且QS越小qS波振幅衰減越嚴(yán)重．

6 結(jié) 論

圖7 對(duì)表3中介質(zhì)1、2、3、4分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的x分量記錄Fig．7 Records of x－component in point（1500m，500m）simulated by media 1，2，3and 4according to table 3

采用旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分法模擬黏彈TTI介質(zhì)中地震波場(chǎng)的傳播過(guò)程，模擬精度較高，且易于實(shí)現(xiàn)．同時(shí)，將完全匹配層吸收邊界條件應(yīng)用到黏彈TTI介質(zhì)旋轉(zhuǎn)交錯(cuò)網(wǎng)格高階有限差分模擬中，能有效地消除人工邊界所產(chǎn)生的反射．

對(duì)模擬得到的波場(chǎng)進(jìn)行分析表明：

（1）黏彈TTI介質(zhì)的各向異性主要影響波前面的形態(tài)，并且波前面形態(tài)復(fù)雜；qS波波前面和同相軸的三分叉現(xiàn)象普遍，且其同相軸一般不再是雙曲線型；當(dāng)TI介質(zhì)傾斜時(shí)，三個(gè)分量上均能夠觀測(cè)到橫波分裂現(xiàn)象，而且各類波形的同相軸變得不對(duì)稱；并且TI介質(zhì)傾斜時(shí)，地震波通過(guò)界面反射后，其反射特征變得更為復(fù)雜，能夠觀測(cè)到反射橫波的分裂現(xiàn)象．

圖8 對(duì)表3中介質(zhì)1、2、3、4分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的y分量記錄Fig．8 Records of y－component in point（1500m，500m）simulated by media 1，2，3and 4according to table 3

圖9 對(duì)表3中介質(zhì)1、2、3、4分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的z分量記錄Fig．9 Records of z－component in point（1500m，500m）simulated by media 1，2，3and 4according to table 3

圖10 對(duì)表3中介質(zhì)1、5、4分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的x分量記錄Fig．10 Records of x－component in point（1500m，500m）caused by media 1，5and 4according to table 3

圖11 對(duì)表3中介質(zhì)1、5、3分別模擬得到的質(zhì)點(diǎn)（1500m，500m）處的x分量記錄Fig．11 Records of x－component in point（1500m，500m）caused by media 1，5and 3according to table 3

（2）黏彈TTI介質(zhì)的黏彈性主要影響地震波的振幅衰減和波形畸變，對(duì)應(yīng)于膨脹滯彈性形變的縱波品質(zhì)因子主要影響qP波的振幅衰減和波形畸變，對(duì)應(yīng)于剪切滯彈性形變的橫波品質(zhì)因子主要影響qS波的振幅衰減和波形畸變；而且，當(dāng)品質(zhì)因子越小時(shí)，地震波振幅衰減的越多，接收到的子波的峰值時(shí)間越往后延遲，且其波形變寬、畸變也越嚴(yán)重．

（References）

［1］ Carcione J M．Wave propagation in anisotropic linear viscoelastic media：Theory and simulated wavefields．Geophysical Journal International，1990，101（3）：739－750．

［2］ Carcione J M．Constitutive model and wave equations for linear，viscoelastic，anisotropic media．Geophysics，1995，60（2）：537－548．

［3］ Carcione J M．Staggered mesh for the anisotropic and viscoelastic wave equation．Geophysics，1999，64（6）：1863－1866．

［4］ Zhang Z J，Wang G J，Harris J M．Multi－component wavefield simulation in viscous extensively dilatancy anisotropic media．Physics of the Earth and Planetary Interiors，1999，114（1－2）：25－38．

［5］ 張忠杰，滕吉文，賀振華．EDA介質(zhì)中地震波速度、衰減與品質(zhì)因子方位異性研究．中國(guó)科學(xué)（E輯），1999，29（6）：569－574．Zhang Z J，Teng J W，He Z H．Azimuthal anisotropy of seismic velocity，attenuation and Q value in viscous EDA media．Science in China（Series E）（in Chinese），1999，29（6）：569－574

［6］ 杜啟振，楊慧珠．方位各向異性黏彈性介質(zhì)波場(chǎng)有限元模擬．物理學(xué)報(bào)，2003，52（8）：2010－2014．Du Q Z，Yang H Z．Finite－element methods for viscoelastic and azimuthally anisotropic media．Acta Physica Sinica（in Chinese），2003，52（8）：2010－2014．

［7］ 杜啟振．各向異性黏彈性介質(zhì)偽譜法波場(chǎng)模擬．物理學(xué)報(bào)，2004，53（12）：4428－4434．Du Q Z．Wavefield forward modeling with the pseudospectral method in viscoelastic and azimuthally anisotropic media．Acta Physica Sinica（in Chinese），2004，53（12）：4428－4434．

［8］ 杜啟振，王延光，付水華．方位各向異性黏彈性介質(zhì)波場(chǎng)數(shù)值模擬．地球物理學(xué)進(jìn)展，2006，21（2）：502－504．Du Q Z，Wang Y G，F(xiàn)u S H．Wavefield forward modeling with the pseudo－spectral method in viscoelastic and azimuthally anisotropic media．Progress in Geophysics（in Chinese），2006，21（2）：502－504．

［9］ 郭智奇，劉財(cái)，楊寶俊等．黏彈各向異性介質(zhì)中地震波場(chǎng)模擬與特征．地球物理學(xué)進(jìn)展，2007，22（3）：804－810．Guo Z Q，Liu C，Yang B J．Seismic wave fields modeling and feature in viscoelastic anisotropic media．Progress in Geophysics（in Chinese），2007，22（3）：804－810．

［10］ 譚未一，趙兵，張中杰．黏彈性VTI介質(zhì)地震波模擬特征分析．地球物理學(xué)進(jìn)展，2007，22（4）：1305－1311．Tan W Y，Zhao B，Zhang Z J．The analysis of the wave field characteristics in viscoelastic VTI medium．Progress in Geophysics（in Chinese），2007，22（4）：1305－1311．

［11］ 李軍，蘇云，李錄明．各向異性黏彈性介質(zhì)波場(chǎng)數(shù)值模擬．大慶石油學(xué)院學(xué)報(bào)，2009，33（6）：38－42．Li J，Su Y，Li L M．Wavefield numeric simulation in anisotropic viscoelastic media．Journal of Daqing Petroleum Institute（in Chinese），2009，33（6）：38－42．

［12］ 劉洋，李承楚，牟永光．具有傾斜對(duì)稱軸的橫向各向同性介質(zhì)中的彈性波．石油地球物理勘探，1998，33（2）：161－169．Liu Y，Li C C，Mou Y G．Elastic wave in transverse isotropic media with dipping symmetric axis．Oil Geophysical Prospecting（in Chinese），1998，33（2）：161－169．

［13］ Qin Y L，Zhang Z J，Li S L．CDP mapping in tilted transversely isotropic（TTI）media．Part I：Method and effectiveness．Geophysical Prospecting，2003，51（4）：315－324．

［14］ Qin Y L，Zhang Z J，Xu S H．CDP mapping in tilted transversely isotropic（TTI）media．Part II：Velocity analysis by combining CDP mapping with a genetic algorithm．Geophysical Prospecting，2003，51（4）：325－332．

［15］ Virieux J．P－SV wave propagation in heterogeneous media：velocity－stress finite－difference method．Geophysics，1986，51（4）：889－901．

［16］ 馮英杰，楊長(zhǎng)春，吳萍．地震波有限差分模擬綜述．地球物理學(xué)進(jìn)展，2007，22（2）：487－491．Feng Y J，Yang C C，Wu P．The review of the finitedifference elastic wave motion modeling．Progress in Geophysics（in Chinese），2007，22（2）：487－491．

［17］ Saenger E H，Gold N，Shapiro S A．Modeling the propagation of elastic waves using a modified finite－difference grid．Wave Motion，2000，31（1）：77－92．

［18］ Saenger E H，Bohlen T．Finite－difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid．Geophysics，2004，69（2）：583－591．

［19］ 陳浩，王秀明，趙海波．旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分及其完全匹配層吸收邊界條件．科學(xué)通報(bào)，2006，51（17）：1985－1994．Chen H，Wang X M，Zhao H B．Rotated staggered grid finite－difference and the perfectly matched layer．Chinese Science Bulletin（in Chinese），2006，51（17）：1985－1994．

［20］ 張魯新，符力耘，裴正林．不分裂卷積完全匹配層與旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分在孔隙彈性介質(zhì)模擬中的應(yīng)用．地球物理學(xué)報(bào)，2010，53（10）：2470－2483．Zhang L X，F(xiàn)u L Y，Pei Z L．Finite difference modeling of Biot＇s poroelastic equations with unsplit convolutional PML and rotated staggered grid．Chinese J．Geophys．（in Chinese），2010，53（10）：2470－2483．

［21］ 杜啟振，孫瑞艷，張強(qiáng)．組合邊界條件下二維三分量TTI介質(zhì)波場(chǎng)數(shù)值模擬．石油地球物理勘探，2011，46（2）：187－195．Du Q Z，Sun R Y，Zhang Q．Numerical simulation of threecomponent elastic wavefield in 2DTTI media in the condition of the combined boundary．Oil Geophysical Prospecting（in Chinese），2011，46（2）：187－195．

［22］ Liu Y，Sen M K．An implicit staggered－grid finite－difference method for seismic modelling．Geophysical Journal International，2009，179（1）：459－474．

［23］ Igel H，Mora P，Riollet B．Anisotropic wave propagation through finite－difference grids．Geophysics，1995，60（4）：1203－1216．

［24］ Berenger J P．A perfectly matched layer for the absorption of electromagnetic waves．Journal of Computational Physics，1994，114（2）：185－200．

［25］ Collino F，Tsogka C．Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media．Geophysics，2001，66（1）：294－307

［26］ Zhang Z J，Teng J W，Wan Z C．Establishment of parameteric equations of seismic wave fronts in 2－D anisotropic media．Chinese Science Bulletin，1994，39（22）：1890－1894．

［27］ 裴正林，王尚旭．任意傾斜各向異性介質(zhì)中彈性波波場(chǎng)交錯(cuò)網(wǎng)格高階有限差分法模擬．地震學(xué)報(bào)，2005，27（4）：441～451．Pei Z L，Wang S X．A staggered－grid high－order finite－difference modeling for elastic wave field in arbitrary title aniotropic media．Acta Seismologica Sinica（in Chinese），2005，27（4）：441～451．

Rotated staggered grid high－order finite－difference numerical modeling for wave propagation in viscoelastic TTI media

YAN Hong－Yong1，2，LIU Yang1，2＊

1 State Key Laboratory of Petroleum Resource and Prospecting，China University of Petroleum，Beijing102249，China 2 CNPC Key Laboratory of Geophysical Exploration，China University of Petroleum，Beijing102249，China

On the basis of Carcione′s theories of viscoelasticity and anisotropy，two－dimensional，three－component，first－order velocity－stress wave equations of viscoelastic tilted transversely isotropic（viscoelastic TTI）media are presented and a rotated staggered grid any－order finitedifference scheme is used to numerically solve the equations．Equations of the perfectly matched layer（PML）are derived for the wave equations in viscoelastic TTI media and the rotated staggered grid any－order finite－difference scheme is also used to solve these equations．Results of numerical modeling indicate that the modeling precision is high and the absorbing boundary condition is good in the viscoelastic TTI media，and high－precision snapshots of wave field and synthetic seismograms can be obtained，and they can reflect the characteristic of viscoelasticity and anisotropy．

Viscoelasticity，Anisotropy，Rotated staggered－grid，F(xiàn)inite difference，Perfectly match layer

P631收修定稿2011－03－17，2012－03－16收修定稿

國(guó)家自然科學(xué)基金項(xiàng)目（41074100）、國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863）項(xiàng)目（2007AA06Z218）和教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃（NCET－10－0812）聯(lián)合資助．

嚴(yán)紅勇，男，1982年生，中國(guó)石油大學(xué)（北京）博士研究生，主要從事地震波傳播與成像、地震波吸收衰減及地震資料高分辨率處理等方面的研究工作．E－mail：yanhongyong＠163．com

＊通信作者劉洋，E－mail：wliuyang＠vip．sina．com

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

10．6038／j．issn．0001－5733．2012．04．031．

Yan H Y，Liu Y．Rotated staggered grid high－order finite－difference numerical modeling for wave propagation in viscoelastic TTI media．Chinese J．Geophys．（in Chinese），2012，55（4）：1354－1365，doi：10．6038／j．issn．0001－5733．2012．04．031．

10．6038／j．issn．0001－5733．2012．04．031

（本文編輯 汪海英）