王洪華, 吳祺銘, 龔俊波

桂林理工大學地球科學學院, 廣西桂林 541004

0 引言

探地雷達(Ground Penetrating Radar, GPR)數(shù)值模擬是研究高頻電磁波探測問題的重要手段,在分析復雜介質(zhì)中電磁波傳播規(guī)律和實測資料處理與解釋中均發(fā)揮著重要作用(Uduwawala et al., 2005; 劉四新和曾昭發(fā), 2007; 馮德山和王珣, 2017; Lei et al., 2022).研究表明,實際地下介質(zhì)大都具有頻散特征,相對介電常數(shù)是關(guān)于頻率的復數(shù);電磁波在其中傳播時會發(fā)生畸變(Liu and Fan, 1999; 王洪華和戴前偉, 2014; Lu et al., 2020).因此,在GPR數(shù)值模擬中考慮介質(zhì)的頻散特征可有效匹配實測GPR信號的振幅、頻率、走時、相位等參數(shù),提高實測數(shù)據(jù)偏移、層析成像、全波形反演的精度和可靠性(Sena et al., 2005; Oden et al., 2007; Bitri and Grandjean, 1998; Deparis and Garambois, 2009; Wu et al., 2022).

由于頻率域中電位移矢量可表示為電場強度與復相對介電常數(shù)的乘積,因而在頻率域開展頻散介質(zhì)GPR模擬具有原理簡單和易實現(xiàn)的優(yōu)點(Bitri and Grandjean, 1998; Doyon and Giroux, 2018; Feng et al., 2019),但其涉及到許多頻點的多次重復計算,降低計算效率.雖然時間域中電位移矢量表示為電場強度與復相對介電常數(shù)的卷積,但借助輔助微分方程法可避免卷積計算,而且時間域模擬可利用當前時刻前多個時間步的電磁場來迭代求解當前時刻的電磁場,可顯著提高計算效率.因而,利用時間域數(shù)值模擬方法如時域有限差分法(Teixeira et al., 1998; Bergmann et al., 1998; Cheng et al., 2019; ElMahgoub et al., 2012)、時域偽譜法(Liu and Fan, 1999; Fan et al., 2001; Gao et al., 2004)、時域辛算法(Sha et al., 2007; Yang et al., 2021; Lei et al., 2022)等分析高頻電磁波在Debye、Dude和Cole-Cole頻散介質(zhì)中傳播規(guī)律和異常體響應特征得到廣泛應用(Jung et al., 1999; Petropoulos, 2005; Giannakis et al., 2012).然而,時域有限差分法和時域偽譜法大都直接離散麥克斯韋方程(一階電磁波動方程),求解變量多,計算內(nèi)存需求較大;且Yee網(wǎng)格離散復雜介質(zhì)結(jié)構(gòu),階梯狀網(wǎng)格會產(chǎn)生較大的幾何離散誤差,數(shù)值頻散嚴重,在模擬強非均勻介質(zhì)時會出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象(Gurel and Oguz, 2000; Giannakis et al., 2016; Wei et al., 2017).而時域有限元法大都求解二階電磁波動方程,求解變量較少,降低計算內(nèi)存,提高計算效率;而且可采用非結(jié)構(gòu)化網(wǎng)格離散復雜地質(zhì)結(jié)構(gòu),減少模型離散誤差,提高模擬精度(Di and Wang, 2004; Di et al., 2006; Angulo et al., 2011; 馮德山和王珣, 2016;王洪華等, 2019).目前,時域有限元法及其改進方法被廣泛應用于頻散介質(zhì)GPR模擬,分析高頻電磁波在頻散介質(zhì)中的傳播特征(Lu et al., 2005; Diaz Angulo et al., 2011; Gedney et al., 2012; Liu and Fan, 2019).

計算機內(nèi)存空間有限,利用時域有限元法求解GPR波動方程,人為截斷計算空間時,易在截斷位置附近產(chǎn)生明顯的非物理性反射電磁波(張彬等, 2017),需在截斷位置處設(shè)置無反射邊界條件,以使計算區(qū)域中的電磁波傳播至截斷位置處時能全部被充分吸收,不再返回到計算區(qū)域.目前,吸收性能最佳的無反射邊界條件是完全匹配層(Perfectly Matched Layer, PML)(Berenger, 1994),在非頻散介質(zhì)GPR時域有限元模擬中得到廣泛應用(馮德山和王珣, 2017; 王洪華等, 2019).電磁波在非頻散和頻散介質(zhì)傳播時滿足的方程不同,針對非頻散介質(zhì)模擬提出的PML邊界條件實現(xiàn)方法直接應用于頻散介質(zhì)模擬時,會導致數(shù)值不穩(wěn)定或模擬精度差等問題(Liu and Fan, 1999).而且,PML層內(nèi)滿足的波動方程數(shù)值離散后的反射系數(shù)不嚴格為零,特別是當激勵源靠近邊界,波場以近水平入射時,遠偏移距位置會產(chǎn)生明顯的虛假低頻反射,嚴重降低模擬精度(Komatitsch and Martin, 2007).為此,Kuzuoglu和Mittra(1996)在常規(guī)PML復拉伸算子基礎(chǔ)上,引入復頻移拉伸算子,以改進坐標變換關(guān)系,提出了復頻移完全匹配層邊界條件(Complex Frequency Shifted PML,CFS-PML),以有效吸收掠射波和低頻波(Berenger, 2002; Drossaert and Giannopoulos, 2007).CFS-PML實現(xiàn)方法主要有波場分裂和非分裂法兩種,波場分裂法需引入輔助變量,并構(gòu)建輔助微分方程,數(shù)值實現(xiàn)難度大,占用內(nèi)存較多(Liu et al., 2014; 覃發(fā)兵等, 2019;Liu et al., 2019).為此,在彈性波和電磁波時域有限元模擬的非分裂形式的CFS-PML邊界條件實現(xiàn)方法相繼被提出,如卷積法(Roden and Gedney, 2000)、遞歸積分法(Drossaert and Giannopoulos, 2007)和輔助微分方程法(趙建國等, 2014).卷積法通過遞推公式計算時域PML方程中出現(xiàn)的卷積項,可避免直接計算卷積.遞歸積分法利用遞歸積分離散時域PML方程中出現(xiàn)的卷積項,以提高計算效率(田坤等, 2013).上述兩種方法都屬于線性近似法.輔助微分方程PML通過引入輔助變量來構(gòu)建輔助微分方程,以避免時間域卷積項的計算,具有計算精度高的優(yōu)點(馮德山和王珣, 2017).目前,基于輔助微分方程法的非分裂CFS-PML實現(xiàn)方法在彈性波時域有限元模擬中得到廣泛應用(Martin and Komatitsch, 2009; Basu, 2009; Matzen, 2011; Ma et al., 2019).

本文在上述研究基礎(chǔ)上,為避免波場分裂及卷積計算,通過合理構(gòu)造輔助微分方程,提出了一種基于非分裂CFS-PML邊界條件的Debye頻散介質(zhì)GPR時域有限元模擬算法,并通過2個典型GPR模型的數(shù)值試驗,以分析構(gòu)建的非分裂CFS-PML邊界條件的吸收效果及GPR高頻電磁波在頻散介質(zhì)中的傳播特征.

1 Debye頻散介質(zhì)中的二階電磁波動方程

假定地電模型的走向為y軸,則無源條件下的二階時間域電磁波動方程為(王洪華等, 2019):

(1)

其中,Ey為y方向上的電場強度;μ,σ,ε分別是介質(zhì)的磁導率、電導率和復介電常數(shù);t為時間.

式(1)兩邊進行傅里葉變換,整理可得(Bitri and Grandjean, 1998):

(2)

假定介質(zhì)的復介電常數(shù)隨頻率變化滿足Debye頻散關(guān)系,可表示為(劉四新和曾昭發(fā), 2007)

(3)

其中,ε0和ε∞分別表示頻率趨于0和無窮大的介電常數(shù);τ表示馳豫時間.

將式(3)代入式(2)并整理,可得Debye頻散介質(zhì)滿足的GPR二階電磁波動方程為

(4)

2 基于輔助微分方程的非分裂CFS-PML邊界條件實現(xiàn)方法

PML邊界條件的基本思想是利用復拉伸坐標變換重構(gòu)GPR波動方程,以使電磁波能量在PML層呈指數(shù)衰減(Berenger, 1994).PML層內(nèi)的復拉伸坐標變換公式可表示為(Chew and Liu, 1996)

(5)

(6)

(7)

其中,kη≥1和αη≥0分別是η方向上的收縮因子和頻移因子.收縮因子的作用是將大角度入射波矯正為垂直入射波,進入PML層能夠得到充分衰減(Drossaert and Giannopoulos, 2007).頻移因子的引入可使得復平面的極點從實軸移動到虛負半平面,從而有效吸收切入射波引起的低頻波和隱失波(趙建國等, 2014).

實際計算中,復拉伸坐標函數(shù)(5)式中的各參數(shù)一般要求沿PML區(qū)域和求解區(qū)域的交界面外法線方向漸進變化,其表達式可設(shè)置為(Collino and Tsogka, 2001)

(8)

式中,η表示PML層中的節(jié)點到PML層內(nèi)界面的距離;L為PML層厚度;m為調(diào)解因子,dmax,kmax,αmax為常數(shù).

將式(5)代入式(4),可得:

(9)

式(9)兩邊同乘以sxsz,并整理可得:

(10)

(11)

將式(7)代入式(11)并整理,可得:

(12)

利用傅里葉反變換將式(12)轉(zhuǎn)化到時間域,并整理可得:

(13)

(14)

當收縮因子k=1,α=0時,式(13)退化為常規(guī)PML層內(nèi)滿足的電磁波方程.可見,式(13)為施加不同PML邊界條件的統(tǒng)一電磁波方程形式,在基于常規(guī)PML邊界條件的GPR數(shù)值模擬程序上,只需增加描述CFS-PML邊界條件中收縮因子及衰減因子數(shù)組和相應的系數(shù)矩陣即可實現(xiàn)CFS-PML邊界條件的加載,有易于推廣應用.

3 時域有限元法求解的GPR波動方程

假定Ω和Γ分別為計算模型及其PML外邊界;nx和nz分別為計算模型外邊界外法向量的x和z分量;φ為電場Ey的變分試函數(shù);利用Galerkin法可推導式(13)的積分弱形式為(徐世浙, 1994; 劉有山等, 2013):

(15)

式(15)形成的時域有限元方程為

(16)

其中,系數(shù)矩陣的計算公式為

(17)

(18)

其中,N是形函數(shù),Nx和Nz分別是形函數(shù)關(guān)于x和z的導數(shù).本文采用四變形單元和線性基函數(shù)進行空間離散.

由于Newmark-β差分法相比于中心差分法計算精度和穩(wěn)定性更好(Newmark, 1959),因此本文采用Newmark-β差分法對式(16)進行時間離散.式(16)中各式的一般形式為

(19)

其通解可表示為

(M+γΔtC+βΔt2K)Ut+Δt=

(20)

(21)

將式(21)的通解代入式(16)第一式,并整理可得:

與一般式不同的是,式(16)中其余四式中右端項同時包含場量關(guān)于時間的一階和二階導數(shù)項,無法直接進行Newmark-β差分.為此,引入中間變量ω,并分別等于方程左端和右端項(張玉強,2017),則式(16)第二式可表示為

(23)

式(23)中的兩式分別以相同的Δt進行Newmark-β差分離散,可得:

(24)

(25)

式(24)和式(25)中右端項相同,則其左端項也應相同,可得:

(26)

式(26)整理可得:

(27)

同理,式(16)第三、四、五式的Newmark時間差分離散可分別表示為

(28)

(30)

式(22)、(27)、(28)、(29)、(30)為式(16)的Newmark時間差分迭代公式,程序?qū)崿F(xiàn)過程中,依次利用式(22)、(27)、(29)、(30)、(28)分別計算Qt+Δt,Rt+Δt,Px,t+Δt,Pz,t+Δt,Ey,t+Δt,從而實現(xiàn)Debye頻散介質(zhì)中GPR高頻電磁波場的迭代計算.

4 數(shù)值計算

為測試本文構(gòu)建的Debye頻散介質(zhì)GPR時域有限元模擬的非分裂CFS-PML邊界條件的正確性和有效性,建立了2種Debye頻散介質(zhì)I和II,其中介質(zhì)I的復相對介電常數(shù)為6+4i,電導率為0.005 S·m-1,弛豫時間為0.1 ns;介質(zhì)II的復相對介電常數(shù)為10+8i,電導率為0.01 S·m-1,弛豫時間為0.2 ns.圖1展示了頻散介質(zhì)I和頻散介質(zhì)II的復相對介電常數(shù)隨頻率變化曲線.由圖可見,在0.1～10 GHz范圍內(nèi),頻散介質(zhì)的復相對介電常數(shù)隨頻率變化而變化;頻率趨于零或無窮大時,復相對介電常數(shù)都趨于穩(wěn)定.

圖1 Debye頻散介質(zhì)I (a)和介質(zhì)II (b)的復相關(guān)介電常數(shù)隨頻率變化曲線

4.1 非分裂CFS-PML邊界條件的吸收效果分析

圖2為一個大小為4.0 m×1.0 m的狹長型均勻頻散介質(zhì)模型,模型中充填Debye頻散介質(zhì)I.利用時域有限元法進行模擬計算時,計算區(qū)域被四邊形網(wǎng)格剖分成800×200個單元,網(wǎng)格單元大小0.005 m×0.005 m,模型內(nèi)PML層厚度為0.15 m.為保證時間迭代的穩(wěn)定性,時間步長取為0.01 ns.發(fā)射源是中心頻率為0.6 GHz的雷克子波,位于左上角(0.175 m, 0.175 m),距離左側(cè)和上側(cè)PML層都為5個網(wǎng)格,如圖中五角星所示;三個接收點坐標位置分別為(2.175 m, 0.175 m)、(2.175 m, 0.525 m)、(0.275 m, 0.825 m),分別靠近上側(cè)PML層、模型中間和下側(cè)PML層,如圖中黑色圓圈所示.

圖2 狹長型均勻頻散介質(zhì)模型示意圖

利用基于非分裂PML和CFS-PML邊界條件的Debye頻散介質(zhì)與非頻散介質(zhì)GPR時域有限元法進行計算的17 ns和28 ns時刻的波場快照如圖3所示,其中圖3a、3b和3c、3d分別是PML和CFS-PML邊界條件下頻散介質(zhì)中的波場快照;圖3e、3f和3g、3h分別為PML和CFS-PML邊界條件下非頻散介質(zhì)中的波場快照;黑色虛線表示PML層與計算區(qū)域邊界.由圖可知,由于發(fā)射源靠近左側(cè)和上側(cè)PML層,易產(chǎn)掠射波和低頻反射波,常規(guī)PML邊界條件對掠射波和低頻反射波吸收不完全,在上側(cè)和左側(cè)PML層附近出現(xiàn)較為明顯的掠射波和低頻反射波,如圖3a、3b、3e、3f中白色箭頭所示;且低頻反射波能量隨入射角和傳播距離的增大而增強,出現(xiàn)位置延伸到計算區(qū)域,嚴重影響波前面的形態(tài).而施加CFS-PML邊界條件后,遠偏移距處的掠射波和低頻反射波都幾乎被吸收完全,未產(chǎn)生可見的虛假反射,波前面形態(tài)非常完整,如圖3c、3d、3g、3h所示.上述現(xiàn)象也可從圖4所示的PML和CFS-PML邊界條件下地表接收到的GPR剖面中清晰可見,PML邊界條件下,遠偏移距處的直達波下方出現(xiàn)較為明顯的低頻虛假反射波,如圖4a和圖4c中PML和CFS-PML邊界條件下的接收點1處接收到單道波形與參考解進行對比,分別如圖5a和圖5b所示.參考解是在無任何邊界條件下計算模型擴大6倍,發(fā)射源與接收點相對位置不變,以保證接收點在計算時間內(nèi)接收不到模型邊界反射的計算得到.由圖可知:PML邊界條件下,由于接收點1處接收到了近水平入射產(chǎn)生的掠射波和低頻虛假反射,導致在頻散介質(zhì)中20 ns附近的波形與參考解吻合較差,如局部放大圖所示;與PML邊界條件相比,施加CFS-PML邊界條件后,接收點1處接收到的單道波形與參考解高度吻合,誤差非常微小.對比圖5a和5b可知,PML和CFS-PML邊界條件下頻散介質(zhì)中產(chǎn)生的振幅差異分別約為6.32×10-5V·m-1和3.70×10-6V·m-1;而PML和CFS-PML邊界條件下非頻散介質(zhì)中產(chǎn)生的振幅差異約為1.10×10-3V·m-1和8.57×10-5V·m-1;可見,頻散介質(zhì)中產(chǎn)生的低頻反射波振幅要比非頻散介質(zhì)更弱.

圖3 PML和CFS-PML邊界條件下均勻頻散和非頻散介質(zhì)模型中17 ns(左側(cè))和28 ns時刻(右側(cè))的波場快照(a)、(b)和(c)、(d)分別為PML和CFS-PML邊界條件下的頻散介質(zhì)模型中的波場快照;(e)、(f)和(g)、(h)分別為PML和CFS-PML邊界條件下非頻散介質(zhì)模型中的波場快照.

圖4 PML邊界條件和CFS-PML邊界條件下地表接收到的GPR剖面(a)和(b)分別為PML和CFS-PML邊界條件下頻散介質(zhì)模型中地表接收到的GPR剖面; (c)和(d)分別為PML和CFS-PML邊界條件下非頻散介質(zhì)模型中地表接收到的GPR剖面.

圖5 PML和CFS-PML邊界條件下頻散介質(zhì)(a)和非頻散介質(zhì)(b)模型中接收點1的單道波形與參考解對比

圖6與圖7為接收點2和接收點3處的單道波形與參考解對比.由于發(fā)射源位于模型左上角,接收點2和接收點3只能接收到近垂直入射和小角度入射波,因而PML和CFS-PML邊界條件下的單道波形均與參考解擬合較好.由此可見,PML和CFS-PML邊界條件都可實現(xiàn)外行波的有效吸收,CFS-PML邊界條件對近垂直入射和小角度入射波的吸收效果并無明顯改善.

圖6 PML和CFS-PML邊界條件下頻散介質(zhì)(a)和非頻散介質(zhì)(b)模型中接收點2的單道波形與參考解對比

圖7 PML和CFS-PML邊界條件下頻散介質(zhì)(a)和非頻散介質(zhì)(b)模型中接收點3的單道波形與參考解對比

圖8為PML和CFS-PML邊界條件下的頻散介質(zhì)和非頻散介質(zhì)模型(不包含PML層)內(nèi)波場能量隨時間衰減曲線.由圖可見:當發(fā)射源開始激發(fā)時,在約0～2.5 ns之間,能量被注入計算模型中,能量急劇增加;2.5 ns之后,由于波場傳播左側(cè)和上側(cè)PML層,部分能量被PML邊界條件吸收,開始出現(xiàn)緩慢的下降趨勢;電磁波在非頻散介質(zhì)和頻散介質(zhì)分別傳播14 ns 和16 ns后,到達下側(cè)PML層,能量進一步被吸收;與總能量相比下降了約7個數(shù)量級.與PML邊界條件相比,CFS-PML邊界條件下的波場能量在25 ns以后更被有效吸收;與非的黑色箭頭所示;而圖4b和4d所示的CFS-PML邊界條件下的GPR剖面中的直達波下方未出現(xiàn)可見的低頻反射波,直達波更為清晰可見.分別對比圖3b和3f、圖4a和4c可知,由于高頻電磁波在非頻散介質(zhì)中傳播時要比在頻散介質(zhì)中衰減更弱,因而在非頻散介質(zhì)中遠偏移距處產(chǎn)生的低頻反射波能量更強.

圖8 PML和CFS-PML邊界條件下頻散介質(zhì)和非頻散介質(zhì)模型內(nèi)波場能量衰減曲線

頻散介質(zhì)相比,波場在頻散介質(zhì)中傳播時衰減更快.由上述分析結(jié)果可知,CFS-PML邊界條件的吸收效果要優(yōu)于常規(guī)PML邊界.

4.2 高頻電磁波在頻散介質(zhì)中的傳播特征分析

為分析介質(zhì)頻散特性對GPR信號的影響規(guī)律,建立了一個2.5 m×2.0 m的層狀介質(zhì)模型,如圖9所示.上層介質(zhì)為頻散介質(zhì)I,厚度為1.0 m;下層介質(zhì)為頻散介質(zhì)II.GPR測量采用共發(fā)射源測量方式,發(fā)射天線T和接收天線R都位于地表,其中發(fā)射天線固定于0.1 m位置處,而接收天線首先位于0.3 m,并以0.1 m的步長等距離移動20次.其他計算參數(shù)與狹長型模型相同.

圖9 層狀GPR模型示意圖

圖10分別為頻散層狀介質(zhì)和非頻散層狀介質(zhì)的GPR正演剖面.其中,頻散層狀介質(zhì)的相對介電常數(shù)為圖1所示的頻散介質(zhì)I和頻散介質(zhì)II的復相對介電常數(shù)隨頻率變化曲線,非頻散介質(zhì)的相對介電常數(shù)為中心頻率為600 MHz時圖1所對應的相對介電常數(shù)實部,分別為5.75和9.27,電導率與頻散介質(zhì)中的電導率相同.為了更好地顯示遠偏移距處的反射波特征,每道GPR信號都利用直達波的最大振幅值進行了歸一化.圖10b所示的非頻散介質(zhì)中的GPR子波幾乎相同,而圖10a中頻散介質(zhì)的GPR子波被拉寬,分辨率降低;且子波寬度隨偏移距增大而增大,這是由于電磁波在頻散介質(zhì)中傳播時,高頻成分衰減更快所致.

圖10 層狀模型的GPR正演剖面(a) 頻散介質(zhì); (b) 非頻散介質(zhì).

利用不同中心頻率的雷克子波作為發(fā)射源進行計算獲得的水平位置1.5 m處的單道波形對比如圖11所示,由圖可知,發(fā)射源中心頻率越大,1.5 m處接收到的波形振幅越弱,特別是當中心頻率為600 MHz和800 MHz時,反射波幾乎衰減完全.由此可見中心頻率越高,在頻散介質(zhì)中傳播時衰減更強,反射波更不易被識別.

圖11 不同中心頻率的雷克子波模擬的1.5 m處的波形對比

為分析介質(zhì)頻散對GPR信號振幅和傳播速度的影響,提取圖10中頻散和非頻散介質(zhì)的直達波和反射波的振幅隨偏移距變化曲線如圖12所示.由圖可知:頻散和非頻散介質(zhì)中的直達波和反射波的能量都隨偏移距的增大而減小;相比非頻散介質(zhì),頻散介質(zhì)中的直達波和反射波衰減更快;相比反射波能量的衰減,不管是頻散介質(zhì)還是非頻散介質(zhì),直達波衰減的更快.圖13為利用圖10中GPR數(shù)據(jù)計算得到的速度譜,由圖所示:相比圖13b中非頻散介質(zhì)的速度譜,圖13a中頻散介質(zhì)的速度譜中的能量更發(fā)散,通過提取速度譜中能量最強位置對應的速度,估計的頻散介質(zhì)中的電磁波速度為0.117 m·ns-1,而非頻散介質(zhì)中的電磁波傳播速度為0.123 m·ns-1,與真實值0.1251 m·ns-1相比誤差更小.由此可見,電磁波在頻散介質(zhì)中的傳播速度要比非頻散介質(zhì)更低.

圖12 圖10中直達波和反射波信號的振幅隨偏移距變化曲線

圖13 圖10中的GPR數(shù)據(jù)計算得到的速度譜(a) 頻散介質(zhì); (b) 非頻散介質(zhì).

圖14a和圖14b分別為圖10a所示的頻散介質(zhì)中水平位置0.8 m和1.4 m的單道波形,圖14c和圖14d為利用短時傅里葉變換計算的時頻分析剖面.圖14e和14f分別為圖10b所示的非頻散介質(zhì)中水平位置0.8 m和1.4 m的單道波形,圖14g和圖14h為相應的時頻分析剖面.從圖14c提取的直達波和反射波的主頻分別為390 MHz和48 MHz,與子波的中心頻率600 MHz相比,高頻成分衰減特別嚴重;且當傳播到1.4 m時,直達波和反射波的主頻降低為268 MHz和8 MHz.而圖14g和圖14h所示的非頻散介質(zhì)的時頻分析剖面中,直達波和反射波主頻均為561 MHz,只有少量高頻成分得到衰減.可見,與電磁波在非頻散介質(zhì)中傳播相比,電磁波在頻散介質(zhì)中傳播時高頻成分衰減更強,直達波和反射波的主頻更小,且主頻隨傳播距離的增大而減小.上述分析結(jié)果表明:相比于非頻散介質(zhì),高頻電磁波在頻散介質(zhì)中傳播衰減更強,子波持續(xù)時間增大,分辨率降低,傳播速度降低,直達波和反射波的主頻更小.

圖14 圖10中水平位置0.8 m和1.4 m處的單道波形及其時頻分析剖面 (a)、(e)和(c)、(g)分別為頻散介質(zhì)和非頻散介質(zhì)中水平位置0.8 m處的單道波形及其時頻分析剖面; (b)、(f)和(d)、(h)分別為頻散介質(zhì)和非頻散介質(zhì)中水平位置1.4 m處的單道波形及其時頻分析剖面.

5 結(jié)論

本文通過合理構(gòu)造輔助微分方程,推導了Debye頻散介質(zhì)GPR時域有限元模擬的非分裂CFS-PML邊界條件實現(xiàn)公式,避免了電磁波場分裂和卷積計算.在此基礎(chǔ)上,用Galerkin法和Newmark-β差分法推導了相應的時域有限元方程及其時間差分離散格式,構(gòu)建一種基于非分裂CFS-PML邊界條件的Debye頻散介質(zhì)GPR時域有限元模擬算法.兩個GPR模型的模擬結(jié)果表明:本文構(gòu)建的基于輔助微分方程的非分裂CFS-PML邊界條件實現(xiàn)方法可有效地吸收大角度入射的低頻虛假反射波,提高模擬精度;相比于非頻散介質(zhì),高頻電磁波在頻散介質(zhì)中傳播衰減更強、子波持續(xù)時間增大、分辨率和傳播速度降低、直達波和反射波的主頻更小,分析結(jié)果有助于提高實測資料的解譯精度.