大定源回線瞬變電磁場數值濾波算法

劉　云， 宋　滔， 王　赟

(中國科學院　地球化學研究所　礦床地球化學國家重點實驗室，貴陽　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