14-物探化探计算技术-井下瞬变电磁法中全期视电阻率的牛顿求解法

第３６卷　第４期２０１４年７月
物探化探计算技术
ＣＯＭＰＵＴＩＮＧ　ＴＥＣＨＮＩＱＵＥＳ　ＦＯＲ　ＧＥＯＰＨＹＳＩＣＡＬ　ＡＮＤ　ＧＥＯＣＨＥＭＩＣＡＬ　
ＥＸＰＬＯＲＡＴＩＯＮＶｏｌ．３６　Ｎｏ．４
Ｊｕｌ．２０１４
　收稿日期：２０１３－１０－２８ 改回日期：２０１４－０６－
０５基金项目：广西自然科学基金（２０１２０４ＡＡ２３３５）；桂林电子科技大学科研基金（ＵＦ１１０１２
）作者简介：曾庆宁（１９６３－），男，博士，教授，研究方向为现代信号处理及优化算法，Ｅ－ｍａｉｌ：ｑｉｎｇｎｉｎｇｚｅｎｇ
＠１２６．ｃｏｍ。

文章编号：１００１－１７４９（２０１４）０４－０４０１－
０４井下瞬变电磁法中全期视电阻率的牛顿求解法
曾庆宁，张法全，张海如
（桂林电子科技大学　信息与通信学院，桂林　５４１００４
）摘　要：全期视电阻率的求解是瞬变电磁法进行地质反演的关键环节。

针对井下瞬变电磁法中全期视电阻率的隐函数特点，
引进非线性函数求根的牛顿法，给出了井下全期视电阻率的牛顿求解法及其算法步骤。

同时还与二分法做了比较，说明牛顿法与二分法一样，可以满足求解全期视电阻率的任意精度要求，
并且牛顿法比二分法具有更快的收敛速度。

关键词：瞬变电磁法；井下；反演；视电阻率
中图分类号：Ｐ　６３１．３＋
２５　文献标志码：Ａ　ＤＯＩ：１０．３９６９／ｊ
．ｉｓｓｎ．１００１－１７４９．２０１４．０４．０４ 近年发展起来的井下瞬变电磁法，
因其相关设备体积小、成本低、地层分辨率高等优势，成为巷道
超前探测地质异常体的重要方法［
１，７］。

井下瞬变电磁法原理与地面瞬变电磁法原理基本相同，但由于矿井瞬变电磁探测是在井下巷道内进行的，电磁场为全空间而非半空间。

井下与地面一样，主要通过反演地质构造体的视电阻率来发现异常地质，而视电阻率的求解，可通过后期近似公式或前期近似公式进行计算
［１，７］
，但难免存在误差，而
且在某些采样时刻，这种误差可能是无法容忍的。

为此，可考虑使用后全期法和前全期法进行精确求
解［１－８
］，但这种求解法需要使用复杂非线性方程求
根技术。

二分法是一种行之有效的非线性方程求根方法，已成功地用于对地面瞬变电磁法视电阻率的
任意精度求解［２－３
］，该方法亦可用于井下瞬变电磁法视电阻率的任意精度求解［１］。

然而，二分法通常
需要迭代的次数很多，收敛速度较慢。

本次研究引进非线性方程求根的牛顿方法，并将其用于井下瞬变电磁法视电阻率的任意精度求解。

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

１　全期视电阻率
在井下瞬变电磁法中，通常采用的是重叠回线或中心回线，其接收线圈于时刻ｔ对二次场产生的
感应电压为［
１］
Ｖ（ｔ）＝μ０
ＩＳｒ槡πａｔＦ（ｕ
）（１
）其中
Ｆ（ｕ）＝ｕ３ｅ－
ｕ
２
（２
）ｕ＝
ａ
２
μ０ρ
槡
ｔ（３
）其中μ０＝４π×１０－
７为真空的磁导率；Ｓｒ为接收天线的有效面积；Ｉ为发射电流，ａ为发射线框的有效
半径，可通过ａ＝Ｓ／槡π计算；
Ｓ为发射线框的有效面积。

由感应电压Ｖ（ｔ）通过式（１）求解视电阻率ρ的方法称为全期法。

式（２）所示的Ｆ（ｕ）称为核函数，
对其求导数可得
Ｆ′（ｕ）＝（３－２ｕ２）ｕ２ｅ－
ｕ
２
（４
）在（０，＋∞）上Ｆ
′（ｕ）有唯一的零点ｕ０＝
槡
３２
≈１．
２２４７（５）容易证明：Ｆ（ｕ）在０＜ｕ≤ｕ０时，
是ｕ的单调递增函数，Ｆ（ｕ）在ｕ≥ｕ０时，是ｕ的单调递减函数。

因此，对给定的时间ｔ，由式（１）求解视电阻率ρ通常可获得两个解：①可在（０，ｕ０］之间获得，
由此获得的解称之为“后全期解”，并称相应的方法为“后全期法”；②可在［ｕ０，＋∞）之间获得，由此获得的解称之为“
前全期解”，并称相应的方法为“前全期法”。

对实测的感应电压Ｖ０（ｔ），由于各方面的噪声与误差效应，有时难免导致式（１）求解视电阻率ρ可能无解，这时可使用αＶ０（ｔ）代替Ｖ０（ｔ）的方法（α＜
１
），通过重新定义视电阻率而获得解决［４］。

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

对此，如果令
ｇ（ρ）＝
μ０
ＩＳｒ槡
πａｔＦ（ｕ）－Ｖ０（ｔ）（６）则ρ的求解问题转化为函数ｇ（ρ）的求根问题。

由
式（３）和Ｆ（ｕ）的单调性，ｇ（ρ）将在（０，ρ
０］内单调递增，在［ρ
０，＋∞）内单调递减，其中ρ０＝μ０ａ２４ｕ２０
ｔ＝μ０ａ２
６ｔ（７）因此视电阻率ρ的后全期解与前全期解可分别在［ρ０，＋∞）与（０，ρ０］中求ｇ（ρ）的根而获得。

二分法和牛顿法都是对单调函数进行任意精确求根的常用方法。

文献［３］分别详细描述了二分法求解半空间全期视电阻率的具体步骤，而全空间全期视电阻率的二分法算法步骤与半空间时是完全类似的
［１］。

牛顿法与二分法相比，通常具有所需迭代次数少、收敛速度更快的特点。

下面先介绍牛顿求根法，然后给出全期视电阻率的牛顿求解法。

２　牛顿求根法
设函数ｆ（ｘ）在区间［ｂ，ｃ
］单调并且有根。

不失一般性，不妨设ｆ（ｘ）为单调递减函数。

如图（１）所示，任取初始值ｘ０∈［ｂ，ｃ］，对ｉ＝１，２，…，逐步迭代计算ｆ（ｘ）于点（ｘｉ－１，ｆ（ｘｉ－
１））处的切线ＡＢ在ｘ轴的截距ｘｉ，则ｘｉ将收敛于ｆ（
ｘ）在区间［ｂ，ｃ
］的唯一根ｘ＊，
这就是牛顿求根法的基本思想。

图１　牛顿法示意图
Ｆｉｇ．１　Ｄｉａｇ
ｒａｍｍａｔｉｃ　ｓｋｅｔｃｈ　ｏｆ　Ｎｅｗｔｏｎ　ｍｅｔｈｏｄ切线ＡＢ的斜率为ｆ（ｘ）在ｘｉ－
１的导数ｆ′（ｘｉ－
１），因此，切线ＡＢ的方程为ｙ＝ｆ′（ｘｉ－１）（ｘ－ｘｉ－１）＋ｆ（ｘｉ－
１）（８）设该直线与ｘ轴的截距为ｘｉ，
则ｘｉ＝ｘｉ－
１－ｆ（ｘｉ－
１）ｆ′ｘｉ－
１（９
）式（９
）即为牛顿求根法的迭代公式。

３　全期视电阻率的牛顿求解法
根据上述第一节，对给定的时间，视电阻率ρ的全期解有两个，分别为后全期解与前全期解。

后全期解可在［ρ０，ρｂ］内用牛顿法求解，而前全期解可在［ρａ，ρ
０］内用牛顿法求解，其中ρｂ为视电阻率上界的一个估值，通常可取ρｂ＝１０５
Ω·ｍ，ρ
ａ为式电阻率下界的一个估值，通常可取ρａ＝
１０－５Ω·ｍ。

由式（９）可知，牛顿求解法中需用到函数的导数。

因此由式（６）求解ｇ（ρ）的根时，需用到ｇ（ρ）的导数，容易得出ｇ′（ρ）为
ｇ′（ρ）＝
μ０ＩＳｒ槡πａｔＦ′（ｕ）·ｕ′ρ（１０
）其中
Ｆ′（ｕ）＝（３－２ｕ２）ｕ２ｅ－
ｕ
２
（１１）ｕ＝
ａ
２
μ０ρ
槡
ｔ（１２
）ｕ′ρ＝－
ａ４
μ
０槡
ｔ
ρ－３２（１３
）于是视电阻率ρ的后（前）全期解的牛顿解法步
骤可描述如下：
２
０４ 物探化探计算技术３６卷
步骤１：任意给定视电阻率ρ的精度误差ε＞０，任取ρ（０）∈（ρ０，ρｂ），令ｉ＝０。

步骤２：ｉ＝ｉ＋１，按式（６）、式（２）、式（３）计算ｇ［ρ（ｉ－１）］，并按式（１０）、式（１１）、式（１２）、式（１３）计算ｇ′［ρ（ｉ－１）］。

步骤３：按下式计算ρ（ｉ）：
ρ（ｉ）＝ρ（ｉ－１）－ｇ［ρ（ｉ－１）］
ｇ′ρｉ－１。

步骤４：如果｜ρ（ｉ）－ρ（ｉ－１）｜＜ε，则迭代停止，ρ（ｉ）即为所求之解，否则转步骤２继续进行迭代。

４　牛顿法与二分法的比较
二分法与牛顿法均为对井下全期视电阻率进行任意精确求解的方法。

二分法原理简单，实现容易，但收敛速度较慢；而牛顿法通常具有比二分法快得多的收敛速度。

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

假设井下全空间为均匀介质空间，介质的视电阻率８０Ωｍ。

采用同点探测方式，发射天线有效面积２００ｍ２，接收天线有效面积６０ｍ２，发射电流２Ａ，关断时间为１３０μｓ，时间均匀采样的采样间隔为２０μｓ，测道数为１２８。

于是，通过式（１）－式（３）可获得正演数据Ｖ（ｔｉ），ｉ＝１、２、…、１２８。

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

表１　二分法与牛顿法平均迭代次数
Ｔａｂ．１　Ｉｔｅｒａｔｉｏｎｓ　ｏｆ　ｂｉｓｅｃｔｉｏｎ　ａｎｄ　Ｎｅｗｔｏｎ　ｍｅｔｈｏｄｓ
精度，次１０－１　１０－２　１０－３　１０－５　１０－６　１０－８
二分法２４　２７　３０　３７　４０　４７
牛顿法４　４　５　５　６　６在表１数据的计算过程中，假设视电阻率ρ在１０－８Ω·ｍ～１０６Ω·ｍ之间，牛顿法的初始值取为ρ（０）＝５０。

计算结果表明：二分法和牛顿法均能很好地收敛于给定的视电阻率。

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

应当注意：牛顿法对函数导数及初值较为敏感，而且前述的ｇ′（ρ）的绝对值在很多地方很小，因此，实际计算后（前）全期视电阻率时，应在前述的步骤４中增加对ρ（ｉ）是否仍落在区间［ρ０，ρｂ］（［ρａ，ρ０］）的判断，如果超出区间范围，则应停止迭代。

此外，牛顿法与二分法一样，对特别后期（或特别前期）的采样时刻，由于感应电压的值太小且实际测量中又难免随机噪声的影响，容易导致解的较大误差，这时最好用后期近似公式法（或前期近似公式法）代替。

５　应用实例
实测数据来源于对陕西陈家山煤矿某井下巷道掘进头处的ＴＥＭ探测，目的是探测采煤掘进前方是否有积水，以避免透水事故发生。

探测时，分别对前方上倾１５°扇面、水平方向扇面和下倾１５°扇面进行探测。

每个扇面的探测范围从左边－４５°至右边４５°，每１５°一个测点，共７个测点。

发射天线有效面积８ｍ２，接收天线有效面积４０ｍ２，采用共轴方式，发射与接收天线相距７ｍ，发射电流２Ａ，采样频率２５Ｈｚ，每测点３０测道，关断时间为１４０μｓ。

图２为上述水平方向扇面反演剖面图。

反演中，视电阻率采用后全期解，且均使用牛顿法求解，初值均取ρ（０）＝５０Ω·ｍ，而ρａ＝１０－５Ω·ｍ，ρｂ＝１０５Ω·ｍ，ε＝１０－３。

由图２可见，在掘进头水平扇形探测面的右偏２５°～３１°、距离８０ｍ～１５０ｍ存在低阻异常区域。

事后证明，此区域正是该煤矿已经充水的一个采空区，与反演结果完全吻合。

图２　水平扇形反演图
Ｆｉｇ．２　Ｈｏｒｉｚｏｎｔａｌ　ｆａｎｌｉｇｈｔ　ｉｎｖｅｒｓｉｏｎ　ｐｉｃｔｕｒｅ
６　结束语
对井下全期视电阻率的求解，给出了牛顿求解法及其具体的算法步骤。

牛顿求解法与二分法一样，可用于求后全期和前全期任意精度要求的解，但
３
０
４
４期曾庆宁等：井下瞬变电磁法中全期视电阻率的牛顿求解法
牛顿法所需的迭代次数比二分法少得多，具有更快的收敛速度。

参考文献：
［１］　杨海燕，邓居智，张华，
等．矿井瞬变电磁法全空间视电阻率解释方法研究［Ｊ］．地球物理学报，２０１０，５３（３）：６５１－６５６．
［２］　张成范，
翁爱华，孙世栋，等．计算矩形大定源回线瞬变电磁测深全区视电阻率［Ｊ］．吉林大学学报：地球科学版，２００９，３９（４）：７５５－７５８．
［３］　李文尧，
晏冲为．中心回线瞬变电磁法全期视电阻率的二分法求解［Ｊ］．昆明理工大学学报：自然科学版，２０１３，３８（２）：２６－３３．
［４］　白登海，ＭＡＸＷＥＬＬ　
Ａ　ＭＥＪＵ，卢健，等．时间域瞬变电磁法中心方式全区视电阻率的数值计算［Ｊ］．地球物理学报，２００３，４６（９）：６９７－７０４．
［５］　翁爱华，
陆冬华，刘国兴．利用连分式定义瞬变电磁法全区视电阻率研究［Ｊ］．煤田地质与勘探，２００３，３１（３）：５６－５９．
［６］　苏朱刘，
胡文宝．中心回线方式瞬变电磁测深虚拟全区式电阻率和一维反演方法［Ｊ］．石油物探，２００２，４１（２）：２１６－２２１．
［７］　李好，
胡运兵，吴燕清．应用矿井瞬变电磁法超前探测煤矿井下含水体［Ｊ］．矿业安全与环保，２０１２，３９（５）：４９－５２．
［８］　牛之琏．
时域电磁法原理［Ｍ］．长沙：中南大学出版社，２００７．
Ｎｅｗｔｏｎ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｓｏｌｖｉｎｇ　ｆｕｌｌ－ｔｉｍｅ　ａｐｐａｒｅｎｔ　ｒｅｓｉｔｉｖｉｔｙ　
ｏｆ　ｕｎｄｅｒｇｒｏｕｎｄ－ｍｉｎｅｔｒａｎｓｉｅｎｔ　ｅｌｅｃｔｒｏｍａｇ
ｎｅｔｉｃ　ｍｅｔｈｏｄＺＥＮＧ　Ｑｉｎｇ－ｎｉｎｇ，
ＺＨＡＮＧ　Ｆａ－ｑｕａｎ，ＺＨＡＮＧ　Ｈａｉ－ｒｕ（Ｃｏｌｌｅｇｅ　ｏｆ　Ｉｎｆｏｒｍａｔｉｏｎ　ａｎｄ　Ｔｅｌｅｃｏｍｍｕｎｉｃａｔｉｏｎｓ，Ｇｕｉｌｉｎ　Ｕｎｉｖｅｒｓｉｔｙ　
ｏｆ　Ｅｌｅｃｔｒｏｎｉｃ　Ｔｅｃｈｎｏｌｏｇｙ，Ｇｕｉｌｉｎ　５４１００４）Ａｂｓｔｒａｃｔ：Ｓｏｌｖｉｎｇ　ｆｕｌｌ－ｔｉｍｅ　ａｐｐａｒｅｎｔ　ｒｅｓｉｔｉｖｉｔｙ　ｉｓ　ｔｈｅ　ｋｅｙ　ｓｔｅｐ　ｆｏｒ　ｇｅｏｌｏｇｉｃａｌ　ｉｎｖｅｒｓｉｏｎ　ｏｆ　Ｔｒａｎｓｉｅｎｔ　Ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ　Ｍｅｔｈｄ（ＴＥＭ）．Ｉｎ　ｔｈｉｓ　ｐａｐｅｒ，Ｎｅｗｔｏｎ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｓｏｌｖｉｎｇ　ｔｈｅ　ｒｏｏｔ　ｏｆ　ｎｏｎ－ｌｉｎｅａｒ　ｆｕｎｃｔｉｏｎ　ｉｓ　ｉｎｔｒｏｄｕｃｅｄ　ａｎｄ　ｔｈｅｎ　ｕｓｅｄ　ｆｏｒ　ｓｏｌｖｉｎｇｔｈｅ　ｆｕｌｌ－ｔｉｍｅ　ａｐｐａｒｅｎｔ　ｒｅｓｉｔｉｖｉｔｙ　ｏｆ　ｕｎｄｅｒｇｒｏｕｎｄ－ｍｉｎｅ　ＴＥＭ，ａｃｃｏｒｄｉｎｇ　ｔｈｅ　ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ　ｏｆ　ｔｈｅ　ｉｍｐｌｉｃｉｔ　ｆｕｎｃｔｉｏｎ　ｏｆ　ｔｈｅ　ａｐ－ｐａｒｅｎｔ　ｒｅｓｉｔｉｖｉｔｙ．Ｍｅａｎｗｈｉｌｅ，ｔｈｅ　ｃｏｍｐａｒｉｓｏｎ　ｏｆ　Ｎｅｗｔｏｎ　ａｌｇｏｒｉｔｈｍ　ａｎｄ　ｔｈｅ　ｂｉｓｅｃｔｉｏｎ　ａｌｇｒｉｔｈｍ　ｉｎｄｉｃａｔｅｓ　ｔｈａｔ　ｔｈｅ　ｐｒｅｓｅｎｔ　Ｎｅｗｔｏｎａｌｇｏｒｉｔｈｍ　ｃａｎ　ａｐｐｒｏｘｉｍａｔｅ　ｔｈｅ　ｆｕｌｌ－ｔｉｍｅ　ａｐｐａｒｅｎｔ　ｒｅｓｉｔｉｖｉｔｙ　ｗｉｔｈ　ａｎｙ　ａｃｃｕｒａｃｙ　ｊｕｓｔ　ｌｉｋｅ　ｔｈｅ　ｂｉｓｅｃｔｉｏｎ　ａｌｇｏｒｉｔｈｍ，ｂｕｔ　ｔｈｅ　Ｎｅｗｔｏｎａｌｇｏｒｉｔｈｍ　ｍａｙ　
ｃｏｎｖｅｒｇｅ　ｍｕｃｈ　ｆａｓｔｅｒ　ｔｈａｎ　ｔｈｅ　ｂｉｓｅｃｔｉｏｎ　ａｌｇｏｒｉｔｈｍ．Ｋｅｙ　
ｗｏｒｄｓ：ｔｒａｎｓｉｅｎｔ　ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ　ｍｅｔｈｏｄ；ｕｎｄｅｒｇｒｏｕｄ－ｍｉｎｅ；ｉｎｖｅｒｓｉｏｎ；ａｐｐａｒｅｎｔ　ｒｅｓｉｔｉｖｉｔｙ４
０４ 物探化探计算技术３６卷。