瞬变电磁法全区视电阻率的二分搜索算法

　[收稿日期]2009202212

　[基金项目]国家自然科学基金项目(40874094)。

　[作者简介]陈清礼(19652),男,1987年大学毕业,博士(后),副教授,现主要从事电磁勘探方面的研究工作。

瞬变电磁法全区视电阻率的二分搜索算法

陈清礼　(长江大学地球物理与石油资源学院,湖北荆州434023)

[摘要]为了快速且精确地计算瞬变电磁法的全区视电阻率,研究了中心回线观测方式感应电动势的特

性,发现均匀半空间中心回线观测方式瞬变电磁测深法的感应电动势随电阻率的增大而单调下降的特性。

据此设计出了计算全区视电阻率的一种快速且精确的算法:二分搜索算法。该算法对电阻率的可能范围

不断进行二等分空间,直到找到一个电阻率,其对应感应电动势与实际观测的电动势相符。理论数据的

计算表明,采用二分搜索算法计算一个测点的全区视电阻率只需要1min 左右,相对误差小于011%,证

明了该算法快速且精确。

[关键词]瞬变电磁法;全区视电阻率;二分搜索算法;中心方式

[中图分类号]P631134[文献标识码]A [文章编号]100029752(2009)022*******

瞬变电磁法(TEM )最早于20世纪30年代提出,到50～60年代,成功地完成了瞬变电磁法的一维正、反演,建立了瞬变电磁法的解释理论和野外工作方法之后,瞬变电磁法才开始进入实用阶段[1]。此后,瞬变电磁法在世界范围内得到了广泛的研究和应用[2～10]。我国于70年代初开始引进和研究TEM ,在石油勘探[11,12]、换流站接地极电阻率测试[13]、地下水寻找[14,15]、地热勘探[16]、工程勘探[17]、矿产调查[18,19]等诸多领域得到了广泛的应用。

由于瞬变电磁场与地下电阻率之间的关系式非常复杂且呈非线性,无法获得全区视电阻率的解析表达式,目前广泛采用的都是早、晚期近似的方案[20]。全区视电阻率的计算一直是研究人员追求的一个目标,但又是困扰研究人员的一个疑难问题。许多学者探讨了如何定义和计算全程视电阻率。白登海给出了一种时间域瞬变电磁法视电阻率的数值计算方法,该方法根据中心方式磁场垂直分量时间变化率的核函数的表现特征,把整个瞬变过程分为早期阶段、早期到晚期的转折点和晚期阶段。分别得到早期视电阻率和晚期视电阻率,然后通过转折点构成一条完整的全程视电阻率曲线[21]。杨生给出了中心回线装置发射电流为斜阶跃波形条件下全区视电阻率迭代反演计算方法[22]。李建平把回线分解为水平电偶极子,然后给出利用电偶极子求取全区视电阻率的方法[23]。王华军依据均匀半空间瞬变响应曲线随地下电导率、发射回线边长与观测时间具有平移伸缩特性,提出了一种直接计算全区视电阻率的方法[24]。付志红等人研究了斜阶跃场源瞬变电磁法的全程视电阻率数值计算[25]。虽然研究的人员很多,但直到现在,还没有一种比较完善的计算全区视电阻率的方法,主要体现在速度、精度和应用条件等方面。

笔者设计出了一种快速、精确的计算中心回线TEM 全区视电阻率算法即二分搜索法,从实测的感应电动势数据直接计算不同时刻的视电阻率。理论计算表明,该算法能更精确地计算全区视电阻率,其相对误差小于011%,同时计算速度非常快,计算一个测点需要的时间不到1min 。

1　算法的理论基础

美国地球物理学家Raab 和Frisehkneecht 推导的中心回线装置的感应电动势表达式[26]为:

V (t )=qI ρπ3/2L 33<(z )-(3z +2z 3)<(z )u (t )(1)

・

54・石油天然气学报(江汉石油学院学报)　

2009年4月　第31卷　第2期Journal of Oil and G as T echnology (J 1J PI )　Apr 12009　Vol 131　No 12

式中,z =2πL/τ;扩散参数τ=2πρt ×107;概率积分<(z )=∫z 02πe -t 2d t ; <(z )=2πe -z 2

;L 为发射回线边长;I 为发射电流;q 为接收线框有效面积;ρ为电阻率。u (t )为阶跃发射电流的电压响应函数,取u (t )=1。

由式(1)可知,感应电动势曲线V (t )是电阻率ρ的一个非线性的复杂函数:V (t )=f (t ,ρ)。理论上说,

给定一条感应电动势曲线,可以由式(1)计算出与之对应的电阻率ρ。均匀半空间的情况下,计算出的是半空间的真电阻率,在不同时刻应该是一个常数;而在非均匀半空间的情况下,则是视电阻率。由于V (t )=

f (t ,ρ

)非常复杂,不能得到其反函数ρ(t )=f -1(t ,V ),因而无法直接由感应电动势计算出视电阻率曲线。因此,不同的学者从多个不同的角度出发提出了各种各样的计算视电阻率的方案。目前常用的是采用早、晚期近似的方法,其最大的缺点是无法获得过渡区的视电阻率。

从抽象关系式V (t )=f (t ,ρ

)不难看出,由实测感应电动势曲线求取全区视电阻率,实质上归结为在可能的电阻率值中搜索出一个合理的值。搜索得到的电阻率值由式(1)计算出的感应电动势与实测感应电动势相等。这样,核心问题就转化为由电阻率值通过式(1)计算出感应电动势。因此,设计良好的数值计算算法需要分析研究感应电动势的特性,以便优化搜索策略。

如果感应电动势随电阻率的增大而具有单调变化的特性,那么就可以设计出一种求取全区视电阻率的算法:二分搜索算法。由于式(1)非常复杂,不能通过理论分析来认识这个特性,只能通过模拟计算来验证是否具有单调特性。为了确保计算的精度,利用MA TL AB 数学工具来完成这个任务。依据式(1),编写了相应的MA TL AB 计算程序,并对时间t 分别为01001、0101、0108、011、1、10、100、1000ms 的感应电动势进行了计算,计算结果如图1所示(图中感应电动势已取对数)。

计算表明,除了在1

μs 和10μs 这2种情况外,从0108ms 到100ms ,感应电动势随电阻率的增加而单调下降。说明了在极早期(t <0108ms )和极晚期(t ≥1s )感应电动势随电阻率不具有单调性外,在0108ms

2　计算全区视电阻率的二分搜索算法

对于在时刻t 的实测感应电动势V 0(t ),需要找到该时刻的一个视电阻率值ρa (t )。该视电阻率值ρa (t )通过式(1)计算的理论感应电动势V (t )与实测感应电动势V 0(t )应该相等。由于在任意时刻,感应电动势随电阻率的增大而单调下降。因此,利用单调下降的特性,设计出了求取全区视电阻率的一种算法,其基本

思想是在电阻率的可能数值范围[ρa ,ρb ]内,二等分区间[ρa ,ρb ],以其中点ρm =(ρa +ρb )/2的电阻率值代

入式(1)计算中点的理论感应电动势V m (t )。如果实测感应电动势V 0(t )大于理论中点的感应电动势

V m (t ),说明视电阻率值应在[ρa ,ρm ]区间而不在[ρm ,ρb ]区间。

这样就把搜索范围减少了一半,提高了计算速度。因此命名该算法为二分搜索算法。具体步骤如下:①设在时刻t 的实测电动势为V 0(t ),设其对应时

间t 的视电阻率为ρa (t ),为了求出ρa (t ),依据工区地层电阻率的最小值ρa 和最大值ρb 可以确定视电阻率

的可能范围是[ρa ,ρb ],如无法确定,一般取为[011,100000]。②计算区间[ρa ,ρb ]的中点ρm =(ρa +ρb )/2,

将ρm 作为电阻率的值代入式(1)中,计算出时刻t 处的理论感应电动势值V (t )。

③用理论计算的感应电动势值V (t )与实测的感应电动势V 0(T )进行比较,如果两者的相对误差小于给定的参数ε,则此时的ρm 就是与实测电动势为V 0(T )相对应的真实视电阻率。计算结束。否则,进入下一步。④如果V 0(t )>V (t ),依据感应电动势随电阻率的增大而单调下降的事实,说明ρm 值比真实视电阻率偏小,真实视电阻率应在

[ρa ,ρm ],把ρm 赋给ρb ;否则,真实视电阻率应在[ρm ,ρb ]区间,把ρm 赋给ρa 。⑤计算区间[ρa ,ρm ]的长度,如果小于给定的精度,则结束;否则,进入下一步。⑥重复步骤②、③、④和⑤。

现在来分析该算法的复杂性。假设工区的电阻率在[0,10000]之间,那么1次分割,区间数变为21,2次分割变为22,N 次分割变为2N 。此时区间的长度为10000/2N ,当N =20时,区间的长度小于0101。这就

・

64・　石油天然气学报(江汉石油学院学报)2009年4月