王樹奇 黃茸茸 吳楠

摘 要：為了獲取準(zhǔn)確的井下采空區(qū)的地理信息，利用時(shí)域不連續(xù)伽遼金法（DGTD）和高階疊層矢量基函數(shù)的特點(diǎn)，建立了井下采空區(qū)三維空間模型，研究了高階疊層矢量基函數(shù)階數(shù)對(duì)井下采空區(qū)電磁正演精度的影響，得到井下采空區(qū)電磁場(chǎng)強(qiáng)度的空間分布規(guī)律。井下采空區(qū)電磁正演模擬結(jié)果表明，采用2.5階疊層矢量基函數(shù)的DGTD方法計(jì)算精度比傳統(tǒng)低階棱邊基函數(shù)提高了74%，證明了該方法對(duì)復(fù)雜目標(biāo)電磁特性研究的實(shí)用性，并能準(zhǔn)確可靠地實(shí)現(xiàn)井下采空區(qū)地質(zhì)結(jié)構(gòu)的電磁探測(cè)。關(guān)鍵詞：井下采空區(qū);時(shí)域不連續(xù)伽遼金法;高階疊層矢量基函數(shù);電磁正演中圖分類號(hào)：TM 153

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1672-9315（2021）02-0348-07

DOI：10.13800/j.cnki.xakjdxxb.2021.0220開(kāi)放科學(xué)（資源服務(wù)）標(biāo)識(shí)碼（OSID）：

Electromagnetic forward modeling based on discontinuous

Galerkin time domain method in underground goaf

WANG Shuqi1，HUANG Rongrong1，WU Nan2

（1.College of Communication and InformationEngineering，Xian University of Science and Technology，Xian 710054，China;

2.School of Information Engineering，Shaanxi Xueqian Normal University，Xian 710100，China）Abstract：In order to obtain accurate geographic information of underground goaf，a threedimensional spatial model of underground goaf is established by using the characteristics of discontinuous Galerkin time domain （DGTD）and highorder hierarchical vector basis functions.The influence of the order of highorder hierarchical vector basis functions on the accuracy of electromagnetic forward modeling of underground goaf is examined and the spatial distribution and response characteristics of electromagnetic field intensity in underground goaf are analyzed.The results of electromagnetic forward modeling of goaf show that the calculation accuracy of DGTD method using hierarchical vector basis functions of order 2.5 is 74% higher than that of the traditional loworder edgebased basis function，indicating the practicability of this method in the study of electromagnetic characteristics of complex targets.And it is able to accurately and reliably realize the electromagnetic detection of geological structure in goaf.Key words：underground goaf;time domain discontinuous Galerkin method;high order hierarchical vector basis functions;electromagnetic forward modeling0 引 言

煤炭是中國(guó)能源支柱產(chǎn)業(yè)，然而大范圍的煤炭開(kāi)采在礦山中留下了大量的采空區(qū)，導(dǎo)致地下巖體原有的力學(xué)平衡被打破，扭曲的巖體隨時(shí)可能發(fā)生位移、巖爆等事故;更嚴(yán)重的是采空區(qū)會(huì)被瓦斯、地下水等充填，給礦區(qū)的開(kāi)采作業(yè)埋下重大的安全隱患，采用電磁數(shù)值計(jì)算對(duì)礦井復(fù)雜地質(zhì)結(jié)構(gòu)進(jìn)行探測(cè)分析具有重要的意義。

國(guó)內(nèi)外學(xué)者針對(duì)時(shí)域電磁正演的方法做了大量研究，但對(duì)于井下采空區(qū)的時(shí)域電磁正演的研究多采用時(shí)域有限差分（finite difference time domain，F(xiàn)DTD）方法和時(shí)域矢量有限元（finite element time domain，F(xiàn)ETD）[1-3]。于景邨和常江浩采用三維有限差分法對(duì)老空水的礦井瞬變電磁響應(yīng)進(jìn)行了研究[4]。楊道學(xué)提出了基于卷積完全匹配層的交替方向隱式的FDTD方法，并將其應(yīng)用于探地雷達(dá)正演[5]。FDTD方法具有理論難度低，計(jì)算效率高等優(yōu)點(diǎn)，但是該方法中使用的六面體剖分單元在復(fù)雜幾何結(jié)構(gòu)中具有局限性，且不支持高階基函數(shù)等問(wèn)題將直接影響計(jì)算精度。張永超和李宏杰等研究了時(shí)域矢量有限元三維正演[6]，拓展了礦井瞬變電磁正演對(duì)復(fù)雜地質(zhì)模型的適用性，棱邊矢量基函數(shù)的使用很好地解決了節(jié)點(diǎn)方法中的“偽解”，但只是用了低階棱邊基函數(shù)且隱式離散格式，計(jì)算效率低。時(shí)域不連續(xù)伽遼金方法（discontinuous galerkin time domain，DGTD）[7-12] 繼承了FETD 采用非結(jié)構(gòu)網(wǎng)格對(duì)復(fù)雜幾何外形擬合好，便于使用高階基函數(shù)的優(yōu)點(diǎn);該算法對(duì)Maxwell 方程采用伽遼金加權(quán)法得到弱解形式，放松了單元之間的邊界條件，單元之間通過(guò)數(shù)值通量聯(lián)系交換數(shù)據(jù);高階疊層矢量基函數(shù)的引入是一個(gè)突破，它滿足單元邊界上的 場(chǎng)切向連續(xù)性，消除了節(jié)點(diǎn)法中存在的“偽解”并且可顯著提高該方法的計(jì)算精度;在時(shí)間離散方面既支持顯式時(shí)間格式也支持隱式時(shí)間格式;在空間離散方面由于單元間的高度獨(dú)立性，對(duì)于每個(gè)單元可以選用不同階的基函數(shù)進(jìn)行離散，需要高精度時(shí)采用高階基函數(shù)，對(duì)精度要求低時(shí)采用低階基函數(shù)?？梢詫⒋笮投喑叨葐?wèn)題被拆分成為一個(gè)個(gè)相對(duì)獨(dú)立的小問(wèn)題，使得分析問(wèn)題更為簡(jiǎn)單方便，以上原因更加說(shuō)明了研究DGTD算法的價(jià)值與意義。文中首次將基于高階疊層矢量基函數(shù)的DGTD方法應(yīng)用于井下采空區(qū)三維正演研究，建立了井下采空區(qū)三維空間模型，分析了疊層矢量基函數(shù)階數(shù)對(duì)正演結(jié)果的影響，試驗(yàn)驗(yàn)證了井下采空區(qū)模型正演的計(jì)算精度。

1 DGTD方法

三維非均勻各向同性介質(zhì)Maxwell旋度方程為

μ

H

t

=-×E-σm

H

ε

E

t

=×H-σe

E

（1）

式中 ε為介電系數(shù);μ為磁導(dǎo)系數(shù);σe為電導(dǎo)率;σm為導(dǎo)磁率;

E

和H為電場(chǎng)強(qiáng)度與磁場(chǎng)強(qiáng)度。對(duì)公式（1）采用伽遼金加權(quán)法，在四面體單元內(nèi)進(jìn)行體積分，可以得到以下弱解形式

∫T m

μ

H

t

+×E+σm

H

·q′dv=0

∫T m

ε

E

t

-×H+σe

E

·q′dv=0

（2）

q′為加權(quán)函數(shù)，使用以下恒等式對(duì)公式（2）中的

×E和

×H項(xiàng)進(jìn)行代數(shù)變換

∫T m

（×H）

·mq′dv

=∫T m·（H×

mq′）dv

+

∫T m（×

mq′）·Hdv

=T m（m×H）

·

mq′ds

+

∫T m（×

mq′）·Hdv

∫T m

（×E）

·mq′dv

=∫T m·（E×

mq′）dv

+

∫T m（×

mq′）·Edv

=T m（m×E）

·

mq′ds

+

∫T m（×

mq′）·Edv

（3）在相鄰區(qū)域單元交界面上采用數(shù)值通量[13-14]的形式來(lái)保證單元之間切向場(chǎng)的連續(xù)性，數(shù)值通量格式如下

m×Hm*=

m×Hm

+kmh[m×（Hm+-Hm）]

-

vme[m×（m×

（Em+-Em））]

m×Em*=

m×Em

+kme[m×（Em+-Em）]

+

vmh[m×（m×

（Hm+-Hm））]

（4）

式中

是四面體的外法向單位矢量;

Hm*，

Em*稱為數(shù)值通量，

定義在2個(gè)相鄰單元的交界面上，由2個(gè)相接面的電磁場(chǎng)合成，數(shù)值通量在放寬單元之間連續(xù)性邊界條件的基礎(chǔ)上，實(shí)現(xiàn)了相鄰單元之間場(chǎng)量信息的交互;

Em，

Hm

為當(dāng)前四面體單元的電場(chǎng)磁矢量;

Em+，

Hm+

為相鄰四面體單元的電磁場(chǎng)矢量。

將數(shù)值通量代入公式（3）中可得

∫T m（×Hm*）·

mq′dv=

∫T m

（×Hm）·

mq′dv

-

T m

kmh[

m×（Hm-Hm+）]·mq′ds+

T m

vme[

m×（

m×（Em-Em+）]·mq′ds

∫T m（×Em*）·

mq′dv=

∫T m

（×Em）·

mq′dv

-

T m

kme[

m×（Em-Em+）]·mq′ds-

T m

vmh[

m×（

m×（Hm-Hm+）]·mq′ds

（5）公式（5）的空間部分采用基函數(shù)在每個(gè)單元中展開(kāi)

H≈H=∑Qmq=1hmq（t）mq（r）

E≈E=∑Qmq=1emq（t）mq（r）

（6）整理簡(jiǎn)化后得到基于四面體單元的DGTD顯式半離散方程組

μMdtHm

+（σmM-Fvh）Hm

+F+vhHm+

=-（S-Fke）Em-F+keEm+

εMdtEm

+（σeM-Fve）Em

+F+veEm+

=（S-Fkh）Hm+F+khHm+

（7）在得到空間離散格式之后將進(jìn)行時(shí)間離散，最后根據(jù)時(shí)間離散格式完成時(shí)間迭代?？刹扇〔煌绞綄?duì)得到的空間離散方程進(jìn)行時(shí)間偏導(dǎo)數(shù)的離散，若采用隱式時(shí)間離散方案則將降低計(jì)算效率，文中采用已經(jīng)廣泛應(yīng)用于時(shí)域算法中的蛙躍算法進(jìn)行時(shí)間離散[15]，蛙躍算法是一種顯式時(shí)域離散方法，電磁場(chǎng)相隔半個(gè)時(shí)間步交替進(jìn)行迭代計(jì)算。規(guī)定電場(chǎng)E在整數(shù)時(shí)間步tn采樣，磁場(chǎng)H在半整數(shù)步tn+1/2進(jìn)行采樣。在公式（7）中對(duì)時(shí)間的導(dǎo)數(shù)采取中心差分的格式處理，由二階中心差分代替一階時(shí)間導(dǎo)數(shù)[16]

（dtHm）n

=Hmn+1/2-

Hmn-1/2

Δt

+O（Δt2）

（dtEm）n+1/2

=Emn+1-

Emn

Δt

+O（Δt2）

（8）

由于在迭代過(guò)程中會(huì)出現(xiàn)電場(chǎng)的半整數(shù)步、磁場(chǎng)的整數(shù)步，對(duì)其取平均近似[16]有以下形式

Hmn

=Hmn+1/2+

Hmn-1/2

2

+O（Δt2）

Emn+1/2

=Emn+1+

Emn

2

+O（Δt2）

（9）

通過(guò)初始的激勵(lì)源與邊界條件，就可以利用DGTD逐步迭代求得計(jì)算域中電磁場(chǎng)分布隨時(shí)間的變化

Hmn+1/2

=αmHmn-1/2

+βmM-1[-（S-Fke）

Emn

-

F+ke

Em+n

+FvhHmn-1/2

-F+vhHm+n-1/2]Emn+1

=αeEmn

+βeM-1[（S-Fkh）

Hmn+1/2+

F+kh

Hm+n+1/2

+FveEmn

-F+veEm+n

]

（10）

其中系數(shù)為

αm=1-Δtσm2μ

1+Δtσm2μ

，βm=Δt

μ

1+Δtσm2μ

αe=1-Δtσe2ε

1+Δtσe2ε

，βe=Δt

ε

1+Δtσe2ε

2 高階疊層矢量基函數(shù)為了克服低階矢量基函數(shù)精度較低的缺點(diǎn)，NOTAROS等提出了高階插值型矢量基函數(shù)[17-19]。高階插值型矢量基函數(shù)具有線性獨(dú)立性好、物理解釋明確、編程實(shí)現(xiàn)方便等優(yōu)點(diǎn)。但在一個(gè)四面體內(nèi)插值型矢量基函數(shù)具有更多的未知量，為了滿足交界面處切向場(chǎng)連續(xù)，插值矢量基不允許不同階次基函數(shù)的混合，這導(dǎo)致了在原本低階單元即可足夠精確描述的區(qū)域引入了大量冗余自由度，從而增加了存儲(chǔ)量和計(jì)算時(shí)間，降低了算法分析的效率。高階疊層型矢量基函數(shù)[20-21]的提出正是為了解決高階插值型矢量基函數(shù)的這些缺點(diǎn)。首先，高階疊層型基函數(shù)具有疊層嵌套的特點(diǎn)，高階疊層構(gòu)矢量基函數(shù)的構(gòu)建是基于低階基函數(shù)逐階遞增的，即高階的疊層基函數(shù)里面包含有階數(shù)較低的基函數(shù)，這在保證高精度的情況下降低了基函數(shù)構(gòu)造難度以及存儲(chǔ)量。其次，借助高階疊層矢量基函數(shù)，時(shí)域不連續(xù)伽遼金算法可以實(shí)現(xiàn)計(jì)算域內(nèi)高階和低階單元的混合計(jì)算，數(shù)值結(jié)果表明高階疊層矢量基函數(shù)的引用使得對(duì)目標(biāo)離散時(shí)采用較大尺寸單元得到的結(jié)果與采用較小尺寸單元剖分的低階基函數(shù)方案有相當(dāng)?shù)木?。疊層矢量基函數(shù)是一種基于棱邊、面和體積的基函數(shù)，它將自由度賦予棱邊、面和體積而不是單元節(jié)點(diǎn)。它隱含了散度為零的邊界條件，消除了偽解，非常適合用來(lái)表示矢量場(chǎng)。綜上可知，高階疊層矢量基函數(shù)在處理大規(guī)模電磁問(wèn)題和多尺度電磁問(wèn)題上都具有非常大的優(yōu)勢(shì)，具有很高的實(shí)用性與靈活性。表1給出了文中所采用的基于四面體單元的0.5階到2.5階疊層矢量基函數(shù)形式。當(dāng)然，疊層型矢量基函數(shù)還有更高的階的形式，但隨著基函數(shù)階數(shù)的增加，自由度會(huì)以數(shù)倍增加，這將增加計(jì)算難度和時(shí)間成本。并且當(dāng)基函數(shù)階數(shù)為2.5階時(shí)，該算法就可以獲得較好的計(jì)算精度。因此文中就0.5階至2.5階基函數(shù)展開(kāi)研究。

3 矩形諧振腔模擬計(jì)算為了分析不同高階疊層矢量基函數(shù)的DGTD方法的計(jì)算精度，建立了尺寸為2.0 m×0.3 m×0.2 m的無(wú)損耗矩形PEC邊界諧振腔模型，如圖2、圖3所示。

圖2中腔體被劃分為2 261個(gè)均勻網(wǎng)格尺寸的四面體，圖3所示計(jì)算域被劃分為2 292個(gè)網(wǎng)格大

小不同的四面體單元。在點(diǎn)（0，0，0）m處激發(fā)

一個(gè)中心頻率為310 MHz，極化方向?yàn)椋?，1，1）的布萊克曼-哈里斯（blackman harris window，BHW）電偶極子源。

分別選取0.5階、1.5階、2.5階以及混合階疊層矢量基函數(shù)時(shí)，DGTD方法在觀測(cè)點(diǎn)（0.7，0，0）m處得到的時(shí)域波形分別如圖4（a）～圖4（d）所示。

從圖4可以看出，當(dāng)疊層基函數(shù)為0.5階時(shí)，其結(jié)果的計(jì)算精度較差。隨著高階疊層矢量基函數(shù)的階數(shù)增加，數(shù)值計(jì)算結(jié)果的精度有了顯著提升，經(jīng)分析見(jiàn)表2所示結(jié)果。當(dāng)基函數(shù)階數(shù)分別為0.5階、1.5階、2.5階以及混合階疊層矢量基函數(shù)時(shí)，計(jì)算結(jié)果的相對(duì)誤差為43.94%、2244%、229%和2.42%。從0.5階到2.5階，DGTD方法的計(jì)算相對(duì)誤差降低了4165%，計(jì)算精度得到了明顯的改善;混合階基函數(shù)的計(jì)算中基函數(shù)的階數(shù)根據(jù)四面體單元的尺寸進(jìn)行選擇，其結(jié)果的相對(duì)誤差與2.5階基函數(shù)結(jié)果相近，該結(jié)果表明基于

高階疊層矢量基函數(shù)的DGTD方法在保證計(jì)算精度的同時(shí)可實(shí)現(xiàn)計(jì)算域內(nèi)高階和低階單元的混合計(jì)算，對(duì)復(fù)雜目標(biāo)的電磁特性研究具有很高的實(shí)用性。

4 井下采空區(qū)電磁正演建立煤層區(qū)域?yàn)?0 m×10 m×10 m，井下采空區(qū)尺寸為4 m×4 m×2 m的三維正演地質(zhì)模型，如圖5，圖6所示。背景設(shè)置為均勻介質(zhì)，取相對(duì)介電常數(shù)為4.09褐煤進(jìn)行模擬驗(yàn)證。以空間模型的中心位置作為原點(diǎn)，掘進(jìn)工作面位于z方向上。在（0 m，0 m，4 m）處設(shè)置中心頻率為31 MHz的BHW電偶極子激勵(lì)源。接收點(diǎn)放置于掘進(jìn)工作面上，其位置坐標(biāo)為（0 m，1 m，5 m）。

為了研究所選取的基函數(shù)對(duì)井下電磁響應(yīng)特性的影響，設(shè)計(jì)了地質(zhì)模型，分別將0.5階、1.5階、2.5階疊層矢量基函數(shù)應(yīng)用于電磁正演中，結(jié)果如圖7所示。對(duì)比不同階基函數(shù)對(duì)采空區(qū)的電磁響應(yīng)特征的影響，當(dāng)基函數(shù)為0.5階、1.5階、25階3種不同階數(shù)時(shí)均可分辨反射波信號(hào)。從圖7所示的基于DGTD方法井下采空區(qū)電場(chǎng)響應(yīng)結(jié)果可看出，隨著基函數(shù)階數(shù)的增加，直達(dá)波和反射波的幅值精確度有了顯著提高，且采空區(qū)反射波曲線拖延現(xiàn)象得到了顯著改善。對(duì)比DGTD方法中不同階疊層基函數(shù)正演波形，誤差分析結(jié)果見(jiàn)表3。計(jì)算中將接收點(diǎn)的場(chǎng)值記錄并與仿真軟件結(jié)果進(jìn)行比較，其相對(duì)誤差定義為

Error=20log（|E-Eref|/Eref max）（dB）

式中 E為時(shí)域波形數(shù)值解;Eref為仿真軟件計(jì)算結(jié)果。

經(jīng)分析得：井下采空區(qū)正演的相對(duì)誤差從0.5階的-19.12 dB縮小至2.5階的-73.19 dB，計(jì)算精度提高了74%。由以上分析可得高階疊層矢量基函數(shù)的階數(shù)將直接影響DGTD方法的正演計(jì)算精度，采空區(qū)正演時(shí)域結(jié)果隨基函數(shù)階數(shù)的增加，幅值的精確度大幅提升。

5 結(jié) 論

1）利用基于高階疊層矢量基函數(shù)的DGTD方法對(duì)井下采空區(qū)地質(zhì)結(jié)構(gòu)進(jìn)行三維電磁正演，解決了傳統(tǒng)時(shí)域電磁方法在模擬曲線邊界時(shí)存在階梯效應(yīng)，對(duì)復(fù)雜、曲面目標(biāo)不能準(zhǔn)確建模、不支持高階基函數(shù)、隱式離散以及“偽解”等問(wèn)題。

2）高階疊層矢量基函數(shù)的引入提高了DGTD方法的計(jì)算精度，隨著高階疊層矢量基函數(shù)階數(shù)由0.5階、1.5階到2.5階增加，基于DGTD方法的井下采空區(qū)電磁正演計(jì)算結(jié)果的精確度分別提高了69%和74%;當(dāng)基函數(shù)為2.5階時(shí)，基于DGTD方法的井下采空區(qū)電磁正演計(jì)算結(jié)果可精準(zhǔn)得出直達(dá)波絕對(duì)幅值的大小、反射波絕對(duì)幅值的大小、反射波形的拖尾程度等表征，達(dá)到較好的電磁正演的效果。該方法為井下采空區(qū)的電磁正演提供了一種新的思路。

參考文獻(xiàn)（References）：

[1]

孫懷鳳，程銘，吳啟龍，等.瞬變電磁三維FDTD正演多分辨網(wǎng)格方法[J].地球物理學(xué)報(bào)，2018，61（12）：5096-5104.

SUN Huaifeng，CHENG Ming，WU Qilong，et al.A multiscale grid scheme in threedimensional transient electromagnetic modeling using FDTD[J].Chinese Journal of Geophysics，2018，61（12）：5096-5104.

[2]王樹奇，寧箭，王振.基于FDTD的煤層工作面上保護(hù)層檢測(cè)正演分析[J].采礦技術(shù)，2020，20（1）：145-149.

WANG Shuqi，NING Jian，WANG Zhen.Forward analysis of protective layer detection on coal seam working surface based on FDTD method[J].Mining Technology，2020，20（1）：145-149.

[3]郭文峰，曹志勇，衛(wèi)紅學(xué)，等.塌陷采空區(qū)的正演模擬及波場(chǎng)分析[J].地球物理學(xué)進(jìn)展，2015，30（2）：847-852.

GUO Wenfeng，CAO Zhiyong，WEI Hongxue，et al.The forward modeling and the analysis of wave field in the collapse of gob area[J].Progress in Geophysics，2015，30（2）：847-852.

[4]于景邨，常江浩，蘇本玉，等.老空水全空間瞬變電磁法探測(cè)三維數(shù)值模擬研究[J].煤炭科學(xué)技術(shù)，2015，43（1）：95-99，103.YU Jingcun，CHANG Jianghao，SU Benyu，et al.Study on whole space transient electromagnetic method prospect threedimensional numerical modeling of gob water[J].Coal Science and Technology，2015，43（1）：95-99，103.

[5]楊道學(xué)，趙奎，曾鵬，等.基于露天開(kāi)采地下盲空區(qū)探地雷達(dá)正演模擬及應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào)，2019，38（11）：2288-2298.

YANG Daoxue，ZHAO Kui，ZENG Peng，et al.Forward modeling and application of ground penetrating radar in blind underground cavities of opencast mining[J].Chinese Journal of Rock Mechanics and Engineering，2019，38（11）：2288-2298.

[6]張永超，李宏杰，邱浩，等.礦井瞬變電磁法的時(shí)域矢量有限元三維正演[J].煤炭學(xué)報(bào)，2019，44（8）：2361-2368.

ZHANG Yongchao，LI Hongjie，QIU Hao，et al.3D forward modeling of mine transient electromagnetic by timedomain vector finite element[J].Journal of China Coal Society，2019，44（8）：2361-2368.

[7]REN Q，SUN Q，TOBN L，et al.EB schemebased hybrid SE-FE DGTD method for multiscale EM simulations[J].IEEE Transactions on Antennas and Propagation

，2016，64（9）：4088-4091.

[8]REN Q，TOBN L E，SUN Q，et al.A new 3-D nonspurious discontinuous galerkin spectral element timedomain （DG-SETD）method for maxwells equations[J].

IEEE Transactions on Antennas and Propagation，

2015，63（6）：2585-2594.

[9]楊謙.三維DGTD若干關(guān)鍵技術(shù)研究[D].西安：西安電子科技大學(xué)，2018.

YANG Qian.Study on some key technique of 3D DGTD[D].Xian：Xidian University，2018.

[10]BAN Z G，SHI Y，YANG Q，et al.GPUaccelerated hybrid discontinuous galerkin time domain algorithm with universal matrices and local time stepping method[J].IEEE Transactions on Antennas and Propagation，

2020，68（6）：4738-4752.

[11]李林茜，魏兵，楊謙，等.二維TM波時(shí)域非連續(xù)伽略金算法理論數(shù)值通量研究[J].電波科學(xué)學(xué)報(bào)，2016，31（5）：877-882.

LI Linqian，WEI Bing，YANG Qian，et al.Study on numberical flux of node DGTD method：TM case[J].Chinese Journal of Radio Science，2016，31（5）：877-882.

[12]ZHAO L，CHEN G，YU W H.GPU accelerated discontinuous Galerkin time domain algorithm for electromagnetic problem of electrically large objects[J].

Progress Electromagnetic Rearch，2016，67：137-151.

[13]TIAN C，SHI Y，CHAN C H.An improved vector wave equationbased discontinuous galerkin time domain method and its hybridization with maxwells equationbased discontinuous galerkin time domain method[J].

IEEE Transactions on Antennas and Propagation，2018，66（11）：6170-6178.

[14]CHEN J，LIU Q H.Discontinuous galerkin timedomain methods for multiscale electromagnetic simulations：A review[J].Proceedings of the IEEE，2013，101（2）：242-253.

[15]LIU M，SIRENKO K，BAGCI H.An efficient discontinuous galerkin finite element method for highly accurate solution of maxwell equations[J].IEEE Transactions on Antennas and Propagation，2012，60（8）：3992-3998.

[16]ALVAREZ J，ANGULO L D，CABELLO M R，et al.An analysis of the leapfrog discontinuous galerkin method for maxwells equations[J].IEEE Transactions on Antennas and Propagation，2014，62（2）：197-207.

[17]NOTAROS B M.Higher order frequencydomain computational electromagnetics[J].IEEE Transactions on Antennas and Propagation，2008，56（8）：2251-2276.

[18]尹文祿.高階矢量有限元方法在電磁領(lǐng)域中的研究及應(yīng)用[D].北京：國(guó)防科學(xué)技術(shù)大學(xué)，2010.

YIN Wenlu.Study and application on higher order vector FEM in electromagnetic field[D].Beijing：National University of Defense Technology，2010.

[19]LVAREZ GONZLEZ，JESS.A discontinuous Galerkin finite element method for the timedomain solution of maxwell equations[D].Granada，Ceutay Melilla：Universidad De Granada，2014.

[20]WEBB J P.Hierarchal vector basis functions of arbitrary order for triangular and tetrahedral finite elements[J].IEEE Transactions on Antennas and Propagation，

1999，47（8）：1244-1253.

[21]WANG J，WEBB J P.Hierarchal vector boundary elements and padaption for 3-D electromagnetic scattering[J].IEEE Transactions on Antennas and Propagation，

1997，45（12）：1869-1879.