王樹奇 黃茸茸 吳楠

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

文獻標志碼：A

文章編號：1672-9315（2021）02-0348-07

DOI：10.13800/j.cnki.xakjdxxb.2021.0220開放科學（資源服務(wù)）標識碼（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 引 言

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

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

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為電場強度與磁場強度。對公式（1）采用伽遼金加權(quán)法，在四面體單元內(nèi)進行體積分，可以得到以下弱解形式

∫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ù)，使用以下恒等式對公式（2）中的

×E和

×H項進行代數(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]的形式來保證單元之間切向場的連續(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個相鄰單元的交界面上，由2個相接面的電磁場合成，數(shù)值通量在放寬單元之間連續(xù)性邊界條件的基礎(chǔ)上，實現(xiàn)了相鄰單元之間場量信息的交互;

Em，

Hm

為當前四面體單元的電場磁矢量;

Em+，

Hm+

為相鄰四面體單元的電磁場矢量。

將數(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ù)在每個單元中展開

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

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

（6）整理簡化后得到基于四面體單元的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ù)時間離散格式完成時間迭代。可采取不同方式對得到的空間離散方程進行時間偏導(dǎo)數(shù)的離散，若采用隱式時間離散方案則將降低計算效率，文中采用已經(jīng)廣泛應(yīng)用于時域算法中的蛙躍算法進行時間離散[15]，蛙躍算法是一種顯式時域離散方法，電磁場相隔半個時間步交替進行迭代計算。規(guī)定電場E在整數(shù)時間步tn采樣，磁場H在半整數(shù)步tn+1/2進行采樣。在公式（7）中對時間的導(dǎo)數(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）

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

Hmn

=Hmn+1/2+

Hmn-1/2

2

+O（Δt2）

Emn+1/2

=Emn+1+

Emn

2

+O（Δt2）

（9）

通過初始的激勵源與邊界條件，就可以利用DGTD逐步迭代求得計算域中電磁場分布隨時間的變化

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

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

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

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

一個中心頻率為310 MHz，極化方向為（1，1，1）的布萊克曼-哈里斯（blackman harris window，BHW）電偶極子源。

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

參考文獻（References）：

[1]

孫懷鳳，程銘，吳啟龍，等.瞬變電磁三維FDTD正演多分辨網(wǎng)格方法[J].地球物理學報，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的煤層工作面上保護層檢測正演分析[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)紅學，等.塌陷采空區(qū)的正演模擬及波場分析[J].地球物理學進展，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]于景邨，常江浩，蘇本玉，等.老空水全空間瞬變電磁法探測三維數(shù)值模擬研究[J].煤炭科學技術(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]楊道學，趙奎，曾鵬，等.基于露天開采地下盲空區(qū)探地雷達正演模擬及應(yīng)用[J].巖石力學與工程學報，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]張永超，李宏杰，邱浩，等.礦井瞬變電磁法的時域矢量有限元三維正演[J].煤炭學報，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].西安：西安電子科技大學，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波時域非連續(xù)伽略金算法理論數(shù)值通量研究[J].電波科學學報，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].北京：國防科學技術(shù)大學，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.