基于MCPML邊界條件的頻率域聲波方程高精度模擬

董士琦 韓立國 胡 勇 羅玉欽

(吉林大學地球探測科學與技術(shù)學院，吉林長春130026)

1 引言

頻率域波場數(shù)值模擬以數(shù)據(jù)選擇靈活、適合多炮同時模擬、沒有時間累積誤差以及計算效率高等優(yōu)勢，近年來越發(fā)受到人們的重視。Lysmer等[1]最早使用有限元方法實現(xiàn)了頻率域正演模擬；Marfurt等[2，3]對比了時間域和頻率域不同正演算法的精度和計算效率；Pratt等[4]給出了聲波方程的頻率域5點差分格式，但是數(shù)值頻散嚴重；Jo等[5]提出了頻域聲波方程的9點差分格式，將45°坐標系和0°坐標系結(jié)合引入加權(quán)參數(shù)來構(gòu)造Laplace算子的差分格式，通過構(gòu)造阻抗矩陣并加以適當?shù)恼鹪春瘮?shù)求解Helmhotz方程實現(xiàn)頻率域的波場模擬，將一個波長所需最小網(wǎng)格數(shù)減小到3.3；Stekl等[6]將9點差分格式應(yīng)用于彈性波數(shù)值模擬，并研究了彈性波在黏彈性介質(zhì)中的傳播； Shin等[7]構(gòu)造了聲波方程25點差分格式，但是仍然只有二階差分精度，且計算效率不高。近年來，業(yè)界廣泛研究了如何提高頻率域正演精度和效率[8-21]。曹紅書等[22]提出了由0°坐標系和45°坐標系四階9點差分格式結(jié)合而成的四階17點差分格式，提高了正演精度。目前，頻率域正演普遍應(yīng)用9點差分格式，面對復(fù)雜介質(zhì)構(gòu)造時精度較低，且存在較大的數(shù)值頻散。

在吸收邊界研究方面，Bérenger[23]首先在電磁學中提出了完全匹配層(PML)吸收邊界；Chew等[24]實現(xiàn)了將復(fù)拉伸坐標系引入PML邊界條件。大量研究證明PML邊界條件較傳統(tǒng)的吸收邊界條件(如指數(shù)衰減吸收邊界、旁軸近似吸收邊界等)具有更好的衰減效果[25-28]。Kuzuoglu等[29]對坐標進行復(fù)頻移變換提出了CPML邊界條件，改善了大入射角的吸收效果和低頻成分產(chǎn)生的奇異值現(xiàn)象；Meza-Fajardo等[30]提出了在衰減區(qū)域中加入互相正交的衰減因子實現(xiàn)了多軸完全匹配層(MPML)，加強了各向異性介質(zhì)中波場能量衰減的穩(wěn)定性。近年來，業(yè)界還研究了PML邊界條件的改進[31-41]，田坤等[42]提出了非分裂的多軸卷積完全匹配層(MCPML)邊界條件，同時具有CPML和MPML邊界條件的優(yōu)點，在改進邊界衰減效果的同時保證了計算效率。目前，頻率域正演計算主要應(yīng)用傳統(tǒng)PML邊界條件，可以發(fā)揮頻率域正演多炮同時模擬的優(yōu)勢，但對不同炮點位置的衰減效果不同，尤其是對大入射角的波場衰減效果有限。

為進一步提高正演模擬精度、改進邊界吸收效果，本文推導了頻域17點差分格式并應(yīng)用了MCPML邊界條件。在頻率域應(yīng)用MCPML邊界條件只需要對常規(guī)PML邊界公式進行簡單變換，不需要逆Fourier變換后進行卷積計算，在改善衰減效果的同時可以提高計算效率。在頻率域?qū)Ρ攘薓CPML與PML和CPML的衰減效果，并應(yīng)用Marmousi進行試算，驗證了方法的有效性。

2 頻率—空間域聲波方程17點差分格式的構(gòu)造

二維各向同性頻率—空間域聲波方程由時間域聲波方程做Fourier變換可得

(1)

A(x,z,ω)P(x,z,ω)=-S(x,z,ω)

(2)

式中A為阻抗矩陣，是一大型復(fù)數(shù)稀疏矩陣(圖1)。通過計算阻抗矩陣并加載震源函數(shù)，即可得到該頻率下的單頻波場。

圖1 17點差分阻抗矩陣結(jié)構(gòu)示意圖 波場計算區(qū)域網(wǎng)格邊長為6，阻抗矩陣大小為36×36

為了構(gòu)造離散化的波動方程，需要對Laplace算子做差分近似。Laplace算子的0°坐標系和45°坐標系的四階9點差分格式分別為

(3)

(4)

式中：Pm,n表示[xm,zn]=[x0+(m-1)Δx,z0+(n-1)Δz]處的波場值，本文取Δx=Δz=Δh。圖2給出了0°坐標系、45°坐標系以及混合坐標系Laplace算子的四階差分示意圖。

將式(3)和式(4)兩個Laplace算子加權(quán)平均，可得

(5)

式中a是加權(quán)系數(shù)。因此， 頻域聲波方程可以表示為

(6)

圖2 Laplace算子的四階差分示意圖 (a)0°坐標系四階9點差分格式； (b)45°坐標系 四階9點差分格式； (c)四階17點差分格式

波場中任一點的波場值可以由它本身以及周圍16個點的加權(quán)平均表示

P(xm,zn)=bPm,n+c(Pm-1,n+Pm+1,n+Pm,n-1+

Pm,n+1)+d(Pm-1,n-1+Pm+1,n-1+Pm-1,n+1+

Pm+1,n+1)+e(Pm-2,n+Pm+2,n+Pm,n-2+Pm,n+2)+

f(Pm-2,n-2+Pm+2,n-2+Pm-2,n+2+Pm+2,n+2)

(7)

式中：b、c、d、e、f是加權(quán)系數(shù)， 且有b+4c+4d+4e+4f=1。將上式展開并整理，即可得到頻率域波動方程的17點差分格式

T1Pm,n+T2(Pm-1,n+Pm+1,n+Pm,n-1+Pm,n+1)+

T3(Pm-2,n+Pm+2,n+Pm,n-2+Pm,n+2)+

T4(Pm-1,n-1+Pm+1,n-1+Pm-1,n+1+Pm+1,n+1)+

T5(Pm-2,n-2+Pm+2,n-2+Pm-2,n+2+Pm+2,n+2)

=-S(ω)

(8)

圖3為傳統(tǒng)的二階混合網(wǎng)格9點差分格式與四階17點差分格式頻散曲線對比，其中縱坐標vph/v表示數(shù)值速度與真實速度之比，比值越接近1表示數(shù)值頻散越小。由圖可見：傳統(tǒng)的二階混合網(wǎng)格9點差分格式最小波長所需要的網(wǎng)格點數(shù)H為3.3

圖3 兩種差分格式不同傳播角度頻散曲線對比 (a)二階混合網(wǎng)格9點差分格式； (b)四階17點差分格式

(圖3a)；在H=2.56時二階9點差分格式的速度誤差達到±14%(圖3a)，而四階17點差分格式速度誤差在±1%以內(nèi)(圖3b)。

不同差分格式阻抗矩陣的構(gòu)造不同，求解Helmhotz方程的耗時也不同。17點差分阻抗矩陣對角線條數(shù)以及分解稀疏矩陣求解Helmhotz方程的計算復(fù)雜程度都大于9點差分格式，因而會帶來更多的計算耗時。表1比較了網(wǎng)格數(shù)為400×400的兩種差分格式求解Helmhotz方程的耗時。

表1 兩種差分格式的計算耗時

3 頻率域PML邊界條件

PML邊界條件[23]是一種對邊界反射衰減效果非常好的處理方法，基本思想是在波場模擬計算的區(qū)域外邊界加入若干個層做吸收，核心為足夠的邊界厚度以及衰減因子。

對笛卡爾坐標進行了復(fù)拉伸坐標變換(以x方向為例)

(10)

式中：ω是角頻率；G(x)是x方向的衰減因子；sx是x方向復(fù)拉伸系數(shù)。對于二階偏導，有

(11)

則含有PML邊界條件的二維頻域聲波方程變?yōu)?/p>

k2P(x,z,ω)=-S(ω)

(12)

衰減因子G有兩種定義方式，分別為

(13)

和

(14)

上兩式中：m是常數(shù)，一般取2或3；vmax是介質(zhì)中的最大縱波速度；R是邊界反射系數(shù)，一般取值范圍為10-4～10-6；L為PML層厚度；l為PML區(qū)域內(nèi)部計算點到PML邊界的距離；f0為子波主頻；a0是經(jīng)驗值，取1.79。

圖4為17點差分格式兩種衰減因子G1與G2衰減效果對比，其計算區(qū)域為2400m×2400m，震源頻率為20Hz，均勻介質(zhì)的縱波速度為1000m/s，PML層厚度為60m。由圖4可知，使用衰減因子G1的頻域波場有邊界反射干擾(圖4a)，使用衰減因子G2的波場圖像幾乎沒有邊界反射(圖4b)，因此本文在進行頻域正演時衰減因子采用G2。

圖4 采用不同衰減因子的頻域波場 (a)G1； (b)G2

4 非分裂的MCPML邊界條件

4.1 非分裂的CPML邊界條件

為了改進大入射角時常規(guī)PML方法衰減效果，Drossaert等[43]引入復(fù)頻移因子α和收縮因子χ，提出了基于復(fù)頻移變換的非分裂卷積完全匹配層(CPML)邊界條件，用遞歸積分的方法避免了直接計算卷積，從而提高了計算效率。以x方向為例，式(9b)變?yōu)?/p>

(15)

式中：χx≥1；αx≥0；復(fù)頻移因子α可以避免當頻率過小時衰減因子出現(xiàn)奇異值，收縮因子可以使大角度的入射波在進入PML區(qū)域后向該區(qū)域的衰減法向彎曲，從而提高衰減效果[44]。通常設(shè)

(16)

通常取αmax=f0，χmax=1。

4.2 分裂的MPML邊界條件

二維的一階速度—應(yīng)力方程為

(17)

式中：W是由速度分量vx、vz和應(yīng)力分量σxx、σzz、σxz構(gòu)成的列向量；B和C是系數(shù)矩陣。傳統(tǒng)分裂的PML將速度分量分為正交的兩個方向，以vx為例，將其分成x方向分量和z方向分量之和，即

(18)

則在計算區(qū)域內(nèi)滿足方程

(19)

相應(yīng)的匹配層內(nèi)的方程為

(20)

MPML邊界條件在傳統(tǒng)分裂的PML邊界條件基礎(chǔ)上加入了正交的衰減因子。以垂直于x軸的衰減邊界中的vx為例，有

(21)

(22)

p(z/x)是不同方向衰減因子比例系數(shù)，也叫穩(wěn)定性因子，通過調(diào)整該參數(shù)可使邊界衰減效果穩(wěn)定。穩(wěn)定性因子可通過經(jīng)驗或試算獲得。

4.3 MCPML邊界條件

MCPML邊界條件將分裂的MPML和非分裂的CPML邊界條件相結(jié)合，同時保留復(fù)頻移坐標變換以及在衰減區(qū)域內(nèi)加入相互正交的衰減因子。以垂直于x方向的衰減區(qū)域為例，該方向上復(fù)拉伸系數(shù)與式(15)相同，在z方向上的復(fù)拉伸系數(shù)為

(23)

對于與衰減區(qū)域平行的衰減因子取值范圍不包括四個邊角。四個邊角的加載方式與PML邊界一樣，只做簡單的疊加。由于在各個衰減邊界區(qū)域內(nèi)加入了相互正交的衰減因子，在模擬大角度入射邊界時，衰減效果會顯著提升。通過改變穩(wěn)定性因子的取值，MCPML可轉(zhuǎn)化為前三種邊界條件，因此MCPML邊界條件可以看成常規(guī)PML方法的拓展。

根據(jù)式(11)可以導出含有MCPML邊界條件的17點差分格式為

R1Pm,n+R2(Pm-1,n+Pm+1,n)+R3(Pm,n-1+Pm,n+1)+

R4(Pm-2,n+Pm+2,n)+R5(Pm,n-2+Pm,n+2)+

R6(Pm-1,n-1+Pm+1,n+1+Pm+1,n-1+Pm-1,n+1)+

R7(Pm-2,n-2+Pm+2,n+2+Pm+2,n-2+Pm-2,n+2)

=-S(ω)

(24)

應(yīng)用式(24)實現(xiàn)了17點有限差分格式與MCPML邊界條件的結(jié)合(圖5)，可以通過編程實現(xiàn)高精度頻率域正演模擬

。

5 頻率域PML邊界條件對比

為驗證MCPML邊界條件在大入射角情況下的有效性，對各向同性均勻介質(zhì)進行17點差分格式數(shù)值模擬。將震源放在計算區(qū)域的左上角模擬大入射角的情況。定義：模型網(wǎng)格數(shù)為400×400，網(wǎng)格間距為6m，子波主頻為20Hz，聲波速度為1000m/s，PML厚度為60m。由于震源位于計算區(qū)域的左上角，震源激發(fā)后波將以大角度入射計算區(qū)域的上邊界和左邊界。圖6a顯示PML邊界在上邊界和左邊界吸收效果不好，反射波與直達波疊加在一起，加強了波形振幅；在下方的匹配層中出現(xiàn)了較強的邊界反射，嚴重干擾了正常波場。圖6b顯示CPML邊界稍微改善了對大入射角的吸收效果；圖6c顯示MCPML邊界條件，由于在吸收邊界區(qū)域加入了相互正交的衰減因子，對于大入射角時衰減效果有了顯著提升。在計算效率方面，應(yīng)用三種邊界條件計算耗時只有毫秒級別的差別(表2)；對于多炮正演模擬的計算，由于炮點位置不同產(chǎn)生的邊界反射可以忽略。

圖6 大入射角情況下三種邊界條件衰減效果對比

(a)PML邊界； (b)CPML邊界； (c)MCPML邊界(p(x/z)=0.3，p(z/x)=0.3)。左為頻率域波場； 右為z=36m的波場強度

表2 三種邊界條件的計算耗時對比

6 Marmousi模型試算

Marmousi模型(圖7a)網(wǎng)格數(shù)為576×183，網(wǎng)格間距為16m，吸收邊界厚度L=20，雷克子波主頻為20Hz，震源置于模型左上角。

對比圖7b與圖7c可見，應(yīng)用CPML邊界條件的正演結(jié)果上邊界衰減效果有限，反射波與直達波疊加形成干擾(圖7b 紅圈處)；而應(yīng)用MCPML邊界條件的正演結(jié)果邊界衰減完全，可提高模擬精度。試驗頻率間隔為0.01Hz，共計算了4000個頻率切片。將頻率域正演結(jié)果進行逆Fourier變換為時域波場。對比圖7d與圖7e可見，在加載MCPML邊界條件時，頻率域9點差分格式正演波場快照轉(zhuǎn)換有明顯的數(shù)值頻散現(xiàn)象(圖7d 紅圈處)，而頻率域17點差分格式正演結(jié)果波形清晰，數(shù)值頻散微弱。圖7f為17點差分格式MCPML邊界條件正演的6.25s波場快照，波前和反射波清晰，幾乎無數(shù)值頻散和邊界反射干擾，證明了本文方法的優(yōu)越性。

圖7 Marmousi模型試算結(jié)果 (a)Marmousi速度模型； (b)17點差分格式CPML邊界單頻波場； (c)17點差分格式MCPML邊界單頻波場； (d)9點差分 格式2.5s波場快照； (e)17點差分格式2.5s波場快照； (f)17點差分格式MCPML邊界條件正演結(jié)果6.25s波場快照

7 結(jié)論

本文推導了頻率域17點有限差分正演格式，與傳統(tǒng)的9點差分格式相比，精度高、數(shù)值頻散弱，更適用于復(fù)雜地質(zhì)情況下的地震波模擬。應(yīng)用MCPML邊界條件，改善了邊界對大角度入射波的衰減效果。在實際應(yīng)用方面：用本文方法進行多炮正演時，震源在不同位置同時激發(fā)產(chǎn)生的邊界反射波振幅很小，可忽略不計；在頻率域應(yīng)用MCPML邊界條件只需要對復(fù)拉伸系數(shù)加以簡單變換，不需要對波場分裂計算，易于編程實現(xiàn)，與時間域相比，不需要做逆Fourier變換和時域卷積計算，計算效率更高。

本文方法的缺點是17點差分格式構(gòu)造阻抗矩陣并求解Helmhotz方程的計算量非常大，計算效率較9點差分格式低，今后需要從算法和計算機硬件兩個方面提高本文方法的計算效率，并進一步改善邊界條件，使波場在邊界處能夠衰減得更徹底。

[1] Lysmer J，Drake L A.A finite element method for seismology.Methods in Computational Physics：Advances in Research & Applications，1972，11：181-216.

[2] Marfurt K J.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics，1984，49(4)：533-549.

[3] Marfurt K J，Shin C S.The future of iterative modeling in geophysical exploration∥Supercomputers in Seismic Exploration(Handbook of Geophysical Exploration：Seismic Exploration，Volume 21).Pergamon Press，New York，1989，203-228.

[4] Pratt R G.Frequency-domain elastic wave modeling by finite difference：A tool for crosshole seismic imaging.Geophysics，1990，55(5)：626-632.

[5] Jo C H，Shin C S,Suh J H.An optimal 9-point,finite difference,frequency-space,2-D scalar wave extrapolator.Geophysics，1996，61(2)：529-537.

[6] Stekl I，Pratt R G.Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators.Geophysics，1998，63(5)：1779-1794.

[7] Shin C，Sohn H.A frequency-space 2-D scalar wave extrapolator using extended 25-point finite-difference operator.Geophysics，1998，63(1)：289-296

[8] Brossier R，Virieux J，Operto S.Parsimonious finite-volume frequency-domain method for 2-D P-SV wave modelling.Geophysical Journal International，2008，175(2)：541-559.

[9] OpertoS，VirieuxJ，RibodettiAetal.Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic(TTI) media.Geophysics，2009，74(5)：T75-T95.

[10] Dupuy B，De Barros L，Garambois S et al.Wave propa-gation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous Galerkin method.Geophysics，2011，76(4)：N13-N28.

[11] Hustedt B，Operto S，Virieux J.Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling.Geophysical Journal International，2004，157(3)：1269-1296.

[12] Min D J，Shin C，Kwon B D et al.Improved frequency-domain elastic wave modeling using weighted-averaging difference operators.Geophysics，2000，65(3)：884-895.

[13] 馬召貴，王尚旭，宋建勇.頻率域波動方程正演中的多網(wǎng)格迭代算法.石油地球物理勘探，2010，45(1)：1-5. Ma Zhaogui，Wang Shangxu，Song Jianyong.Multi-grid iterative algorithm in frequency domain wave equation modeling.OGP，2010，45(1)：1-5.

[14] 韓利，韓立國，李翔等.二階聲波方程頻域PML邊界條件及頻域變網(wǎng)格步長并行計算.吉林大學學報(地球科學版)，2011，41(4)：1226-1232. Han Li，Han Liguo，Li Xiang et al.PML boundary conditions for second-order acoustic wave equations and variable grid parallel computation in frequency-domain modeling.Journal of Jilin University(Earth Science Edition)，2011，41(4)：1226-1232.

[15] 高正輝，孫建國，孫章慶等．基于完全匹配層構(gòu)建新方法的2.5維聲波近似方程數(shù)值模擬．石油地球物理勘探，2016，51(6)：1128-1133. Gao Zhenghui，Sun Jianguo，Sun Zhangqing et al.2.5D acoustic approximate equation numerical simulation with a new construction method of the perfectly matched layer.OGP，2016，51(6)：1128-1133.

[16] 蔡志成，顧漢明，成景旺等．大規(guī)模變網(wǎng)格三維地震正演MPI并行策略與實現(xiàn)．石油地球物理勘探，2017，52(3)：468-476. Cai Zhicheng，Gu Hanming，Cheng Jingwang et al.Parallel large-scale variable grid 3D seismic wave modeling.OGP，2017，52(3)：468-476.

[17] 殷文，印興耀，吳國忱等.高精度頻率域彈性波方程有限差分方法及波場模擬.地球物理學報，2006，49(2)：561-568. Yin Wen，Yin Xingyao，Wu Guochen et al.The me-thod of finite difference of high precision elastic wave equations in the frequency domain and wave field simulation.Chinese Journal of Geophysics，2006，49(2)：561-568.

[18] 劉璐，劉洪，劉紅偉.優(yōu)化15點頻率—空間域有限差分正演模擬.地球物理學報，2013，56(2)：644-652. Liu Lu，Liu Hong，Liu Hongwei.Optimal 15-point finite difference forward modeling in frequency-space domain.Chinese Journal of Geophysics，2013，56(2)：644-652.

[19] 梁文全，王彥飛，楊長春．聲波方程數(shù)值模擬矩形網(wǎng)格有限差分系數(shù)確定法.石油地球物理勘探，2017，52(1)：56-62. Liang Wenquan，Wang Yanfei，Yang Changchun.Acoustic wave equation modeling with rectangle grid finite difference operator and its linear time space domain solution.OGP，2017，52(1)：56-62.

[20] 周聰，劉江平，羅銀河等．二維頻率域全波場有限差分數(shù)值模擬方法．石油地球物理勘探，2014，49(2)：278-287. Zhou Cong，Liu Jiangping，Luo Yinhe et al.2D full-wavefield modeling in frequency domain using finite-difference.OGP，2014，49(2)：278-287.

[21] 馬曉娜，李志遠，谷丙洛等.2D聲波頻率域數(shù)值模擬中幾種有限差分方法的對比分析.地球物理學進展，2015，30(2)：878-888. Ma Xiaona，Li Zhiyuan，Gu Bingluo et al.Comparisons and analysis of several optimization finite-differencing schemes in 2D acoustic frequency-domain numerical modeling.Progress in Geophysics，2015，30(2)：878-888.

[22] 曹紅書，陳景波.聲波方程頻率域高精度正演的17點格式及數(shù)值實現(xiàn).地球物理學報，2012，55(10)：3440-3449. Cao Hongshu，Chen Jingbo.A 17-point scheme and its numerical implementation for high-accuracy modeling of frequency-domain acoustic equation.Chinese Journal of Geophysics，2012，55(10)：3440-3449.

[23] Bérenger J P.A perfectly matched layer for the absorp-tion of electromagnetic waves.Journal of Computational Physics，1994，114(9)：185-200.

[24] Chew W C，Weedon W H.A 3-D perfectly matched medium from modified Maxwell's equations with stretched coordinates.Microwave & Optical Techno-logy Letters，1994，7(13)：599-604.

[25] Engquist B，Majda A.Absorbing boundary conditions for numerical simulation of waves.Mathematics of Computation， 1977，31(139)：629-651.

[26] Shin C.Sponge boundary condition for frequency-domain modeling.Geophysics，1995，60(6)：1870-1874.

[27] Higdon R L.Absorbing boundary conditions for elastic wave.Geophysics，1991，56(2)：231-241.

[28] Liao Z P，Wang H L，Yang B P et al.A transmitting boundary for transient wave analysis.Science China(Series A)，1984，27(10)：1063-1076.

[29] Kuzuoglu M，Mittra R.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave & Guided Wave Letters，2002，6(12)：447-449.

[30] Meza-Fajardo K C，Papageorgiou A S.A non-convolutional split-field，perfectly matched layer for wave propagation in isotropic and anisotropic elastic media：stability analysis.Bulletin of the Seismological Society of America，2008，98(4)：1811-1836.

[31] Hu W，Abubakar A，Habashy T.Application of the nearly perfectly matched layer in acoustic wave mo-deling.Geophysics，2006，72(5)：SM169-SM175.

[32] Gao H，Zhang J.Implementation of perfectly matched layers in an arbitrary geometrical boundary for elastic wave modeling.Geophysical Journal International，2008，174(3)：1029-1036.

[33] NilssonS，PeterssonN，Sj?greenBetal.Stablediffe-rence approximations for the elastic wave equation in second order formulation.Siam Journal on Numerical Analysis，2007，45(5)：1902-1936.

[34] Martin R，Komatitsch D，Bruthiaux E et al.Convolution and nonconvolution perfectly matched layer techniques optimized at grazing incidence for high-order wave propagation modelling.Geophysical Research Abstracts，2010，12：EGU2010-4268-2.

[35] Meza-Fajardo K C，Papageorgiou A S.On the stability of a non-convolutional perfectly matched layer for isotropic elastic media.Soil Dynamics & Earthquake Engineering，2010，30(3)：68-81.

[36] Zeng C，Xia J H，Miller R D et al.Application of the multiaxial perfectly matched layer (MPML) to near-surface seismic modeling with Rayleigh waves.Geophysics，2011，76(3)：T43-T52.

[37] 熊章強，毛承英.聲波數(shù)值模擬中改進的非分裂式

PML邊界條件.石油地球物理勘探，2011，46(1)：35-39. Xiong Zhangqiang，Mao Chengying.Improved unsplit PML boundary conditions in acoustic wave numerical simulation.OGP，2011，46(1)：35-39.

[38] 張顯文，韓立國，黃玲等.基于遞歸積分的復(fù)頻移PML彈性波方程交錯網(wǎng)格高階差分.地球物理學報，2009，52(7)：1800-1807. Zhang Xianwen，Han Liguo，Huang Ling et al.A staggered-grid high-order difference method of complex frequency-shifted PML based on recursive intergration for elastic wave equation.Chinese Journal of Geophysics，2009，52(7)：1800-1807.

[39] 劉洋.波動方程時空域有限差分數(shù)值解及吸收邊界條件研究進展.石油地球物理勘探，2014，49(1)：35-46. Liu Yang.Research progress on time-space domain finite difference numerical solution and absorption boundary conditions of wave equation.OGP，2014，49(1)：35-46.

[40] 王守東.聲波方程完全匹配層吸收邊界.石油地球物理勘探，2003，38(1)：31-34. Wang Shoudong.Acoustic wave equation PML absorbing boundary.OGP，2003，38(1)：31-34.

[41] 李佩笑，林偉軍，張秀梅等.常規(guī)分裂和非分裂完全匹配層吸收邊界比較研究.聲學學報，2015，40(1)：44-53. Li Peixiao，Lin Weijun，Zhang Xiumei et al.Comparisons for regular spitting and non-splitting perfectly matched layer absorbing boundary conditions.Acta Acustica，2015，40(1)：44-53

[42] 田坤，黃建平，李振春等.多軸卷積完全匹配層吸收邊界條件.石油地球物理勘探，2014，49(1)：143-152. Tian Kun，Huang Jianping，Li Zhenchun et al.Multi-axial convolution perfectly matched layer absorption boundary condition.OGP，2014，49(1)：143-152.

[43] Drossaert F H，Giannopoulos A.A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves.Geophysics，2007，72(2)：T9-T17.

[44] Zhang W，Shen Y.Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling.Geophysics，2010，75(4)：T141-T154.