最優(yōu)系數(shù)有限單元法色散介質(zhì)GPR頻率域正演模擬

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

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

0 引言

探地雷達(GPR)是一種利用高頻電磁波探測地下內(nèi)部結(jié)構(gòu)及地質(zhì)體物性變化的地球物理方法,它被廣泛應(yīng)用于城市道路病害、管線探測等淺地表勘探領(lǐng)域(Daniels,2004).土壤、混凝土、巖石等地下淺層介質(zhì)大多具有色散特性,其介電常數(shù)是關(guān)于頻率的復(fù)函數(shù),導致電磁波的相速度和衰減系數(shù)不斷變化,波形產(chǎn)生衰減與畸變(Wang and Oristaglio,2000;Oden et al.,2007).在GPR數(shù)值模擬中,簡單地將地下介質(zhì)設(shè)置為非色散介質(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)等被眾多學者應(yīng)用到Debye、Drude和Cole-Cole等色散介質(zhì)的研究中.時域正演方法應(yīng)用較為廣泛,但對色散特性的處理會涉及卷積項,需要存儲前一時刻的波場,且針對不同的色散模型需調(diào)整時間離散格式(Pratt and Worthington,1990).頻率域數(shù)值模擬中電位移矢量可表示為電場強度與復(fù)相對介電常數(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方法基于變分原理和加權(quán)余量法,求解變量較少、占用內(nèi)存小,適應(yīng)求解復(fù)雜電磁問題,因此本文采用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)方法是將權(quán)函數(shù)等同于形函數(shù),該方法下就能得到一致質(zhì)量矩陣(Mullen and Belytschko,1982);或是將單元的質(zhì)量集中在單元節(jié)點上(王月英,2007),得到集中質(zhì)量矩陣(Marfurt,1984).為了獲取更精確的數(shù)值解,Marfurt(1984)以及Christon(1999)提出折衷質(zhì)量矩陣,改善了積分單元的頻散特性.此后,Min等(2003)在FEFD中針對不同剖分單元進行加權(quán)平均得到優(yōu)化后的總體剛度、質(zhì)量矩陣;Guddati和Yue(2004)通過改變雙線性矩形單元積分點位置,構(gòu)造新的質(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)的研究基礎(chǔ)上,結(jié)合有限元質(zhì)量、剛度矩陣自身的約束條件,利用最小二乘方法得到僅需三個參數(shù)組合下的精確優(yōu)化系數(shù),所需的優(yōu)化參數(shù)更少,避免了對源項進行歸一化處理,具有更好的普適性.同時引入基于無限吸收函數(shù)的精確完全匹配層(EPML)(Feng et al.,2019)進一步提高了正演效率.通過頻散曲線、算法收斂性、吸收邊界以及復(fù)雜色散模型正演等對比實驗,驗證了本文算法相較于一致、集中和折衷矩陣等常規(guī)FEFD算法在壓制數(shù)值頻散方面更具優(yōu)勢,且能有效提高吸收邊界性能.

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

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

根據(jù)Maxwell方程組,可由時間域標量波動方程推導出GPR二維TM波頻率域傳播方程.假設(shè)電磁波在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法應(yīng)用于GPR正演時,需要對模擬區(qū)域進行人為截斷.目前在截斷邊界處對反射波吸收效果較好的邊界條件為完全匹配層(PML)(Berenger,1994;王洪華等,2019),結(jié)合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)有關(guān)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加權(quán)余量法,并代入完美導體邊界條件,得到式(5)的弱解形式:

=?wjωμ0JzdΩ,

(8)

其中w為權(quán)函數(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是與頻率相關(guān)的所有源共用的系數(shù)矩陣,f是離散化激勵源項,E是對應(yīng)節(jié)點的離散電場值.

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

數(shù)值頻散的程度取決于剖分單元大小、波的傳播角度以及單位波長內(nèi)的網(wǎng)格點數(shù)(梁昊,2016),選擇合適的參數(shù)可以使數(shù)值頻散最小化.根據(jù)Galerkin法,設(shè)四邊形單元剖分下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)電場設(shè)其表達式為

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

(14)

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

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

A7E7+A12E12+A17E17+A8E8+A13E13+A18E18

+A9E9+A14E14+A19E19=0,

(15)

設(shè)節(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.設(shè)真實相速度為V,則歸一化相速度Vph/V可表示為

Vph/V的值越接近1,表示數(shù)值解與真實解的誤差越小.Vph/V與a、b、c、d、e等參數(shù)相關(guān),可通過選取最優(yōu)化的參數(shù)使其接近1,設(shè)置目標函數(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ù),且改變相應(yīng)的單元系數(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)

相應(yīng)的矩陣形式為

(24)

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

本文采用最小二乘算法求解超定方程式(24).結(jié)合約束條件(22),可獲得優(yōu)化系數(shù)a、b、c、d、e的值.為了方便對比,分別求出了bde組合、abde組合以及abcde組合下的優(yōu)化系數(shù),結(jié)果如表1所示.分析表1可知,ade與bde組合求出的優(yōu)化系數(shù)相同,abde組合下的優(yōu)化系數(shù)值大小與前兩組組合基本相同.五個參數(shù)的優(yōu)化參數(shù)abcde組合存在多解性,且需要對源項權(quán)重進行修改才能得到與真實解相同的場值.因此本文采用約束條件式(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的頻散關(guān)系如圖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)化矩陣的四條歸一化相速度曲線會相交于一點.此點的意義十分重要:首先此點對應(yīng)的歸一化相速度與1的誤差較小,說明壓制頻散的效果較強.其次,此點是不同傳播角度下歸一化相速度曲線相交的點,意味著滿足此點的網(wǎng)格剖分條件下,即使在不同傳播角度,優(yōu)化矩陣都能達到相同的壓制數(shù)值頻散的效果.因此該點對應(yīng)的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īng)優(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ū)域進行頻率域正演并利用傅里葉逆變換至時間域,得到相應(yīng)波場快照及單道波形,四種方法的計算時間分別都約為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ù)值頻散效應(yīng)在一致矩陣和集中矩陣中比較明顯:一致矩陣方法的頻散效應(yīng)超前,集中矩陣方法的頻散現(xiàn)象滯后.折衷矩陣的頻散效應(yīng)既有超前也有滯后,但相比前兩種方法,其頻散程度較小.而優(yōu)化矩陣相較于其他方法未見明顯的數(shù)值頻散現(xiàn)象.

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

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

圖5a為接收點處四種方法數(shù)值解與解析解的單道波形圖,圖5b則是四種方法對應(yīng)的單道波形與解析解的絕對誤差.波在135 ns左右到達,一致與集中矩陣方法的頻散效應(yīng)分別在波形到達前和到達后較為明顯,且二者的波形與解析解相差較大.折衷矩陣方法的單道波形整體與解析解擬合較好,但在小圖以及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ù)值頻散,也能適應(yīng)常規(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所示.左側(cè)一半設(shè)置為Debye介質(zhì)模型,對應(yīng)G約為22.3;右側(cè)一半介質(zhì)設(shè)置為Cole-Cole模型,對應(yīng)G約為20.7.將源置于模擬區(qū)域中心,左右接收點坐標分別為(2.12 m,2.04 m)和(5.96 m,2.04 m).對模擬區(qū)域進行頻率域正演并利用傅里葉逆變換至時間域得到相應(yīng)波場快照及其單道波形.

表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) 對應(yīng)吸收邊界殘差對比圖

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

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

3.3 復(fù)雜色散介質(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,左側(cè)PVC管中心位置為(0.5 m,0.6 m),管內(nèi)充盈液體;右側(cè)為PVC空管,中心位置為(1.5 m,0.6 m).模型具體參數(shù)如表6所示.對模擬區(qū)域進行頻率域正演,并采用傅里葉逆變換至時間域,可得到相應(yīng)不同時刻波場快照及剖面圖.

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

圖8 城市道路病害模型圖

將激勵源置于空氣層中心,不同算法正演結(jié)果如圖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)中,一致矩陣下的正演結(jié)果出現(xiàn)了超前的數(shù)值頻散,這驗證了前文的分析,而優(yōu)化矩陣下的正演結(jié)果未出現(xiàn)明顯的數(shù)值頻散現(xiàn)象.由于加載了EPML吸收邊界,在7.8 ns時刻的優(yōu)化矩陣算法正演結(jié)果未出現(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ì)層中,左右兩側(cè)PVC管的雙曲線反射波,其弧頂分別在水平位置0.5 m、1.5 m處.右側(cè)的PVC空管的雙曲線反射波并不完整,這是由于空管上方的色散矩形體對雷達波的衰減較強,導致波到達空管處的能量較弱.

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

4 結(jié)論

針對常規(guī)FEFD算法在色散介質(zhì)中存在的數(shù)值頻散問題,提出了一種基于EPML的最優(yōu)系數(shù)FEFD算法.該算法采取歸一化相速度與1的誤差最小策略,通過最小二乘法求解目標函數(shù),實現(xiàn)了僅需三個優(yōu)化參數(shù)組合下對有限元剛度矩陣與質(zhì)量矩陣的優(yōu)化.實驗結(jié)果表明:優(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頻率域正演模擬精度與效率.

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