楊 鑫，姜彥南，魏 兵

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

?

·數(shù)理科學·

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

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

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

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

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

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

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

1 瞬變電磁法基本理論

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

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

(1)

(2)

(3)

(4)

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

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

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

(5)

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

1.2 邊界條件

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

(6)

(7)

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

2 時域有限差分方法

2.1 二階差分法

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

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

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

2.2 高階差分法

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

以1維可導函數(shù)f(x)為例,對f(x)進行關(guān)于x網(wǎng)格劃分,每個網(wǎng)格編號為Δxi(i=1，2，3，…)。則f(x)在每個網(wǎng)格Δxi中點處Δxi-Δxi/2的一階導數(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為

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

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

3 源的驗證和分析

3.1 頻域激勵源

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

(8)

Er(ω，ρ)=

(9)

(10)

(11)

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

3.2 源的數(shù)值計算

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

(12)

其中，x> 0，Jr為r階貝塞爾函數(shù)。對于這種類型的積分計算，可采用高斯數(shù)值積分方法實現(xiàn)。而對于頻域源中復雜的拉普拉斯變換，文中選用Gaver-Stehfest[15]變換法進行數(shù)值變換。由于Gaver-Stehfest變換方法是純實數(shù)域的運算，且只需要少量的離散拉普拉斯變量值S(一般取12～16個)就可實現(xiàn)反變換的計算，所以該方法廣泛應用于求解拉普拉斯逆變換問題。其變換思路大致為：設所求逆變換的時刻為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域的離散點數(shù)。由式(13)可見,Gaver-Stehfest變換可極大地降低拉普拉斯變換計算的復雜性。這樣,通過Hankel積分變換求出頻率域場表達式,然后再對頻率域解作逆拉普拉斯變換,即可求得時空域的初始場分布,計算結(jié)果如圖2所示。圖中顯示了當T=6 ms時電磁場部分分量在地下的縱向分布情況。

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

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

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

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

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

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

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

4 結(jié) 語

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

[1] 丁艷飛, 白登海, 許誠. 均勻半空間表面大定源瞬變電磁響應的快速算法[J].地球物理學報,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)斷時間的回線源激發(fā)TEM三維時域有限差分正演[J].地球物理學報, 2013, 56(3): 1049-1064.

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

[8] 殷長春, 任秀艷, 劉云鶴. 航空瞬變電磁法對地下典型目標體的探測能力研究[J]. 地球物理學報, 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ū)多金屬找礦中的應用研究[J].地球物理學進展, 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] 宋汐瑾, 黨瑞榮, 郭寶龍,等. 井中磁源瞬變電磁響應特征研究[J]. 地球物理學報, 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

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

楊鑫，男，重慶人，博士生，從事計算電磁學研究。

P631

A

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