2. 中国科学院研究生院, 北京 100049
2. Graduate University, Chinese Academy of Sciences, Beijing 100049, China
瞬变电磁法以其探测深度广、分辨率高、装置轻便、工作效率高等优点广泛应用于资源勘探与工程勘查领域[1].其基本原理是,导电介质在一次电流脉冲场激励下会产生涡流,在发射电流间断期间涡流不会立即消失,在其周围空间形成随时间衰减的二次磁场.二次磁场随时间衰减的规律主要取决于介质的导电性、规模、埋深,以及发射电流的形态和频率.因此,我们可以通过接收线圈测量的二次场的衰减特性了解地下的电性结构[2].
利用时间域瞬变电磁法(TDEM/TEM)进行大深度面积探测时最常用的一种观测方式是大定源(大回线源)方式[3-5].大回线源一般指边长达数百米甚至上千米的矩形或正方形发射回线,测量时回线源固定不动,接收装置在回线内外不同位置移动.这种装置的优点是可以使用较大的发射功率,抗干扰能力强,信噪比高,探测深度大.另外,这种装置一旦铺好回线后,可采用多台接收机同时工作[6],或者把接收线圈排成阵列,形成阵列式观测系统,适合于大范围的面积测量,效率高.这种方式的缺点是瞬变场响应的空间分布复杂,数据处理比较困难.测点偏离发射框中心越远,计算速度越慢.这对于大面积测量的海量数据来说是无法接受的,所以目前大多数处理方法是采用中心方式的晚期近似,误差比较大.尤其是早期和中期信息不能得到有效利用,影响了浅层探测效果.计算速度一直是影响TEM 大定源方式全程精确解的瓶颈,本文的目的就是希望寻找一种可行的办法,既能保证计算精度,又能提高计算效率.
计算时间域瞬变电磁响应常用的手段是频-时转换[7].最早的转换方法是余弦多项式近似算法[2, 8-9],该方法要求对核函数抽样次数太多,计算速度很慢[7].20世纪80年代之后,Gaver-Stehfest逆拉普拉斯变换法[10-12]被用于TDEM 响应的计算,但该方法计算量大、计算速度仍然较慢[7].此外,还有Guptasarma算法(数字滤波法)[13],余弦变换数字滤波法[14],网格嵌套插值法[12]及有限项的级数展开方法[15-19]等.数字滤波法对应的装置模型简单,适用性不强.本文提出一种快速计算大定源方式时间域瞬变电磁全程响应的算法,可有效提高计算效率.该方法将全程瞬变过程分为早中晚三个阶段,对于早期和晚期阶段,将积分核函数进行渐近展开[20];在中期阶段,对内层积分(即误差函数erf)进行有理Chebyshev 渐近展开[21-22],对外层定积分采用Romberg 积分[23].实际计算表明,该方法在保证计算精度的同时,可大大提高计算速度.
2 理论基础设一个矩形发射线框置于均匀半空间表面,如图 1所示,采用直角坐标系,坐标原点位于矩形线框的中心,Z轴垂直于发射框所在平面指向内.矩形框的长、宽分别为2a和2b.我们的目的就是计算当发射框中电流突然中断以后地面任意观测点的二次瞬变响应场.
设发射框中载有直流电I,当供电电流中断后,矩形线框在均匀半空间表面产生的时间域二次瞬变响应场(bz及其时间变化率$\frac{\partial {{b}_{z}}}{\partial t}$)的表达式为[3-6]:
(1) |
和
(2) |
其中μ0 =4π×10-7h/m 为真空中的磁导率,I为回线中所载的电流,t为时间,从电流关断时算起.A′、B′、C′和D′分别为四个边所产生的磁场的积分函数,A″、B″、C″和D″分别为四个边所产生的磁场变化率的积分函数,其具体表达式为
其中,
(3) |
其中
将(1)和(2)式分别改写为
其中
(4) |
(5) |
Y(Z)和Y′(Z)分别为二次瞬变磁场垂直分量bz及其时间变化率∂bz/∂t的归一化响应核函数.计算瞬变响应实际上就是计算(4)和(5)式所表达的核函数.
3 核函数的计算方法核函数Y(Z)和Y′(Z)没有简单的解析解,只能通过数值手段求近似解.下面把整个瞬变过程分为早期、中期和晚期三个阶段进行讨论.
3.1 早期和晚期阶段的数值算法核函数的数值表达式可由误差函数erf(w)的级数展开式得到.erf(w)有两种展开形式,一种为晚期表达式(对应小w),另一种为早期表达式(对应大w).在实际情况中,空间变量r总是有限的,所以由(3)式可知,w的极限性质与(ρt)-1/2 和Z的极限性质保持一致.
3.1.1 核函数的晚期表达式当w较小时(即ρt较大,Z较小时)[18]:
(6) |
(7) |
将(3)、(6)和(7)式分别代入(4)和(5)式,求解积分后便可得到Y(Z)和Y′(Z)的晚期渐近式:
(8) |
(9) |
其中ζ(i)和γ(i)为多项式系数,表 1给出了前40项的值.φ(i)=-φ1 -φ2 +φ3 +φ4,
级数(6)和(7)虽然是在小变量条件下的展开式,但(8)和(9)式的收敛域为(- ∞,+ ∞).从理论上讲,对于一个已知的Z值,只要(8)和(9)式中求和数i足够大,便可以求得比较精确的Y(Z)和Y′(Z)值.但是,计算机的计算能力总是有限的,计算(8)和(9)式时i只能取有限值,否则就会发生溢出错误.另外,当i值较大时,计算耗时将大幅度增加.经试验验证,实际计算时i取有限项便可满足速度和精度要求((8)和(9)式的最佳截断阶数和使用区间参照表 3).在晚期阶段,由于二次涡流场已经扩散到地下较远的地方,对于地面测点来说,相当于一个偶极子产生的场,所以不论测点位于发射框内还是框外,其表达式都是相同的.
当w很大(即ρt很小,Z很大时)时,误差函数可以表示为[24-25]:
(10) |
(11) |
其中
将(10)式分别代入(4)和(5)式,逐项积分后得到Y(Z)和Y′(Z)的早期渐近式.在早期阶段,二次涡流场位于发射框附近,其表现特性严重依赖于发射框的形状和测点相对于发射框的位置.为了便于理解和计算,我们把观测区域分为三类(见图 2),早期场在三类区域的表现特性具有很大差别.
I:测点位于发射框内(|x| <a,|y|<b)
(12) |
(13) |
χ(i)和ξ(i)为多项式系数,表 2给出了χ(i)和ξ(i)前4项的值.
II:测点位于发射框外非对角区域(|x|>a,|y|<b;|x|<a,|y|>b)
① 当x>a,|y|<b时:
(14) |
(15) |
当x<-a,|y|<b时,(14)和(15)式中令x= |x|即可.
② 当|x|<a,y>b时:
(16) |
(17) |
当|x|<a,y<-b时,(16)和(17)式中令y= |y|即可.
III:测点位于发射框外对角区域(|x|>a,|y|>b)
当x>a,y>b时:
(18) |
(19) |
其他情况下,(18)和(19)式中只需令x=|x|,y= |y|即可.
(11) 和(12)式的收敛条件为(具体计算方法见附录A).
(20) |
其中(x,y)为测点坐标.当Z小于该临界值时(12)和(13)式将出现发散(或震荡)现象.例如,当矩形框为100m×100m, 测点坐标为(20,0),i=2时,Z>2.208.实际上,(12)和(13)式求和只需要取前边有限几项就可以满足精度要求((12)和(13)式的最佳截断阶数和使用区间参照表 3).
3.2 中期阶段的数值算法对于时间域瞬变电磁全程响应的计算,晚期和早期都相对比较容易处理,最难的是对二者之间的过渡阶段(即中期阶段)的处理.当测点与框边之间的距离大于边长的25%时,早期曲线和晚期曲线可以很好地衔接起来.然而,随着测点与框边之间距离的缩小,晚期表达式的收敛速度明显减慢,早期表达式的收敛半径亦随之减小.增加展开式的项数,可以实现早期曲线和晚期曲线的衔接,但是计算所需的CPU 时间将大幅度增加.因此,在中期阶段时,我们采用数值积分方法来求解.
(4) 和(5)式的内层积分为误差函数erf(w),计算时采用Cody[21-22]提出的有理Chebyshev渐近展开式:
其中pi和qi是多项式系数,具体数值见表 4和表 5(这里虽然只给出了前几项的值,但足够用于实际计算).
采用该方法,可将误差函数erf(w)的计算误差控制在6×10-19以内[21-22],并且能够大幅度提高计算速度.对于外层的定积分则采用常规的Romberg数值积分法[23]进行计算.
4 理论模型验证 4.1 近似解的误差评估设矩形发射框置于10Ωm 的水平均匀半空间表面,线框大小为400m×400m, 发送电流为1A.
由于瞬变场的对称性,在进行理论模型验证时,我们只考虑图 1中第一象限的测点.测线选择两个不同方向(见图 3).测线1通过发射框中心沿Y轴覆盖Ⅰ类和Ⅱ类区域,测线2通过发射框中心和一个顶点,覆盖Ⅰ类和Ⅲ类区域.
设P1 和Q1 分别为核函数bz和∂bz/∂t的精确解,精确解由常规的积分方法求得.P2 和Q2 分别为bz和∂bz/∂t根据本文算法求得的近似解,令
其中,RDb 为bz的近似解相对于其精确解的百分比误差,RDd 为∂bz/∂t的近似解相对于其精确解的百分比误差.求解P1 和Q1 时,(8)、(9)、(12)和(13)式的展开项数和使用区间参照表 3.
图 4 给出了两条测线上不同测点处∂bz/∂t的近似解与精确解的相对误差随时间的变化,图 5给出了两条测线上不同测点处bz的近似解与精确解的相对误差随时间的变化.从图 4,5可以看到,在计算的整个时间区间上(一般的TEM 设备的记录时间范围为几微妙(nμs)到几百毫秒),无论是∂bz/∂t还是bz,其近似曲线与精确解的相对误差非常小.对于∂bz/∂t来说,不同测点的误差可以控制在0.0001%以内;对于bz来说,不同测点的误差可以控制在0.0002% 以内.从图 4,5 还可以看到:(1)晚期渐近值与精确解之间的误差均小于2×10-8%.实际上,随着晚期展开项数的增加,晚期渐近值与精确解之间的误差会进一步减小.但是考虑到计算速度,在进行数值计算时,只展开了前40项;(2)早期阶段近似值与精确解的相对误差较大,在1×10-7%附近徘徊,但已经满足计算的精度要求.当取展开式的前几项时,即可得到很好的渐近效果.由于早期阶段存在震荡现象,随着早期渐近式展开项数的增加,渐近曲线并没有进一步逼近精确值,反而出现了衰退现象.这就涉及到了渐近级数的最佳阶数问题;(3)中期阶段的误差在1×10-7%~0.0002%之间.在计算式(4)和(5)的外层积分时,采用Romberg数值积分法,截断误差设为1×10-6.实际上,随着所选误差值的减小,整体误差也将减小,但计算时间却会成倍增长;(4)随着测点到矩形框中心距离的增加,近似值与精确值之间的误差突变点在向右推移.这是因为,随着测点到矩形框中心距离的增加,晚期渐近级数的收敛速度逐渐减慢造成的.
在其它条件相同的情况下(即10Ωm 均匀半空间,400m×400m 发射框,同一台计算机.其中,所用计算机配置参数如下:CPU 为AMD AthlonX24400+,双核处理器,主频为2.3 GHz, 内存为1GB,使用年限为4年),分别采用常规算法和本文的近似算法计算不同位置的核函数bz和∂bz/∂t,然后比较两种算法所花费的CUP 时间(表 6和表 7).实际测试表明,两种算法在不同测点所花费的时间都不同,这与理论预测是一致的.在发射框中心,由于对称性,当测点距离边框的距离大于边长的25%时,单点计算所需的CPU 时间可以控制在1s以内,比常规计算快7倍;其他测点处单点计算所需的CPU 时间亦可控制在3s以内,比常规计算快4倍.
常规算法:晚期阶段,对(4)式和(5)式中的内层积分采用(6)式,早期阶段采用(11)式。对(4)式和(5)式中的外层积分采用传统的辛普森积分方法.
5 结论对置于均匀导电半空间表面上的大回线框所产生的瞬变磁场响应,本文给出了一种新的数值计算方法.该算法的思想是将核函数分为早期、中期和晚期三个阶段,对早晚期采用渐近展开的方法求解积分,对于中期采用有理展开和数值积分求解相结合的方法.理论模型证明与传统算法相比,该算法在保证精度的同时,大大提高了计算速度.
(1) 当测点位于回线线框内部且测点距框边的距离大于边长的25%时,瞬变响应的渐近解同精确解的相对误差小于0.0001%,并且计算每个点所需的CPU 时间可以控制在1s以内,比常规计算快7倍.
(2) 当测点位于距离框边小于边长的25% 时,瞬变响应的渐近解同精确解的相对计算误差可以控制在0.0002%以内,且计算每个点的CPU 时间亦可控制在3s以内,比常规计算快4倍.
(3) 本文所提出算法的误差极大值存在于早期、中期和晚期曲线结合点的位置上.因此,如果找出一种更合理的方法来选取早期、中期和晚期曲线的结合点,便可以进一步提高该算法的计算精度;对于核函数的展开,针对早期和晚期阶段,我们采取普通的多项式展开.而多项式展开普遍存在的问题是收敛速度和收敛域的问题,如果我们能够找到它们的有理函数展开,有可能能够进一步加快算法的计算速度.
附录A Y(Z)和Y′(Z)早期渐近式收敛区间由(10)式收敛的条件为
(A1) |
经整理:
(A2) |
由(3)式可得
(A3) |
其中ri为观测点距离四个边框的距离.
由于
(A4) |
(A5) |
(A6) |
(A7) |
由均值不等式可得:
所以
(A8) |
由于S=4ab,代入(A8)式即可得到(20)式.
致谢感谢UmitAvsar博士在TEM 理论和绘图方面的指导和帮助.
[1] | 薛国强, 李貅, 宋建平, 等. 回线源瞬变电磁成像的理论分析及数值计算. 地球物理学报 , 2004, 47(2): 338–343. Xue G Q, Li X, Song J P, et al. Theoretical analysis and numerical calculation of loop-source transient electromagnetic imaging. Chinese J. Geophys. (in Chinese) , 2004, 47(2): 338-343. |
[2] | Kaufman A A, Keller G V. Frequency and Transient Sounding. New York: Elsevier Science Publishers B V , 1983: 1-38. |
[3] | Poddar M. A rectangular loop source of current on a two-layered earth. Geophysical Prospecting , 1982, 30(1): 101-114. DOI:10.1111/gpr.1982.30.issue-1 |
[4] | Rao K P, Saraf P D, Mallick K. Electromagnetic depth sounding on a transitional earth using large rectangular loop source // Deep Electromagnetic Exploration. Lecture Notes in Earth Sciences. New Delhi, India: Narosa Publishing House, 1998, 83: 628-638. |
[5] | Rao K P, Babu G A. Microsoft C. NET program and electromagnetic depth sounding for large loop source. Computers & Geosciences , 2009, 35(7): 1369-1378. |
[6] | Poddar M. A rectangular loop source of current on multilayered earth. Geophysics , 1983, 48(1): 107-109. DOI:10.1190/1.1441398 |
[7] | 陈向斌, 胡社荣, 张超. 瞬变电磁场响应计算的频-时域转换方法综述. 工程地球物理学报 , 2008, 5(2): 242–246. Chen X B, Hu S R, Zhang C. Review of frequency-domain and time-domain transformation method in transient electromagnetic field response computation. Chinese Journal of Engineering Geophysics (in Chinese) , 2008, 5(2): 242-246. |
[8] | 朴化荣. 电磁测深法原理. 北京: 地质出版社, 1990 : 112 -149. Piao H R. Principles of Electromagnetic Sounding (in Chinese). Beijing: Geological Publishing House, 1990 : 112 -149. |
[9] | 翁爱华, 王雪秋. 长偏移距瞬变电磁测深甚晚期响应及视电阻率的数值计算. 地震地质 , 2003, 25(4): 664–668. Weng A H, Wang X Q. Numerical simulations of very-late time response and apparent resistivity in long-offset TEM sounding. Seismology and Geology (in Chinese) , 2003, 25(4): 664-668. |
[10] | Knight J H, Raiche A P. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method. Geophysics , 1982, 47(1): 47-50. DOI:10.1190/1.1441280 |
[11] | 罗延钟, 昌彦君. G-S变换的快速算法. 地球物理学报 , 2000, 43(5): 684–690. Luo Y Z, Chang Y J. A rapid algorithm for G-S transform. Chinese J. Geophys. (in Chinese) , 2000, 43(5): 684-690. |
[12] | Raiche A P. Transient electromagnetic field computations for polygonal loops on layered earths. Geophysics , 1987, 52(6): 785-793. DOI:10.1190/1.1442345 |
[13] | Guptasarma D. Computation of the time-domain response of a polarizable ground. Geophysics , 1982, 47(11): 1574-1576. DOI:10.1190/1.1441307 |
[14] | 王华军. 正余弦变换的数值滤波算法. 工程地球物理学报 , 2004, 1(4): 329–333. Wang H J. Digital filter algorithm of the sine and cosine transform. Chinese Journal of Engineering Geophysics (in Chinese) , 2004, 1(4): 329-333. |
[15] | Spies B R, Raiche A P. Calculation of apparent conductivity for the transient electromagnetic (coincident loop) method using an HP-67 calculator. Geophysics , 1980, 45(7): 1197-1204. DOI:10.1190/1.1441117 |
[16] | Raiche A P, Spies B R. Coincident loop transient electromagnetic master curves for interpretation of two-layer earths. Geophysics , 1981, 46(1): 53-64. DOI:10.1190/1.1441139 |
[17] | Hasegawa K. On the step response and the apparent conductivity for a stratified earth by a horizontal electric dipole. Geophysical Exploration , 1985, 38(3): 21-31. |
[18] | 白登海, MejuM A, 卢健, 等. 时间域瞬变电磁法中心方式全程视电阻率的数值计算. 地球物理学报 , 2003, 46(5): 697–704. Bai D H, Meju M A, Lu J, et al. Numerical calculation of all-time apparent resistivity for the central loop transient electromagnetic method. Chinese J. Geophys. (in Chinese) , 2003, 46(5): 697-704. |
[19] | 蒋邦远. 实用近区磁源瞬变电磁法勘探. 北京: 地质出版社, 1998 : 1 -60. Jiang B Y. Near Zone Magnetic Source Transient EM Prospect (in Chinese). Beijing: Geological Publishing House, 1998 : 1 -60. |
[20] | 李家春, 周显初. 数学物理中的渐近方法. 北京: 科学出版社, 2002 . Li J C, Zhou X C. Asymptotic Methods in Mathematical Physics (in Chinese). Beijing: Science Press, 2002 . |
[21] | Cody W J. Rational chebyshev approximations for the error function. Mathematics of Computation , 1969, 23(107): 631-637. DOI:10.1090/S0025-5718-1969-0247736-4 |
[22] | Cody W J. Performance evaluation of programs for the error and complementary error functions. ACM Transactions on Mathematical Software , 1990, 16(1): 29-37. DOI:10.1145/77626.77628 |
[23] | Burden R L, Faires J D. Numerical Analysis (Seventh Edition). Beijing: Higher Education Press, 2003 . |
[24] | Dingle R B. Asymptotic Expansions: Their Derivation and Interpretation. London: Academic Press, 1973 . |
[25] | Abramowitz M, Stegun I A. Handbook of Mathematical Functions. The USA: Dover Publications , 1964: 295-308. |