陈 辉，杜兴忠，孙永清

(中国电建集团 贵阳勘测设计研究院有限公司，贵阳 550081)

0 引言

瞬变电磁法(Time domain electromagnetic methods，TEM)是利用线圈或接地电极观测二次场(响应场)，通过对这些响应信息提取和分析，从而达到探测地下地质体的方法[1]。国内、外许多学者对于瞬变电磁法从1维正反演，到2.5维正反演，再到3维正反演做了大量工作，均取得了很多有价值的研究成果，但这些研究大部分是基于低阻异常体的响应特征出发，很少研究纯高阻体地质状况的响应特性。笔者从建立1维高阻层状模型和2.5维高阻矩形柱模型出发，分别研究其正演响应特征，并依据正演响应结果，对反演做了初步分析。

1 TEM 2.5维正反演概述

1.1 正演模拟

2.5维正演模拟先利用拉氏变换，将2.5维电磁问题从时间域转换到拉氏域，再利用傅氏变换从拉氏域转换到傅氏域中，并构建拉氏傅氏域中电磁场分量沿走向方向的对偶二阶微分方程，由对偶二阶微分方程和边界条件形成边值问题；然后，采用有限元法求解该边值问题，得到电磁场的解。将求解的方程结果经过逆拉氏傅氏变换，可以得到各节点在时间域三维空间的电磁场分量。

对于瞬变电磁2.5维的对偶二阶方程的推导，前人已做了大量研究，这里直接引用，中心回线TEM 2.5维沿走向方向的电磁场二阶对偶偏微分方程形式为[2-3]：

a4u+b4up

(1)

a5v

(2)

其中，

(3)

(4)

其中，系数ηu和ηv计算式为：

(5)

(6)

电磁场二阶对偶微分式(1)和式(2)，边界条件式(3)和式(4)四个方程一起构成了在(s,m)域中沿走向方向(这里为y方向)电磁分量的边值问题，采用有限元法进行求解。

(7)

式中：求和下限t=i1是(m+1)/2的整数部分。

采用余弦变换数值滤波算法做逆傅氏变换，根据正余弦变换数值滤波算法[5]，得到式(8)。

(8)

1.2 反演分析

1.2.1 阻尼最小二乘法

瞬变电磁常采用的是阻尼最小二乘法进行反演解释，它是1963年马奎特(Marquardt)[6-7]通过改进高斯-牛顿法的法方程组而得到的一种新的反演方法。对于目标函数式(9)的极小值问题，阻尼最小二乘法的处理方式为式(10)。

(9)

(JTJ+λI)δ=ε

(10)

其中：J表示雅可比矩阵；λ表示阻尼因子(它是一个正数，用以控制校正的方向和步长)；δ是修正量即δ=p(k+1)-p(k)，从式(10)中可以明显看出它是λ的函数即δ=δ(λ)，ε是残差向量。所以阻尼最小二乘法的迭代公式为：

p(k+1)=p(k)+(JTJ+λI)-1ε

(11)

在反演算法中，雅克比矩阵J的计算精度与反演结果的准确性有直接关系，设从模型空间映射到数据空间的函数为φ，模型参数用p表示。对雅克比矩阵行列中的每个元素的一阶偏导数的计算，采用有限差分方法中的一阶向前差分为式(12)。

(12)

阻尼最小二乘反演法的反演是一个循环迭代的过程，迭代一次修改一次模型参数，直到由模型计算出的电磁响应值与实测或者理论数据间的差值满足设定的要求。

1.2.2 阻尼特征参数法

阻尼特征参数反演法(Damped Eigenparameter Inversion approach)，是由Jupp[8]和Vozoff[9]提出并用于地球物理数据联合反演解释中，它结合了马奎特方法和奇异值截断的优点，将经奇异值分解后得到的特征值分为重要参数和不重要参数，截断小的奇异值，通过改变奇异值来修正阻尼因子。相对于阻尼最小二乘方法从两次计算结果来改变阻尼因子，该方法综合所有的数据特性，从整体上指导反演过程的进行，具有快速稳定的特点。

对于时间域瞬变电磁勘探，假设M个测点感应电动势的数据表示成向量为d=(d1,d2,…,dM)T，表示地下介质电导率等电性参数的模型参数可以表示成P=(P1,P2,…,PN)T。

大多数地球物理问题并不是线性的，对于非线性反演问题的线性化用雅可比矩阵进行描述，它是一种从N维参数空间到M维数据空间的线性映射。引入误差向量ε，方程可写成：

(13)

当采用式(13)计算时，如果奇异值Sj较小，计算会出现不稳定甚至不收敛的情况，为解决该问题，可以通过阻尼因子tj改善反演过程，于是引入一个阻尼修正步长：

(14)

每迭代一次就修改一次模型参数，由式(12)可得出，新模型参数的计算式为式(15)。

P(new)=P(old)+VδQ

(15)

新模型参数生成后，重新计算新模型的响应值和误差向量。反演过程就是一个循环迭代的过程，直到满足预先设置的结束准则迭代终止。

1.2.3 反演计算步骤

本文实现反演的计算步骤如下：

1)设置反演参数，读取反演观测数据。

2)设置初始模型，计算其正演理论响应值。

3)计算模型理论响应值和反演观测数据的误差。

4)判断误差是否满足要求，不满足要求跳到步骤5)，满足要求直接跳到步骤9)。

5)计算雅克比矩阵。

6)计算新的模型参数及其理论响应值。

7)计算新模型的理论响应值和反演观测数据的误差。

8)再次判断误差是否满足要求或者达到最大迭代次数，不满足要求或未达到最大迭代次数则跳到步骤5)，满足要求或已达到最大迭代次数则跳到步骤9)。

9)输出得到迭代结束后的模型参数。

2 正演模拟

2.1 高阻层状模型1维TEM响应特征

根据TEM一维正演理论[10-12]，构建一维层状K型地电模型如图1，具体模型参数见表1。

图1 层状模型

表1 层状模型参数

高阻层状模型1～高阻层状模型4和两种均匀半空间模型(电阻率为10 Ω·m和103Ω·m)的瞬变电磁响应曲线共同绘入图2，从图2中可以看到，对于相同电阻率的高阻层：同一层厚时(模型1和模型2、模型3和模型4)，围岩电阻率相对越高，高阻层响应曲线进入衰减的时间相对越早。对于高阻层状模型整体来看，在关断电流后的某个时刻(设这个时刻为t1)后响应曲线衰减速度相对较快，然后在某一时刻(设这个时刻为t2)达到一个最大值，其后衰减速度逐渐减慢。t1、t2值可通过模型1～模型4与对应均匀半空间模型的相对比值来确定(图3)，比值越远离1线表示模型响应值和均匀半空间响应值差异相对越大，反之亦然。图3中比值远离1线的起始位置(可选择与1线的相对误差等于0.01%作为起始位置判定点)对应的时刻即为t1的值，离1线垂直距离最远的点对应的时刻即为t2值。从图2、图3中还可以看到，同一围岩电阻率时(模型1和模型3、模型2和模型4)，高阻层相对越厚，进入衰减过程t1时刻后正演响应随时间衰减相对越快，达到衰减最大t2时刻相对较晚，反之亦然。

图2 高阻层状模型TEM响应曲线

图3 模型和均匀半空间TEM响应的相对比值

层状模型5～层状模型8的高阻层电阻率差异不同倍数时TEM响应见图4、图6、图8、图10。采用相邻倍数响应的比值曲线对比高倍数与前一倍数响应的差异情况，见图5、图7、图9、图11，并计算各模型相邻倍数响应相对误差均值见表2。从图3～图10和表2中可以看到，地层层厚相同时，同一模型在高阻层电阻率超过围岩电阻率103倍后，相邻倍数TEM响应的相对误差均值均小于0.03%，说明超过103倍后，即使高阻层电阻率较围岩电阻率变化很大，但TEM响应值变化却很小(几乎无变化)。

表2 相邻倍数TEM响应相对误差平均值表

图4 电阻率差异不同倍数时TEM响应曲线(模型5)

图5 相邻倍数响应的比值曲线(模型5)

图6 电阻率差异不同倍数时TEM响应曲线(模型6)

图7 相邻倍数响应的比值曲线(模型6)

图8 电阻率差异不同倍数时TEM响应曲线(模型7)

图9 相邻倍数响应的比值曲线(模型7)

图11 相邻倍数响应的比值曲线(模型8)

2.2 高阻矩形柱模型2.5维TEM响应特征

本次研究采用有限元法，将所要求解的连续场离散为很多相互连接的较小单元，把连续的求解域离散化，在每个单元内用插值函数进行插值，建立起每个单元网格各自的公式，最后再将这些单元集合起来求解，以各单元节点函数值为未知量的线性方程组，得出原来的场值问题的解[13]。建立离散网格，向下前5层每层层厚设置为10 m，其后的每层厚度比上层按6%递增，深度范围0 m～400 m；横向采用等间距为20 m的均匀网格，宽度为-650 m～650 m；测点范围为-590 m～590 m，测点间距为20 m。发射电流最大值为10 A，中心回线发射线圈半径为13 m，匝数为4匝；发射电流波形为方波，电流平稳时间为7.6 ms，7.6 ms关断电流；二次场采样时间为9.60 ms，26个时间采样道。均匀半空间中建立矩形柱体模型，模型顶部离地面深度为50 m，宽度为120 m，模型深度为115 m(图12)。设围岩的电阻率为103Ω·m，高阻体的电阻率为105Ω·m，低阻体的电阻率为10 Ω·m。

图12 地下单一矩形异常体模型(模型9)

矩形柱模型的2.5维TEM响应曲线见图13，从图13中可以看到：均匀半空间中高阻、低阻矩形柱体的2.5维TEM响应曲线，在异常体的正上方位置(测点0值位置)均存在明显的异常形态，两者的异常响应形态相反[14]。

图13 矩形柱体2.5维TEM响应曲线

参考模型9各个参数，构建地下高阻、低阻同时存在的2.5维柱状体模型(图14)，测点范围为-590 m～590 m，测点间距为20 m，围岩的电阻率为103Ω·m，高阻体的电阻率为105Ω·m，低阻体的电阻率为10 Ω·m。得到其瞬变电磁响应曲线如图15所示。从图14中可以看到，当同时存在低阻和高阻模型时，低阻体TEM响应曲线非常明显，而高阻体响应曲线尽管也有一定异常形态，但相对低阻体响应曲线而言并不明显。

图14 地下两个矩形异常体模型(模型10)

图15 地下两个矩形异常体TEM响应曲线

3 反演浅析

根据高阻体模型的2.5维正演响应特征，采用阻尼最小二乘反演方法和阻尼特征参数反演法，分别对模型10进行反演分析，两种反演方法结果数据对比见表3，反演结果图见图16。从表3中可以看到，阻尼特征参数反演法均方根相对误差(RMS)和噪信比(RSR)较低，迭代次数较少，收敛较快，反演耗时较少；从图16中可以看出，两种反演方法的效果均能较准确地反演出低阻体的特征，也能在一定程度上反演出高阻体的形态，但阻尼最小二乘法对高阻体的反演结果信息并没有阻尼特征参数法丰富，说明阻尼特征参数反演法对高阻体的反演信息量及效果相对较好。

图16 高阻-低阻矩形柱体模型反演电导率结果

表3 两种反演方法结果数据对比

4 结果与认识

通过建立1维高阻层状模型和2.5维高阻矩形柱模型，模拟了地下高阻体地电模型的瞬变电磁正演响应，讨论了围岩电阻率差异、高阻体电阻率差异等情况下的1维响应特征，分析了高阻矩形柱体的2.5维响应特征，并对其响应结果进行反演分析，得到了一些认识如下：

1)对1维高阻层状模型正演模拟可知：同一层厚时围岩电阻率相对越高，高阻层响应曲线进入衰减的时间相对越早；高阻层越厚，进入衰减过程t1时刻后TEM响应随时间衰减相对越快，达到衰减最大t2时刻较晚，反之亦然。地层层厚相同时，同一模型在高阻层电阻率超过围岩电阻率103倍后，其TEM响应值变化很小(几乎无变化)。

2)对2.5维高阻矩形柱状体模型正演模拟可知：当地下单独存在高阻或低阻异常体时，在异常体位置TEM响应明显，两者的异常响应形态相反。当地下同时存在低阻和高阻异常体时，低阻体TEM响应明显，而高阻体响应并不明显。

3)分别采用阻尼最小二乘反演方法和阻尼特征参数反演法,对2.5维高阻矩形柱状体模型进行反演分析可知：阻尼特征参数反演法均方根相对误差和噪信比较低，迭代次数较少，收敛较快，反演耗时较少，对高阻体的反演结果信息相对丰富、效果较明显。

4)由于影响瞬变电磁响应的因素复杂，笔者在各向同性介质中所建模型存在局限性，在实际情况下，即使相同条件，也很难达到相同效果，这也是今后需要进一步研究的方向。