戴世坤 曾 鈴* 周印明③ 李 昆 陳輕蕊 凌嘉宣

(①中南大學(xué)地球科學(xué)與信息物理學(xué)院,湖南長沙 410083; ②中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室，湖南長沙 410083； ③東方地球物理公司綜合物化探處，河北涿州 072751)

0 引言

在電磁場地球物理方法中，均勻?qū)訝罱橘|(zhì)電磁場解析表達(dá)式含有漢克爾積分[1-2]，其核函數(shù)為0階和1階Bessel函數(shù)。隨著其宗量的增大，Bessel函數(shù)呈現(xiàn)快速振蕩和慢衰減的特征，造成漢克爾積分難以實(shí)現(xiàn)高效、高精度的計(jì)算，特別是高頻和大收發(fā)距情況下難度更大。為了解決這個(gè)問題，學(xué)者們提出了許多方法，主要分為兩大類： 第一類為數(shù)字濾波方法； 第二類為直接積分法(Direct integration method，DIM)。Ghosh[3]首先將數(shù)字濾波法引入地球物理數(shù)值計(jì)算； 針對(duì)電磁法，Johansen等[4]提出了漢克爾數(shù)字濾波系數(shù)的解析計(jì)算方法。之后，很多研究者給出了不同優(yōu)化濾波系數(shù)[5-6]。蔡盛[7]對(duì)比了5組高精度快速漢克爾長濾波系數(shù)的計(jì)算精度，研究結(jié)果表明在一定的收發(fā)距和頻率范圍內(nèi)，計(jì)算精度較高，但耗時(shí)太大，隨著濾波系數(shù)的增大而增長?？傮w來說，數(shù)字濾波法適用于中、低頻和適中的收發(fā)距情況，在高頻和大收發(fā)距的情況下計(jì)算精度很低。考慮到數(shù)字濾波法的局限性，眾多學(xué)者利用多項(xiàng)式精確擬合核函數(shù)的思想，達(dá)到求解解析解的目的，提出了直接積分類方法。代表性方法包括高斯積分法[8]、基于快速正弦余弦算法的漢克爾變換算法[9]、直接法[10]、離散復(fù)鏡像法[11-12]、譜域法[13]、基于高階窗函數(shù)的結(jié)合非線性變換的算法[14]、外推積分方法(QWE)[15]、三次樣條插值法[16]等。其中最后一種方法采用三次樣條函數(shù)對(duì)核函數(shù)進(jìn)行插值，得到一系列核函數(shù)為多項(xiàng)式的簡單Sommerfeld積分，再利用Bessel函數(shù)的遞推公式和Lommel公式的漸進(jìn)展開求解。這些方法雖能較好地應(yīng)用于低頻和適中收發(fā)距的均勻?qū)訝罱橘|(zhì)電磁場的計(jì)算，但隨著頻率和收發(fā)距的增大，這些方法的計(jì)算精度逐漸降低，甚至失效。針對(duì)高頻電磁場的計(jì)算，鄭圣談等[17]提出了高密度采樣數(shù)字濾波法，但當(dāng)介質(zhì)導(dǎo)電率變小、介電常數(shù)變大、收發(fā)距變大、頻率增加時(shí)，該方法計(jì)算精度會(huì)大幅度降低，不具普適性。

針對(duì)電磁場漢克爾積分核函數(shù)復(fù)雜的特點(diǎn)，本文提出一種高效、高精度直接積分方法，用于頻率和收發(fā)距廣普適用的電磁場高效、高精度的數(shù)值計(jì)算。其基本思想是Bessel函數(shù)可以在兩個(gè)區(qū)間分別用不同的多項(xiàng)式展開，每一區(qū)間的漢克爾積分可離散為多個(gè)單元積分之和，每個(gè)單元被積函數(shù)采用三次樣條插值函數(shù)表示，由此可求得積分的解析解，并通過求和獲得漢克爾積分的數(shù)值解。在此基礎(chǔ)上，利用均勻全空間電偶極子電磁場的解析解，正確選取積分范圍，并合理剖分積分單元。數(shù)值解與解析解對(duì)比表明，本文算法正確、可靠。與數(shù)字濾波類算法的比較表明，本文算法廣泛適用于不同頻率和不同收發(fā)距電磁場的計(jì)算，具有較強(qiáng)的普適性。

1 方法原理

均勻?qū)訝罱橘|(zhì)電磁場的解析解可用漢克爾積分表示

(1)

式中：g(m)為大地響應(yīng)函數(shù)，m表示積分參數(shù)；r為發(fā)射源到接收點(diǎn)的水平距離，mr為宗量；Jn(·)為n階Bessel函數(shù)。

式(1)中，Bessel函數(shù)為0階和1階，根據(jù)特殊函數(shù)手冊[18]，可以分別在兩個(gè)區(qū)間用不同的多項(xiàng)式展開。第一區(qū)間為0≤mr≤4，Bessel函數(shù)可展開為簡單的代數(shù)多項(xiàng)式，每個(gè)積分單元內(nèi)代數(shù)多項(xiàng)式與響應(yīng)函數(shù)的乘積可用三次樣條插值函數(shù)表示。第二區(qū)間為4

在第一個(gè)區(qū)間(0≤mr≤4)有

J0(mr)=-0.0005014415t7+0.0076771853t6-

0.0709253492t5+0.4443584263t4-

1.7777560599t3+3.9999973021t2-

3.9999998721t+1.0+e0(mr)

(2)

式中截?cái)嗾`差e0≤1×10-8。

在第二個(gè)區(qū)間(4

(3)

式中

(4)

其中截?cái)嗾`差e1、e2≤1×10-8。

對(duì)第一區(qū)間0≤mr≤4，將積分范圍均勻剖分為N個(gè)單元。由式(2)可知，Bessel函數(shù)的展開式為簡單的代數(shù)多項(xiàng)式，因此設(shè)f(m)=g(m)J0(mr)，這個(gè)區(qū)間的積分F1可表示為

(5)

式中mi(i=1,2,…,N,N+1)為離散采樣點(diǎn)。每個(gè)單元內(nèi)被積函數(shù)fi用三次樣條插值函數(shù)[18-19]表示為

fi=di(m-mi)3+ci(m-mi)2+

bi(m-mi)+ai

(6)

式中ai、bi、ci、di為單元i的三次樣條插值函數(shù)的系數(shù)。設(shè)Li=mi+1-mi為單元i的長度，單元i的積分解析表達(dá)式為

(7)

因此，第一區(qū)間各單元積分累加的結(jié)果為

(8)

對(duì)于第二區(qū)間4

(9)

設(shè)

(10)

(11)

(12)

式中

(13)

因此，第二區(qū)間各單元積分累加求和結(jié)果為

(14)

最終漢克爾積分為

F=F1+F2

(15)

2 積分參數(shù)的選取

為了提高計(jì)算精度和效率，需合理選取積分范圍，恰當(dāng)?shù)剡M(jìn)行單元剖分，而這些參數(shù)的選取依賴于響應(yīng)函數(shù)的特征和變化規(guī)律[19-20]。本文以均勻全空間電偶極子源電磁場為例，通過研究不同響應(yīng)函數(shù)的變化規(guī)律，確定積分范圍，并選定單元剖分方案。

在均勻全空間中，將x方向上電偶極子源產(chǎn)生的電、磁場用漢克爾積分表示，其解析解的漢克爾積分表達(dá)式分別為

(16a)

(17)

(18)

取電導(dǎo)率σ=0.0001S/m，相對(duì)介電常數(shù)ε=1，電偶極子源中心坐標(biāo)為(0,0,0)，偶極矩為1C·m，頻率v=1×108Hz， 積分變量m在[1×10-12,1×108]區(qū)間按對(duì)數(shù)等間距均勻采樣，每一個(gè)對(duì)數(shù)單位內(nèi)采樣點(diǎn)為5000個(gè)。計(jì)算上述三種響應(yīng)函數(shù)的值，結(jié)果如圖1所示。

圖1 三種響應(yīng)函數(shù)隨m的變化曲線

由圖1可以看出，m從1×10-6增大到1×10-1時(shí)，三種響應(yīng)函數(shù)的對(duì)數(shù)值呈線性增大； 當(dāng)m為1×10-1～1×101時(shí)，三種響應(yīng)函數(shù)值隨m變化較快，尤其在1×100附近出現(xiàn)尖脈沖；m>102時(shí)三種響應(yīng)函數(shù)隨著m的增大按照對(duì)數(shù)以近似線性規(guī)律迅速衰減； 在m接近1×103時(shí)，三種函數(shù)的值均衰減至低于10-20。根據(jù)上述分析結(jié)果，m的積分范圍可選取為[0,1×104]，并整體上按照對(duì)數(shù)等間距規(guī)律進(jìn)行單元剖分。特別地，在m為1×10-1～1×102時(shí)，因響應(yīng)函數(shù)值變化較快，單元剖分可適當(dāng)加密，尤其在m接近1×100時(shí)，響應(yīng)函數(shù)變化劇烈，單元剖分需格外加密。

根據(jù)以上原則，對(duì)均勻全空間電偶極子源產(chǎn)生電磁場用漢克爾積分表達(dá)式進(jìn)行數(shù)值計(jì)算，并與解析解進(jìn)行比較，以驗(yàn)證算法的正確性。下文分別計(jì)算均勻全空間中1×104、1×108Hz的電場響應(yīng)。

對(duì)v=1×104Hz，計(jì)算zr=40m所在平面的電磁場，計(jì)算范圍：x方向?yàn)?500～500m，y方向?yàn)?500～500m，m的采樣區(qū)間為0～1×104。在區(qū)間1×10-3～1×101加密積分單元，總采樣點(diǎn)數(shù)為373。圖2為電場分量的數(shù)值解、解析解及相對(duì)誤差平面圖。由圖可見，數(shù)值解與解析解曲線形態(tài)基本一致，相對(duì)誤差都很小，最大為8×10-6(Ex分量)，即便在源附近，其誤差也不大于1%，表明本文算法的正確性和可靠性，且精度較高。

圖2 v=1×104Hz時(shí)電場分量Ex(a)、Ey(b)、Ez(c)的解析解(左)、數(shù)值解(中)及其相對(duì)誤差(右)

對(duì)v=1×108Hz，計(jì)算zr=10m所在平面的電場。計(jì)算的平面范圍x方向?yàn)?15～15m，y方向?yàn)?15～15m，m的采樣區(qū)間為0～1×105，在區(qū)間1×10-2～1×100加密采樣點(diǎn)，總采樣點(diǎn)數(shù)為811。圖3為該模型的電場分量數(shù)值解、解析解及相對(duì)誤差平面圖。從圖3可看出，數(shù)值解與解析解曲線形態(tài)基本一致，相對(duì)誤差很小，其中Ey分量的誤差最大，最大為4×10-6，即便如此，在源附近的誤差也都不大于1%，這表明算法在高頻時(shí)計(jì)算的場值是正確、可靠的，且精度較高。

圖3 v=1×108Hz時(shí)電場分量Ex(a)、Ey(b)、Ez(c)的解析解(左)、數(shù)值解(中)及其相對(duì)誤差(右)

上述兩個(gè)頻率的計(jì)算結(jié)果證明了本文算法正確、可靠，且能適應(yīng)高頻電磁場的計(jì)算。

3 模型算例

設(shè)計(jì)均勻全空間模型和均勻?qū)訝罱橘|(zhì)模型，分別使用本文DIM算法和常用的數(shù)字濾波類算法模擬電磁場，以驗(yàn)證本文算法對(duì)不同頻率、不同收發(fā)距電磁場模擬計(jì)算的廣普適用性。

3.1 均勻全空間模型

首先驗(yàn)證本文算法對(duì)于頻率的廣普適用性。

在均勻全空間中，設(shè)發(fā)射源到接收點(diǎn)的水平距離R=100m，計(jì)算頻率范圍為1×10-2～1×1010Hz，電導(dǎo)率σ=0.0001S/m，相對(duì)介電常數(shù)ε=1，電偶極子源中心坐標(biāo)為(0,0,0)，偶極矩為1C·m。分別使用本文算法與四種經(jīng)典的數(shù)字濾波算法計(jì)算該模型的數(shù)值解與解析解的相對(duì)誤差，其中四種數(shù)字濾波算法的濾波系數(shù)選取見表1。數(shù)字濾波方法①和②參見文獻(xiàn)[21]，方法③和方法④分別參見文獻(xiàn)[22]和文獻(xiàn)[23]。這四種方法及本文提出的直接積分法計(jì)算結(jié)果如圖4所示，其中黑色直線表示相對(duì)誤差為1%。從圖4可看出，在頻率v<1×106Hz的情況下，這四種數(shù)字濾波類算法的計(jì)算誤差均小于1%；隨著頻率的升高，四種濾波系數(shù)類算法的誤差越來越大，而本文直接積分法的計(jì)算誤差小于1%。表明本文提出的直接積分法計(jì)算精度高，在廣泛的頻率范圍內(nèi)適用。

表1 四種常用的數(shù)字濾波算法的濾波系數(shù)長度取值

為驗(yàn)證該方法對(duì)收發(fā)距的廣普適用性，在均勻全空間中，取頻率v=1×107Hz， 發(fā)射源到接收點(diǎn)的水平距離r范圍為1×10-3～1×103m，其他參數(shù)不變。使用本文算法與四種經(jīng)典的數(shù)字濾波算法分別計(jì)算不同收發(fā)距時(shí)的數(shù)值解與解析解的相對(duì)誤差，其中四種數(shù)字濾波算法的濾波系數(shù)選取同上(表1)。計(jì)算結(jié)果見圖5?？梢钥闯?，在收發(fā)距r<2m時(shí)，直接積分法和四種數(shù)字濾波類算法的計(jì)算誤差均小于1%； 但隨著收發(fā)距的增大，四種濾波系數(shù)類算法的誤差越來越大，而直接積分法在任意r時(shí)，計(jì)算誤差均不大于1%。表明本文提出的直接積分法計(jì)算精度高，在廣泛的收發(fā)距范圍內(nèi)適用。

圖4 均勻全空間模型不同頻率下直接積分法與四種

圖5 均勻全空間模型不同收發(fā)距下直接積分法與四種

3.2 層狀介質(zhì)模型電磁場不同算法對(duì)比

設(shè)計(jì)一個(gè)三層層狀介質(zhì)模型，模型參數(shù)見圖6。電偶極子源中心坐標(biāo)為(0,0,0)，偶極矩為1C·m。

圖6 三層層狀模型斷面示意圖

計(jì)算頻率v=1×108Hz、z=10m所在平面的電場分量，計(jì)算平面范圍x方向和y方向均為-15～15m。

圖7和圖8分別為表1中第④種數(shù)字濾波法與本文直接積分法計(jì)算的Ex、Ey、Ez實(shí)部與虛部平面等值線圖。由圖可見，數(shù)字濾波法計(jì)算的Ex、Ey、Ez等值線圓滑度較低，表明計(jì)算結(jié)果有一定誤差，而直接積分算法計(jì)算的電場等值線圓滑度較高，表明計(jì)算結(jié)果更準(zhǔn)確、精度更高。

圖7 三層層狀介質(zhì)模型直接積分法(上)與數(shù)字濾波法(下)計(jì)算的電場分量Ex(a)、Ey(b)和Ez(c)實(shí)部平面等值線圖

圖8 三層層狀介質(zhì)模型直接積分法(上)與數(shù)字濾波法(下)計(jì)算的電場分量Ex(a)、Ey(b)和Ez(c)虛部平面等值線圖

4 結(jié)論

本文提出一種高效、高精度的直接積分方法，適用于不同頻率和不同收發(fā)距的電磁場模擬計(jì)算，尤其適合高頻電磁場的計(jì)算，具體取得以下結(jié)論。

(1)利用Bessel函數(shù)可以在兩個(gè)積分區(qū)間分別用不同的多項(xiàng)式展開。將每一區(qū)間的漢克爾積分離散成多個(gè)單元積分之和，每個(gè)單元被積函數(shù)采用三次樣條插值函數(shù)表示，由此可求得積分的解析解。通過疊加，求得漢克爾積分的數(shù)值解。

(2)響應(yīng)函數(shù)的特征和變化規(guī)律是積分參數(shù)選取的關(guān)鍵因素。為此，以均勻全空間電偶極子源電磁場為例，研究了三種典型的響應(yīng)函數(shù)隨積分變量的變化規(guī)律。計(jì)算結(jié)果表明，三種響應(yīng)函數(shù)隨積分變量呈對(duì)數(shù)規(guī)律變化，據(jù)此按對(duì)數(shù)規(guī)律確定積分范圍和單元剖分：在響應(yīng)函數(shù)值變化快的地方，單元剖分可適當(dāng)加密，尤其在尖脈沖附近，單元剖分需格外加密；在響應(yīng)函數(shù)值變化慢的地方，單元剖分可適當(dāng)稀疏。這樣合理剖分積分單元可提高直接積分法的計(jì)算精度和效率。均勻全空間電偶極子源電磁場數(shù)值解與解析解對(duì)比表明本文算法更正確、可靠。

(3)分別設(shè)計(jì)均勻全空間模型和層狀介質(zhì)模型，對(duì)本文提出的直接積分法與數(shù)字濾波方法電磁場數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果表明本文積分法對(duì)頻率和收發(fā)距具有廣普適用性，尤其對(duì)高頻電磁場，其優(yōu)勢更加突出，為高頻電磁場的計(jì)算提供了重要的方法基礎(chǔ)。

附錄A 1階Bessel函數(shù)漢克爾積分的解析表達(dá)式推導(dǎo)過程

根據(jù)特殊函數(shù)手冊[18]，對(duì)J1(mr)函數(shù)分兩個(gè)區(qū)間進(jìn)行展開。

在第一個(gè)區(qū)間(0≤mr≤4)

0.0236616773t5+0.1777582922t4-

0.8888839648t3+2.6666660544t2-

3.9999999710t+1.9999999998)+e0(mr)

(A-1)

其中截?cái)嗾`差e0≤1×10-8。

在第二個(gè)區(qū)間(4

(A-2)

式中

(A-3)

0.000999941t-3+0.000266891t-2-

0.001601836t-1+0.093749994)+e2(mr)

(A-4)

其中截?cái)嗾`差e1、e2≤1×10-8。

對(duì)于第一區(qū)間(0≤mr≤4)，將其剖分為N個(gè)單元，由式(A-2)可知，Bessel函數(shù)的展開式為簡單的代數(shù)多項(xiàng)式，因此設(shè)f(m)=g(m)J1(mr)，這部分的積分為

(A-5)

設(shè)mi(i=1,2,…,N,N+1)為離散采樣點(diǎn)，每個(gè)單元內(nèi)被積函數(shù)fi可采用三次樣條插值函數(shù)表示為

fi=di(m-mi)3+ci(m-mi)2+

bi(m-mi)+ai

(A-6)

(A-7)

第一區(qū)間各單元積分累加表達(dá)式為

(A-8)

對(duì)第二區(qū)間(4

(A-9)

設(shè)

(A-10)

在每個(gè)單元內(nèi)，F(xiàn)p、Fq變化規(guī)律簡單，可采用三次樣條插值函數(shù)近似表示為

(A-11)

(A-12)

式中

(A-13)

可求得第二區(qū)間各單元的積分為

(A-14)

據(jù)此，漢克爾積分為

F=F1+F2

(A-15)