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

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

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

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

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

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

TTI介質(zhì)；分裂式完全匹配層；縱橫波耦合方程；一階波動方程；高階交錯網(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ù)值模擬是各向異性逆時偏移和各向異性全波形反演的基礎(chǔ)［1］。1986年，THOMSEN［2］提出了在大多數(shù)情況下地球介質(zhì)為弱各向異性（各向異性參數(shù)為10%～20%）的觀點，并在弱各向異性假設(shè)下提出弱各向異性橫向各向同性（TI）理論，即用3個各向異性參數(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波波動方程，該方程雖然采用了令橫波速度為0的TTI聲波近似，但仍然是一個P－SV波耦合的波動方程，因為各向異性波場傳播過程中偽SV波仍然存在［5］，而且在角度劇變（即各向異性的對稱軸在空間變化較快）情況下波場傳播存在嚴重的不穩(wěn)定性問題［6］。為了解決TTI介質(zhì)聲波近似各向異性波場傳播存在的不穩(wěn)定性問題，F(xiàn)LETCHER等［7］引入適當(dāng)?shù)臋M波分量來減小角度劇變時產(chǎn)生的不穩(wěn)定性，但是同時也帶來了較強的偽SV波假象。DUVENECK等［8－9］從Hooke定律和運動學(xué)方程出發(fā)推導(dǎo)得到新的TTI介質(zhì)方程；ZHANG等［10－11］從VTI介質(zhì)彈性波方程出發(fā)采用自共軛旋轉(zhuǎn)方式推導(dǎo)出一個自共軛3DTTI波動方程。但上述方程仍然存在偽SV波假象且在角度劇變時不穩(wěn)定。FOWLER等［12］從特征值矩陣分析的角度對各種TTI介質(zhì)耦合方程的精度和穩(wěn)定性進行了綜合分析。BUBE等［13－14］討論了二階TI介質(zhì)波動方程在TTI介質(zhì)聲波近似情況下不穩(wěn)定性產(chǎn)生的原因并推導(dǎo)了新的聲波和彈性波TI介質(zhì)一階波動方程。

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

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

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

1 方法原理

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

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］，此時推導(dǎo)得到的2DTTI介質(zhì)方程為［4］：

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

式中：φ為三維方位角。當(dāng)θ＝0時，（3）式退化為3D VTI方程（垂直對稱軸的TI方程）；當(dāng)θ＝90°時，（3）式退化為3DHTI方程（水平對稱軸的TI方程）；當(dāng)三維方位角φ＝0時，（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ì)方程。但是這種方程在實際計算中往往會產(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波波動方程為：

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

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

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

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

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

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

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

從t＝1s時刻的波場快照（圖2a，圖2b）可以看出，ZHOU等［4］方程基本上還能保持穩(wěn)定傳播，F(xiàn)LETCHER等［7］方程由于引入了橫波分量，產(chǎn)生了偽橫波（SV波）噪聲（圖2b中紅色箭頭所示）。而從t＝1．5s時刻的波場快照（圖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］各向異性波動方程是一種穩(wěn)定的TTI介質(zhì)波動方程，我們討論的二維和三維TTI介質(zhì)正演模擬將采用該方程。

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

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

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

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

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

其中，輔助變量形式為：

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

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

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

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

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

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

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

2 數(shù)值計算實例

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

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

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

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

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

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

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

對一系列公式推導(dǎo)研究可以發(fā)現(xiàn)，所有3DTTI介質(zhì)方程的SPML方程都具有相同的形式，差分離散形式也相同，只是PML輔助變量表達形式不一致，其具體形式根據(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ì)地震波場正演數(shù)值模擬和全波形反演研究工作。

國家自然科學(xué)基金項目（41604110）、國土資源部海底礦產(chǎn)資源重點實驗室開放基金項目（KLMMR－2015－A－14）、國家高技術(shù)研究發(fā)展計劃（863計劃）項目（2013AA092501）和國家科技重大專項（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）．