0 引言

瞬变电磁法(Transient Electromagnetic Me-thod，TEM)是一种重要的可控源电法勘探方法，该方法采用不接地回线或接地导线向地下发射脉冲电流，在一次场的间歇期，观测大地中的感应电磁场(通常观测垂直方向的感应电动势)，也称为二次场。通过研究二次场随时间和空间的变化规律，分析大地的电性分布特征。因采用人工源，并且观测二次场，该方法具有数据精度高、不受一次场影响等优点，已被广泛地应用于工程、水文、环境、资源等勘探领域^[1-5]。

尽管当前TEM数据的二维和三维正反演解释取得了巨大进步^[6-11]，然而由于反演计算量大，需要大规模的计算平台，难以普遍应用于实际资料^[12-16]，故一维反演仍是目前实际资料处理解释主要手段之一^[17-21]。然而，TEM数据的一维反演存在反演剖面连续性差的问题，且数据记录周期的晚期信号弱，单点反演的深部数据受噪声影响大。

基于大地介质电性连续变化的特点，Auken等^[22]首先将横向约束反演方法用于直流电电磁数据反演；Vallée等^[23]将横向约束方法应用于时间域航空电磁数据的反演，有效改善了反演效果；蔡晶等^[15]、殷长春等^[16]分别将该方法用于航空电磁数据频域和时域反演，明显改善了整个反演剖面的连续性，提高了解释结果的准确性。

传统的横向约束反演方法只沿测线方向进行横向约束。对于三维测网的观测数据，若只沿测线方向进行横向约束而不考虑相邻测线之间构造的关联性，会大大降低地质解释的可靠性。为此，Viezzoli等^[24]提出了基于空间约束的航空电磁数据反演方法，即不仅沿测线方向，也垂直于测线方向施加横向约束。殷长春等^[25]将空间约束用于航空电磁拟三维反演，通过理论和实测数据测试证明了该方法适合求解多测线航空电磁数据的最优化反演。

本文从横向约束原理出发，将横向约束引入地面三维测网的TEM数据的一维反演解释，给出了一种距离加权进行相邻点横向约束的拟三维反演方法。该方法能有效地实现横向约束，保证相邻测点模型的连续性，且不要求测网是严格规则的，实现了一般测网条件下的横向约束拟三维反演。模型正演数据和实际数据测试结果验证了该方法的正确性、有效性和实用性。

1 横向约束拟三维反演方法

传统的横向约束拟二维反演方法首先集成剖面数据，再沿测线施加横向约束，同时反演测线上多测点的地电参数^[16]。该方法的实质是把相邻测点地电模型的差异作为约束项加入目标函数，对模型进行横向平滑处理，以保证相邻测点间模型的连续性，进而压制因为噪声导致的单个测点地电参数突变。TEM勘探中，测点设计常采用测网方式，且通常测点分布是非完全规则的。基于横向约束的原理，本文采用一种距离加权相邻点进行横向约束的拟三维反演方法，即：对所有测点数据进行集成，对所有测点均采用半径r之内、距其最近的p个测点进行距离加权横向约束，反演测点处的地电模型，实现一种广泛适用的TEM数据横向约束拟三维反演。

距离加权相邻点横向约束的原理如图 1所示。图中测点按顺序编号，各测点的反演由距其最近的p个点进行横向约束。例如，33号测点由距其最近的24、34、25、32号测点(此时p=4)进行横向约束，这4个测点分别称为33号点的第1~4号相邻点。

设三维测网有l个测点，按一定顺序进行编号，并设每个测点有k个衰减延时(实际各个测点的延时数可以不一致，此处为了表述方便，都设为k)。这l个测点的电磁数据可表示为

$ \begin{array}{l} \mathit{\boldsymbol{d}} = [{d_{11}},{d_{12}}, \cdots ,{d_{1k}}, \cdots ,{d_{i1}},{d_{i2}}, \cdots ,\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {d_{ik}}, \cdots ,{d_{l1}},{d_{l2}}, \cdots ,{d_{lk}}{]^{\rm{T}}} \end{array} $

(1)

假设反演模型为n层，只考虑各层的电导率σ，则l个测点的总模型可表示为

$ \begin{array}{l} \mathit{\boldsymbol{m}} = [{\sigma _{11}},{\sigma _{12}}, \cdots ,{\sigma _{1n}}, \cdots ,{\sigma _{i1}},{\sigma _{i2}}, \cdots ,\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\sigma _{in}}, \cdots ,{\sigma _{l1}},{\sigma _{l2}}, \cdots ,{\sigma _{\ln }}{]^{\rm{T}}} \end{array} $

(2)

以各测点为中心，在横向约束半径r内搜索距其最近的p个测点(如果数据点不足，允许少于p个测点)，并分别记录这p个点的序号及其距离。

基于吉洪诺夫正则化反演思想，反演目标函数中包含数据拟合项及正则化项。将横向约束作为正则化项加入反演目标函数，依次将各测点与其第1~p个相邻点的模型差异加入目标函数，再将各测点模型的纵向粗糙度加入目标函数，由此得到横向及纵向上均施加了平滑约束的反演目标函数。求解该反演目标函数，即可得到拟三维反演结果。距离加权横向约束的拟三维反演目标函数可写为

$ \begin{array}{l} \phi = {[F(\mathit{\boldsymbol{m}}) - \mathit{\boldsymbol{d}}]^{\rm{T}}}{\mathit{\boldsymbol{C}}^{ - 1}}[F(\mathit{\boldsymbol{m}}) - \mathit{\boldsymbol{d}}] + {k_{\rm{v}}}{({\mathit{\boldsymbol{E}}_{\rm{v}}}m)^{\rm{T}}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({\mathit{\boldsymbol{E}}_{\rm{v}}}\mathit{\boldsymbol{m}}) + {k_{\rm{h}}}\sum\limits_{i = 1}^p {{{({\mathit{\boldsymbol{E}}_{{\rm{h}},i}}\mathit{\boldsymbol{m}})}^{\rm{T}}}} {\mathit{\boldsymbol{D}}_i}({\mathit{\boldsymbol{E}}_{{\rm{h}},i}}\mathit{\boldsymbol{m}}) \end{array} $

(3)

式中：F为正演算子；C为协方差矩阵，用于基于数据误差对观测数据拟合进行加权；k_v、k_h分别为纵向和横向约束因子；E_v为各测点模型的纵向一阶差分矩阵；E_{h, i}为各测点第i相邻点横向约束一阶差分矩阵(若某测点的第i相邻点不存在时，E_{h, i}中对应于该测点的行元素全为零)；D_i为第i相邻点横向约束距率加权矩阵；E_vm和E_{h, i}m分别表示所有测点对应模型的纵向粗糙度及其第i个相邻点的横向粗糙度。E_v、E_{h, i}及D_i可分别表示为

$ \begin{array}{l} {\mathit{\boldsymbol{E}}_{\rm{v}}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left[ {\begin{array}{*{20}{c}} { - 1}&1&0& \cdots &{}&{}&0\\ 0&{ - 1}&1&0& \cdots &{}&0\\ {}&{}&{}& \vdots &{}&{}&{}\\ 0& \cdots &{}& \cdots &0&{ - 1}&1 \end{array}} \right]_{[l \times (n - 1)] \times (l \times n)}}\\ \end{array} $

(4)

$ \begin{array}{l} {\mathit{\boldsymbol{E}}_{{\rm{h}},i}} = {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left[ {\begin{array}{*{20}{c}} { - 1}&0& \cdots &1&0& \cdots &0\\ 0&{ - 1}&0& \cdots &1& \cdots &0\\ {}&{}&{}& \vdots &{}&{}&{}\\ 0& \cdots &0&1&0& \cdots &{ - 1} \end{array}} \right]_{(l \times n) \times (l \times n)}} \end{array} $

(5)

$ {\mathit{\boldsymbol{D}}_i} = {\left[ {\begin{array}{*{20}{c}} {\frac{1}{{{r_{1,i}}}}}&0& \cdots &{}&{}&0\\ 0&{\frac{1}{{{r_{1,i}}}}}&0& \cdots &{}&0\\ {}&{}& \vdots &{}&{}&{}\\ {}&{}&{}&0&{\frac{1}{{{r_{l,i}}}}}&0\\ 0& \cdots &{}&{}&0&{\frac{1}{{{r_{l,i}}}}} \end{array}} \right]_{(l \times n) \times (l \times n)}} $

(6)

式中r_{l, i}表示第l个测点与其第i个相邻点的距离。在实际数据处理过程中，需要设定一个最小距离r₀，当r_{l, i}<r₀时，用r₀代替。式(3)中依次加入了各测点与其相邻p个测点地电模型的差异对模型进行横向约束，求该目标函数的极小，即能得到横向平滑约束的最小二乘解，从而实现了横向约束拟三维反演。

式(3)中的纵向约束及横向约束项为反演目标函数的平滑正则化项，根据吉洪诺夫正则化的思想，纵向约束因子k_v及横向约束因子k_h可依据工区内地层的纵向和横向变化特征初选，在保证数据拟合的前提下，以后验方式选择合理的值。

由于各测点都选取相邻测点进行横向约束，该方法不要求测网严格规则，同时对各测点的排序也不做要求，具有广泛的适应性。此外，该方法同样适用于二维数据的反演，此时该方法即为常规的横向约束拟二维反演(LCI)。

2 横向约束拟三维反演目标函数的求解

高斯—牛顿法具有迭代稳定、收敛快的特点，因此本文采用高斯—牛顿迭代法求解式(3)中横向约束拟三维反演目标函数的极小。

设m₀为模型初值，Δm为模型修正量，将F在m₀处一阶展开，目标函数对Δm求导，并令导数为零，得到高斯—牛顿法模型修正量的求解式为

$ \begin{array}{*{20}{l}} {\left[ {{\mathit{\boldsymbol{J}}^{\rm{T}}}{\mathit{\boldsymbol{C}}^{ - 1}}\mathit{\boldsymbol{J}} + {k_{\rm{v}}}\mathit{\boldsymbol{E}}_{\rm{v}}^{\rm{T}}{\mathit{\boldsymbol{E}}_{\rm{v}}} + {k_{\rm{h}}}\sum\limits_{i = 1}^p {\left( {\mathit{\boldsymbol{E}}_{{\rm{h}},i}^{\rm{T}}{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{E}}_{{\rm{h}},i}}} \right)} } \right]\Delta \mathit{\boldsymbol{m}}}\\ {{\kern 1pt} = {\mathit{\boldsymbol{J}}^{\rm{T}}}{\mathit{\boldsymbol{C}}^{ - 1}}\left[ {\mathit{\boldsymbol{d}} - F\left( {{\mathit{\boldsymbol{m}}_0}} \right)} \right] - {k_{\rm{v}}}\mathit{\boldsymbol{E}}_{\rm{v}}^{\rm{T}}{\mathit{\boldsymbol{E}}_{\rm{v}}}{\mathit{\boldsymbol{m}}_0} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {k_{\rm{h}}}\sum\limits_{i = 1}^p {\left( {\mathit{\boldsymbol{E}}_{{\rm{h}},i}^{\rm{T}}{\mathit{\boldsymbol{D}}_i}{\mathit{\boldsymbol{E}}_{{\rm{h}},i}}{\mathit{\boldsymbol{m}}_0}} \right)} } \end{array} $

(7)

式中J为各测点正演函数与相应模型雅可比矩阵的集成，即顺序地以各测点一维雅可比矩阵为块的分块对角矩阵，该矩阵仅包含l×k×n个非零元素。

由式(7)可求解Δm, 得到修正后的模型m=m₀+Δm，再把m赋给m₀，这样反复迭代，直到达到预设的拟合误差或最大迭代次数，即得到最终反演结果。

式(7)反演迭代的关键在于求取正演函数相对于模型的雅可比矩阵。本文重点研究最常用的回线源装置瞬变电磁法。一维正演采用首先在频率域进行、再通过正弦变换转换到时间域的经典算法，雅可比矩阵的求取则采用拟正演的方法。

水平矩形回线源布设于地面，测点可位于地面上除发射导线处的任何位置，记录数据为垂直磁场的感应电动势。对于回线源频率域的响应，可将矩形线圈的四边看作是4个导线源，则回线源的响应可表示为四个导线源的叠加场，地面上的磁场垂直分量^[26]为

$ {H_{\rm{z}}} = \sum\limits_{j = 1}^4 {\frac{I}{{4{\rm{ \mathit{ π} }}}}} \int_{ - {L_j}/2}^{{L_j}/2} {\frac{y}{R}} \int_0^\infty {(1 + {r_{{\rm{TE}}}})} \frac{{{\lambda ^2}}}{{{\mu _0}}}{J_1}(\lambda R){\rm{d}}\lambda {\rm{d}}{x^\prime } $

(8)

式中：J₁(λR)为一阶贝塞尓函数，λ为内侧积分变量，$R=\sqrt{\left(x-x^{\prime}\right)^{2}+y^{2}}$，x为测点到该边中垂线的距离，x′为沿该边进行积分的积分变量，y为测点与该边的距离；L_j为发射回线某边的长度；r_TE为地面TE模式的反射系数；I为发射电流；μ₀为空气中的磁导率。

通过正弦变换

$ {{\dot H}_{\rm{z}}}=\frac{2}{{\rm{ \mathit{ π} }}}\int_0^\infty {{\rm{Im}}[{H_{\rm{z}}}(\omega )]\sin (\omega t){\rm{d}}\omega } $

(9)

可得到H_z磁场时间域的响应$\dot H$_z(磁场的时间变化率)。

根据式(8)和式(9)，只有反射系数r_TE中包含模型参数，因此只需要对r_TE求导。本文借鉴MT Occam反演^[27]中的正演函数对地层电导率的求导方法，基于反射系数递推，采用拟正演方法计算瞬变响应对于各层的偏导数。

^r_TE的表达式^[26]为

$ {r_{{\rm{TE}}}} = \frac{{{Y_0} - {{\hat Y}_1}}}{{{Y_0} + {{\hat Y}_1}}} $

(10)

$ {{\rm{Y}}_0} = \frac{{{{\rm{u}}_0}}}{{{{{\rm{\hat z}}}_0}}} $

(11)

其中

$ {{\hat Y}_1} = {Y_1}\frac{{{{\hat Y}_2} + {Y_1}\tanh ({u_1}{h_1})}}{{{Y_1} + {{\hat Y}_2}\tanh ({u_1}{h_1})}} $

(12)

$ {{\hat Y}_i} = {Y_i}\frac{{{{\hat Y}_{i + 1}} + {Y_i}\tanh ({u_i}{h_i})}}{{{Y_i} + {{\hat Y}_{i + 1}}\tanh ({u_i}{h_i})}} $

(13)

$ {{\hat Y}_N} = {Y_N} $

(14)

$ {Y_i} = \frac{{{u_i}}}{{{{\hat z}_i}}} $

(15)

式中：$\hat{z}_{i}$=iμ_iω，μ_i为各地层的磁导率，本文中都取为μ₀(空气磁导率)；${u_i} = {\left( {{\lambda ^2} - k_i^2} \right)^{\frac{1}{2}}}$，k_i²=-iωμ_iσ_i，特别地，u₀表示空气中的参数，其值可根据空气中的磁导率μ₀和电导率σ₀计算得到；N表示大地模型的层数；h_i表示各地层厚度。反射系数r_TE对各层电导率σ_i的导数，可利用$\hat{Y}_{1}$的递推快速计算得到

$ \frac{{\partial {r_{{\rm{TE}}}}}}{{\partial {\rm{lg}}{\sigma _i}}} = \ln 10 \times {\sigma _i}\frac{{\partial {r_{{\rm{TE}}}}}}{{\partial {\sigma _i}}} = \frac{{ - 2\frac{{\partial {{\hat Y}_1}}}{{\partial {\sigma _i}}}{Y_0}}}{{{{({Y_0} + {{\hat Y}_1})}^2}}} $

(16)

$\hat{Y}_{1}$由$\hat{Y}_{\rm N}$递推而来，因此$\frac{\partial \hat{Y}_{1}}{\partial \sigma_{i}}$可展开为

$ \begin{array}{l} \frac{{\partial {{\hat Y}_1}}}{{\partial {\sigma _i}}} = \frac{{\partial {{\hat Y}_1}}}{{\partial {{\hat Y}_2}}} \times \frac{{\partial {{\hat Y}_2}}}{{\partial {{\hat Y}_3}}} \times \cdots \times \\ \quad \quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\partial {{\hat Y}_i}}}{{\partial {{\hat Y}_{i + 1}}}} \times \cdots \times \frac{{\partial {{\hat Y}_i}}}{{\partial {\sigma _i}}}\quad j + 1 \le i \end{array} $

(17)

式(16)表明，计算下一层电导率的导数可以利用上一层计算的中间结果，从而减少计算量，提高计算效率。

需要注意的是，计算$\frac{{\partial {{\hat Y}_j}}}{{\partial {{\hat Y}_{j + 1}}}}$涉及到多个ω值及多个λ值，因此需要存储多个递推中间值；由于${{{\hat Y}_j}}$是从最下层往上递推，因此需要预先递推计算。

3 模型试算

为验证上述横向约束拟三维反演方法的效果，设计包含一个高阻体及一个低阻体的模型(图 2a和图 2b)进行正演模拟，再用正演合成数据进行反演测试。模拟采用中心回线装置，设发射回线为100m×100m的正方形，电流为10A，接收线圈面积为200m²，观测时间范围为0.00001~0.00500s，测网点线距均为100m，共计21×21个测点。正演场中加入典型的1 nV/m²随机背景噪声。由于本实验的目的是验证横向约束拟三维反演的有效性，因此对各测点对应的模型采用一维正演。

图 2c为过高阻体及低阻体中心测线(y=950m)的加噪感应电动势多测道曲线。由于TEM信号早期强、晚期弱，所以图 2c中早期段中没有明显噪声，而晚期段中噪声明显。

采用本文横向约束拟三维反演方法对正演数据(图 2c)进行反演，同时进行单点一维反演，以便对比分析两种方法的效果。两种方法反演的拟合误差均设为小于5%。反演结果见图 3。从图 3a可以看到，单点一维反演结果中，高阻体及低阻体的形态大致能够恢复，但高阻体内存在明显的横向不连续特征；反演剖面的下部则出现明显不合理的横向不连续现象及电阻率突变，其原因是晚期噪声的影响。从图 3b可以看出，横向约束拟三维反演对模型的恢复效果很好，低阻体及高阻体形态基本完整，边界清晰；异常体以下的区域形态也恢复较好。可见横向约束拟三维反演通过增强横向连续性，能够有效地压制单点反演造成的不合理模型突变。此外，在横向上，拟三维反演的地层界面很清晰，并没有因为横向约束导致界面模糊，原因在于地层界面产生的瞬变异常(图 2c)边界清晰，反演函数中数据拟合项的影响强于横向约束项。

4 实际数据应用

在中东某油田的采油区，地震资料显示在浅部(深度约200m)可能存在凹陷，但不能完全确定。这类浅部凹陷严重影响地震资料的深部成像质量。为查清该浅部结构，改善地震剖面成像效果，在该区部署了TEM勘探，调查400m深度以上地层结构。

TEM采用200m×200m的中心回线装置，观测数据为垂直磁场分量。设计线距为400m，点距为250m，但由于工区位于采油区，油田设施众多，为避开油田设施的干扰，实际测点分布很不规则(图 4)。

由于测点分布很不规则，因此采用横向约束拟三维反演方法。为了对比分析反演方法效果，对其中一条测线(N1测线)同时进行常规单点一维反演。图 5为N1测线经面积、电流归一化后的观测感应电动势等值线剖面图，图 6为该测线单点一维反演和横向约束拟三维反演电阻率剖面及其对应的归一化感应电动势正演剖面。N2、N3测线反演电阻率剖面见图 7，图 8为全工区拟三维电阻率反演结果在不同海拔上的水平切片。

可以看出，单点一维电阻率反演剖面能大致体现沿测线的构造特征，但各测点间的连续性不好，横向上孤立异常突出，地层电性特征显示也不甚清楚；而横向约束拟三维反演的电性层横向连续性较好，地层电性特征清楚，很大程度上避免了单点反演中的不合理电性突变，有效提高了反演剖面的精度。还可以看出，相邻的N1、N2及N3测线反演剖面上电性层横向连续，电性特征清楚，电性层在各剖面上能够很好地进行追踪，中部的目标凹陷(横坐标2~4km)被清晰地勾绘了出来。图 7中反演电阻率剖面叠加了后期微测井视电阻率曲线，可见反演剖面与电测井曲线吻合较好，进一步说明本文横向约束拟三维反演结果的可靠性。图 8中，地层横向变化特征清楚，清晰地显示了中部浅层凹陷的形态及走向。

实际数据应用结果说明本文方法能够有效地实现一般测网条件下的TEM数据横向约束拟三维电阻率反演。

5 结论

针对TEM数据一维电阻率反演中存在反演剖面连续性差的问题，本文将横向约束引入TEM数据反演中；考虑到实际测点分布并不严格规则的特点，构建了距离加权相邻点横向约束的拟三维反演方法，并采用高斯—牛顿法进行求解。

模型合成数据及实际数据的计算结果说明，距离加权相邻点横向约束拟三维反演可有效实现横向约束，保证相邻测点模型的连续性，提高反演解释效果，且对一般非规则测网TEM数据也是适用的，具有广泛的适用性。