地–巷可控源电磁三维反演算法研究

3-D Inversion Algorithm for Surface-Tunnel Controlled-Source Electromagnetic Data

摘要: 可控源电磁法是一种重要的地球物理勘探方法，也是金属矿等战略性矿产资源的勘查的重要方法之一。常规的可控源电磁法通常采用地面发射–地面接收的观测装置，在一定程度上限制了其勘探深度。若使用地面发射–巷道接收的观测装置，在离地下目标体更近的巷道采集数据，有望提升对深部结构的探测能力。然而，当前地–巷可控源电磁法鲜有研究，为研究地–巷可控源电磁法的可行性和勘探效果，文章首先开发了地–巷可控源电磁三维反演算法，然后设计了理论地–巷勘探模型，对该模型进行了反演测试，结果表明，联合反演地面–巷道数据可有效提升可控源电磁法的深部探测分辨率，为地–巷可控源电磁法的实际应用提供了理论参考和依据。

Abstract: Controlled-source electromagnetic (CSEM) method is a crucial geophysical exploration technique and plays a significant role in the exploration of strategic mineral resources such as metal ores. Conventional CSEM typically employs a surface transmitter-surface receiver configuration, which, to some extent, limits its exploration depth. By adopting a surface transmitter-underground tunnel receiver configuration, where data are collected in tunnels closer to underground targets, the detection capability for deep structures can potentially be enhanced. However, research on surface-tunnel CSEM remains scarce. To investigate the feasibility and exploration effectiveness of surface-tunnel CSEM, this paper first develops a 3D inversion algorithm for surface-tunnel CSEM. Subsequently, a theoretical surface-tunnel exploration model is designed, and inversion tests are conducted on this model. The results demonstrate that joint inversion of surface-tunnel data can effectively improve the resolution of CSEM for deep exploration, providing theoretical references and a foundation for the practical application of surface-tunnel CSEM.

1. 引言

可控源电磁法具有观测范围广、工作效率高、勘探深度大等特点，被广泛应用于矿产资源勘探、工程勘查、地下水勘查等领域[1]-[3]，取得了良好的应用效果，尤其是针对金属矿等战略性矿产资源的勘探，效果显著。然而，随着当前浅部矿产资源的逐渐枯竭，开展深部矿产资源勘查已经是必须的战略选择，这也对可控源电磁法的探测深度和精度提出了更高的要求。老矿山深边部找矿是当前深部矿产资源探测的重要方向，而开发生产多年的老矿山通常建有大量的地下采矿巷道，若能在巷道内布置测点进行数据采集，则有望获得更强的深部异常感应信号，从而提升深部探测能力和效果。然而，当前地–巷可控源电磁法的研究相对较少，可控源电磁法反演研究也主要是针对地面数据。

反演成像是可控源电磁法数据解释的基础，反演问题实际上是针对目标函数的最优化问题。对于可控源电磁问题，由于观测数据与地下模型之间的关系是非线性的，其反演问题是一种非线性最优化问题。对于非线性最优化问题，使用无需任何近似的全局类的最优化算法是最为合适的，但全局类的方法通常需要大量的正演计算，对于大规模三维反演，其时间成本是难以接受的，因此主要被应用在一维、二维反演中[4]-[7]。当前，电磁法三维反演中使用的主流最优化算法仍然是局部搜索的梯度类算法，尽管梯度类的方法不能保证搜索到全局最优解，但相对于全局类的方法在处理三维问题时在效率上具有显著的优势，在电磁法反演中得到广泛应用，主要包括高斯牛顿法(GN) [8] [9]，非线性共轭梯度法(NLCG) [10]-[12]，和有限内存拟牛顿法(L-BFGS) [13] [14]。三种方法各有优劣，本文选择无需显式计算和存储灵敏度矩阵的L-BFGS方法作为反演最优化算法。

传统的反演方法大多基于结构化的六面体网格，难以实现局部加密。地下巷道通常是狭长的结构，截面通常只有数米，需要非常小的网格单元进行剖分，而巷道延伸可达几十米至上百米，使用结构化的六面体网格对巷道进行剖分时，会在巷道以外的区域产生大量的冗余的网格，增加计算成本。为此，本文引入八叉树网格，开发了基于八叉树网格的地–巷可控源电磁三维反演算法，八叉树网格具备良好的局部加密能力，可以有效地离散巷道模型，并且节省网格数量，从而节省计算资源消耗。

2. 反演方法

2.1. 目标函数

反演采用如下正则化的目标函数[15]

$\varphi \left(m\right)={\varphi }_d\left(m\right)+\lambda {\varphi }_m\left(m\right)$ (1)

其中$m$ 为模型参数向量，${\varphi }_d\left(m\right)$ 为数据拟合项，${\varphi }_m\left(m\right)$ 为模型约束项，$\lambda$ 为正则化因子，用于平衡数据拟合与模型约束之间的权重。本文在反演算法中采用冷却法更新正则化因子，即当前后两次迭代的拟合差小于预先设定的阈值时，就将正则化因子减小到原来的十分之一。数据拟合项的具体表达式为

${\varphi }_d\left(m\right)={\Vert W_d\left(d_{pre}-d_{obs}\right)\Vert }_2^2$ (2)

其中$W_d$ 为数据协方差矩阵，是由数据误差的倒数组成的对角阵，$d_{pre}$ 为预测数据向量，$d_{obs}$ 为观测数据向量。模型约束项的表达式为

${\varphi }_m\left(m\right)={\Vert W_m\left(m-m_{\text{ref }}\right)\Vert }_2^2$ (3)

其中$W_m$ 为模型约束矩阵，$m_{\text{ref }}$ 为参考模型向量。

反演中采用均方根误差(RMS)定义数据拟合差，表达式为

$\text{RMS}=\sqrt{\frac{{\varphi }_d\left(m\right)}{N_d}}$ (4)

其中$N_d$ 为观测数据个数。理想情况下，RMS下降到1时，可认为反演收敛，此时可认为预测数据在噪声水平下拟合了观测数据。

2.2. L-BFGS最优化算法

本文采用L-BFGS算法求解目标函数最小化问题，在第k次迭代，L-BFGS算法的模型更新公式为

$m_{k+1}=m_k+{\alpha }_k{p}_k$ (5)

其中${\alpha }_k$ 是步长，本文通过基于Wolfe条件的非精确线搜索得到，$p_k$ 为搜索方向，其表达式为

$p_k=-H_k^{-1}{g}_k$ (6)

其中$g_k$ 为目标函数梯度向量，$H_k$ 为目标函数二阶导数矩阵，即海森矩阵，表达式分别为

$g_k=2\mathrm{Re}\left\{J_k^{\text{T}}{W}_d^{\text{T}}{W}_d{\left[f\left(m_k\right)-d_{obs}\right]}^{*}\right\}+2\lambda W_m^{\text{T}}{W}_m\left(m_k-m_{ref}\right)$ (7)

$H_k=2\mathrm{Re}\left\{J_k^{\text{T}}{W}_d^{\text{T}}{W}_d{J}_k\right\}+2\lambda W_m^{\text{T}}{W}_m$ (8)

本文采用采用谱元法结合直接求解器求解反演中的正演问题[16]，然后使用伴随正演的方法高效求解式(7)中所表示的目标函数梯度。式(8)所表示的海森矩阵是一个大型密实矩阵，显式计算和存储以及求逆都是非常困难的，L-BFGS算法无需精确计算海森矩阵，而是通过构建一个海森矩阵逆矩阵的近似矩阵$B_k$ 来获得搜索方向，则式(6)可重写为

$p_k=-B_k{g}_k$ (9)

为了减小内存消耗，L-BFGS算法通常只利用最近几次迭代的信息构建$B_k$ ，其表达式为

$B_k=\left(v_{k-1}^{\text{T}}\cdots v_{k-m}^{\text{T}}\right)B_k^0\left(v_{k-m}\cdots v_{k-1}\right)+\left(v_{k-1}^{\text{T}}\cdots v_{k-m+1}^{\text{T}}\right){\rho }_{k-m}{s}_{k-m}{s}_{k-m}^{\text{T}}\left(v_{k-m+1}\cdots v_{k-1}\right)+\cdots +{\rho }_{k-1}{s}_{k-1}{s}_{k-1}^{\text{T}}$ (10)

其中

${\rho }_k=\frac{1}{y_k^T{s}_k},v_k=I-{\rho }_k{y}_k{s}_k^{\text{T}}$ (11)

$s_k=m_k-m_{k-1},y_k=g_k-g_{k-1}$ (12)

L-BFGS算法不需要通过显式计算式(10)然后代入式(9)的方式得到搜索方向，可以利用基于上述表达式推导得到的双循环递归算法[17]直接计算得到搜索方向，该算法仅涉及到一系列向量运算，其计算速度和内存消耗都是非常小的，非常适用于三维反演。

3. 模型算例

设计如图1所示的地–巷可控源电磁法勘探模型，电导率为0.01 S/m的均匀半空间下含两个异常体，几何尺寸为1 km × 1 km × 0.5 km，浅部异常体的顶部埋深为300 m，电导率为0.1 S/m，深部异常体顶部埋深为1200 m，电导率为0.2 S/m。在两个异常体之间，深度为1 km的位置设置9个巷道，巷道沿x轴方向延伸，长度为1.8 km，横截面尺寸为6 m × 6 m，巷道的间距为200 m，巷道内电导率设为空气电导率10⁻⁸ S/m。在地面和巷道内同时布置测点，点距为100 m，地面共计441个测点，巷道内共计153个测点。场源布置在地面距测区中心5 km处，长度为100 m，电流为1 A。利用正演计算1~8192 Hz的14个频率的Ex分量数据，并加入3%的高斯噪声，作为反演观测数据。反演区域为−1.2 km < x < 1.2 km，−6 km < y < 1.2 km，0 km < z < 3 km，巷道处的网格剖分如图2所示，利用八叉树网格可以实现对巷道区域的局部加密，而不影响周边的网格，避免了增加冗余的网格单元。

对于该模型，本文采用两种反演方案进行对比，第一种方案只反演地面数据，第二种反演方案同时反演地面数据和巷道数据，考虑到高频信号探测深度有限，其中巷道数据仅使用100 Hz 以下的数据。反演中，巷道中的电导率固定(空气电导率)，不随反演过程变化。

Figure 1. Diagram of the of the surface-tunnel model

图1. 地面–巷道模型示意图

Figure 2. Cross-sectional view of octree mesh discretization at the tunnel (Top right corner is a localized zoomed-in 5x display)

图2. 巷道处八叉树网格剖分截面(右上角为局部放大5倍显示)

图3展示了两种反演的数据拟合差曲线，反演中收敛阈值均设为1.001，仅反演地面数据时，初始拟合差为7.05，经过65 次迭代后收敛，联合反演地面数据和巷道数据时，初始拟合差为11.67，经过63 次迭代后收敛，两种反演都较为稳定，收敛性也较为接近。

Figure 3. Convergence curve of inversion iteration for the surface-tunnel model

图3. 地–巷模型反演迭代收敛曲线

Figure 4. Sectional view of inversion results for the surface-tunnel model

图4. 地–巷模型反演结果切片

图4展示了两种反演的结果。左侧一列为仅反演地面数据的结果，右侧一列为联合反演地面和巷道数据的结果，自上而下分别为y = 0 m，x = 0 m，z = 550 m，z = 1300 m处的切片，黑色方框表示真实异常体的位置。分析对比图中结果可以发现，对于浅部异常体，两种反演都能较好的恢复出来，验证了算法的正确性，反演恢复的异常位置基本与真实异常位置重合，形态也接近真实异常体，仅在异常底部边界略有差异，这是因为反演中使用了光滑约束；对于深部异常体，只反演地面数据时几乎没有反映，而联合反演地面和巷道数据则可以恢复出深部异常体，证明巷道数据的加入有效地提升了深部分辨率。需要说明的是，反演恢复出的深部异常体在位置和形态上基本能反映真实异常体的信息，但恢复的精度有限，整体上看，深部异常体的反演效果相较于浅部异常体要差，这是因为场源在地面激发，深部区域无论是一次场还是感应二次场强度都较弱，尽管在巷道内置测点相对于地面测点可以采集到更强的信号，但强度仍然有限，也限制了其分辨率。此外，本文所设计的巷道模型仅为简化的理论模型，实际矿区地下情况更加复杂，并且巷道内通常建有轨道、支护钢架、机电设备等基础设施，也会对电磁信号造成较强的干扰。因此，对于地巷可控源电磁勘探，需要考虑信号强度及强噪声干扰的问题，若能通过合理的技术手段在巷道内采集到一定强度的有效信号，则可以在一定程度提升可控源电磁法的探测能力和深部分辨率。

4. 结论

本文成功开发了地–巷可控源电磁三维反演算法，采用具备局部加密能力的八叉树网格进行模型离散，利用无需显式计算和存储海森矩阵的L-BFGS算法求解反演目标最小化问题。设计了地面–巷道可控源勘探模型，对该模型进行了反演，并对比了仅反演地面数据和联合反演地面数据及巷道数据两种采集方案的反演效果，反演结果表明，相对于传统的地面可控源电磁法，地–巷可控源电磁法可以在一定程度上改善深部探测的分辨率，具备一定的实际应用前景。本文仅从理论层面进行了初步研究，未来将考虑建立更加复杂、更加接近真实地质情况的模型进行深入研究。

基金项目

地球深部探测与矿产资源勘查国家科技重大专项(2024ZD1002202)。

NOTES

^*通讯作者。