王珣, 朱磊, 馮德山*, 許德如, 丁思元, 劉碩

1 中南大學 地球科學與信息物理學院, 長沙 4100832 東華理工大學 地球科學學院, 南昌 330013

0 引言

探地雷達(GPR)是一種利用高頻電磁波探測地下內(nèi)部結構及地質(zhì)體物性變化的地球物理方法,它被廣泛應用于城市道路病害、管線探測等淺地表勘探領域(Daniels,2004).土壤、混凝土、巖石等地下淺層介質(zhì)大多具有色散特性,其介電常數(shù)是關于頻率的復函數(shù),導致電磁波的相速度和衰減系數(shù)不斷變化,波形產(chǎn)生衰減與畸變(Wang and Oristaglio,2000;Oden et al.,2007).在GPR數(shù)值模擬中,簡單地將地下介質(zhì)設置為非色散介質(zhì),會影響GPR信號的精確解譯.因此,需開展色散介質(zhì)的GPR正演研究(Liu et al.,2007),以便更有效地指導雷達工程探測實踐(劉四新和曾昭發(fā),2007;朱尉強和黃清華,2016).

目前,GPR針對色散介質(zhì)的數(shù)值模擬方法可分為時間域、頻率域兩大類.時間域數(shù)值模擬方法如時域有限差分法(Teixeira,2008;馮德山等,2014;Giannakis and Giannopoulos,2014;李靜等,2010,2016)、時域有限單元法(馮德山等,2013;王洪華和戴前偉,2014;馮德山和王珣,2017;王洪華等,2018;Feng et al.,2018;Liu et al.,2019)、時域偽譜法(Liu and Fan,1999;李展輝等,2009;Huang et al.,2010;Fang and Lin,2012)、時域辛算法(方宏遠和林皋,2013;Fang et al.,2013,2019;雷建偉等,2020;Lei et al.,2022)等被眾多學者應用到Debye、Drude和Cole-Cole等色散介質(zhì)的研究中.時域正演方法應用較為廣泛,但對色散特性的處理會涉及卷積項,需要存儲前一時刻的波場,且針對不同的色散模型需調(diào)整時間離散格式(Pratt and Worthington,1990).頻率域數(shù)值模擬中電位移矢量可表示為電場強度與復相對介電常數(shù)的乘積,對色散項的處理更為簡單直接,同時不受時間推進算法的限制,從而避免了求解所有時刻波場記錄.因此GPR色散介質(zhì)的頻率域正演方法天然地具有原理清晰、易實現(xiàn)的優(yōu)點(Bitri and Grandjean,1998).頻率域正演方法主要有頻域有限差分法(Shin and Fan,2012)與頻域有限單元法(FEFD)(Jin,2015;Feng et al.,2019).頻域有限差分法大都直接離散一階麥克斯韋方程,求解變量多,計算內(nèi)存需求較大.FEFD方法基于變分原理和加權余量法,求解變量較少、占用內(nèi)存小,適應求解復雜電磁問題,因此本文采用FEFD進行色散介質(zhì)的GPR正演.

利用FEFD方法對介質(zhì)進行空間離散時,傳播波場的采樣點數(shù)量不足會導致GPR波動方程的數(shù)值解與真實解之間誤差越來越大,產(chǎn)生數(shù)值頻散(Ihlenburg and Babu?ka,1995;Guddati and Yue,2004).通常采用更細的網(wǎng)格尺寸或更高階的插值基函數(shù)去壓制數(shù)值頻散,但這兩種策略的計算成本較大,降低了正演效率.此外網(wǎng)格尺寸的細化伴隨著一定的數(shù)值偽影,而高階的方法在高頻帶會出現(xiàn)龍格現(xiàn)象(Marfurt,1984).為避免網(wǎng)格加密與高階近似所帶來的問題,優(yōu)化低階有限元的質(zhì)量、剛度矩陣常被用于壓制數(shù)值頻散:傳統(tǒng)方法是將權函數(shù)等同于形函數(shù),該方法下就能得到一致質(zhì)量矩陣(Mullen and Belytschko,1982);或是將單元的質(zhì)量集中在單元節(jié)點上(王月英,2007),得到集中質(zhì)量矩陣(Marfurt,1984).為了獲取更精確的數(shù)值解,Marfurt(1984)以及Christon(1999)提出折衷質(zhì)量矩陣,改善了積分單元的頻散特性.此后,Min等(2003)在FEFD中針對不同剖分單元進行加權平均得到優(yōu)化后的總體剛度、質(zhì)量矩陣;Guddati和Yue(2004)通過改變雙線性矩形單元積分點位置,構造新的質(zhì)量矩陣與剛度矩陣.以上方法大都單獨針對一致矩陣、集中矩陣或二者組合系數(shù)進行研究,或在犧牲計算資源的情況下對有限元的剛度矩陣、質(zhì)量矩陣同時優(yōu)化,難以獲得較高的正演效率與精度.Ogunbo和Shin(2022a,2022b)受優(yōu)化的九點格式有限差分法啟發(fā),采用最速下降法求解出五個參數(shù)組合下的優(yōu)化系數(shù),驗證了優(yōu)化矩陣在不同傳播角度及寬頻帶下,比集中、一致和折衷矩陣具有更優(yōu)越的壓制數(shù)值頻散的性能.

本文在Ogunbo和Shin(2022a,2022b)的研究基礎上,結合有限元質(zhì)量、剛度矩陣自身的約束條件,利用最小二乘方法得到僅需三個參數(shù)組合下的精確優(yōu)化系數(shù),所需的優(yōu)化參數(shù)更少,避免了對源項進行歸一化處理,具有更好的普適性.同時引入基于無限吸收函數(shù)的精確完全匹配層(EPML)(Feng et al.,2019)進一步提高了正演效率.通過頻散曲線、算法收斂性、吸收邊界以及復雜色散模型正演等對比實驗,驗證了本文算法相較于一致、集中和折衷矩陣等常規(guī)FEFD算法在壓制數(shù)值頻散方面更具優(yōu)勢,且能有效提高吸收邊界性能.

1 基于EPML的色散介質(zhì)GPR二維頻率域控制方程

1.1 色散介質(zhì)模型GPR頻率域波動方程

根據(jù)Maxwell方程組,可由時間域標量波動方程推導出GPR二維TM波頻率域傳播方程.假設電磁波在z方向沒有變化,引入時諧因子ejωt得到:

(1)

(2)

代入(1)式得到Helmholtz方程:

(3)

在GPR中常見的色散模型有Cole-Cole模型和Debye模型,許多土壤和巖石可以使用Cole-Cole函數(shù)去模擬(Giannakis et al.,2016).若不考慮電導率的影響,相對介電常數(shù)可表示為(Cole and Cole,1941,1942):

(4)

式中ε∞為介質(zhì)的光頻介電常數(shù),ε0為靜態(tài)介電常數(shù),τ為偶極子的弛豫時間,參數(shù)β(0≤β≤1)表示弛豫時間分散的程度,參數(shù)β=0時,式(4)簡化為Debye色散介質(zhì)(Atteia and Hussein,2010).

1.2 EPML吸收邊界條件

FEFD法應用于GPR正演時,需要對模擬區(qū)域進行人為截斷.目前在截斷邊界處對反射波吸收效果較好的邊界條件為完全匹配層(PML)(Berenger,1994;王洪華等,2019),結合GPR控制方程:

(5)

其中,sx與sy為坐標拉伸變量.常規(guī)PML取得最佳吸收效果的過程本質(zhì)上是求解一個多參數(shù)最優(yōu)化問題,不僅計算繁瑣,且得到的參數(shù)不具有普適性(Gedney,1996).EPML吸收邊界的提出,簡化了邊界條件的加載過程(Bermúdez et al.,2004;Cimpeanu et al.,2015).參考Feng等(2019)有關EPML邊界條件的討論,其sx與sy計算公式為

sx=1-iσx/k,sy=1-iσy/k,

(6)

電導率σx和σy表示EPML中不同方向的吸收函數(shù),σx具體計算公式為

(7)

其中,LPML為PML層厚,dx表示PML內(nèi)部單元中心到模擬區(qū)域邊界的距離.式(7)所表達的吸收函數(shù)形式僅包含最基本的PML厚度及節(jié)點位置,可以極大簡化參數(shù)優(yōu)化過程,提高正演效率.

2 最優(yōu)系數(shù)的頻率域有限元法

2.1 有限元法

根據(jù)Galerkin加權余量法,并代入完美導體邊界條件,得到式(5)的弱解形式:

=?wjωμ0JzdΩ,

(8)

其中w為權函數(shù),利用雙線性插值基函數(shù)對場值Ez進行展開:Ez=∑NiEe,Ni為形函數(shù).單元場向量Ee可表示為

(9)

(10)

其中Ke為單元剛度矩陣,Me為單元質(zhì)量矩陣,fe為單元源向量.其具體表達式:

Me=?Ωesxsyk2NiNjdΩ,

fe=?Ωejωμ0JzNidΩ.

將單元場向量、源向量及單元剛度、質(zhì)量矩陣擴展成整體矩陣,形成多源系統(tǒng)的線性方程組:

AE=f,

(11)

式中A是與頻率相關的所有源共用的系數(shù)矩陣,f是離散化激勵源項,E是對應節(jié)點的離散電場值.

2.2 最優(yōu)參數(shù)選取

數(shù)值頻散的程度取決于剖分單元大小、波的傳播角度以及單位波長內(nèi)的網(wǎng)格點數(shù)(梁昊,2016),選擇合適的參數(shù)可以使數(shù)值頻散最小化.根據(jù)Galerkin法,設四邊形單元剖分下x與y方向上單元邊長為lx與ly,Ke與Me具體表示為

(12)

(13)

對于傳統(tǒng)的有限元方法:a=2,b=1,c=4,d=2,e=1.利用16個四邊形單元對二維無源區(qū)域進行離散,剖分單元邊長為h,如圖1所示.

圖1 計算區(qū)域剖分示意圖

二維頻率域電磁TM波,在xy平面內(nèi)電場設其表達式為

E=E0e-jk(xcosθ+ysinθ),

(14)

其中,θ為傳播方向與x軸的夾角.設波長為λ,單位波長內(nèi)的網(wǎng)格點數(shù)為G,電磁波的傳播速度為v,可知有:λ=2πv/ω,k=ω/v,G=λ/h.

剖分區(qū)域內(nèi)部的節(jié)點周圍都有八個節(jié)點,以圖1中編號13的節(jié)點為例,節(jié)點的離散場值Ei與每個節(jié)點對應系數(shù)Ai滿足關系(Chen et al.,2013):

A7E7+A12E12+A17E17+A8E8+A13E13+A18E18

+A9E9+A14E14+A19E19=0,

(15)

設節(jié)點13處的電場值為u0,式(15)可表示為

(16)

根據(jù)歐拉公式ejx=cosx+jsinx,式(16)可表示為

(17)

由波的傳播規(guī)律可知:kh=2π/G.定義:

(18)

將(18)式代入式(17),由于u0不為零,化簡得:

(19)

此處k是通過有限元方法得到的數(shù)值解波數(shù),數(shù)值解相速度為Vph=ω/k.設真實相速度為V,則歸一化相速度Vph/V可表示為

Vph/V的值越接近1,表示數(shù)值解與真實解的誤差越小.Vph/V與a、b、c、d、e等參數(shù)相關,可通過選取最優(yōu)化的參數(shù)使其接近1,設置目標函數(shù):

J(a,b,c,d,e,G,θ)=

(21)

通常θ與G的取值范圍(Shin and Sohn,1998)分別為:θ∈[0,π/2],G∈[Gmin,Gmax]=[4,400].

觀察參數(shù)a、b、c、d、e在式(21)中分布特征,同比放大縮小一組參數(shù)會存在無數(shù)組優(yōu)化參數(shù),且改變相應的單元系數(shù)矩陣,需要對源向量進行歸一化處理.因此需要對參數(shù)施加約束條件,根據(jù)質(zhì)量守恒定理,令a、b、c、d、e滿足:

a+b=3,c+2d+e=9.

(22)

根據(jù)上式約束條件,目標函數(shù)確定好G與θ后,實際上只有三個獨立參數(shù).若單獨采用一個約束條件a+b=3或c+2d+e=9,則會有四個獨立參數(shù)組合.不同參數(shù)組合下的目標函數(shù)形式不同.以參數(shù)組合ade為例,將式(22)代入式(21)中得到:

6G2(1-P-Q+PQ)a+2π2(2-P-Q)d

+2π2(1-PQ)e=18π2+9G2(2PQ-P-Q),

(23)

相應的矩陣形式為

(24)

其中l(wèi)與r分別為G、θ的離散個數(shù):m=1,2,…,l;n=1,2,…,r.令G、θ的離散值分別為Gn和θm,矩陣中元素S的表達式如下:

本文采用最小二乘算法求解超定方程式(24).結合約束條件(22),可獲得優(yōu)化系數(shù)a、b、c、d、e的值.為了方便對比,分別求出了bde組合、abde組合以及abcde組合下的優(yōu)化系數(shù),結果如表1所示.分析表1可知,ade與bde組合求出的優(yōu)化系數(shù)相同,abde組合下的優(yōu)化系數(shù)值大小與前兩組組合基本相同.五個參數(shù)的優(yōu)化參數(shù)abcde組合存在多解性,且需要對源項權重進行修改才能得到與真實解相同的場值.因此本文采用約束條件式(22)所得到三參數(shù)組合優(yōu)化系數(shù)更加的準確,普適性更好.

表1 不同參數(shù)組合下得到的最優(yōu)化參數(shù)取值

選取ade組合下的最優(yōu)化參數(shù),并得到常規(guī)、集中以及折衷類型(Ogunbo and Shin,2022a)下的a,b,c,d,e不同參數(shù)取值,如表2所示.

表2 不同類型參數(shù)取值

2.3 頻散及精度分析

根據(jù)公式(20)可得到歸一化相速度Vph/V與單位波長內(nèi)網(wǎng)格點數(shù)G的頻散關系如圖2所示.隨著G的減少與傳播角度的改變,歸一化相速度曲線逐漸偏離1,數(shù)值頻散隨之產(chǎn)生.一致、集中矩陣的曲線偏離程度較大,其次是折衷矩陣;而優(yōu)化矩陣的歸一化相速度曲線與1的偏離程度相較于前三種方法極小.歸一化相速度與1的偏差在小于1%時,一致矩陣和集中矩陣的G至少大于12.9,折衷矩陣的G至少大于6.2.而優(yōu)化系數(shù)矩陣的G僅需大于4.8的情況下,其歸一化相速度與1的最大絕對誤差即可小于0.2%.通過頻散曲線分析,優(yōu)化系數(shù)矩陣與一致、集中和折衷矩陣相比,在壓制數(shù)值頻散方面的性能極好,可顯著降低歸一化相速度的誤差.

圖2 (a) 一致、集中、折衷和優(yōu)化矩陣的歸一化相速度曲線; (b) 優(yōu)化矩陣歸一化相速度曲線

值得注意的是,當G為4.8時,優(yōu)化矩陣的四條歸一化相速度曲線會相交于一點.此點的意義十分重要:首先此點對應的歸一化相速度與1的誤差較小,說明壓制頻散的效果較強.其次,此點是不同傳播角度下歸一化相速度曲線相交的點,意味著滿足此點的網(wǎng)格剖分條件下,即使在不同傳播角度,優(yōu)化矩陣都能達到相同的壓制數(shù)值頻散的效果.因此該點對應的G是最優(yōu)化的單位波長下的網(wǎng)格點數(shù).

3 探地雷達有限元色散介質(zhì)模型正演模擬

3.1 網(wǎng)格精度對比

3.1.1 不同網(wǎng)格尺寸精度對比試驗

模擬區(qū)域為10.0 m×10.0 m的Debye色散模型:ε∞=10.0,ε0=15.0,σ=0.001 S·m-1,τ為3 ns,將中心頻率為100 MHz的雷達波脈沖激勵源置于模擬區(qū)域中心.利用四種不同方法對模擬區(qū)域進行頻率域正演,改變正方形剖分單元邊長h,測試在不同網(wǎng)格尺寸下優(yōu)化矩陣算法的收斂性.四種方法在滿足相同剖分條件下的運行時間基本相同,h分別為0.04 m、0.05 m、0.08 m、0.1 m、0.2 m,五種剖分條件下的計算時間分別為0.8 s、0.5 s、0.22 s、0.15 s、0.06 s.

表3為不同方法在不同剖分網(wǎng)格大小下數(shù)值解與解析解的誤差的二范數(shù).同等自由度下,優(yōu)化矩陣下的網(wǎng)格精度相較于集中、一致、折衷矩陣分別提高了89.48%、87.91%、53.57%.圖3是四種方法誤差二范數(shù)與剖分單元邊長h的收斂曲線,隨著h的減小,優(yōu)化矩陣下誤差二范數(shù)總是最小,其次是折衷、一致、集中矩陣.該實驗證明優(yōu)化矩陣在保障正演效率的同時,也能提高正演模擬的精度.

表3 不同網(wǎng)格大小下四種方法的誤差L2范數(shù)

圖3 不同方法下數(shù)值解誤差2范數(shù)收斂曲線

3.1.2 最優(yōu)網(wǎng)格點精度對比實驗

模擬區(qū)域激勵源與介電參數(shù)不變,取單元邊長h=0.186 m,根據(jù)歸一化相速度理論曲線,對應優(yōu)化矩陣最優(yōu)單位波長網(wǎng)格點數(shù)(G=4.8),剖分單元數(shù)量200×200.將源置于模擬區(qū)域中心,接收點坐標(27.714 m,9.486 m).采樣時間為0.1 ns,采樣點數(shù)為3000,加載10層常規(guī)PML,對模擬區(qū)域進行頻率域正演并利用傅里葉逆變換至時間域,得到相應波場快照及單道波形,四種方法的計算時間分別都約為190 s.

圖4為100 ns時刻的波場快照圖,黑色箭頭表示數(shù)值頻散發(fā)生的區(qū)域.其中圖4a黑色實心圓表示激勵源的位置,黑色倒三角表示接收點位置.為方便對比四種方法的頻散特征,本文將四種方法在100 ns時刻的1/4波形拼接為一個波場快照.圖4b表示為四種方法與解析解的殘差對比圖.此時波未到邊界處不會產(chǎn)生反射,且圖4a中的頻散現(xiàn)象屬于正演方法所帶來的數(shù)值頻散而非色散介質(zhì)產(chǎn)生的色散現(xiàn)象.數(shù)值頻散效應在一致矩陣和集中矩陣中比較明顯:一致矩陣方法的頻散效應超前,集中矩陣方法的頻散現(xiàn)象滯后.折衷矩陣的頻散效應既有超前也有滯后,但相比前兩種方法,其頻散程度較小.而優(yōu)化矩陣相較于其他方法未見明顯的數(shù)值頻散現(xiàn)象.

圖4 (a) 在均勻頻散介質(zhì)下100ns時刻不同方法的波場快照圖;(b) 不同方法與解析解的絕對誤差

實際上,上述特征是由歸一化相速度與1的關系所決定的:歸一化相速度大于1時會出現(xiàn)數(shù)值頻散超前的現(xiàn)象,而小于1時會出現(xiàn)頻散滯后的現(xiàn)象.從圖2中可知,一致矩陣和集中矩陣的曲線分別在大于1與小于1的區(qū)域分布,折衷矩陣為前兩者的線性組合,其歸一化相速度曲線在1兩側都有分布;而優(yōu)化矩陣的歸一化相速度曲線與1的誤差極小,數(shù)值頻散現(xiàn)象并不明顯.

圖5a為接收點處四種方法數(shù)值解與解析解的單道波形圖,圖5b則是四種方法對應的單道波形與解析解的絕對誤差.波在135 ns左右到達,一致與集中矩陣方法的頻散效應分別在波形到達前和到達后較為明顯,且二者的波形與解析解相差較大.折衷矩陣方法的單道波形整體與解析解擬合較好,但在小圖以及135 ns前的波形中可以看出,仍然存在一定的頻散現(xiàn)象.紅色虛線所代表的優(yōu)化矩陣方法,無論是整體還是局部,其波形都與解析解最為接近,從誤差圖中可以明顯看出其壓制數(shù)值頻散的優(yōu)勢.

圖5 (a) G=4.8網(wǎng)格剖分下均勻頻散介質(zhì)下不同方法單道波形圖; (b) 不同方法下數(shù)值解與解析解絕對誤差

值得注意的是,從圖5b中可知,一致、集中和折衷矩陣方法的波形在220 ns后還有反射波能量殘留.實際上該反射波是由數(shù)值頻散產(chǎn)生的,可見傳統(tǒng)方法加載10層的常規(guī)PML的吸收效果并不理想,而優(yōu)化矩陣方法不僅能夠壓制數(shù)值頻散,也能適應常規(guī)PML吸收邊界.

3.2 吸收邊界對比

受模擬區(qū)域介電參數(shù)以及剖分單元數(shù)量的限制,正演模擬時難以將G控制在最優(yōu)網(wǎng)格點,需要一種更具普適性的吸收邊界.本節(jié)考慮G不等于最優(yōu)網(wǎng)格點的剖分下,討論EPML與常規(guī)PML的吸收效果對比.選擇中心頻率為100 MHz的雷達波脈沖激勵源,模擬區(qū)域8 m×8 m,剖分單元大小h=0.04 m,模擬區(qū)域介電參數(shù)如表4所示.左側一半設置為Debye介質(zhì)模型,對應G約為22.3;右側一半介質(zhì)設置為Cole-Cole模型,對應G約為20.7.將源置于模擬區(qū)域中心,左右接收點坐標分別為(2.12 m,2.04 m)和(5.96 m,2.04 m).對模擬區(qū)域進行頻率域正演并利用傅里葉逆變換至時間域得到相應波場快照及其單道波形.

表4 模型介電參數(shù)

圖6為PML與EPML吸收邊界在58 ns時刻的吸收效果對比圖,黑色實心圓表示激勵源的位置,黑色倒三角表示接收點位置,黑色矩形框表示模擬區(qū)域邊界,黑色虛線表示介質(zhì)分界面.加載常規(guī)PML層厚度為0.2 m時,波形在邊界處有明顯反射波.0.4 m厚度的常規(guī)PML與0.2 m的EPML下的波形,其反射波能量明顯減弱.從表5可以看出,增加常規(guī)PML的厚度,會犧牲正演運行時間以及內(nèi)存占用;而EPML只需加載0.2 m便可達到0.4 m厚度常規(guī)PML邊界的吸收效果.

表5 不同吸收邊界運行參數(shù)

圖6 (a)—(c) 不同吸收邊界厚度下的波場快照; (d)—(f) 對應吸收邊界殘差對比圖

圖7為不同吸收邊界條件下在接收點處的單道波形.圖中可見加載3種不同PML的直達波都能與參考解較好地吻合.為了分析細微的差別,將圖7中加黑色框線67～148 ns的單道雷達波形進行放大得到小圖,在不同的色散介質(zhì)中,0.2 m的EPML與0.4 m的常規(guī)PML二者吸收效果相近,且都能與參考解較好地擬合;而0.2 m的常規(guī)PML下的反射波能量擾動十分明顯.上述結果表明EPML在能夠在保證吸收效果的同時,也能提高模擬效率.

圖7 R1 (a) 和R2 (b) 接收點處不同吸收邊界厚度下的單道波形

3.3 復雜色散介質(zhì)模型正演

圖8為Debye色散介質(zhì)的城市道路病害模型,模擬區(qū)域2 m×0.8 m,采用邊長0.0125 m的規(guī)則四邊形剖分單元,剖分網(wǎng)格160×60.激勵源主頻900 MHz,放置在空氣層中,激發(fā)接收位置如圖中所示.第一層為隨機雙相介質(zhì)混凝土層(馮德山和王珣,2016),背景介質(zhì)為水泥砂漿,其中骨料散射體占比45%,最大和最小骨料粒徑分別為0.02 m和0.0125 m,該層水平位置0.75 m處存在一個寬度為0.0125 m貫通空氣裂隙.第二層為色散土壤,此層存在一個大小0.2 m×0.1 m的矩形色散異常體,中心位置為(1.6 m,0.3 m).第三層為埋有兩根PVC管的均勻介質(zhì),PVC管壁厚度0.0125 m,內(nèi)徑0.1 m,左側PVC管中心位置為(0.5 m,0.6 m),管內(nèi)充盈液體;右側為PVC空管,中心位置為(1.5 m,0.6 m).模型具體參數(shù)如表6所示.對模擬區(qū)域進行頻率域正演,并采用傅里葉逆變換至時間域,可得到相應不同時刻波場快照及剖面圖.

表6 模型介電參數(shù)表

圖8 城市道路病害模型圖

將激勵源置于空氣層中心,不同算法正演結果如圖9所示.圖9(a、c、e)為加載常規(guī)PML下一致矩陣的FEFD算法下波場快照,圖9(b、d、f)為加載EPML下優(yōu)化矩陣的FEFD算法下的波場快照,二者吸收邊界厚度相同.電磁波在5.2 ns經(jīng)過色散矩形異常體,在6.5 ns時到達PVC空氣管并開始產(chǎn)生反射波,7.8 ns時刻經(jīng)過PVC左管產(chǎn)生反射波.如圖中黑色箭頭所示,圖9(a、c、e)中,一致矩陣下的正演結果出現(xiàn)了超前的數(shù)值頻散,這驗證了前文的分析,而優(yōu)化矩陣下的正演結果未出現(xiàn)明顯的數(shù)值頻散現(xiàn)象.由于加載了EPML吸收邊界,在7.8 ns時刻的優(yōu)化矩陣算法正演結果未出現(xiàn)反射波,而一致矩陣算法波場快照出現(xiàn)了反射波,驗證了EPML良好的吸收性能.

圖9 (a)、(c)、(e) 常規(guī)算法下不同時刻波場快照; (b)、(d)、(e) 優(yōu)化算法下不同時刻波場快照

將發(fā)射天線和接收天線均放在空氣層,采用自激自收的方式,偏移距0.0125 m,每次同步移動0.0125 m,共記錄160道數(shù)據(jù),得到傳統(tǒng)算法與優(yōu)化算法下正演剖面如圖10所示.

圖10 (a) 常規(guī)算法正演剖面圖; (b) 優(yōu)化算法下正演剖面圖

從左右剖面圖都可看出,混凝土層和色散土壤層的邊界是清晰的,而色散土壤層與均勻介質(zhì)層的邊界較為模糊.在混凝土層,由于隨機雙相介質(zhì)中骨料的影響,雷達波散射較嚴重,波形發(fā)生扭曲,出現(xiàn)較多的干擾雜波,同時在該層水平位置0.75 m處可以看到由空氣裂隙產(chǎn)生的繞射波.色散土壤層中可觀察到弧頂在水平位置1.6 m處的矩形異常體的雙曲線反射波.均勻介質(zhì)層中,左右兩側PVC管的雙曲線反射波,其弧頂分別在水平位置0.5 m、1.5 m處.右側的PVC空管的雙曲線反射波并不完整,這是由于空管上方的色散矩形體對雷達波的衰減較強,導致波到達空管處的能量較弱.

一致矩陣算法下的剖面圖10a中黑色箭頭表示數(shù)值頻散較為明顯的地方,且從第二三層界面處可以看出圖中的相位、幅值也發(fā)生了一定變化.從紅色箭頭可以看出,在常規(guī)PML下可以看到明顯的角點反射,而優(yōu)化矩陣下算法由于加載了EPML,在邊緣處未出現(xiàn)明顯的反射波.

4 結論

針對常規(guī)FEFD算法在色散介質(zhì)中存在的數(shù)值頻散問題,提出了一種基于EPML的最優(yōu)系數(shù)FEFD算法.該算法采取歸一化相速度與1的誤差最小策略,通過最小二乘法求解目標函數(shù),實現(xiàn)了僅需三個優(yōu)化參數(shù)組合下對有限元剛度矩陣與質(zhì)量矩陣的優(yōu)化.實驗結果表明:優(yōu)化矩陣的歸一化相速度在單位波長僅需4.8個網(wǎng)格點下便可達到誤差小于0.2%的精度,而傳統(tǒng)方法不僅需要更多的網(wǎng)格點,且精度更低.通過與傳統(tǒng)方法的頻散特征對比分析,驗證了優(yōu)化矩陣在壓制頻散方面優(yōu)良特性.EPML吸收邊界條件的引入,簡化了吸收參數(shù)優(yōu)化過程,5層EPML即可達到10層常規(guī)PML的吸收效果,能夠有效節(jié)約計算成本,提高色散介質(zhì)GPR頻率域正演模擬精度與效率.

致謝衷心感謝編輯以及審稿專家提出的建設性意見,對文章質(zhì)量的提高有很大幫助.