1 引言

矩形大定源是目前瞬变电磁法(TEM)中重要的野外工作装置(Fitterman and Stewart, 1986; Danielsen et al., 2003).但该观测方法资料反演依然是制约其应用的一个难题，因此目前在工程应用中, 主要基于视电阻率(Spies and Eggers, 1986; 白登海等, 2003; 张成范等, 2009)和一维反演的模型进行资料解释(Xue et al., 2004, 2012; 翁爱华等, 2008; Zhou and Xue, 2014; Li et al., 2016).在三维研究方面，大多数研究围绕三维正演模拟进行，SanFilipo和Hohmann(1985)首先采用积分方程法，后来Newman和Hohmann(1988)基于频率域电磁场模拟和时频转换技术，完成三维正演工作, 并且与SanFilipo和Hohmann(1985)的结果对比，验证算法的正确性; Wang和Hohmann(1993)采用时间域有限差分法并考虑模型空间磁导率的变化.在国内，孙怀凤等(2013)应用时间域有限差分法计算了矩形回线瞬变电磁响应，并将发射机关断效应考虑到模拟过程中; 李建慧等(2013)开展了利用G-S变换将矢量有限元得到的频率域三维响应变换为瞬变响应的研究.这些工作为瞬变电磁数据三维反演奠定了基础.

但直接用于大定源瞬变电磁数据反演的技术尚处于研究中，因此矩形大定源瞬变电磁三维反演解释研究同样远远落后于应用需要.虽然其他瞬变电磁反演的技术，从原理上也可以应用到本方法的资料反演中，如Wang等(1994)采用共轭梯度法完成井中瞬变电磁数据的三维时间域反演；Newman和Commer(2005)改进了Wang等(1994)和Zhdanov和Portniaguine(1997)的算法，在大型并行计算平台上实现了长偏移距瞬变电磁(LOTEM)数据的三维反演; 基于Haber等(2007)的工作，Oldenburg等(2013)利用线性方程组直接求解软件(MUMPS)在并行计算平台上解决大尺度多源时间域电磁数据三维反演问题.但直接将上述时间域反演方法应用于大定源TEM方法中存在困难，最主要的困难是三维大定源瞬变正演的实现及与反演技术的精巧结合.然而Mitsuhata等(2001)成功地从LOTEM数据中恢复出频率域数据，并应用频率域三维反演技术进行反演，获得良好效果(Sasaki et al., 2015).最近，Weng等(2017)完成了时间域到频率域的反演变换，并应用到大定源瞬变电磁数据到频率域响应转换.因此, 若能进一步完成大定源频率域电磁数据的三维反演技术，将可从频率域角度完成大定源瞬变电磁数据的三维解释.

为此，本文首先概述非线性共轭梯度反演算法的基本原理，然后采用交错网格有限差分技术计算三维模型频率域正演响应，且与Newman和Hohmann(1988)给出的结果对比，验证算法的正确性.最后分别给出的单个低阻异常体和不同埋深的双低阻异常体模型合成的数据进行的三维频率域反演，验证了算法的可行性和有效性.因此，基于本文算法，利用时频转换技术，可以对野外实测的矩形大定源瞬变电磁数据进行三维频率域反演，为瞬变电磁数据的三维解释奠定了基础.

2 三维频率域正演 2.1 基本原理

为了克服地面矩形大定源电磁法在三维数值模拟过程中源描述的困难以及源奇异性对数值稳定性的影响，将其作用用背景场E^p代替.这样，三维正演数值模拟将针对异常二次电场E^s进行.取正时谐，忽略位移电流，二次场E^s满足(Commer, 2003)

(1)

式中k_p²=-iωμ_p(σ_p+iωε)，σ为电导率；μ和ε分别为背景模型的磁导率和介电常数，在模拟中近似等于真空中磁导率和介电常数.

式(1)采用交错网格有限差分技术离散(Yee, 1966; Alumbaugh et al., 1996)，最后形成如下矩阵方程：

(2)

其中K为大型系数矩阵，s为含有背景场和边界条件的右手项.由于剖分时，边界网格足够大，模拟区域中局部异常在最外边界处产生的二次电场的切向分量近似为零.该方程的求解采用最小残差方法(Mackie et al., 1994)，利用散度校正提高计算精度和收敛速度.

大定源产生的背景场E^p的计算采用如下的思路：首先采用面积更小的方形虚框覆盖矩形大定源回线(翁爱华等，2010)；然后，借助虚界面方法(Weng et al., 2016)计算方形虚源等效的各个圆回线在全空间中任意点的电场，叠加得到大定源回线产生的全空间背景场.

一旦计算出电场E^s，有计算网格边界中点总场E=E^s+E^p，进而可以通过

(3)

求解网格单元表面中心处磁场(Kelbert et al., 2014)，而具体计算点处的电磁场则由这些电场或者磁场通过空间插值得到(翁爱华等, 2012).

2.2 正确性验证

三维大定源频率响应模拟精度验证采用如下的验证思路：将数值模拟的频率响应利用余弦变换(Patra and Mallick, 1979)进一步变换为时间响应，并将其与前人瞬变响应的精度比较反推频率响应的模拟精度.精度验证的模型来自Newman和Hohmann(1988)，即均匀半空间嵌入一个三维低阻异常体模型，具体的异常体大小和观测装置如图 1所示.数值模拟时，以异常体中心在地面的投影为坐标原点，网格剖分数量为25×19×20个，最小网格大小为20 m×20 m×20 m.首先计算[-200, 200]区间内的测点的频率域响应，然后利用余弦变换获得感应电动势(感应电动势数据被1 m²的单匝接收线圈归一化).选取1 ms和30 ms两个时间道的计算结果与Newman和Hohmann(1988)的结果对比，结果如图 2所示.比较可见, 两种算法的结果具有较好的一致性，进而我们认为本文的频率响应是正确的.

2.3 近源的频率响应

大定源瞬变电磁是近源电磁法，从观测时间上说属于晚期瞬变电磁测量；因此，将其瞬变响应变换到频率域，其频率响应低频成分明显，在一级近似条件下，场不随频率变化，可能不具有交流测深的功能.为了探讨这个问题，分别计算了图 1中异常体上方和外侧两个位置处有无异常体时的频率响应曲线，如图 3所示.从图可见，在准静态近似条件的频率范围内，响应实部确实与频率无关，但虚部与频率相关；并且在没有异常体时，两个测点的频率响应特征基本一样，其频率响应符合均匀半空间变化规律.当存在异常体时，两个测点频率响应的实部和对应的均匀半空间基本一样，但虚部在低频时完全偏离均匀半空间响应，高频基本没有变化.这个特征表明，图 1中模型异常体的存在只反映在相对低频的测深响应上；在高频时，探测深度没有达到异常体位置，故没有响应异常.因此，从频率响应可以得到如下结论，近区瞬变电磁法具有频率测深功能；只是这种响应特征主要通过虚部响应反映出来.传统意义上的近区不具有测深功能是从场强度，在近似条件下得到的，而测深感应被掩盖在比实部小很多的虚部响应中.在这个意义上，如果在频率域测量，则要求非常高的相位测量精度；而瞬变电磁通过测量感应电动势，以仪器的灵敏度换取测量精度.因此，虽然大定源瞬变电磁是近区测量，但从时间域变换到频率域，是可以获得地下深部信息的.

3 反演方法 3.1 基本流程

地面矩形大定源三维频率反演采用非线性共轭梯度方法, 该方法通过计算非线性目标函数

(4)

的梯度解决正则化反演问题(Avdeev, 2005; Egbert and Kelbert, 2012; Siripunvaraporn, 2012).式中m为模型参数，m_ref为参考模型矢量，为正演算子；d为数据矢量，而C_d和C_m分别为数据方差矩阵和模型方差矩阵；λ为正则化因子，平衡数据预测误差和模型约束.一个典型的反演流程(翁爱华等, 2012; 刘云鹤和殷长春, 2013)经过简化后示意于表 1中.

3.2 梯度的计算

根据Egbert和Kelbert (2012)，目标函数的梯度g可以表示为

(5)

式中，为灵敏度矩阵.可以看出，计算梯度g，只需要计算出灵敏度J后直接代入(5)式即可.但显式计算灵敏度矩阵，需要大量正演计算和内存需求(Avdeev, 2005；Siripunvaraporn, 2012).因此根据Egbert和Kelbert(2012)，进一步将灵敏度矩阵J及J^T表示为

(6)

及

(7)

这里，K为正演系数矩阵.，.

这样，从将梯度计算中J或者J^T与向量乘积转换为系数矩阵K或者K^T与相关矢量的乘积，即所谓的正演或者伴随正演(Egbert and Kelbert, 2012)，从而只需要2次正演的计算量，并且不需要显式计算并存储灵敏度矩阵J.这是非线性共轭梯度方法成功应用于三维电磁数据的反演的最大特点.

4 反演结果 4.1 模型1

图 4所示为100 Ωm的均匀半空间嵌入一个1 Ωm的低阻异常体模型.异常体大小为80 m×280 m×60 m，顶部埋深60 m.发射回线为300 m×600 m，位于异常体左侧，发射电流1 A.异常体中心在地面上的投影点为坐标原点.回线内布设14条测线，每条测点上14个测点，点距和线距均为40 m.在0.125~8192 Hz频率区间按照2的幂次方布置17个频点.模型剖分个数为18×19×20，最小网格的大小为40 m×40 m×10 m，x、y正负方向均有5个扩边网格，z正方向有5个扩边网格，扩边系数分别为2、2和1.5.

反演在PC上进行串行运算，Intel(R) Core(TM) i5-2320处理器，主频3.0 GHz，4 G内存.数据中加入5%的高斯随机噪声，均匀半空间模型作为反演的初始模型，迭代120次，用时4个小时.图 5给出了具体的反演结果.图 5a为深度75 m处的平面切片，该平面横切低阻异常体，具有较高的横向分辨率.过异常的x和y方向的断面图分别如图 5b和5c.从图可以看出，异常体的位置与理论模型一致，证明了矩形大定源非线性共轭梯度三维频率域反演的可行性和有效性.图 5d的迭代过程中均方根误差(RMS)的变化情况显示，随着迭代的进行，数据拟合误差逐渐减小，迭代稳定收敛，但在一定迭代次数后，误差减小变慢，收敛变慢.传统的大定源方法一般只处理中间一定区域的数据.在本例中，接收点位于矩形回线内外的任意位置，体现了本文的算法可以处理任意位置数据的优点.

4.2 模型2

图 6所示为100 Ωm的均匀半空间嵌入两个埋深不同、电阻率为5 Ωm的低阻异常体模型(Sasaki et al., 2015).异常体大小为400 m×800 m×500 m，埋深分别为200 m和400 m，工作布置与Sasaki等(2015)的一致，点距和线距都为400 m，构成9×5的测网.矩形大定源发射回线的边长为1000 m×1000 m，供电电流1 A.反演数据中加入5%的高斯随机噪声，从电阻率为100 Ωm的均匀半空间开始，经过110次迭代，得到的三维频率域反演的结果如图 7.从图可以看出，地面矩形大定源非线性共轭梯度三维频率域反演较准确地揭示出异常体的电性特征和空间位置.图 8的RMS曲线显示NLCG反演算法在前期收敛速度较快，但后期收敛速度变慢.显然，要想到达更高的拟合精度，需要较多的迭代次数，这也是非线性共轭梯度算法的弱点.值得指出的是，在本例中，对接收点位置没有要求，观测可以在矩形回线内外的任意位置上进行.在这点上，本文提出的正演和反演算法，可以实现矩形大定源全区观测数据的反演.

5 结论

本文通过数值模拟，研究了矩形大定源频率数据的非线性共轭梯度反演技术，并用理论模型验证了该算法的可行性和反演效果.通过数值分析，可以获得如下的结论：

(1) 采用交错网格有限差分针对二次电场满足的Helmholtz方程实现矩形大定源频率响应的三维数值模拟.模拟过程中，利用虚界面法结合大定源虚框等效原理获得背景场，克服了大定源数值模拟过程中源描述的困难以及源奇异性；与前人结果的一致性验证了本文正演算法准确性；

(2) 采用非线性共轭梯度反演方法，能由大定源频率电磁响应，获得地下电阻率的三维分布，从数值上验证了本文反演算法的可行性；

(3) 在数值计算过程中，参与反演计算的测点分布范围突破传统的回线中间1/3范围的要求，而是分布在整个回线内，甚至回线外.因此，本文提出的反演算法，具有更大的灵活性，并且可以提高大定源瞬变电磁的野外工作效率.