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

1 中南大學地球科學與信息物理學院, 長沙 410083

2 東華理工大學地球科學學院, 南昌 330032

0 引言

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

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

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

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

1 DGTD算法求解GPR波動方程

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

(1)

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

為了簡化方便推導,將(1)式化簡為(Lu et al.,2004)

(2)

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

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

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

(4)

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

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

(5)

(6)

(7)

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

(8)

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

(9)

2 DGTD 算法影響因素

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

2.1 數(shù)值通量

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

(10)

(11)

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

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

2.2 時間離散

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

(1)標準四級四階顯式RK方法(explicit RK,ERK)

(12)

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

(13)

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

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

(14)

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

(15)

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

(16)

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

(17)

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

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

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

i+j≤N,

(18)

其中,

(19)

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

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

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

3 DGTD精度與效率分析

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

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

(20)

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

(21)

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

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

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

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

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

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

3.2 時間離散格式的影響

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

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

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

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

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

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

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

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

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

(22)

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

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

表5 誤差為0.5‰時,單元大小與基函數(shù)階次之間的關系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)格剖分方式對于計算精度的影響.因此,采用了如圖1所示的4種網(wǎng)格來進行剖分,分別為等腰直角三角形、等邊三角形、頂角為120°的等腰三角形及頂角為30°的等腰三角形,將它們與Delaunay剖分、前沿推進剖分方法進行對比.為了更加清楚地探究網(wǎng)格剖分方式的影響,只改變剖分三角形的形狀,其他模擬條件保持一致.設置模擬區(qū)域的網(wǎng)格數(shù)K大致相等,約為500個網(wǎng)格.采用不同的基函數(shù)階次N(1,2,3,4,5,6)進行正演模擬,所得到的誤差值如圖6所示.由圖可見,隨著基函數(shù)階次的提升,四個網(wǎng)格的誤差均在減小;而基函數(shù)階次相同時,等邊三角形的誤差是最小的,其次是Delaunay剖分三角形和前沿推進剖分,而頂角為120°的等腰三角形誤差最大.當基函數(shù)的階次變大時,因網(wǎng)格剖分方式所引起的誤差也在減小.總體上,它們都處于一個誤差級別,相差不大.由此可得,高階DGTD算法對于網(wǎng)格具有較好的適應性,網(wǎng)格剖分方式上的不足可以由更高階次的基函數(shù)來補充.

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

4 數(shù)值案例

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

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

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

5 結(jié)論

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

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

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