井下采空区时域不连续伽辽金法电磁正演

王树奇，黄茸茸，吴 楠

(1.西安科技大学 通信与信息工程学院，陕西 西安 710054;2.陕西学前师范学院 信息工程学院，陕西 西安 710100)

0 引 言

煤炭是中国能源支柱产业，然而大范围的煤炭开采在矿山中留下了大量的采空区，导致地下岩体原有的力学平衡被打破，扭曲的岩体随时可能发生位移、岩爆等事故；更严重的是采空区会被瓦斯、地下水等充填，给矿区的开采作业埋下重大的安全隐患，采用电磁数值计算对矿井复杂地质结构进行探测分析具有重要的意义。

国内外学者针对时域电磁正演的方法做了大量研究，但对于井下采空区的时域电磁正演的研究多采用时域有限差分(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旋度方程为

(1)

式中ε为介电系数；μ为磁导系数；σe为电导率；σm为导磁率；E和H为电场强度与磁场强度。

对公式(1)采用伽辽金加权法，在四面体单元内进行体积分，可以得到以下弱解形式

(2)

φq′为加权函数，使用以下恒等式对公式(2)中的×E和×H项进行代数变换

(3)

在相邻区域单元交界面上采用数值通量[13-14]的形式来保证单元之间切向场的连续性，数值通量格式如下

对上述数学关系式 y=134.1x2-182.8x+385.3进行一阶求导，确定出在定义域范围内日耗电最低时的抽油机悬点载荷利用率为68.2%，称之为第一抽油机悬点载荷利用率。按照类似的方法对所得的40组数据分别进行分析，得出40组在各自电动机负载率下的最优抽油机悬点载荷利用率，具体数据如表1。

(4)

图1 数值通量Fig.1 Numerical flux

将数值通量代入公式(3)中可得

(5)

公式(5)的空间部分采用基函数在每个单元中展开

(6)

整理简化后得到基于四面体单元的DGTD显式半离散方程组

(7)

在得到空间离散格式之后将进行时间离散，最后根据时间离散格式完成时间迭代。可采取不同方式对得到的空间离散方程进行时间偏导数的离散，若采用隐式时间离散方案则将降低计算效率，文中采用已经广泛应用于时域算法中的蛙跃算法进行时间离散[15]，蛙跃算法是一种显式时域离散方法，电磁场相隔半个时间步交替进行迭代计算。规定电场E在整数时间步tn采样，磁场H在半整数步tn+1/2进行采样。

在公式(7)中对时间的导数采取中心差分的格式处理，由二阶中心差分代替一阶时间导数[16]

(8)

由于在迭代过程中会出现电场的半整数步、磁场的整数步，对其取平均近似[16]有以下形式

(9)

(10)

其中系数为

(11)

2 高阶叠层矢量基函数

为了克服低阶矢量基函数精度较低的缺点，NOTAROS等提出了高阶插值型矢量基函数[17-19]。高阶插值型矢量基函数具有线性独立性好、物理解释明确、编程实现方便等优点。但在一个四面体内插值型矢量基函数具有更多的未知量，为了满足交界面处切向场连续，插值矢量基不允许不同阶次基函数的混合，这导致了在原本低阶单元即可足够精确描述的区域引入了大量冗余自由度，从而增加了存储量和计算时间，降低了算法分析的效率。

高阶叠层型矢量基函数[20-21]的提出正是为了解决高阶插值型矢量基函数的这些缺点。首先，高阶叠层型基函数具有叠层嵌套的特点，高阶叠层构矢量基函数的构建是基于低阶基函数逐阶递增的，即高阶的叠层基函数里面包含有阶数较低的基函数，这在保证高精度的情况下降低了基函数构造难度以及存储量。其次，借助高阶叠层矢量基函数，时域不连续伽辽金算法可以实现计算域内高阶和低阶单元的混合计算，数值结果表明高阶叠层矢量基函数的引用使得对目标离散时采用较大尺寸单元得到的结果与采用较小尺寸单元剖分的低阶基函数方案有相当的精度。叠层矢量基函数是一种基于棱边、面和体积的基函数，它将自由度赋予棱边、面和体积而不是单元节点。它隐含了散度为零的边界条件，消除了伪解，非常适合用来表示矢量场。

综上可知,高阶叠层矢量基函数在处理大规模电磁问题和多尺度电磁问题上都具有非常大的优势，具有很高的实用性与灵活性。表1给出了文中所采用的基于四面体单元的0.5阶到2.5阶叠层矢量基函数形式。当然，叠层型矢量基函数还有更高的阶的形式，但随着基函数阶数的增加，自由度会以数倍增加，这将增加计算难度和时间成本。并且当基函数阶数为2.5阶时，该算法就可以获得较好的计算精度。因此文中就0.5阶至2.5阶基函数展开研究。

表1 叠层矢量基函数Table 1 Hierarchical vector basis functions

3 矩形谐振腔模拟计算

为了分析不同高阶叠层矢量基函数的DGTD方法的计算精度，建立了尺寸为2.0 m×0.3 m×0.2 m的无损耗矩形PEC边界谐振腔模型，如图2、图3所示。

图2 矩形PEC边界均匀单元尺寸谐振腔Fig.2 Uniformly meshed rectangular resonator with PEC boundary

图2中腔体被划分为2 261个均匀网格尺寸的四面体，图3所示计算域被划分为2 292个网格大小不同的四面体单元。在点(0，0，0)m处激发一个中心频率为310 MHz，极化方向为(1，1，1)的布莱克曼-哈里斯(blackman harris window，BHW)电偶极子源。

图3 矩形PEC边界不同单元尺寸边界谐振腔Fig.3 Non-uniformly meshed rectangular resonator with PEC boundary

分别选取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%、22.44%、2.29%和2.42%。从0.5阶到2.5阶，DGTD方法的计算相对误差降低了41.65%，计算精度得到了明显的改善；混合阶基函数的计算中基函数的阶数根据四面体单元的尺寸进行选择，其结果的相对误差与2.5阶基函数结果相近，该结果表明基于高阶叠层矢量基函数的DGTD方法在保证计算精度的同时可实现计算域内高阶和低阶单元的混合计算，对复杂目标的电磁特性研究具有很高的实用性。

图4 基于叠层基函数的DGTD方法接收点数值结果Fig.4 Numerical results of the receiving point of the DGTD method based on the hierarchical vector basis functions

表2 矩形谐振腔模拟计算误差分析Table 2 Error analysis of simulation of rectangular resonator

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)。

图5 井下采空区三维模型Fig.5 Three-dimensional model of underground goaf

图6 井下采空区三维模型剖分Fig.6 Meshes of three-dimensional model of underground goaf

为了研究所选取的基函数对井下电磁响应特性的影响，设计了地质模型，分别将0.5阶、1.5阶、2.5阶叠层矢量基函数应用于电磁正演中，结果如图7所示。对比不同阶基函数对采空区的电磁响应特征的影响，当基函数为0.5阶、1.5阶、2.5阶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方法的正演计算精度，采空区正演时域结果随基函数阶数的增加，幅值的精确度大幅提升。

图7 基于DGTD方法煤层采空区正演结果Fig.7 Forward modeling results of coal seam goaf based on DGTD method

表3 井下采空区电磁正演误差分析Table 3 Error analysis of electromagnetic forward modeling results in underground goaf

5 结 论

1)利用基于高阶叠层矢量基函数的DGTD方法对井下采空区地质结构进行三维电磁正演，解决了传统时域电磁方法在模拟曲线边界时存在阶梯效应，对复杂、曲面目标不能准确建模、不支持高阶基函数、隐式离散以及“伪解”等问题。

2)高阶叠层矢量基函数的引入提高了DGTD方法的计算精度，随着高阶叠层矢量基函数阶数由0.5阶、1.5阶到2.5阶增加，基于DGTD方法的井下采空区电磁正演计算结果的精确度分别提高了69%和74%；当基函数为2.5阶时，基于DGTD方法的井下采空区电磁正演计算结果可精准得出直达波绝对幅值的大小、反射波绝对幅值的大小、反射波形的拖尾程度等表征，达到较好的电磁正演的效果。该方法为井下采空区的电磁正演提供了一种新的思路。