井下瞬变电磁法中全期视电阻率的牛顿求解法

曾庆宁 ，张法全， 张海如

(桂林电子科技大学 信息与通信学院，桂林 541004)

近年发展起来的井下瞬变电磁法，因其相关设备体积小、成本低、地层分辨率高等优势，成为巷道超前探测地质异常体的重要方法[1，7]。

井下瞬变电磁法原理与地面瞬变电磁法原理基本相同，但由于矿井瞬变电磁探测是在井下巷道内进行的，电磁场为全空间而非半空间。井下与地面一样，主要通过反演地质构造体的视电阻率来发现异常地质，而视电阻率的求解，可通过后期近似公式或前期近似公式进行计算[1，7]，但难免存在误差，而且在某些采样时刻，这种误差可能是无法容忍的。为此，可考虑使用后全期法和前全期法进行精确求解[1-8]，但这种求解法需要使用复杂非线性方程求根技术。二分法是一种行之有效的非线性方程求根方法，已成功地用于对地面瞬变电磁法视电阻率的任意精度求解[2-3]，该方法亦可用于井下瞬变电磁法视电阻率的任意精度求解[1]。然而，二分法通常需要迭代的次数很多，收敛速度较慢。

本次研究引进非线性方程求根的牛顿方法，并将其用于井下瞬变电磁法视电阻率的任意精度求解。通过实例，说明了牛顿法与二分法一样，用于求解井下全期视电阻率是可行的，并且 在同等精度要求下，牛顿法比二分法所需的迭代次数要少许多，明显具有更快的收敛速度，为井下瞬变电磁法基于视电阻率的地质反演，提供了一种满足任意精度要求且速度更快的运算方法。

1 全期视电阻率

在井下瞬变电磁法中，通常采用的是重叠回线或中心回线，其接收线圈于时刻t对二次场产生的感应电压为[1]

(1)

其中

F(u)=u3e-u2

(2)

(3)

由感应电压V(t)通过式(1)求解视电阻率ρ的方法称为全期法。式(2)所示的F(u)称为核函数，对其求导数可得

F′(u)=(3-2u2)u2e-u2

(4)

在(0,+∞)上F′(u)有唯一的零点

(5)

容易证明：F(u)在 0

对实测的感应电压V0(t)，由于各方面的噪声与误差效应，有时难免导致式(1)求解视电阻率ρ可能无解，这时可使用αV0(t)代替V0(t)的方法(α<1)，通过重新定义视电阻率而获得解决[4]。

求解视电阻率ρ的后全期解与前全期解，是瞬变电磁法进行正确反演和资料解释的关键环节，获得视电阻率ρ的后全期精确值和前全期精确值是我们所期望的。对此，如果令

(6)

则ρ的求解问题转化为函数g(ρ)的求根问题。由式(3)和F(u)的单调性，g(ρ)将在(0,ρ0]内单调递增，在[ρ0,+∞)内单调递减，其中

(7)

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

二分法和牛顿法都是对单调函数进行任意精确求根的常用方法。文献[3]分别详细描述了二分法求解半空间全期视电阻率的具体步骤，而全空间全期视电阻率的二分法算法步骤与半空间时是完全类似的[1]。

牛顿法与二分法相比，通常具有所需迭代次数少、收敛速度更快的特点。下面先介绍牛顿求根法，然后给出全期视电阻率的牛顿求解法。

2 牛顿求根法

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

图1 牛顿法示意图Fig．1 Diagrammatic sketch of Newton method

切线AB的斜率为f(x)在xi-1的导数f′(xi-1)，因此，切线AB的方程为

y=f′(xi-1)(x-xi-1)+f(xi-1)

(8)

设该直线与x轴的截距为xi，则

(9)

式(9)即为牛顿求根法的迭代公式。

3 全期视电阻率的牛顿求解法

根据上述第一节，对给定的时间，视电阻率ρ的全期解有两个，分别为后全期解与前全期解。后全期解可在[ρ0,ρb]内用牛顿法求解，而前全期解可在[ρa,ρ0]内用牛顿法求解，其中ρb为视电阻率上界的一个估值，通常可取ρb=105Ω·m，ρa为式电阻率下界的一个估值，通常可取ρa=10-5Ω·m。

由式(9)可知，牛顿求解法中需用到函数的导数。因此由式(6)求解g(ρ)的根时，需用到g(ρ)的导数，容易得出g′(ρ)为

(10)

其中

F′(u)=(3-2u2)u2e-u2

(11)

(12)

(13)

于是视电阻率ρ的后(前)全期解的牛顿解法步骤可描述如下：

步骤1：任意给定视电阻率ρ的精度误差ε>0，任取ρ(0)∈(ρ0,ρb)，令i=0。

步骤2：i=i+1，按式(6)、式(2)、式(3)计算g[ρ(i-1)]，并按式(10)、式(11)、式(12)、式(13)计算g′[ρ(i-1)]。

步骤3：按下式计算ρ(i):

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

4 牛顿法与二分法的比较

二分法与牛顿法均为对井下全期视电阻率进行任意精确求解的方法。二分法原理简单，实现容易，但收敛速度较慢；而牛顿法通常具有比二分法快得多的收敛速度。

通过试例，比较二分法与牛顿法两者的收敛速度。

假设井下全空间为均匀介质空间，介质的视电阻率80 Ωm。采用同点探测方式，发射天线有效面积200 m2，接收天线有效面积60 m2，发射电流2A，关断时间为130 μs，时间均匀采样的采样间隔为20 μs，测道数为128。于是，通过式(1)-式(3)可获得正演数据V(ti),i=1、2、…、128。

表1列出了对全期视电阻率不同精度要求下，通过V(ti)，用二分法和牛顿法在所有测道求解后全期视电阻率所需的平均迭代次数。

表1 二分法与牛顿法平均迭代次数

在表1数据的计算过程中，假设视电阻率ρ在10-8Ω·m～106Ω·m之间，牛顿法的初始值取为ρ(0)=50。计算结果表明：二分法和牛顿法均能很好地收敛于给定的视电阻率。

从表1中容易看出，牛顿法仅需几次迭代即可达到很高的精度，在同等精度要求下，牛顿法所需的迭代次数比二分法要少许多。

应当注意：牛顿法对函数导数及初值较为敏感，而且前述的g′(ρ)的绝对值在很多地方很小，因此，实际计算后(前)全期视电阻率时，应在前述的步骤4中增加对ρ(i)是否仍落在区间[ρ0,ρb]([ρa,ρ0])的判断，如果超出区间范围，则应停止迭代。此外，牛顿法与二分法一样，对特别后期(或特别前期)的采样时刻，由于感应电压的值太小且实际测量中又难免随机噪声的影响，容易导致解的较大误差，这时最好用后期近似公式法(或前期近似公式法)代替。

5 应用实例

实测数据来源于对陕西陈家山煤矿某井下巷道掘进头处的TEM探测，目的是探测采煤掘进前方是否有积水，以避免透水事故发生。探测时，分别对前方上倾15°扇面、水平方向扇面和下倾15°扇面进行探测。每个扇面的探测范围从左边-45°至右边45°，每15°一个测点，共7个测点。发射天线有效面积8 m2，接收天线有效面积40 m2，采用共轴方式，发射与接收天线相距7 m，发射电流2 A，采样频率25 Hz，每测点30测道，关断时间为140 μs。

图2为上述水平方向扇面反演剖面图。反演中，视电阻率采用后全期解，且均使用牛顿法求解，初值均取ρ(0)=50 Ω·m，而ρa=10-5Ω·m，ρb=105Ω·m，ε=10-3。由图2可见，在掘进头水平扇形探测面的右偏25°～31°、距离80 m～150 m存在低阻异常区域。事后证明，此区域正是该煤矿已经充水的一个采空区，与反演结果完全吻合。

图2 水平扇形反演图Fig．2 Horizontal fanlight inversion picture

6 结束语

对井下全期视电阻率的求解，给出了牛顿求解法及其具体的算法步骤。牛顿求解法与二分法一样，可用于求后全期和前全期任意精度要求的解，但牛顿法所需的迭代次数比二分法少得多，具有更快的收敛速度。

参考文献：

[1] 杨海燕, 邓居智, 张华,等. 矿井瞬变电磁法全空间视电阻率解释方法研究[J]. 地球物理学报, 2010, 53(3): 651-656.

[2] 张成范,翁爱华,孙世栋,等.计算矩形大定源回线瞬变电磁测深全区视电阻率[J].吉林大学学报：地球科学版, 2009,39(4):755-758.

[3] 李文尧，晏冲为. 中心回线瞬变电磁法全期视电阻率的二分法求解[J]. 昆明理工大学学报：自然科学版, 2013, 38(2):26-33.

[4] 白登海，MAXWELL A MEJU, 卢健，等. 时间域瞬变电磁法中心方式全区视电阻率的数值计算[J]. 地球物理学报，2003,46(9)： 697-704.

[5] 翁爱华，陆冬华，刘国兴. 利用连分式定义瞬变电磁法全区视电阻率研究[J]. 煤田地质与勘探，2003，31(3): 56-59.

[6] 苏朱刘，胡文宝. 中心回线方式瞬变电磁测深虚拟全区式电阻率和一维反演方法[J]. 石油物探，2002，41(2): 216-221.

[7] 李好，胡运兵，吴燕清. 应用矿井瞬变电磁法超前探测煤矿井下含水体[J]. 矿业安全与环保，2012，39(5): 49-52.

[8] 牛之琏.时域电磁法原理[M].长沙：中南大学出版社，2007.