CFS－PML在探地雷達(dá)FDTD法模擬中應(yīng)用

石 明，馮德山，吳 奇

（1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083；2.貴州省有色金屬和核工業(yè)地質(zhì)勘查局 物化探總隊(duì)，貴州 都勻 558004；3.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410083）

0 引 言

探地雷達(dá)（ground penetrating radar，GPR）數(shù)值模擬是地球物理重要研究?jī)?nèi)容，對(duì)建立的模型開(kāi)展數(shù)值模擬能有效指導(dǎo)地質(zhì)解釋.應(yīng)用時(shí)域有限差分（finite difference time domain，F(xiàn)DTD）法進(jìn)行GPR 數(shù)值模擬時(shí)，因計(jì)算機(jī)內(nèi)存空間的有限，總是在某處把網(wǎng)格空間截?cái)?，使之成為有限區(qū)域.這樣，在網(wǎng)格空間截?cái)嗵幘蜁?huì)產(chǎn)生強(qiáng)烈的非物理的電磁反射干擾波.為了在截?cái)噙吔缣幉灰鹛摷俜瓷?，吸收邊界條件對(duì)模擬精度起著關(guān)鍵作用.目前在截?cái)噙吔缣幊Ｍㄟ^(guò)加載完全匹配層（perfectly matched layer，PML）來(lái) 改 善 吸 收 效果，但無(wú)論是PML 還是單軸各向異性完全匹配層（uniaxial anisotropic perfectly matched layer，UPML）對(duì)隱失波、低頻波、大角度入射的掠角波都不能很好地吸收.為此，本文引入了復(fù)頻移卷積完全匹配層（complex frequency shifted convolution perfectly matched layer，CFS－PML）以改善截?cái)噙吔缣帉?duì)隱失波、低頻波、掠角波的吸收效果.

早期提出的 Taflove－Brodwin 邊界［1］、Engquist－Majda 邊界［2］、Mur 邊界［3］、Liao 邊界［4］，通常在模擬區(qū)域的外邊界仍有0.5% ～5.0%的反射系數(shù)，它們逐漸被PML 邊界條件所取 代.PML邊界條件最初是由Berenger 提出的［5－6］，它在FDTD網(wǎng)格外邊界加載一種非物理的吸收媒質(zhì)，使進(jìn)入PML 的透射波迅速衰減，將邊界條件的研究向前推進(jìn)了一大步；但Berenger PML的理論體系是非Maxwell方程的，物理機(jī)制模糊，同時(shí)，其電磁場(chǎng)分量分裂技術(shù)增加了計(jì)算內(nèi)存與數(shù)值實(shí)現(xiàn)的難度［7］，而且只對(duì)行波有吸收效果，對(duì)空域中衰減的隱失波、低頻波以及小角度入射的掠角波無(wú)吸收效果.為了改善邊界條件的吸收效果，Sacks等［8］和Gedney［9］提出了UPML.UPML不需要對(duì)電磁場(chǎng)分裂，計(jì)算更簡(jiǎn)潔、易編程實(shí)現(xiàn)，因其在吸收系數(shù)中引入了線性的吸收因子，UPML邊界對(duì)隱失波、低頻波等干擾波有一定的吸收效果.詹應(yīng)林等［10］將UPML 應(yīng)用于二維FDTD探地雷達(dá)正演中；Feng等［11－12］將UPML應(yīng)用到交替方向隱式時(shí)域有限差分（ADI－FDTD）法中實(shí)現(xiàn)了探地雷達(dá)的二、三維數(shù)值模擬；李靜等［13］開(kāi)展了三維UPML邊界條件的高階FDTD探地雷達(dá)正演.盡管UPML對(duì)于隱失波、低頻波等干擾波有一定的吸收效果，但仍然不是特別完美.

考慮到復(fù)頻移拉伸函數(shù)的PML 技術(shù)可以吸收低頻信號(hào)在邊界表面的反射，同時(shí)對(duì)隱失波等大角度掠射波也有很好的吸收效果，Kuzuoglu等［14］在復(fù)頻率偏移完全匹配層的變量Sk中引入低頻分量吸收參數(shù)αk，將復(fù)平面的極點(diǎn)從實(shí)軸移動(dòng)到負(fù)虛半平面上，加強(qiáng)了對(duì)低頻波、隱失波的吸收；Roden等［15］將復(fù)頻移拉伸函數(shù)加入到電磁波模擬的PML 方程，有效壓制了信號(hào)的假反射；Wang等［16］將CFS－PML應(yīng)用到ADI－FDTD 法中；Drossaert等［17－18］和Komatitsh 等［19］提出了基于遞歸積分的波場(chǎng)非分裂的復(fù)頻移PML（recursive integration CFS－PML，簡(jiǎn)稱(chēng)CFSRIPML）；Zhang 等［20］和Martin 等［21］將 基 于 輔助微分方程的CFS－PML 應(yīng)用到彈性波數(shù)值模擬中，該方法易于擴(kuò)展到高階差分形式；張顯文等［22］將遞歸復(fù)頻移PML應(yīng)用到交錯(cuò)網(wǎng)格高階差分彈性波波動(dòng)方程正演中；張魯新等［23］將不分裂的完全匹配層（CPML）結(jié)合旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分技術(shù)應(yīng)用到孔隙彈性介質(zhì)模擬中；Giannopoulos［24］和Li等［25］將基于CFS－RIPML 吸收邊界技術(shù)應(yīng)用在電磁波波動(dòng)方程的有限差分正演中.相比常規(guī)PML、UPML 吸收邊界，CFS－PML 對(duì)于隱失波、低頻波等具有更好的吸收效果，所以CFS－PML有非常廣闊的發(fā)展前景.

1 零極點(diǎn)匹配Z變換與逆Laplace變換

根據(jù)Kuzuoglu等［14］的研究，定義復(fù)坐標(biāo)拉伸坐標(biāo)系中的算子與拉伸變量分別為

式中：κk的引入是為了改善PML 對(duì)表面波的吸收特性；σk為PML內(nèi)k方向電導(dǎo)率參數(shù)；αk則是為了改善PML 對(duì)低頻分量的吸收特性.σk和αk為大于0的正實(shí)數(shù)，κk為大于等于1的正實(shí)數(shù).根據(jù)式（2）中的變量（k＝x，y，z）的Z變換，利用關(guān)系式j(luò)ω→s，將式（2）變換到S域，則有

其中

然后利用零極點(diǎn)匹配Z變換（matchedZtransform）關(guān)系：

將式（3）變換到Z域

其中

求式（3）的倒數(shù)為

其中

定義倒數(shù)的逆Laplace變換為則式（8）的拉普拉斯逆變換為［26］

式中：δ（t）是單位沖激函數(shù)，u（t）是單位階躍函數(shù)，ζk（t）為

2 二維TM 波的CFS－PML 邊界條件推導(dǎo)

由電磁波傳播理論可知［27］，GPR遵循的Maxwell方程的6個(gè)耦合偏微分方程在二維情況下變?yōu)?個(gè)偏微分方程，共有Ez、Hx、Hy3個(gè)分量，稱(chēng)之為T(mén)M 波：

其中E為電場(chǎng)強(qiáng)度，V／m；H為磁場(chǎng)強(qiáng)度，A／m；ε為介電常數(shù)，F(xiàn)／m；μ為磁導(dǎo)率，H／m；σe為電導(dǎo)率，S／m；σm為等效磁阻率，Ω／m.

將頻率域方程（13）～（15）進(jìn)行反Fourier變換，并將式（10）代入，可得到時(shí)域方程

式（16）～（18）中涉及時(shí)域卷積的計(jì)算，直接計(jì)算效率較低，為了有效地計(jì)算方程中的卷積，根據(jù)離散脈沖ζk（t）的定義，可得到

其中

將式（19）和（20）代入式（16）～（18）可得到：En＋1z（i，j）＝CA（m）·Enz（i，j）＋CCx（m）·

其中ψEzx、ψEzy、ψHxy和ψHyx為引入的輔助變量，它們可采用遞歸卷積方法計(jì)算：

符號(hào)標(biāo)記CA（m）～CRy（m）中的電磁參數(shù)標(biāo)號(hào)m的取值與式（21）～（23）中右端電場(chǎng)或磁場(chǎng)分量的空間位置相同，其表達(dá)式如下：

在二維情況下，模型內(nèi)部κx＝κy＝1，σx＝σy＝0，αx＝αy＝0，Sx＝Sy＝1.在CFS－PML 邊界區(qū)域，這些參數(shù)被賦予不同的值，設(shè)σx和σy大于0，允許傳播的波被吸收，另一方面，使κx和κy大于1，有利于CFS－PML區(qū)域吸收衰減波.αx和αy大于0，能提高對(duì)低頻波的吸收.為使CFS－PML 達(dá)到最優(yōu)吸收效果，設(shè)置吸收參數(shù)從CFS－PML 外邊界到內(nèi)邊界逐步減小到0.對(duì)κx和κy有

其中d（0≤d≤δ）為CFS－PML 內(nèi)的計(jì)算點(diǎn)到CFS－PML區(qū)域內(nèi)邊界的距離；δ為CFS－PML 區(qū)域的厚度，計(jì)算過(guò)程中常取8～10個(gè)空間網(wǎng)格步長(zhǎng).m被稱(chēng)為CFS－PML 的指數(shù)參數(shù)，κkmax為參數(shù)κk的最大值.對(duì)于參數(shù)σx和σy有

因?yàn)樽鴺?biāo)擴(kuò)展變量是一維函數(shù)，這里κk和σk只沿k（x或y）方向變化.為了衰減波的有效吸收，在常見(jiàn)的程序計(jì)算中常設(shè)置：

3 CFS－PML吸收邊界條件效果分析

3.1 均勻介質(zhì)中邊界條件對(duì)比

建立圖1所示3.0 m×3.0 m 的模型，模型中均勻介質(zhì)相對(duì)介電常數(shù)為10，電導(dǎo)率為0.002 S／m，在模型的中心加入500 MHz的脈沖Ricker子波，應(yīng)用不帶邊界條件及加入不同邊界條件的FDTD 法對(duì)該模型進(jìn)行正演，以說(shuō)明不同邊界條件的優(yōu)劣，正演中取網(wǎng)格空間步長(zhǎng)為dx＝dy＝0.005m，dt＝11.79×10－12s，時(shí)窗長(zhǎng)度為3.5×10－8s.

圖1 均勻介質(zhì)模型示意圖Fig.1 The sketch map of homogeneous medium model

圖2為不加載任何邊界條件正演所得的15、22、30ns雷達(dá)波場(chǎng)快照；由圖2（a）的15ns快照中可見(jiàn)，雷達(dá)波以圓形方式向外傳播；圖2（b）雷達(dá)波在截?cái)噙吔绲?條邊、4個(gè)角點(diǎn)引起了較強(qiáng)的虛假反射波；圖2（c）截?cái)喾瓷洳ㄓ忠栽ㄇ跋喾捶较騻鞑?，從左邊的電?chǎng)數(shù)量級(jí)上可以看出，該截?cái)喾瓷洳ǖ哪芰枯^強(qiáng)，取值在－200～200V／m變化，幾乎與一次波前的能量在一個(gè)數(shù)量級(jí).

圖3（a）為采用二階Mur吸收邊界處理截?cái)噙吔绾蟮?0ns波場(chǎng)快照，4條邊的人為截?cái)鄰?qiáng)反射波得到了很好的吸收，僅4個(gè)角點(diǎn)處的邊界反射波仍然可見(jiàn)；但從左邊的電場(chǎng)取值范圍0.5～1.0V／m 不難看出，Mur邊界條件對(duì)邊界處反射波起了較好的吸收作用，但是效果仍然有待提高.圖3（b）、（c）分別為加入U(xiǎn)PML 吸收邊界、CFSPML吸收邊界30ns的波場(chǎng)快照，內(nèi)部區(qū)域仍清晰可見(jiàn)4條截?cái)噙吔缣幍姆瓷洌?個(gè)角點(diǎn)處的反射消失不見(jiàn)，再看左邊的電場(chǎng)取值范圍，顯然4條邊界的反射波能量是非常弱的，其能量在－0.1～0.1 V／m 變化，說(shuō)明這兩種邊界條件較二階Mur吸收邊界條件有數(shù)量級(jí)的提高.

3.2 UPML與CFS－PML吸收效果對(duì)比驗(yàn)證

圖2 無(wú)邊界條件時(shí)的不同時(shí)刻波場(chǎng)快照三維顯示Fig.2 The 3Dvisualization map of wave field snapshots at different time without boundary condition

圖3 加載不同邊界條件的30ns波場(chǎng)快照三維顯示Fig.3 The 3Dvisualization map of wave field snapshots at 30ns with different boundary conditions

為了對(duì)比UPML 與CFS－PML 邊界條件對(duì)隱失波、低頻波以及小角度入射的掠角波等干擾波的吸收效果，設(shè)置圖4所示雙層介質(zhì)模型，模型尺寸為4.0m×1.0m，模型中上層介質(zhì)1的介電常數(shù)為6.0，電導(dǎo)率為0.020S／m，介質(zhì)2的介電常數(shù)為3.0，電導(dǎo)率為0.001S／m，在S（x＝2.0 m，y＝0.1m）位置加入一個(gè)頻率為900 MHz的Ricker子波脈沖，分別應(yīng)用8個(gè)網(wǎng)格的UPML與CFS－PML邊界條件的FDTD 法對(duì)該模型進(jìn)行正演，正演中取網(wǎng)格空間步長(zhǎng)為dx＝dy＝0.005m，dt＝11.79×10－12s，指數(shù)參數(shù)m＝4，CFS－PML的κkmax＝5，時(shí)窗長(zhǎng)度為4.0×10－8s.

圖4 雙層介質(zhì)模型Fig.4 Model of double－layer medium

圖5、6為8、10、12、15、20、25ns的電場(chǎng)分量Ez的波場(chǎng)快照，其中圖5為UPML 吸收邊界，圖6為CFS－PML吸收邊界.首先取最大值的0.5%設(shè)為閾值，將小于該閾值的波場(chǎng)值置零，為了加強(qiáng)小振幅的顯示，以0.3為冪進(jìn)行歸一化.圖5、6為經(jīng)過(guò)上述處理，歸一化后的波場(chǎng)快照.從圖5可以看出，對(duì)于狹長(zhǎng)模型或井筒模型，以掠射角入射的雷達(dá)波在模型邊界上產(chǎn)生虛假的反射，作為有效波的干擾回到模型區(qū)域內(nèi)，UPML上邊界上產(chǎn)生一系列的呈有規(guī)律分布的虛假反射波.而Kuzuoglu等［14］提出的CFS－PML 拉伸函數(shù)將奇異性 由實(shí)軸移到復(fù)平面的負(fù)虛半軸，改善了吸收邊界性能.從圖6可以看出，基于CFS－PML 吸收邊界法對(duì)入射到邊界層內(nèi)波的能量進(jìn)行了充分的吸收，即使對(duì)于以掠射角入射到邊界層內(nèi)的狹長(zhǎng)模型吸收效果也較UPML理想.

3.3 吸收邊界的全局反射誤差分析

圖5 UPML邊界條件下的電場(chǎng)分量Ez 不同時(shí)刻波場(chǎng)快照Fig.5 The wave field snapshots of electric field component Ezat different time under UPML boundary condition

圖6 CFS－PML邊界條件下的電場(chǎng)分量Ez 不同時(shí)刻波場(chǎng)快照Fig.6 The wave field snapshots of electric field component Ezat different time under CFS－PML boundary condition

圖7 P1處UPML和CFS－PML與解析解對(duì)比圖Fig.7 The comparison map of UPML，CFS－PML and analytical solution at P1position

全局反射誤差是反映吸收邊界吸收性能的關(guān)鍵參數(shù)，它能反映一種邊界條件在整個(gè)區(qū)域模擬計(jì)算中的誤差精度.為了說(shuō)明CFS－PML 與UPML的吸收效果，仍以圖4所示的模型為例，觀測(cè)點(diǎn)P1的坐標(biāo)為（0.08m，0.08m）.圖7（a）、（b）分 別 為 加 載UPML 和CFS－PML 吸 收 邊 界FDTD 法正演所得信號(hào)與精確信號(hào)的擬合曲線，其中的精確信號(hào)為將模型區(qū)域擴(kuò)大至原模型4倍而獲得的信號(hào).從圖7中可以看出，觀測(cè)點(diǎn)P1處CFS－PML與解析解擬合更好，而UPML 在邊界處因掠射波、低頻波不能被較好吸收，在19ns左右與解析解擬合得并不太好.通過(guò)模擬計(jì)算對(duì)比驗(yàn)證CFS－PML、UPML 兩種吸收邊界在觀測(cè)點(diǎn)與激勵(lì)源位置角度最大的P1 點(diǎn)的全局反射誤差，以更直觀地充分反映吸收邊界的整體吸收效果.全局反射誤差的獲得是通過(guò)在設(shè)置激勵(lì)源與觀測(cè)點(diǎn)位置關(guān)系的基礎(chǔ)上，分別采用正常模型和擴(kuò)大4倍模型獲得觀測(cè)點(diǎn)處的實(shí)測(cè)接收信號(hào)Es和參考模型信號(hào)Eref，并求取參考模型信號(hào)中振幅最大的值Erefmax，通過(guò)以下計(jì)算公式獲得全局反射誤差：

通過(guò)設(shè)置相同的觀測(cè)參數(shù)及吸收邊界層數(shù)獲得了兩種不同吸收邊界的反射誤差對(duì)比，如圖8所示，CFS－PML吸收邊界的反射誤差相比UPML 吸收邊界其整體誤差更小，說(shuō)明CFS－PML 吸收邊界較UPML有一定幅度的提高與改善.

圖8 P1處兩種吸收邊界反射誤差對(duì)比圖Fig.8 The comparison map of reflection error under UPML and CFS－PML at P1position

4 結(jié) 論

（1）推導(dǎo)了CFS－PML邊界條件的GPR 正演FDTD 差分公式，應(yīng)用均勻介質(zhì)中GPR 波場(chǎng)邊界處的傳播快照，對(duì)比沒(méi)有加載邊界條件、加載Mur、UPML、CFS－PML 邊界處的波場(chǎng)反射波的強(qiáng)弱，實(shí)例結(jié)果表明UPML 與CFS－PML 能較好地處理GPR 截?cái)噙吔?，具有較好的吸收效果.

（2）設(shè)置雙層介質(zhì)狹長(zhǎng)模型對(duì)比了UPML與CFS－PML的吸收效果，波場(chǎng)快照與全局反射誤差都 說(shuō) 明，CFS－PML 較UPML 對(duì) 隱 失 波、低 頻波、掠角波有更高效的吸收，非常適合于寬頻帶脈沖GPR 正演中FDTD 法的吸收邊界.由于CFSPML邊界條件的優(yōu)良吸收特性，其能廣泛地應(yīng)用于電磁場(chǎng)仿真、無(wú)線傳播、電磁兼容、微波技術(shù)等領(lǐng)域的FDTD 計(jì)算中.

［1］ Taflove A，Brodwin M E.Numerical solution of steady－state EM scattering problems using the timedependent Maxwell′s equations ［J］.IEEE Transactions on Microwave Theory and Techniques，1975，23（8）：623－630.

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

［3］ Mur G.Absorbing boundary conditions for the finite－difference approximation of the time－domain electromagnetic－field equations ［J］.IEEE Transactions on Electromagnetic Compatibility，1981，EMC－23（4）：377－382.

［4］ LIAO Zhen－peng，Wong H L，YANG Bai－po，et al.A transmitting boundary for transient wave analyses ［J］.Science in China，Ser.A，1984，27（10）：1063－1076.

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

［6］ Berenger J P.Three－dimensional perfectly matched layer for the absorption of electromagnetic waves［J］.Journal of Computational Physics，1996，127（2）：363－379.

［7］ 葛德彪，閆玉波.電磁波時(shí)域有限差分方法［M］.西安：西安電子科技大學(xué)出版社，2005：67－81.GE De－biao，YAN Yu－bo.Finite－Difference Time－Domain Method for Electromagnetic Waves ［M］.Xi′an：Xidian University Press，2005：67－81.（in Chinese）

［8］ Sacks Z S，Kingsland D M，Lee R，etal.Perfectly matched anisotropic absorber for use as an absorbing boundary condition ［J］.IEEE Transactions on Antennas and Propagation，1995，43（12）：1460－1463.

［9］ Gedney S.An anisotropic perfectly matched layer absorbing media for the truncation of FDTD lattice［J］.IEEE Transactions on Antennas and Propagation，1996，44（12）：1630－1639.

［10］ 詹應(yīng)林，昌彥君，曹中林.基于UPML 吸收邊界的探地雷達(dá)數(shù)值模擬研究［J］.資源環(huán)境與工程，2008，22（2）：235－238.ZHAN Ying－lin，CHANG Yan－jun，CAO Zhonglin.Modeling of ground penetrating radar based on UPML boundary conditions ［J］.Resources Environment and Engineering，2008，22（2）：235－238.（in Chinese）

［11］ FENG De－shan，DAI Qian－wei.GPR numerical simulation of full wave field based on UPML boundary condition of ADI－FDTD［J］.NDT ＆ E International，2011，44（6）：495－504.

［12］ 馮德山，謝 源.基于單軸各向異性完全匹配層邊界條件的ADI－FDTD三維GPR 全波場(chǎng)正演［J］.中南大學(xué)學(xué)報(bào)（自然科學(xué)版），2011，42（8）：2363－2371.FENG De－shan，XIE Yuan.Three dimensional GPR numerical simulation of full wave field based on UPML boundary condition of ADI－FDTD ［J］.Journal of Central South University（Science and Technology），2011，42（8）：2363－2371.（in Chinese）

［13］ 李 靜，曾昭發(fā)，黃 玲，等.三維探地雷達(dá)數(shù)值模擬中UPML 邊界研究［J］.物探化探計(jì)算技術(shù)，2010，32（1）：6－12.LI Jing，ZENG Zhao－fa，HUANG Ling，etal.Study of UPML boundary for three dimensional GPR simulation ［J］.Computing Techniques for Geophysical and Geochemical Exploration，2010，32（1）：6－12.（in Chinese）

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

［15］ Roden J A，Gedney S.Convolutional PML（CPML）：an efficient FDTD implementation of the CFS－PML for arbitrary media［J］.Microwave and Optical Technology Letters，2000，27（5）：334－339.

［16］ WANG Lin－nian，LIANG Chang－h(huán)ong.A new implementation of CFS－PML for ADI－FDTD method［J］.Microwave and Optical Technology Letters，2006，48（10）：1924－1928.

［17］ Drossaert F H，Giannopoulos A.A nonsplit complex frequency－shifted PML based on recursive integration for FDTD modeling of elastic waves［J］.Geophysics，2007，72（2）：9－17.

［18］ Drossaert F H，Giannopoulos A.Complex frequency shifted convolution PML for FDTD modeling of elastic waves［J］.Wave Motion，2007，44（7－8）：593－604.

［19］ Komatitsch D，Martin R.An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation ［J］.Geophysics，2007，72（5）：155－167.

［20］ ZHANG Wei，SHEN Yang.Unsplit complex frequency－shifted PML implementation using auxiliary differential equations for seismic wave modeling［J］.Geophysics，2010，75（4）：141－154.

［21］ Martin R，Komatitsch D，Gedney S，etal.A highorder time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations（ADE－PML）［J］.Computer Modeling in Engineering and Sciences，2010，56（1）：17－40.

［22］ 張顯文，韓立國(guó)，黃 玲，等.基于遞歸積分的復(fù)頻移PML彈性波方程交錯(cuò)網(wǎng)格高階差分法［J］.地球物理學(xué)報(bào)，2009，52（7）：1800－1807.ZHANG Xian－wen，HAN Li－guo，HUANG Ling，etal.A staggered－grid high－order difference method of complex frequency－shifted PML based on recursive integration for elastic wave equation［J］.Chinese Journal of Geophysics，2009，52（7）：1800－1807.（in Chinese）

［23］ 張魯新，符力耘，裴正林.不分裂卷積完全匹配層與旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分在孔隙彈性介質(zhì)模擬中的應(yīng)用［J］.地 球 物 理 學(xué) 報(bào)，2010，53（10）：2470－2483.ZHANG Lu－xin，F(xiàn)U Li－yun，PEI Zheng－lin.Finite difference modeling of Biot′s poroelastic equations with unsplit convolutional PML and rotated staggered grid［J］.Chinese Journal of Geophysics，2010，53（10）：2470－2483.（in Chinese）

［24］ Giannopoulos A.An improved new implementation of complex frequency shifted PML for the FDTD method ［J］.IEEE Transactions on Antennas and Propagation，2008，56（9）：2995－3000.

［25］ LI Jing，ZENG Zhao－fa，HUANG Ling，etal.GPR simulation based on complex frequency shifted recursive integration PML boundary of 3D high order FDTD［J］.Computers and Geosciences，2012，49：121－130.

［26］ 李建雄.時(shí)域有限差分法中完全匹配層的實(shí)現(xiàn)算法研究［D］.天津：天津大學(xué)，2007.LI Jian－xiong.Research on algorithms for implementing perfectly matched layers in the finite difference time domain method［D］.Tianjin：Tianjin University，2007.（in Chinese）

［27］ Yee K S.Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media［J］.IEEE Transactions on Antennas and Propagation，1996，14（3）：302－307.