馮德山, 楊良勇, 王珣

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院, 長沙 410083 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 有色金屬成礦預(yù)測教育部重點(diǎn)實(shí)驗(yàn)室, 長沙 410083

?

探地雷達(dá)FDTD數(shù)值模擬中不分裂卷積完全匹配層對倏逝波的吸收效果研究

馮德山1,3, 楊良勇1,2, 王珣1,3

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院, 長沙 410083 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 有色金屬成礦預(yù)測教育部重點(diǎn)實(shí)驗(yàn)室, 長沙 410083

介紹了CPML邊界條件的原理，推導(dǎo)了CPML的GPR正演FDTD差分公式，對比分析了Berenger PML、UPML、CPML三種PML對倏逝波的吸收性能.開展了PML邊界中關(guān)鍵參數(shù)κ和α的選取實(shí)驗(yàn)，確定了參數(shù)的取值范圍與選取原則.然后，以二維TM波為例，研究了倏逝波產(chǎn)生的機(jī)理，分析了決定逝波性吸收性能的影響因素.均勻介質(zhì)的波場快照、檢測點(diǎn)的反射誤差及全局反射誤差對比，說明了3種邊界條件對傳輸波都具有較好的吸收能力，而對低頻倏逝波的吸收表現(xiàn)迥異，其中CPML因?yàn)橐肓藚?shù)α，對倏逝波的吸收效果最佳，但離散化造成的全域誤差也最大.最后，應(yīng)用加載UPML和CPML邊界條件的FDTD程序，開展了GPR二維剖面法、寬角法矩狀地電模型及三維復(fù)雜模型的正演，展示了倏逝波反射對雷達(dá)正演剖面及波場快照的影響.進(jìn)一步對比了UPML與CPML對倏逝波的吸收表現(xiàn)優(yōu)劣，結(jié)果顯示，CPML可有效減少邊界反射誤差，并能取得滿意的精度，綜合考慮對倏逝波的吸收、全域誤差、編程難易程度等因素，在GPR正演中推薦使用CPML.

探地雷達(dá)； 不分裂卷積完全匹配層； 單軸各向異性完全匹配層； 倏逝波； 時(shí)域有限差分法

1 引言

探地雷達(dá)(Ground Penetrating Radar，GPR)是一種對地下或物體內(nèi)部不可見目標(biāo)體進(jìn)行定位的電磁無損探測技術(shù)，它根據(jù)地下介質(zhì)的電性差異，在雷達(dá)剖面上進(jìn)行成像，可分析地下介質(zhì)的電性結(jié)構(gòu)及空間形態(tài)(曾昭發(fā)等，2010；馮德山等，2014).為了更好地刻畫雷達(dá)波動現(xiàn)象、理解其傳播規(guī)律與物理過程，有必要開展GPR波動方程正演.GPR正演的模擬方法有許多，包括李展輝等(2009)將基于Runge Kutta方法的時(shí)間步迭代求解錯格時(shí)域偽譜法應(yīng)用于三維井中GPR模擬，有效解決了時(shí)域偽譜法的Gibbs現(xiàn)象；方宏遠(yuǎn)等(2012,2013)利用Hamilton系統(tǒng)的辛分塊龍格庫塔方法開展了道路結(jié)構(gòu)層雷達(dá)波傳播特征研究；底青云和王妙月(1999)推導(dǎo)了含衰減項(xiàng)的雷達(dá)波有限元方程，實(shí)現(xiàn)了GPR波的有限單元法正演；馮德山等(2012)提出了一種結(jié)合透射邊界與Sarma邊界的混合邊界條件，改善了有限單元法GPR正演的邊界反射.盡管GPR數(shù)值模擬方法種類許多，但時(shí)域有限差分法(Finite Difference Time Domain，F(xiàn)DTD)因具有直接時(shí)域計(jì)算、編程容易、節(jié)約存儲空間和計(jì)算時(shí)間、適合并行計(jì)算等優(yōu)點(diǎn)仍占據(jù)最重要的地位，得到了廣泛的應(yīng)用.Xu(1997)、Giannopoulos(2005)、Irving和Knight(2006)、Cassidy和Millington(2009)、李靜等(2010)、馮德山等(2010)應(yīng)用FDTD法進(jìn)行GPR數(shù)值模擬時(shí)，由于計(jì)算機(jī)內(nèi)存有限，模擬網(wǎng)格區(qū)域必須截?cái)酁橛邢迏^(qū)域，為了在截?cái)噙吔缣幉灰鹛摷俜瓷洌ǔＴ诰W(wǎng)格外圍加上吸收邊界條件(Absorbing Boundary Condition，ABC)來削弱邊界的反射(Berenger，1999).早期提出的Taflove & Brodwin(1975)邊界、Engquist-Majda邊界(1977)、Liao邊界(1984)有效地減弱了邊界反射；Mur(1981)在單行波的基礎(chǔ)上提出了Mur二階吸收邊界，得到了廣泛的應(yīng)用；Mei和Fang (1992)提出了超吸收邊界條件，它是對單行波方程為基礎(chǔ)的吸收邊界條件的一種有效改進(jìn)，但是在模擬區(qū)域的角點(diǎn)處未能達(dá)到很好的改善.這些邊界條件通常在模擬區(qū)域的外邊界仍有0.5%～5%的反射系數(shù)，他們逐漸被完全匹配層(Perfectly Matched Layer，PML)所取代.

PML邊界條件最初是由Berenger(1994)提出的，它在FDTD網(wǎng)格外邊界加載一種非物理的吸收媒質(zhì)，可以使電磁波在匹配層中無反射地通過并按指數(shù)規(guī)律衰減，相比于Mur二階吸收邊界，完全匹配層的吸收效果可以提高幾個(gè)數(shù)量級，將邊界條件的研究向前推進(jìn)了一大步.Katz等(1994)將Berenger PML的使用范圍延伸到三維模擬區(qū)域；Moerloose和Stuchly(1995)通過分析和數(shù)值實(shí)驗(yàn)指出Berenger PML對倏逝波的吸收毫無效果，且在低頻情況下，會存在很大的數(shù)值誤差；Berenger(1997)提出由倏逝波產(chǎn)生的數(shù)值反射主要存在于與PML電導(dǎo)率有關(guān)的截止頻率之下，并且指出了Moerloose和Stuchly論文中僅考慮了傳輸波平行于PML界面的情況，理論上，倏逝波越強(qiáng)，PML對其的吸收效果也越好，但是由于時(shí)間上和空間上的離散，效果不盡人意；Berenger(1999)對波導(dǎo)應(yīng)用中的倏逝波的反射進(jìn)行了系統(tǒng)的討論，并且對電導(dǎo)率遞增的PML層的倏逝波反射系數(shù)進(jìn)行了研究.Berenger PML的理論體系是非Maxwell方程的，物理機(jī)制模糊，同時(shí)，其電磁場分量分裂技術(shù)增加了計(jì)算內(nèi)存與數(shù)值實(shí)現(xiàn)的難度(葛德彪和閆玉波，2011；2005)，而且只對行波有吸收效果，對空域中衰減的隱失波、低頻波以及入射角度較小的掠角波無吸收效果.為了改善邊界條件的吸收效果.Sacks(1995)和Gedeny(1996)提出了單軸各向異性完全匹配層(Uniaxisal PML，UPML)，UPML邊界無須對電磁場分裂，計(jì)算更簡潔、易編程實(shí)現(xiàn)，因其在吸收系數(shù)中引入了線性的吸收因子，UPML邊界對隱失波、凋落波等干擾波有一定的吸收效果，能適用于有耗、無耗、色散等介質(zhì)中.詹應(yīng)林(2008)將UPML應(yīng)用于二維FDTD探地雷達(dá)正演中；馮德山等(2011)將UPML邊界應(yīng)用于ADI-FDTD差分中實(shí)現(xiàn)了GPR的三維數(shù)值模擬.盡管UPML邊界對于凋落波、隱失波、低頻波的干擾波有一定的吸收效果，但仍然不是特別完美.為了克服UPML低頻區(qū)性能變差的缺點(diǎn)，Kuzuoglu和Mittra(1996)考慮到復(fù)頻移拉伸函數(shù)的PML技術(shù)可以吸收低頻信號在邊界表面的反射，同時(shí)對隱失波等大角度掠射波也具有好的吸收效果，為此，提出了復(fù)頻偏移技術(shù)(Complex Frequency Shifted，CFS)加強(qiáng)了對低頻、隱失波的吸收；Roden和Gedney(2000)將復(fù)頻移拉伸函數(shù)加入到電磁波模擬的PML方程，有效壓制了信號的假反射；Berenger(2002)驗(yàn)證了CFS PML能夠?qū)Φ皖l區(qū)倏逝波取得非常優(yōu)越的吸收效果；Wang和Liang(2006)將復(fù)頻移PML邊界(Complex Frequency Shift Perfectly Matched Layer, CFS-PML)應(yīng)用到ADI-FDTD中；李展輝和黃清華(2014)將CFS-PML應(yīng)用在瞬變電磁法正演中，以替代傳統(tǒng)狄利克雷邊界條件；Drossaert和Giannopoulos(2007)以及 Dimitri和Roland(2007)提出了基于遞歸積分的波場非分裂的復(fù)頻移PML邊界(Recursive Integration CFS-PML, 簡稱CFS-RIPML)；Zhang和Shen(2010)及Martin等(2010)將基于輔助微分方程的CFS-PML應(yīng)用到彈性波數(shù)值模擬中，該方法易于擴(kuò)展到高階差分形式；張顯文等(2009)將遞歸復(fù)頻移PML應(yīng)用到交錯網(wǎng)格高階差分彈性波波動方程正演中；Giannopoulos(2008)及Li等(2012)將基于CFS-RIPML吸收邊界技術(shù)應(yīng)用在電磁波波動方程的有限差分正演中.Roden和Gedney(2000)在CFS-PML基礎(chǔ)上應(yīng)用卷積技術(shù)提出了卷積完全匹配層(Convolutional PML，CPML)，大大提高了運(yùn)算效率.張魯新等(2010)將不分裂的完全匹配層(CPML)結(jié)合旋轉(zhuǎn)交錯網(wǎng)格有限差分技術(shù)應(yīng)用到孔隙彈性介質(zhì)模擬中；相比常規(guī)PML、UPML吸收邊界，CPML對于隱失波、低頻波等具有更好的吸收效果，所以 CPML具有非常廣闊的發(fā)展前景.

本文分析了PML邊界中各參數(shù)的設(shè)置，采用波場快照著重探討了二、三維UPML、CPML對倏逝波的吸收效果對比.然后，采用GPR二維剖面法和寬角法正演剖面直觀地展示了UPML和CPML兩種邊界條件下倏逝波的反射對正演剖面的影響，同時(shí)在三維情況下展示了復(fù)雜模型的波場快照，直觀地顯示了倏逝波的形態(tài).結(jié)果顯示，CPML可有效減少邊界反射誤差，并能取得滿意的精度，為GPR的正演提供了合適的吸收邊界.

2 卷積完全匹配層原理

CPML吸收邊界條件是建立在伸縮坐標(biāo)PML基礎(chǔ)上.在伸縮坐標(biāo)中，z方向的安培環(huán)路定律為

(1)

其中si為坐標(biāo)伸縮因子：

(2)

把公式(1)從頻域轉(zhuǎn)換到時(shí)域，由于坐標(biāo)伸縮因子與頻率有關(guān)，因此轉(zhuǎn)換到時(shí)域后公式右邊會存在卷積形式，即

(3)

(4)

(5)

其中，δ(t)和u(t)分別是單位沖擊函數(shù)和單位階躍函數(shù)，將(5)式代入到(3)式中，有

(6)

直接在遞推方程(6)中計(jì)算時(shí)域卷積的效率是很低的，為了提高卷積計(jì)算效率，將ξi(t)離散沖擊響應(yīng)定義為

(7)

式(7)中

).

(8)

根據(jù)公式(7)和(8)，按照Yee氏網(wǎng)格將(6)式進(jìn)行時(shí)間和空間上的離散，可得

(9)

公式(9)中的離散卷積計(jì)算起來十分復(fù)雜，但是，由于Z0i(m)是簡單的指數(shù)形式，從而它們的和可以使用遞歸卷積來遞歸得到(葛德彪和閆玉波，2011；2005).因此，引入一組新的輔助表達(dá)式ψi.

(10)

式(10)中

bi=e， (i=x,y,z)

(13)

其中ai已在式(8)中給出，進(jìn)一步整理公式(10)，可得

(14)

式(14)中標(biāo)號m=(i,j,k+1/2)，CA、CB與常規(guī)FDTD中表示的內(nèi)容相同.其余場量可類似推出.由式(14)可看出，除κ項(xiàng)以及后面的ψ項(xiàng)外，公式(14)前面部分與常規(guī)FDTD公式十分相似.故編程時(shí)可令模擬區(qū)域內(nèi)κ=1，而PML層內(nèi)令κ漸進(jìn)變化，從而將模擬區(qū)域與PML邊界層中采用統(tǒng)一的差分公式，然后，再在PML區(qū)域中加上后面的ψ項(xiàng)，這種編程方式可使程序更簡潔,有效減少計(jì)算量.為了減少反射誤差，在PML層內(nèi)σ、κ、α參數(shù)分布為非均勻的，各參數(shù)單調(diào)變化，以z方向?yàn)槔?，通常?/p>

(15)

其中d為PML層網(wǎng)格數(shù)，z0為PML與模擬區(qū)域分界處的網(wǎng)格編號.由于電場和磁場在Yee網(wǎng)格中空間位置相差半個(gè)網(wǎng)格，故σ、κ、α等參數(shù)也需在PML網(wǎng)格上根據(jù)電場和磁場的對應(yīng)位置來分別設(shè)置，同時(shí)，z0代表的分界處對于電場和磁場也有所不同，需要相差半個(gè)空間網(wǎng)格.

3 完全匹配層對倏逝波的吸收效果

在PML邊界條件中，BerengerPML、UPML、CPML這3種完全匹配層應(yīng)用較廣泛，對于傳輸波都具有很好的吸收效果，但對倏逝波的吸收，尤其是低頻倏逝波的吸收，卻有著截然不同的效果.倏逝波(evanescentwave)又稱凋落波、隱失波(Berenger，1998).以二維TM波為例對比這3種邊界條件對倏逝波的吸收能力.假設(shè)真空中雷達(dá)傳輸波的傳播方向?yàn)閄，倏逝波的方向?yàn)閅，X與坐標(biāo)軸x軸夾角為θ，其中θ的選取為負(fù)數(shù)，任一場量的場值都可寫為(Berenger，2000)

ψvacuum=ψ0ejω t-j(ω cosh/c)Xe-(ωsinh/c)Y,

(16)

指數(shù)項(xiàng)e-(ω sinh/c)Y為倏逝波的自然衰減項(xiàng).在BerengerPML中，假設(shè)PML層垂直于x軸，即σy=0，任一分量的場值可寫為

(17)

|ψvacuum|e.

(18)

分析衰減項(xiàng)表達(dá)式，結(jié)合倏逝波cosh>1可知，理論上，PML對倏逝波的吸收能力較傳輸波更強(qiáng)，倏逝波越強(qiáng)，衰減越快.同時(shí)衰減項(xiàng)與入射角θ存在聯(lián)系，入射角越大，對倏逝波的吸收能力越弱.當(dāng)倏逝波足夠強(qiáng)時(shí)，即cosh足夠大，倏逝波可以在極小的距離x內(nèi)充分吸收.但是在FDTD離散網(wǎng)格中，在cosh足夠大時(shí)，倏逝波甚至可在小于一個(gè)網(wǎng)格距離內(nèi)被吸收.顯然，有限的FDTD網(wǎng)格無法完成變化巨大的倏逝波的采樣工作，所以必然會造成劇烈的反射誤差，尤其是在低頻區(qū).Berenger(2000)通過數(shù)值計(jì)算驗(yàn)證指出，在低頻區(qū)，真空與均勻PML分界面上，傳輸波相對于倏逝波的反射系數(shù)更小，而倏逝波越強(qiáng)，其反射系數(shù)也越大.所以倏逝波強(qiáng)烈反射的關(guān)鍵原因在于PML層的離散化.在BerengerPML中，可通過增加PML層的厚度來壓制倏逝波.

對于UPML，其關(guān)鍵參數(shù)為s=κ+σ/jωε0，當(dāng)κmax=1時(shí)，也即κ=1，可以證明，與BerengerPML完全等價(jià).在離散FDTD中計(jì)算時(shí)，UPML采用的是半隱式差分，而BerengerPML在PML中采用的是指數(shù)差分，而且編程平臺的精度有限，故在數(shù)值計(jì)算時(shí)，難免會存在舍入誤差.但兩者的數(shù)值誤差可以忽略不計(jì)，所以在數(shù)值上，κmax=1時(shí)的UPML與BerengerPML也等價(jià).在UPML中，當(dāng)κmax增大時(shí)，對倏逝波的吸收能力也加強(qiáng)，但隨之而來的代價(jià)是，由于κ的離散化，導(dǎo)致它在PML網(wǎng)格間遞增值也逐漸變大，UPML的全域誤差也劇烈增大.

CPML中，坐標(biāo)伸縮系數(shù)為si=κi+σi/(αi+jωε0)，很顯然，fα=αi/2πε0為一關(guān)鍵的分界參數(shù)，當(dāng)f?fα?xí)r，si趨近于si=κi+σi/jωε0，此時(shí)CPML與UPML無異.特別地，當(dāng)κi=1時(shí)，|ψPML|變?yōu)?/p>

(19)

可以看到公式(19)與公式(17)相同，即在高頻時(shí)，CPML與BerengerPML、UPML的作用相當(dāng).當(dāng)f?fα?xí)r，si趨近于si=κi+σi/αi，即si變?yōu)閷?shí)坐標(biāo)伸縮因子.CPML中任一場分量的值可寫為

(20)

特別地，當(dāng)κi=1時(shí)，將fα的值代入，|ψPML|可變?yōu)?/p>

(21)

倘若將ψvacuum中的自然衰減項(xiàng)合并進(jìn)去，則式(21)中總的衰減相當(dāng)于倏逝波在距離(1+σx/αx)x上的自然衰減.分析式(21)可知，當(dāng)f?fα?xí)r，CPML并不能對傳輸波(sinh=0)起到吸收作用，但是對倏逝波能起到出色的吸收作用(sinh>0)，這是其他PML所沒有的優(yōu)點(diǎn).公式(21)中的指數(shù)衰減項(xiàng)要小于公式(19)中的指數(shù)衰減項(xiàng).對于BerengerPML來說，正是因?yàn)橘渴挪ㄋp太快而不能被網(wǎng)格合理采樣，所以在低頻區(qū)吸收性能較差，而CPML中的指數(shù)衰減項(xiàng)的衰減能力適中，雖然還是存在一定的反射，但是基本上能夠被網(wǎng)格合理采樣，所以CPML對倏逝波的吸收能力要強(qiáng)于BerengerPML.注意到式(21)中，指數(shù)系數(shù)與頻率相關(guān)，即f和sinh存在一定的變動，但實(shí)際上，在頻率在fα之下的區(qū)域，fsinh接近于常數(shù)，故CPML在低頻區(qū)能取得穩(wěn)定的倏逝波吸收效果.

4 倏逝波的吸收數(shù)值實(shí)驗(yàn)

從BerengerPML、UPML、CPML在高頻區(qū)時(shí)的波場幅值公式中可看出，3種PML對倏逝波的吸收性能都與cosθ相關(guān).當(dāng)θ趨近于π/2時(shí)，cosθ趨近于0，此時(shí)傳輸波的方向與PML表面接近于平行時(shí)，PML對倏逝波的吸收效果最弱.為了清楚地反應(yīng)倏逝波的反射，分析和對比倏逝波吸收的相關(guān)參數(shù)，開展大角度入射下的TM波實(shí)驗(yàn).將激勵源置于模擬區(qū)域的左上角距左分界面和上分界面一個(gè)網(wǎng)格處，模擬區(qū)域統(tǒng)一設(shè)為真空.激勵源為主頻900MHz的blackman-harris脈沖，取空間步長為0.012m，時(shí)間步長為0.012ns，UPML中參數(shù)κmax分別等于1、3、5、7、9、11，CPML中κmax=1、αmax=0.008.為了能夠更好地呈現(xiàn)出低頻倏逝波異常反射，圖中右邊的顏色棒(colorbar)統(tǒng)一采用lines顯示方式.圖1為不同PML參數(shù)條件下400時(shí)間步的電場Ez波場快照，在圖1(a—g)圖中，左分界面和上分界面處都出現(xiàn)比較強(qiáng)烈的倏逝波的反射.從圖1(b—g)中看出，隨著κmax的增大，UPML對倏逝波的吸收效果逐漸改善，但并不是κmax越大越好，如果κmax過大，會因?yàn)樵赑ML中分布時(shí)增加速率太大而造成強(qiáng)烈的后期虛假反射，通常取κmax=5～11為宜.注意到，在κmax=1時(shí)，UPML與BerengerPML在數(shù)學(xué)與物理是完全等價(jià)的，所以在圖1a和圖1b中Ez波形完全一致.

圖1h為邊界條件采用CPML時(shí)同時(shí)刻的雷達(dá)波場快照，其中主要參數(shù)設(shè)置為κmax=1、αmax=0.008.可以看到CPML對倏逝波取得了非常出色的吸收效果.通過分析這個(gè)實(shí)驗(yàn)可知，BerengerPML是UPML參數(shù)κmax=1時(shí)的特例，通過選取合適的κmax值可以改善UPML邊界對倏逝波吸收效果，而CPML中參數(shù)α對倏逝波吸收效果具有非常重要的作用，通過選取合適的κmax及αmax的參數(shù)組合，可以達(dá)到對倏逝波的充分吸收.

圖1 Berenger PML、UPML、CPML三種吸收邊界在400時(shí)間步時(shí)Ez的波場快照(a) Berenger PML；(b) κmax=1 的UPML；(c) κmax=3的UPML； (d) κmax=5的UPML；(e) κmax=7的UPML； (f) κmax=9的UPML； (g) κmax=11的UPML； (h) κmax=1、αmax=0.008的CPML.Fig.1 The Ez snapshot of Berenger PML、UPML、CPML in the timestep of 400

圖2 三種PML對倏逝波吸收性能的檢驗(yàn)?zāi)Ｐ图皺z測點(diǎn)分布示意圖Fig.2 The testing model and distribution of detection point for absorption properties of evanescent wave in three kinds of PML

為了進(jìn)一步對比BerengerPML、UPML、CPML對倏逝波的吸收效果，設(shè)置圖2所示的簡單模型，模型網(wǎng)格數(shù)量為200×200，通過設(shè)置一個(gè)檢測點(diǎn)來分析3種不同邊界條件下的反射誤差.脈沖源點(diǎn)A距上側(cè)和左側(cè)PML層僅一個(gè)網(wǎng)格，檢測(接收)點(diǎn)B在A點(diǎn)右側(cè)100個(gè)網(wǎng)格處，激勵源為900MHz的Blackman-harris脈沖，空間步長設(shè)為0.012m，時(shí)間步長為0.025ns.

(22)應(yīng)用式(22)計(jì)算3種PML所得B點(diǎn)反射誤差如圖3所示，其中UPML中參數(shù)κmax=5，CPML中κmax=5，αmax=0.008.由圖3可看出，CPML對倏逝波的吸收效果有很大的改善，大約提高了一個(gè)數(shù)量級.至于后期反射誤差劇烈抖動，推斷是κ和α離散化所致.對比于BerengerPML和UPML可知，κmax=5時(shí)，UPML對倏逝波的吸收有一定的改善，但是仍達(dá)不到CPML邊界條件的吸收效果.

(23)

BerengerPML、κmax=5的UPML和κmax=5，αmax=0.008的CPML的全域誤差如圖4所示.圖中可見，BerengerPML的全域誤差最小.在相同的κ情況下，CPML和UPML的全域誤差在同一個(gè)數(shù)量級，CPML的全域誤差要比UPML的大，但都在合理的范圍之內(nèi).因此在CPML與UPML之間權(quán)衡時(shí)，考慮到對倏逝波的吸收效果，CPML作為優(yōu)先考慮的邊界.

圖3 B點(diǎn)不同邊界條件的反射誤差對比Fig.3 Relative error at point B of different boundary conditions

圖4 不同邊界條件的全域誤差對比Fig.4 Contrast of global error under different boundary conditions

5 CPML邊界在GPR正演中的應(yīng)用及其對倏逝波的吸收效果分析

5.1 CPML在二維GPR正演中的應(yīng)用實(shí)例

為了更好地說明CPML邊界條件在GPR正演截?cái)噙吔缣帉Ψ瓷洳ǖ奈招Ч贛atlab平臺編寫了CPML邊界的二維雷達(dá)FDTD正演程序.設(shè)置模擬區(qū)域網(wǎng)格數(shù)量為200×200空間步長，空間步長為0.006 m，時(shí)間步長為0.015 ns，CPML邊界層厚度設(shè)為10個(gè)網(wǎng)格.為了驗(yàn)證程序的正確性，采用最簡單的均勻介質(zhì)模型，模擬區(qū)域內(nèi)介質(zhì)的相對介電常數(shù)為3.5，電導(dǎo)率為2.5 S·m-1.激勵源為主頻900 MHz的Blackman-harris脈沖，位于模擬區(qū)域正中心，CPML參數(shù)取κmax=5，αmax=0.008.

圖5為應(yīng)用FDTD算法正演所得的不同時(shí)間步電場Ez傳播快照.由圖5(a—b)可見，當(dāng)激勵源位于模擬區(qū)域中心時(shí)，電磁波以圓形方式向外傳播，這與電磁波傳播規(guī)律是一致的，從而證明了程序的正確性.從圖5c中可見，當(dāng)雷達(dá)波傳播至邊界處時(shí)，由于CPML邊界條件的加載，使截?cái)噙吔绲玫搅溯^好的處理，波形得以傳播出模擬區(qū)域而沒有產(chǎn)生強(qiáng)反射.圖5d說明雷達(dá)波形傳播出去后，區(qū)域內(nèi)部后期也未見截?cái)噙吔缣幍姆瓷洳?，說明了CPML邊界條件在正演中取得了滿意的吸收效果.

由于二維GPR探測中常采用剖面方式顯示雷達(dá)數(shù)據(jù)，而雷達(dá)數(shù)據(jù)采集方式主要分為兩類：剖面法與寬角法.為了加深CPML對二維雷達(dá)正演剖面中倏逝波的吸收效果的理解，建立圖6所示的矩狀地電模型，背景介質(zhì)介電常數(shù)為11.0，電導(dǎo)率為0.1 mS·m-1.矩形模型上頂面埋深為0.52 m，寬度為0.16 m，厚度為0.16 m，矩狀小異常體的介電常數(shù)6.0，電導(dǎo)率0.5 mS·m-1.應(yīng)用FDTD算法對該模型進(jìn)行剖面法和寬角法兩種采樣方式開展正演模擬.

為此，采用發(fā)射天線(激勵源)和接收天線(記錄點(diǎn))距為0.24 m，自左向右同步移動.網(wǎng)格步長為0.004 m，時(shí)間步長為0.02 ns，時(shí)窗長度為24 ns.激勵源仍為900 MHz的Blackman-harris脈沖，模擬區(qū)域的外邊界分別采用κmax=1的UPML和κmax=5，αmax=0.008的CPML，激勵源的位置均距上側(cè)PML內(nèi)表面1個(gè)網(wǎng)格.圖7分別呈現(xiàn)了UPML與CPML波場快照，圖中白色線框即為異常體.從圖7a中可看出，在邊界處存在著明顯的倏逝波反射，而在圖7b中，倏逝波的反射幾乎不可見.

圖5 加載CPML邊界的不同時(shí)刻FDTD法正演雷達(dá)波場快照Fig.5 The snapshot of radar wavefeild in difference moment by FDTD with CPML boundary

二維情況下，倏逝波的反射開始出現(xiàn)的位置是在PML內(nèi)表面附近且距離激勵源約1個(gè)波長處，故在剖面法中，若自激自收或收發(fā)天線位置靠近時(shí)，則記錄剖面上并沒有倏逝波的反射出現(xiàn).但倘若收發(fā)天線水平方向上存在著一定距離，則雷達(dá)剖面上會存在著明顯的倏逝波反射.

圖6 矩狀地電模型Fig.6 Rectangular geoelectric model

圖8a為采用UPML的雷達(dá)正演剖面圖，圖中可見，在直達(dá)波之后存在著連續(xù)的倏逝波反射干擾，如圖8a中黑色虛線框所示.而CPML對倏逝波的壓制效果比較理想，在正演剖面圖8b中未出現(xiàn)倏逝波反射的干擾.對比圖8(a—b)可知，由于倏逝波的反射的存在會一定程度上影響正演剖面的正確性，所以在大收發(fā)距時(shí)CPML邊界要優(yōu)于UPML邊界.

再考慮該矩狀模型寬角法的正演結(jié)果，將發(fā)射天線置于矩狀模型的正上方，接收天線置于發(fā)射天線兩側(cè)，共記錄200道波形.設(shè)置網(wǎng)格步長為0.004 m，時(shí)間步長為0.02 ns，模擬時(shí)窗長度為24 ns，激勵源仍為900 MHz的Blackman-harris脈沖，激勵源的位置距離上側(cè)PML內(nèi)表面1個(gè)網(wǎng)格.圖9(a—c)分別為加載κmax=1與κmax=5的UPML和κmax=5，αmax=0.008的CPML邊界條件雷達(dá)正演剖面.圖9a中UPML內(nèi)表面處的倏逝波的反射依次到達(dá)記錄點(diǎn)，在正演剖面中形成黑色虛線所圈的“掃帚形”的干擾波，即使采用圖9b中κmax=5的UPML倏逝波反射有一定的減弱，但在正演剖面中仍然存在；分析圖9c可見，加載了CPML邊界的寬角法雷達(dá)正演剖面中未出現(xiàn)倏逝波反射的干擾，說明CPML邊界條件對倏逝波具有良好的吸收特性.

圖7 矩狀模型正演波場快照(a) UPML邊界Ez波場快照; (b) CPML邊界Ez波場快照.Fig.7 The snapshot of Ez with rectangular geoelectric model(a) Ez snapshot with UPML boundary; (b) Ez snapshot with CPMLboundary.

圖8 矩狀模型剖面法探地雷達(dá)正演剖面(a) UPML邊界; (b) CPML邊界.Fig.8 The forward modelling profile of GPR using profile method to rectangular geoelectric model(a) UPML boundary; (b) CPMLboundary.

圖9 矩狀模型寬角法探地雷達(dá)正演剖面(a) κmax=1的UPML邊界; (b) κmax=5的UPML邊界; (c) κmax=5，α=0.008的CPML邊界.Fig.9 The forward modelling profile of GPR using wide angle method to rectangular geoelectric model

圖10 加載UPML、CPML邊界在不同時(shí)間步時(shí)電場Ez的三維雷達(dá)波場快照(a) 180個(gè)時(shí)間步的UPML波場快照; (b) 260個(gè)時(shí)間步的UPML波場快照; (c) 340個(gè)時(shí)間步的UPML波場快照;(d) 180個(gè)時(shí)間步的CPML波場快照; (e) 260個(gè)時(shí)間步的CPML波場快照; (f) 340個(gè)時(shí)間步的CPML波場快照.Fig.10 The Ez wavefield snapshot in different moment with UPML ,CPML boundary(a) The UPML wavefield snapshot in the moment of timestep 180; (b) The UPML wavefield snapshot in the moment of timestep 260;(c) The UPML wavefield snapshot in the moment of timestep 340; (d) The CPML wavefield snapshot in the moment of timestep 180;(e) The CPML wavefield snapshot in the moment of timestep 260; (f) The CPML wavefield snapshot in the moment of timestep 340.

5.2 CPML在三維GPR正演中的應(yīng)用實(shí)例

實(shí)例分析三維GPR正演中CPML對倏逝波的吸收效果，設(shè)置三維模擬區(qū)域大小為200×200×200個(gè)空間步長，介質(zhì)為真空介質(zhì)，空間步長為0.01 m，時(shí)間步長為0.0167 ns.激勵源與上例一致，為了較好地體現(xiàn)倏逝波的反射情況，方便分析邊界條件對倏逝波的吸收效果，將源置于x、y平面中心，且距上側(cè)PML層1個(gè)網(wǎng)格，便于雷達(dá)波在大角度情況下入射到PML邊界.PML吸收邊界分別加載κmax=1的UPML和κmax=5，αmax=0.008的CPML.

圖10(a—c)為加載UPML邊界條件下180、260、340三個(gè)時(shí)間步的Ez波場快照；圖10(d—f)為加載CPML邊界條件下180、260、340三個(gè)時(shí)間步的Ez波場快照.分析圖10(a—c)可知，隨著雷達(dá)波前向外擴(kuò)散，加載UPML邊界的三維雷達(dá)波形中出現(xiàn)了明顯的倏逝波反射現(xiàn)象，并且倏逝波的反射沿著UPML的內(nèi)表面向外擴(kuò)散.而圖10(d—f)加載了CPML邊界的三維波形中倏逝波能量非常弱、幾乎不可見.兩者對比說明：相比于UPML，CPML對倏逝波具有更好的吸收效果.

圖11 三維情形下正方體與球體異常的復(fù)雜模型Fig.11 Cube and sphere objection in three-dimension

為突出在復(fù)雜介質(zhì)下CPML吸收邊界對倏逝波的吸收效果，建立如圖11所示的三維復(fù)雜模型.背景介質(zhì)的介電常數(shù)為5，電導(dǎo)率為0，其大小為1 m×1 m×1 m.模型中存在兩個(gè)異常體，分別為邊長0.03 m的正方體和直徑0.06 m的球體.異常正方體與異常球體的介電常數(shù)為3.0，電導(dǎo)率為3 mS·m-1.激勵源與上例一致，同樣將源置于x、y平面中心.PML吸收邊界分別加載κmax=1的UPML和κmax=5，αmax=0.05的CPML.圖12a為加載UPML邊界條件下的Ez波場快照；圖12b為加載CPML邊界條件的Ez波場快照；圖12c為z方向上的兩道切片，可以看到第一道切片中倏逝波的反射呈規(guī)則的同心圓狀，在第二道切片中可以看到兩個(gè)異常體的明顯散射；圖12d中同樣有兩道切片，在第一道切片中倏逝波的反射明顯減弱.由此可見，雷達(dá)波遇到正方體異常和球體異常時(shí)會產(chǎn)生散射.UPML邊界的波場快照中，上部存在著明顯的倏逝波反射.而CPML邊界的波場快照中，倏逝波的反射明顯要小.所以在GPR正演中，倏逝波的反射會嚴(yán)重影響波場快照和正演剖面的正確性，有必要采取更好的吸收邊界來壓制倏逝波的反射.

6 結(jié)論

(1) 探討了Berenger PML、UPML、CPML三種PML邊界條件的內(nèi)在關(guān)系，對比了它們對倏逝波的吸收效果.Berenger PML由于阻抗匹配條件的特殊性，主要應(yīng)用于真空介質(zhì).理論上，Berenger PML對倏逝波的吸收能力要強(qiáng)于雷達(dá)傳輸波，但實(shí)際中由于FDTD的離散，特別是在低頻區(qū)，Berenger PML對倏逝波的吸收能力要遠(yuǎn)差于傳輸波，但可以通過加大PML層的厚度來壓制倏逝波的反射.相比于Berenger PML，UPML不需分裂場量，而且編程相對簡單，其適用范圍更廣.CPML為基于復(fù)頻移和卷積計(jì)算技術(shù)的完全匹配層，對倏逝波，特別是低頻倏逝波，有著出色的吸收效果.

(2) 詳細(xì)探討了PML邊界條件中關(guān)鍵參數(shù)κ和α對倏逝波吸收效果的影響.實(shí)驗(yàn)表明：UPML中，參數(shù)κmax通常在1～11范圍內(nèi)取值較為合適，當(dāng)參數(shù)κmax=1時(shí)的UPML與Berenger PML在數(shù)學(xué)上和物理上完全等價(jià).隨著κmax的增加，UPML對倏逝波的吸收能力加強(qiáng)，但是當(dāng)κmax過大超過取值范圍時(shí)，其離散所引起的全域誤差會影響計(jì)算精度.而CPML中有兩個(gè)關(guān)鍵參數(shù)κ和α，其中αmax通常在0～0.05范圍內(nèi)，當(dāng)αmax過大時(shí)，也會增大模擬區(qū)域中的全域誤差，與UPML相同的κ值情況下，CPML的全域誤差要大于UPML，但兩者全域誤差在同一個(gè)數(shù)量級.通常二維情況下CPML在κmax=1，αmax=0.008左右時(shí)可取得不錯的倏逝波吸收效果.

圖12 加載UPML、CPML邊界下復(fù)雜模型的三維雷達(dá)波場快照(a) UPML波場快照; (b) CPML波場快照; (c) UPML波場快照，z方向的2個(gè)切片; (d) CPML波場快照，z方向的2個(gè)切片.Fig.12 Three-dimensional radar snapshots of complex models under UPML and CPML boundary(a) The UPML wavefield snapshot; (b) The CPML wavefield snapshot; (c) The UPML wavefield snapshot; (d) The CPML wavefield snapshot.

(3) 分別將UPML與CPML加載于GPR正演的邊界條件處理中，了解到倏逝波的反射常出現(xiàn)在PML內(nèi)表面的特性，驗(yàn)證了倏逝波反射的存在，而且倏逝波的反射通常與雷達(dá)傳輸波入射到PML層的角度有關(guān)，入射角越大，倏逝波的反射越強(qiáng).倏逝波的存在會對雷達(dá)正演剖面的正確性造成干擾，因此，有必要壓制倏逝波的反射.實(shí)驗(yàn)證明：CPML不僅對雷達(dá)傳輸波有著良好的吸收能力，對倏逝波也能取得不錯的吸收效果.故在GPR正演的截?cái)嗵?，推薦加載CPML邊界條件.

Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves.JournalofComputationalPhysics, 114(2): 185-200.

Berenger J P. 1996. Perfectly matched layer for the FDTD solution of wave-structure interaction problems.IEEETransactionsonAntennasandPropagation, 44(1): 110-117.

Berenger J P. 1997. Numerical reflection of evanescent waves from perfectly matched layers.IEEEAntennasandPropagationSocietyInternationalSymposium,Montreal,Canada,July, (2): 1888-1891.Berenger J P. 1998.An effective PML for the absorption of evanescent waves in waveguides.IEEEMicrowaveandGuidedWaveLetters, 8(5): 188-190.Berenger J P. 1999. Evanescent waves in PML′s: Origin of the numerical reflection in wave-structure interaction problems.IEEETrans.AntennasPropag, 47(10): 1497-1503.

Berenger J P. 2000. Numerical reflection of evanescent waves by PMLs: Origin and interpretation in the FDTD case. Expected consequences to other finite methods.Int.J.Numer.Model, 13(2): 2-3.Berenger J P. 2002. Numerical reflection from FDTD-PMLs: a comparison of the split PML with the unsplit and CFS PMLs.IEEETransactionsonAntennasandPropagation, 50(3): 258-265.

Berenger J P. 2002. Application of the CFS PML to the absorption of evanescent waves in wave guides.IEEEMicrowaveandWirelessComponentsLetters, 12(6): 218-220.

Cassidy N J, Millington T M. 2009. The application of finite-difference time-domain modelling for the assessment of GPR in magnetically lossy materials.JournalofAppliedGeophysics, 67(4): 296-308.

Di Q Y, Wang M Y.1999. 2D finite element modeling for radar wave.ChineseJ.Geophys. (in Chinese), 42(6): 818-825.

Dimitri K，Roland M. 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation.Geophysics, 72(5): 255-167.Drossaert F H，Giannopoulos A. 2007. A nonsplit complex frequency-shifted PML based on recursive integration for FDTD modeling of elastic waves.Geophysics, 72(2): 9-17.Drossaert F H，Giannopoulos A. 2007. Complex frequency shifted convolution PML for FDTD modeling of elastic waves.WaveMotion, 44: 593-604.

Engquist B, Majda A. 1977. Absorbing boundary conditions for the numerical simulation of waves.Math.Comput., 31(139): 629-651.

Fang H Y, Lin G. 2013.Simulation of GPR wave propagation in complicated geoelectric model using symplectic method.ChineseJ.Geophys.(in Chinese), 56(2): 653-659.

Fang H Y, Lin G. 2012.Symplectic partitioned Runge-Kutta methods for two-dimensional numerical model of ground penetrating radar.Computers&Geosciences, 49: 323-329.Feng D S, Chen C S, Dai Q W. 2010. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD.ChineseJ.Geophys. (in Chinese), 53(10): 2484-2496.

Feng D S, Chen C S, Wang H H.2012. Finite element method GPR forward simulation based on mixed boundary condition.ChineseJ.Geophys.(in Chinese), 55(11): 3774-3785.Feng D S,Xie Y. 2011. Three dimensional GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD radar.JournalofCentralSouthUniversity:ScienceandTechnology, 42(8): 2363-2372.Feng D S, Chen J W, Wu Q. 2014. A hybrid ADI-FDTD subgridding scheme for efficient GPR simulation of dispersion media.ChineseJ.Geophys. (in Chinese), 57(4): 1322-1334.

Giannopoulos A. 2005. Modelling ground penetrating radar by GprMax.ConstructionandBuildingMaterials, 19(10): 755-762.

Gedney S D. 1996. An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices.IEEETransactionsonAntennasandPropagation,44 (12): 1630-1639.

Gedney S D. 1996. An Anisotropic PML absorbing media for the FDTD simulation of fields in lossy and dispersive media.Electromagnetics, 16(4): 399-415.

Ge D B，Yan Y B. 2005. Finitedifference Time Domain Method for Electromagnetic Waves(in Chinese). Xi′an：Xidian University Press.Giannopoulos A. 2008. An improved new implementation of complex frequency shifted PML for the FDTD method.IEEETransactionsonAntennasandPropagation, 56(9): 2995-3000.Irving J, Knight R. 2006. Numerical modeling of ground-penetrating radar in 2-D using MATLAB.Computers&Geosciences, 32(9): 1247-1258.Katz D S,Thiele E T,Taflove A. 1994. Validation and extension to three dimensions of the Berenger PML absorbing boundary condition for FD-TD meshes.IEEEMicrowaveandGuidedWaveLetters, 4(8): 268-270.

Kuzuoglu M, Mittra R. 1996. Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEEMicrowaveandGuidedWaveLetters, 6(12): 447-449.

Liao Z P, Wong H L, Yang B P, et al. 1984.A transmitting boundary for transient wave analysis.ScientiaSinica(SeriesA),27(10): 1063-1076.Li J, Zeng Z F, Huang L, Liu F S. 2012. GPR simulation based on complex frequency shifted recursive integration PML boundary of 3D high order FDTD.Computers&Geosciences, 49:121-130.

Li J, Zeng Z F, Wu F S, et al. 2010.Study of three dimension high order FDTD simulation of GPR.ChineseJ.Geophys.(inChinese), 53(4)：974-981.Li J X. 2007. Research on algorithms for implementing perfectly matched layers in the finite difference time domain method [Doctor′s thesis]. School of Electronic Information Engineering, Tianjin university.Li Z H, Huang Q H, Wang Y B. 2009. A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media.ChineseJ.Geophys. (in Chinese),52(7): 1915-1922.

Li Z H, Huang Q H. 2014. Application of the complex frequency shifted perfectly matched layer absorbing boundary conditions in transient electromagnetic method modeling.ChineseJ.Geophys. (in Chinese), 57(4): 1292-1299.Martin R, Komatitsch D, Gedney S D, et al. 2010. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML).CMES, 56(1): 17-40.Mei K K, Fang J. 1992. Superabsorption—a method to improve absorbing boundary conditons.IEEETrans.onAntennasandPropagation, 40:1001-1010.Mur G. 1981.Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations.IEEETrans.ElectomagnCompat, EMC-23(4): 377-382.Moerloose J D,Stuchly M. 1995. Behavior of Berenger′s ABC for evanescent waves.IEEEMicrowaveandGuidedWaveLetters, 5(10): 344-346.Roden J, Gedney S D. 2000. Convolution PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media.Microwaveandopticaltechnologyletters, 27: 334-339.

Sacks Z S, Kingsland D M, Lee R, et al. 1995. A perfectly matched anisotropic absorber for use as an absorbing boundary condition.IEEETransactionsonAntennasandPropagation,43 (8): 1460-1463.

Taflove A, Brodwin M E. 1975. Numerieal solution of steady-state EM scattering Problems using the time-dependent Maxwell′s equations.IEEETrans.MicrowaveTheoryTech, 23: 623-630.

Wang L N, Liang C H. 2006. A new implementation of CFS-PML for ADI-FDTD method.MicrowaveandOpticalTechnologyLetters, 48(10): 1924-1928.

Xu T, McMechan G A. 1997. GPR attenuation and its numerical simulation in 2.5 dimensions.Geophysics, 62: 403-414.

Zeng Z F, Liu S X, Feng X. 2010.Theory and Application of Ground Penetrating Radar(in Chinese), Beijing: Electronic Industry Press.

Zhang L X, Fu LY, 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.Zhang W, Shen Y. 2010. Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling.Geophysics, 75(4): 141-154.

Zhang X W, Han L G, Huang L, et al. 2009. A staggered-grid high-order difference method of complex frequency-shifted PML based on recursive integration for elastic wave equation.ChineseJ.Geophys. (in Chinese),52(7): 1800-1807.Zhan Y L, Chang Y J, Cao Z L. 2008. Modeling of groundpenetrating radar based on UPML boundary conditions.ResourcesEnvironmentandEngineering, 22(2): 235-238.

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

底青云，王妙月. 1999. 雷達(dá)波有限元仿真模擬. 地球物理學(xué)報(bào)，42(6)：818-825.

方宏遠(yuǎn)，林皋. 2013. 基于辛算法模擬探地雷達(dá)在復(fù)雜地電模型中的傳播.地球物理學(xué)報(bào)，56(2)：653-659.

馮德山，陳佳維，吳奇. 2014. 混合ADI-FDTD亞網(wǎng)格技術(shù)在探地雷達(dá)頻散媒質(zhì)中的高效正演．地球物理學(xué)報(bào)，57(4)：1322-1334.

馮德山，陳承申，戴前偉. 2010.基于UPML邊界條件的交替方向隱式有限差分法GPR全波場數(shù)值模擬. 地球物理學(xué)報(bào)，53(10)：2484-2496.

馮德山，陳承申，王洪華. 2012. 基于混合邊界條件的有限單元法正演模擬. 地球物理學(xué)報(bào)，55(11): 3774?3785.

馮德山，謝源. 2011. 基于單軸各向異性完全匹配層邊界條件的ADI-FDTD三維GPR全波場正演. 中南大學(xué)學(xué)報(bào)(自然科學(xué)版)，42(8)：2363-2372.

葛德彪，閆玉波. 2005. 電磁波時(shí)域有限差分法. 西安：西安電子科技大學(xué)出版社，67-81.

李建雄. 2007.時(shí)域有限差分法中完全匹配層的實(shí)現(xiàn)算法研究[博士論文]. 天津：天津大學(xué)電子信息工程學(xué)院.

李靜，曾昭發(fā)，吳豐收等. 2010. 探地雷達(dá)三維高階時(shí)域有限差分法模擬研究.地球物理學(xué)報(bào)，53(4)：974- 981.

李展輝，黃清華，王彥賓. 2009. 三維錯格時(shí)域偽譜法在頻散介質(zhì)井中雷達(dá)模擬中的應(yīng)用.地球物理學(xué)報(bào)，52(7)：1915-1922.

李展輝，黃清華. 2014.復(fù)頻率參數(shù)完全匹配層吸收邊界在瞬變電磁法正演中的應(yīng)用. 地球物理學(xué)報(bào)，57(4)：1292-1299.

曾昭發(fā)，劉四新，馮晅. 2010.探地雷達(dá)原理與應(yīng)用. 北京：電子工業(yè)出版社.

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

張顯文，韓立國，黃玲等. 2009.基于遞歸積分的復(fù)頻移PML彈性波方程交錯網(wǎng)格高階差分法. 地球物理學(xué)報(bào)，52(7)：1800-1807.

詹應(yīng)林，昌彥君，曹中林.2008. 基于UPML 吸收邊界的探地雷達(dá)數(shù)值模擬研究.資源環(huán)境與工程，22(2)：235-238.

(本文編輯 胡素芳)

The unsplit convolutional perfectly matched layer absorption performance analysis of evanescent wave in GPR FDTD forward modeling

FENG De-Shan1,3， YANG Liang-Yong1,2， WANG Xun1,3

1SchoolofGeosciencesandInfo-Physics，CentralSouthUniversity，Changsha410083，China2InstituteofGeologyandGeophysics,ChineseAcademyofSciences,Beijing100029,China3KeyLaboratoryofMetallogenicPredictionofNonferrousMetals,MinistryofEducation,Changsha410083，China

This article introduces the principle of CPML boundary conditions, and deduces its FDTD formula for GPR forward modeling, and then compares and analyzes the absorption properties of evanescent wave in Berenger PML, UPML, CPML. The numerical experiment of key parametersκandαin PML were performed so as to determine its range and selection principles. Then, using the two-dimensional TM wave as an example, we studied the mechanism of evanescent wave generation and the influence factors of absorption properties. Wavefield snapshot of uniform medium, reflection error at detection point and global reflection error comparison shows that the three kinds of PML have fabulous absorption capacity to travelling-wave, but quite different to evanescent wave, especially in the low-frequency. Because the parameterαwas introduced, the absorption of evanescent wave in CPML is exceptionally excellent; however, the global error is the largest caused by discretization. Finally, we used the FDTD program with UPML and CPML boundary conditions to carry out 2D GPR forward modeling using profile method and wide-angle method on rectangular geoelectric model and 3D uniform medium model. From the experiment, we observed that there is a huge influence of evanescent wave on both the radar forward modeling profile and wavefield snapshot. Furthermore, taking the absorption performance of evanescent wave by UPML and CPML into consideration, we can determine that CPML can effectively reduce the reflection error in truncation boundary, and obtain a satisfying accuracy. Considering the absorption performance of evanescent wave, the global error, programming difficulty and other factors, in the GPR forward modeling, we highly recommend the use of CPML boundary condition in GPR studies.

Ground Penetrating Radar；Convolutional perfectly matched layer; Uniaxial perfectly matched layer; Evanescent wave; Finite difference time domain method

10.6038/cjg20161232.

國家自然科學(xué)基金資助項(xiàng)目(41574116)，中南大學(xué)創(chuàng)新驅(qū)動項(xiàng)目(2015CX008)，教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃(NCET-12-0551)，中南大學(xué)教師研究基金(2014JSJJ001)，中南大學(xué)升華育英人才計(jì)劃，湖湘青年創(chuàng)新創(chuàng)業(yè)平臺培養(yǎng)對象項(xiàng)目共同資助.

馮德山，男，1978年生，博士，教授，從事地球物理數(shù)據(jù)處理與正反演研究. E-mail：fengdeshan@126.com

10.6038/cjg20161232

P631

2015-11-15，2016-09-12收修定稿

馮德山, 楊良勇, 王珣. 2016. 探地雷達(dá)FDTD數(shù)值模擬中不分裂卷積完全匹配層對倏逝波的吸收效果研究. 地球物理學(xué)報(bào)，59(12):4733-4746,

Feng D S， Yang L Y， Wang X. 2016. The unsplit convolutional perfectly matched layer absorption performance analysis of evanescent wave in GPR FDTD forward modeling.ChineseJ.Geophys. (in Chinese),59(12):4733-4746,doi:10.6038/cjg20161232.