航空時(shí)間域電磁法回線源有限差分初始場(chǎng)計(jì)算

許洋鋮 林 君 嵇艷鞠 鞏源峰 關(guān)珊珊

(吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院/地球信息探測(cè)儀器教育部重點(diǎn)實(shí)驗(yàn)室,吉林 長(zhǎng)春130025)

再根據(jù)垂直磁場(chǎng)和水平磁場(chǎng)的關(guān)系[15],求出磁場(chǎng)水平分量。圖1給出了采用回線積分法和納比積安中心回線解析法計(jì)算的圓形線圈均勻半空間模型,回線中心點(diǎn)垂直二次磁場(chǎng)響應(yīng)曲線。為了對(duì)比各種方法,文中各種方法均采用均勻半空間模型,同樣的計(jì)算參數(shù),如表1所示。

由圖1中可見,回線積分法和解析表達(dá)式計(jì)算結(jié)果一致,誤差精度能達(dá)到0.000 1%,計(jì)算精度完全滿足有限差分初始場(chǎng)的要求。

1.引 言

直升機(jī)時(shí)間域電磁系統(tǒng)(airborne time domain electromagnetic)具有效率高、成本低、分辨率高、操控靈活等優(yōu)點(diǎn),被廣泛應(yīng)用于礦產(chǎn)資源、地質(zhì)調(diào)查和環(huán)境監(jiān)測(cè)等領(lǐng)域。世界上發(fā)達(dá)國(guó)家已經(jīng)將直升機(jī)航空電磁測(cè)量系統(tǒng)用于商業(yè)化服務(wù)[1],我國(guó)立項(xiàng)開展吊艙式時(shí)間域直升機(jī)航空電磁勘查系統(tǒng)的研發(fā)。

目前,時(shí)間域航空電磁野外數(shù)據(jù),多采用層狀均勻大地模型進(jìn)行反演解釋,當(dāng)在電阻率橫向變化比較明顯的區(qū)域,層狀均勻大地不能反映真實(shí)的地質(zhì)情況,只有二維、三維異常體模型更符合客觀實(shí)際地質(zhì)條件。時(shí)間域有限差分方法(a finite-difference approach for time-domain solution,簡(jiǎn)稱 FDTD),比較適合模擬地下復(fù)雜形狀的三維導(dǎo)電體以及電阻率劇烈變化的模型,具有計(jì)算簡(jiǎn)單、快速等特點(diǎn)。

對(duì)于時(shí)間域瞬變電磁響應(yīng)的數(shù)值計(jì)算,1985年,HOHMANN和ADHIDJAJA[2]在YEE[3]提出交錯(cuò)網(wǎng)格的基礎(chǔ)上,實(shí)現(xiàn)了二維導(dǎo)電體的顯式有限差分瞬變電磁響應(yīng)正演,1992年,LEPPIN M[4]利用隱格式有限差分,計(jì)算了方形線圈三維導(dǎo)電體的瞬變電磁響應(yīng),1993年,WANG T和HOHMANN G W[5]采用FDTD方法實(shí)現(xiàn)了回線源的三維瞬變電磁正演問(wèn)題,近年來(lái),魏文博、譚捍東[6]和宋維琪[7]等研究了回線源三維瞬變電磁有限差分正演計(jì)算,閆述[8]、李桐林[9]等也研究了二維長(zhǎng)導(dǎo)線源有限差分正演,都取得了一定的成果。

采用有限差分法計(jì)算時(shí)間域電磁響應(yīng)時(shí),初始場(chǎng)值的計(jì)算和邊界條件的選擇直接決定了計(jì)算精度和速度。HOHMAN和ADHIDJAJA等人的計(jì)算都是將回線源初始場(chǎng)等效成有限個(gè)磁偶極子初始場(chǎng)的和。在吊艙式時(shí)間域航空電磁系統(tǒng),采用圓形回線發(fā)射、中心裝置,初始場(chǎng)計(jì)算含有雙重貝塞爾函數(shù),若采用等效成磁偶極子進(jìn)行近似,會(huì)給迭代帶來(lái)較大誤差[10],而且對(duì)每個(gè)剖分網(wǎng)格都要計(jì)算,將導(dǎo)致計(jì)算量變得非常龐大,為此,本文研究了圓形回線源的初始場(chǎng)數(shù)值計(jì)算方法,提出了混合法計(jì)算初始電磁場(chǎng)。

2.初始場(chǎng)數(shù)值計(jì)算

有限差分法中,對(duì)于初始時(shí)刻的瞬變電磁場(chǎng)的計(jì)算,根據(jù)電磁場(chǎng)理論[11],在足夠短的時(shí)間內(nèi),假設(shè)大地的上部分是均勻半空間,采用均勻半空間電磁響應(yīng)的積分方程進(jìn)行計(jì)算初始電磁場(chǎng)。選擇初始時(shí)間[12]為

式中：μ1是大地表層磁導(dǎo)率;σ1是大地表層電導(dǎo)率;Δmin是最小網(wǎng)格步長(zhǎng)。在地面及空氣中圓回線的二次場(chǎng)電磁響應(yīng)表達(dá)式為

式中 ：核函數(shù) F0=rTEeu0(z-h)λ,u0=(λ2-ω2μ0ε0+iω μ0σ0)1/2,為空氣中磁導(dǎo)率,I為發(fā)射電流,a 為發(fā)射回線半徑,h為發(fā)射線圈距地面高度,z為接收線圈距地面距離,rTE為地表反射系數(shù),J1為一階貝塞爾函數(shù),為零階貝塞爾函數(shù)。

上述式(2)～(4)中都含有雙重貝塞爾函數(shù),無(wú)法用解析式表達(dá)出來(lái)。由于有限差分方程迭代的計(jì)算量大,為此,要求初始場(chǎng)值的計(jì)算方法不僅精度高,而且速度快。針對(duì)式(2)～(4)采用了回線積分法、宗量近似法、貝塞爾函數(shù)展開法計(jì)算初始電磁場(chǎng)值,利用漢克爾變換[13]、高斯積分和Guptasarma數(shù)字濾波方法[14],實(shí)現(xiàn)了時(shí)間域電磁響應(yīng)計(jì)算。

2.1 回線積分法

回線積分法是將電偶極子產(chǎn)生的切向電場(chǎng)分量(不含接地項(xiàng))沿發(fā)射回線進(jìn)行積分,先得到發(fā)射回線的感應(yīng)電壓,再根據(jù)互易等效原理,計(jì)算接收線圈的感應(yīng)電壓,最后根據(jù)法拉第電磁感應(yīng)定律,以及圓形線圈與方形回線面積等效,得到頻率域二次磁場(chǎng)表達(dá)式,為

再根據(jù)垂直磁場(chǎng)和水平磁場(chǎng)的關(guān)系[15],求出磁場(chǎng)水平分量。圖1給出了采用回線積分法和納比積安中心回線解析法計(jì)算的圓形線圈均勻半空間模型,回線中心點(diǎn)垂直二次磁場(chǎng)響應(yīng)曲線。為了對(duì)比各種方法,文中各種方法均采用均勻半空間模型,同樣的計(jì)算參數(shù),如表1所示。

由圖1中可見,回線積分法和解析表達(dá)式計(jì)算結(jié)果一致,誤差精度能達(dá)到0.000 1%,計(jì)算精度完全滿足有限差分初始場(chǎng)的要求。

圖2為采用回線積分法計(jì)算地面垂直二次磁場(chǎng)分布。計(jì)算時(shí)間為71 898.547 32 s。雖然沿回線積分的精度很高,可是它的計(jì)算速度太慢,不能滿足FDTD對(duì)速度的要求。

圖1 回線積分法中心點(diǎn)垂直磁場(chǎng)

表1 幾種方法的計(jì)算參數(shù)

圖2 回線積分法地面垂直磁場(chǎng)

由于發(fā)射回線在空間各處電磁場(chǎng)響應(yīng)沒(méi)有解析式(地面中心點(diǎn)除外),而回線積分法在中心點(diǎn)的計(jì)算值,同與其他方法相比,更接近理論值,所以下面在地面其他點(diǎn)將采用回線積分法作為比較標(biāo)準(zhǔn)。

2.2 宗量近似法

宗量近似法是將貝塞爾函數(shù)根據(jù)求解區(qū)域,進(jìn)行大、小宗量近似。在回線外,當(dāng)a?r時(shí),采用小宗量近似,有以垂直磁場(chǎng)為例,磁場(chǎng)響應(yīng)表達(dá)式(2)寫為

電磁響應(yīng)表達(dá)式簡(jiǎn)化為含一重貝塞爾函數(shù)積分,計(jì)算變得簡(jiǎn)單。圖3為采用小宗量近似法計(jì)算的地面各點(diǎn)垂直磁場(chǎng)響應(yīng)。

圖3 小宗量近似計(jì)算的地面垂直磁場(chǎng)

計(jì)算所用時(shí)間為0.323 043 s,從圖3中可看出,小宗量近似只有在a?r和a?r時(shí)才與回線積分法計(jì)算的垂直磁場(chǎng)有很好近似,在a和r相差不大的情況下有較大誤差,可達(dá)5%,此時(shí)不適合有限差分迭代計(jì)算。

采用大宗量近似時(shí)[16],利用Hankel函數(shù)的漸進(jìn)特性,將磁場(chǎng)積分式分為兩部分,即

由Hankel函數(shù)的大宗量近似性質(zhì)和折線逼近算法有

利用貝塞爾函數(shù)的導(dǎo)數(shù)關(guān)系以及差商代替導(dǎo)數(shù),有

式中：

此時(shí)將包含雙重貝塞爾函數(shù)的無(wú)限積分轉(zhuǎn)換為只含有一個(gè)貝塞爾函數(shù)積分的有限個(gè)有限區(qū)間的和,在每個(gè)有限區(qū)間上采用7點(diǎn)高斯積分。圖4為采用大宗量近似計(jì)算的地面電磁響應(yīng)分布。

圖4 大宗量近似法地面垂直磁場(chǎng)

計(jì)算時(shí)間為223.977549 s。由圖4中可以看出,大宗量近似與回線積分的誤差精度能達(dá)到0.01%,但是由于有有限個(gè)區(qū)間的積分,雖然計(jì)算速度比回線積分法快,但是還是不適合計(jì)算大量的初始場(chǎng),所以不適合FDTD。

2.3 貝塞爾函數(shù)展開法

貝塞爾函數(shù)展開法是將含有發(fā)射半徑的貝塞爾函數(shù)J0(λ a)項(xiàng)寫入核函數(shù)中,核函數(shù)寫為：F1=rTE eu0(z-h)λ J1(λ a)。將 J1(λ a)進(jìn)行級(jí)數(shù)展開 ,以垂直磁場(chǎng)為例：

將貝塞爾函數(shù)展式(17)代入表達(dá)式(16)中,直接計(jì)算即可得到電磁場(chǎng)響應(yīng)。圖5為貝塞爾函數(shù)展開法計(jì)算的電磁場(chǎng)響應(yīng)分布。

圖5 貝塞爾函數(shù)展開法地面垂直磁場(chǎng)

計(jì)算時(shí)間為1.941 0311 s,從圖5中可看出,用貝塞爾函數(shù)展開直接計(jì)算的垂直磁場(chǎng)在地面各點(diǎn)的響應(yīng)曲線,在r≤30 a時(shí)同回線積分法計(jì)算誤差精度能達(dá)到0.01%,但是在r＞30 a時(shí)開始出現(xiàn)震蕩,這是由于核函數(shù)中引入了貝塞爾函數(shù)的緣故,所以只能計(jì)算r≤30 a的初始電磁場(chǎng)。

2.4 幾種算法對(duì)比分析

將回線積分法、宗量近似法、貝塞爾函數(shù)展開法進(jìn)行對(duì)比、分析,各種方法的計(jì)算精度、速度、適用區(qū)間,如表2所示。

表2 幾種方法比較

回線積分法精度較高,可以計(jì)算全空間的電磁場(chǎng),但是計(jì)算速度慢,不適合大規(guī)模計(jì)算電磁場(chǎng)。

小宗量近似法計(jì)算速度較快,但是只能是在極端的情況下,精度才能滿足要求;在收發(fā)距與圓線圈半徑相差不大的情況下,誤差比較大,不能計(jì)算全部空間的電磁場(chǎng);大宗量近似雖然精度高,可以計(jì)算全部空間的電磁場(chǎng),但是也受到計(jì)算速度的限制。

貝塞爾函數(shù)展開法計(jì)算較為簡(jiǎn)單,計(jì)算時(shí)間快,但是只能在r≤30 a時(shí)精度滿足要求,在r＞30 a時(shí)出現(xiàn)震蕩。

3.混合法計(jì)算初始電磁場(chǎng)

采用有限差分計(jì)算時(shí),需要計(jì)算大地表層的初始電磁場(chǎng),網(wǎng)格劃分越細(xì),計(jì)算精度越高,計(jì)算量也越大,為此,提出了采用混合法進(jìn)行計(jì)算初始電磁場(chǎng),在r≤30 a采用貝塞爾函數(shù)展開法、在r＞30 a時(shí)采用小宗量近似法計(jì)算FDTD的初始場(chǎng)。

圖6為采用混合法計(jì)算均勻大地表面的時(shí)電場(chǎng)和時(shí)磁場(chǎng)分布。由圖中可見,初始場(chǎng)基本沒(méi)有截?cái)嗾`差。

4.結(jié) 論

在計(jì)算有限差分法初始電磁場(chǎng)時(shí),需要計(jì)算大量的電磁場(chǎng),對(duì)計(jì)算精度和速度都有要求。通過(guò)對(duì)回線積分法、宗量近似法、貝塞爾函數(shù)展開法三種方法計(jì)算初始場(chǎng)進(jìn)行對(duì)比,研究表明,為了計(jì)算全空間的電磁場(chǎng)初值,采用任何單一的方法,都存在缺陷,只有采用混合法,將貝塞爾函數(shù)展開與小宗量近似法結(jié)合,不僅計(jì)算速度快,而且精度高,滿足FDTD對(duì)計(jì)算時(shí)間和計(jì)算精度的要求,適合FDTD的初始場(chǎng)計(jì)算。

[1]雷 棟,胡祥云,張素芳.航空電磁法的發(fā)展現(xiàn)狀[J].地質(zhì)找礦論叢,2006,21(1)：40-53.LEI Dong,HU Xiangyun,ZHANG Sufang.Development status quo of airborne electromagnetic[J].Contributions to Geology and Mineral Resources Research,2006,21(1)：40-53.(in Chinese)

[2]ADHIDJIAJA J I,HOHM ANN G W,ORISTAGLIO M L.Two-dimensional transient electrom-agnetic responses[J].Geophsics,1985,50(12)：2849-2861.

[3]YEE K S.Numerical solution of initial boundary value problems involvingMaxwell equations in isotropic media[J].IEEE Trans.Antennas propa-gat,May 1966,AP-14(3)：302-307.

[4]LEPPIN M.Electromagnetic modeling of 3-D Sources over 2-D inhomogeneities in the time do-Main[J].Geophsics,1992,57(8)：994-1003.

[5]WANG T and HOHMANN G W.A finite-difference,time-domain solution for three-dimensiona lelectromagnetic modeling[J].Geophsics,1993,58(6)：797-809.

[6]譚捍東,余欽范,John Booker,等.大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)報(bào),2003,46(5)：705-711.TAN Handong,YU Qinfan,John Booker,et al.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J].Chinese Journal of Geophysics,2003,46(5)：705-711.(in Chinese)

[7]宋維琪,仝兆歧.3D瞬變電磁場(chǎng)的有限差分正演計(jì)算[J].石油地球物理勘探,2000,35(6)：751-757.SONG Weiqi,TONG Zhaoqi.Forward finite differential calculation for 3-D transient electromagnetic field[J].Oil Geophysical Prospecting,2000,35(6)：751-757.(in Chinese)

[8]閆 述,陳明生,傅君眉.瞬變電磁場(chǎng)的直接時(shí)域數(shù)值分析[J].地球物理學(xué)報(bào),2002,45(2)：275-283.YAN Shu,CHEN Mingsheng,FU Junmei.Direct time-domain numerical analysis of transient electromagnetic fields[J].Chinese Journal of Geophysics,2002,45(2)：275-283.(in Chinese)

[9]徐凱軍,李桐林.時(shí)域瞬變場(chǎng)電磁場(chǎng)有限差分法[J].世界地質(zhì),2004,23(3)：301-305.XU KAijun,LI Tonglin.A finite-difference time-domain solution for transient electromagnetic fields[J].World Geology,2004,23(3)：301-305.(in Chinese)

[10]朱凱光,林 君,程德福.基于MAT LAB的高頻諧變電磁場(chǎng)的數(shù)值計(jì)算[J].電波科學(xué)學(xué)報(bào),2002,17(3)：224-228.ZHU Kaiguang,LIN Jun,Cheng Defu.Numerical computation of high frequency electromagnetic fields under MA TLAB[J].Chinese Journal of Radio Science,2002,17(3)：224-228.(in Chinese)

[11]NABIGHIAN M N.勘查地球物理-電磁法[M].北京：地質(zhì)出版社,1992：204-207.

[12]ORISTASGLIO M L and HOHMANN W.Diffusion of electromagnetic field into a two-dimentional earth：a finite-difference approach[J].Geophsics,1984,49(7)：870-894.

[13]GUPT ASARMA D,SINGH B.New digital linear filters for Hankel J0 and J1 transforms[J].Geophysical Prospecting,1997,45(5)：745-762.

[14]GUPTASA RM A D.Computation of the time-domain response of a polarizable ground[J].Geophysics,1982,47(11)：953-963.

[15]NABIGHIAN M N.Toward a three-dimensional automatic interpreta-tion of potential field data generalized Hilbert transforms：fundamental relations[J].Geophsics,1984,49(6)：780-786.

[16]華 軍,蔣延生,汪文秉.雙重貝塞爾函數(shù)積分的數(shù)值計(jì)算[J].煤田地質(zhì)與勘探,2001,29(3)：58-62.HUA Jun,JIANG Yansheng,WANG Wenbing.The numerical integration of dual hankel transformation[J].Coal Geology&Exploration,2001,29(3)：58-62.(in Chinese)