瞬变电磁全期视电阻率的弦截无限逼近算法

曾庆宁 罗 瀛 刘 帅 龙 超

(①桂林电子科技大学认知无线电与信息处理教育部重点实验室，广西桂林 541004；②桂林电子科技大学信息与通信学院，广西桂林 541004)

0 引言

瞬变电磁法(TEM)在石油勘探 、矿产调查、地下水勘察、工程地质勘探等诸多领域得到了广泛应用[1-6]。传统的瞬变电磁反演是将测量的感应电动势转化为地质体的视电阻率[7-10]，从而得到探区真实的地质图像。

可通过近似公式法[1,3](包括后期近似公式与前期近似公式)或者改进的近似公式法[11-13]，由感应电动势求解视电阻率，这些方法计算相对简单，但通常存在较大误差，而且在某些时段，尤其是在反映中浅地层地质的早期时段，这种误差可能导致错误的地质反演结果。特别是当目标体积较小或目标体的电阻率与围岩的电阻率相近时，这种误差更可能导致错误反演结果。为此，可使用后全期法和前全期法进行无限精确求解[14-17]，以克服这种误差的影响，使反演成像更为精细和准确，但这种无限精确求解法需要使用复杂非线性函数的求根技术。

二分法和牛顿法是两种行之有效的非线性函数无限精确求根方法[18-21]，前者已成功应用于地面瞬变电磁数据的视电阻率任意精度求解[14-17]，后者成功应用于井下瞬变电磁数据的视电阻率任意精度求解[18]。然而，二分法通常需要较多的迭代次数，收敛速度较慢；牛顿法需要首先求出瞬变电磁感应电动势函数的导数，对于井下瞬变电磁信号，该导数的表达式并不复杂，但对于地面瞬变电磁信号，其表达式则非常复杂，尤其是其计算误差很大，往往导致算法失效。

本文引进非线性函数求根的弦截方法[21]，并将其用于地面重叠回线或中心回线瞬变电磁视电阻率的任意精度求解，给出了弦截无限逼近算法。文中还进一步说明了在同等精度要求下，弦截法比二分法所需的迭代次数更少，具有更快的收敛速度。同时，弦截法可避免类似牛顿法那样需要对瞬变电磁信号求导的问题，为瞬变电磁数据的视电阻率反演提供了另一种满足任意精度要求且收敛速度较快的计算方法。最后针对地下水仓的实测瞬变电磁数据，应用弦截法求解其全期视电阻率的无限精确解，并据此进行了地质反演成像，效果较好。

1 全期视电阻率

瞬变电磁法中，经常采用重叠回线或中心回线法采集数据，接收线圈于时刻t测得的感应电动势为[1]

(1)

式中核函数

(2)

(3)

通过式(1)求解视电阻率ρ的方法称为全期法。可以证明，F(z)是z的上凸函数[20-21]，因此可以找到唯一的点z0，使F(z)在z≤z0时是z的单调递增函数，在z≥z0时是单调递减函数[17-18]。通过计算，得到z0=1.613630000000744。因此，对给定的时间t，由式(1)求解视电阻率ρ通常可获得两个解: 在区间(0,z0]获得的解称为“后全期解”，并称相应的方法为“后全期法”；在[z0,+∞)获得的解称为“前全期解”，并称相应的方法为“前全期法”。

从式(3)可以看出，z与ρt有关。如果t较大，则ρt较大，因而当z≤z0时用“后全期解”或“后全期法”较合理；反之，如果t较小，则应该用“前全期解”或“前全期法”。

对实测的感应电动势V0(t)，由于噪声和误差的影响，可能导致利用式(1)求解视电阻率ρ时出现无解，这时可使用αV0(t)(α<1)代替V0(t)，重新求解视电阻率，参数α的定义和计算参见文献[19]。

求解视电阻率ρ的后全期解与前全期解，获得其精确解是瞬变电磁数据反演和解释的关键环节。令

(4)

则ρ的求解问题转化为求取非线性函数g(ρ)的解。由式(3)及F(z)的上凸特性，函数g(ρ)在(0,ρ0]内单调递增，在[ρ0,+∞)内单调递减，其中

(5)

因此视电阻率ρ的后全期解与前全期解可分别在[ρ0,+∞)与(0,ρ0]区间通过求g(ρ)的解获得。

二分法、牛顿法是对非线性单调函数进行任意精确求根的方法，文献[14-18]详细描述了利用这两种方法求解全期视电阻率的具体做法。弦截法也是对单调函数进行任意精确求根的方法，但弦截法与二分法相比，具有迭代次数少、收敛快的特点，而且无须象牛顿法那样对函数求导。下面先介绍弦截求根法，然后给出全期视电阻率的弦截求解方法和步骤。

2 弦截求根法

设函数f(x)在区间[b,c]内单调并且有根。不失一般性，不妨设f(x)为单调递减函数。如图1所示，任取初始值x1,x2∈[b,c]，逐步迭代计算f(x)过点(xi-1,f(xi-1))和(xi-2,f(xi-2))的弦截线AB在x轴的截距xi(i=3,4,…)，则xi将收敛于f(x)在区间[b,c]的唯一根x*，这就是弦截求根法的基本思想。

弦截线AB的表达式为

图1 弦截法示意图

设直线AB与x轴的截距为xi，则弦截法无限逼近求根的迭代公式为

(6)

3 全期视电阻率的弦截逼近算法

根据上述原理，对每个给定的时间点，都有两个视电阻率ρ的全期解，分别为后全期解与前全期解。后全期解可在[ρ0,ρb]内用弦截法求解，而前全期解可在[ρa,ρ0]内用弦截法求解，其中ρb为视电阻率上界的一个估值，通常取ρb=1.0×106Ω·m;ρa为视电阻率下界的一个估值，通常取ρa=1.0×10-6Ω·m。

视电阻率ρ的后全期解的弦截解法步骤如下：

(1)任意给定视电阻率ρ的精度ε>0，任取ρ(1)、ρ(2)∈[ρ0,ρb]；令i=2，运用式(2)～式(4)计算g[ρ(i)]和g[ρ(i-1)]。

(2)令i=∶i+1，计算g[ρ(i)]。

(3)计算

ρ(i)=ρ(i-2)-g[ρ(i-2)]×

(4)若|ρ(i)-ρ(i-1)|<ε，则迭代停止，ρ(i)即为所求之解； 否则，转步骤(2)继续进行迭代计算。

前全期解计算步骤据此类推，不再赘述。

4 理论模型计算效果分析

视电阻率的计算方法可分为近似公式法和无限逼近法。近似公式法计算量小且算法稳定，但得到的视电阻率误差较大，尤其是在反映中浅层地质状况的早期采样时间段，相对误差有时可达5%以上。而无限逼近法可获得误差任意小的视电阻率，缺点是计算量大，稳定性也不如近似公式法。二分法[14-17]、牛顿法[18]及本文提出的弦截法均为无限逼近法，理论上而言，牛顿法收敛最快，弦截法次之，二分法最慢[20-24]。牛顿法需要利用瞬变电磁信号的导数函数，对于式(1)～式(3)表示的地面瞬变电磁信号而言，其导数表达式过于复杂，其计算误差影响也很大，因此，牛顿法并不适用于地面瞬变电磁法的视电阻率计算，而二分法和弦截法无需导数计算。下面通过两个理论模型，比较二分法与弦截法的计算效果。

4.1 均匀半空间

假设均匀半空间的电阻率为80Ω·m。采用中心回线探测方式，发射天线有效面积为200m2，接收天线的有效面积为60m2，发射电流为2A，关断时间为130μs，采样间隔为20μs，共采集128个数据。通过式(1)～式(3)可获得正演数据V(ti),i=1,2,…,128。

表1列出了在均匀半空间全期视电阻率计算过程中不同精度条件下，基于V(ti)分别用二分法和弦截法对所有点求解全期视电阻率所需的平均迭代次数。这里假设视电阻率ρ的变化范围为1.0×10-8～1.0×106Ω·m， 弦截法的初始值ρ(1)=40Ω·m,ρ(2)=50Ω·m。

4.2 层状模型

假设地下包含三个电性层，各层电阻率分别为50、8、70Ω·m，层厚分别为150、50m、+∞，其他各项参数同均匀半空间模型。表2为不同方法在不同精度要求下的迭代次数统计。

表1 均匀半空间模型二分法与弦截法计算视电阻率的平均迭代次数统计

表2 层状模型二分法与弦截法计算视电阻率的平均迭代次数统计

对比表1与表2可以看出，在同等精度要求时，弦截法所需的迭代次数比二分法小得多。由于每次迭代两种算法的运算量基本相当，因此弦截法所需的总运算量也比二分法少得多。

需特别指出的是，弦截法与牛顿法对初值都较敏感，实际计算后/前全期视电阻率时，应在前述的步骤4中增加对ρ(i)是否仍落在区间[ρ0,ρb]或[ρa,ρ0]进行判断； 如果超出范围，则应停止迭代，改用其他方法(比如近似公式法)。此外，在极后期或极前期，由于感应电动势太小且随机噪声较大，与二分法和牛顿法一样，如果利用弦截法计算视电阻率，解的误差会较大，这种情况下，采用近似公式法计算视电阻率，误差会小一些且计算过程更稳健。

5 应用实例

实测数据采集自广西州景矿区，地下300m处有一水仓。使用地面大定源回线TEM探测法，大回线设计在水仓的正上方，呈矩形状，矩形边长分别为500m和400m；在回线内部设计4条测线，间距为30m，每条测线有8个测点，间距为25m；每个测点的二次场测道为62道，测道先密后疏；接收天线有效面积为100m2，发射电流为15A，采样频率为6.25Hz，关断时间为250μs，数据叠加次数为16。任一测点与大回线任意一边的最近距离均大于120m，因此可使用瞬变电磁二次感应电动势的中心回线公式(式(1))进行反演。

反演时对浅层盲区进行填充，视电阻率采用后全期解，且均使用弦截法进行求解，初值均取ρ(1)=30Ω·m,ρ(2)=40Ω·m,ρa=1.0×10-5Ω·m,ρb=1.0×105Ω·m，逼近精度ε=1.0×10-3。图2为其中一条过水仓测线的视电阻率反演剖面，可见反演结果很好地展现了地下300m深度水仓的准确位置，甚至水仓的大致形状及大小(水仓为低阻，即图中深蓝色区域)也清晰可见。

图2 过水仓测线的视电阻率反演剖面

6 结束语

对瞬变电磁法中全期视电阻率的求解，给出了弦截无限逼近求解法及其具体的算法步骤。弦截求解法与二分求解法、牛顿求解法一样，可用于求解后全期和前全期视电阻率任意精度解，但弦截法无需像牛顿法那样计算瞬变电磁信号的导数，避免了复杂的导数计算，且所需的迭代次数比二分法少得多，具有更快的收敛速度。