馮德山, 劉碩, 王珣*, 丁思元, 張華, 蘇玄, 陳磊, 顏照坤

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院, 長(zhǎng)沙 410083

2 東華理工大學(xué)地球科學(xué)學(xué)院, 南昌 330032

0 引言

探地雷達(dá)(Ground penetrating radar,GPR)是一種利用高頻電磁波查明地下界面或者地質(zhì)體空間位置、結(jié)構(gòu)的勘探方法.由于它操作簡(jiǎn)單、效率快、結(jié)果直觀等優(yōu)點(diǎn),現(xiàn)已被應(yīng)用于多個(gè)領(lǐng)域,例如：民用基礎(chǔ)設(shè)施檢測(cè)、地震災(zāi)害生命探測(cè)與救援等(劉瀾波和錢榮毅,2015;胡群芳等,2020).但因地下介質(zhì)較為復(fù)雜,又存在各種自然干擾和人為噪聲,導(dǎo)致雷達(dá)剖面各種波形混雜,難以區(qū)分.因此需要開展GPR正演模擬,了解雷達(dá)波的傳播過(guò)程和特征,提高探地雷達(dá)的解譯精度,為逆時(shí)偏移和反演奠定理論基礎(chǔ)(劉四新和曾昭發(fā),2007;朱尉強(qiáng)和黃清華,2016;葛德彪和魏兵,2019).

探地雷達(dá)數(shù)值模擬方法很多,包括時(shí)域有限差分法(finite-difference time-domain,FDTD)(Yee,1966;李靜等,2016)、時(shí)域有限單元法(finite element time-domain,FETD)(馮德山和王珣,2017;王洪華等,2018)、時(shí)域間斷Galerkin算法(discontinuous Galerkin time-domain,DGTD)(K?nig et al.,2010)、時(shí)域偽譜法(pseudo spectral time-domain,PSTD)(李展輝等,2009;Huang et al.,2010;Fang and Lin,2012)、有限體積法(finite-volume time-domain,FVTD)(Munz et al.,2000)、辛算法(方宏遠(yuǎn)和林皋,2013;雷建偉等,2020)等,這些算法各有特色.FDTD算法雖簡(jiǎn)單易懂,計(jì)算速度快,但頻散較嚴(yán)重,通常用于直交網(wǎng)格,對(duì)于幾何不規(guī)則的地質(zhì)體,并不能緊密貼合其邊界(葛德彪和閆玉波,2005;李靜等,2010).FETD算法能利用非規(guī)則的三角形對(duì)模擬區(qū)域進(jìn)行剖分,可處理地表起伏問(wèn)題,但需要求解大型矩陣,計(jì)算效率低(馮德山等,2013;Jin,2014;王洪華等,2019).PSTD算法計(jì)算精度高,占用內(nèi)存相對(duì)較小,但不適用于復(fù)雜的地質(zhì)體模型,不利于并行計(jì)算.FVTD算法本身包含幾何信息,易處理復(fù)雜網(wǎng)格問(wèn)題,但算法本身較為復(fù)雜,不易提高計(jì)算精度.而DGTD算法結(jié)合了FETD與FVTD算法的優(yōu)點(diǎn),可與非規(guī)則網(wǎng)格結(jié)合,易于實(shí)現(xiàn)網(wǎng)格自適應(yīng)和并行,計(jì)算精度高,受到眾多研究學(xué)者的青睞(Hesthaven and Warburton,2002).1973年,Reed和Hill(1973)在微分三角形離散坐標(biāo)方程中提出了DG的思想.Hu等(1999)在二維波動(dòng)方程的半離散背景下研究了DGTD算法的色散和耗散特性.Lu等(2005)推導(dǎo)了Debye線性色散介質(zhì)和完全匹配層區(qū)域中Maxwell方程的統(tǒng)一公式.Diaz Angulo等(2011)將DGTD算法應(yīng)用于GPR三維模擬中,并與FDTD算法進(jìn)行對(duì)比.Yang等(2017)采用了非結(jié)構(gòu)化四面體的正交基函數(shù),求解了三維非均勻介質(zhì)的DGTD迭代公式,說(shuō)明該算法具有較高的精度和寬頻模擬能力.這些論文都證明了DGTD算法具有精度高、收斂性高的優(yōu)點(diǎn),表明了DGTD算法在GPR正演模擬中的優(yōu)勢(shì).

隨著DGTD算法在計(jì)算電磁學(xué)領(lǐng)域的深入,很多學(xué)者也對(duì)DGTD算法進(jìn)行了各種改進(jìn).DGTD算法最大的一個(gè)優(yōu)勢(shì)是易于實(shí)現(xiàn)并行計(jì)算,作者在之前的論文中已實(shí)現(xiàn)GPU-DGTD算法,在此不再進(jìn)行贅述(Feng et al.,2022).劉梅林和劉少斌(2008)將高階Runge-Kutta DGTD方法應(yīng)用于電磁場(chǎng)諧振腔的研究中.Xu等(2020)從數(shù)學(xué)角度研究了常系數(shù)雙曲線性方程的RKDG方法,詳細(xì)研究了能量方程的矩陣傳遞過(guò)程,給出了在標(biāo)準(zhǔn)CFL條件下保證L2范數(shù)穩(wěn)定的充分條件.Xu和Zhang(2022)討論了用基于偏迎風(fēng)數(shù)值通量的DG方法求解一維變系數(shù)雙曲方程,分別給出了其在半離散和全離散情況下的穩(wěn)定性分析和hp誤差估計(jì).李林茜等(2016)從二維TM情形弱解方程出發(fā),討論了當(dāng)前三角單元和相鄰單元進(jìn)行場(chǎng)量交互時(shí)數(shù)值通量的物理意義和不同形式.雖然眾多學(xué)者在通量、時(shí)間離散方面進(jìn)行了不同程度的研究,但對(duì)于探地雷達(dá)正演模擬,關(guān)于DGTD算法計(jì)算效率和精度的影響因素分析,仍有待進(jìn)一步深入研究.

基于上述DGTD算法的研究現(xiàn)狀和其他學(xué)者的研究基礎(chǔ),本文將DGTD算法應(yīng)用到了探地雷達(dá)正演模擬中.從數(shù)值通量、時(shí)間離散格式、單元大小與基函數(shù)階次、網(wǎng)格剖分方式4個(gè)方面,分析了它們對(duì)于DGTD算法的計(jì)算精度和效率的影響,并利用火星烏托邦平原模型證明了DGTD算法的有效性和實(shí)用性.

1 DGTD算法求解GPR波動(dòng)方程

探地雷達(dá)在地下的傳播符合電磁波理論,即Maxwell方程組.二維時(shí)間域TM波的標(biāo)量方程為(Taflove and Brodwin,1975)

(1)

其中,ε、σ和μ為介質(zhì)的物性參數(shù),分別為介電常數(shù)、電導(dǎo)率和磁導(dǎo)率;t為電磁場(chǎng)的傳播時(shí)間;Hx、Hy和Ez分別表示x、y、z方向上的電磁場(chǎng)分量;Jz為加載在z方向的電流源.

為了簡(jiǎn)化方便推導(dǎo),將(1)式化簡(jiǎn)為(Lu et al.,2004)

(2)

-?Ωm(D(Um)φ)dΩm=-∮?ΩmD(Um)φ·nmdl,(3)

如果不對(duì)局部解或者基函數(shù)加任何條件,那么對(duì)所有單元均成立的(3)式將會(huì)在單元邊界上具有多解,因此,必須引入數(shù)值通量統(tǒng)一兩個(gè)相鄰單元的邊界值,記為Um*.則(3)式變?yōu)?/p>

-∮?ΩmD(Um*)φ·nmdl.

(4)

將(4)式的左邊再做一次分部積分,可得到其強(qiáng)格式方程(Hesthaven and Warburton,2008)：

∮?Ωm[D(Um)-D(Um*)]φ·nmdl,

(5)

(6)

(7)

利用Runge-Kutta(RK)算法進(jìn)行時(shí)間離散.為了方便理解和推導(dǎo),進(jìn)行如下變換：

(8)

采用s級(jí)顯式RK格式進(jìn)行離散,在時(shí)間步長(zhǎng)Δt內(nèi)可以表示為(Butcher,1976)

(9)

2 DGTD 算法影響因素

利用DGTD算法進(jìn)行GPR正演時(shí),影響其正演速度和精度的因素有很多.本文中作者選擇了幾個(gè)相對(duì)重要的因素展開分析與討論,主要包括：數(shù)值通量、時(shí)間離散格式、單元大小與基函數(shù)階次、網(wǎng)格剖分方式.

2.1 數(shù)值通量

數(shù)值通量具有很多種形式,為了GPR正演模擬的高精度,并且在計(jì)算過(guò)程中保持穩(wěn)定性和收斂性,應(yīng)當(dāng)首先選擇具有守恒性、單調(diào)性和相容性的數(shù)值通量.局部Lax-Friedrichs數(shù)值通量簡(jiǎn)單易懂,便于推導(dǎo),是目前最快速、有效的方法之一.在第m個(gè)三角單元的邊界上,二維TM模式下局部Lax-Friedrichs數(shù)值通量的表達(dá)式為

(10)

(11)

其中,εm為當(dāng)前單元的介電常數(shù);μm為當(dāng)前單元的磁導(dǎo)率.

表1 DGTD中局部Lax-Friedrichs數(shù)值通量的三種形式Table 1 Three forms of local Lax-Friedrichs numerical fluxes in DGTD

2.2 時(shí)間離散

對(duì)顯式的RK方案來(lái)說(shuō),級(jí)數(shù)s與階數(shù)p的關(guān)系為p≤s.當(dāng)p≤4時(shí),p=s是可能存在的;而p≥5時(shí),p級(jí)p階的RK方案則不存在.這表明如若RK方案在5階以上,則需要更多的中間變量和步驟來(lái)進(jìn)行計(jì)算,這無(wú)疑會(huì)增大計(jì)算量,并且需要大型的儲(chǔ)存空間.因此,本文為了在保證計(jì)算精度的同時(shí),盡量提高計(jì)算速率,選擇了4階RK方案來(lái)進(jìn)行時(shí)間上的離散,主要有以下三種(Hesthaven and Warburton,2008)：

(1)標(biāo)準(zhǔn)四級(jí)四階顯式RK方法(explicit RK,ERK)

(12)

其穩(wěn)定性條件為

(13)

其中CN表示當(dāng)前網(wǎng)格下,N階DGTD算法的穩(wěn)定性條件.

(2)五級(jí)四階保持強(qiáng)穩(wěn)定的顯式RK方法(strong-stability-preserving RK,SSPRK)

(14)

其中,ai,bi,ci,di均為SSPRK系數(shù)(Hesthaven and Warburton, 2008).從上式中可以看出,此式需要儲(chǔ)存5個(gè)中間變量,且比ERK多一個(gè)函數(shù)運(yùn)算,其穩(wěn)定性條件為

(15)

(3)低存儲(chǔ)五級(jí)四階顯式RK方法(low-storage ERK,LSERK)

(16)

式中pi和ki為中間變量;αi,βi,γi為L(zhǎng)SERK的系數(shù)(Hesthaven and Warburton,2008).其穩(wěn)定性條件為

(17)

由(16)式可以看出,LSERK只需要儲(chǔ)存兩個(gè)中間變量即可.逐級(jí)進(jìn)行求解時(shí),只保存此級(jí)的p和k值,因此對(duì)內(nèi)存的需求小.SSPRK和LSERK雖然都需要進(jìn)行5級(jí)迭代,但允許更大的時(shí)間步長(zhǎng),這相當(dāng)于抵消了多余的計(jì)算量.因此綜合來(lái)看,LSERK這種方法比其他兩種RK方法更具有優(yōu)勢(shì).

2.3 單元大小與基函數(shù)階次

在采用DGTD算法進(jìn)行GPR正演模擬時(shí),單元大小和基函數(shù)的階次對(duì)正演模擬的計(jì)算效率和計(jì)算精度有著至關(guān)重要的影響.在直角等腰三角形I中,二維N階基函數(shù)的定義為(Hesthaven and Warburton,2008)：

i+j≤N,

(18)

其中,

(19)

當(dāng)基函數(shù)階次較低或者網(wǎng)格數(shù)較少時(shí),計(jì)算量小,計(jì)算效率高,但計(jì)算精度低.隨著基函數(shù)階次的提高或網(wǎng)格數(shù)的增大,將會(huì)提高計(jì)算精度,但也會(huì)增大計(jì)算量,降低計(jì)算效率.因此,需要對(duì)單元大小與基函數(shù)階次進(jìn)行影響因素分析,達(dá)到計(jì)算效率與計(jì)算精度之間的平衡.

2.4 網(wǎng)格剖分方式

在GPR正演模擬時(shí),需要了解DGTD算法對(duì)不同網(wǎng)格的敏感程度,在網(wǎng)格剖分時(shí)避免因網(wǎng)格畸變所引起的計(jì)算精度降低.Delaunay剖分和前沿推進(jìn)剖分都是GPR正演模擬中常用的網(wǎng)格剖分方式.所有三角形的外接圓均滿足空?qǐng)A性質(zhì)的三角剖分,稱為Delaunay三角剖分,即每個(gè)三角形的外接圓范圍內(nèi)(邊界除外),不包含其他頂點(diǎn)(張嶺和郝天珧,2006).而前沿推進(jìn)法以相鄰三點(diǎn)為基礎(chǔ),檢驗(yàn)夾角的大小.如若為鈍角,則在角的平分線上生成一個(gè)新的節(jié)點(diǎn);如若為銳角,則將兩條前沿的非公共點(diǎn)直接相鄰,生成一個(gè)新的單元.前沿推進(jìn)法生成的網(wǎng)格具有適應(yīng)能力強(qiáng),網(wǎng)格質(zhì)量高等優(yōu)勢(shì)(馬鈞霆等,2015).因此,本文為了探究DGTD算法對(duì)于網(wǎng)格的依賴程度,選取了四種網(wǎng)格與常規(guī)用的非規(guī)則Delaunay三角剖分、前沿推進(jìn)剖分進(jìn)行對(duì)比.四種網(wǎng)格分別為：(a)等腰直角三角形;(b)等邊三角形;(c)頂角為120°的等腰三角形;(d)頂角為30°的等腰三角形.其示意圖如圖1所示.

3 DGTD精度與效率分析

為了分析數(shù)值通量、時(shí)間離散格式、單元大小與基函數(shù)階次、網(wǎng)格剖分方式對(duì)DGTD算法精度和效率的影響,選取了具有解析解的均勻模型來(lái)進(jìn)行正演模擬,以便分析誤差.此模型區(qū)域大小為1 m×1 m,其相對(duì)介電常數(shù)εr=4,電導(dǎo)率σ=6 mS·m-1,相對(duì)磁導(dǎo)率μr=1.激勵(lì)源為主頻900 MHz的雷克子波,在(0.5 m, 0.5 m)處激發(fā)(子圖中的三角形)(圖2).接收點(diǎn)位于(0.6 m, 0.6 m)(子圖中的圓點(diǎn)).誤差公式為

圖1 不同網(wǎng)格形狀示意圖(a) 等腰直角三角形; (b) 等邊三角形; (c) 頂角為120°的等腰三角形; (d) 頂角為30°的等腰三角形.

(20)

其中,utest為采用不同參數(shù)的DGTD算法正演模擬的結(jié)果;uref為二維均勻模型的解析解：

(21)

圖2 二維均勻模型及其解析解灰色子圖為模型示意圖,三角形為激勵(lì)源,圓形為接收點(diǎn).

3.1 數(shù)值通量的影響

利用圖2所示的二維均勻模型,分析迎風(fēng)通量、中心通量和補(bǔ)償通量對(duì)GPR正演模擬的影響.利用三角形對(duì)整個(gè)模擬區(qū)域進(jìn)行剖分,基函數(shù)階次N=3,時(shí)間采樣間隔dt=1×10-11s,總模擬時(shí)間為10 ns.表2給出了不同數(shù)值通量的誤差.為了更加直觀地感受不同τ值的數(shù)值流量與解析解之間的誤差,選擇了幾個(gè)典型流量做了A-Scan對(duì)比圖,如圖3所示.

從表2中可以看出,中心通量的誤差是最大的,而在τ大于等于1/8后,誤差基本在(3～6)×10-4之間,相差無(wú)幾.這是因?yàn)樵谑?10)中,υ所乘的那一項(xiàng)稱為耗散項(xiàng).中心通量的υ=0,是沒(méi)有耗散的,但是會(huì)出現(xiàn)偽解.在這個(gè)案例中,它的相對(duì)誤差達(dá)到了31.7236,誤差較大.補(bǔ)償通量和迎風(fēng)通量雖是耗散的,但是可以解決偽解的問(wèn)題.因此,耗散項(xiàng)對(duì)于GPR正演模擬是非常有必要的.并且從表2中也可以看到,τ的值不能過(guò)小,否則相對(duì)誤差會(huì)很大,給之后的GPR數(shù)據(jù)解譯帶來(lái)困難.同樣的結(jié)論也可以由圖3得到.當(dāng)τ=1/50和τ=1/20時(shí),波形與解析解存在很大的差異.而隨著τ值的增大,波形與解析解基本一致,誤差減小.在該模型中,τ=1/4時(shí),誤差最小.為了尋找到最優(yōu)的τ值,對(duì)物性參數(shù)不同的模型采取了不同的基函數(shù)階次和網(wǎng)格進(jìn)行模擬,結(jié)果表明,τ的最優(yōu)區(qū)間為[1/8,1].最優(yōu)的τ值是根據(jù)模型變化的,不同的單元大小、基函數(shù)階次也會(huì)改變其取值.為了避免出現(xiàn)偽解,又提高計(jì)算精度,推薦采用τ=1/2的補(bǔ)償通量來(lái)進(jìn)行高階DGTD正演模擬.若想進(jìn)一步提高精度,則需要對(duì)具體的模型案例做額外的分析實(shí)驗(yàn).

表2 不同數(shù)值通量的相對(duì)誤差值Table 2 Relative errors of different numerical fluxes

圖3 不同τ值的數(shù)值通量與解析解之間的A-Scan對(duì)比圖

3.2 時(shí)間離散格式的影響

為了對(duì)比ERK、SSPRK和LSERK的優(yōu)勢(shì),同樣利用圖2所示的二維均勻模型來(lái)進(jìn)行正演模擬.因?yàn)镋RK允許的時(shí)間步長(zhǎng)較小,而SSPRK和LSERK允許的時(shí)間步長(zhǎng)較大,因此先采用一個(gè)臨界的時(shí)間步長(zhǎng)值對(duì)比其穩(wěn)定性條件.設(shè)置dt=3.9×10-11s,時(shí)間模擬總時(shí)長(zhǎng)為10 ns,采用補(bǔ)償通量τ=1/2來(lái)進(jìn)行空間離散.三種離散方式下不同時(shí)刻的波場(chǎng)快照如圖4所示.由圖可見(jiàn),在設(shè)置的時(shí)間步長(zhǎng)條件下,2 ns時(shí),三者的波場(chǎng)快照都是穩(wěn)定的狀態(tài),并沒(méi)有什么區(qū)別;而隨著時(shí)間的推移,3 ns時(shí),ERK已經(jīng)開始出現(xiàn)色散的跡象,而基于SSPRK和LSERK時(shí)間離散格式的波場(chǎng)快照很光滑,不存在雜波;4 ns時(shí),ERK的波場(chǎng)快照已經(jīng)被污染,而SSPRK和LSERK還是很穩(wěn)定的,并沒(méi)有產(chǎn)生色散.因此可以看出,ERK的穩(wěn)定性條件比SSPRK和LSERK的穩(wěn)定性條件低,SSPRK和LSERK時(shí)間離散方案可以允許更大的時(shí)間步長(zhǎng).

圖4 三種時(shí)間離散方案下不同時(shí)刻的波場(chǎng)快照(a), (b), (c) 分別為ERK、SSPRK和LSERK在2 ns時(shí)刻的波場(chǎng)快照;(d), (e), (f)為3 ns;(g), (h), (i)為4 ns.

此外,為了更加了解這三個(gè)時(shí)間離散方案的優(yōu)勢(shì)和劣勢(shì),在滿足各自穩(wěn)定性條件下,采用極限時(shí)間步長(zhǎng)進(jìn)行離散,對(duì)比其占用內(nèi)存空間、正演模擬所需時(shí)間和精度,對(duì)比結(jié)果如表3所示.ERK為四級(jí)四階,而SSPRK和LSERK為五級(jí)四階,但當(dāng)三種時(shí)間離散格式分別采取了不同的時(shí)間步長(zhǎng)時(shí),其正演模擬時(shí)間幾乎相同.這是因?yàn)镾SPRK和LSERK雖然每一個(gè)時(shí)間步會(huì)多一個(gè)函數(shù)計(jì)算,但是它們?cè)试S更大的時(shí)間步長(zhǎng),這相當(dāng)于把額外的工作量抵消;在精度上,因?yàn)槿N方法均為4階精度,所以誤差保持在同一水平,約等于4.31×10-4;在儲(chǔ)存量上則有很大的差別,每一次啟動(dòng)時(shí)間離散方案時(shí),ERK需要額外儲(chǔ)存4個(gè)矩陣,SSPRK需要儲(chǔ)存5個(gè)矩陣,而LSERK只需要額外儲(chǔ)存2個(gè)矩陣.對(duì)于大型的復(fù)雜模型和三維正演模擬,LSERK有著巨大的內(nèi)存優(yōu)勢(shì),因此本文在接下來(lái)的影響因素分析中,將采取LSERK方案進(jìn)行時(shí)間上的離散.

表3 不同離散時(shí)間格式的計(jì)算條件及特征Table 3 Characteristics and calculation conditions for different temporal integrations

3.3 單元大小與基函數(shù)階次的影響

通過(guò)以上的實(shí)驗(yàn)可知,補(bǔ)償通量τ=1/2和LSERK時(shí)間離散方案在GPR正演模擬中有較大的優(yōu)勢(shì).除此之外,網(wǎng)格單元的大小與基函數(shù)的階次對(duì)收斂速度也有很大的影響.同樣利用圖2所示的二維均勻模型進(jìn)行正演模擬,采用不同的網(wǎng)格尺寸對(duì)模型進(jìn)行空間離散,即網(wǎng)格數(shù)K分別為128、200、512、800、1250和2048;此外,在不同單元大小的情況下,分別采用基函數(shù)階次N=1, 2, 3, 4, 5, 6來(lái)進(jìn)行正演模擬,其相對(duì)誤差值如表4所示.從中可以看出,網(wǎng)格數(shù)的增大和基函數(shù)階次的提高均可以減小誤差.而誤差到達(dá)1.8×10-4時(shí),不再向下減小,這是因?yàn)檎菽M的誤差由空間離散和時(shí)間離散共同決定,此時(shí)已達(dá)到其時(shí)間離散的最高精度.因此當(dāng)再減小網(wǎng)格尺寸和提高基函數(shù)階次時(shí),誤差將不會(huì)減小.為了對(duì)其收斂性進(jìn)一步分析,繪制了如圖5所示的誤差趨勢(shì)圖.圖5a為單元數(shù)K對(duì)誤差的影響,從中可以看出,基函數(shù)階次相同的情況下,網(wǎng)格數(shù)越大,網(wǎng)格尺寸越小,其誤差越小;此外,分析不同階次下誤差曲線的斜率可以看出,階次越高,其誤差收斂性越好.圖5b為階次N對(duì)誤差的影響.同一網(wǎng)格下,階次越高,誤差越小;此外,網(wǎng)格數(shù)越多,其誤差收斂性越高.分析兩圖可知,提高基函數(shù)階次N或增大網(wǎng)格數(shù)K,均可以提高其誤差的收斂性.

表4 不同單元大小與基函數(shù)階次的相對(duì)誤差值Table 4 Relative errors of different grid sizes and orders of basis function

圖5 誤差曲線趨勢(shì)圖(a) 誤差隨單元K的變化; (b) 誤差隨階次N的變化.

從圖5可以發(fā)現(xiàn),階次較小時(shí),需要的網(wǎng)格數(shù)很多.而當(dāng)階次變大時(shí),因時(shí)間離散的限制原因,增加網(wǎng)格數(shù)并不能提高計(jì)算精度,反而會(huì)增加計(jì)算量,降低計(jì)算效率.因此,需要在單元大小與基函數(shù)階次之間尋找一個(gè)合適的關(guān)系,平衡計(jì)算效率與計(jì)算精度.電磁波波長(zhǎng)λ計(jì)算公式為

(22)

其中,c為電磁波在介質(zhì)中傳播的速度;f為激勵(lì)源的中心頻率;c0為電磁波在真空中的傳播速度,約為3.0×108m·s-1;εr為相對(duì)介電常數(shù);μr為相對(duì)磁導(dǎo)率.由(22)式計(jì)算出,圖2所示的二維均勻模型背景介質(zhì)的電磁波波長(zhǎng)為16.6 cm.圖5a中紅色虛線所示的誤差為5×10-4.根據(jù)斜率計(jì)算出3,4,5,6階在誤差為0.5‰時(shí)所需要的網(wǎng)格數(shù)量,進(jìn)而計(jì)算出單元大小.計(jì)算結(jié)果如表5所示.

從表中,可以得到d/N≈λ/15這個(gè)關(guān)系式.為測(cè)試此公式的普遍性,采用不同的基函數(shù)階次、不同大小的網(wǎng)格進(jìn)行正演模擬,計(jì)算其單元大小與誤差,得到相同的結(jié)果.

表5 誤差為0.5‰時(shí),單元大小與基函數(shù)階次之間的關(guān)系Table 5 The relationship between the grid size and the basis function order when the error is 0.5‰

3.4 網(wǎng)格剖分方式的影響

在GPR正演模擬中,為了避免因網(wǎng)格畸形引起的誤差,需要了解不同的網(wǎng)格剖分方式對(duì)于計(jì)算精度的影響.因此,采用了如圖1所示的4種網(wǎng)格來(lái)進(jìn)行剖分,分別為等腰直角三角形、等邊三角形、頂角為120°的等腰三角形及頂角為30°的等腰三角形,將它們與Delaunay剖分、前沿推進(jìn)剖分方法進(jìn)行對(duì)比.為了更加清楚地探究網(wǎng)格剖分方式的影響,只改變剖分三角形的形狀,其他模擬條件保持一致.設(shè)置模擬區(qū)域的網(wǎng)格數(shù)K大致相等,約為500個(gè)網(wǎng)格.采用不同的基函數(shù)階次N(1,2,3,4,5,6)進(jìn)行正演模擬,所得到的誤差值如圖6所示.由圖可見(jiàn),隨著基函數(shù)階次的提升,四個(gè)網(wǎng)格的誤差均在減小;而基函數(shù)階次相同時(shí),等邊三角形的誤差是最小的,其次是Delaunay剖分三角形和前沿推進(jìn)剖分,而頂角為120°的等腰三角形誤差最大.當(dāng)基函數(shù)的階次變大時(shí),因網(wǎng)格剖分方式所引起的誤差也在減小.總體上,它們都處于一個(gè)誤差級(jí)別,相差不大.由此可得,高階DGTD算法對(duì)于網(wǎng)格具有較好的適應(yīng)性,網(wǎng)格剖分方式上的不足可以由更高階次的基函數(shù)來(lái)補(bǔ)充.

圖6 不同網(wǎng)格剖分方式的誤差圖

4 數(shù)值案例

中國(guó)首次火星探測(cè)任務(wù)“天問(wèn)一號(hào)”配備了探地雷達(dá)系統(tǒng),于2020年7月23日成功發(fā)射,于2021年5月15日成功著陸于火星烏托邦平原南部(Dong et al.,2022).選取最優(yōu)的DGTD參數(shù)組合,應(yīng)用到了火星烏托邦平原正演模擬中,表明所分析的DGTD算法對(duì)于復(fù)雜模型的適應(yīng)性和準(zhǔn)確性.根據(jù)烏托邦平原可能的地下結(jié)構(gòu)(Li et al.,2022),建立圖7a所示的模型進(jìn)行數(shù)值模擬.整個(gè)模擬區(qū)域?yàn)?20 m×60 m,相對(duì)磁導(dǎo)率為1,各層的介質(zhì)和相對(duì)介電常數(shù)、電導(dǎo)率在表6中給出.編號(hào)②的地層中,有一些灰色的小巖塊,為含冰玄武巖,用來(lái)模擬火山噴發(fā)時(shí)上涌的巖塊,后又被沉積物所覆蓋,并且離地面越近,巖塊越小.采用DGTD算法進(jìn)行正演模擬時(shí),網(wǎng)格總數(shù)為16295,基函數(shù)階次N=3,時(shí)間采樣間隔dt=3×10-11s,總模擬時(shí)間為800 ns.采用主頻為20 MHz的雷克子波作為激勵(lì)源,激發(fā)點(diǎn)和接收點(diǎn)都設(shè)置在了地表處,如圖7a中的紅色三角形所示,總共接收195道數(shù)據(jù).

圖7b為DGTD算法模擬的地質(zhì)雷達(dá)剖面.為了更好地解譯地下介質(zhì),對(duì)該剖面做高度校正,讓其更加貼合實(shí)際地形.由圖7b可見(jiàn),大約在170～200 ns、0～96 m處有一條反射波,那是①和②的地層分界面;400ns的位置也存在一條明顯的反射波,那是②的下界面,經(jīng)過(guò)高度校正后的反射波波形準(zhǔn)確地反映了含冰風(fēng)成沉積物的起伏.在0～60 m的水平位置,該反射波的幅值較小,而60～120 m處,地層分界面較明顯.這是由于③磁鐵礦的相對(duì)介電常數(shù)與②差異較小,因此反射波較弱.圖7b中紅色橢圓形A所圈出的位置,為上涌到②號(hào)地層(含冰風(fēng)成沉積物)處的小含冰玄武巖塊的雙曲線反射波;虛線矩形B所表示的位置為③號(hào)地層的下界面;實(shí)線矩形C處為④與⑤的地層分界面.這些反射波都與圖7a的模型相對(duì)應(yīng).這說(shuō)明了最優(yōu)參數(shù)組合的DGTD算法在探地雷達(dá)正演模擬中的準(zhǔn)確性,可指導(dǎo)火星實(shí)測(cè)數(shù)據(jù)的解譯.

圖7 復(fù)雜模型示意圖及其剖面圖(a) 火星烏托邦平原二維模型示意圖; (b) 剖面圖.

5 結(jié)論

(1) 作者將DGTD算法引入到GPR正演模擬中,推導(dǎo)了DGTD算法的弱解形式,介紹了三種數(shù)值通量和時(shí)間離散方案.從數(shù)值通量、時(shí)間離散格式、單元大小與基函數(shù)階次、網(wǎng)格剖分方式4個(gè)方面,分析了DGTD算法的計(jì)算精度和效率.

(2) 數(shù)值模擬實(shí)驗(yàn)表明,Lax-Friedrichs數(shù)值通量中τ=1/2時(shí)的補(bǔ)償通量,既可以消除偽解,又可以提高計(jì)算精度;低存儲(chǔ)顯式Runge-Kutta方案在GPR模擬中更具有優(yōu)勢(shì).而提高基函數(shù)的階次或增大網(wǎng)格數(shù),均可以提高算法精度;誤差為0.5‰時(shí),網(wǎng)格與階次的一個(gè)適用關(guān)系為d/N≈λ/15.不同的網(wǎng)格剖分方式表明,網(wǎng)格的形狀對(duì)于高階DGTD算法的影響較小,也證明DGTD算法對(duì)于網(wǎng)格的適應(yīng)性強(qiáng).

(3) 本文中作者采用DGTD算法僅開展了二維TM模式下的GPR正演.下一步的研究方向?yàn)閷GTD算法應(yīng)用到大尺度三維GPR正演模擬中,并實(shí)現(xiàn)hp型自適應(yīng),拓寬DGTD算法的適用性,為GPR反演奠定理論基礎(chǔ).