黏性聲波方程數(shù)值正演中的匹配Z變換完全匹配層吸收邊界

趙建國(guó)，史瑞其，陳競(jìng)一，潘建國(guó)，王宏斌

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

2中海油研究總院，北京 100027

3Tulsa大學(xué)地球科學(xué)學(xué)院，Oklahoma，74104USA

4中石油勘探開(kāi)發(fā)研究院西北分院，蘭州 730000

1 引言

地震波數(shù)值正演是研究地震波在地下介質(zhì)傳播規(guī)律的重要手段（金勝汶等，1994；劉禮農(nóng)等，2004；謝桂生等，2005）.在使用數(shù)值方法模擬無(wú)限或半無(wú)限大介質(zhì)中的地震波傳播會(huì)因計(jì)算機(jī)內(nèi)存產(chǎn)生截?cái)噙吔?，從而產(chǎn)生虛假數(shù)值反射，吸收邊界條件技術(shù)可以用于消除截?cái)噙吔缧?yīng)產(chǎn)生的虛假反射.經(jīng)過(guò)三十多年的發(fā)展，在地震波數(shù)值領(lǐng)域中主要產(chǎn)生了兩大類吸收邊界技術(shù)：海綿邊界條件（Cerjan etal.，1985；Sochacki etal.，1987）和旁軸近似邊界條件（Clayton and Engquist，1977；Stacey，1988；Higdon，1991；Quarteroni etal.，1988），但是這些吸收邊界技術(shù)在很多情況下表現(xiàn)不佳，比如對(duì)于大角度入射波和低頻能量不能很好的吸收.Bérenger（1994）在Maxwell電磁波方程中提出了完全匹配層（Perfectly Matched Layer，PML），用于消除時(shí)域有限差分法電磁波模擬中的虛假反射.PML對(duì)于較大范圍入射角和頻率成分的波具有高效的吸收作用，理論上對(duì)于入射波具有零反射系數(shù)，目前已經(jīng)成為常用的吸收邊界技術(shù).之后，PML很快被證明在地震波場(chǎng)數(shù)值模擬同樣高效（Chew and Liu，1996；Hasting etal.，1996）.

后來(lái)，PML從數(shù)學(xué)上證明可以同復(fù)拉伸坐標(biāo)的概念聯(lián)系起來(lái)（Chew and Weedon，1994）.但Kuzuoglu和Mittra（1996）指出，使用傳統(tǒng)復(fù)拉伸坐標(biāo)算子的PML對(duì)于隱矢波和極低頻波成分吸收效能會(huì)出現(xiàn)下降，并且在大角度特別是切入射情況下（如震源靠近PML吸收邊界的情況下）會(huì)產(chǎn)生較大程度的虛假反射，因此為解決該問(wèn)題，Kuzuoglu和Mittra（1996）在電磁波模擬中引入了復(fù)頻移拉伸算子用于提高PML的吸收效能.在傳統(tǒng)PML算法中，為在時(shí)間域中求解復(fù)拉伸坐標(biāo)下的波動(dòng)方程，需要對(duì)波場(chǎng)變量進(jìn)行非物理分裂，但經(jīng)典的分裂場(chǎng)形式的PML（Chew and Weedon，1994）不能直接使用復(fù)頻移拉伸算子.在PML實(shí)現(xiàn)形式上，一些非分裂場(chǎng)格式被提出（Cheng and Zhao，2011）.Roden和Gedney（2000）在Maxwell方程模擬中提出了卷積PML來(lái)實(shí)施復(fù)頻移拉伸算子，在PML方程中出現(xiàn)卷積項(xiàng)從而得到非分裂場(chǎng)的形式.之后，使用復(fù)頻移拉伸算子的卷積PML也被引入地震波模擬中（Drossaert and Giannopoulos，2007；Komatitsch and Martin，2007）.

Li和Dai（2008）在電磁波模擬中提出了匹配Z變換PML（Matched Z-transform PML，MZT-PML），同卷積PML是基于Fourier變換理論導(dǎo)出不同，MZT-PML是基于匹配Z變換理論得到非分裂場(chǎng)PML格式，也避免了卷積PML中的片斷常數(shù)近似（Roden and Gedney，2000），其推導(dǎo)過(guò)程也相對(duì)簡(jiǎn)單（Li and Dai，2008）.具有非分裂場(chǎng)格式的MZTPML技術(shù)可以在地震波方程模擬中直接應(yīng)用復(fù)頻移拉伸算子（Shi etal.，2012），本文推導(dǎo)了適用于速度-壓力形式黏聲波方程的MZT-PML吸收邊界控制方程，并對(duì)相應(yīng)的數(shù)值正演結(jié)果進(jìn)行了吸收效果分析.

2 標(biāo)準(zhǔn)線性固體黏聲波方程

在實(shí)際介質(zhì)中，由于地層的吸收衰減作用引起的速度耗散，地震波頻率會(huì)隨著傳播距離的增大而降低，因此近年來(lái)地震勘探研究越對(duì)黏彈性介質(zhì)的關(guān)注不斷提高（Carcione etal.，1988；Liu and Chang，1996）.常用的黏彈性介質(zhì)模型有Kelvin-Voigt模型、Maxwell模型和標(biāo)準(zhǔn)線性固體（Standard Linear Solid，SLS）模型，但Maxwell模型不能描述彈性蠕變，Kelvin-Voigt不能描述應(yīng)力弛豫.SLS模型（Carcione，2007）可以同時(shí)描述上述兩種現(xiàn)象，因此相對(duì)符合實(shí)際黏彈性介質(zhì)中的地震波傳播規(guī)律.本文主要目的在于說(shuō)明PML吸收邊界的應(yīng)用，因此為簡(jiǎn)便起見(jiàn)使用基于單SLS模型的黏聲波方程.在二維Cartesian直角坐標(biāo)系中，速度-壓力格式的

SLS模型黏聲波方程如下（Carcione，2007）：

其中，ω0為震源的中心角頻率.由（1）式可見(jiàn)，當(dāng)弛豫時(shí)間常數(shù)τε＝τσ時(shí)，（1）式即退化為聲波方程.

3 匹配Z變換完全匹配層

將（1）式變換到頻率域，有

其中，i為虛數(shù)單位，ω為角頻率，根據(jù)復(fù)拉伸坐標(biāo)變換（Chew and Weedon，1994）的概念，將空間偏導(dǎo)數(shù)?η（η＝x，z）替換為復(fù)拉伸坐標(biāo)?~η＝S-1η?η.其中Sη為拉伸坐標(biāo)算子，對(duì)PML的吸收特性有重要影響.傳統(tǒng)拉伸坐標(biāo)算子為

其中，實(shí)系數(shù)dη≥0，使波場(chǎng)振幅在PML內(nèi)部呈指數(shù)級(jí)別衰減.使用傳統(tǒng)拉伸算子（4）式的PML存在大角度入射波的吸收效能不足的問(wèn)題，這是因?yàn)閭鹘y(tǒng)拉伸算子只對(duì)行波具有吸收作用，而對(duì)于沿求解區(qū)域和PML區(qū)域界面平行方向傳播的隱矢波不起作用，另外對(duì)于極低頻成分也存在奇異性因而也不能較好地消除極低頻反射（Chew and Liu，1996）.為解決該問(wèn)題，復(fù)頻移拉伸坐標(biāo)算子被提出：

式（5）中通過(guò)引入輔助衰減系數(shù)κη≥1和頻移系數(shù)αη≥0，將復(fù)平面的極點(diǎn)從從實(shí)軸移動(dòng)到負(fù)虛半平面上，以提高對(duì)隱矢波和低頻波的吸收.易見(jiàn)（5）式比（4）式更具有一般性，令κη＝1和αη＝0即為傳統(tǒng)拉伸算子（4）式.在復(fù)拉伸坐標(biāo)下，（3）式改寫(xiě)為因?yàn)槔熳鴺?biāo)算子Sη為頻率ω的函數(shù)，因此將（6）式變換回時(shí)間域會(huì)產(chǎn)生卷積項(xiàng)（Roden and Gedney，2000）.考慮到頻率域的乘法運(yùn)算和時(shí)間域的卷積運(yùn)算都等價(jià)于Z域的乘法運(yùn)算，因此將（6）式變換到Z域可以使推導(dǎo)過(guò)程相對(duì)簡(jiǎn)單.

MZT-PML的主要思想在于將拉伸坐標(biāo)算子（5）式的倒數(shù)分裂為兩部分的乘積，且令iω僅存在于其中一項(xiàng)中：

利用對(duì)應(yīng)關(guān)系iω?s將（6）式變換到S域，有

其中，算符ZT［］·表示匹配Z變換運(yùn)算，Δt表示時(shí)間采樣間隔.為簡(jiǎn)短起見(jiàn)，這里只以式（6a）的MZTPML格式為例進(jìn)行說(shuō)明；另外這里假設(shè)PML吸收邊界由η＝0（η＝x，z）起始.將（6a）由頻率域變換到Z域有：

其中，（1-z-1）／Δt是iω的Z變換，考慮映射關(guān)系iω??t和?t?（1-z-1）／Δt，將（9）式代入（10）式可得

將（12）轉(zhuǎn)化為

式（15）即為關(guān)于聲壓場(chǎng)p連同輔助變量Ψvx，x和Ψvz，z的有限差分時(shí)間更新格式.Ψvx，x和Ψvz，z需要在更新速度場(chǎng)v和壓力場(chǎng)p前計(jì)算，注意到式（15c）中同時(shí)含有不同離散時(shí)刻（n和n+1時(shí)間步）的Ψ，因此Ψn應(yīng)當(dāng)在更新前暫存于一個(gè)臨時(shí)變量中.（6）式中關(guān)于速度場(chǎng)vx、vz和記憶變量R的MZT-PML格式可以同理導(dǎo)出.

4 數(shù)值正演

在數(shù)值模擬時(shí)，空間上采用八階交錯(cuò)網(wǎng)格有限差分格式（Virieux，1986），并采用比Taylor展開(kāi)系數(shù)精度更高的Holberg差分系數(shù)（Holberg，1987）來(lái)計(jì)算所有空間導(dǎo)數(shù).Holberg系數(shù)通過(guò)在給定階數(shù)截?cái)嘤?jì)算得到，降低了群速度產(chǎn)生的誤差從而整體上減小了數(shù)值頻散.例如，x方向空間偏導(dǎo)數(shù)在i+1／2網(wǎng)格點(diǎn)處的差分離散格式為

其中，八階精度Holberg差分系數(shù)C1＝1.231666，C2＝-0.1041182，C3＝0.02063707，C4＝-0.003570998.模型計(jì)算區(qū)域大小為5250m×1500m的長(zhǎng)方形區(qū)域，由701×201個(gè)網(wǎng)格點(diǎn)構(gòu)成，有限差分元胞邊長(zhǎng)為Δx＝Δz＝7.5m.介質(zhì)密度ρ＝2000kg·m-3，縱波速度cp＝3000m／s，品質(zhì)因子Q＝30.時(shí)間上，采用二階顯式差分格式，因此時(shí)間更新步長(zhǎng)Δt由需滿足Courant穩(wěn)定性條件：

取Δt＝1.2ms，模擬持續(xù)時(shí)長(zhǎng)4.8s，共計(jì)4000時(shí)間步.采用點(diǎn)聲壓震源激發(fā)，震源子波為一導(dǎo)數(shù)高斯函數(shù)其中震源中心頻率f0＝20Hz，子波時(shí)延t0＝1.2／f0＝0.06s，激發(fā)點(diǎn)位于（517.5m，180m）處；這里設(shè)置2個(gè)檢波點(diǎn)，檢波點(diǎn)1和檢波點(diǎn)2分別位于（742.5m，1380m）和（5017.5m，142.5m）處.計(jì)算區(qū)域四周由PML吸收邊界包圍，厚度為75m，占10個(gè)有限差分元胞.故震源靠近頂部吸收邊界，檢波點(diǎn)1和2分別靠近底部和頂部吸收邊界（圖1）.一般在PML區(qū)域內(nèi)，要將復(fù)頻移算子（5）式中各參數(shù)設(shè)置為沿PML區(qū)域和求解區(qū)域的交界面外法線方向漸進(jìn)變化：

其中，η為PML內(nèi)一點(diǎn)到與交界面的距離，L為PML的厚度.各衰減參數(shù)取值如下（Li and Matar，2010；Festa and Vilotte，2005）：d0＝-（Pd+1）cplog10（Rc）／2L，κ0＝-（Pκ+1）blog10（Rc）／2L，α0＝πf0；其中Rc＝10-5，Pd＝Pκ＝2，Pα＝1，b＝10Δh，Rc為理論反射系數(shù)，Δh為空間元胞邊長(zhǎng).這里也給出使用復(fù)拉伸坐標(biāo)算子（3）式的傳統(tǒng)分裂場(chǎng)PML吸收邊界條件下的黏聲波方程模擬（Liu and Tao，1997）結(jié)果，用以分析對(duì)比.另外，為檢驗(yàn)吸收邊界效果需要同參考解對(duì)比，本文參考解是在不含吸收邊界的足夠大區(qū)域（模擬時(shí)間內(nèi)邊界反射不會(huì)到達(dá)檢波器）使用相同的數(shù)值離散格式計(jì)算得到，震源和檢波器相對(duì)位置不變.

圖1 計(jì)算區(qū)域示意圖圖中黑色內(nèi)框線表示PML與求解區(qū)域的交界，圓圈表示震源位置，三角表示檢波點(diǎn)位置Fig.1 Diagrammatic figure of the computational model The inner black lines denote the interface between PML and solution domain，the circle denotes the source location and the triangles denote receiver locations

圖2 傳統(tǒng)PML邊界條件下（a）0.48s，（b）0.96s和（c）1.68s的黏聲波聲壓場(chǎng)p的波場(chǎng)快照（振幅已加10倍增益）Fig.2 Snapshots of the pressure field pof viscoacoustic wave at（a）0.48s，（b）0.96sand（c）1.68s in the classical PML condition，and the amplitudes has been scaled by a factor of 10

圖3 MZT-PML邊界條件下（a）0.48s，（b）0.96s和（c）1.68s的黏聲波聲壓場(chǎng)p的波場(chǎng)快照（振幅已加10倍增益）Fig.3 Snapshots of the pressure field pof viscoacoustic wave at（a）0.48s，（b）0.96sand（c）1.68s in the MZT-PML condition，and the amplitudes has been scaled by a factor of 10

該算例中，由于震源靠近頂部吸收邊界，形成了大角度入射PML的情況.從波場(chǎng)快照?qǐng)D2可見(jiàn)，傳統(tǒng)PML對(duì)大角度入射波吸收不完全，產(chǎn)生了低頻虛假反射，且虛假反射能量隨著入射角度和傳播距離的增大而增強(qiáng).從1.68s的波場(chǎng)可見(jiàn)，在遠(yuǎn)偏移距低頻反射尤為明顯.而圖3中MZT-PML邊界由于使用復(fù)頻移拉伸算子，對(duì)大角度入射波吸收很好，沒(méi)有出現(xiàn)圖1中的低頻反射；在圖3中1.68s波場(chǎng)中，即使在遠(yuǎn)偏移距也沒(méi)有可見(jiàn)的虛假反射.

從圖4a可見(jiàn)，在檢波點(diǎn)1處無(wú)論是傳統(tǒng)PML還是MZT-PML的計(jì)算結(jié)果都和參考解吻合良好，說(shuō)明在小角度入射情況下，兩種吸收邊界都可以高效吸收外行波場(chǎng)能量.而從圖4b中檢波點(diǎn)2的地震記錄可見(jiàn)，PML情況下的地震記錄上產(chǎn)生了明顯的虛假反射抖動(dòng)，與參考解吻合不佳，這是因?yàn)闄z波器2記錄到了大角度入射波產(chǎn)生的低頻虛假反射，這些虛假反射在遠(yuǎn)偏移距更強(qiáng).而在檢波器2處，MZT-PML情況下的記錄與參考解吻合良好，說(shuō)明MZT-PML相對(duì)于傳統(tǒng)PML，對(duì)切入射波的吸收有明顯改善.這里也給出完全彈性介質(zhì)情況下（令弛豫時(shí)間常數(shù)τε＝τσ）聲波方程在MZT-PML邊界條件下的模擬結(jié)果：圖5中，對(duì)比聲波方程和黏聲波方程的地震記錄，由于介質(zhì)的黏滯性，導(dǎo)致聲波振幅衰減強(qiáng)烈（如圖5b）.但MZT-PML邊界條件下，無(wú)論是聲波方程還是黏聲波方程的計(jì)算結(jié)果都與各自參考解吻合良好.

另外，長(zhǎng)時(shí)間數(shù)值正演中的穩(wěn)定性也是吸收邊界效能的一個(gè)重要指標(biāo)（Komatitsch and Martin，2007；Chen and Bording，2010；Chen，2011），本文將波場(chǎng)傳播時(shí)間延長(zhǎng)至120s（共105時(shí)間步），并記錄求解區(qū)域Ω中（不包括PML吸收邊界）的波場(chǎng)總能量E為

圖4 傳統(tǒng)PML和MZT-PML邊界條件下得到的（a）檢波點(diǎn)1和（b）檢波點(diǎn)2的黏聲波方程地震記錄Fig.4 Seismograms calculated in classical PML and MZT-PML conditions at（a）receiver 1and（b）receiver 2

圖5 MZT-PML邊界條件下得到的（a）檢波點(diǎn)1和（b）檢波點(diǎn)2的聲波、黏聲波方程地震記錄Fig.5 Seismograms of acoustic and viscoacoustic wave calculated in MZT-PML condition at（a）receiver 1and（b）receiver 2

圖6 吸收邊界條件下黏聲波場(chǎng)的總能量衰減（a）0～6s（b）0～120s，縱坐標(biāo)軸為半對(duì)數(shù)坐標(biāo)顯示Fig.6 Total energy decay of the viscoacoustic wave field on a semi-logarithmic scale for y-axis in the figure（a）0～6s（b）0～120s

從圖6中展示了彈性和黏彈性介質(zhì)中兩種吸收邊界條件下的波場(chǎng)總能量衰減情況.由圖6a可見(jiàn)，當(dāng)數(shù)值正演開(kāi)始時(shí)（約0.1s左右），震源將能量注入介質(zhì)使得總能量急劇上升，之后由于黏彈性介質(zhì)本身的衰減，使得黏聲波情況下的總量明顯低于彈性聲波.由于PML吸收邊界的作用，也導(dǎo)致兩種介質(zhì)中的波場(chǎng)能量開(kāi)始逐漸下降.在約1.7s時(shí)能量出現(xiàn)一個(gè)陡降，因?yàn)檫@時(shí)直達(dá)波完全離開(kāi)求解區(qū)域，理論上講該時(shí)刻之后求解區(qū)域中的能量全部是虛假能量，這時(shí)的能量衰減情況可以較為準(zhǔn)確地反映吸收邊界的效能，從圖6a可見(jiàn)MZT-PML情況下波場(chǎng)剩余能量低于PML情況（紅藍(lán)色曲線之間，以及黑色實(shí)虛線之間的對(duì)比）.由圖6b可見(jiàn)，在數(shù)值正演后期（6s之后），能量仍舊持續(xù)下降，說(shuō)明MZTPML在彈性聲波、黏聲波方程模擬中具有長(zhǎng)時(shí)間數(shù)值穩(wěn)定性.在黏聲波介質(zhì)中，MZT-PML情況下總能量持續(xù)下降至10-30J數(shù)量級(jí)，低于PML情況下的10-13J數(shù)量級(jí).

5 討論和結(jié)論

本文基于匹配Z變換的思想建立了用于黏聲波方程正演的非分裂MZT-PML，這種非分裂PML技術(shù)可以直接使用復(fù)頻移拉伸算子.在有限差分黏聲波方程的數(shù)值正演實(shí)驗(yàn)中，同傳統(tǒng)分裂場(chǎng)PML邊界條件的應(yīng)用效果進(jìn)行了比較.數(shù)值實(shí)驗(yàn)結(jié)果說(shuō)明，同傳統(tǒng)PML相比，當(dāng)震源靠近吸收邊界時(shí)，MZT-PML可以明顯改善對(duì)大角度入射波的吸收.另外，MZT-PML在完全彈性介質(zhì)和黏彈性介質(zhì)情況下均表現(xiàn)良好.從105時(shí)間步的能量衰減計(jì)算可見(jiàn)，MZT-PML邊界條件下波場(chǎng)總能量在數(shù)值正演后期總體上呈持續(xù)下降趨勢(shì)，說(shuō)明MZT-PML在黏聲波方程模擬中具有長(zhǎng)時(shí)間穩(wěn)定性.該方法具有比傳統(tǒng)PML更好的吸收衰減效能.

Bérenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.J.Comput.Phys.，114（2）：185-200.

Carcione J M，Kosloff D，Kosloff R.1988.Wave propagation simulation in a linear viscoelastic medium.Geophys.J.Int.，95（3）：597-611.

Carcione J M.2007.Wave Fields in Real Media：Theory and Numerical Simulation of Wave Propagation in Anisotropic，Anelastic，Porous and Electromagnetic Media.2nd edn.Amsterdam，The Netherlands：Elsevier Science.

Cerjan C，Kosloff D，Kosloff R，etal.1985.A nonreflecting boundary condition for discrete acoustic and elastic wave equations.Geophysics，50（4）：705-708.

Chen J Y，Bording R P.2010.Application of the nearly perfectly matched layer to the propagation of low-frequency acoustic waves.J.Geophys.Eng.，7（3）：277-283.

Chen J Y.2011.Application of the nearly perfectly matched layer for seismic wave propagation in 2Dhomogeneous isotropic media.Geophysical Prospecting，59（4）：662-672.

Cheng J Y，Zhao J G.2011.Application of the nearly perfectly matched layer to seismic-wave propagation modeling in elastic anisotropic media.B.Seismol.Soc.Am.，101（6）：2866-2871.

Chew W C，Liu Q H.1996.Perfectly matched layers for elastodynamics：A new absorbing boundary condition.J.Comput.Acoust.，4（4），341-359.

Chew W C，Weedon W H.1994.A 3Dperfectly matched medium from modified Maxwell′s equations with stretched coordinates.Microw.Opt.Tech.Lett.，7（13）：599-604.

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

Drossaert F H，Giannopoulos A.2007.Complex frequency shifted convolution PML for FDTD modelling of elastic waves.Wave Motion，44（7-8）：593-604.

Engquist B，Majda A.1977.Absorbing boundary conditions for the numerical simulation of waves.Mathematics of Computation，31（139）：629-651.

Festa G，Vilotte J P.2005.The Newmark scheme as velocity–stress time-staggering：an efficient PML implementation for spectral element simulations of elastodynamics.Geophys.J.Int.，161（3）：789-812.

Hasting F D，Schneider L B，Broschat S L.1996.Application of the perfectly matched layer（PML）absorbing boundary condition to elastic wave propagation.J.Acoust.Soc.Am.，100（5）：3061-3069.

Higdon R L.1991.Absorbing boundary conditions for elastic waves.Geophysics，56（2）：231-241.

Holberg O.1987.Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena.Geophysical Prospecting，35（6）：629-655.

Jin S W，Chen B Y，Ma Z T.1994.Three-dimensional wave equation forward modeling by the finite-difference methods.Acta Geophysica Sinica（in Chinese），37（6）：804-810.

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.

Kuzuoglu M，Mittra R.1996.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave Guided Wave Letters，6（12）：447-449.

Li J X，Dai J F.2008.Modified Z-transform-based FDTD algorithm for the anisotropic perfectly matched layer.International Journal of Numerical Modelling：Electronic Networks，Devices and Fields，21（5）：279-286.

Li Y F，Matar O B.2010.Convolutional perfectly matched layer for elastic second-order wave equation.J.Acoust.Soc.Am.，127（3）：1318-1327.

Liu L N，Cui F L，Zhang J F.2004.Seismic modeling with one-way wave equation in 3-D complex structures.Chinese J.Geophys.（in Chinese），47（3）：514-520.

Liu Q H，Chang C.1996.Compressional head waves in attenuative formations：forward modeling and inversion.Geophysics，61（6）：1908-1920.

Liu Q H，Tao J P.1997.The perfectly matched layer for acoustic waves in absorptive media.J.Acoust.Soc.Am.，102（4）：2072-2082.

Quarteroni A，Tagliani A，Zampieri E.1988.Generalized Galerkin approximations of elastic waves with absorbing boundary conditions.Comput.Method.Appl.M.，163（1-4）：323-341.

Roden J A，Gedney S D.2000.Convolution PML（CPML）：An efficient FDTD implementation of the CFS-PML for arbitrary media.Microw.Opt.Tech.Lett.，27（5）：334-339.

Shi RQ，Wang SX，Zhao JG.2012.An unsplit complex-frequencyshifted PML based on matched Z-transform for FDTD modelling of seismic wave equations.J.Geophys.Eng.，9（2）：218-229.Sochacki J，Kubichek R，George J，etal.1987.Absorbing boundary conditions and surface waves.Geophysics，52（1）：60-71.

Stacey R.1988.Improved transparent boundary formulations for the elastic-wave equation.B.Seismol.Soc.Am.，78（6）：2089-2097.

Virieux J.1986.P-SV wave propagation in heterogeneous media：velocity-stress finite-difference method.Geophysics，51（4）：889-901.

Xie G S，Liu H，Li Y M，etal.2005.Seismic modeling by reflection／transmission operator and one-way wave equation under the condition of fluctuating reflectors.Chinese J.Geophys.（in Chinese），48（5）：1172-1178.

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

金勝汶，陳必遠(yuǎn)，馬在田.1994.三維波動(dòng)方程有限差分正演方法.地球物理學(xué)報(bào)，37（6）：804-810.

劉禮農(nóng)，崔風(fēng)林，張劍鋒.2004.三維復(fù)雜構(gòu)造中地震波模擬的單程波方法.地球物理學(xué)報(bào)，47（3）：514-520.

謝桂生，劉洪，李幼銘等.2005.界面起伏條件下反射／透射算子+單程波方程的地震波模擬方法.地球物理學(xué)報(bào)，48（5）：1172-1178.