CFS－PML在探地雷達FDTD法模擬中應用

石 明，馮德山，吳 奇

（1.中南大學 地球科學與信息物理學院，湖南 長沙 410083；2.貴州省有色金屬和核工業地質勘查局 物化探總隊，貴州 都勻 558004；3.中南大學 有色資源與地質災害探查湖南省重點實驗室，湖南 長沙 410083）

0 引 言

探地雷達（ground penetrating radar，GPR）數值模擬是地球物理重要研究內容，對建立的模型開展數值模擬能有效指導地質解釋.應用時域有限差分（finite difference time domain，FDTD）法進行GPR 數值模擬時，因計算機內存空間的有限，總是在某處把網格空間截斷，使之成為有限區域.這樣，在網格空間截斷處就會產生強烈的非物理的電磁反射干擾波.為了在截斷邊界處不引起虛假反射，吸收邊界條件對模擬精度起著關鍵作用.目前在截斷邊界處常通過加載完全匹配層（perfectly matched layer，PML）來 改 善 吸 收 效果，但無論是PML 還是單軸各向異性完全匹配層（uniaxial anisotropic perfectly matched layer，UPML）對隱失波、低頻波、大角度入射的掠角波都不能很好地吸收.為此，本文引入了復頻移卷積完全匹配層（complex frequency shifted convolution perfectly matched layer，CFS－PML）以改善截斷邊界處對隱失波、低頻波、掠角波的吸收效果.

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

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

1 零極點匹配Z變換與逆Laplace變換

根據Kuzuoglu等［14］的研究，定義復坐標拉伸坐標系中的算子與拉伸變量分別為

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

其中

然后利用零極點匹配Z變換（matchedZtransform）關系：

將式（3）變換到Z域

其中

求式（3）的倒數為

其中

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

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

2 二維TM 波的CFS－PML 邊界條件推導

由電磁波傳播理論可知［27］，GPR遵循的Maxwell方程的6個耦合偏微分方程在二維情況下變為3個偏微分方程，共有Ez、Hx、Hy3個分量，稱之為TM 波：

其中E為電場強度，V／m；H為磁場強度，A／m；ε為介電常數，F／m；μ為磁導率，H／m；σe為電導率，S／m；σm為等效磁阻率，Ω／m.

將頻率域方程（13）～（15）進行反Fourier變換，并將式（10）代入，可得到時域方程

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

其中

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

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

符號標記CA（m）～CRy（m）中的電磁參數標號m的取值與式（21）～（23）中右端電場或磁場分量的空間位置相同，其表達式如下：

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

其中d（0≤d≤δ）為CFS－PML 內的計算點到CFS－PML區域內邊界的距離；δ為CFS－PML 區域的厚度，計算過程中常取8～10個空間網格步長.m被稱為CFS－PML 的指數參數，κkmax為參數κk的最大值.對于參數σx和σy有

因為坐標擴展變量是一維函數，這里κk和σk只沿k（x或y）方向變化.為了衰減波的有效吸收，在常見的程序計算中常設置：

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

3.1 均勻介質中邊界條件對比

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

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

圖2為不加載任何邊界條件正演所得的15、22、30ns雷達波場快照；由圖2（a）的15ns快照中可見，雷達波以圓形方式向外傳播；圖2（b）雷達波在截斷邊界的4條邊、4個角點引起了較強的虛假反射波；圖2（c）截斷反射波又以原波前相反方向傳播，從左邊的電場數量級上可以看出，該截斷反射波的能量較強，取值在－200～200V／m變化，幾乎與一次波前的能量在一個數量級.

圖3（a）為采用二階Mur吸收邊界處理截斷邊界后的30ns波場快照，4條邊的人為截斷強反射波得到了很好的吸收，僅4個角點處的邊界反射波仍然可見；但從左邊的電場取值范圍0.5～1.0V／m 不難看出，Mur邊界條件對邊界處反射波起了較好的吸收作用，但是效果仍然有待提高.圖3（b）、（c）分別為加入UPML 吸收邊界、CFSPML吸收邊界30ns的波場快照，內部區域仍清晰可見4條截斷邊界處的反射，但是4個角點處的反射消失不見，再看左邊的電場取值范圍，顯然4條邊界的反射波能量是非常弱的，其能量在－0.1～0.1 V／m 變化，說明這兩種邊界條件較二階Mur吸收邊界條件有數量級的提高.

3.2 UPML與CFS－PML吸收效果對比驗證

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

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

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

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

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

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

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

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

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

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

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

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

4 結 論

（1）推導了CFS－PML邊界條件的GPR 正演FDTD 差分公式，應用均勻介質中GPR 波場邊界處的傳播快照，對比沒有加載邊界條件、加載Mur、UPML、CFS－PML 邊界處的波場反射波的強弱，實例結果表明UPML 與CFS－PML 能較好地處理GPR 截斷邊界，具有較好的吸收效果.

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

［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］ 葛德彪，閆玉波.電磁波時域有限差分方法［M］.西安：西安電子科技大學出版社，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］ 詹應林，昌彥君，曹中林.基于UPML 吸收邊界的探地雷達數值模擬研究［J］.資源環境與工程，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 全波場正演［J］.中南大學學報（自然科學版），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］ 李 靜，曾昭發，黃 玲，等.三維探地雷達數值模擬中UPML 邊界研究［J］.物探化探計算技術，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－hong.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］ 張顯文，韓立國，黃 玲，等.基于遞歸積分的復頻移PML彈性波方程交錯網格高階差分法［J］.地球物理學報，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］ 張魯新，符力耘，裴正林.不分裂卷積完全匹配層與旋轉交錯網格有限差分在孔隙彈性介質模擬中的應用［J］.地 球 物 理 學 報，2010，53（10）：2470－2483.ZHANG Lu－xin，FU 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］ 李建雄.時域有限差分法中完全匹配層的實現算法研究［D］.天津：天津大學，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.