地面与半航空瞬变电磁法三维联合反演

刘明宏， 蔡红柱,2*， 杨浩， 熊咏春, 胡祥云,2

1 中国地质大学(武汉)地球物理与空间信息学院， 武汉 430074 2 地质过程与矿产资源国家重点实验室， 武汉 430074

0 引言

瞬变电磁法(Transient electromagnetic method, TEM)是一种电磁感应类地球物理探测方法，通过观测激励源断电后的二次场响应，对观测信号进行成像或反演分析可获得地下空间的电性分布信息.目前，瞬变电磁法已经被广泛应用于矿产勘查(武欣等, 2016; Yang and Oldenburg, 2012)、油气勘探(Li and Constable, 2010)和地下水调查(Auken et al., 2017; Trabelsi et al., 2013)等领域.

瞬变电磁法通常采用磁性回线源或电性接地导线源作为激励，可在地面或空中等平台进行观测.传统的地面中心回线源瞬变电磁观测装置具有探测深度大、分辨率高、抗干扰能力强等特点，应用非常广泛，但在山地、丛林、沼泽、冰川等艰险环境，地面工作难度极大.航空电磁法是探测地面环境复杂区域的有效手段，然而航空电磁测量对仪器设备技术要求高，成本昂贵、飞行风险大.将发射源设置在地面，使用无人机在空中观测的SATEM系统成为了一种新兴的地球物理探测技术(Meng et al., 2020; Wu et al., 2019).SATEM系统相对于地面系统，采用空中接收工作方式，具有效率高、不受地表环境限制的优点；相对于航空系统，采用地面发射方式，具有发射功率大、数据信噪比高、勘探深度大的优势(Ito et al., 2014; Smith et al., 2001).图1展示了地面与半航空瞬变电磁装置系统的示意图.地面中心线源装置通常在回线框内的一定范围内采集数据，半航空瞬变电磁一般采用接地长导线作为发射源，在偏离源的区域飞行采集数据.

图1 地面瞬变电磁与半航空瞬变电磁装置示意图Fig.1 Schematic diagram of the ground-based TEM (GTEM) and semi-airborne TEM

半航空电磁方法的开发旨在于增加航空电磁探测的深度.Nabighian(1988)提出了SATEM的概念.Elliott(1996, 1998)开发了基于回线源的FLAIRTEM系统，并成功应用于巴布亚新几内亚的矿床调查.加拿大Fugro公司开发了回线源TerraAir系统(Smith et al., 2001).Mogi等(1998)开发了首套采用接地导线发射源SATEM系统(GREATEM)，成功应用于地热调查、火山构造探测和海水入侵等多个领域，展示了电性激励源SATEM系统具有高信噪比、探测深度大的优势(Abd Allah et al., 2013; Ito et al., 2014; Mogi et al., 2009).我国在SATEM系统的研发与应用上也取得了不错的进展.嵇艳鞠等(2013)和Wang等(2013)研制了一套以无人飞艇为搭载平台的时域电磁接收系统，并开展了海水入侵调查和水资源探测试验.刘富波等(2017)基于无人直升机搭载平台，在山东昌邑莲花山铁矿区开展了电磁探测飞行试验.随着无人机技术的迅速发展，以无人机为搭载平台的SATEM成为了近年的研究热点.方涛等(2015)、Xue等(2018)、吴启龙(2019)和谢小国等(2021)将无人机SATEM系统应用于地下采空区、地下巷道、隧道工程调查、古河道探测等领域.

在解释方面，一些学者采用视电阻率成像的方法进行瞬变电磁数据解释(Christiansen et al., 2006; 马振军等, 2021; 李貅等, 2015; 张莹莹等, 2015; 赵越等, 2015).此外，采用一维反演获得层状电导率模型是瞬变电磁数据解释的主要方法(Farquharson and Oldenburg, 1993; 李永兴等, 2010; 李展辉和黄清华, 2018)，但根据一维反演结果构建的拟二维剖面通常会产生电阻率突变以至于难以进行地质解释.随着横向约束和空间约束方法的引入，瞬变电磁一维反演结果的空间连续性得到了增强(Auken et al., 2008; Viezzoli et al., 2008).然而在复杂的地质情况下，传统的视电阻率成像方法和一维反演很难恢复精确的地下结构.

三维反演是对TEM数据精确解释的有效方法.正演是反演的基础，前人在TEM三维正演上进行了大量研究.积分方程法(Wannamaker et al., 1984; Zhdanov et al., 2006; Zhu et al., 2019; 殷长春和刘斌, 1994)、有限差分法(Commer et al., 2015; Commer and Newman, 2004; Maaø, 2007; Wang and Hohmann, 1993; Yu et al., 2017; 孙怀凤等, 2013; 余翔等, 2017)被应用于电磁三维正演模拟.随着计算硬件的快速发展，尤其是高性能集群的应用，TEM三维反演已成为可能.基于积分方程法，结合时频转换技术和移动足迹(“moving footprint”)方法，Cox等(2012)开发了TEM三维反演方法，Wang等(2021)开发了适用于多道瞬变电磁数据的三维并行反演算法.基于有限差分法和时频转换方法，Liu和Yin(2016)实现了航空TEM三维反演.随后，Su等(2021)将其推广到TEM各向异性三维反演.基于规则网格的有限体积法的TEM三维反演也有少量报道(Oldenburg et al., 2013; Ren et al., 2018; Yang et al., 2014).

在先前报道的基于积分方程、有限差分或有限体积法的TEM三维反演，采用了规则网格对模型进行离散，难以考虑复杂的地质构造(如地形)的影响.近年来，基于非结构化四面体网格的有限元法在地球物理电磁模拟中的应用受到广泛关注(Badea et al., 2001; Cai et al., 2017; Fu et al., 2015; Jahandari et al., 2017; Li et al., 2018; Um et al., 2010, 2012; Yin et al., 2019).Um等(2010)使用基于非结构化四面体网格的有限元方法来计算时域电磁场.此外，他们还提出了一种时间步长加倍的方法减少所需计算的时间步进(Um et al., 2012)，随后引入了并行策略加速计算(Fu et al., 2015).Cai等(2017)基于非结构化四面体网格的时域有限元方法，结合自适应时间步进法和混合边界条件策略，提高了正演的准确性和计算效率.殷长春等(2019)提出了基于自适应有限元方法的时间域三维正演算法，提高计算效率.Jiang等(2019)采用全场控制方程的有限元法进行了地-井时域电磁模拟，并利用煤矿含水层模型进行了验证.Li等(2020)开发了考虑地形影响的接地线源瞬变电磁法有限元三维正演算法，并将其应用于Ovoid地区硫化物矿床的数值模拟.

近年来，基于矢量有限元的三维TEM反演的研究也取得进展.Liu等(2019)基于无约束四面体网格矢量有限元法实现了三维TEM反演，反演采用BFGS优化方法.Qi等(2022)开发了一种结合局部网格技术和解耦正反演网格技术的航空TEM数据三维反演方法.Zhang等(2021)提出了一种三重网格方法的TEM三维反演策略，正演计算使用细的非结构化网格以保证精度，灵敏度计算使用粗的非结构化网格以减少计算，反演采用规则的结构化网格以便于约束反演参数.基于有限元的TEM三维反演多报道于地面装置和航空装置，但少有关于SATEM的三维反演报道.

可靠性和分辨率是地球物理反演最关注的问题.单一的电磁探测装置对地下结构的灵敏度具有局限性，难以获得地质结构的全面信息.不同TEM装置对地下的电性结构的灵敏度不同，联合地面与半航空TEM数据将形成互补优势，比单独观测方式对地下的电性结构具有更高的灵敏度.另一方面，地球物理反演普遍存在多解性问题，多种地球物理方法的数据联合反演是提高反演结果可靠性的有效方法.联合反演已被广泛应用于地球物理数据解释(Gallardo and Meju, 2004; Lelièvre et al., 2012; Zhdanov et al., 2012)，但关于地面TEM与SATEM的三维联合反演未见报道.

为了实现可靠的瞬变电磁三维反演，本文开发了并行的瞬变电磁三维正演与反演算法.正演采用基于非结构化四面体网格的矢量有限元法，可以精确模拟复杂结构.为了提高计算效率，在求解正演和灵敏度计算时，本文采用自适应时间步长的方法减少时间步进，并利用Intel MKL PARDISO并行直接求解器求解大型线性方程组.在反演中，为了提高反演收敛速度，采用具有拟二阶收敛性的高斯-牛顿法反演策略，通过最小二乘QR分解(LSQR)方法获得模型更新方向.为结合不同装置的探测优势以及降低反演的多解性影响，基于构建统一目标函数和数据权重调整的反演策略，实现了地面TEM和SATEM数据的联合反演.最后通过理论模型研究，分析了地面TEM和SATEM的探测能力，并对比单独反演与联合反演结果，验证了联合反演方法的优势.

1 瞬变电磁三维正演

三维正演模拟是认识电磁规律与开展电磁反演的基础.从时域的准静态Maxwell方程组出发，在忽略位移电流和考虑真空磁导率的情况下得到(Ward and Hohmann, 1988)：

(1)

(2)

(3)

本文采用基于非结构化四面体网格的一阶矢量有限元方法求解方程(3)(Cai et al., 2017; Jin, 2014; Liu et al., 2019; Um et al., 2012).在不同的时刻，电场被赋值到单元边缘上，四面体单元内的电场可用沿四面体边界的电场的线性组合表示：

(4)

方程(3)的时间导数项需要使用有限差分近似表示.在时域有限元方法中，后退欧拉法因其无条件稳定而常被采用(Haber, 2014; Jin, 2014; Um et al., 2012).在相同步进次数条件下，二阶后退欧拉法的时间差分比一阶方法具有更高精度.殷长春等(2019)和Liu等(2019)使用了具有相同步长的二阶后退欧拉方法近似时间导数项.为了减少所需模拟的时间步进次数，本文采用了灵活的自适应时间步进方法，在一定的步数后，尝试增加时间步长(Um et al., 2012; Cai et al., 2017).

采用二阶后退欧拉方法，在tn时刻的电场时间导数可以表示为

(5)

通过伽辽金有限元分析获得如下的有限元线性方程组(Jin, 2014)：

AnEn=bn，

(6)

其中

(7)

(8)

式中An∈(Ned×Ned)，是一个大型的稀疏矩阵，bn∈(Ned×1)，K和M是刚度矩阵(Jin, 2014).

求解方程6需要给定适当边界条件和初始条件.本文采用Dirichlet边界条件，假设电场在建模区域的外边界上趋近于零.为了匹配Dirichlet边界条件，将模拟区域拓展至一个很大的边界.在测点和发射源的附近加密网格以确保计算精度.在核心区域以外逐渐增大网格的尺寸，以减少网格数量.对于初始条件的选择，通常可以采用两种方式.当计算全波形响应时采用电场为零的初始条件.当忽略上升电流与持续电流只考虑关断电流的响应时，通常求解直流电场的拉普拉斯方程获得初始电场.本文采用了如图2所示的梯形发射波，其脉宽时间较短，不能忽略上升电流(I)、持续电流的影响，因此需要进行全波形计算，考虑电场为零的初始条件：E0=0.本文通过将离散的发射电流密度加入公式(8)实现激励源的加载，可以模拟任意发射波形.电流关断后启用自适应时间步进，初始时间步长大小取决于发射波形的离散.

图2 梯形波Fig.2 Trapezoidal waveform

确定边界条件和初始条件后，可通过迭代法或直接法求解方程(6)得到某时刻的电场En.由于每个时间步进都必须求解有限元方程组，考虑到时间步长数量很大，采用迭代法并不实用(Jin, 2014; Um et al., 2012).采用直接法求解更具优势，当时间步长相同时，刚度矩阵保持不变，An矩阵分解的结果可以被重复使用(Jin, 2014; Oldenburg et al., 2013).在自适应步进方法中，矩阵An只需改变几次，可以显著减少矩阵An的分解次数.为了提高计算效率，采用 Intel MKL PARDISO并行直接求解器求解方程(6)(Schenk et al., 2001).

通过正向步进的方式逐步求解方程(6)得到电场，测点的数据可以通过以下插值矩阵获得：

d=QE，

(9)

其中，Q是与时间和空间相关的插值矩阵，E为所有时刻的电场.

2 高斯-牛顿法反演

TEM反演是不适定问题.为了克服这个问题以得到具有地球物理意义的解，Tikhonov正则化被引入到反演目标泛函中(Tikhonov and Arsenin, 1977)：

φ(m)=φd(m)+λφm(m)，

(10)

其中，φd是数据拟合项，φm是模型稳定项，λ是用于平衡数据拟合和模型约束正则化参数.φd和φm分别定义为

(11)

(12)

其中，m是模型参数，为了保证反演中模型参数具有物理意义，取对数域的电导率作为模型参数，m=log(σ).Wd为数据加权矩阵，瞬变电磁数据跨越多个数量级，为了平衡不同时间通道的数据权重，采用观测数据的倒数进行构建，并在数据变号处数据权重调整为零，以避免无穷大值，类似的方法已成功应用于瞬变电磁反演(Cox et al., 2012).F(m)表示模型m的正演响应，dobs为观测数据，mref为含先验信息的参考模型.Wm为模型粗糙度矩阵，由反演单元与邻近单元的体积和距离的权重构建(Cai et al., 2021).

对目标函数进行二阶泰勒级数展开，并忽略高阶项可以获得以下法方程：

式中δmn为模型更新量，数据误差rn-1=F(mn-1)-dobs，n为迭代次数.法方程通常采用共轭梯度法(Conjugate gradient，CG)求解，该方法需要求解近似海森矩阵以获得有效的预条件因子.在数值计算上，法方程等价于如下最小二乘形式(Hansen, 1998; Aster et al., 2018)：

(14)

LSQR算法是利用Lanczos子空间投影方法求解最小二乘问题的有效方法(Paige and Saunders, 1982)，相对于PCG算法，LSQR算法具有更高的稳定性(Grayver et al., 2013).求解方程组(14)前需要获得灵敏度矩阵Jd，Jd的准确性影响着反演的稳定性和准确性.为了求取数据的灵敏度矩阵，需要先求取电场对模型参数的导数Je,n，式(6)对m取导数和整理得

(15)

第n步的数据对模型参数的灵敏度(Jacobian)可以表示为

(16)

式中，Nd和Ne分别表示数据的数量和模型空间四面体单元数量.通常不需要求解整个模型空间的灵敏度矩阵，只考虑反演区域的灵敏度矩阵，Jd,n的大小变为Nd×Nm，Nm为反演区域的模型参数的数量.数据的数量Nd通常远少于模型参数的数量Nm，因此本文求取灵敏度矩阵的转置.因为系数矩阵An是对称矩阵，则有

(17)

获得模型更新方向δm后，通过线搜索策略找到最优的模型更新步长l，模型更新可以表示为(Nocedal and Wright, 2006)：

mn=mn-1+lδm.

(18)

正则化参数在权衡数据拟合项和模型约束项中起着重要作用，关系着反演过程的稳定.在初始反演迭代中，本文采用近似谱分析来获得正则化参数(Long et al., 2020)，在后续迭代中使用冷却方法减小正则化参数.此外，观测数据与预测数据之间的归一化拟合差可用以下公式表示(Zhdanov, 2009)：

(19)

当归一化拟合差达到设定阈值或收敛停止时，反演终止迭代.

3 联合反演

联合反演可以结合不同探测方法的优势以降低反演的多解性.对于不同物性的联合反演方法，通常需要构建不同物性之间的相关性约束.在地球物理反演中构建相关性约束的方法主要有：交叉梯度(Gallardo and Meju, 2004)、Gramian约束(Zhdanov et al., 2012)、模糊c均值聚类法(Lelièvre et al., 2012)等.地面TEM与SATEM是对同一物性(电导率)进行反演，因此，在联合反演策略中，直接采用与单独反演相同的目标函数，而不需要构建相关性约束关系(Commer and Newman, 2009).

在联合反演中，观测数据dobs为地面TEM的数据和SATEM的数据组合成的列向量：

(20)

式中，dGTEM和dSATEM分别表示地面TEM的数据子集和SATEM的数据子集.

与单独反演不同，在联合反演中数据的权重需要进行适当调整.地面TEM与SATEM的数据采集密度不同，在相同测量面积上，半航空的数据量可能远大于地面数据.如果简单采用与单独反演相同的数据加权矩阵，当地面TEM数据量远小于SATEM数据时，地面数据在反演中难以发挥作用.我们希望在联合反演中，不同数据子集能发挥相应的作用，以结合两种装置的探测能力.因此，本文采用改变不同数据子集的数据加权矩阵的方式，来调整地面TEM数据和SATEM数据的权重(Commer and Newman, 2009).

本文以相同观测区域上不同地球物理数据的数量为依据确定权重系数α1和α2.假设子集dGTEM和dSATEM的数据量分别为N1和N2，则

α1×N1=α2×N2，

(21)

当指定α1=1时，α2=N1/N2.相应的Wd则为

(22)

上式中WGTEM、WSATEM分别表示地面TEM和SATEM的数据加权矩阵.在联合反演中，只需通过改变α2，就可以调节不同数据子集的权重，确保不同数据子集权重的平衡.联合反演中的灵敏度计算、法方程求解等流程与单独反演完全相同.

4 理论模型研究

为了验证所提出的反演算法的有效性，本文设计了两个理论模型.在数值实验中，地面TEM与SATEM均采用如图2所示的梯形发射波.观测时间取等对数间隔10-4.5～10-1.5s，共25个时间道.在所有的模型案例中，使用添加2%的高斯噪声的垂直分量dHz/dt合成数据进行反演.

4.1 模型1：无地形模型

模型1电性结构和观测装置如图3所示.导线源长1500 m，第1号源在X=-1500 m处，第2号源在X=1500 m处.设置6个大小均为500 m×500 m的大回线源，仅在回线圈内观测，采用两种观测点距，50 m点距的数据量为6×81×25=12150个，200 m点距的数据量为6×9×25=1350个.SATEM测区范围在X:-700～700 m，Y:-450～450 m之间，点距为50 m，数据量为551×25=13775个.浅部两个异常体大小为300 m×600 m×60 m，中心点深度为60 m，电导率分别为0.1 S·m-1和0.001 S·m-1；深部异常体大小为400 m×600 m×200 m，中心点深度为350 m，电导率为0.2 S·m-1；背景电导率为0.01 S·m-1.

图3 模型1的切片图(a) Z=60 m切片，红色线段、白色线框分别为SATEM发射源、地面TEM回线框发射源，黑色、白色测点的距离分别为50 m、200 m； (b) Y =0 m切片.Fig.3 Slice view of the model 1(a) Slice at Z=60 m, red line and white loops represent the SATEM electrical sources and the ground TEM loop sources respectively, the spacing of the black and white dots are 50m and 200 m, respectively; (b) Slice at Y = 0.

在环境复杂地区难以进行密集数据采集.为了测试地面TEM的中心回线源装置的探测能力与不同测点数量对地面瞬变电磁三维反演结果的影响，对不同点距的数据进行反演.正演与反演使用不同的网格，反演网格取消对异常体处加密，在模型1的所有反演中均使用同一套网格，整个模型区域大小为20 km×20 km×20 km.反演区域大小为3 km×2 km×1 km，反演区域四面体单元数量为393940，整个模型的四面体单元数量为1037619，反演初始模型电导率为0.01 S·m-1的均匀半空间.

本文的数值实验在中国地质大学(武汉)高性能计算集群上完成，每个节点包含2个Intel(R) Xeon(R) 2.5 GHz的CPU，每个CPU包含20个核，单节点内存为384GB.模型1的反演使用1个节点计算，50 m点距和200 m点距的每次反演迭代分别需要342 min和346 min，测点数量对反演时间影响较小.反演收敛情况如图7a所示，最终反演结果的数据归一化拟合差在2%左右.

图4为地面TEM单独反演的结果.浅层两个薄板异常体的位置和电导率恢复较好，但深部异常体的反演结果类似一个“环状”结构，只能恢复异常体的边界，在异常体中心形成了一个空洞.随着测点数量的减少，反演结果分辨率降低，假异常增加.对于实际探测，难以进行合理的地质解释.

图4 模型1的地面TEM反演结果(a,b) 50 m点距数据反演结果; (c,d) 200 m点距数据反演结果； (a,c) 三维视图； (b,d) Z=350 m切片.Fig.4 The ground TEM inversion results for model 1(a,b) Inversion results with station spacing of 50 m; (c,d) Inversion results with station spacing of 200 m; (a,c) Three-dimensional view; (b,d) Slice view at Z=350 m.

为了测试SATEM装置的探测能力，对两个不同发射源的模拟数据分别进行反演.网格、初始模型和其他参数与上相同，每次反演迭代约4 h.反演收敛情况如图7b所示，最终反演结果的数据归一化拟合差在5%左右.因为在较小的偏移距内进行观测，衰减曲线会出现“变号”，数据拟合差稍大于噪声水平.

SATEM单独反演结果如图5所示.两个反演结果均恢复了大致的地下电性分布情况.但与地面TEM装置反演结果相比，SATEM对异常体边界的分辨能力较差.SATEM单独反演的结果也存在明显的假异常，且在不同发射源的情况下，反演结果的假异常不同，这给数据解释带来巨大挑战.

图5 模型1的SATEM反演结果(a,b) 第1号源数据反演结果； (c,d) 第2号源数据反演结果； (a,c) 三维视图； (b,d) Z=350切片.Fig.5 SATEM inversion results for model 1(a,b) Inversion results with the data excited by the first source; (c,d) Inversion results with the data excited by the first source; (a,c) Three-dimensional view; (b,d) Slice view at Z=350 m.

单一观测方式探测能力的局限性和反演的多解性严重影响反演结果的准确性.联合多种测量装置进行反演是降低多解性的有效手段.我们分别将不同采集密度的地面TEM数据与SATEM第2号源的数据进行联合反演.联合反演采用与单独反演相同的网格、初始模型和反演参数.联合反演耗时与地面TEM单独反演相当，每次反演迭代约341 min.反演收敛情况如图7c所示，地面TEM点距50 m、200 m与SATEM的联合反演最终结果的数据归一化拟合差分别为5.6%和3.9%.

图6为SATEM与地面TEM联合反演结果.联合反演对三个异常体的边界和电导率均得到有效重建.地面TEM单独反演的深部异常体的“环状”结构以及SATEM单独反演的假异常，在联合反演中均得到有效压制.尤其采用50 m点距的地面TEM数据时，联合反演结果基本上没有明显假异常.当地面TEM数据为200 m点距时，虽然数据量只有50 m点距的1/9，但联合反演结果也明显好于单独反演的结果.

图6 模型1的地面TEM与SATEM联合反演结果(a,b) 地面TEM点距50 m数据与SATEM数据联合反演； (c,d) 地面TEM点距200 m数据与SATEM数据联合反演； (a,c) 三维视图； (b,d) Z=350切片.Fig.6 Results of ground TEM and SATEM joint inversion for model 1(a,b) Joint inversion between the ground TEM data with 50 m station spacing and the SATEM data; (c,d) Joint inversion between the ground TEM data with 200 m station spacing and the SATEM data; (a,c) Three-dimensional view; (b,d) Slice view at Z=350 m.

图7 模型1反演的收敛曲线(a) 地面TEM反演； (b) SATEM反演； (c) 地面TEM与SATEM联合反演.Fig.7 Convergence plot for different inversions of model 1(a) The ground TEM inversion; (b) SATEM inversion; (c) Joint inversion between the ground TEM data and the SATEM data.

4.2 模型2：带地形模型

为了更接近真实地质情况，建立了一个带有地形和多个异常体的模型.模型2的地形、电性结构和观测装置如图8所示.该模型地形起伏高差约100 m，X轴负方向较为平坦，X轴正方向地形起伏，分别模拟喀斯特地貌的平原与峰丛.导线源长1000 m，在X=-1500 m处.设置8个大小为500 m×500 m的大回线源，仅在回线内观测，点距50 m，观测数据量为8×81×25=16200个.SATEM测区范围在X:-950～950 m，Y:-450～450 m之间，点距50 m，数据量为741×25=18525个.异常体1的大小为300 m×600 m×100 m，中心点埋深约100 m，电导率为0.001 S·m-1；异常体2的大小为600 m×600 m×100 m，倾斜放置，电导率为0.2 S·m-1；异常体3和异常体4的大小均为600 m×200 m×100 m，中心点埋深约100～200 m，电导率分别为0.2 S·m-1和0.001 S·m-1；背景电导率为0.01 S·m-1.

反演网格取消对异常体处加密，模型2的所有反演中使用同一套网格.反演区域X方向范围：-2000～2000 m，Y方向范围：-1000～1000 m，Z方向厚度约1000 m.反演区域四面体单元数量为 367077，整个模型的四面体单元数量为713753，反演初始模型为0.01 S·m-1的均匀半空间.反演收敛情况如图12所示，地面TEM、SATEM单独反演和联合反演最终结果的数据归一化拟合差分别为8.3%、5.4%和10.2%.

图9、10和11为带地形模型反演结果的不同方向切片图.在图9b、10b和11b的地面TEM单独反演结果中，地面中心回线观测方式对该模型的反演结果很不理想，两个高导异常体连在一起，无法恢复高导体的准确信息.两个低导异常体的恢复效果较差，只能看到模糊的大致情况.我们认为是因为中心回线装置的观测数据中缺少不同偏移距离的信息导致难以区分该模型的两个高导异常体.在图9c、10c和11c的SATEM单独反演结果中，四个异常体的恢复结果较好，对于底部埋深约有600 m的倾斜高导体也恢复了异常体的总体形态，但在异常体附近存在一些小的假异常.图9d、10d和11d的联合反演结果中，四个异常体的形态恢复较好，SATEM单独反演的局部假异常也被有效弱化.在X=-1500 m附近出现了新的假异常，其产生原因可能是由于此处不在观测数据的覆盖范围内，数据对该区域灵敏度较小.此外，该假异常较之单独反演的假异常更为微弱，对数据解释影响较小.

图9 模型2在Z=100 m的模型与反演结果切片(a) 理论模型； (b) 地面TEM反演； (c) SATEM反演； (d) 联合反演.Fig.9 Slice view of the true model and the inversion results at Z=100 m for model 2(a) Synthetic model; (b) The ground TEM inversion; (c) SATEM inversion; (d) Joint inversion.

图10 模型2在Y=-200 m的模型与反演结果切片(a) 理论模型； (b) 地面TEM反演； (c) SATEM反演； (d) 联合反演.Fig.10 Slice view of the true model and the inversion results at Y=-200 m for model 2(a) Synthetic model; (b) The ground TEM inversion; (c) SATEM inversion; (d) Joint inversion.

图11 模型2在Y =200 m的模型与反演结果切片(a) 理论模型； (b) 地面TEM反演； (c) SATEM反演； (d) 联合反演.Fig.11 Slice view of the true model and the inversion results at Y = 200 m for model 2(a) Synthetic model; (b) The ground TEM inversion; (c) SATEM inversion; (d) Joint inversion.

图12 模型2反演的收敛曲线Fig.12 Convergence plot for different inversions of model 2

5 结论

本文提出了一套三维TEM反演算法.正演算法采用基于非结构化四面体网格的矢量有限元法，可以对起伏地形、复杂地质结构、任意形状发射源和任意波形进行精确模拟.在三维正演算法的基础上开发出了基于高斯-牛顿法的高效并行反演算法；构建了地面TEM和SATEM的联合反演方案.该算法可适用于不同TEM系统的单独反演及不同系统的联合反演.通过理论模型研究验证了算法的有效性.在模型1中，对比了地面TEM点距50 m数据和点距200 m数据的单独反演，点距200 m数据反演结果分辨率较高、假异常较少，说明了增加TEM数据量可以一定程度上提高反演结果质量.对比两个模型的地面TEM与SATEM的单独反演可以发现，单一的TEM探测方式的探测能力具有局限性，单独反演通常伴随一定的假异常，且不同的探测装置的单独反演结果产生不同的假异常，给解释带来巨大的挑战，为获得可靠结果有必要开展联合反演.联合地面TEM和SATEM数据联合反演，可以结合不同观测方式的探测优势，以及压制由于反演多解性而形成的假异常，反演结果的可靠性和分辨率得到提高.本文所提出的方法算法可以为复杂地质条件下开展地面TEM与SATEM探测提供理论支撑.