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

[收稿日期]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 <t <1s 的范围内,感应电动势随电阻率的增大而单调下降。

大于1s 出现锯齿是因为概率积分原因所致。

目前最先进的观测系统V8的有效时间段为0108～25ms ,在该时间段内感应电动势随电阻率的增大而单调下降。

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月图1　感应电动势随电阻率的变化 图2　概率积分精度对计算结果的影响是说当分割20次时,电阻率的误差不超过0101Ω・m 。

这表明通过20次计算,计算出的视电阻率的精度就可以达到0101Ω・m 。

即使实测感应电动势曲线在时间轴上有10000个点,总计算量也不超过20×104次,对现在的计算机来说,这样的计算量很小。

虽然20×104次的计算量很小,但每一次都要进行正演计算。

因此,正演计算决定了二分搜索算法的性能。

而正演的关键是概率积分,采用梯度法计算概率积分时,为了控制计算的精度,采用步长收缩10倍来提高数值积分计算精度。

控制参・74・第31卷第2期陈清礼:瞬变电磁法全区视电阻率的二分搜索算法　数goal 表示2次相邻计算的相对误差,经过大量计算表明,该参数应小于10-6;否则,正演计算得到的电阻率不可靠,表现为感应电动势随电阻率非单调下降。

如图2所示,分别是goal 取010001、0100001和01000001时的情况,前2种情况下,曲线不符合单调递减的特性,与MA TLAB 程序计算结果不一致,明显不正确;而后一种情况则正确。

概率积分是最耗时间的部分,因为对概率积分的要求非常高,所以积分步长非常小,计算量因而很大,整个计算量的绝大部分都在此。

实际计算表明,采用梯度法数值方法计算概率积分时,计算一个测点的全区视电阻率需要的时间约4h 。

为了改进概率积分,研究采用了一种改进的算法,该算法在保障计算精度的前提下,可以把时间由4h 提高到1min ,极大地提高了计算的效率。

3　理论数据的计算理论数据因为有标准的参考结果,最能检验算法的正确性和精度,因此这里只用理论数据来证实该算法的可靠性。

层状介质模型与均匀半空间计算视电阻率的公式都是式(1),只需要讨论均匀半空间的情形即可。