0 引言

时间域航空电磁法具有速度快、成本低、通行性好、勘探深度大等优点，特别适合我国地形和地质条件复杂的地区，比如沙漠、沼泽、高山以及植被覆盖区^[1]，目前已广泛应用于金属矿产、地下水勘查、地质填图及环境工程等领域^[2]。由于数据采样密集，在具有很高的横向分辨率的同时，其海量数据也给电磁二、三维反演方法带来了挑战，因此目前实用的数据处理方法主要是电导率深度成像和一维反演。

电导率深度成像并不属于反演方法。它是把观测数据转换为如视电导率的中间参数。这些参数表征地下电性结构的主要特征。通过将这些中间参数和转换深度联系起来，就可以对地下结构进行成像。目前时间域航空电磁法主要应用的成像方法有EMFlow^[3]、电导率深度转换^[4]、基于半空间模型^[5]和假层半空间模型^{[6, 7, 8]}的查表法等。嵇艳鞠等^[9]还研究了基于多元线性回归的直升机时间域航空电磁数据三维异常体电导率深度识别算法。

成像方法速度快，然而精度较低，比较适合于海量数据的现场、快速处理。相比之下，反演方法能获得更加可靠的地下电阻率分布结构。黄皓平和王维中^[10]将广义逆反演方法应用到二层大地时间域航空电磁数据的反演；罗勇等^[11]研究了时间域航空电磁一维阻尼特征参数反演方法。这两种反演方法都属于梯度方法，反演速度快，但是对初始模型要求高，容易陷入局部极小。强建科等^[12]将Occam反演引入到航空瞬变电磁一维反演中。Occam方法^[13]应用了正则化反演的思想，由于在目标函数中加入了模型光滑约束项，造成反演结果中地层分界面不清晰。Sattel^[14]还尝试了Zohdy反演方法。

上述成像和反演方法都假设地下为一维层状介质模型。由于航空电磁系统数据量巨大，大规模二、三维反演耗时太长，难以实现^[15]，时间域航空电磁数据反演尤其困难。考虑到航空电磁系统空间采样非常密集(点距通常为几米)、具有很高的横向分辨率，相邻测点每一层的电阻率和层厚横向变化应该不太大，特别是在一些沉积地区。当目标层电阻率和层厚横向分布比较连续时，我们可采取横向约束反演方法。横向约束反演方法将测线数据作总体反演，并通过施加横向约束保证反演断面的连续性。它是一种拟二维反演方法。

Auken等^[16]首先将横向约束反演方法引入到直流电法数据反演。Auken和Christiansen^[17]、Santos^[18]和Auken等^[19]分别将横向约束方法运用到二维电阻率数据、EM34剖面数据和三维瞬变电磁数据的反演。Simon等^[20]和蔡晶等^[21]将横向约束反演应用到频率域航空电磁数据中，取得了不错的效果。Vallée和Smith首次使用横向约束反演处理时间域航空电磁数据^[22]。然而，他们仅考虑了简单的两层模型。笔者将加权横向约束反演(weighted laterally-constrained inversion，WLCI)方法^[21]应用到时间域航空电磁数据反演中。我们首先介绍时间域航空电磁法一维正演模拟方法，进而详细阐述时间域航空电磁数据加权横向约束反演理论，结合理论模型和实测数据反演验证加权横向约束反演方法在时间域航空电磁数据处理中的有效性，并阐述加权因子的选取对反演结果的影响。 1 时间域航空电磁法一维正演理论

1.1 垂直磁偶极子的电磁响应公式

定义直角坐标系的原点位于磁偶极子源正下方的地表，z轴向下方向为正。在一维层状各向同性介质中，高度为h的垂直磁偶极子源产生的磁场z分量表达式^[23]为

点击浏览公式

点击浏览公式

1.2 频率域电磁响应到时间域的转换

通过傅里叶变换，我们可以将频率域电磁响应转换到时间域，即

点击浏览公式

点击浏览公式

我们利用Yin等^[24]给出的160点汉克尔滤波系数计算式(4)中的半整数阶汉克尔积分。

对于常规航空电磁系统发射波形，航空电磁响应可通过褶积方法计算。根据褶积定理及阶跃响应与脉冲响应之间的关系^[23]，我们得到

点击浏览公式

2 WLCI反演理论及其在时间域航空电磁数据处理中的应用 2.1 反演数据和参数

时间域航空电磁观测数据为各个时间道的磁感应强度或者磁感应强度随时间的变化率(感应电动势)，以矢量形式表示为

点击浏览公式

点击浏览公式

为了 保证反演参数即各层的电阻率和层厚度为正，对模型参数取对数。将第i个测点的模型参数写成矢量形式：

点击浏览公式

点击浏览公式

灵敏度矩阵又称雅可比矩阵，它由电磁响应数据对模型参数的偏导数组成，表示模型参数对正演响应结果的影响程度。对于时间域航空电磁系统，沿测线方向第k(k=1,2,…,N _s)个测点的灵敏度矩阵可以表示为

点击浏览公式

为方便，我们将正演方程(1) 、(4) 和(5)简化为

点击浏览公式

时间域航空电磁响应与地电参数(电阻率和层厚度)之间存在复杂的非线性关系，为使计算简便，在初始模型 m ₀处对F 做泰勒展开并忽略高阶项，得到

点击浏览公式

点击浏览公式

2.4 电阻率和层厚度横向约束

如前所述，单点反演方法有时会出现反演界面不连续、相邻测点反演参数横向变化剧烈的现象；为此，我们采用拟二维的横向约束反演方法。横向约束反演把相邻测点地电参数的差异作为约束项加入目标函数，以达到光滑目的。为此，我们设

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

在实际反演中，可以通过加权来调整对各层电阻率和层厚度光滑度的横向约束强度。为此，将式(17)改写成

点击浏览公式

点击浏览公式

2.5 深度横向约束

为了得到横向光滑的反演模型，在对电阻率和层厚进行约束的基础上有时还可以加入深度约束。为此，我们设

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

集成数据拟合方程、电阻率和层厚度横向约束方程及深度横向约束方程，得到考虑横向约束的总体反演方程：

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

反演中首先选定一个初始模型，利用方程(37)进行迭代更新，直到拟合误差不再减小或者拟合误差小于一个设定的阀值。

3 反演算例

3.1 理论模型反演

为了验证加权横向约束反演方法在时间域航空电磁数据处理中的有效性，首先对理论模型正演模拟数据进行横向约束反演，并与传统的单点最小二乘反演结果进行对比。航空电磁系统为CGG公司HELITEM MULTIPULSE系统^[27]。发射线圈半径长17.5 m，匝数为1，发射电流为一个4 ms宽的半正弦波和1 ms宽的梯形波，发射基频为30 Hz，接收线圈位于发射线圈侧上方0.5 m，等效面积为1 m²，接收时间道为56道。

设计的理论模型如图 1a所示。表层为电阻率10 Ω·m、厚度20 m的低阻盖层，底层是电阻率为20 Ω·m的黏土层，第二层山谷状结构是电阻率为100 Ω·m的砂砾层。从水文地质的角度考虑，砂砾层很重要，因为它可能是一个保护很好的含水层。采样点间隔为5 m，共101个测点，飞行高度如图 1b所示。每个测点均使用一维模型正演得到dB_z/dt响应数据。为了模拟实测数据，我们对飞行高度和每个测点每个时间道的dB_z/dt响应加了3%的高斯白噪声。

a.理论模型参数； b.飞行高度。 图 1 理论模型 Fig. 1 Synthetic model

图选项

图 2为理论模型单点最小二乘和WLCI反演结果。反演初始模型均选择电阻率为100 Ω·m的均匀半空间(假设三层模型具有相同电阻率，第一层和第二层的层厚为50 m)。实际反演中对电阻率、层厚和深度(3N-2)个参数采用不同的加权因子，研究加权因子对地电参数横向连续性的影响特征。图 2b和图 2c的加权因子分别为(1,1,1,1,0,0,1)，(2.5,2.5,2.5,2.5,0.0,0.0,2.5)。

a.单点最小二乘反演结果；b.WLCI反演结果（加权因子为1,1,1,1,0,0,1）；c. WLCI反演结果（加权因子为2.5,2.5,2.5,2.5,0.0,0.0,2.5）。 图 2 理论模型反演结果 Fig. 2 Inversion results for synthetic model

图选项

从图 2可以看出：1) 单点最小二乘反演和WLCI反演均能反映地下三层地电结构。然而，单点最小二乘反演结果中，层界面比较粗糙，各个测点电阻率的横向连续性差，电阻率呈条带状分布。相比之下，由于施加了横向约束，WLCI方法得到的电阻率横向分布均匀，层界面光滑连续。2) 图 2b和图 2c由于采取了不同的横向约束加权因子，反演层界面的光滑程度和电阻率的连续性存在差异。图 2b由于采用的加权因子太小，在一定程度上改善了反演结果，但电阻率条带仍较明显。相比之下，图 2c由于采用了合适的加权因子，反演结果的电性分界面连续，非常接近真实地电模型。由此我们得出结论：加权因子的选择对改善反演模型的光滑度和连续性非常重要。对于横向连续性较好的地层，可以选取较大加权因子以改善反演断面的连续性。

3.2 实测数据反演

为了进一步验证横向约束反演处理时间域航空电磁数据的有效性，我们对CGG公司使用HELITEM MULTIPULSE系统在加拿大阿尔伯塔省Fort Mcmurray市北西60 km某测区采集的实测数据进行反演。选取10320测线496750——498000段共100个测点的数据 (点距约12 m) 进行WLCI反演。

图 3为实测数据单点最小二乘及WLCI反演结果。初始模型选择电阻率均为100 Ω·m的五层模型，第一层到第四层的层厚均为50 m。图 3a为单点最小二乘反演结果；图 3b的电阻率、层厚和深度的加权因子都为0 (没有横向约束)；图 3c的电阻率、层厚和深度的加权因子都为0.25；图 3d的电阻率、层厚和深度的加权因子都为1.5。

a.单点最小二乘反演结果；b.WLCI反演结果（加权因子均为0）；c.WLCI反演结果（加权因子均为0.25）；d.WLCI反演结果（加权因子均为1.5）。 图 3 实测数据反演结果 Fig. 3 Inversion results for survey data

图选项

从图 3可以看出：1)单点最小二乘和WLCI反演都反映出地下五层电阻率结构。表层电阻率为10~60 Ω·m，层厚为15~40 m；第二层的电阻率为20~50 Ω·m，厚度为30~50 m；第三层的电阻率为2~5 Ω·m，厚度为35~50 m，第四层的电阻率为10~20 Ω·m，厚度为35~50 m；基底的电阻率大于50 Ω·m。由文献^[27]给出的地质资料可知，测区主要是Athabasca砂岩层，且成层性较好。表层是约35 m厚相对良导的冰渍层，内嵌导电性更好的黏土、淤泥薄层。冰渍层下是35~40 m厚的相对高阻的Grand Rapids砂岩。砂岩层下是60~80 m厚的页岩，是测区主要的良导体。图 3反演结果中第三层和第四层即为低阻页岩层。在页岩层之下是存储着油砂矿的Fort Mcmurray层，深处依次是Devonian石灰石和盐层。2)图 3a中单点最小二乘反演结果第二层和第四层电阻率横向不连续，呈明显条带状。3)图 3b各层电阻率、层厚和深度加权因子均为0(没有横向约束)的反演结果和单点最小二乘反演结果相似却不完全一样。这是由于单点最小二乘反演允许各个测点迭代次数不相同，而WLCI反演各个测点迭代次数相同所致。4)由于施加了横向约束，明显改善了图 3c和图 3d反演结果的横向连续性，反演的界面更加光滑，电阻率分布比较连续。5)由于加权因子较小，图 3c中第四层在497600——498000段电阻率和层厚出现了不连续；相比之下，图 3d由于施加较强约束，反演的电阻率和层厚横向连续性均较好，真实地反映了地下电性分布，较好地吻合该测区地质信息^[27]。由此得出结论：在横向约束条件下，正确的地电参数会在不同测点间传播。对于电性结构已知的测点，可以结合已知地质信息(比如钻孔资料)给定正确的初始地电模型，在反演中保持电阻率、层厚和深度不变，从而通过施加横向约束改善整个反演断面的连续性和光滑性。

4 结论

1)与单点最小二乘反演相比，横向约束反演结果中相邻测点的电阻率和层厚光滑连续，电阻率虚假异常得到有效压制，层界面清晰。

2)通过引入参数加权实现对不同参数的光滑度和连续性的制约，可进一步提高反演解的稳定性。

3)加权因子对反演结果有较大影响，加权因子选择太小，横向约束强度低，达不到光滑效果；反之，加权因子太大，模型参数被过度光滑，重要地质特征被平均和圆滑，严重时造成整体反演失败。在实际应用中，应该综合考虑地质资料选择合适的加权因子。

4)理论和实测数据的反演结果表明，拟二维加权横向约束反演是处理时间域航空电磁数据的有效方法。

然而，必须指出的是，对于构造发育地区，反演区段选择应考虑到构造运动造成的地下电性的不连续性。同时，本文中也没有详细研究加权因子的自适应选取等问题，这将是我们未来的研究重点。

CGG公司、Greg Hodges和Tianyou Chen博士为我们提供了HELITEM MULTIPULSE的野外实测数据；同时，殷长春“千人计划”研究团队的其他成员也给予了帮助，在此一并表示感谢。