大定源回線瞬變電磁場(chǎng)數(shù)值濾波算法

劉　云， 宋　滔， 王　赟

(中國(guó)科學(xué)院　地球化學(xué)研究所　礦床地球化學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，貴陽(yáng)　550081)

?

大定源回線瞬變電磁場(chǎng)數(shù)值濾波算法

劉云， 宋滔*， 王赟

(中國(guó)科學(xué)院地球化學(xué)研究所礦床地球化學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，貴陽(yáng)550081)

在前人的工作基礎(chǔ)上，給出一種新的快速計(jì)算框線源激發(fā)的瞬變電磁場(chǎng)的數(shù)值濾波算法，其中內(nèi)層積分的Hankel積分式采用47點(diǎn)J1型線性濾波系數(shù)，外層積分采用“n點(diǎn)式”Gauss-Legendre數(shù)值積分法。通過(guò)線性迭加技術(shù)，將多個(gè)單線源組合成任意多邊形的瞬變電磁的發(fā)射框源。算例分析表明，數(shù)值計(jì)算結(jié)果和解析解的最大相對(duì)誤差小于0.007%，單測(cè)點(diǎn)、單時(shí)間道的計(jì)算耗時(shí)為0.1 s。該方法可應(yīng)用于瞬變電磁法生產(chǎn)實(shí)際中確定出最小延遲時(shí)間，以避免“框線效應(yīng)”對(duì)瞬變電磁場(chǎng)的畸變影響。

大定源回線； 瞬變電磁場(chǎng)； 數(shù)值濾波算法； 框線效應(yīng)

0　引言

大定源回線瞬變電磁法(Fixed Loop TEM)是時(shí)間域電磁法勘探中重要的測(cè)量裝置，它是在地面上布置一個(gè)邊長(zhǎng)為幾百米甚至上千米的回線框，以階躍波或其他波形電流場(chǎng)源為激勵(lì)源，在斷電瞬間觀測(cè)地下介質(zhì)產(chǎn)生的二次感應(yīng)電磁場(chǎng)隨時(shí)間變化的衰減特性，從而探測(cè)出地下介質(zhì)的電性結(jié)構(gòu)分布[1-2]。它具有工作效率高、勘探深度范圍大、資料信噪比高等優(yōu)點(diǎn)，已廣泛應(yīng)用于地質(zhì)構(gòu)造、金屬礦、煤田等地球物理勘探領(lǐng)域[3-6]。但是由框線源引起的場(chǎng)源附近早期瞬變場(chǎng)的“框線效應(yīng)”問(wèn)題，給瞬變電磁法的資料處理和解釋帶來(lái)了極大困難[7]。因此，有必要研究大定源回線裝置瞬變電磁法的場(chǎng)源模擬技術(shù)。

在Nabighian等[8]提出的有限長(zhǎng)電流線源的瞬變電磁場(chǎng)“解析法”，以及丁艷飛等[9]的框源結(jié)構(gòu)瞬變電磁場(chǎng)“分階段法”的基礎(chǔ)上，給出另外一種簡(jiǎn)單、快速、精度高的計(jì)算大定源回線激發(fā)瞬變電磁場(chǎng)的數(shù)值濾波算法，并分析場(chǎng)源附近早期瞬變場(chǎng)的“框線效應(yīng)”現(xiàn)象，從而有效計(jì)算出生產(chǎn)實(shí)際中的最小延遲時(shí)間。

1　均勻半空間線電流源瞬變電磁場(chǎng)表達(dá)式

如圖1(a)所示，在均勻半空間地電模型下，當(dāng)測(cè)線垂直穿過(guò)線源的任意位置時(shí)， Nabighian等[8]給出了電流線源地表瞬變電磁場(chǎng)歸一化感應(yīng)電動(dòng)勢(shì)f(單位，v·m-2)的解析表達(dá)式。即：

當(dāng)測(cè)線位于線源的中線時(shí)，如圖1(b)所示，令線源的長(zhǎng)度為2L，則式(1)可以進(jìn)一步簡(jiǎn)化為[8]:

圖1　電流線源瞬變電磁場(chǎng)Fig.1　Transient electromagnetic field of current line source(a)測(cè)線位于線源任何位置；(b)測(cè)線位于線源中線

(2)

因此，式(2)為式(1)特殊情況下的解析表達(dá)式，式(2)可以用來(lái)驗(yàn)證式(1)數(shù)值解的精度。

2　積分式的數(shù)值濾波計(jì)算方法

在式(1)中有兩層積分式，分別為定積分式和J1型Hankel積分式，需要用數(shù)值濾波方法計(jì)算。

2.1J1型Hankel積分式的數(shù)值濾波計(jì)算方法

在式(1)中，

(3)

稱為J1型Hankel積分式，數(shù)值濾波計(jì)算公式為[11-12]式(4)。

i=1,2,…,n

(4)

其中：a為初始步長(zhǎng)；s為位移系數(shù)。

Guptasarma等[13]及阮百堯[11]給出了47點(diǎn)、140點(diǎn)J1型Hankel數(shù)值濾波系數(shù)。通過(guò)與解析解驗(yàn)算[14-16]，這兩套數(shù)值濾波系數(shù)的計(jì)算結(jié)果相對(duì)誤差極小，均能滿足工程計(jì)算中的精度要求。表1給出47點(diǎn)J1型Hankel數(shù)值濾波系數(shù)。

表1　47點(diǎn)J1型Hankel數(shù)值濾波系數(shù)

2.2定積分式的數(shù)值濾波計(jì)算方法

在式(1)中，要計(jì)算定積分式(式(5))，

(5)

常規(guī)的數(shù)值積分法為“5點(diǎn)式”或“10點(diǎn)式”Gauss數(shù)值積分法[17]。由于在野外實(shí)際中，線源尺寸往往會(huì)比較大，“5點(diǎn)式”或“10點(diǎn)式”Gauss數(shù)值積分法很難達(dá)到較高的計(jì)算精度，這里采用“n點(diǎn)式”Gauss-Legendre數(shù)值求積法[10]，其數(shù)值濾波計(jì)算公式見(jiàn)式(6)。

(6)

其中：Wi為n個(gè)數(shù)值濾波系數(shù)，可以根據(jù)線源長(zhǎng)度任意設(shè)定n值大小。選取的n值越大，計(jì)算精度較高，則計(jì)算耗時(shí)較長(zhǎng)。具體詳見(jiàn)參考文獻(xiàn)[10]中的數(shù)值計(jì)算理論和計(jì)算程序。

3　大定源回線的場(chǎng)源結(jié)構(gòu)

圖2為瞬變電磁法大定源回線的矩形場(chǎng)源結(jié)構(gòu)，分別是由4個(gè)單線源組合成的矩形回線框源。在式(1)單線源的理論基礎(chǔ)上，可以組合得到四邊形大定源回線瞬變電磁場(chǎng)的歸一化感應(yīng)電動(dòng)勢(shì)為[18]

f=f1+f2+f3+f4。

(7)

同理可推，任意多邊形大定源回線的場(chǎng)源結(jié)構(gòu)也同樣滿足這樣的構(gòu)成規(guī)律。

圖2　大定源回線Fig.2　The fixed loop(a)框內(nèi)測(cè)量方式；(b)框外測(cè)量方式

4　算例分析

為驗(yàn)證方法的可靠性和有效性，設(shè)計(jì)兩個(gè)理論模型進(jìn)行數(shù)值模擬和分析。

4.1數(shù)值解的精度分析

地電模型如圖3所示，均勻半空間電阻率為100 Ω·m，電流線源長(zhǎng)度為200 m，測(cè)線位于線源的中線，測(cè)線x為-2 000 m～2 000 m，發(fā)射電流強(qiáng)度為1 A，延遲時(shí)間為0.3 ms。

圖3　均勻半空間單線源地電模型Fig.3　The homogeneous half space geoelectric model with single line source

由式(1)計(jì)算得到的數(shù)值解為f1，式(2)計(jì)算得到的解析解為f2，令

(8)

其中，RE為數(shù)值解相對(duì)于解析解的百分比誤差。

計(jì)算單個(gè)延遲時(shí)間的計(jì)算機(jī)耗時(shí)為0.1 s，模擬結(jié)果如圖4所示。在整條測(cè)線上測(cè)點(diǎn)的數(shù)值解和解析解最大相對(duì)誤差為0.006 6%，其余測(cè)點(diǎn)的相對(duì)誤差均小于0.003%，兩個(gè)誤差極值點(diǎn)位于二次渦流向外擴(kuò)散位置處。計(jì)算表明，這里所采用的數(shù)值濾波算法具有較高的模擬精度和計(jì)算速度。

圖4　單線源瞬變電磁場(chǎng)數(shù)值解和解析解比較圖Fig.4　Comparison of numerical and analytical solution of transient electromagnetic field by single line source(a)0.3 ms歸一化感應(yīng)電動(dòng)勢(shì)；(b)數(shù)值解相對(duì)解析解的百分誤差

4.2大定源回線瞬變電磁場(chǎng)分析

地電模型圖5，均勻半空間電阻率為100 Ω·m，測(cè)量區(qū)域x為-2 000 m～2 000 m，y為-2 000 m～2 000 m。點(diǎn)距為10 m，大回線發(fā)射框尺寸為1 000×1 000 m2，發(fā)射電流強(qiáng)度為1 A，延遲時(shí)間為1 μs～1 s，以10為底對(duì)數(shù)間隔采樣，共計(jì)61個(gè)時(shí)間道。

圖6為1 μs時(shí)瞬變電磁場(chǎng)歸一化感應(yīng)電動(dòng)勢(shì)模擬結(jié)果圖。從圖6可知，在早期的瞬變電磁場(chǎng)仍然集中于框源附近，其場(chǎng)的分布形狀與框源形狀相似。框源附近的早期瞬變電磁場(chǎng)的畸變較大，即所謂的“框線效應(yīng)”，而處于框源中心的瞬變電磁場(chǎng)比較均勻，“框線效應(yīng)”影響較弱。框源附近早期瞬變電磁場(chǎng)的“框線效應(yīng)”對(duì)測(cè)量數(shù)據(jù)的影響較大，給數(shù)據(jù)處理和資料解釋帶來(lái)極大困難。

在生產(chǎn)實(shí)際中，為了有效避免“框線效應(yīng)”影響，可以采用大定源回線瞬變電磁場(chǎng)數(shù)值模擬的辦法來(lái)選取實(shí)際測(cè)量中的最小延遲時(shí)間。圖7為選取y=0測(cè)線、框內(nèi)數(shù)據(jù)(-500 m～500 m)的多時(shí)間道歸一化感應(yīng)電動(dòng)勢(shì)組成的綜合剖面圖。

圖5　均勻半空間大定源回線地電模型Fig.5　The homogeneous half space geoelectric model with fixed loop

圖6　1 μs瞬變電磁場(chǎng)模擬結(jié)果Fig.6　Transient electromagnetic field modeling results at 1 μs(a)1 μs歸一化感應(yīng)電動(dòng)勢(shì)水平剖面圖;(b)y=0測(cè)線1 μs歸一化感應(yīng)電動(dòng)勢(shì)曲線圖

圖7　多時(shí)間道歸一化感應(yīng)電動(dòng)勢(shì)綜合剖面圖Fig.7　Multi-time gates normalized induced electromotive force comprehensive profile

顯然，從1×10-3ms～7.9×10-2ms的綜合剖面圖呈“U”字形分布形態(tài)，表明“框線效應(yīng)”對(duì)數(shù)據(jù)的畸變影響比較大，這一延遲時(shí)間范圍內(nèi)的數(shù)據(jù)不能用于反演和資料解釋；從2.5×10-1ms后，“框線效應(yīng)”影響消失，曲線形態(tài)正常。因此，選擇的最小延遲時(shí)間，即第一個(gè)時(shí)間道為2.5×10-1ms為宜。同理可知，當(dāng)同時(shí)存在框內(nèi)、外測(cè)點(diǎn)數(shù)據(jù)時(shí)， 選取的最小延遲時(shí)間要比較框內(nèi)測(cè)量時(shí)選取的最小延遲時(shí)間大。這說(shuō)明當(dāng)大定源回線瞬變電磁法采用尺寸較大的發(fā)射框源時(shí)，一部分的早期瞬變電磁數(shù)據(jù)由于受到“框線效應(yīng)”的畸變影響，不得不采取這種剔除數(shù)據(jù)的辦法，這對(duì)于地表的淺部信息反映是一個(gè)較大損失，這也是作者方法在實(shí)際應(yīng)用方面的參考和局限之處。

5　結(jié)論

給出了一種新的計(jì)算均勻半空間大定源回線瞬變電磁場(chǎng)的數(shù)值濾波算法，分別采用了“n點(diǎn)式”Gauss-Legendre數(shù)值求積法，以及47點(diǎn)J1型Hankel變換線性濾波器數(shù)值求積法，單線源組合成大定源回線的線性迭加技術(shù)。與常規(guī)的算法相比較，大大簡(jiǎn)化了計(jì)算步驟，具有較高計(jì)算速度和數(shù)值模擬精度，使編程更容易實(shí)現(xiàn)。通過(guò)對(duì)大定源回線瞬變電磁場(chǎng)的數(shù)值模擬，可以確定出生產(chǎn)實(shí)際中最小延遲時(shí)間，從而有效避免瞬變電磁法的“框線效應(yīng)”現(xiàn)象。

[1]米薩克 N．納比吉安.勘查地球物理電磁法(第一卷)理論[M]．北京：地質(zhì)出版社，1992．

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

[2]蔣邦遠(yuǎn)．實(shí)用近區(qū)磁源瞬變電磁勘探[M]．北京：地質(zhì)出版社，1998．

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

[3]陳貴生．瞬變電磁法在金屬礦產(chǎn)勘查上的應(yīng)用效果及存在問(wèn)題探討[J]．礦產(chǎn)與地質(zhì)，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]薛國(guó)強(qiáng)，李貅．瞬變電磁隧道超前預(yù)報(bào)成像技術(shù)[J]．地球物理學(xué)報(bào)，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號(hào)巖體勘查中的應(yīng)用與研究[J]．吉林地質(zhì)，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]韓自強(qiáng)，羅姣，劉濤，等．定源回線瞬變電磁法在煤礦富水區(qū)調(diào)查中的應(yīng)用[J]．地球物理學(xué)進(jìn)展，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]包乃利，劉鴻福，余傳濤．大定源瞬變電磁法激勵(lì)場(chǎng)及邊框效應(yīng)研究[J]．煤田地質(zhì)與勘探，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]丁艷飛,白登海,許誠(chéng)．均勻半空間表面大定源瞬變電磁響應(yīng)的快速算法[J]．地球物理學(xué)報(bào),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]．北京:科學(xué)出版社,2002．

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

[11]阮百堯．均勻水平大地上頻率域垂直磁偶源電磁場(chǎng)數(shù)值濾波解法[J]．桂林工學(xué)院學(xué)報(bào),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]．北京:冶金工業(yè)出版社,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]張偉，王緒本，覃慶炎．漢克爾變換的數(shù)值計(jì)算與精度的對(duì)比[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]徐世良．計(jì)算機(jī)常用算法(第二版)[M]．北京: 清華大學(xué)出版社，1995．

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

[18]劉云，王緒本，段長(zhǎng)生．大回線瞬變電磁正演模擬在工程實(shí)際中的應(yīng)用[C]．第十屆中國(guó)國(guó)際地球電磁學(xué)術(shù)討論會(huì)論文集，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

國(guó)家“973計(jì)劃”項(xiàng)目(2014CB440905)；國(guó)家自然科學(xué)基金(41440031)；貴州省科學(xué)技術(shù)基金(2014GZ93278)

劉云(1973-)，男，博士，副研究員，研究方向?yàn)榈厍蛭锢頂?shù)值模擬與反演成像，E-mail： liu_yun@mail.gyig.ac.cn。

宋滔(1988-)，男，博士，研究方向?yàn)榈厍蛭锢頂?shù)值模擬及與演成像，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