城市地下空洞等值反磁通瞬变电磁三维正演模拟及分析

Simulation and Analysis of 3D Forward Modeling of Opposing Coils Transient Electromagnetic Field in Urban Underground Cavity

摘要: 等值反磁通瞬变电磁法(Opposing Coils TEM，简称OCTEM)由于具有浅层盲区小、施工效率高等优点，目前在城市浅层地球物理勘探领域受到了广泛的应用。为研究OCTEM在城市地下空洞探测中的应用，基于等值反磁通理论与矢量有限元瞬变电磁三维正演算法，构建低阻异常模型、高阻异常模型，分析了不同填充状态情况下的空洞OCTEM响应规律。并遵循相似性原则设计进行物理实验，对比数值计算和物理实验结果论证了数值模拟结果的正确性。结果表明：OCTEM对高阻异常与低阻异常都具有良好的识别能力，但对于低阻异常体OCTEM更为敏感。研究结果可为OCTEM在实际空洞探测应用中提供参考。

Abstract: The opposing coils transient electromagnetic method (OCTEM) has been widely used in the field of urban shallow geophysical exploration due to its advantages of small shallow blind area and high efficiency. In order to study the application of OCTEM in urban underground cavity detection, based on the opposing coils theory and vector finite element transient electromagnetic 3D forward algorithm, low resistance anomaly model and high resistance anomaly model were constructed, and the OCTEM response law of cavity under different filling conditions was analyzed. The physical experiment was designed according to the similarity principle, and the correctness of the numerical simulation results was verified by comparing the numerical calculation and the physical experiment results. The results show that OCTEM has good recognition ability for both high and low resistance anomalies, but is more sensitive to low resistance anomalies. The results can provide a reference for the application of OCTEM in cavity detection.

1. 前言

城市地下空洞是常见的地质灾害之一，潜伏地下的空洞会造成地基土软弱，承载力不够，进而导致城市道路地面发生沉降、塌陷事故的发生，未经发现的地下空洞如定时炸弹一般严重危害着人民群众人民财产安全。因此，对城市地下空洞进行探测，查明空洞的分布范围、规模并进行及时的治理意义重大。等值反磁通瞬变电磁法(Opposing Coils TEM，简称OCTEM)相比于传统瞬变电磁法具有浅层盲区小、横向探测精度高、施工简便等优点，目前在城市浅层地质地球物理勘探领域受到了广泛应用 [1] [2] [3]。在实际的城市地下空洞探测应用中，由于地下介质的复杂性，空洞可能表现为空气填充的高阻空洞或者水–泥填充的低阻溶洞，如何区分地质体的准确电性情况而给地质勘探造成极大的困难。为提高OCTEM空洞探测的研究与应用水平，对城市地下空洞开展OCTEM三维正演模拟。

目前三维瞬变电磁数值模拟主要采用有限差分法、有限元法或积分方程法 [4] [5] [6]。早在1980年，Kuo和Cho [7] 在二维有限元瞬变电磁正演算法的基础上，对三维矢量有限元瞬变电磁正演算法开展了研究。1987年，Gupta [8] [9] 等首次将积分方程法引入矢量有限元瞬变电磁3D数值模拟中，降低稀疏矩阵维数，减少稀疏矩阵的计算量，提高了计算速度，并对三维低阻异常模型电磁响应进行了模拟。1994年，殷长春 [10] 结合并矢格林函数理论和积分方程法研究了三维地质体的频率电磁响应。1999年唐新功、胡文宝 [11] 等人通过在体积分方程法引入张量格林函数对多种三维薄板的瞬变电磁响应开展了研究。2008年，Borner [12] 等人提出了一种适用于瞬态偶极子源电磁场模拟的三维有限元瞬变电磁算法。2012年，李健慧 [13] 等人利用矢量有限元法实现了基于回线源的瞬变电磁3D数值模拟。2014年，Chung [14] 等人提出了一种全域三维有限元算法，并计算了三维瞬变电磁响应。2016年，李瑞雪 [15] 等人为解决积分方程–矢量有限元瞬变电磁三维正演模拟中存在的伪解问题，增强了边界条件的连续性，同时采用Flow-through Hankel transform快速计算技术，提高计算速度，并分析了三维异常体和三维薄板模型的电磁响应特征。2018年，马炳镇 [16] 等人研究了起伏地形下地面瞬变电磁响应规律。2019年，张永超 [17] 等人基于矢量有限元理论提出了矿井瞬变电磁法三维正演算法，探究了关断时间与关断波形的影响，并分析总结了立方体形采空区的电磁响应规律特征。

本文在矢量有限元瞬变电磁三维正演算法的基础上构建反向对偶磁源，建立OCTEM三维正演算法，通过构建三维长方体空洞模型，改变空洞异常体的电阻率，研究在不同填充状态下空洞的OCTEM电磁响应规律，总结响应特征。并遵循相似性原则建立物理实验，通过物理实验验证了城市地下空洞OCTEM三维数值模拟结果的正确性。

2. 三维正演算法理论

2.1. 控制方程

在各向同性的三维均匀介质中，在准静态条件(忽略位移电流)下，Maxwell方程组有如下微分形式 [18]：

$\nabla \times E+\mu \frac{\partial H}{\partial t}=0$ (1-1)

$\nabla \times H-\sigma E=J_s$ (1-2)

$\nabla \cdot B=0$ (1-3)

$\nabla \cdot E=\rho$ (1-4)

公式中：E为电场强度；B为磁感应强度；H为磁场强度；D为电位移矢量；J为电流密度； $J_s$ 为源电流密度； $\rho$ 为自由电荷密度， $\mu$ 为磁导率， $\sigma$ 为电导率。

则当取时间因子为 ${\text{e}}^{i\omega t}$ 时，公式(1-1)可表示为：

$\nabla \times E+i\omega \mu H=0$ (2)

然后将公式(2)的两边同时取旋度可得

$\nabla \times \nabla \times E+i\omega \mu \nabla \times H=0$ (3)

再将式(3)代入式(1-2)中得电场的矢量波动方程为

$\nabla \times \nabla \times E+i\omega \mu \sigma E=-i\omega \mu J_s$ (4)

上式中的电场是包含一次电场(均匀半空间电场)与二次电场(异常体激发的二次电场)的总场，然而在直接求解总场的过程中，由于场源具有空间奇异性，在模拟场源附近电场快速变化时需要对场源附近区域做超细化的剖分，增加计算量。所以本文选择先求取一次电场与二次电场，再相加得到总场。将总场分解为一次电场 $E^P$ 与二次电场 $E^s$，则有

$E=E^P+E^s$ (5)

将式(5)带入式(4)，电场矢量波动方程分解整理后得

$\nabla \times \nabla \times E^P+i\omega \mu \sigma E^P=-i\omega \mu J_s$ (6)

$\nabla \times \nabla \times E^s+i\omega \mu \sigma E^s=-i\omega \mu \Delta \sigma E^P$ (7)

式(7)为二次电场的扩散方程，其中 $\Delta \sigma$ 表示电导率异常，是异常体电导率与背景电导率的差值，则二次电场扩散方程的残数可表示为：

$R=\nabla \times \nabla \times E^s+i\omega \mu \sigma E^s+i\omega \mu \Delta \sigma E^P$ (8)

而对于一次场 $E^P$，直接采用Nabighian [19] 的推导结果对均匀半空间模型求解：

$E\left(\rho ,z\right)=\frac{-Ia}{\sigma }{\displaystyle {\int }_0^{\infty }\left(u_1-\lambda \right){\text{e}}^{-u_{1}z}{\text{e}}^{-\lambda h}\lambda J_1\left(\lambda a\right)J_1\left(\lambda \rho \right)\text{d}\lambda }$ (9)

其中，h为线圈距离地面的高度， $u_1=\sqrt{{\lambda }^2+i\omega {\mu }_0\sigma }$。

由于本文要模拟的是等值反磁通装置， $E^P$ 可视为上下两个反向磁源单独作用所引起响应的矢量叠加所得，即

$E^P=E_{Tx(+)}+E_{Tx(-)}$ (10)

当h = 2d时，即单独上方反向电流环发射时：

$E_{Tx(-)}=\frac{-Ia}{\sigma }{\displaystyle {\int }_0^{\infty }\left(u_1-\lambda \right){\text{e}}^{-u_{1}z}{\text{e}}^{-2\lambda d}\lambda J_1\left(\lambda a\right)J_1\left(\lambda \rho \right)\text{d}\lambda }$ (11)

当h = 0时，即地表单独正向电流环发射时：

$E_{Tx(+)}=\frac{-Ia}{\sigma }{\displaystyle {\int }_0^{\infty }\left(u_1-\lambda \right){\text{e}}^{-u_{1}z}\lambda J_1\left(\lambda a\right)J_1\left(\lambda \rho \right)\text{d}\lambda }$ (12)

在无源区域得介质分界面上，电场E的切向分量具有连续性，即：

$n\times \left(E_1-E_2\right)=0$ (13)

其中， $E_1$ 与 $E_2$ 分别为介质分界面两边的电场。

而在无穷远边界上，电场满足Dirichlet边界条件，即电场在无穷远边界上的切向分量为零：

$n\times E=0$ (14)

2.2. 矢量有限单元分析

本文采用六面体单元进行矢量单元分析。单元各棱边的编号规则如图1所示。

Figure 1. Hexahedron vector element

图1. 六面体矢量单元

则在每个六面体单元内的二次电场可由十二条棱边的矢量基函数表示为

$E^s={\displaystyle \underset{i=1}{\overset{12}{\sum }}{N}_i{E}_i^s}$ (15)

$E_i^s$ 代表棱边上的二次电场值，由于矢量基函数的散度为零，即满足：

$\nabla \cdot N_i=\frac{\alpha N_{ix}}{\alpha x}+\frac{\alpha N_{iy}}{\alpha y}+\frac{\alpha N_{iz}}{\alpha z}=0$ (16)

将式(16)两边同时取旋度，并将式(17)带入可得：

$\nabla \cdot E^s={\displaystyle \underset{i=1}{\overset{12}{\sum }}\left(\nabla \cdot N_i\right)E_i^s}=0$ (17)

式(17)可证明由矢量基函数所表示的二次电场的散度为0，保证二次电场解的有效性。

将式(16)代入到式(8)并引入Galerkin加权余量法，单元内的加权余量积分为

$R^e=\left[A^e\right]\left\{E^s\right\}+\left[B^e\right]\left\{E^s\right\}+\left[C^e\right]\left\{E^P\right\}$ (18)

其中， $\left[A^e\right]$ 、 $\left[B^e\right]$ 、 $\left[C^e\right]$ 均为12 × 12阶的单元矩阵：

$\left[A^e\right]={\displaystyle {\iiint }_{V^e}\left(\nabla \times N_i^e\right)\cdot \left(\nabla \times N_j^e\right)\text{d}V}$ (19)

$\left[B^e\right]={\displaystyle {\iiint }_{V^e}i\omega \mu \sigma N_i^e\cdot N_j^e\text{d}V}$ (20)

$\left[C^e\right]={\displaystyle {\iiint }_{V^e}i\omega \mu \Delta \sigma N_i^e\cdot N_j^e\text{d}V}$ (21)

将上述所有的单元矩阵按照编码规则进行组合，并令所有单元内加权余量积分为零，得到总体矩阵方程组为

$K{E}^s=-S$ (22)

其中，K为 $N\times N$ 阶的稀疏矩阵， $E^s$ 、S为N阶向量，N是模型剖分后的总棱边数。

采用CG算法求解方程组(22)得到二次电场，带入式(5)获得总场，再通过法拉第定律求得总磁场：

$H=-\frac{1}{i\omega \mu }{\displaystyle \underset{i=1}{\overset{12}{\sum }}{E}_i\nabla \times N_i}$ (23)

在解出频率域内的总磁场后，然后通过Fourier变换将总磁场转换到时间域。

3. 数值模拟

为研究无充填与全充填空洞的等值反磁通瞬变电磁响应特征，建立如图2所示的三维地电模型，模型区域大小为1000 m × 1000 m × 800 m (长 × 宽 × 高)，其中长、宽、高分别对应x、y、z方向，模型采用矩形单元进行网格剖分，并对异常边界加密处理，矩形单元尺寸随着与发射源距离的增加而逐渐增大，共剖分41 × 41 × 26个单元。空洞表示为大小100 m × 100 m × 50 m (长 × 宽 × 高)的块状异常体，空洞的中心投影点位于地表坐标(500，500)处，埋深为30 m。正向发射线圈TX(+)与反向发射线圈TX(−)均采用半径为a = 10 m的圆形回线，正向发射线圈TX(+)水平置于地表中心h = 0 m处，反向发射线圈TX(−)置于地表中心上方h = 4 m处，发射电流大小均为10 A；在发射的间歇期，均在h = 2 m处采用等效面积为100 m²的中心回线装置进行观测，观测时间为10⁻⁵至10⁻³秒，对数等距分为20个时间门。测线沿东西方向，线距25 m，共13条测线，每条测线13个测点，点距25 m，共计169个测点，主测线中心测点位于异常体中心，共计81个测点，低阻模型总计算耗时103043.32 s，平均耗时609.7 s。

为更立体、更直观的观测低阻空洞与高阻空洞的OCTEM电磁响应特征，将各个测点相同时间门的电磁响应绘制成电磁响应时间切片。低阻模型的不同时间道电磁响应切片分别如图3所示。由图3可知，在低阻空洞模型，早期在远离异常体的区域，观测点相同时间道的响应幅值均保持在较低的水平，整体无较大变化，在异常体边界区域，响应幅值开始出现明显的增大，进入异常内部区域后，响应幅值继续增加，并在达到一定水平后趋于平缓，形成类似梯形的异常响应峰值平台，非常好的反映了异常体的形态。而在进入中晚期后，随着观测信号深度的加深，空洞异常的影响减小，信号幅值衰减，响应幅值峰值平台消失，整体形态渐渐趋于向圆形变化。

Figure 2. Diagram of three-dimensional cube cavity model

图2. 三维立方体空洞模型示意图

Figure 3. OCTEM electromagnetic response time slice of low resistance cavity model. (a) 0.008 ms; (b) 0.165 ms; (c) 0.338 ms

图3. 低阻空洞模型数值模拟OCTEM 电磁响应时间切片。(a) 0.008 ms；(b) 0.165 ms；(c) 0.338 ms

高阻模型的不同时间道电磁响应切片分别如图4所示。由图4可知，与低阻空洞模型一样，在远离异常体的区域，观测点相同时间道的响应幅值不大，都近似于背景响应，在异常体边界至内部，由于高阻空洞对涡流热损耗较强，电磁响应信号衰减更加强烈，异常体范围内的响应赋值低于周围的水平，形成凹槽，与低阻空洞模型形成鲜明的对比。

Figure 4. OCTEM electromagnetic response time slice of high resistance cavity model. (a) 0.01 ms; (b) 0.025 ms; (c) 0.06 ms

图4. 高阻空洞模型数值模拟OCTEM电磁响应时间切片。(a) 0.01 ms；(b) 0.025 ms；(c) 0.06 ms

4. 物理模拟

为了验证数值模拟结果的准确性，遵循物理模拟的相似性远原则，设计建立相似比为1:1000的空洞物理模型并进行实验。实验采用大小为100 cm × 100 cm × 80 cm水槽模型模拟道路下的半空间，为更好的模拟地下的真实情况，注入沙土充当背景介质。本次实验仪器采用湖南五维地质科技有限公司为坑道探测所研发的ADTEM-18高精度瞬变电磁系统。发射线圈与接收线圈均采用半径为1 cm的圆形线圈，发射电流为4.2 A，发射频率25 Hz。采用实验共设置13条观测线，线间距2.5 cm，每条线设置13个观测点，点距为2.5 cm。低阻模型与高阻模型的物理模拟的主测线OCTEM电磁响应多测道剖面图如图5、图6所示，电磁响应时间切片分别由图7、图8所示，由图7、图8可知，采集数据存在一定干扰，与数值模拟结果存在一定差别，但总体上物理模拟结果是符合数值模拟的响应规律的。

Figure 5. Multi-channel section of OCTEM response on main survey line of low resistance model

图5. 低阻模型主测线OCTEM响应多测道剖面图

Figure 6. Multi-channel section of OCTEM response on main survey line of high resistance model

图6. 高阻模型主测线OCTEM响应多测道剖面图

Figure 7. OCTEM electromagnetic response time slice of low resistance model. (a) 0.162 ms; (b) 0.206 ms; (c) 0.335 ms

图7. 低阻模型物理模拟OCTEM 电磁响应时间切片。(a) 0.162 ms；(b) 0.206 ms；(c) 0.335 ms

Figure 8. OCTEM electromagnetic response time slice of high resistance model. (a) 0.263 ms; (b) 0.885 ms; (c) 8.576 ms

图8. 高阻模型物理OCTEM电磁响应时间切片。(a) 0.263 ms；(b) 0.885 ms；(c) 8.576 ms

5. 结论

1) 基于矢量有限单元法，采用六面体网格对计算区域进行剖分，结合等值反磁通理论对反向对偶磁源进行分析，应用伽辽金法推导有限元方程，通过傅里叶变换将电磁响应从频率域变换到时间域，实现了空洞的时间域OCTEM电磁响应数值模拟计算，并通过与物理模拟结果对比验证了正确性。

2) 对比数值模拟和物理实验发现，采用矢量有限元法数值模拟计算的城市地下空洞OCTEM响应异常与物理模拟响应趋势一致，而且异常幅值明显，边界清晰，无论是充水–泥的岩溶和纯空的空洞，OCTEM都有较好的识别能力，因此，OCTEM应用于城市地下空洞探测是可行的。