均勻半空間回線源非中心點頻率域電磁場計算

劉 磊，范 濤，李 貅，智慶全

（1．中國煤炭科工集團 西安研究院有限公司，西安 710077；2．長安大學地質工程與測繪學院，西安 710054）

均勻半空間回線源非中心點頻率域電磁場計算

劉 磊1，范 濤1，李 貅2，智慶全2

（1．中國煤炭科工集團 西安研究院有限公司，西安 710077；2．長安大學地質工程與測繪學院，西安 710054）

在回線源頻率域非中心點磁場計算過程中，會遇到雙貝塞爾函數積分問題。這里給出了一種基于靜磁場畢奧－薩伐定律計算雙貝塞爾函數積分的數值方法，將頻率域電磁場寫成由回線源直接激發的一次場和由地下導電介質感應出的二次場的和的形式，根據一次場和二次場被積函數的形態區別，對兩部分的積分采用不同的計算方法。為了驗證該方法的準確性，將回線源剖分成若干個電偶極子，用電偶極子疊加的電磁場與本文提出方法計算的結果對比。對比結果表明，二者的相對誤差很小，這有助于研究回線源電磁場的分布規律。

瞬變電磁法 ；雙貝塞爾函數 ；畢奧－薩伐定律

0 引言

在回線源頻率域電磁場正演模擬研究中，中心點磁場響應特征研究比較詳盡，而非中心點電磁場的研究相對較少。相比于中心點計算表達式，非中心點電磁場的表達式更為復雜，是含有雙貝塞爾函數的積分形式，該類積分由于強震蕩性和積分區間無限大，要求得精確的積分結果十分困難。

對于震蕩函數的積分，Filon［1］采用折線逼近的方法，折線逼近的辦法在積分區間有限時，計算效果較好，當積分區間無限時，折線逼近方法會遇到很大的困難；Iserles等［2］用構造新的函數表達式的方法對截斷點到∞區間的積分進行逼近；陳入云等［3］在Iserles的基礎上，給出了一種改進的漸近算法，并根據誤差理論給出了數值積分的誤差，該算法對于單個貝塞爾函數積分效果較好，但不能完全適用于雙重貝塞爾函數的積分；張爽等［4］利用貝塞爾函數展開形式進行積分計算并計算了誤差；華軍［5］對折線逼近法進行了改進，他將無限的積分區間分為兩個積分區間，在有限的積分區間內利用一般的數值積分方法，在截斷的無限區間內，使用正余弦變換的方法計算積分；Chave［6］利用連分式的辦法解決了有限區間內計算速度和計算精度之間的矛盾，但對于積分區間無限長的情況，效果不佳；薛國強等［7］利用華軍給出的計算方法研究了方形回線瞬變電磁的瞬變電磁場響應特征。無論使用何種算法，積分區間都是無限的，要取得精確的計算結果，需要的計算時間都較長。本研究根據回線源頻率域電磁場的推導過程，結合一次場對應的物理意義，將無限區間的雙貝塞爾函數積分等價成沿發射回線的回路積分，簡化了積分計算過程并提高了計算速度。

1 回線源產生的頻率域電磁場

均勻半空間情況下大回線頻率域垂直磁場表達式為［8］

在磁場表達式的積分中，含有0階和1階貝塞爾函數，被積函數震蕩性很強且積分區間無限，不能直接求取。在頻率域，磁場是由激發電流產生的一次場和由導電介質感應產生的二次場共同構成。從標量赫茲位出發，分別對一次場和感應場表達式進行分析，可得計算總場的簡便方法。均勻半空間回線源頻率域磁場標量赫茲位的表達式為

根據標量赫茲位與磁場的關系

回線源頻率域磁場的垂直分量可以寫成如下的形式

其中：

式中 ：f1是一次場的表達式，由激發電流直接產生；f2是感應場的表達式，由地下渦旋電流產生；a是發射回線的半徑；r是計算點與回線圓心偏離的距離。

要精確計算積分表達式的值，必須了解被積分函數f1、f2的形態，在計算過程中，選取回線半徑a的值為100 m，偏移距離r的值為60 m。一次場和二次場核函數形態如圖1及圖2所示。

一次場核函數是高震蕩函數，不能用常規的方法得到積分結果，而二次場核函數在出現極大值之后迅速衰減，用一般的積分方法可以計算得到。

對定義在區間［0，1］上的連續函數積分，可以用函數在某些點的值與其對應的權值的積之和代替積分結果：

ai、xi分別為對應的高斯系數和高斯點，對于積

圖1 f1被積函數形態Fig．1 The pattern of integration expression f 1

圖2 f2被積函數形態Fig．2 The pattern of integration expression f2

分區間在［a，b］上的積分，通過積分變換將積分區間轉換到［0，1］區間：

在本文的計算過程中，將二次場積分區間截斷，截斷區間為［0，1］，將截斷區間剖分成若干子區間，在每個子區間內采用七點高斯積分，可以比較精確地得到二次場積分結果。

一次場是由發射電流直接激發，而回線源激發的磁場可通過畢奧－薩伐定律計算，雙貝塞爾函數無限區間的積分可以等價地表示成沿線圈的回路積分：

對公式（9）進行進一步化簡得到關于一次場的表達式

從表達式（10）可以看出，積分在a＝r（接收點在線圈所在的位置）時奇異，這是由于不考慮發射回線截面積大小時，激發電流密度無限大，該點處的磁場強度無法確定。

關于一次磁場的高震蕩無窮積分，根據其對應的物理含義，可以用一個比較簡單的定積分公式表示，該積分又可以利用本文所提到的高斯積分公式求取，線圈平面內任意一點的場值都很容易得到，這樣大大簡化了非中心點磁場計算過程，節省計算量。

2 驗證實例

均勻半空間回線源中心點存在解析解，為了驗證本文方法的正確性，選取發射回線半徑為100 m，半空間電阻率為100Ω·m，選定計算頻率范圍為10－2Hz～106Hz，中心點處解析解與本文所提出的數值解對比結果如圖3及圖4所示。

圖3 中心點計算結果對比Fig．3 The contrast of result in the center point

從計算的結果看出，將總場分解成由回線產生的一次場由地下介質產生的二次場的計算結果和解析結果相差很小，在頻率為10－2Hz～106Hz范圍內，解析解和數值解實部、虛部最大相對誤差均小于0．4%。

圖4 中心點場值計算相對誤差Fig．4 The relative error of the result in the center point

當測點位于中心點時，偏移距為0 m，被積函數中只出現一個貝塞爾函數，為進一步驗證本文給出的方法在任何偏移距下都成立，需給出一個比較準確的對比結果。采用文獻［9］中的方法，將回線源剖分成若干段，任意點電磁場的值等于所有偶極子在該點產生的場值的疊加。任意點磁場可以用式（11）表示。

選取發射回線半徑為100 m，偏移距 r＝60 m，半空間電阻率為100Ω·m，兩種方法計算得到的非中心點場值及誤差對比結果見圖5及圖6。

圖5 非中心點計算結果對比Fig．5 The contrast of result in the non－central point

由非中心點的計算結果表明，當計算頻率在10－2Hz～106Hz范圍內時，隨著頻率的增加，基于偶極子疊加計算的結果與本文給出的方法相比，實部和虛部的相對誤差都有所增加，相對誤差最大值小于4%。

圖6 非中心點計算結果實、虛部誤差分布Fig．6 Therelative error distribution of the real and imaginary part

由于與激發源的相對位置不同，磁場在回線源不同位置有不同的變化規律，選取回線半徑為100 m，半空間電阻率值為100Ω·m，激發頻率分別為10－2Hz～106Hz，接收點分別位于0 m、30 m、60 m、90 m、120 m、150 m，利用本文給出的方法，計算了不同位置磁場結果。

圖7 不同偏移距場實部與頻率關系Fig．7 The relationship of the frequency and the real part of magnetic field at different offset

在發射線圈內部，隨著偏移距離的變大，磁場實、虛部隨頻率的變化規律一致，對于實部，磁場幅值逐漸變大，虛部最大值對應的頻率增高（圖7）；在線框的外部，磁場的實部與虛部符號均與線圈內相反，這與激發場在線圈位置處方向發生變化相對應（圖8）。

3 結論

圖8 不同偏移距場值虛部與頻率關系Fig．8 The relationship of the frequency and the imaginary part of magnetic field at different offset

1）回線源中心點與非中心點磁場計算結果表明，可以將回線源頻率域的場分開寫成一次場和二次場的和，根據各自代表的物理含義，頻率域磁場關于一次場的無限區間雙貝塞爾函數積分，可以等價地寫成一種沿發射回線的閉合回路積分形式。這避免了計算震蕩積分和積分區間無限大的問題，大大簡化計算過程，節省了計算時間。

2）利用畢奧－薩伐定律求得一次場后，二次場利用一般的數值積分方法就可以完成，利用本文給出方法與解析解結果及文獻［9］給出的方法對比表明，本文提出的方法是有效的。

3）對非中心點磁場的計算結果表明，非中心點磁場隨激發頻率的變化與中心點有相似的規律，在發射線圈內，非中心點場值與中心點場值實部表現為幅值的差異，虛部表現為最大值對應的頻率發生移動，在發射線圈外部，實部與虛部的符號均與線圈內相反。

［1］ FILON L N G．On a quadrature formula for trigonometric integrals［J］．Proc Royal Soc Edinburgh，1928（49）：38－47．

［2］ ISERLES A，NORSETT S P．On the numerical quadrature methods for highly oscillatory integrals and their implementation［J］．BIT Numerical mathematic，2004（44）：755－772．

［3］ 陳入云，向淑晃．一類含貝塞爾函數積分的數值算法［J］．重慶工學院學報：自然科學版，2008，22（11）：83－88．

［4］ 張爽，郭欣，宋立軍．利用貝塞爾函數的級數形式進行數值計算的誤差分析［J］．吉林大學學報，2004，14（2）：57－59．

［5］ 華軍，蔣延生，汪文秉．雙重貝塞爾函數積分的數值計算［J］．煤田地質與勘探，2001，29（3）：58－61．

［6］ ALAN D．Chave，numerical integration of related Hankel transforms by quadrature and continued fraction expansion［J］．Geophysics，1983，48（12）：1671－ 1686．

［7］ 薛國強，李貅，郭文波，等．大回線源瞬變電磁場響應特征［J］．石油地球物理勘探，2007，42（5）：586－590．

［8］ 李貅．瞬變電磁測深的理論與應用［M］．陜西：陜西科學技術出版社，2002．

［9］ 李建平，李桐林，張亞東．層狀介質任意形狀回線源瞬變電磁場正反演［J］．物探與化探，2012，36（2）：256－258．

A new calculation method for double bessel function integration

LIU Lei1，FAN Tao1，LI Xiu2，ZHI Qin－quan2
（1．Xi'an Research Institute of China Coal Technology＆Engineering Group Corp，Xi'an 710077，China；2．College of Geology Engineering and Geomatics，Chang'an University，Xi'an 710054，China）

According to the double Bessel function integration account in the process of calculation non－center point transient electromagnetic field，a new numerical method based on Biot Savart Law was offered in this paper．We describe the field as the summation of primary field stimulate by the loop and secondary field generated by the media，after analysis the pattern of the kernel function，calculated the two integral function with respective method．In order to confirm the validity，subdivision the loop into some electric dipole，comparing the result acquire by superposition every electric dipole to that gain from the method offered，it prove the effectiveness of the new method．

TEM；double Bessel function；Biot Savart law

P 631．3＋25

A

10．3969／j．issn．1001-1749．2014．05．07

1001-1749（2014）05-0555-05

2014-02-24 改回日期：2014-05-23

“十二五”國家科技重大專項（2011ZX05040-002）；國家自然科學基金項目（41104087）

劉磊（1988－），男，碩士，從事煤田地質勘探方面研究，Ej-mail：chdliulei＠163．com。