航空电磁作为一种快速有效的勘查手段，在矿产资源、地下水和环境工程领域得到越来越广泛的应用.三维航空电磁正演技术目前已相当成熟，如采用规则网格的有限差分法(Newman and Alumbaugh, 1995;陈辉等，2011)，非结构网格有限元法(Yin et al., 2016a, 2016b; Ren et al., 2013)等.近年由于“Moving Footprint”概念(Cox et al., 2010)和矩阵直接求解技术(Oldenburg et al., 2013)的引入，三维航空电磁反演也取得了相当大的进展，如基于积分方程(Cox et al., 2012)和有限差分法(Liu and Yin, 2016)的三维高斯牛顿反演已应用到航空电磁实测数据的反演解释之中.然而，航空电磁数据量巨大，地球物理电磁二、三维反演由于计算硬件要求高和计算效率低，目前尚没有完全实用化.航空电磁采样密集，同时系统的“Footprint”范围较小(Yin et al., 2014)，相邻测点的地下电性差异不大，特别是在一些沉积岩地区，一维模型能很好地模拟地下电性分布特征，因此航空电磁数据处理和解释目前仍然普遍采用一维模型成像和反演方法.

航空电磁数据一维反演常用的是对每个测点单独进行反演.周道卿等(2010)将马奎特反演应用到频率域航空电磁数据反演中，罗勇等(2014)将一维阻尼特征参数反演应用于航空电磁数据反演中.然而，由于反演解的不稳定性，单个测点一维反演结果中经常出现电阻率或层厚度沿横向方向发生突变的现象，造成反演的地电断面粗糙和不连续.因此，在地下介质横向连续性较好的地区(比如沉积岩地区)，施加横向约束进行多测点联合反演是一个很好的解决方案.

Viezzoli等(2008)指出，为获取连续性较好的反演结果可以从数据域和模型域两方面入手.从数据域就是在反演前对观测数据进行光滑.该方法过去被广泛应用于航空电磁数据处理中.数据经光滑后信噪比增加，但分辨率降低.从模型域的角度看，就是在反演过程中对相邻测点的模型参数施加约束，不需要对数据进行光滑，有效地保留了数据中详实的地下导电率信息(Viezzoli et al., 2008).本文研究针对模型域的约束反演方法.

横向约束方法已被用于不同的地球物理数据反演之中(Auken et al., 2000).Auken等(2002)系统总结了横向约束反演思想，并将其成功运用于直流电法数据反演.对于二维电阻率和三维瞬变电磁数据同样可利用横向约束进行反演(Auken and Christiansen, 2004；Auken et al., 2004).Siemon等(2009)和Tartaras和Beamish(2006)将横向约束反演推广到航空电磁领域，分别对直升机和固定翼航空电磁数据进行反演.蔡晶等(2014)在频率域航空电磁数据中引入横向加权约束，获得很好的反演结果.殷长春等(2016)将加权横向约束应用到时间域航空电磁数据反演中，并通过对CGG航空电磁系统HELITEM MULTIPULSE的实测数据反演验证了方法的有效性.

传统的横向约束反演方法并没有考虑相邻测线之间的构造关联性，因此往往只能给出沿测线方向连续性较好的反演结果.由于测线之间没有约束，垂直于测线的地质信息难以在数据反演中应用，反演得到的电阻率沿测线方向的关联性没有得到体现，给地下构造解释带来不确定性(Viezzoli et al., 2009).此外，必须注意的是，垂直测线方向往往构造连续性较好，因此施加测线间连续性约束可大大改善反演效果.Viezzoli等(2008)提出了空间约束反演(spatially constrained inversion，SCI)，即沿测线和垂直测线方向同时施加约束，利用航空电磁数据验证算法有效性，并将该算法成功运用于澳大利亚墨累河下游地区的环境评估.近年，Santos等(2011)将拟三维电导率成像运用于DUALEM-421数据反演，而Vignoli等(2015)提出了一种基于梯度正则化的锐空间约束反演方法(sharp spatially constraint inversion).

传统空间约束反演通常采用阻尼最小二乘法.这种基于最小二乘准则的反演方法，是通过最小化由数据残差构成的目标函数求取模型“最优解”.反演过程只追求模拟数据与原始观测数据的最大拟合(周道卿等，2010).当利用最小二乘进行空间约束反演时，为了克服奇异值分解以及搜索阻尼因子时伴随正演造成的计算效率问题，通常将大区域的数据集划分为许多小子集，然后以每个子集作为一个单元进行空间约束反演(Viezzoli et al., 2008).反演过程中各测点的约束按照Delaunay三角化准则进行.当数据量庞大时这种反演策略耗时巨大.为提高拟三维空间约束反演效率，减少内存需求，本文采用有限内存L-BFGS最优化方法，由于无需显式存储Hessian逆矩阵和灵敏度矩阵，缩减了内存需求，同时保证目标函数有更好的收敛性.本文约束反演算法与前人采用的纯数学意义上的约束反演(比如基于Denaulay准则)存在本质差异：本文算法考虑的是地质意义上的约束.它既考虑沿测线方向相邻测点间的约束(基于稠密采样)，又考虑垂直测线方向构造连续性约束.通过理论和实测数据实验测试表明，本文设计的反演算法计算效率高、稳定可靠，非常适合求解多测线、大数据量的航空电磁最优化反演问题.

航空电磁横向约束反演采用的是水平层状模型.假设时间域航空电磁系统使用圆形回线作为发射源，在回线中心或正上方测量随时间变化的磁感应强度或感应电动势.定义直角坐标系的原点为圆形回线中心在地表的投影，x, y平面位于大地表面，z轴向下为正.发射回线中的电流强度为I，半径为a，高度为h.对于一维层状各向同性介质，根据Nabighian(1988)，准静态条件下回线中心或正上方频率域磁场垂直分量表达式为

(1)

式中，J₁(λa)为一阶贝塞尔函数，r_TE是地下电性相关的反射系数，其递推计算参考(Nabighian，1988).公式(1)中包含的贝塞尔函数在(0，∞)区间上的积分可以用快速汉克尔变换计算.

根据频率域和时间域电磁信号的关系，可以通过傅里叶变换将宽频范围内的频率域电磁响应转换到时间域(殷长春等，2013)：

(2)

其中，ω是角频率，t是时间，F(ω)为频率域响应，B_S(t)是时间域阶跃响应.由频率域电磁场实虚部的奇偶性并结合余弦函数和0.5阶贝塞尔函数的关系可得

(3)

对于时间域航空电磁任意发射波形，使用由殷长春等(2013)给出的褶积公式计算时间域电磁响应，即

(4)

式中，*表示褶积，I(t)是随时间变化的电流脉冲，B(t)是磁感应强度，dB(t)/dt是磁感应的时间变化率(感应电动势).由此，可以通过阶跃响应和发射波形褶积计算任意发射波形的时间域航空电磁响应.

为了描述沿测线方向和垂直测线方向模型参数之间的约束关系，本文首先建立一个多层介质模型.进而，基于横向约束反演思想设计两种模型方案：一种方案类似于Occam反演的思想——假设层数和层厚度固定，仅反演电阻率；另一方案是层数固定，反演地下介质电阻率与厚度.通过在目标函数中引入约束项，实现航空电磁数据的拟三维模型空间约束反演.

对于一个工区包含多条航空电磁测线，首先依据测点的地理坐标建立一个统一的拟三维的模型网格，如图 1所示，其中x轴方向表示测线方向，y轴垂直于测线方向，z轴垂直向下.测点和网格中心一一对应(当测线数据不是规则网格时，可先进行网格化处理).每个测点都使用N_z层模型，参考瞬变电磁场的扩散规律我们将剖分的层厚度按指数方式递增.假设沿测线方向有N_x个网格(N_x为测点个数)；垂直于测线方向有N_y个网格(N_y为测线条数)，m_{i, j}^k表示第(i, j)个网格节点的第k个模型参数.必须指出的是，虽然我们定义三维网格模型，然而本文所有正演模拟和雅克比矩阵的计算均使用一维模型，即各测点电磁响应仅与其对应网格节点的垂向电性分布有关，不受周围网格电性的影响.相邻测点电性之间的关联通过双向空间约束实现，因此我们称之为“拟三维”模型空间约束反演.当一次只反演一条测线数据时，本文的算法简化成传统的横向约束反演(Laterally Constrained Inversion, LCI).

图 1

Fig. 1

图 1 拟三维模型 Fig. 1 Quasi 3D model

对于本文的拟三维模型空间约束反演，将正则化目标函数定义为

(5)

其中d为观测数据，m为模型参数向量，依据约束方案分别由各层电阻率或电阻率和层厚度组成.C_d是数据方差矩阵，f(m)为模型响应算子.(5)式中第二项为本文采用的反演模型约束项，引入正则化因子λ可在数据拟合和模型约束之间取得平衡.

对(5)式求导可得目标函数的梯度：

(6)

其中，约束矩阵ψ(m)由各测点的模型参数约束项构成ψ_{i, j}^k (i=1, 2, …，N_x, j=1, 2, …, N_y, k=1, 2, …，2N_z-1(不固定层厚)或k=1, 2, …，N_z(固定层厚)).J为雅克比矩阵，本文采用差分方法计算，通过和半解析解对比试验发现，采用3%模型改变量可满足数值计算精度要求.

对于公式(5)和(6)中的正则化项，采用对每个网格节点的模型参数及其周围网格(前、后、左、右)的模型参数施加约束，则第(i, j)个网格的第k个模型参数的约束项可表示为

(7)

式中， 分别表示沿着x、y方向施加约束的相对权重.针对边界测点无法使用全部约束项的情况，采用该测点周围可用的部分测点作空间约束，此时(7)式中的约束项进行相应的取舍.在实际进行空间约束反演时，为计算(6)式中的目标函数梯度，需要计算模型约束项对模型参数m_{i, j}^k的导数.由(7)式可得

(8)

通过选择合适的权重，使相邻测点的模型参数差的加权平均最小化，实现相邻测点模型参数相互约束.由于每个测点又与其他相邻测点进行约束，使得模型参数信息沿测线和垂直测线两个方向按照距离权重依次向外传递直至覆盖全局，从而整个测区成为一个关联的整体，实现了整个测区拟三维模型空间约束反演.

在完成模型约束后，我们可以建立航空电磁拟三维模型空间约束反演的技术流程(刘云鹤等，2013).本文采用Nocedal和Wright(1999)提出的L-BFGS步长搜索技术进行反演.假设目标函数存在二阶导数，对其在已知点m_l处进行泰勒展开并忽略高阶项可得

(9)

其中，H_l为Hessian矩阵，p_l为模型搜索方向，φ_l=φ(m_l).令，可得

(10)

由此，模型参数迭代由如下形式给出：

(11)

式中α_l为迭代步长.从一个起始模型出发按照公式(10)进行迭代，直到数据拟合差小于事先设定的阈值.Hessian矩阵的逆矩阵迭代公式为

(12)

(13)

(14)

由上式可知，通过存储n个最近的向量组s_i, y_i, i=l－n, …, l－1对H_l^－1进行校正，即可生成Hessian矩阵的逆矩阵.该向量组包含了n个最新迭代的曲率信息，实践证明n值通常取3~20即可迭代得到满意解.

由此，可得L-BFGS方法实现拟三维模型空间约束反演的技术流程如下：

(1) 针对初始模型利用(7)式计算约束项；

(2) 利用(6)式计算目标函数的梯度；

(3) 假设一个初始对称正定矩阵(通常为单位矩阵)，利用(11)—(13)式递推计算Hessian矩阵的逆矩阵；

(4) 利用公式(9)获取搜索方向p_l，再利用(10)式进行模型更新，同时利用(5)式计算目标函数值并判断是否满足要求.如满足要求，迭代终止，否则利用新模型替代初始模型，转入第(1)步.

应用L-BFGS方法进行迭代时，线性搜索的关键在于如何确定搜索方向p_l及搜索步长α_l.对于搜索步长α_l的选择，我们通常采用如下充分下降条件和曲率条件，即

(15)

(16)

其中0 < c₁ < c₂ < 1，c₁通常取值很小，可取为c₁=10^－4，c₂通常取值较大，本文取为0.9.公式(14)和(15)合称为Wolfe条件.

有限内存的拟牛顿法只需计算目标函数的一阶导数，无须计算二阶导数，因而计算量大大减小.同时由于计算Hessian矩阵的逆矩阵时只需存储有限的早期迭代信息，大大节省内存需求，特别适用于大数据量的航空电磁数据反演.

本节讨论拟三维空间约束反演对理论和实测数据的应用效果.我们考虑Fugro/CGG公司的HELITEM MULTIPULSE系统.HELITEM系统发射双脉冲，即先发射一个大磁矩的半正弦波，然后再发射一个小磁矩的梯形波，如图 2所示.HELITEM系统采用中心回线，发射线圈高度30 m, 接收线圈高度30.5 m，发射磁矩约39×10⁴ Am².

图 2

Fig. 2

图 2 航空电磁HELITEM MULTIPULSE系统发射波形 Fig. 2 Transmitting waveform for HELITEM MULTIPULSE system

本文假设电阻率为100 Ωm的半空间中埋藏一个缓起伏三维导电层，电阻率为20 Ωm，其下底面水平，上底面埋深从20 m到80 m变化，模型主剖面参见图 3a.由于山坡状良导层变化非常缓慢，相邻测点之间电性变化较小，可采用拟三维空间约束方法进行航空电磁数据反演.正演模拟采用时域有限元法(殷长春, 2016c)，累计设计21条测线，线距50 m，点距10 m，共计2121个测点.我们对每个测点每个时间道的电磁响应添加3%的高斯噪声.

图 3

Fig. 3

图 3 加入3%高斯噪声的不同测线HELITEM理论数据反演结果 (a)—(c)理论模型; (d)—(f)单点反演; (g)—(i)固定层厚SCI反演; (j)—(l)可变层厚SCI反演 Fig. 3 Inverted sections of synthetic HELITEM data with 3% Gaussian noise added (a)—(c) Theoretical model; (d)—(f) single-station inversions; (g)—(i) SCI inversion with fixed layer-thickness; (j)—(l) SCI inversion with variable layer-thickness.

反演使用的初始模型层数设为20层，第一层的层厚度为5 m，以后各层按1.12倍递增.我们首先取w_x=w_y=1沿测线方向前后和垂直测线方向左右各取一个点，对当前测点实施空间约束.对于固定层厚的拟三维反演方案，在利用L-BFGS方法进行迭代反演时，仅更新电阻率，其他模型参数不变；而对于可变层厚的拟三维反演方案，则同时更新电阻率和层厚度.由于两种空间约束反演模式的参数约束范围不同，正则化因子的影响也不同.本文首先通过试验(见后面讨论)确定如下最优正则化因子：对于层厚可变的空间约束反演，正则化因子取为50；而对固定层厚的空间约束反演，正则化因子取为10.图 3给出不同测线理论模型数据进行单点反演、固定层厚和可变层厚SCI的反演结果.

由图 3可以看出：(1)无论对于主剖面L1还是旁侧剖面L4和L7，单点反演结果受噪声影响严重，存在明显的电阻率突变、界面粗糙.与理论模型对比发现，拟合界面效果差；(2)由于施加空间约束，两种SCI反演方案给出的反演断面的横向连续性均有很大改善，与理论模型吻合较好.其中，采用可变层厚的反演策略给出的结果更光滑、更接近真实模型.因此，相比于单点反演跳跃严重的现象，两种空间约束反演策略对噪声有较明显的抑制作用，反演的界面均比较光滑，其中可变层厚度反演结果最理想.

为了更清楚地了解拟三维空间约束反演的效果，图 4展示上述三维模型反演结果在z=45 m处的水平切面.图 4a的理论模型中低阻区为一个标准圆.对比三种反演结果可以看出：(1)受噪声影响，单点反演结果切面形状极不规则，低阻分布零散，边界严重畸变；(2)两种拟三维反演给出的切面基本反映理论模型的形状及范围，边界连续性较好.其中可变层厚的空间约束反演最接近真实模型.

对于本模型的21条测线和2121个测点，假设模型剖分层数20层，在计算机(3.40 GHzCPU和8 GB内存)上利用单点反演耗时10小时47分，利用固定层厚SCI反演耗时4小时18分，而利用可变层厚的SCI反演耗时7小时13分.可以看出SCI反演明显具有精度高和速度快的优点.

图 4

Fig. 4

图 4 HELITEM理论数据反演结果切面图 (a)理论模型；(b)单点反演；(c)固定层厚SCI反演；(d)可变层厚SCI反演. Fig. 4 Planar view of inversions of synthetic HELITEM data (a) Theoretical model; (b) Single-station inversions; (c) SCI inversion with fixed layer-thickness; (d) SCI inversion with variable layer-thickness.

为分析权重对空间约束反演效果的影响，我们以可变层厚反演策略为例，给出两个方向施加不同权重时HELITEM理论数据的反演结果.对比图 5中的结果发现，当w_x与w_y的比值变大时，层界面起伏变大，特别是底界面形状发生明显畸变.比值越大，界面畸变越严重.在实际应用时，可根据点线距相对大小调整w_x与w_y的比值，以达到最佳反演效果.本文中我们对理论和实测数据均选取w_x=w_y=1.从图 5可以看出，反演出的中间层顶、底界面均比较光滑.

图 5

Fig. 5

图 5 不同加权因子理论数据反演结果 Fig. 5 Inversions of synthetic data for different weightings

正则化因子λ的选择是成功实现SCI反演的关键因素之一.λ选的过大导致数据拟合较差，反演结果严重失真；反之，λ选的太小空间约束效果不明显.为获得最优λ值，我们以可变层厚SCI反演为例说明如何确定最优正则化因子.图 6给出不同λ时模型主剖面和旁侧剖面的反演结果.由图可以看出：当λ=0时，即无约束条件下对所有数据同时反演时(实质上是测区所有数据联合反演)，模型剖面的边界相对粗糙；当λ=50时，边界光滑连续，旁侧剖面效果更明显；当λ继续增加到>100时，边界连续性仍然较好，但低阻层底界面发生畸变，旁侧剖面效果更差.因此，本文针对可变层厚的SCI反演选择λ=50作为正则化因子.类似地，针对固定层厚的SCI反演，我们选定λ=10为正则化因子.

图 6

Fig. 6

图 6 不同正则化因子理论数据反演结果 Fig. 6 Inversions of synthetic data for different regularization factors

为进一步比较各种不同反演方法的有效性，我们对加拿大某区利用HELITEM MULTIPULSE系统采集的航空电磁数据进行反演.如前所述，HELITEM MULTIPULSE系统在前端发射一个高能量的半正弦波，在末端再发射一个快速关断的方波或梯形波，具有勘探深度大、近地表分辨率高的特点(殷长春等，2015).测区累计飞行观测6条测线，线距250 m，每条测线有248个测点，点距40 m.假设反演层数为12层.图 7给出单点SVD反演、沿测线方向的横向约束LCI反演以及固定层厚和可变层厚条件下拟三维空间约束SCI反演结果.由图可以看出：(1)单点反演结果电阻率发生突变，横向不连续性现象严重；(2)沿测线方向的LCI反演结果中，中间良导层的连续性有一定改善，但高阻基底电阻率的连续性较差，特别是下半空间电阻率呈明显条带状；(3)固定层厚的SCI算法对于地下主要构造反映效果明显，对整个测区的电阻率层状分布有较清晰的展示，然而电阻率横向连续性仍然不理想；(4)可变层厚的SCI算法相比于单点反演和固定层厚SCI反演，界面更光滑，电阻率横向连续性更好，与测区的地质特征吻合很好(Chen et al., 2014).

图 7

Fig. 7

图 7 HELITEM MULTIPULSE系统实测数据反演结果 (a)单点反演；(b)横向约束LCI反演；(c)固定层厚SCI反演；(d)可变层厚SCI反演. Fig. 7 Inversion results of survey data from HELITEM MULTIPULSE system (a) Single-station inversions; (b) LCI inversion; (c) SCI inversion with fixed layer-thickness; (d) SCI inversion with variable layer-thickness.

本文基于空间约束思想设计两种约束模式，并利用L-BFGS成功实施航空电磁拟三维模型空间约束反演，得出如下结论：

(1) 空间约束反演沿着测线和垂直测线两个方向对测点电性参数施加约束，通过邻近测点传递信息，从而将所有测点联结成为一个整体，构建了一种精确和高效的航空电磁数据反演方法.较之于单点和LCI反演，两种空间约束方案均能获得更好的反演结果.

(2) 本文提出的两种反演约束反演方案中，固定层厚约束反演仅能对电阻率施加横向约束.相比之下，可变层厚约束反演既可对层厚又可对电阻率施加横向约束，反演效果更好.

(3) 通过对比反演时间可以看出，当数据量较大时，基于L-BFGS算法的反演策略在保证反演精度的前提下具有速度快、计算效率高的优势.

(4) 正则化因子和加权系数对约束反演效果有较大影响，实际数据反演是应通过大量模型试验确定.

考虑到本文模拟飞行测区和实测数据量均较小，将整个测区数据同时进行反演.对于测区较大的情况，受计算条件的限制，可能无法完成庞大数据的一次性反演.此时，可将测区划分成若干个小区块，之间保持适当重叠，在各小数据区块完成约束反演后，通过加权平均进行参数融合，实现整个测区的大数据反演.另外，在本文中仅讨论了比较简单的利用测点周围4个相邻测点进行约束的情况，对于更复杂的约束(比如多点约束)或许能更好地改善反演效果，这些将是我们未来的研究方向.

感谢CGG公司的陈天友博士为本研究提供HELITEM MULTIPULSE实测数据.作者向吉林大学电磁“千人计划”研究团队成员，特别是刘云鹤副教授在本文算法和后续文章准备过程中提供的帮助表示感谢.