地面回线源瞬变电磁法一维反演系统及其应用

张文波，张莹，李建慧

(中国地质大学(武汉) 地球物理与空间信息学院，湖北 武汉 430074)

0 引言

瞬变电磁法在地下水调查、金属矿产勘探、煤矿防治水等领域发挥着重要作用。瞬变电磁法三维反演是当前的研究热点，但实用化的资料处理解释手段仍以一维反演为主[1-2]。

电磁法一维反演的常用方法有最小构造反演、Occam反演、横向约束反演(laterally constrained inversion，LCI)、空间约束反演(spatially constrained inversion，SCI)，以及模拟退火法、深度学习等多种人工智能优化算法。一维最小构造反演是求取多层介质模型的最光滑解，即在一定拟合差条件下，使得模型粗糙度最小，因此，目标函数不仅包含数据拟合项，也包含模型正则化项[3-4]。Occam反演与最小构造反演类似，不同之处在于正则化因子的选取：Occam反演通过一元函数优化来选取正则化因子，而最小构造反演中正则化因子通常固定不变或者通过一些非常简单的算法选取[5]，因此，Occam反演更加稳定，对正则化初始值依赖程度低，但其计算量偏大。LCI方法同时反演电阻率和层厚度，并在目标函数中加入相邻测点之间电阻率、层厚度和界面深度约束项[6]，因此单条测线所有测点数据需整体反演。SCI与LCI类似，不仅考虑同一条测线相邻测点之间的约束，还需考虑其他测线相邻测点之间的约束[7]，因此测区内所有测点数据需整体反演。

目前，国内外多个课题组均已开发了瞬变电磁法一维反演软件。如：加拿大英属哥伦比亚大学Oldenburg课题组在国际上率先实现了用于回线源瞬变电磁法数据处理的一维最小构造反演[8]，开发了实用化的软件包EM1DTM，并成功应用于多个野外实例[9-10]；德国科隆大学Tezkan课题组采用阻尼最小二乘法和Occam算法开展了回线源瞬变电磁法、长偏移距瞬变电磁法、射频大地电磁法等方法的单一反演或2种以上方法联合反演[11-14]，并开发了EMUPLUS软件包；丹麦奥胡斯大学Auken课题组基于阻尼最小二乘法开发了瞬变电磁法一维LCI程序，用于解决反演电阻率和层界面横向连续性差等问题[15]，还开发了SCI程序，用于解决反演电阻率和层界面水平方向连续性差的问题[7]；Kirkegaard和Auken在并行计算、迭代法求解方程组、灵敏度矩阵近似等方面优化了算法，进而提升了反演效率[16]。上述反演算法和策略均集成于AarhusInv软件之中。此外，IX1D、Maxwell等商业软件也都包含瞬变电磁法一维反演模块。

国内多个课题组也开展了相关研究。吉林大学殷长春课题组先后开发了适用于航空瞬变电磁法资料处理解释的一维LCI和SCI程序[17-18]，采用Occam算法对m序列发射波形多道瞬变电磁法数据开展了反演研究[19]，还采用深度学习方法实现了航空瞬变电磁法的一维反演快速成像[20]；北京大学黄清华课题组开展了地面回线源瞬变电磁法的一维最小构造反演研究，该算法适用于水平方向发生形变的复杂发射回线[21]；中国矿业大学(北京)程久龙课题组开展了粒子群优化——阻尼最小二乘混合算法研究，并应用于矿井瞬变电磁法探测[22]；山东大学孙怀凤等采用模拟退火法开展了回线源瞬变电磁法一维反演，并取得了良好效果[23]。

综上所述，瞬变电磁法一维反演技术已趋于成熟，并在资料处理中发挥着重要作用。近年来，笔者团队初步开发了一套地面回线源瞬变电磁法一维反演系统。该一维反演系统基于Fortran语言编写，并融合了以下2个开源程序：Occam反演程序[24]与电磁场从频率域向时间域转换的正弦和余弦变换程序[25]。为了提高计算效率，我们对汉克尔变换、矩阵与矢量乘积、一般矩阵乘积等环节开展了基于OpenMP的并行计算。

1 反演算法原理

上述地面回线源瞬变电磁法一维反演系统包括基于高斯牛顿法的最小构造反演和Occam反演，基于阻尼最小二乘法的横向约束反演和空间约束反演(图1)。其中，一维正演基于自研算法开发[26-27]，灵敏度矩阵采用扰动法计算。

图1 一维反演框架Fig.1 The framework of the 1D inversion system

1.1 最小构造反演与Occam反演

基于Tikhonov正则化方式构建反演目标函数φ，有：

φ=‖Wd(dobs-dpre)‖2+λ‖Wm(m-m*)‖2。

(1)

式中:右边第一项为实测数据与正演模型响应的拟合差函数，第二项为模型正则化函数，‖·‖表示l2范数；dobs为实测数据向量；dpre为模型参数向量m相应的正演模型响应向量，由瞬变电磁法一维正演计算；Wd为数据加权矩阵，此处为一个对角矩阵，其元素为实测数据标准差的倒数；Wm为模型粗糙度算子；λ为正则化因子；m*为参考模型，先验信息包含其中。这两种反演算法中，m为各层介质电导率。

对于目标函数式(1)，采用高斯牛顿法[28]最优化时，第k+1次迭代的正规方程为：

(2)

式中：Jk为正演函数在模型参数mk处的一阶导数，即雅克比矩阵；δmk+1为模型更新量。最小构造反演中，通常采用“冷却法”更新正则化因子[29]。算法流程见图2。

图2 一维最小构造反演流程Fig.2 The flow chart for the 1D minimum-structure inversion

Occam反演中，上述目标函数和正规方程可保持不变，其与最小构造反演的主要区别在于如何更新正则化因子，具体策略见文献[5]和[28]。

1.2 LCI和SCI

LCI和SCI中，构建如下目标函数φ：

φ=‖Wd(dobs-dpre)‖2+α‖Rp(m-m*)‖2+

β‖Rh(m-m*)‖2。

(3)

式中:右边第一项为实测数据与正演模型响应的拟合差函数；第二项为模型电导率和层厚约束项，Rp和α分别为其约束矩阵和权重；第三项为模型层界面深度约束项，Rh和β分别为其约束矩阵和权重。其中，Rp和Rh具体形式可参考文献[6],其他参数与上文一致。

对于目标函数式(3)，采用阻尼最小二乘法[28]最优化时，第k+1次迭代的正规方程为：

(4)

式中：I为单位矩阵，λ为阻尼因子，仍采用“冷却法”更新。图3为LCI和SCI流程图，合成正规方程后，其流程与最小构造反演一致。

图3 一维LCI和SCI流程Fig.3 The flow chart for the 1D LCI and SCI

LCI和SCI的区别在于参与约束的测点不同。如图4所示,LCI中对于某一测点(红色)，参与约束的为该测点所在测线中的2个相邻点(绿色)，而SCI中对于某一测点(红色)，参与约束的为该测点的所有邻近点(绿色)。本研究中，以某一测点为圆心，按特定搜索半径自动寻找各个方向与该测点距离最近的测点，并将这些测点作为约束测点。

图4 LCI和SCI参与约束测点示意Fig.4 The diagram showing the observation points used as constrained points in the LCI and SCI

2 应用实例

2.1 实例概述

内蒙古自治区那仁宝力格煤田勘探区内地势平坦，新近系火山喷出岩比较发育，火山岩浆喷溢时的通道穿越煤层时不但对煤层产生破坏，而且对煤的变质也会有不同程度的影响。为了保证今后煤田的顺利开采，瞬变电磁法被用于确定岩浆上涌通道，以及通道周边玄武岩与围岩的界面。

从火山口和火山锥的存在和分布形态看，岩浆以中心式喷发为主，裂隙式喷发次之；岩性组合主要有橄榄拉斑玄武岩和玄武岩。勘探区内，玄武岩体电阻率可达上千欧姆米，而砂岩和泥岩地层电阻率仅为几十欧姆米，二者电性差异明显[30]。这也是采用瞬变电磁法圈定玄武岩岩体在沉积岩中分布范围的前提。

如图5所示，野外工作采用了大回线源装置，发射回线边长为600m×300m，其中600m边长沿测线方向布置。本次野外工作布设了11条测线，线距为100m，点距为40m，共计486个测点。160测线为控制测线，穿越了整个玄武岩露头，其他测线均位于玄武岩露头中心区域。对于160测线，铺设1次回线源只采集6个测点数据；对于其他测线，铺设1次回线源可采集2条相邻测线的12个测点数据。采用加拿大Phoenix公司生产的V8多功能电法仪，观测时间序列为分布于0.5～20ms之间的30个时刻。

图5 那仁宝力格煤田测点分布示意Fig.5 The distribution diagram of measuring points in the Narenbaolige coalfield

2.2 单点反演

为了验证算法的正确性，首先将本研究开发的最小构造反演和Occam反演程序与商业软件IX1D对比。对于最小构造反演和Occam反演，初始模型为36层层状介质，每层电阻率均为101.5Ω·m，目标拟合差为3。IX1D软件使用上述初始模型时，在多个测点无法拟合实测数据，因此使用了软件自动估算的初始模型。由于IX1D软件无法设置目标拟合差，其在多数测点的最终均方根拟合差(RMS)大于3(图6所示)；而本研究开发的最小构造反演收敛性良好，所有测点均方根拟合差均在3以内。

图6 IX1D软件和本研究最小构造反演结果的RMS对比Fig.6 The comparison of RMS for IX1D software and the minimum-structure inversion method presented here

图7为160线由单点反演地电模型拼接而成的电阻率二维断面。图中，高阻异常体(玄武岩岩体)形态对称，呈“锅底状”分布，并且3种方法的反演电阻率分布范围基本一致。其中，IX1D反演结果在相邻测点之间跳跃性较强，玄武岩岩体与沉积岩界面连续性较差。这一现象与部分测点最终均方根拟合差偏大有关。此外，3种反演结果与钻孔揭露的玄武岩厚度吻合良好，仅在钻孔ZK3处差异明显。该钻孔终孔深度为716.6m，未揭露下伏沉积岩，结合该钻孔位于高阻异常体中心，推断其极有可能位于火山通道分布范围之内。因此，一维反演能够获得玄武岩岩体除岩浆上涌通道区域外的分布形态，却无法直接反映岩浆上涌通道信息。这一现象与电磁场的传播规律有关，即随着探测深度的增大，电磁成像分辨率随之下降，也与玄武岩岩体“倒锥形”的真实形态有关。

图8为测区最小构造反演结果俯视图，色标与图7保持一致。图中，测线中心区域电阻率最大，色标为桔红色至红色，测线两端电阻率偏低，显示为绿色甚至蓝色。颜色的这种分布规律一定程度上反映了玄武岩厚度。另外，尽管多数测点与其四周相邻测点的电阻率在同一个色标范围内，但它们之间电阻率差异仍然较大。

图7 160线单点反演电阻率断面Fig.7 The section view stitched from the single-point 1D inversion models for survey line 160

图8 测区最小构造反演电阻率平面(z=0 m)Fig.8 The plane view of the 1D minimum-structure inversion results (z=0 m)

2.3 LCI反演

根据最小构造反演和Occam反演结果，研究区域电性结构相对简单，界面连续性良好，适合开展LCI反演。反演中，初始模型设置为三层层状介质，每层电阻率均为101.5Ω·m，第一、第二层的厚度均为100m；约束权重α和β依次设置为10、100和1000；目标拟合差仍为3。在移动工作站(处理器为4核8线程、2.7GHz主频)、Ubuntu操作系统条件下，单条测线一次反演耗时约为10min，最终拟合差均小于3。

图9为160线LCI反演电阻率断面， 其中玄武岩岩体(红、黄和绿色区域)形态与上一节中单点反演结果基本一致，依然无法直接反映岩浆上涌通道信息。图9中，随着约束权重的增大，相邻测点之间的玄武岩岩体电阻率差异缩小，以及玄武岩与沉积岩界面深度连续性不断增强，并且反演获取的玄武岩岩体厚度也略有增大，但在钻孔ZK3处仍未达到400m。这种横向连续性的增强是LCI反演的特点，也是与单点反演结果的主要区别。

图9 160线LCI反演电阻率断面Fig.9 The section view stitched from LCI inversion models for survey line 160

图10为测区LCI反演结果俯视图(α=β=1000),图中同一测线相邻测点之间的电阻率和层厚连续性良好，但测线之间相邻测点的电阻率差异依然较大。

图10 测区LCI反演结果俯视图(z=0 m，α=β=1 000)Fig.10 The plane view of the LCI results (z=0 m, α=β=1 000)

2.4 SCI反演

SCI反演同样适用于该研究区域。反演中，初始模型设置与上一节LCI反演相同；约束权重α和β依次设置为10、100和1000；目标拟合差仍为3。整个测区一次反演耗时约为60min。图11为160线SCI反演电阻率断面，其与LCI反演结果(图9)非常相似，同一条测线相邻测点之间电阻率和层厚连续性随约束权重α和β变化的规律也一致。

图11 160线SCI反演电阻率断面Fig.11 The section view stitched from SCI inversion models for survey line 160

图12～图14为不同约束权重下测区SCI反演结果俯视图。如图所示,随着约束权重的增大，任一测点与四周相邻测点的电阻率连续性不断增强，目标函数中空间约束项作用逐渐凸显。

图12 测区SCI反演结果俯视图(z=0 m，α=β=10)Fig.12 The plane view of the SCI results (z=0 m, α=β=10)

图13 测区SCI反演结果俯视图(z=0 m，α=β=100)Fig.13 The plane view of the SCI results (z=0 m,α=β=100)

图14 测区SCI反演结果俯视图(z=0 m，α=β=1 000)Fig.14 The plane view of the SCI results (z=0 m, α=β=1 000)

3 结论

本文详细叙述了笔者团队研发的地面回线源瞬变电磁法一维反演系统，并将其应用于内蒙古那仁宝力格煤田玄武岩岩体形态探测实例。该实例中，单点反演、LCI和SCI刻画了玄武岩岩体除岩浆上涌通道区域之外的分布形态。下一步，将以上述一维反演结果为初始模型，开展该实例的三维反演研究，力争获得更加精确的玄武岩岩体形态。