TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究

張 衡，劉 洪，李 博，丁仁偉，李麗青，李福元

（1．國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局，廣東廣州510075；2．中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京100029；3．中國石油化工股份有限公司石油物探技術(shù)研究院，江蘇南京211103；4．山東科技大學(xué)地球科學(xué)與工程學(xué)院，山東青島266590）

TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究

張 衡1，劉 洪2，李 博3，丁仁偉4，李麗青1，李福元1

（1．國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室，中國地質(zhì)調(diào)查局廣州海洋地質(zhì)調(diào)查局，廣東廣州510075；2．中國科學(xué)院地質(zhì)與地球物理研究所，中國科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京100029；3．中國石油化工股份有限公司石油物探技術(shù)研究院，江蘇南京211103；4．山東科技大學(xué)地球科學(xué)與工程學(xué)院，山東青島266590）

針對(duì)傾斜橫向各向同性（TTI）介質(zhì)聲波波動(dòng)方程的特點(diǎn)，研究了TTI介質(zhì)分裂式完全匹配層（Split perfectly matched layer，SPML）吸收邊界條件。首先對(duì)常見的幾種TTI介質(zhì)聲波波動(dòng)方程進(jìn)行了歸納，并從TTI介質(zhì)波場(chǎng)傳播穩(wěn)定性的角度進(jìn)行對(duì)比分析，結(jié)果表明，引入橫波分量的TTI介質(zhì)縱橫波耦合方程適用于TTI介質(zhì)。然后從TTI介質(zhì)縱橫波耦合二階波動(dòng)方程出發(fā)，推導(dǎo)得到其一階波動(dòng)方程的形式，進(jìn)而推導(dǎo)出一階波動(dòng)方程形式的SPML波動(dòng)方程，并給出了高階交錯(cuò)網(wǎng)格有限差分算法的具體實(shí)現(xiàn)過程。數(shù)值模擬結(jié)果表明，SPML吸收邊界條件能達(dá)到很好的人工邊界反射吸收效果，相比優(yōu)化海綿吸收邊界條件，其人工邊界反射吸收效果更好。

TTI介質(zhì)；分裂式完全匹配層；縱橫波耦合方程；一階波動(dòng)方程；高階交錯(cuò)網(wǎng)格有限差分

Keywords：TTI media，split perfectly matched layer，P－wave and SV－wave coupled equation，first－order wave equation，high－order staggered－grid finite－difference

各向異性正演數(shù)值模擬是各向異性逆時(shí)偏移和各向異性全波形反演的基礎(chǔ)［1］。1986年，THOMSEN［2］提出了在大多數(shù)情況下地球介質(zhì)為弱各向異性（各向異性參數(shù)為10%～20%）的觀點(diǎn)，并在弱各向異性假設(shè)下提出弱各向異性橫向各向同性（TI）理論，即用3個(gè)各向異性參數(shù)——Thomsen各向異性參數(shù)ε（縱波各向異性參數(shù)）、δ（變異系數(shù)）和γ（橫波各向異性參數(shù)）表述P波、SV波及SH波特征。此后30年，關(guān)于各向異性介質(zhì)的研究基本都是基于弱各向異性理論。

1996年，TSVANKIN［3］基于THOMSEN的弱各向異性理論推導(dǎo)得到P－SV波VTI介質(zhì)精確相速度頻散方程。2006年，ZHOU等［4］從TSVANKIN相速度頻散方程出發(fā)，推導(dǎo)得到TTI介質(zhì)qP波波動(dòng)方程，該方程雖然采用了令橫波速度為0的TTI聲波近似，但仍然是一個(gè)P－SV波耦合的波動(dòng)方程，因?yàn)楦飨虍愋圆▓?chǎng)傳播過程中偽SV波仍然存在［5］，而且在角度劇變（即各向異性的對(duì)稱軸在空間變化較快）情況下波場(chǎng)傳播存在嚴(yán)重的不穩(wěn)定性問題［6］。為了解決TTI介質(zhì)聲波近似各向異性波場(chǎng)傳播存在的不穩(wěn)定性問題，F(xiàn)LETCHER等［7］引入適當(dāng)?shù)臋M波分量來減小角度劇變時(shí)產(chǎn)生的不穩(wěn)定性，但是同時(shí)也帶來了較強(qiáng)的偽SV波假象。DUVENECK等［8－9］從Hooke定律和運(yùn)動(dòng)學(xué)方程出發(fā)推導(dǎo)得到新的TTI介質(zhì)方程；ZHANG等［10－11］從VTI介質(zhì)彈性波方程出發(fā)采用自共軛旋轉(zhuǎn)方式推導(dǎo)出一個(gè)自共軛3DTTI波動(dòng)方程。但上述方程仍然存在偽SV波假象且在角度劇變時(shí)不穩(wěn)定。FOWLER等［12］從特征值矩陣分析的角度對(duì)各種TTI介質(zhì)耦合方程的精度和穩(wěn)定性進(jìn)行了綜合分析。BUBE等［13－14］討論了二階TI介質(zhì)波動(dòng)方程在TTI介質(zhì)聲波近似情況下不穩(wěn)定性產(chǎn)生的原因并推導(dǎo)了新的聲波和彈性波TI介質(zhì)一階波動(dòng)方程。

TTI介質(zhì)正演時(shí)，吸收邊界條件至關(guān)重要，因?yàn)槲者吔鐥l件極大地影響正演模擬效果。地下介質(zhì)被視為無限半空間，而計(jì)算區(qū)域有限，因此在人為截取的有限空間內(nèi)求解波動(dòng)方程會(huì)導(dǎo)致很強(qiáng)的人工邊界反射，必須采用吸收邊界條件對(duì)人工邊界反射進(jìn)行有效吸收。目前常用的吸收邊界條件包括單程波旁軸近似吸收邊界條件［15］和衰減吸收邊界條件［16－17］兩大類。CLAYTON等［15］提出的單程波旁軸近似吸收邊界條件應(yīng)用廣泛，在小角度反射時(shí)吸收效果較好，但對(duì)大角度入射的波吸收效果較差。CERJAN等［16］提出了添加阻尼層的海綿吸收邊界條件，該邊界條件不受波動(dòng)方程形式的限制，很容易實(shí)現(xiàn)，但是邊界吸收的效果仍較差。BORDING［18］對(duì)CERJAN等［16］提出的阻尼系數(shù)進(jìn)行了優(yōu)化，計(jì)算得到優(yōu)化的海綿吸收系數(shù)，但是計(jì)算效率相對(duì)CERJAN等的方法有所降低。BERENGER［17］最早在電磁學(xué)領(lǐng)域提出了完全匹配層（PML）吸收邊界條件，能達(dá)到完美的邊界吸收效果。CHEW等［19］將PML吸收邊界條件用于地震波數(shù)值模擬中，取得了很好的應(yīng)用效果。BéCACHE等［20］對(duì)各向異性介質(zhì)的PML吸收邊界條件進(jìn)行了深入研究，指出各向異性介質(zhì)應(yīng)用PML吸收邊界條件時(shí)存在不穩(wěn)定性。PML吸收邊界條件是目前應(yīng)用效果最好的吸收邊界條件，代表吸收邊界技術(shù)研究的前沿發(fā)展方向，但是其公式推導(dǎo)和編程實(shí)現(xiàn)復(fù)雜，需要對(duì)不同的波動(dòng)方程推導(dǎo)不同的PML方程形式。

作為PML吸收邊界條件技術(shù)研究的一個(gè)重要方面，分裂式PML（SPML）吸收邊界條件已經(jīng)得到廣泛應(yīng)用。COLLINO等［21］提出二維彈性波介質(zhì)一階速度 應(yīng)力方程的SPML算法；裴正林［22－23］將該算法推廣到三維彈性波介質(zhì)，實(shí)現(xiàn)了基于SPML吸收邊界條件的三維彈性波介質(zhì)波動(dòng)方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬。

本文研究的TTI介質(zhì)時(shí)間域SPML吸收邊界條件能很好地吸收邊界反射，實(shí)現(xiàn)TTI介質(zhì)高精度數(shù)值模擬。

1 方法原理

1．1 TTI介質(zhì)耦合qP波波動(dòng)方程

P－SV波VTI介質(zhì)精確相速度頻散方程為［3］：

式中：f＝1－v2SZ／v2PZ，vSZ表示橫波速度，vPZ表示縱波速度；vP（θph）表示相速度，θph表示相速度傳播角度；ε和δ表示Thomsen各向異性參數(shù)［2］。

在推導(dǎo)過程中直接引入TTI介質(zhì)聲波近似（TTI介質(zhì)聲波近似情況下直接令橫波速度vSZ＝0）［24］，此時(shí)推導(dǎo)得到的2DTTI介質(zhì)方程為［4］：

式中：θ為對(duì)稱軸傾角；p為地震波場(chǎng)；q為輔助波場(chǎng)；H1和H2為輔助變量，其3DTTI方程形式為：

式中：φ為三維方位角。當(dāng)θ＝0時(shí)，（3）式退化為3D VTI方程（垂直對(duì)稱軸的TI方程）；當(dāng)θ＝90°時(shí)，（3）式退化為3DHTI方程（水平對(duì)稱軸的TI方程）；當(dāng)三維方位角φ＝0時(shí)，（3）式簡化為2DTTI方程。VTI介質(zhì)和HTI介質(zhì)都可以看作是TTI介質(zhì)的某種特殊情形。

本文將ZHOU等［4］、DUVENECK等［8－9］和ZHANG等［10－11］直接采用vSZ＝0的TTI介質(zhì)聲波近似思想推導(dǎo)得到的TTI介質(zhì)方程歸納為基于TTI介質(zhì)聲波近似的TTI介質(zhì)方程。但是這種方程在實(shí)際計(jì)算中往往會(huì)產(chǎn)生數(shù)值不穩(wěn)定問題，尤其是在角度劇變的情況下［6］。FLETCHER等［7］完全放棄了采用TTI介質(zhì)聲波近似的思路，提出在方程推導(dǎo)中保留橫波分量。引入橫波分量的方式能有效消除TTI介質(zhì)聲波近似產(chǎn)生的不穩(wěn)定問題，F(xiàn)LETCHER等［7］引入橫波分量推導(dǎo)得到的2D TTI介質(zhì)縱橫波耦合qP波波動(dòng)方程為：

TTI介質(zhì)波場(chǎng)模擬首先要滿足物理穩(wěn)定性條件，即P－SV波VTI介質(zhì)精確相速度頻散方程（（1）式）中的相速度的平方不能小于0［26］：

由（5）式推導(dǎo)得到TI介質(zhì)波場(chǎng)模擬中的物理穩(wěn)定性條件為：

在TTI介質(zhì)聲波近似情況下，vSZ＝0，則f＝1，由（6）式推導(dǎo)得到TTI介質(zhì)聲波近似情況下的物理穩(wěn)定性條件為：

在變速介質(zhì)情況下，由于對(duì)變系數(shù)求導(dǎo)會(huì)帶來微分算子的附加項(xiàng)，影響算子的穩(wěn)定性，因此，我們借助對(duì)稱算子理論，將原始偏微分方程適當(dāng)對(duì)稱化，得到與常系數(shù)偏微分方程一樣穩(wěn)定的微分算子［11，27－28］。

經(jīng)實(shí)際測(cè)試，TTI介質(zhì)波場(chǎng)傳播不穩(wěn)定的情況主要由角度劇變產(chǎn)生，包括對(duì)稱軸傾角θ或者三維方位角φ快速變化的區(qū)域［6］。本文設(shè)計(jì)了一個(gè)楔形2DTTI介質(zhì)模型用于算法穩(wěn)定性測(cè)試（圖1），模型的主要特點(diǎn)是角度劇變（圖1d），因此能很好地測(cè)試各向異性波動(dòng)方程算法的穩(wěn)定性。模型橫向網(wǎng)格數(shù)和縱向網(wǎng)格數(shù)都為300，橫向網(wǎng)格間距和縱向網(wǎng)格間距都為10m。

圖1 楔形2DTTI介質(zhì)模型各向異性參數(shù)變化情況

數(shù)值模擬時(shí)，雷克震源子波主頻取25Hz，傳播時(shí)間為3s，時(shí)間采樣間隔為0．000 5s，吸收邊界條件采用優(yōu)化海綿吸收邊界條件［18］，分別采用TTI介質(zhì)聲波近似方程（ZHOU等［4］方程）和引入橫波分量的TTI介質(zhì)縱橫波耦合方程（FLETCHER等［7］方程）來進(jìn)行數(shù)值模擬。

從t＝1s時(shí)刻的波場(chǎng)快照（圖2a，圖2b）可以看出，ZHOU等［4］方程基本上還能保持穩(wěn)定傳播，F(xiàn)LETCHER等［7］方程由于引入了橫波分量，產(chǎn)生了偽橫波（SV波）噪聲（圖2b中紅色箭頭所示）。而從t＝1．5s時(shí)刻的波場(chǎng)快照（圖2c，圖2d）可以看出，ZHOU等［4］方程由于楔形界面處角度劇變已經(jīng)產(chǎn)生不穩(wěn)定現(xiàn)象，F(xiàn)LETCHER等［7］方程由于引入了橫波分量，能減小不穩(wěn)定情況的產(chǎn)生，仍然能保持穩(wěn)定傳播。因此FLETCHER等［7］各向異性波動(dòng)方程是一種穩(wěn)定的TTI介質(zhì)波動(dòng)方程，我們討論的二維和三維TTI介質(zhì)正演模擬將采用該方程。

本文對(duì)常見的TTI介質(zhì)波動(dòng)方程進(jìn)行了歸納，認(rèn)為基于TTI介質(zhì)精確相速度頻散方程推導(dǎo)得到的TTI介質(zhì)波動(dòng)方程可以分為兩類：一類為TTI介質(zhì)聲波近似方程；另一類為TTI介質(zhì)縱橫波耦合方程。其中TTI介質(zhì)聲波近似方程可以視為TTI介質(zhì)縱橫波耦合方程在橫波速度為0時(shí)的某種特例。從TTI介質(zhì)波場(chǎng)傳播穩(wěn)定性的角度對(duì)這兩類方程進(jìn)行對(duì)比分析可知，引入橫波分量的TTI介質(zhì)縱橫 波耦合方程適用于TTI介質(zhì)。

圖2 楔形2DTTI介質(zhì)模型波場(chǎng)快照

1．2 2DTTI介質(zhì)時(shí)間域分裂式完全匹配層（SPML）波動(dòng)方程推導(dǎo)

基于FLETCHER等［7］的二階2DTTI介質(zhì)縱橫波耦合波動(dòng)方程（（4）式）來推導(dǎo)2DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程。TTI介質(zhì)波動(dòng)方程的特點(diǎn)是含有比較復(fù)雜的交叉導(dǎo)數(shù)項(xiàng)，如何有效處理TTI介質(zhì)方程中的交叉導(dǎo)數(shù)項(xiàng)是PML吸收邊界條件研究的難點(diǎn)。本文的思路是針對(duì)TTI介質(zhì)波動(dòng)方程的特點(diǎn)對(duì)二階方程中的交叉導(dǎo)數(shù)項(xiàng)進(jìn)行分裂，得到其一階方程的形式，進(jìn)而推導(dǎo)得到TTI介質(zhì)一階SPML波動(dòng)方程。

引入4個(gè)輔助波場(chǎng)項(xiàng)ψx，ψz，Ωx，Ωz，從FLETCHER等［7］的二階2DTTI介質(zhì)波動(dòng)方程出發(fā)，推導(dǎo)得到一階2DTTI介質(zhì)波動(dòng)方程的形式，對(duì)方程進(jìn)行傅里葉變換并引入X和Z方向的頻率域PML拉伸函數(shù)，推導(dǎo)得到2DTTI介質(zhì)頻率域SPML波動(dòng)方程，將此方程反傅里葉變換到時(shí)間域，即可得到2DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程：

其中，輔助變量形式為：

其中，引入了PML衰減因子dx，dz。本文采用的PML衰減因子表達(dá)式為［21］：

式中：Lpml為PML邊界層的厚度；l為計(jì)算點(diǎn)距PML邊界的距離；Rcoeff為理論反射系數(shù)，本文取Rcoeff＝0．000 1。

一系列公式推導(dǎo)研究可以發(fā)現(xiàn)，所有2DTTI介質(zhì)的SPML方程都具有相同的形式，差分離散形式也相同，只是PML輔助變量表達(dá)形式不一致，其具體形式根據(jù)不同的TTI方程形式而定。本文只討論二維TTI介質(zhì)情形，但是這種方法也可以擴(kuò)展到三維TTI介質(zhì)（詳細(xì)推導(dǎo)見附錄A），在此不再展開討論。

1．3 2DTTI介質(zhì)SPML波動(dòng)方程數(shù)值實(shí)現(xiàn)

采用高階交錯(cuò)網(wǎng)格有限差分法來實(shí)現(xiàn)1．2節(jié)推導(dǎo)的2DTTI介質(zhì)SPML波動(dòng)方程。對(duì)（8）式進(jìn)行高階交錯(cuò)網(wǎng)格有限差分離散，具體實(shí)現(xiàn)公式如下：

從二維SPML數(shù)值實(shí)現(xiàn)示意圖（圖3）可以看出，二維PML區(qū)域分為4個(gè)邊和4個(gè)角共8塊區(qū)域，具體實(shí)現(xiàn)時(shí)，我們需要分別考慮各PML區(qū)域波場(chǎng)的衰減。理論上，無吸收邊界條件具有最強(qiáng)的邊界反射，而PML吸收層數(shù)越多，邊界反射吸收效果越好。但是增加的PML吸收層數(shù)也帶來了更多的計(jì)算量，從而在一定程度上降低了計(jì)算效率，因此在實(shí)際應(yīng)用中需要在PML吸收層數(shù)與計(jì)算效率之間尋求平衡。

圖3 二維SPML數(shù)值實(shí)現(xiàn)示意圖解

2 數(shù)值計(jì)算實(shí)例

2．1 簡單均勻2DTTI介質(zhì)模型

為了驗(yàn)證本文提出的TTI介質(zhì)SPML吸收邊界條件的有效性，首先采用簡單均勻2DTTI介質(zhì)模型進(jìn)行SPML點(diǎn)源響應(yīng)測(cè)試。模型橫向網(wǎng)格數(shù)和縱向網(wǎng)格數(shù)都為256，橫向網(wǎng)格間距和縱向網(wǎng)格間距都為10m，速度為3 000m／s。TTI介質(zhì)的各向異性參數(shù)ε＝0．3，δ＝0．1，對(duì)稱軸傾角為45°。

均勻TTI介質(zhì)模型的單道記錄能很好地測(cè)試不同PML邊界層條件下的邊界反射吸收效果。數(shù)值模擬時(shí)，雷克震源子波主頻取25Hz，震源位置設(shè)置在模型的中心，單檢波點(diǎn)取在震源位置左側(cè)50個(gè)網(wǎng)格點(diǎn)處，傳播時(shí)間為1．5s，時(shí)間采樣間隔為0．001 5s，采用12階高階交錯(cuò)網(wǎng)格有限差分格式進(jìn)行模擬。采用SPML吸收邊界條件，PML邊界吸收層數(shù)分別設(shè)為10層、20層和30層。取3種情形的單道記錄進(jìn)行對(duì)比（圖4）。由圖4可見，PML吸收層數(shù)為10層時(shí)邊界反射吸收效果不理想；而PML吸收層數(shù)為20層時(shí)已經(jīng)取得了比較好的邊界反射吸收效果（圖4a），但從局部放大圖可以看出，仍然存在著微弱的邊界反射；而PML吸收層數(shù)為30層時(shí)邊界反射導(dǎo)致的波場(chǎng)誤差降到僅為10－5，此時(shí)邊界反射基本被消除（圖4b）。從PML吸收層數(shù)為20層時(shí)SPML波動(dòng)方程在t＝0．72s和t＝0．90s的波場(chǎng)快照（圖5）也可以看到微弱的邊界反射。因此本例選定PML邊界吸收層數(shù)為30層。圖6a到圖6f分別給出了SPML波動(dòng)方程在PML吸收層數(shù)為30層時(shí)各個(gè)時(shí)刻的波場(chǎng)快照。由圖6可見，SPML對(duì)邊界反射吸收效果明顯。

2．2 復(fù)雜BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型

以復(fù)雜BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型為例，對(duì)比TTI介質(zhì)SPML吸收邊界條件與優(yōu)化海綿吸收邊界條件對(duì)邊界反射的吸收效果［18］。本文數(shù)值模擬時(shí)采用截取的部分BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型（圖7），該模型主要模擬了近海岸的地質(zhì)構(gòu)造，特點(diǎn)是左部有一個(gè)被TTI各向異性地層包圍的高速各向同性鹽丘，模型參數(shù)包括縱波速度vPZ，Thomsen縱波各向異性參數(shù)ε，Thomsen變異系數(shù)δ，對(duì)稱軸角度θ，模型橫向網(wǎng)格數(shù)nx＝601，縱向網(wǎng)格數(shù)nz＝361，橫向網(wǎng)格間距和縱向網(wǎng)格間距都為6．25m，其中海水層為各向同性介質(zhì)。數(shù)值模擬時(shí)，雷克震源子波主頻取25Hz，采用SPML吸收邊界條件，PML邊界設(shè)為30層，傳播時(shí)間為2s，時(shí)間采樣間隔為0．000 4s，采用12階高階交錯(cuò)網(wǎng)格有限差分格式進(jìn)行模擬。因?yàn)檎鹪粗糜诟飨蛲缘暮Ｋ畬又校琓TI數(shù)值模擬時(shí)不會(huì)出現(xiàn)影響qP波場(chǎng)傳播的偽橫波（SV波）噪聲，此種情況下不需要采取在震源處添加各向同性薄層消除偽SV波噪聲的策略［5，29］。從BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型的正演單炮記錄（圖8）來看，SPML消除邊界反射效果良好（圖8a），而優(yōu)化海綿吸收邊界條件仍然有著比較強(qiáng)的邊界反射（圖8b中紅色箭頭標(biāo)示），SPML吸收邊界條件明顯優(yōu)于優(yōu)化海綿吸收邊界條件的邊界反射吸收效果。

圖4 均勻2DTTI介質(zhì)模型不同PML吸收層數(shù)單道記錄

圖5 均勻2DTTI介質(zhì)模型采用2DTTI介質(zhì)SPML波動(dòng)方程在PML吸收層數(shù)為20層時(shí)不同時(shí)刻的正演波場(chǎng)快照

圖6 均勻2DTTI介質(zhì)模型采用2DTTI介質(zhì)SPML波動(dòng)方程在PML吸收層數(shù)為30層時(shí)不同時(shí)刻的正演波場(chǎng)快照

圖7 截取的部分BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型各向異性參數(shù)變化情況

圖8 BP2007 2DTTI介質(zhì)海洋標(biāo)準(zhǔn)模型采用不同邊界條件正演的單炮記錄

3 結(jié)束語

本文對(duì)常見的幾種TTI介質(zhì)波動(dòng)方程進(jìn)行了歸納，并從波場(chǎng)傳播穩(wěn)定性的角度進(jìn)行分析，結(jié)果表明引入橫波分量的TTI介質(zhì)縱橫波耦合波動(dòng)方程適用于TTI介質(zhì)?；赥TI介質(zhì)縱橫波耦合二階波動(dòng)方程，推導(dǎo)得到一階方程形式的TTI介質(zhì)SPML波動(dòng)方程，研究發(fā)現(xiàn)，所有2DTTI介質(zhì)的SPML方程都具有相同的形式，差分離散形式也相同，只是PML輔助變量表達(dá)形式不一致，其根據(jù)不同的TTI方程形式而定。數(shù)值模擬結(jié)果表明，本文研究的SPML吸收邊界條件能達(dá)到很好的人工邊界反射吸收效果。本文方法可以進(jìn)一步推廣到三維TTI介質(zhì)的高精度模擬中。至于2DVTI或3DVTI介質(zhì)情形，由于對(duì)稱軸角度為0，波場(chǎng)傳播較穩(wěn)定，此時(shí)適合采用噪聲相對(duì)較少的VTI介質(zhì)聲波近似方程，因?yàn)門TI介質(zhì)縱橫波耦合方程雖然由于引入了橫波分量解決了地震波在各向異性介質(zhì)中傳播的不穩(wěn)定性問題，但是也帶來了更多的橫波噪聲。VTI介質(zhì)波動(dòng)方程的特點(diǎn)是不含比較復(fù)雜的交叉導(dǎo)數(shù)項(xiàng)，可以采用計(jì)算效率更高的非SPML吸收邊界條件來實(shí)現(xiàn)。

［1］ WARNER M，RATCLIFE A，NANGOO T，et al．Ani－sotropic 3Dfull－waveform inversion［J］．Geophysics，2013，78（2）：R59－R80

［2］ THOMSEN L．Weak elastic anisotropy［J］．Geophysics，1986，51（10）：1954－1966

［3］ TSVANKIN I．P－wave signatures and notation for transversely isotropic media：an overview［J］．Geophysics，1996，61（2）：467－483

［4］ ZHOU H B，ZHANG G Q，BLOOR R．An anisotropic acoustic wave equation for modeling and migration in 2DTTI media［J］．Expanded Abstracts of 76thAnnual Internat SEG Mtg，2006：194－198

［5］ GRECHKA V，ZHANG L，RECTOR J W．Shear waves in acoustic anisotropic media［J］．Geophysics，2004，69（2）：576－582

［6］ BAKKER P M，DUVENECK E．Stability analysis for acoustic wave propagation in tilted TI media by finite differences［J］．Geophysical Journal International，2011，185（2）：911－921

［7］ FLETCHER R P，DU X，F(xiàn)OWLER P J．Reverse time migration in tilted transversely isotropic（TTI）media［J］．Geophysics，2009，74（6）：179－187

［8］ DUVENECK E，MILCIK P，BAKKER P M．Acoustic VTI wave equations and their application for anisotropic reverse－time migration［J］．Expanded Abstracts of 78thAnnual Internat SEG Mtg，2008：2186－2190

［9］ DUVENECK E，BAKKER P M．Stable P－wave modeling for reverse－time migration in tilted TI media［J］．Geophysics，2011，76（2）：65－75

［10］ ZHANG Y，ZHANG H Z．A stable TTI reverse time migration and its implementation［J］．Expanded Abstracts of 79thAnnual Internat SEG Mtg，2009：2794－2798

［11］ ZHANG Y，ZHANG H Z，ZHANG G Q．A stable TTI reverse time migration and its implementation［J］．Geophysics，2011，76（3）：WA3－WA11

［12］ FOWLER P J，DU X，F(xiàn)LETCHER R P．Coupled equations for reverse time migration in transversely isotropic media［J］．Geophysics，2010，75（1）：S11－S22

［13］ BUBE K P，NEMETH T，STEFANI J P，et al．Firstorder systems for elastic and acoustic variable－tilt TI media［J］．Geophysics，2012，77（5）：T157－T170

［14］ BUBE K P，NEMETH T，STEFANI J P，et al．On the instability in second－order systems for acoustic VTI and TTI media［J］．Geophysics，2012，77（5）：T171－T186

［15］ CLAYTON R，ENQUIST B．Absorbing boundary conditions for acoustic and elastic wave equations［J］．Bulletin of the Seismological Society of America，1977，67（6）：1529－1540

［16］ 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

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

［18］ BORDING R P．Finite difference modeling－nearly optimal sponge boundary conditions［J］．Expanded Abstracts of 74thAnnual Internat SEG Mtg，2004：1921－1924

［19］ CHEW W C，LIU Q H．Perfectly matched layers for elastodynamics：a new absorbing boundary condition［J］．Journal of Computational Acoustics，1996，4（4）：341－359

［20］ BéCACHE E，F(xiàn)AUQUEUX S，JOLY P．Stability of perfectly matched layers，group velocities and anisotropic waves［J］．Journal of Computational Physics，2003，188（2）：399－433

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

［22］ 裴正林．三維各向異性介質(zhì)中彈性波方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬［J］．石油大學(xué)學(xué)報(bào)（自然科學(xué)版），2004，28（5）：23－29

PEI Z L．Three－dimensional numerical simulation of elastic wave propagation in 3－D anisotropic media withstaggered－grid high－order difference method［J］．Journal of the University of Petroleum，2004，28（5）：23－29

［23］ 裴正林．三維各向同性介質(zhì)彈性波方程交錯(cuò)網(wǎng)格高階有限差分法模擬［J］．石油物探，2005，44（4）：308－315

PEI Z L．Numerical simulation of elastic wave equation in 3－D isotropic media with staggered－grid high－order difference method［J］．Geophysical Prospecting for Petroleum，2005，44（4）：308－315

［24］ ALKHALIFAH T．Acoustic approximations for processing in transversely isotropic media［J］．Geophysics，1998，63（2）：623－631

［25］ 李博，李敏，劉紅偉，等．TTI介質(zhì)有限差分逆時(shí)偏移的穩(wěn)定性探討［J］．地球物理學(xué)報(bào)，2012，55（4）：1366－1375

LI B，LI M，LIU H W，et al．Stability of reverse time migration in TTI media［J］．Chinese Journal of Geophysics，2012，55（4）：1366－1375

［26］ 李博．地震偏移成像和GPU超算技術(shù)［D］．北京：中國科學(xué)院大學(xué)，2011

LI B．Seismic migration and GPU supercomputing technology［D］．Beijing：University of Chinese Academy of Sciences，2011

［27］ CLAERBOUT J F．Imaging the earth’s interior［M］．Palo Alto，California：Blackwell Scientific Publications，1985：1－398

［28］ 劉洪，王秀閩，曾銳，等．單程波算子積分解的象征表示［J］．地球物理學(xué)進(jìn)展，2007，22（2）：463－471

LIU H，WANG X M，ZENG R，et al．Symbol description to integral solution of one－way wave operator［J］．Progress in Geophysics，2007，22（2）：463－471

［29］ 張衡，劉洪，唐祥德，等．基于平均導(dǎo)數(shù)優(yōu)化方法的VTI介質(zhì)頻率空間域正演［J］．地球物理學(xué)報(bào)，2015，58（9）：3306－3316

ZHANG H，LIU H，TANG X D，et al．Forward modeling of VTI media in the frequency－space domain based on an average－derivative optimal method［J］．Chinese Journal of Geophysics，2015，58（9）：3306－3316

附錄A 3DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程

3DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程的推導(dǎo)基于FLETCHER等［7］的二階3DTTI介質(zhì)縱橫波耦合波動(dòng)方程（（4）式）。與二維TTI介質(zhì)相比，三維TTI介質(zhì)時(shí)間域SPML波動(dòng)方程的推導(dǎo)首先需要從FLETCHER等［7］的二階3DTTI介質(zhì)波動(dòng)方程出發(fā)，推導(dǎo)得到一階3DTTI介質(zhì)波動(dòng)方程的形式，然后對(duì)方程進(jìn)行傅里葉變換并引入x，y和z 3個(gè)方向的頻率域PML拉伸函數(shù)，將推導(dǎo)得到的3DTTI介質(zhì)頻率域SPML波動(dòng)方程反傅里葉變換，即可得到3DTTI介質(zhì)時(shí)間域SPML波動(dòng)方程：

推導(dǎo)的輔助變量形式為：

PML波動(dòng)方程推導(dǎo)過程中引入的頻率域PML拉伸函數(shù)ξx，ξy，ξz表示為：

引入的PML衰減因子dx，dy，dz表達(dá)式為：

式中：Lpml是PML邊界層的厚度；l是計(jì)算點(diǎn)距PML邊界的距離；Rcoeff是理論反射系數(shù)值，取0．000 1。

對(duì)一系列公式推導(dǎo)研究可以發(fā)現(xiàn)，所有3DTTI介質(zhì)方程的SPML方程都具有相同的形式，差分離散形式也相同，只是PML輔助變量表達(dá)形式不一致，其具體形式根據(jù)不同的TTI方程形式而定。

（編輯：陳 杰）

The research on split PML absorbing boundary conditions of acoustic equation for TTI media

ZHANG Heng1，LIU Hong2，LI Bo3，DING Renwei4，LI Liqing1，LI Fuyuan1

（1．MLR Key Laboratory of Marine Mineral Resources，Guangzhou Marine Geological Survey，China Geological Survey，Guangzhou510075，China；2．Key Laboratory of Petroleum Resources Research，Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing100029，China；3．Sinopec Geophysical Research Institute，Nanjing211103，China；4．College of Earth Science and Engineering，Shandong University of Science and Technology，Qingdao 266590，China）

We study split perfectly matched layer（SPML）absorbing boundary conditions for TTI media based on the characteristics of TTI media acoustic wave equation．We firstly summarize several common TTI media acoustic wave equations and compare these wave equations in detail from the perspective of wavefield propagation stability for TTI media．We conclude that the TTI acoustic approximation equation introduced SV－wave component is applicable for TTI media．Then we derive the TTI media first－order acoustic wave equation based on the SPML boundary condition from the TTI P－wave and SV－wave coupled second－order wave equation．Afterwards we give the concrete implementation method for high－order staggered－grid finite－difference algorithm．The numerical modeling examples show that the SPML boundary condition can absorb the artificial boundary reflection very well．The boundary reflection absorbing effect of the proposed SPML algorithm is much better compared with the optimal sponge boundary condition．

P631

A

1000－1441（2017）03－0349－13

10．3969／j．issn．1000－1441．2017．03．005

張衡，劉洪，李博，等．TTI介質(zhì)聲波方程分裂式PML吸收邊界條件研究［J］．石油物探，2017，56（3）：349－361

ZHANG Heng，LIU Hong，LI Bo，et al．The research on split PML absorbing boundary conditions of acoustic equation for TTI media［J］．Geophysical Prospecting for Petroleum，2017，56（3）：349－361

2016－02－14；改回日期：2016－04－29。

張衡（1987—），男，博士，工程師，主要從事復(fù)雜各向異性介質(zhì)地震波場(chǎng)正演數(shù)值模擬和全波形反演研究工作。

國家自然科學(xué)基金項(xiàng)目（41604110）、國土資源部海底礦產(chǎn)資源重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目（KLMMR－2015－A－14）、國家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）項(xiàng)目（2013AA092501）和國家科技重大專項(xiàng)（2011ZX05003－003）聯(lián)合資助。

This research is financially supported by the National Natural Science Foundation of China（Grant No．41604110），the Open Fund of MLR Key Laboratory of Marine Mineral Resources（Grant No．KLMMR－2015－A－14），the National High－tech R ＆D Program of China（863Program）（Grant No．2013AA092501）and the National Science and Technology Major Project of China（Grant No．2011ZX05003－003）．