一階彈性波方程數(shù)值模擬中的混合吸收邊界條件

任志明，劉洋＊

1 中國石油大學(xué)（北京）油氣資源與探測國家重點實驗室，北京 102249

2 中國石油大學(xué)（北京）CNPC物探重點實驗室，北京 102249

1 引言

地震數(shù)值模擬是對實際的地質(zhì)、地球物理問題作適當(dāng)?shù)慕?，形成簡化的?shù)學(xué)模型，采用數(shù)值計算方法獲取地震響應(yīng)的過程.地震數(shù)值模擬是了解地震波在地下介質(zhì)中傳播規(guī)律、幫助解釋觀測數(shù)據(jù)的有效手段.然而，在使用計算機進行地震波正演模擬時，由于計算空間的有限性，必然會引入人工截斷邊界.因此，必須對邊界進行適當(dāng)?shù)奶幚?，以消除或減弱這種虛假邊界反射.

數(shù)值模擬中常用的吸收邊界條件有以下三種：（1）單程波邊界條件（Clayton and Engquist，1977；Engquist and Majda，1977；Higdon，1991；Zhou and McMechan，2000；Heidari and Guddati，2006）.這類邊界條件能較好地吸收入射角在一定范圍內(nèi)的反射波，當(dāng)入射角較大時吸收效果較差.（2）衰減邊界條件（Cerjan et al.，1985；Kosloff and Kosloff，1986；Cao and Greenhalgh，1998；Tian et al.，2008）.這類吸收邊界條件在擴邊區(qū)域內(nèi)使波的振幅衰減（通常采用指數(shù)衰減）而被吸收掉，但其具有衰減系數(shù)較難確定、吸收效果較差等缺點.（3）完全匹配層（PML）吸收邊界條件（Bérenger，1994；Zeng et al.，2001；Wang and Tang，2003；Collino and Tsogka，2001；Hu et al.，2007；Gao and Zhang，2008）.這種方法在邊界處加一個匹配層，在匹配層中采用具有衰減效應(yīng)的波動方程來吸收邊界反射，通過改變阻尼因子大小來控制吸收效果.上述三種邊界條件中，PML邊界條件的吸收效果最好，在數(shù)值模擬中得到了廣泛應(yīng)用.

常規(guī) PML 吸收邊界條件 （Bérenger，1994；Collino and Tsogka，2001）是針對一階波動方程設(shè)計的，而且是一種分裂型邊界條件.多年來，PML吸收邊界條件不斷發(fā)展，出現(xiàn)了很多變種.Komatitsch和Tromp（2003）給出了適用于二階波動方程的PML邊界條件.但是由于公式中出現(xiàn)了時間的三階導(dǎo)數(shù)項，計算量及存儲量大大增加.Komatitsch和Martin（2007）指出常規(guī)的PML邊界條件在地震波接近切向入射時失效.因此，一種新的不分裂褶積型PML邊界條件被提出（Komatitsch and Martin，2007；Martin and Komatitsch，2009；Martin et al.，2008；張魯新等，2010；杜啟振等，2011），這種邊界條件在地震波切向入射時也能較好地吸收邊界反射，但計算量更大.

在常規(guī)單程波邊界條件中，內(nèi)部區(qū)域和邊界的波場分別是通過雙程波方程和單程波方程求取的.由于雙程波方程、單程波方程之間的差異，導(dǎo)致內(nèi)部區(qū)域和邊界之間波場存在一個突變，邊界反射是不可避免的.為了克服這一問題，Liu和Sen（2010）提出一種混合吸收邊界條件.該方法在內(nèi)部區(qū)域與邊界之間增加過渡區(qū)域，通過線性加權(quán)單程波波場與雙程波波場實現(xiàn)波場的平滑變化，來減小人工邊界反射，具有計算量小、易于實現(xiàn)和吸收效果好等優(yōu)點.Liu and Sen（2012）將混合吸收邊界條件應(yīng)用于交錯網(wǎng)格彈性波二階位移－應(yīng)力方程正演模擬中.對于二階彈性波方程數(shù)值模擬，混合吸收邊界條件比Komatitsch and Tromp（2003）提出的PML邊界條件更經(jīng)濟可行.但數(shù)值模擬實驗表明，當(dāng)過渡帶點數(shù)較少或計算時間較大時，現(xiàn)有的基于二階位移－應(yīng)力方程的混合吸收邊界條件穩(wěn)定性較差，從而制約了該方法在實際數(shù)值模擬、逆時偏移及全波形反演中的應(yīng)用.與二階位移方程相比，一階速度－應(yīng)力方程具有更高的模擬精度和更好的穩(wěn)定性.在研究地震波傳播規(guī)律中，一階方程更受歡迎（Du et al.，2009；裴正林，2004；王秀明等，2003；楊頂輝，2002；Zhang and Liu，1999；Liu and Sen，2009）.因此，如何將混合吸收邊界條件應(yīng)用到一階彈性波速度－應(yīng)力方程模擬中有待進一步研究.在一階方程數(shù)值模擬中，混合吸收邊界條件與當(dāng)前最受歡迎的分裂PML邊界條件相比，吸收效果、計算量及存儲量有無優(yōu)勢也需要驗證.

本文的主要內(nèi)容包括：簡要回顧了一階彈性波速度－應(yīng)力方程及高階交錯網(wǎng)格有限差分方法；推導(dǎo)了兩種適合于速度－應(yīng)力的單程波方程：二階Higdon單程波方程和一階Higdon單程波方程；提出基于一階彈性波速度－應(yīng)力方程的混合吸收邊界條件實現(xiàn)方法；給出自由邊界的實現(xiàn)方法；分別采用均勻介質(zhì)模型、高縱、橫波速度比模型和Marmousi模型對混合吸收邊界條件和PML邊界條件進行比較；最后對混合吸收邊界條件的穩(wěn)定性做了簡要分析.

2 一階彈性波速度－應(yīng)力方程及高階交錯網(wǎng)格有限差分

二維彈性波速度－應(yīng)力方程如下（Levander，1988）：

式中，（vx，vz）是速度矢量，（τxx，τzz，τxz）是應(yīng)力的三個分量，λ（x，z）和μ（x，z）是拉梅系數(shù)，ρ（x，z）是密度.

采用交錯網(wǎng)格有限差分（Virieux，1986；Luo and Schuster，1990；董良國等，2000）求解方程（1a）－（1e）.2階時間精度和2M階空間精度的差分格式（以式（1a）為例）為

其中，

（裴正林，2004；Liu and Sen，2009）.Δt和h分別為時間和空間采樣間隔，上標(biāo)k為時間節(jié)點索引數(shù)，下標(biāo)（i，j）為空間網(wǎng)格點索引數(shù).

3 Higdon單程波方程

混合吸收邊界條件是通過雙程波波場和單程波波場加權(quán)平均實現(xiàn)的.單程波方程的種類直接決定著混合吸收邊界條件的吸收效果.考慮到計算量和實現(xiàn)復(fù)雜度，本文采用低階Higdon單程波方程.

3.1 二階Higdon單程波方程

二階Higdon位移單程波方程（Higdon，1991）（以左邊界為例）為

β1＝1，β2＝vP／vS，vP和vS分別為縱波和橫波速度.

位移與速度、位移與應(yīng)力的關(guān)系式（Levander，1988）為

聯(lián)立方程（3）和（4），可得二階 Higdon速度－應(yīng)力單程波方程為

為了獲得較高的模擬精度，采用下面的離散格式（Higdon，1991）（以速度的水平分量為例）：

式中的系數(shù)可以通過下面的恒等式得到（Higdon，1991）：

其中，r＝vPΔt／h，b為0.3～0.5之間的常數(shù).

左上角的離散格式可以通過加權(quán)平均左邊界和上邊界的格式得到，形式如下（Liu and Sen，2012）：

3.2 改進的一階Higdon單程波方程

一階Higdon位移單程波方程（Higdon，1991）（以左邊界為例）為

其中，

由上式可知，此方程只對縱波有衰減作用，不能很好地壓制橫波反射.二階Higdon單程波方程（方程（3））中第一項主要衰減縱波，第二項主要衰減橫波.因此，我們構(gòu)造一個新的一階Higdon單程波方程：

聯(lián)立方程（4）和（10），可得改進的一階 Higdon速度－應(yīng)力單程波方程為

為了獲得較高的模擬精度，采用下面的離散格式（Higdon，1991）（以速度的水平分量為例）：

左上角的離散格式仍然可以通過加權(quán)平均左邊界和上邊界的格式得到，形式如下：

4 混合吸收邊界條件

將研究區(qū)域分成邊界（區(qū)域III：A1，B1）、過渡（區(qū)域II：A2，A3，…，AN；B2，B3，…，BN；C1，C2，…，CN；D1，D2，…，DN）和內(nèi)部（區(qū)域I）區(qū)域，如圖1所示.混合吸收邊界條件的實現(xiàn)步驟如下：

（1）對區(qū)域I和區(qū)域II，解雙程波方程（1a）—（1e），得到雙程波波場：

（2）對區(qū)域II和區(qū)域III，解單程波方程（7）—（8）或（12）—（13），得到單程波波場：

（3）對區(qū)域II，將雙程波波場和單程波波場加權(quán)平均得到最終的波場：

ωBi＝ （i－1）／N，ωDi＝ （i－0.5）／N，ωAi＝ωCi＝（i－0.75）／N，（i＝1，2，…，N）分別為速度水平分量、速度垂直分量和應(yīng)力的加權(quán)系數(shù).

混合吸收邊界條件的吸收效果與過渡帶的點數(shù)N密切相關(guān).N＝1時，混合吸收邊界條件退化為單程波邊界條件.本文采用二階Higdon單程波方程和改進的一階Higdon單程波方程構(gòu)建混合吸收邊界條件，分別得到混合二階Higdon吸收邊界條件和混合一階Higdon吸收邊界條件.為了驗證混合吸收邊界條件的有效性，我們將其與數(shù)值模擬中最常用的PML邊界條件進行比較.本文采用分裂PML方法（Bérenger，1994；Collino and Tsogka，2001）.混合吸收邊界條件和PML邊界條件的吸收效果依賴于過渡區(qū)域的點數(shù)（N或LPML），下文的比較都是在LPML＝N＝10條件下進行的.

圖1 混合吸收邊界條件示意圖（Liu and Sen，2012）Fig.1 Illustration of hybrid absorbing boundary conditions（Liu and Sen，2012）

5 自由邊界條件

Lan和Zhang（2011）對數(shù)值模擬中處理自由邊界的方法進行了詳細(xì)研究，主要有直接法和鏡像法兩種.在直接法中，滿足如下方程（Hestholm and Ruud，2000；周竹生等，2007）：

其等價于

由于直接求解方程（16）需要在自由表面附近進行降階處理，可能會導(dǎo)致不穩(wěn)定或模擬精度不夠，因此直接法在實際中很少用.

鏡像法（Levander，1988；Liu and Sen，2012）是另一種處理自由邊界的方法.它通過將自由表面兩側(cè)的應(yīng)力場鏡像反對稱處理，來滿足邊界上應(yīng)力為零的條件.該方法容易實現(xiàn)，但缺少嚴(yán)格的理論支撐.因此，本文采用直接法和鏡像法相結(jié)合的辦法（Robertsson，1996）.

將方程（16b）帶入方程（1c）得

方程（1a）—（1b）中不包括?τxx／?z，因此自由表面上方不需要τxx.則應(yīng)力更新公式如下：

速度更新公式如下：

6 數(shù)值模擬

這里分別采用均勻介質(zhì)，高縱、橫波速度比模型和Marmousi模型來驗證混合吸收邊界條件的有效性.同時，我們也從計算量、存儲量及吸收效果等方面對混合吸收邊界條件和PML邊界條件進行比較.

首先設(shè)計一個均勻介質(zhì)模型.模型大小為2000m×2000m，vP＝4000m／s，vS＝2000m／s，ρ＝2100kg／m3.空間網(wǎng)格大小為10m×10m，時間步長為1ms，記錄長度為2s.震源采用主頻為15Hz的Ricker子波，位于模型中央，施加在速度的水平分量上.

圖2為采用不同吸收邊界條件得到的波場快照（圖（a）、（b）、（c）和（d）中，從上到下依次為x分量和z分量，從左到右依次為370ms、620ms和1000ms時的情況）.由圖2可知，與無吸收邊界條件相比，PML邊界條件、混合一階Higdon吸收邊界條件和混合二階Higdon吸收邊界條件都能較好地吸收邊界反射.為了進一步比較，圖3a給出采用不同吸收邊界條件時點（500m，250m）處的振動圖，其中參考解為模型足夠大時的情況.由圖可知，PML邊界條件、混合一階Higdon吸收邊界條件和混合二階Higdon吸收邊界條件得到的結(jié)果與參考解相近.圖3b是這三種邊界條件模擬結(jié)果與參考解的振動誤差圖.由圖可見，混合二階Higdon吸收邊界條件的吸收效果最好，混合一階Higdon吸收邊界條件次之，PML邊界條件最差.

表1和表2分別為采用不同邊界條件所需要的計算時間和存儲量.本次模擬是在計算機型號為HP p6615cn臺式機（Intel（R）Core（TM）i3處理器、2GB內(nèi)存）上實現(xiàn)的.由表可見，PML邊界條件需要最長的計算時間和最大的存儲空間，混合二階Higdon吸收邊界條件居中，混合一階Higdon吸收邊界條件消耗最少的計算機資源.與常規(guī)的PML方法相比，本文提出的兩種混合Higdon吸收邊界條件在占用較少的計算機資源的情況下能給出更好的吸收結(jié)果，優(yōu)勢明顯.同時考慮計算效率和吸收效果，下文只采用混合一階Higdon吸收邊界條件.

表1 不同吸收邊界條件的計算時間Table 1 The computational time by using different absorbing boundary conditions

圖2 均勻介質(zhì)模型中采用不同邊界條件的波場快照（a）無吸收邊界；（b）PML邊界條件；（c）混合一階Higdon吸收邊界條件；（d）混合二階Higdon吸收邊界條件.圖（a）、（b）、（c）和（d）中，從上到下依次為x分量和z分量，從左到右依次為370ms、620ms和1000ms時的波場快照.Fig.2 Snapshots by using different absorbing boundary conditions for a homogeneous model（a）Non－ABC；（b）PML ABC；（c）Hybrid 1st Higdon ABC；（d）Hybrid 2nd Higdon ABC.x－and z－components are from top to bottom，snapshots at 370ms，620ms and 1000ms from left to right.

表2 不同吸收邊界條件的存儲量Table 2 The memory by using different absorbing boundary conditions

圖3 均勻介質(zhì)模型中采用不同吸收邊界條件時點（500m，250m）處的振動圖（a）和誤差圖（b）圖（a）和（b）中，從上到下依次為x分量和z分量.Fig.3 Records（a）and errors（b）at（500m，250m）by using different absorbing boundary conditions for a homogeneous model（x－and z－components are from top to bottom）

圖4 高縱、橫波速度比模型（Liu and Sen，2012）Fig.4 An elastic model with high vP／vS（Liu and Sen，2012）

在實際彈性波數(shù)值模擬中，縱、橫波速度比太大可能會引起不穩(wěn)定現(xiàn)象（Hestholm，2003；Opr?al and Zahradnik，1999）.因此，采用一個典型的高縱、橫波速度比模型，模型參數(shù)如圖4所示（Opr?al and Zahradnik，1999；Liu and Sen，2012）.空間網(wǎng)格大小為5m×5m，時間步長為0.5ms.震源采用主頻為10Hz的Ricker子波，位于（500m，200m）處，施加在速度的水平分量上.圖5為采用不同邊界條件時得到的波場快照.由圖可見，PML邊界條件和混合一階Higdon吸收邊界條件都能吸收大部分的邊界反射.但在某些部位（箭頭所指的位置），采用混合一階Higdon吸收邊界條件時的邊界反射要更弱.

為了進一步驗證本文提出的基于一階彈性波方程的混合吸收邊界條件的有效性，下面對Marmousi模型（圖6）進行數(shù)值模擬.空間網(wǎng)格大小為10m×10m，時間步長為1ms，記錄長度為4s.震源采用主頻為15Hz的Ricker子波，位于（2500m，0 m）處，施加在速度的水平分量上.本例子中，上邊界采用自由邊界條件.圖7為采用不同邊界條件的地面地震記錄.為了便于比較不同邊界條件的吸收效果，顯示時采用自動增益控制.由圖7b可見，盡管PML邊界條件可以吸收大量邊界反射，但當(dāng)采用增益放大后，可以看到一部分微弱的反射波（箭頭所指的區(qū)域）.而圖7c中，在相同的增益控制下難以看到邊界反射，表明混合吸收邊界條件具有更好的吸收效果.

7 討論

穩(wěn)定性是評價吸收邊界條件好壞的重要指標(biāo).這里將對本文提出的基于速度－應(yīng)力方程的混合吸收邊界條件的穩(wěn)定性進行簡要分析.采用不同時刻內(nèi)部區(qū)域的擾動能量大小來衡量邊界條件的穩(wěn)定性，公式如下：

這里仍然以Marmousi模型為例來對比不同吸收邊界條件的穩(wěn)定性，模型參數(shù)及觀測系統(tǒng)參數(shù)同上.圖8為采用不同邊界條件時內(nèi)部區(qū)域能量大小的變化曲線.由圖可見，當(dāng)時間較小時，混合二階Higdon吸收邊界條件、混合一階Higdon吸收邊界條件及PML邊界條件都是穩(wěn)定的.隨時間增大（大于2000ms）時，混合二階Higdon吸收邊界條件出現(xiàn)了不穩(wěn)定現(xiàn)象，這主要是由于一階雙程波方程和二階單程波方程同時使用導(dǎo)致的.而混合一階Higdon邊界條件仍然是穩(wěn)定的.

圖5 高縱、橫波速度比模型中采用不同邊界條件的波場快照（a）無吸收邊界；（b）PML邊界條件；（c）混合一階 Higdon吸收邊界條件.圖（a）、（b）和（c）中，從上到下依次為x分量和z分量，從左到右依次為2000ms和4000ms的波場快照.Fig.5 Snapshots by using different absorbing boundary conditions for an elastic model with high vP／vS（a）Non－ABC；（b）PML ABC；（c）Hybrid 1st Higdon ABC.x－and z－components are from top to bottom，snapshots at 2000ms and 4000ms from left to right.

為了進一步驗證混合一階Higdon吸收邊界條件的有效性，我們將計算時間增加至8s，并與Liu和Sen的方法（2012）進行了比較，如圖9所示.由于Liu和Sen的方法得到的是位移場，而本文方法求得的是速度場，兩者大小有差異.這里只對比其穩(wěn)定性，結(jié)論如下：

（1）當(dāng)過渡帶點數(shù)較少或計算時間較大時，現(xiàn)有的基于二階位移－應(yīng)力方程的混合吸收邊界條件穩(wěn)定性較差（如：N＝2時，4s以后；N＝5時，7s以后）.

（2）當(dāng)過渡帶點數(shù)較少或計算時間較大時，本文提出的混合一階Higdon吸收邊界條件穩(wěn)定性較好（能量無增大趨勢）.

（3）隨著過渡帶的點數(shù)增加，穩(wěn)定性逐漸變好，但計算量相應(yīng)增加.

圖6 Marmousi vP模型（vS和ρ可以由經(jīng)驗公式求得）Fig.6 The Marmousi vPmodel（vSandρcan be obtained by empirical formulas）

圖7 Marmousi模型中采用不同邊界條件的地面地震記錄（a）無吸收邊界條件；（b）PML邊界條件；（c）混合一階Higdon吸收邊界條件.圖（a）、（b）和（c）中，從左到右依次為x分量和z分量.Fig.7 Seismograms by using different absorbing boundary conditions for the Marmousi model（a）Non－ABC；（b）PML ABC；（c）Hybrid 1st Higdon ABC.x－and z－components are from left to right.

圖8 Marmousi模型中采用不同邊界條件的內(nèi)部區(qū)域能量曲線Fig.8 Variation of energy in the inner area for the Marmousi model by using different absorbing boundary conditions

圖9 Marmousi模型中采用混合一階Higdon吸收邊界條件和Liu和Sen方法的內(nèi)部區(qū)域能量曲線Fig.9 Variation of energy in the inner area for the Marmousi model by using our hybrid 1st Higdon and Liu and Sen′s absorbing boundary conditions

本文提出的混合一階Higdon吸收邊界條件在各種試驗中表現(xiàn)良好，穩(wěn)定可靠，具有較好的應(yīng)用前景.另外，在構(gòu)造混合吸收邊界條件時，單程波方程的種類也很重要，一方面影響邊界條件的吸收效果，另一方面也會對穩(wěn)定性產(chǎn)生影響.我們也嘗試過用Clayton和Engquist（1977）單程波方程構(gòu)造混合吸收邊界條件，但當(dāng)縱、橫波速度比較大時，同樣是不穩(wěn)定的.

8 結(jié)論

本文推導(dǎo)了二階和一階Higdon速度－應(yīng)力單程波方程；提出基于一階彈性波方程的混合吸收邊界條件方法；對多個模型進行了數(shù)值模擬與分析；從計算時間、存儲量和吸收效果等方面，對混合吸收邊界條件和常規(guī)的分裂PML邊界條件進行了比較.主要結(jié)論如下：

（1）在均勻，高縱、橫波速度比及復(fù)雜介質(zhì)的數(shù)值模擬中，混合吸收邊界條件都能夠有效地吸收邊界反射；

（2）與常規(guī)的PML方法相比，混合吸收邊界條件可以在消耗較少的計算時間和存儲量的前提下，給出更好的吸收效果；

（3）與混合二階Higdon吸收邊界條件相比，混合一階Higdon吸收邊界條件具有更好的穩(wěn)定性.

另外，本文方法可以直接推廣到二維或三維復(fù)雜介質(zhì)（黏彈、各向異性）的數(shù)值模擬中，也可以應(yīng)用于涉及到波動方程數(shù)值求解的逆時偏移或全波形反演中.本文提出的吸收邊界條件需要更少的計算機資源，可以提高逆時偏移或全波形反演的計算效率.

Bérenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.J.Comput.Phys.，114（2）：185－200，doi：10.1006／jcph.1994.1159.

Cao S，Greenhalgh S.1998.Attenuating boundary conditions for numerical modeling of acoustic wave propagation.Geophysics，63（1）：231－243，doi：10.1190／1.1444317.

Cerjan C，Kosloff D，Kosloff R，et al.1985.A nonreflecting boundary condition for discrete acoustic and elastic wave equations.Geophysics，50（4）：705－708.doi：10.1190／1.1441945.

Clayton R，Engquist B.1977.Absorbing boundary conditions for acoustic and elastic wave equations.Bull.Seism.Soc.Am.，67（6）：1529－1540.

Collino F，Tsogka C.2001.Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media.Geophysics，66（1）：294－307，doi：10.1190／1.1444908.

Dong L G，Ma Z T，Cao J Z，et al.2000.A staggered－grid highorder difference method of one－order elastic wave equation.ChineseJ.Geophys.（in Chinese），43（3）：411－419.

Du Q Z，Li B，Hou B.2009.Numerical modeling of seismic wavefields in transversely isotropic media with a compact staggered－grid finite difference scheme.Appl.Geophys.，6（1）：42－49，doi：10.1007／s11770－009－0008－z.

Du Q Z，Sun R Y，Zhang Q.2011.Numerical simulation of threecomponent elastic wavefield in 2DTTI media in the condition of the combined boundary.OilGeophys.Prosp.（in Chinese），46（2）：187－195.

Engquist B，Majda A.1977.Absorbing boundary conditions for numerical simulation of waves.Proc.Natl.Acad.Sci.USA，74（5）：1765－1766.

Gao H，Zhang J.2008.Implementation of perfectly matched layers in an arbitrary geometrical boundary for elastic wave modelling.Geophys.J.Int.，174（3）：1029－1036，doi：10.1111／j.1365－246X.2008.03883.x.

Heidari A H，Guddati M N.2006.Highly accurate absorbing boundary conditions for wide－angle wave equations.Geophysics，71（3）：S85－S97，doi：10.1190／1.2192914.

Hestholm S.2003.Elastic wave modeling with free surfaces：Stability of long simulations.Geophysics，68（1）：314－321，doi：10.1190／1.1543217.

Hestholm S O，Ruud B.2000.2Dfinite－difference viscoelastic wave modelling including surface topography.Geophys.Prospect.，48（2）：341－373，doi：10.1046／j.1365－2478.2000.00185.x.

Higdon R L.1991.Absorbing boundary conditions for elastic waves.Geophysics，56（2）：231－241，doi：10.1190／1.1443035.

Hu W Y，Abubakar A，Habashy T M.2007.Application of the nearly perfectly matched layer in acoustic wave modeling.Geophysics，72（5）：SM169－SM175，doi：10.1190／1.2738553.Komatitsch D，Martin R.2007.An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation.Geophysics，72（5）：SM155－SM167，doi：10.1190／1.2757586.

Komatitsch D，Tromp J.2003.A perfectly matched layer absorbing boundary condition for the second－order seismic wave equation.Geophys.J.Int.，154（1）：146－153，doi：10.1046／j.1365－246X.2003.01950.x.

Kosloff R，Kosloff D.1986.Absorbing boundaries for wave propagation problems.J.Comput.Phys.，63（2）：363－376，doi：10.1016／0021－9991（86）90199－3.

Lan H，Zhang Z.2011.Comparative study of the free－surface boundary condition in two－dimensional finite－difference elastic wave field simulation.J.Geophys.Eng.，8（2）：275－286，doi：10.1088／1742－2132／8／2／012.

Levander A R. 1988. Fourth－order finite－difference P－SV seismograms.Geophysics，53（11）：1425－1436，doi：10.1190／1.1442422.

Liu Y，Sen M K.2009.An implicit staggered－grid finite－difference method for seismic modelling.Geophys.J.Int.，179（1）：459－474，DOI：10.1111／j.1365－246X.2009.04305.x.

Liu Y，Sen M K.2010.A hybrid scheme for absorbing edge reflections in numerical modeling of wave propagation.Geophysics，75（2）：A1－A6，doi：10.1190／1.3295447.

Liu Y，Sen M K.2012.A hybrid absorbing boundary condition for elastic staggered－grid modelling.Geophys.Prospect.，60（6）：1114－1132，doi：10.1111／j.1365－2478.2011.01051.x.

Luo Y，Schuster G.1990.Parsimonious staggered grid finitedifference of the wave equation.Geophys.Res.Lett.，17（2）：155－158，doi：10.1029／GL017i002p00155.

Martin R，Komatitsch D.2009.An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation.Geophys.J.Int.，179（1）：333－344，doi：10.1111／j.1365－246X.2009.04278.x.

Martin R，Komatitsch D，Gedney S D.2008.A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seismic wave equation.Comput.Model.Eng.Sci.，37（3）：274－304，doi：10.3970／cmes.2008.037.274.

Opr?al I，Zahradnik J.1999.From unstable to stable seismic modelling by finite－difference method.Phys.Chem.Earth，PartA，24（3）：247－252，doi：10.1016／S1464－1895（99）00026－5.

Pei Z L.2004.Numerical modeling using staggered－grid high－order finite－difference of elastic wave equation on arbitrary relief surface.OilGeophys.Prosp.（in Chinese），39（6）：629－634.

Robertsson J O A.1996.A numerical free－surface condition for elastic／viscoelastic finite－difference modeling in the presence of topography.Geophysics，61（6）：1921－1934，doi：10.1190／1.1444107.

Tian X，Kang I B，Kim G Y，et al.2008.An improvement in the absorbing boundary technique for numerical simulation of elastic wave propagation.J.Geophys.Eng.，5（2）：203－209，doi：10.1088／1742－2132／5／2／007.

Virieux J.1986.P－SV wave propagation in heterogeneous media：Velocity－stress finite－difference method.Geophysics，51（4）：889－901，doi：10.1190／1.1442147.

Wang T，Tang X.2003.Finite－difference modeling of elastic wave propagation：A nonsplitting perfectly matched layer approach.Geophysics，68（5）：1749－1755，doi：10.1190／1.1620648.

Wang X M，Zhang H L，Wang D.2003.Modelling of seismic wave propagation in heterogeneous poroelastic media using a highorder staggered finite－difference method.ChineseJ.Geophys.（in Chinese），46（6）：842－849.

Yang D H.2002.Finite element method of the elastic wave equation and wavefield simulation in two－phase anisotropic media.ChineseJ.Geophys.（in Chinese），45（4）：575－583.

Zeng Y，He J，Liu Q.2001.The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media.Geophysics，66（4）：1258－1266，doi：10.1190／1.1487073.

Zhang L X，F(xiàn)u L Y，Pei Z L.2010.Finite difference modeling of Biot′s poroelastic equations with unsplit convolutional PML and rotated staggered grid.ChineseJ.Geophys.（in Chinese），53（10）：2470－2483，doi：10.3969／j.issn.0001－5733.2010.10.021.

Zhou Z S，Liu X L，Xiong X Y.2007.Finite－difference modelling of Rayleigh surface wave in elastic media.ChineseJ.Geophys.（in Chinese），50（2）：567－573.

Zhang J F，Liu T L.1999.P－SV－wave propagation in heterogeneous media：grid method.Geophys.J.Int.，136（2）：431－438，doi：10.1111／j.1365－246X.1999.tb07129.x.

Zhou H，McMechan G A.2000.Rigorous absorbing boundary conditions for 3－D one－way wave extrapolation.Geophysics，65（2）：638－645，doi：10.1190／1.1444760.

附中文參考文獻(xiàn)

董良國，馬在田，曹景忠等.2000.一階彈性波方程交錯網(wǎng)格高階差分解法.地球物理學(xué)報，43（3）：411－419.

杜啟振，孫瑞艷，張強.2011.組合邊界條件下二維三分量TTI介質(zhì)波場數(shù)值模擬.石油地球物理勘探，46（2）：187－195.

裴正林.2004.任意起伏地表彈性波方程交錯網(wǎng)格高階有限差分法數(shù)值模擬.石油地球物理勘探，39（6）：629－634.

王秀明，張海瀾，王東.2003.利用高階交錯網(wǎng)格有限差分法模擬地震波在非均勻孔隙介質(zhì)中的傳播.地球物理學(xué)報，46（6）：842－849.

楊頂輝.2002.雙相各向異性介質(zhì)中彈性波方程的有限元解法及波場模擬.地球物理學(xué)報，45（4）：575－583.

張魯新，符力耘，裴正林.2010.不分裂卷積完全匹配層與旋轉(zhuǎn)交錯網(wǎng)格有限差分在孔隙彈性介質(zhì)模擬中的應(yīng)用.地球物理學(xué)報，53（10）：2470－2483，doi：10.3969／j.issn.0001－5733.2010.10.021.

周竹生，劉喜亮，熊孝雨.2007.彈性介質(zhì)中瑞雷面波有限差分法正演模擬.地球物理學(xué)報，50（2）：567－573.