基于瞬變磁偶極子的正演算法

楊 鑫，姜彥南，魏 兵

(1.西安電子科技大學(xué) 物理與光電工程學(xué)院, 陜西 西安 710071;2.西安電子科技大學(xué) 信息感知技術(shù)協(xié)同創(chuàng)新中心，陜西 西安 710071；3.桂林電子科技大學(xué) 信息與通信工程學(xué)院, 廣西 桂林 541000)

?

·數(shù)理科學(xué)·

基于瞬變磁偶極子的正演算法

楊 鑫1,2，姜彥南3，魏 兵1,2

(1.西安電子科技大學(xué) 物理與光電工程學(xué)院, 陜西 西安 710071;2.西安電子科技大學(xué) 信息感知技術(shù)協(xié)同創(chuàng)新中心，陜西 西安 710071；3.桂林電子科技大學(xué) 信息與通信工程學(xué)院, 廣西 桂林 541000)

以時(shí)域有限差分方法為技術(shù)手段,研究了瞬態(tài)電磁法中磁偶極子瞬態(tài)場的初始條件實(shí)現(xiàn)方法,提出了一種帶邊界條件的瞬變磁偶極源正演算法。文中先利用數(shù)值積分和Gaver-Stehfest變換等方法,實(shí)現(xiàn)了頻域瞬態(tài)偶極源至?xí)r域空間分布的計(jì)算，采用高階差分方法克服二階差分方程帶來的計(jì)算不穩(wěn)定性。然后結(jié)合地-空邊界和地下邊界條件,獲得瞬態(tài)磁偶極子源的初始場分布,模擬了電波在地下的擴(kuò)散過程,并給出了地表的電動(dòng)勢響應(yīng)曲線。算例的對比結(jié)果驗(yàn)證了文中算法的有效性。

瞬態(tài)場;磁偶極子;時(shí)域有限差分方法;正演計(jì)算

瞬變電磁法(Transient electromagnetic methods)是以電磁感應(yīng)原理為基礎(chǔ)的時(shí)間域人工場源電磁勘探方法。通過對地下響應(yīng)信號的提取和處理,分析異常體的結(jié)構(gòu)形態(tài)、導(dǎo)電性和埋藏深度等特性,達(dá)到找礦和解決地質(zhì)問題等目的。相對于眾多物探方法,瞬變電磁法具有探測深度大、分辨率高和適應(yīng)性強(qiáng)等優(yōu)點(diǎn),廣泛應(yīng)用于工程以及環(huán)境地質(zhì)勘探中[1]。

20世紀(jì)80年代之前，瞬變電磁法的研究主要限于1維和2維模型。之后，2.5維和3維模型成為研究熱點(diǎn)[2-4]。在地電模型的計(jì)算中，根據(jù)地形或環(huán)境的不同常需采用不同形式模擬源，如矩形源、環(huán)形源或者陣列源等[5-6]。而針對不同類型源所激發(fā)的瞬態(tài)場，往往需采用不同的計(jì)算方式以獲得相應(yīng)的數(shù)學(xué)模型。并且，同一物理類型的源因激勵(lì)方式的不同也需進(jìn)行不同的數(shù)學(xué)推導(dǎo)，比如陳曙東等在2012年通過尋找沖激響應(yīng)函數(shù)等方式對自由空間中的回線源在方波和梯形波等不同條件下分別推導(dǎo)出源的瞬變電磁響應(yīng)[7]；2015年殷長春等人利用積分方程等方法推導(dǎo)出了均勻半空間中環(huán)形源在3維空間中的響應(yīng)[8]；2015年Mark Goldman以及Kinza Haider等人研究了瞬變電磁法在海洋探測中的應(yīng)用，并指出利用偶極方法在海洋探測中還處于理論探索階段[9-10]。為了提高探測效率而發(fā)展出的混合探測法，比如瞬變電磁法與直流法(Direct current method)相結(jié)合[11]等。在這些研究中，瞬變電磁法對于一些簡單環(huán)境的計(jì)算是有效的，但對于更復(fù)雜的地電模型如異常非規(guī)則體分布和山區(qū)地形等，這樣的算法就受到了限制[12]。為有效地計(jì)算這些復(fù)雜的模型，就需要解決諸多問題，比如復(fù)雜瞬變源的計(jì)算、邊界條件的確定等，但國內(nèi)外目前關(guān)于瞬變激勵(lì)源與邊界條件的深入研究還不多。理論上講，任意復(fù)雜源都可視為一定數(shù)量偶極子組成，因此，文中對磁偶極子源進(jìn)行深入研究，提出了一種以頻率域磁偶極子作為瞬變激勵(lì)源并結(jié)合吸收邊界條件的新方法。從頻域偶極源出發(fā)，通過一定的數(shù)值計(jì)算手段實(shí)現(xiàn)瞬態(tài)場初始條件計(jì)算，極大地簡化瞬態(tài)源計(jì)算以及程序?qū)崿F(xiàn)的復(fù)雜性，為復(fù)雜源的加入及復(fù)雜環(huán)境下的瞬態(tài)場計(jì)算提供了條件。文中采用時(shí)域有限差分方法(FDTD)實(shí)現(xiàn)地下擴(kuò)散場的計(jì)算。

1 瞬變電磁法基本理論

1.1 瞬態(tài)場擴(kuò)散方程

在用瞬變電磁法處理地球物理問題時(shí),大地可近似為無源、線性和各向同性的媒質(zhì)。瞬變電磁波在地下擴(kuò)散過程中,由于位移電流?D(r,t)/?t十分微弱,特別是在擴(kuò)散后期,位移電流與傳導(dǎo)電流J(r,t)相比可忽略不計(jì)。因此,地下擴(kuò)散電波方程組可忽略位移電流項(xiàng)，其方程組為

(1)

(2)

(3)

(4)

在用FDTD對電磁場進(jìn)行數(shù)值計(jì)算時(shí),一般僅需式(1)和式(2)就可表征電磁波的運(yùn)動(dòng)過程。場在擴(kuò)散晚期表現(xiàn)出似穩(wěn)性,此時(shí)式(1)近似為

這樣，式(3)就不能由式(1)導(dǎo)出,就導(dǎo)致獨(dú)立方程不僅限于式(1)和式(2)。根據(jù)霍姆赫茲定理可知,要確定一個(gè)矢量,必須同時(shí)給定它的散度和旋度。若此時(shí)再用式(1)和式(2)對電磁場進(jìn)行數(shù)值計(jì)算,就會因?yàn)锽(r,t)的不確定性而導(dǎo)致計(jì)算錯(cuò)誤。所以,在擴(kuò)散場的后期,計(jì)算時(shí)就必須將式(3)顯式地包含在內(nèi)，才能保證計(jì)算的正確性。這樣獨(dú)立方程也就變?yōu)槭?1)～(3),構(gòu)成描述地下瞬態(tài)電磁波運(yùn)動(dòng)的方程組。

FDTD處理電磁計(jì)算問題時(shí)，其原理是將場方程轉(zhuǎn)化為關(guān)于計(jì)算時(shí)間的顯示差分迭代方程。因此，方程中必須存在場量對時(shí)間的偏導(dǎo)數(shù)。式(1)能滿足了這一條件，但式(2)卻不然。為解決這一問題，引入了Du Fort-Frankel方法。該方法在式(2)中引入虛擬位移電流[13]Jd(r,t)=γ?E(r,t)/?t后，變?yōu)?/p>

(5)

這樣,結(jié)合式(1)，式(3)和式(5),就可實(shí)現(xiàn)地下電波的模擬計(jì)算。

1.2 邊界條件

地層一般包含多種媒質(zhì)(比如礦藏、水和氣體等),就必定存在諸多不連續(xù)邊界。而利用計(jì)算機(jī)處理的地電模型必定為有限區(qū)域,這種有限的區(qū)域具有的邊界稱為截?cái)噙吔?。所?計(jì)算地電模型時(shí)將存在兩種邊界:不同媒質(zhì)分界面和截?cái)噙吔纭Ｓ晌ㄒ恍远ɡ砜芍?要確定某區(qū)域的電磁場分布,就必須給定邊界處的場分布。對于一些簡單的地球物理模型,常用的地下邊界處理辦法是強(qiáng)制加狄利克雷邊界條件,比如給定計(jì)算域邊界處的磁場或者電場。但對于復(fù)雜的地下環(huán)境,確定地下邊界處的場分布就顯得比較困難,所以文中采用吸收邊界條件方法處理地下截?cái)噙吔?。但對于地空邊界條件的設(shè)置,由于瞬態(tài)源位于地面,若用吸收邊界條件對其進(jìn)行處理,則在地電模型中就相當(dāng)于將源置位于吸收邊界。這種設(shè)置方法將導(dǎo)致地表處場的計(jì)算錯(cuò)誤,所以對于地空邊界還需另行設(shè)置。在直角坐標(biāo)系下,式(6)和式(7) 給出了地空邊界的一種設(shè)置方法。取z軸正向指向自由空間,z=0為地表面,則磁場在地表處沿x和y方向的地空邊界條件分別為

(6)

(7)

式(6)和式(7)給出的物理意義為:若知道波數(shù)域中地表處磁場垂直分量值,就可外推計(jì)算地表上高h(yuǎn)處的磁場水平分量,即由內(nèi)部場確定邊界場，從而實(shí)現(xiàn)地空邊界條件的計(jì)算。

2 時(shí)域有限差分方法

2.1 二階差分法

以三維空間麥克斯韋方程組為例給出時(shí)域有限差分方法的基本計(jì)算思路。首先將電磁波傳播的空間域劃分為一定數(shù)量的方格。通過空間的這種離散劃分,將電場和磁場在空間網(wǎng)格上進(jìn)行交錯(cuò)排列便構(gòu)成了時(shí)域有限差分方法的計(jì)算基礎(chǔ)。電磁場各個(gè)分量的空間離散位置如圖1所示。

圖1 三維空間離散網(wǎng)格及其場分布Fig.1 3D FDTD grid and fields distribution

再將連續(xù)時(shí)間流t離散為一定數(shù)量的時(shí)間間隔(時(shí)間步),Δtn=tn+1-tn,n取正整數(shù)。tn代表時(shí)間的離散點(diǎn),并約定在該整數(shù)離散時(shí)刻點(diǎn)計(jì)算電場值。用tn+Δtn/2代表每個(gè)時(shí)間步的中間時(shí)刻點(diǎn),約定用該時(shí)刻點(diǎn)計(jì)算磁場值。按這種關(guān)系對式(1)、式(2)和式(3)進(jìn)行離散,即可得電磁場的二階差分迭代關(guān)系式。以式(1)的Bx分量為例:

2.2 高階差分法

在電流源關(guān)斷的早期,擴(kuò)散場包含的頻率十分豐富,但集中表現(xiàn)為高頻段。所以對早期擴(kuò)散場進(jìn)行數(shù)值計(jì)算時(shí),應(yīng)該視情況而使用合適的高階差分。文中選用四階差分,由后面的計(jì)算結(jié)果可以表明,這種差分可以近似地計(jì)算早期電磁場的擴(kuò)散過程。下面給出四階差分的實(shí)現(xiàn)方法。

以1維可導(dǎo)函數(shù)f(x)為例,對f(x)進(jìn)行關(guān)于x網(wǎng)格劃分,每個(gè)網(wǎng)格編號為Δxi(i=1，2，3，…)。則f(x)在每個(gè)網(wǎng)格Δxi中點(diǎn)處Δxi-Δxi/2的一階導(dǎo)數(shù)可由以下關(guān)系式求出，

f ′(x)|i+1/2=

a-1f(x)|i-1+a0f(x)|i+

a1f(x)|i+1+a2f(x)|i+2。

其中,參數(shù)ai(i=-1，0，1，2)可由如下方程組求出，

其中，yi為

將地下擴(kuò)散方程按上式離散，可得其四階差分離散形式。以Bx為例給出離散后的四階差分方程，

類似地，可得磁場其他分量以及電場分量的四階差分關(guān)系式。得到差分關(guān)系后,只要離散參數(shù)間滿足空間以及時(shí)間的計(jì)算穩(wěn)定性要求，就可以實(shí)現(xiàn)擴(kuò)散方程的數(shù)值計(jì)算。

3 源的驗(yàn)證和分析

3.1 頻域激勵(lì)源

在圓柱坐標(biāo)系下,位于均勻大地表面的磁偶極子源關(guān)斷瞬間,電磁場在地下的頻域解[14]為

(8)

Er(ω，ρ)=

(9)

(10)

(11)

由式(11)理論上能夠計(jì)算出時(shí)空域的場分布，從而求出場的初始條件。但由于存在emz大積分因子,在計(jì)算機(jī)實(shí)際計(jì)算時(shí),計(jì)算的空間大小一般可達(dá)數(shù)百或數(shù)千米,可見在這種計(jì)算尺度下,該積分因子的計(jì)算結(jié)果將容易溢出。而對一般個(gè)人計(jì)算機(jī)硬件資源而言,計(jì)算區(qū)域的寬度和深度是十分有限的。因此,以時(shí)域解析解為初始條件存在一定的局限性。為克服這一難點(diǎn),文中采用數(shù)值積分法和Gaver-Stehfest方法求解時(shí)域初始場分布。

3.2 源的數(shù)值計(jì)算

源的頻域方程式(8)～(10)為典型的Hankel積分變換形式,因此,可將其簡寫為

(12)

其中，x> 0，Jr為r階貝塞爾函數(shù)。對于這種類型的積分計(jì)算，可采用高斯數(shù)值積分方法實(shí)現(xiàn)。而對于頻域源中復(fù)雜的拉普拉斯變換，文中選用Gaver-Stehfest[15]變換法進(jìn)行數(shù)值變換。由于Gaver-Stehfest變換方法是純實(shí)數(shù)域的運(yùn)算，且只需要少量的離散拉普拉斯變量值S(一般取12～16個(gè))就可實(shí)現(xiàn)反變換的計(jì)算，所以該方法廣泛應(yīng)用于求解拉普拉斯逆變換問題。其變換思路大致為：設(shè)所求逆變換的時(shí)刻為t，先將S域中的拉普拉斯變量離散為Sm= m·ln(2)/t，再把Sm帶回拉普拉斯式中，可得離散的拉普拉斯函數(shù)F(Sm)。最后，用式(13)表示的Gaver-Stehfest變換式求函數(shù)F(s)拉普拉斯逆變換值f(t)，

(13)

其中,Km為Gaver-Stehfest變換系數(shù),n為變換所需S域的離散點(diǎn)數(shù)。由式(13)可見,Gaver-Stehfest變換可極大地降低拉普拉斯變換計(jì)算的復(fù)雜性。這樣,通過Hankel積分變換求出頻率域場表達(dá)式,然后再對頻率域解作逆拉普拉斯變換,即可求得時(shí)空域的初始場分布,計(jì)算結(jié)果如圖2所示。圖中顯示了當(dāng)T=6 ms時(shí)電磁場部分分量在地下的縱向分布情況。

由圖2可見,在源關(guān)斷初期,場的能量集中分布在源附近,計(jì)算結(jié)果與瞬態(tài)源關(guān)斷瞬間的實(shí)際場分布相符合。電場隨著與源點(diǎn)的距離增加而呈明顯的高次衰減分布,且主要分布于地表淺層。此外,電場中所攜帶的能量明顯低于磁場中所攜帶的能量,并且當(dāng)遠(yuǎn)離源時(shí)能量主要存在于磁場中,符合電磁能量在導(dǎo)體中的分布特點(diǎn)。圖2中還清晰地顯示磁場分量分布的“煙圈”效應(yīng),這些“煙圈”中心的位置隨著剖面與源的距離加大而逐漸加深,體現(xiàn)了瞬變場的擴(kuò)散特點(diǎn)。

圖2 地下場縱向剖面圖Fig.2 Longitudinal crossection of underground fields

3.3 激勵(lì)源數(shù)值實(shí)現(xiàn)及驗(yàn)證

圖3～5給出以頻率域磁偶極子源作為初始條件的地下瞬變電磁場的計(jì)算結(jié)果。圖3給出兩個(gè)不同時(shí)刻的地表磁場的垂直分量響應(yīng)電動(dòng)勢數(shù)值解及時(shí)域解析解對比,可見數(shù)值計(jì)算結(jié)果與解析解能很好地吻合。由于時(shí)域解析源是關(guān)于1/ρ(ρ為觀察點(diǎn)至源點(diǎn)的距離)的函數(shù)[15],可見其在源點(diǎn)處存在奇異性,所以場在源點(diǎn)會出現(xiàn)不穩(wěn)定性,這種不穩(wěn)定性隨著場的迭代計(jì)算而逐漸向周圍擴(kuò)散。而數(shù)值解是通過網(wǎng)格離散來實(shí)現(xiàn),就可以避開源點(diǎn)的奇異性,實(shí)現(xiàn)穩(wěn)定的計(jì)算。圖4給出了在地表x軸上50m和200m處的磁場垂直分量的電動(dòng)勢響應(yīng)隨時(shí)間的變化情況。圖5給出了不同時(shí)刻磁場垂直分量在剖面y=0處的地下分布圖,很好地顯示了地下場的擴(kuò)散情況。在擴(kuò)散初期,磁場垂直分量的能量主要集中在地表源附近。隨著時(shí)間的增加,場逐步往地下和水平方向擴(kuò)散,擴(kuò)散過程顯示出了明顯的“煙圈”效應(yīng)。在場計(jì)算的后期,磁場垂直分量的“煙圈”效應(yīng)逐漸消失于整個(gè)計(jì)算域,瞬態(tài)場完成在地下的擴(kuò)散。

圖3 數(shù)值解與解析解對比曲線Fig.3 Comparison of numerical solutions with analytical solutions

圖4 地表磁場垂直分量響應(yīng)曲線Fig.4 Ground surface response curve of the vertical magnetic field

圖5 不同時(shí)刻地下磁場垂直分量分布Fig.5 Vertical component distribution of underground magnetic fields at different time

4 結(jié) 語

從磁偶極子的瞬變場頻域特性出發(fā),計(jì)算電波在地下的時(shí)域場分布,提出了一種結(jié)合吸收邊界條件的瞬變場初始條件計(jì)算方法。先分別利用Gaver-Stehfest變換方法和高斯數(shù)值積分方法求解頻域磁偶極子源在地下的時(shí)域場分布,然后由地表場外推自由空間中的場以確定地空邊界條件,最后利用FDTD實(shí)現(xiàn)了瞬變偶極場的電動(dòng)勢響應(yīng)，計(jì)算并給出了不同時(shí)刻地下場的分布圖,清晰地顯示了地下擴(kuò)散場的“煙圈”效應(yīng)。由驗(yàn)證結(jié)果可知，文中正確地模擬了瞬變磁偶極子源在地下的擴(kuò)散過程。該算法可將源的復(fù)雜計(jì)算過程簡化為基本的數(shù)值計(jì)算,避免了大積分因子的計(jì)算,極大地降低了計(jì)算的復(fù)雜度。文中算法還具有通用性較強(qiáng)的特點(diǎn),可為復(fù)雜瞬變源的模擬提供可靠的參考。

[1] 丁艷飛, 白登海, 許誠. 均勻半空間表面大定源瞬變電磁響應(yīng)的快速算法[J].地球物理學(xué)報(bào),2012,55(6):2087-2096.

[2] SANFILIPO W A, EATON P A, HOHMANN G W. The transient EM response of a prism in a conductive half-space[J].Geophysics, 1985, 50(2): 272.

[3] ENDO M, NOGUCHI K.Three-dimensional modeling considering the topography for the case of the time-domain electromagnetic method[C].2nd International Symposium on 3-Dimensional Electromagnetics, 1999.

[4] MAAO F A. Fast finite-difference time-domain modeling for marine-subsurface electromagnetic problems[J].Geophysics, 2007, 72: 19-23.

[5] UME S, HARRIS J M, ALUMBAUGH D L. A finite element algorithm for 3-D transient electromagnetic modeling[C]∥2009 SEG Annual Meeting. Society of Exploration Geophysicists, 2009.

[6] 孫懷鳳, 李貅, 李術(shù)才，等. 考慮關(guān)斷時(shí)間的回線源激發(fā)TEM三維時(shí)域有限差分正演[J].地球物理學(xué)報(bào), 2013, 56(3): 1049-1064.

[7] 陳曙東, 林君, 張爽. 發(fā)射電流波形對瞬變電磁響應(yīng)的影響[J]. 地球物理學(xué)報(bào), 2012, 55(2): 709-716.

[8] 殷長春, 任秀艷, 劉云鶴. 航空瞬變電磁法對地下典型目標(biāo)體的探測能力研究[J]. 地球物理學(xué)報(bào), 2015, 58(9): 3370-3379.

[9] GOLDMAN M, MOGILATOV V, HAROON A,et al. Signal detectability of marine electromagnetic methods in the exploration of resistive targets[J]. Geophysical Prospecting, 2015, 63(1):192-210.

[10] HAIDER K, ENGESGAARD P, SONNENBORG T O,et al. Numerical modeling of salinity distribution and submarine groundwater discharge to a coastal lagoon in Denmark based on airborne electromagnetic data[J]. Hydrogeology Journal, 2015, 23(2): 217-233.

[11] CHENG Jiulong, LI Fei, PENG Suping, et al. Joint inversion of TEM and DC in roadway advanced detection based on particle swarm optimization[J]. Journal of Applied Geophysics, 2015, 123: 30-35.

[12] 魏明君, 丁云河, 李冰,等. 線源瞬變電磁法在熊耳山區(qū)多金屬找礦中的應(yīng)用研究[J].地球物理學(xué)進(jìn)展, 2015, 30(1): 146-152.

[13] WANG T, HOHMANN G W. A finite-difference time-domain solution for three-dimensional electromagnetic modeling[J]. Geophysics, 1993, 58(6): 797-809.

[14] 米薩克·納比吉安. 勘查地球物理電磁法第一卷理論[M]. 北京: 地質(zhì)出版社, 1992.

[15] 宋汐瑾, 黨瑞榮, 郭寶龍,等. 井中磁源瞬變電磁響應(yīng)特征研究[J]. 地球物理學(xué)報(bào), 2011, 54(4):1122-1129.

(編 輯 李 靜)

The transient properties of magnetic dipole source

YANG Xin1,2, JIANG Yannan3, WEI Bing1,2

(1.School of Physics and Optoelectronic Engineering,Xidian University, Xi′an 710071, China;2.Collaborative Innovation Center of Information Sensing and Understanding, Xidian University, Xi′an 710071， China；3.School of Information and Communication Engineering, Guilin University of Electronic Technology, Guilin 541000, China)

The calculation method of initial condition for transient magnetic dipole is studied by means of FDTD, and a forward algorithm for its transient fields is proposed with a boundary condition. The numerical integration and Gaver-Stehfest transform are applied to achieve the magnetic dipole fields from frequency domain to time domain. And the high-order difference method is used to overcome the calculation instability of second order differential equations. Combinating with the ground-air boundary and subsurface boundary conditions, the initial distributed field of transient magnetic dipole is obtained. The diffusion of the EM wave of the underground is simulated and the EMF response curve at the earth′s surface is presented as well. Finally, the numerical results show the validation of the algorithm.

transient field; magnetic dipole; FDTD method; forward calculation

2015-09-30

國家自然科學(xué)基金資助項(xiàng)目(61231003，61401344，61571348)；廣西自然科學(xué)基金資助項(xiàng)目(2014GXNSFAA118283)

楊鑫，男，重慶人，博士生，從事計(jì)算電磁學(xué)研究。

P631

A

10.16152/j.cnki.xdxbzr.2016-06-005