大定源回线瞬变电磁场数值滤波算法

刘　云， 宋　滔， 王　赟

(中国科学院　地球化学研究所　矿床地球化学国家重点实验室，贵阳　550081)



大定源回线瞬变电磁场数值滤波算法

刘云， 宋滔*， 王赟

(中国科学院地球化学研究所矿床地球化学国家重点实验室，贵阳550081)

在前人的工作基础上，给出一种新的快速计算框线源激发的瞬变电磁场的数值滤波算法，其中内层积分的Hankel积分式采用47点J1型线性滤波系数，外层积分采用“n点式”Gauss-Legendre数值积分法。通过线性迭加技术，将多个单线源组合成任意多边形的瞬变电磁的发射框源。算例分析表明，数值计算结果和解析解的最大相对误差小于0.007%，单测点、单时间道的计算耗时为0.1 s。该方法可应用于瞬变电磁法生产实际中确定出最小延迟时间，以避免“框线效应”对瞬变电磁场的畸变影响。

大定源回线； 瞬变电磁场； 数值滤波算法； 框线效应

0　引言

大定源回线瞬变电磁法(Fixed Loop TEM)是时间域电磁法勘探中重要的测量装置，它是在地面上布置一个边长为几百米甚至上千米的回线框，以阶跃波或其他波形电流场源为激励源，在断电瞬间观测地下介质产生的二次感应电磁场随时间变化的衰减特性，从而探测出地下介质的电性结构分布[1-2]。它具有工作效率高、勘探深度范围大、资料信噪比高等优点，已广泛应用于地质构造、金属矿、煤田等地球物理勘探领域[3-6]。但是由框线源引起的场源附近早期瞬变场的“框线效应”问题，给瞬变电磁法的资料处理和解释带来了极大困难[7]。因此，有必要研究大定源回线装置瞬变电磁法的场源模拟技术。

在Nabighian等[8]提出的有限长电流线源的瞬变电磁场“解析法”，以及丁艳飞等[9]的框源结构瞬变电磁场“分阶段法”的基础上，给出另外一种简单、快速、精度高的计算大定源回线激发瞬变电磁场的数值滤波算法，并分析场源附近早期瞬变场的“框线效应”现象，从而有效计算出生产实际中的最小延迟时间。

1　均匀半空间线电流源瞬变电磁场表达式

如图1(a)所示，在均匀半空间地电模型下，当测线垂直穿过线源的任意位置时， Nabighian等[8]给出了电流线源地表瞬变电磁场归一化感应电动势f(单位，v·m-2)的解析表达式。即：

当测线位于线源的中线时，如图1(b)所示，令线源的长度为2L，则式(1)可以进一步简化为[8]:

图1　电流线源瞬变电磁场Fig.1　Transient electromagnetic field of current line source(a)测线位于线源任何位置；(b)测线位于线源中线

(2)

因此，式(2)为式(1)特殊情况下的解析表达式，式(2)可以用来验证式(1)数值解的精度。

2　积分式的数值滤波计算方法

在式(1)中有两层积分式，分别为定积分式和J1型Hankel积分式，需要用数值滤波方法计算。

2.1J1型Hankel积分式的数值滤波计算方法

在式(1)中，

(3)

称为J1型Hankel积分式，数值滤波计算公式为[11-12]式(4)。

i=1,2,…,n

(4)

其中：a为初始步长；s为位移系数。

Guptasarma等[13]及阮百尧[11]给出了47点、140点J1型Hankel数值滤波系数。通过与解析解验算[14-16]，这两套数值滤波系数的计算结果相对误差极小，均能满足工程计算中的精度要求。表1给出47点J1型Hankel数值滤波系数。

表1　47点J1型Hankel数值滤波系数

2.2定积分式的数值滤波计算方法

在式(1)中，要计算定积分式(式(5))，

(5)

常规的数值积分法为“5点式”或“10点式”Gauss数值积分法[17]。由于在野外实际中，线源尺寸往往会比较大，“5点式”或“10点式”Gauss数值积分法很难达到较高的计算精度，这里采用“n点式”Gauss-Legendre数值求积法[10]，其数值滤波计算公式见式(6)。

(6)

其中：Wi为n个数值滤波系数，可以根据线源长度任意设定n值大小。选取的n值越大，计算精度较高，则计算耗时较长。具体详见参考文献[10]中的数值计算理论和计算程序。

3　大定源回线的场源结构

图2为瞬变电磁法大定源回线的矩形场源结构，分别是由4个单线源组合成的矩形回线框源。在式(1)单线源的理论基础上，可以组合得到四边形大定源回线瞬变电磁场的归一化感应电动势为[18]

f=f1+f2+f3+f4。

(7)

同理可推，任意多边形大定源回线的场源结构也同样满足这样的构成规律。

图2　大定源回线Fig.2　The fixed loop(a)框内测量方式；(b)框外测量方式

4　算例分析

为验证方法的可靠性和有效性，设计两个理论模型进行数值模拟和分析。

4.1数值解的精度分析

地电模型如图3所示，均匀半空间电阻率为100 Ω·m，电流线源长度为200 m，测线位于线源的中线，测线x为-2 000 m～2 000 m，发射电流强度为1 A，延迟时间为0.3 ms。

图3　均匀半空间单线源地电模型Fig.3　The homogeneous half space geoelectric model with single line source

由式(1)计算得到的数值解为f1，式(2)计算得到的解析解为f2，令

(8)

其中，RE为数值解相对于解析解的百分比误差。

计算单个延迟时间的计算机耗时为0.1 s，模拟结果如图4所示。在整条测线上测点的数值解和解析解最大相对误差为0.006 6%，其余测点的相对误差均小于0.003%，两个误差极值点位于二次涡流向外扩散位置处。计算表明，这里所采用的数值滤波算法具有较高的模拟精度和计算速度。

图4　单线源瞬变电磁场数值解和解析解比较图Fig.4　Comparison of numerical and analytical solution of transient electromagnetic field by single line source(a)0.3 ms归一化感应电动势；(b)数值解相对解析解的百分误差

4.2大定源回线瞬变电磁场分析

地电模型图5，均匀半空间电阻率为100 Ω·m，测量区域x为-2 000 m～2 000 m，y为-2 000 m～2 000 m。点距为10 m，大回线发射框尺寸为1 000×1 000 m2，发射电流强度为1 A，延迟时间为1 μs～1 s，以10为底对数间隔采样，共计61个时间道。

图6为1 μs时瞬变电磁场归一化感应电动势模拟结果图。从图6可知，在早期的瞬变电磁场仍然集中于框源附近，其场的分布形状与框源形状相似。框源附近的早期瞬变电磁场的畸变较大，即所谓的“框线效应”，而处于框源中心的瞬变电磁场比较均匀，“框线效应”影响较弱。框源附近早期瞬变电磁场的“框线效应”对测量数据的影响较大，给数据处理和资料解释带来极大困难。

在生产实际中，为了有效避免“框线效应”影响，可以采用大定源回线瞬变电磁场数值模拟的办法来选取实际测量中的最小延迟时间。图7为选取y=0测线、框内数据(-500 m～500 m)的多时间道归一化感应电动势组成的综合剖面图。

图5　均匀半空间大定源回线地电模型Fig.5　The homogeneous half space geoelectric model with fixed loop

图6　1 μs瞬变电磁场模拟结果Fig.6　Transient electromagnetic field modeling results at 1 μs(a)1 μs归一化感应电动势水平剖面图;(b)y=0测线1 μs归一化感应电动势曲线图

图7　多时间道归一化感应电动势综合剖面图Fig.7　Multi-time gates normalized induced electromotive force comprehensive profile

显然，从1×10-3ms～7.9×10-2ms的综合剖面图呈“U”字形分布形态，表明“框线效应”对数据的畸变影响比较大，这一延迟时间范围内的数据不能用于反演和资料解释；从2.5×10-1ms后，“框线效应”影响消失，曲线形态正常。因此，选择的最小延迟时间，即第一个时间道为2.5×10-1ms为宜。同理可知，当同时存在框内、外测点数据时， 选取的最小延迟时间要比较框内测量时选取的最小延迟时间大。这说明当大定源回线瞬变电磁法采用尺寸较大的发射框源时，一部分的早期瞬变电磁数据由于受到“框线效应”的畸变影响，不得不采取这种剔除数据的办法，这对于地表的浅部信息反映是一个较大损失，这也是作者方法在实际应用方面的参考和局限之处。

5　结论

给出了一种新的计算均匀半空间大定源回线瞬变电磁场的数值滤波算法，分别采用了“n点式”Gauss-Legendre数值求积法，以及47点J1型Hankel变换线性滤波器数值求积法，单线源组合成大定源回线的线性迭加技术。与常规的算法相比较，大大简化了计算步骤，具有较高计算速度和数值模拟精度，使编程更容易实现。通过对大定源回线瞬变电磁场的数值模拟，可以确定出生产实际中最小延迟时间，从而有效避免瞬变电磁法的“框线效应”现象。

[1]米萨克 N．纳比吉安.勘查地球物理电磁法(第一卷)理论[M]．北京：地质出版社，1992．

MISAC N．NABIGHIAN. Electromagnetic Method in Applied Geophysics (Volume 1), Theory [M].Beijing: Geological Publishing House, 1992.(In Chinese)

[2]蒋邦远．实用近区磁源瞬变电磁勘探[M]．北京：地质出版社，1998．

JIANG B Y. Near field magnetic source transient electromagnetic exploration [M].Beijing: Geological Publishing House, 1998. (In Chinese)

[3]陈贵生．瞬变电磁法在金属矿产勘查上的应用效果及存在问题探讨[J]．矿产与地质，2006，20(4)：543—547．

CHEN G S.Application result of transient electromagnetic method to metallic mineral resources exploration and existing problems [J].Mineral resources and geology, 2006, 20(4):543-547. (In Chinese)

[4]薛国强，李貅．瞬变电磁隧道超前预报成像技术[J]．地球物理学报，2008，51(3)：894—900．

XUE G Q, LI X. The technology of TEM tunnel prediction imaging [J].Chinese J. Geophysics,2008,51(3):894-900.(In Chinese)

[5]李宁，张文博，刘域田，等．瞬变电磁法在红旗岭铜镍矿3号岩体勘查中的应用与研究[J]．吉林地质，2014，33(2)：69—72．

LI N, ZHANG W B, LIU Y T, et al. Application and research of the transient electromagnetic method in No.3 Rock mass of Hongqiling copper-nickle mine [J].Jilin Geology, 2014, 33(2): 69-72. (In Chinese)

[6]韩自强，罗姣，刘涛，等．定源回线瞬变电磁法在煤矿富水区调查中的应用[J]．地球物理学进展，2015，30(4)：1705—1711．

HAN Z Q, LUO J, LIU T, et al. The application of fixed source loop transient electromagnetic method in the coal mine rich water area survey [J].Progress in Geophysics, 2015, 30(4): 1705-1711. (In Chinese)

[7]包乃利，刘鸿福，余传涛．大定源瞬变电磁法激励场及边框效应研究[J]．煤田地质与勘探，2014,42(2)：80—84．

BAO N L, LIU H F, YU C T. The research on the primary field and rim effect of transient electromagnetic method with large fixed loop [J].Coal Geology & Exploration, 2014, 42(2): 80-84. (In Chinese)

[8]NABIGHIAN M. N，ORISTAGLIO M. L. On the approximation of finite loop sources by two-dimensional line source [J]．Geophysics, 1984, 49(7): 1027—1029．

[9]丁艳飞,白登海,许诚．均匀半空间表面大定源瞬变电磁响应的快速算法[J]．地球物理学报,2012,55(6):2087—2096．

DING Y F, BAI D H, XU C. A rapid algorithm for calculating time domain transient electromagnetic responses of a large fixed loop on the half space [J]. Chinese J. Geophysics, 2012, 55(6):2087-2096. (In Chinese)

[10]何光渝,高永利．Visual Fortran 常用算法集[M]．北京:科学出版社,2002．

HE G Y, GAO Y L. Visual Fortran common algorithm set [M].Beijing: Science Press, 2002. (In Chinese)

[11]阮百尧．均匀水平大地上频率域垂直磁偶源电磁场数值滤波解法[J]．桂林工学院学报,2005,25(1):14—18．

RUAN B Y. Digital filter method of evaluating electromagnetic field from a vertical magnetic dipole above the homogeneous earth [J]. Journal of Guilin University of Technology, 2005, 25(1):14-18. (In Chinese)

[12]程志平．电法勘探教程[M]．北京:冶金工业出版社,2007．

CHENG Z P. Electric exploration course [M].Beijing: Metallurgical industry press, 2007. (In Chinese)

[13]GUPTASARMA D, SINGH B. New digital linear filters for Hankel J0a and J1 transforms [J]．Geophysical Prospecting, 1997,45: 745—762．

[14]ANDERSON W. L.Computer Program Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering[J]．Geophysics, 1979,44(7): 1287—1305．

[15]ALAN D. C.Numerical integration of related Hankel transforms by quadrature and continue fraction expansion[J]．Geophysics, 1983,48(12): 1671—1686．

[16]张伟，王绪本，覃庆炎．汉克尔变换的数值计算与精度的对比[J]．物探与化探,2010,34(6)：753—755．

ZHANG W, WANG X B, QIN Q Y. Research and Application on Numerical Integration of Hankel Transform by Digital Filtering[J]．Geophysical & Geochemical Exploration, 2010,34(6): 753-755．(In Chinese)

[17]徐世良．计算机常用算法(第二版)[M]．北京: 清华大学出版社，1995．

XU S L. Computer algorithms (Second Edition) [M].Beijing: Tsinghua university press, 1995. (In Chinese)

[18]刘云，王绪本，段长生．大回线瞬变电磁正演模拟在工程实际中的应用[C]．第十届中国国际地球电磁学术讨论会论文集，2011:50—54．

LIU Y, WANG X B, DUAN C S. Practical Application of Fixed-loop TEM Forward Modeling[C].The 10thChina International Geo-Electromagnetic Workshop, 2011:50-54. (In Chinese)

A digital filter method for evaluating the fixed-loop transient electromagnetic field

LIU Yun, SONG Tao*, WANG Yun

(State Key Laboratory of Ore Deposit Geochemistry,Institute of Geochemistry Chinese Academy of Sciences,Guiyang550081, China)

Based on the studies of Nabighian (1984) and DING Yan-fei (2012), a new digital filtering algorithm is presented to evaluate the transient electromagnetic (TEM) field induced by a line source. Among them, the inner integral Hankel integral formula using 47-point J1 type linear filtering coefficients, outer integral formula using 'n-point' Gauss-Legendre numerical integration method. Through linear superposition of a plurality of single line source emission source combined into arbitrary polygon transient electromagnetic. Example analysis shows that the maximum relative error of the numerical results and the analytical solution is less than 0.007%, and the computation time of single point and single time channel is 0.1s. The method can be applied to the production of transient electromagnetic method to determine the minimum delay time, so as to avoid the influence of border effects on the transient electromagnetic field.

fixed-loop; transient electromagnetic field; digital filtering algorithm; border effect

2016-04-23改回日期：2016-05-20

国家“973计划”项目(2014CB440905)；国家自然科学基金(41440031)；贵州省科学技术基金(2014GZ93278)

刘云(1973-)，男，博士，副研究员，研究方向为地球物理数值模拟与反演成像，E-mail： liu_yun@mail.gyig.ac.cn。

宋滔(1988-)，男，博士，研究方向为地球物理数值模拟及与演成像，E-mail: shaman_king7@163.com。

1001-1749(2016)04-0437-06

P 631.3

A

10.3969/j.issn.1001-1749.2016.04.01