王琪琪，马 涛，王红兵，谭 鑫，高 鑫，汪强强

(中铁第一勘察设计院集团有限公司 陕西铁道工程勘察有限公司，陕西 西安 710000)

1 引 言

随着大地电磁数据处理解释技术的迅速发展，MT(Magnetotelluric)的应用领域已经从简单的地质构造向复杂的地质构造发展[1-3]。但在这些区域构造离地表很近的地方往往存在着三维局部电性不均匀体，当进行大地电磁法测深，局部不均匀体的边界上会形成积累电荷，使周围的电场突然变大或者减小，从而引起静态效应[4-6]。静态效应往往会使视电阻率拟断面图出现陡立密集的等值线异常，很容易将这种异常解释为地下存在直立的大断裂或纵向延伸的岩脉，同时静态效应也会严重影响反演后的定量解释，如果直接用发生静位移后的视电阻率曲线进行反演，会得到错误的岩层电阻率值和界面深度信息[7]。所以正确认识分析静态效应的响应特征曲线，选择合适的方法识别并进行校正，对大地电磁测深数据的处理与反演解释至关重要[8,9]。

近年来，大地电磁法阻抗张量分解研究工作不断深入，得到了快速发展。赵国泽等[10]将Bahr分解应用到实测电磁数据的处理中，对区域异常分布做了一定的说明；王书明[11]首次将广义逆方法引入到GB分解，并用局部畸变模型进行验证。尹兵祥等[12]提出一种在基于三维构造的正则分解方法，对其参数的物理意义和地质含义进行了解释；杨长福等[13]用视电阻率和相位代替了阻抗张量去实现GB分解，使方程和未知参数的个数减少，提高了计算效率。晋光文等[14]利用Swift和Bahr分解方法对川西-藏东地区的大地电磁阻抗数据进行了校正，分离了区域阻抗和畸变因子，并分析了该地区的构造分布；蔡军涛等[15]通过复数定义提出一种“共轭阻抗变换”分解方法，分解后的阻抗具有和相位张量相似的性质，并得到了区域主轴方位角和畸变因子，同时根据转换后的阻抗定义了新的维性判别参数。谢成良等[16]用相位张量得到电性主轴和维性信息，建立合理的GB分解初始模型，同时将各向异性参数作为未知参数进行计算。

随着大地电磁法勘探规模和计算机处理能力的增大，多测点同步观测将会成为一种趋势，而已有的多频多测点阻抗分解程序需要计算的参数会随测点和频点数呈线性增长。针对传统局部优化算法，如果未知参数没有初值约束，视电阻率的计算结果可能会严重偏离真实值，陷入局部最优，而利用多个测点间的联系来消除静态效应成为了可能[17,18]。

本文利用静态因子在区域分布的统计学性质实现了1D静态校正,并编写了基于粒子群算法和高斯牛顿法的畸变校正程序，通过正演模型验证了该算法的有效性，发现该方法适用于存在多个浅部异常体的区域1D电磁数据。

2 基本原理

2.1 G B分解

1989年，Groom和Bailey[19]对阻抗张量进行了如下分解：

Zm=RCZRT

(1)

其中，R是旋转矩阵;Z是区域未受畸变的一维阻抗张量;C是一个与频率无关的实数畸变矩阵，C=gTSA；g是增益因子，为无量纲参数；T是扭曲矩阵；S是剪切矩阵；A是各向异性因子矩阵，它们的矩阵分别表示如下：

式(2)～式(6)中，θ是主轴方向角，单位为度；Z1为一维区域阻抗张量元素，单位为Ω；t是扭曲因子；e是剪切因子；s是各向异性因子；t,e,s均为无量纲参数。其中t、e组成了畸变矩阵C的可确定部分，g、s组成了不可确定部分，所以实测阻抗又可写为下式：

(7)

2.2 畸变校正

李洋[20]提到对于局部优化算法，初值的确定有很大的影响，本文先通过粒子群全局优化算法得到一个包含畸变参数和区域阻抗的全局最优解空间，然后将全局最优解空间作为高斯牛顿法的初值进行计算。该方法能够更好地避免陷入局部最优，得到更真实的畸变参数，使得校正后区域阻抗更可靠。算法采用的目标函数如下所示：

利用式(7)得到如下参数:

其中,σ=A+B,δ=A-B。建立如下目标函数:

(12)

2.2.1 粒子群算法

粒子群算法是一种全局随机搜索算法[21]。该算法模拟鸟群觅食行为，将鸟看做粒子，通过粒子个体之间共享信息来找寻最优解。每个粒子有速度和位置两个属性，粒子通过自身当前在搜索路径中的最优值和整个粒子群中的最优值来调整速度和位置，其更新方式如下所示：

其中Vi、xi为第i个粒子的速度和位置，ω为惯性权重，c1、c2为学习因子，控制粒子的认知模式，ξ和η为0～1间的随机数，pbest、gbest分别为粒子个体最优值和全局最优值，r为约束因子。当目标函数式满足规定的阈值或者达到最大迭代次数时迭代结束。

2.2.2 高斯牛顿法

目标函数公式(12)二阶泰勒展开如下所示：

(15)

其中,a′是每次迭代得到的畸变参数和区域阻抗;c为a′对应的目标函数值;H为海森矩阵。则γ2(a)梯度为：

(16)

(17)

因为在H(a(k))矩阵中二阶项很小，则H(a(k))≈JTJ，代入上式得：

(18)

当步长足够小或者迭代次数达到阈值时停止迭代，此时a(k+1)为最优解。

2.3 增益因子

上述步骤对GB分解中确定参数进行了拟合，可以得到仅受静态效应影响的阻抗数据，本节主要对不确定参数进行校正，具体步骤如下：

当某一测区的测点数足够多时，地表不均匀体产生的电流畸变在统计计算上将会被光滑[22]，利用这一假设，得到多测点多频点下的SSQ旋转不变量：

(19)

其中,ri为第i个测点;ω为角频率，单位为rad/s。区域构造为1D时，整个测区的SSQ旋转不变量为：

(20)

对应的每个测点视增益表示为：

(18)

gi为第i个测点对应的视增益因子。

3 正演模型分析

建立一个如图1所示的二层模型,第一层厚度为2 km，电阻率为100 Ω·m，第二层电阻率为1 000 Ω·m，浅部存在三个低阻异常体，电阻率为1 Ω·m，x和y方向延伸400 m，z方向延伸500 m。

图1 正演模型

由图2可以看出，当迭代100次之前，畸变参数的拟合误差发生快速下降，之后随着迭代次数的增多误差不再发生变化，约为3.8×10-4，表明该算法收敛速度较快，拟合结果比较精确。

图2 高斯牛顿法迭代过程中均方根误差曲线

图3和图4显示了频率分别为350 Hz和0.07 Hz时不同测点位置上畸变参数和静态因子的横向分布，图中冒号后面的数字表示不同畸变参数随测点变化的标准差。可以发现整体上增益因子g和各向异性因子s随测点的变化比较大，g存在两个极大值和三个极小值，极大值对应异常体间的测点，极小值对应三个异常体正上方，s的分布与之相反。扭曲因子t和剪切因子e随测点变化较小，标准差都在0.005以下，其中频率为320 Hz时静态因子横向分布的标准差小于0.07 Hz，表明当频率较高时,地表异常体对视电阻率的影响不再是静态畸变，而是作为一个局部异常体。

图3 频率为320 Hz时畸变参数和静态因子分布

图4 频率为0.07 Hz时畸变参数和静态因子分布

图5为27号测点校正前后的视电阻率和相位曲线，测点位于异常体一侧，其中xy_correct、yx_correct分别表示校正后的电阻率和相位，xy_distort、yx_distort分别表示存在浅部异常体时的电阻率和相位，xy_nodistort、yx_nodistort分别表示不存在浅部异常体时的电阻率和相位。由图5(a)和图5(b)可以发现,地表异常体的存在使得xy模式和yx模式视电阻率发生了明显分离，xy模式下移，yx模式上移，相位曲线在高频段出现了分离，受到明显静态效应影响，校正后的视电阻率与未含异常体的视电阻率几乎重合，原来分离的视电阻率发生闭合；图5(c)和图5(d)中10 Hz处静态因子存在一个拐点p，与图5(a)和图5(b)对照发现，频率高于10 Hz时，视电阻率和相位曲线包含的是浅部异常体的信息，低于10 Hz时，区域构造受异常体影响曲线发生偏移，此时图5(c)中拐点p所指向的位置，可以作为判断静态效应影响频率范围的一个依据。

图5 27号测点和畸变参数随频率变化曲线

图6为31号测点校正前后的视电阻率和相位曲线，测点位于异常体正上方。xy模式和yx模式受静态效应的影响发生向上偏移，和27号测点一样，校正后的xy模式和yx模式视电阻率曲线相互重合，表现出明显的1D性质，说明该方法很好地消除了畸变的影响。

图6 31号测点校正前后曲线和畸变参数随频率变化曲线

4 结 论

本文针对部分1D区域构造，引入SSQ旋转不变量来计算得到增益因子，结合GB分解实现了区域阻抗的恢复。同时以正演模拟数据为例，获得了能够反映真实地下介质分布的视电阻率曲线和相位以及相关畸变参数，方便后续反演工作的进行，验证了该方法的有效性。最后针对论文内容总结了以下几点：

1)在多测点多频点GB分解的基础上，将粒子群算法和高斯牛顿法结合，避免计算结果陷入局部最优。

2)通过建立的畸变模型发现，相比视电阻率相位受畸变影响更小，并且影响大小随频率发生变化，所以相位可以作为判断区域维性的重要参数。同时针对均匀浅部异常体，不确定畸变参数(增益因子和各向异性因子)呈“正态分布”，在异常体中心最大，逐渐向两侧减小，同时发现不确定参数存在一个拐点，可以用来确定静态效应的影响频段。

3)通过正演模拟数据表明，本文提出的校正方法针对于1D区域构造能够很好地压制多个地表异常体的畸变影响，校正后的视电阻率和相位曲线几乎重合，能够更真实地反映地下介质分布情况。