井下采空区时域不连续伽辽金法电磁正演
2021-04-18王树奇黄茸茸吴楠
王树奇 黄茸茸 吴楠
摘 要：为了获取准确的井下采空区的地理信息，利用时域不连续伽辽金法（DGTD）和高阶叠层矢量基函数的特点，建立了井下采空区三维空间模型，研究了高阶叠层矢量基函数阶数对井下采空区电磁正演精度的影响，得到井下采空区电磁场强度的空间分布规律。井下采空区电磁正演模拟结果表明，采用2.5阶叠层矢量基函数的DGTD方法计算精度比传统低阶棱边基函数提高了74%，证明了该方法对复杂目标电磁特性研究的实用性，并能准确可靠地实现井下采空区地质结构的电磁探测。关键词：井下采空区;时域不连续伽辽金法;高阶叠层矢量基函数;电磁正演中图分类号：TM 153
文献标志码：A
文章编号：1672-9315（2021）02-0348-07
DOI：10.13800/j.cnki.xakjdxxb.2021.0220开放科学（资源服务）标识码（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 引 言
煤炭是中國能源支柱产业，然而大范围的煤炭开采在矿山中留下了大量的采空区，导致地下岩体原有的力学平衡被打破，扭曲的岩体随时可能发生位移、岩爆等事故;更严重的是采空区会被瓦斯、地下水等充填，给矿区的开采作业埋下重大的安全隐患，采用电磁数值计算对矿井复杂地质结构进行探测分析具有重要的意义。
国内外学者针对时域电磁正演的方法做了大量研究，但对于井下采空区的时域电磁正演的研究多采用时域有限差分（finite difference time domain，FDTD）方法和时域矢量有限元（finite element time domain，FETD）[1-3]。于景邨和常江浩采用三维有限差分法对老空水的矿井瞬变电磁响应进行了研究[4]。杨道学提出了基于卷积完全匹配层的交替方向隐式的FDTD方法，并将其应用于探地雷达正演[5]。FDTD方法具有理论难度低，计算效率高等优点，但是该方法中使用的六面体剖分单元在复杂几何结构中具有局限性，且不支持高阶基函数等问题将直接影响计算精度。张永超和李宏杰等研究了时域矢量有限元三维正演[6]，拓展了矿井瞬变电磁正演对复杂地质模型的适用性，棱边矢量基函数的使用很好地解决了节点方法中的“伪解”，但只是用了低阶棱边基函数且隐式离散格式，计算效率低。时域不连续伽辽金方法（discontinuous galerkin time domain，DGTD）[7-12] 继承了FETD 采用非结构网格对复杂几何外形拟合好，便于使用高阶基函数的优点;该算法对Maxwell 方程采用伽辽金加权法得到弱解形式，放松了单元之间的边界条件，单元之间通过数值通量联系交换数据;高阶叠层矢量基函数的引入是一个突破，它满足单元边界上的 场切向连续性，消除了节点法中存在的“伪解”并且可显著提高该方法的计算精度;在时间离散方面既支持显式时间格式也支持隐式时间格式;在空间离散方面由于单元间的高度独立性，对于每个单元可以选用不同阶的基函数进行离散，需要高精度时采用高阶基函数，对精度要求低时采用低阶基函数。可以将大型多尺度问题被拆分成为一个个相对独立的小问题，使得分析问题更为简单方便，以上原因更加说明了研究DGTD算法的价值与意义。文中首次将基于高阶叠层矢量基函数的DGTD方法应用于井下采空区三维正演研究，建立了井下采空区三维空间模型，分析了叠层矢量基函数阶数对正演结果的影响，试验验证了井下采空区模型正演的计算精度。
1 DGTD方法
三维非均匀各向同性介质Maxwell旋度方程为
μ
H
t
=-×E-σm
H
ε
E
t
=×H-σe
E
（1）
式中 ε为介电系数;μ为磁导系数;σe为电导率;σm为导磁率;
E
和H为电场强度与磁场强度。对公式（1）采用伽辽金加权法，在四面体单元内进行体积分，可以得到以下弱解形式
∫T m
μ
H
t
+×E+σm
H
·q′dv=0
∫T m
ε
E
t
-×H+σe
E
·q′dv=0
（2）
q′为加权函数，使用以下恒等式对公式（2）中的
×E和
×H项进行代数变换
∫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）在相邻区域单元交界面上采用数值通量[13-14]的形式来保证单元之间切向场的连续性，数值通量格式如下
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*称为數值通量，
定义在2个相邻单元的交界面上，由2个相接面的电磁场合成，数值通量在放宽单元之间连续性边界条件的基础上，实现了相邻单元之间场量信息的交互;
Em，
Hm
为当前四面体单元的电场磁矢量;
Em+，
Hm+
为相邻四面体单元的电磁场矢量。
将数值通量代入公式（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）的空间部分采用基函数在每个单元中展开
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）在得到空间离散格式之后将进行时间离散，最后根据时间离散格式完成时间迭代。可采取不同方式对得到的空间离散方程进行时间偏导数的离散，若采用隐式时间离散方案则将降低计算效率，文中采用已经广泛应用于时域算法中的蛙跃算法进行时间离散[15]，蛙跃算法是一种显式时域离散方法，电磁场相隔半个时间步交替进行迭代计算。规定电场E在整数时间步tn采样，磁场H在半整数步tn+1/2进行采样。在公式（7）中对时间的导数采取中心差分的格式处理，由二阶中心差分代替一阶时间导数[16]
（dtHm）n
=Hmn+1/2-
Hmn-1/2
Δt
+O（Δt2）
（dtEm）n+1/2
=Emn+1-
Emn
Δt
+O（Δt2）
（8）
由于在迭代过程中会出现电场的半整数步、磁场的整数步，对其取平均近似[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）
其中系数为
α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 高阶叠层矢量基函数为了克服低阶矢量基函数精度较低的缺点，NOTAROS等提出了高阶插值型矢量基函数[17-19]。高阶插值型矢量基函数具有线性独立性好、物理解释明确、编程实现方便等优点。但在一个四面体内插值型矢量基函数具有更多的未知量，为了满足交界面处切向场连续，插值矢量基不允许不同阶次基函数的混合，这导致了在原本低阶单元即可足够精确描述的区域引入了大量冗余自由度，从而增加了存储量和计算时间，降低了算法分析的效率。高阶叠层型矢量基函数[20-21]的提出正是为了解决高阶插值型矢量基函数的这些缺点。首先，高阶叠层型基函数具有叠层嵌套的特点，高阶叠层构矢量基函数的构建是基于低阶基函数逐阶递增的，即高阶的叠层基函数里面包含有阶数较低的基函数，这在保证高精度的情况下降低了基函数构造难度以及存储量。其次，借助高阶叠层矢量基函数，时域不连续伽辽金算法可以实现计算域内高阶和低阶单元的混合计算，数值结果表明高阶叠层矢量基函数的引用使得对目标离散时采用较大尺寸单元得到的結果与采用较小尺寸单元剖分的低阶基函数方案有相当的精度。叠层矢量基函数是一种基于棱边、面和体积的基函数，它将自由度赋予棱边、面和体积而不是单元节点。它隐含了散度为零的边界条件，消除了伪解，非常适合用来表示矢量场。综上可知，高阶叠层矢量基函数在处理大规模电磁问题和多尺度电磁问题上都具有非常大的优势，具有很高的实用性与灵活性。表1给出了文中所采用的基于四面体单元的0.5阶到2.5阶叠层矢量基函数形式。当然，叠层型矢量基函数还有更高的阶的形式，但随着基函数阶数的增加，自由度会以数倍增加，这将增加计算难度和时间成本。并且当基函数阶数为2.5阶时，该算法就可以获得较好的计算精度。因此文中就0.5阶至2.5阶基函数展开研究。
3 矩形谐振腔模拟计算为了分析不同高阶叠层矢量基函数的DGTD方法的计算精度，建立了尺寸为2.0 m×0.3 m×0.2 m的无损耗矩形PEC边界谐振腔模型，如图2、图3所示。
图2中腔体被划分为2 261个均匀网格尺寸的四面体，图3所示计算域被划分为2 292个网格大
小不同的四面体单元。在点（0，0，0）m处激发
一个中心频率为310 MHz，极化方向为（1，1，1）的布莱克曼-哈里斯（blackman harris window，BHW）电偶极子源。
分别选取0.5阶、1.5阶、2.5阶以及混合阶叠层矢量基函数时，DGTD方法在观测点（0.7，0，0）m处得到的时域波形分别如图4（a）～图4（d）所示。
从图4可以看出，当叠层基函数为0.5阶时，其结果的计算精度较差。随着高阶叠层矢量基函数的阶数增加，数值计算结果的精度有了显著提升，经分析见表2所示结果。当基函数阶数分别为0.5阶、1.5阶、2.5阶以及混合阶叠层矢量基函数时，计算结果的相对误差为43.94%、2244%、229%和2.42%。从0.5阶到2.5阶，DGTD方法的计算相对误差降低了4165%，计算精度得到了明显的改善;混合阶基函数的计算中基函数的阶数根据四面体单元的尺寸进行选择，其结果的相对误差与2.5阶基函数结果相近，该结果表明基于
高阶叠层矢量基函数的DGTD方法在保证计算精度的同时可实现计算域内高阶和低阶单元的混合计算，对复杂目标的电磁特性研究具有很高的实用性。
4 井下采空区电磁正演建立煤层区域为10 m×10 m×10 m，井下采空区尺寸为4 m×4 m×2 m的三维正演地质模型，如图5，图6所示。背景设置为均匀介质，取相对介电常数为4.09褐煤进行模拟验证。以空间模型的中心位置作为原点，掘进工作面位于z方向上。在（0 m，0 m，4 m）处设置中心频率为31 MHz的BHW电偶极子激励源。接收点放置于掘进工作面上，其位置坐标为（0 m，1 m，5 m）。
为了研究所选取的基函数对井下电磁响应特性的影响，设计了地质模型，分别将0.5阶、1.5阶、2.5阶叠层矢量基函数应用于电磁正演中，结果如图7所示。对比不同阶基函数对采空区的电磁响应特征的影响，当基函数为0.5阶、1.5阶、25阶3种不同阶数时均可分辨反射波信号。从图7所示的基于DGTD方法井下采空区电场响应结果可看出，随着基函数阶数的增加，直达波和反射波的幅值精确度有了显著提高，且采空区反射波曲线拖延现象得到了显著改善。对比DGTD方法中不同阶叠层基函数正演波形，误差分析结果见表3。计算中将接收点的场值记录并与仿真软件结果进行比较，其相对误差定义为
Error=20log（|E-Eref|/Eref max）（dB）
式中 E为时域波形数值解;Eref为仿真软件计算结果。
经分析得：井下采空区正演的相对误差从0.5阶的-19.12 dB缩小至2.5階的-73.19 dB，计算精度提高了74%。由以上分析可得高阶叠层矢量基函数的阶数将直接影响DGTD方法的正演计算精度，采空区正演时域结果随基函数阶数的增加，幅值的精确度大幅提升。
5 结 论
1）利用基于高阶叠层矢量基函数的DGTD方法对井下采空区地质结构进行三维电磁正演，解决了传统时域电磁方法在模拟曲线边界时存在阶梯效应，对复杂、曲面目标不能准确建模、不支持高阶基函数、隐式离散以及“伪解”等问题。
2）高阶叠层矢量基函数的引入提高了DGTD方法的计算精度，随着高阶叠层矢量基函数阶数由0.5阶、1.5阶到2.5阶增加，基于DGTD方法的井下采空区电磁正演计算结果的精确度分别提高了69%和74%;当基函数为2.5阶时，基于DGTD方法的井下采空区电磁正演计算结果可精准得出直达波绝对幅值的大小、反射波绝对幅值的大小、反射波形的拖尾程度等表征，达到较好的电磁正演的效果。该方法为井下采空区的电磁正演提供了一种新的思路。
参考文献（References）：
[1]
孙怀凤，程铭，吴启龙，等.瞬变电磁三维FDTD正演多分辨网格方法[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].采矿技术，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]郭文峰，曹志勇，卫红学，等.塌陷采空区的正演模拟及波场分析[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]于景邨，常江浩，苏本玉，等.老空水全空间瞬变电磁法探测三维数值模拟研究[J].煤炭科学技术，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]杨道学，赵奎，曾鹏，等.基于露天开采地下盲空区探地雷达正演模拟及应用[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若干关键技术研究[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波时域非连续伽略金算法理论数值通量研究[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]尹文禄.高阶矢量有限元方法在电磁领域中的研究及应用[D].北京：国防科学技术大学，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.
