基于非分裂CFS-PML邊界條件的頻散介質(zhì)GPR時(shí)域有限元模擬

王洪華, 吳祺銘, 龔俊波

桂林理工大學(xué)地球科學(xué)學(xué)院, 廣西桂林 541004

0 引言

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

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

(5)

(6)

(7)

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

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

(8)

式中,η表示PML層中的節(jié)點(diǎn)到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)化到時(shí)間域,并整理可得:

(13)

(14)

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

3 時(shí)域有限元法求解的GPR波動(dòng)方程

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

(15)

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

(16)

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

(17)

(18)

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

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

(19)

其通解可表示為

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

(20)

(21)

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

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

(23)

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

(24)

(25)

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

(26)

式(26)整理可得:

(27)

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

(28)

(30)

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

4 數(shù)值計(jì)算

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

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

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

圖2為一個(gè)大小為4.0 m×1.0 m的狹長型均勻頻散介質(zhì)模型,模型中充填Debye頻散介質(zhì)I.利用時(shí)域有限元法進(jìn)行模擬計(jì)算時(shí),計(jì)算區(qū)域被四邊形網(wǎng)格剖分成800×200個(gè)單元,網(wǎng)格單元大小0.005 m×0.005 m,模型內(nèi)PML層厚度為0.15 m.為保證時(shí)間迭代的穩(wěn)定性,時(shí)間步長取為0.01 ns.發(fā)射源是中心頻率為0.6 GHz的雷克子波,位于左上角(0.175 m, 0.175 m),距離左側(cè)和上側(cè)PML層都為5個(gè)網(wǎng)格,如圖中五角星所示;三個(gè)接收點(diǎn)坐標(biāo)位置分別為(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時(shí)域有限元法進(jìn)行計(jì)算的17 ns和28 ns時(shí)刻的波場快照如圖3所示,其中圖3a、3b和3c、3d分別是PML和CFS-PML邊界條件下頻散介質(zhì)中的波場快照;圖3e、3f和3g、3h分別為PML和CFS-PML邊界條件下非頻散介質(zhì)中的波場快照;黑色虛線表示PML層與計(jì)算區(qū)域邊界.由圖可知,由于發(fā)射源靠近左側(cè)和上側(cè)PML層,易產(chǎn)掠射波和低頻反射波,常規(guī)PML邊界條件對掠射波和低頻反射波吸收不完全,在上側(cè)和左側(cè)PML層附近出現(xiàn)較為明顯的掠射波和低頻反射波,如圖3a、3b、3e、3f中白色箭頭所示;且低頻反射波能量隨入射角和傳播距離的增大而增強(qiáng),出現(xiàn)位置延伸到計(jì)算區(qū)域,嚴(yán)重影響波前面的形態(tài).而施加CFS-PML邊界條件后,遠(yuǎn)偏移距處的掠射波和低頻反射波都幾乎被吸收完全,未產(chǎn)生可見的虛假反射,波前面形態(tài)非常完整,如圖3c、3d、3g、3h所示.上述現(xiàn)象也可從圖4所示的PML和CFS-PML邊界條件下地表接收到的GPR剖面中清晰可見,PML邊界條件下,遠(yuǎn)偏移距處的直達(dá)波下方出現(xiàn)較為明顯的低頻虛假反射波,如圖4a和圖4c中PML和CFS-PML邊界條件下的接收點(diǎn)1處接收到單道波形與參考解進(jìn)行對比,分別如圖5a和圖5b所示.參考解是在無任何邊界條件下計(jì)算模型擴(kuò)大6倍,發(fā)射源與接收點(diǎn)相對位置不變,以保證接收點(diǎn)在計(jì)算時(shí)間內(nèi)接收不到模型邊界反射的計(jì)算得到.由圖可知:PML邊界條件下,由于接收點(diǎn)1處接收到了近水平入射產(chǎn)生的掠射波和低頻虛假反射,導(dǎo)致在頻散介質(zhì)中20 ns附近的波形與參考解吻合較差,如局部放大圖所示;與PML邊界條件相比,施加CFS-PML邊界條件后,接收點(diǎn)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時(shí)刻(右側(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)模型中接收點(diǎn)1的單道波形與參考解對比

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

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

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

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

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

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

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

為分析介質(zhì)頻散特性對GPR信號的影響規(guī)律,建立了一個(gè)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的步長等距離移動(dòng)20次.其他計(jì)算參數(shù)與狹長型模型相同.

圖9 層狀GPR模型示意圖

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

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

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

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

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

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

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

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

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

5 結(jié)論

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