柳建新，張 維*，曹創(chuàng)華，蔡 盛

(1．中南大學 地球科學與信息物理學院，長沙 410083；2．有色資源與地質災害探查湖南省重點實驗室，長沙 410083；3．中南大學 地質調查研究院，長沙 410083)

0 引言

瞬變電磁法的主要工作裝置有中心回線、重疊回線、大定源回線等，其中又以大定源回線裝置應用最為廣泛[1-2]。大回線源裝置的工作方法是采用較長的發(fā)送回線,利用階躍波電流場源激勵,在大地產生過渡過程場,斷電瞬間在大地中形成渦旋交變電磁場。在回線內外一定的范圍內,可觀測到這種由地下介質產生的二次感應電磁場,隨時間變化的衰減特性,進而從測量得到的異常信號中分析出地下不均勻體的導電性能和位置,以達到解決地質問題的目的。大回線裝置與目標體具有最佳耦合、異常幅值大、形態(tài)簡單、受旁側地質體影響小的特點,并對高阻層有很強的穿透能力,對低阻層有較高的分辨能力。這種裝置對鋪設回線的要求不是很嚴格,一旦鋪好回線后,不僅可采用多臺接收機同時工作,而且可以把接收線圈排成陣列,發(fā)展成陣列式接收的觀測系統。這種場源具有發(fā)射磁矩大、場均勻及隨距離衰減慢等特點,適合于密集點距采樣,精細探測。

大回線由于面積較大，不能看作磁偶極子。因此，研究大回線源激發(fā)的電磁場有著重要的意義。國內、外對大回線源產生電磁場的研究并不多，翁愛華[3]等首先求出任意一層電磁波的反射系數和透射系數，然后據此計算了圓形發(fā)射回生的頻率域電磁場；M. Poddar[4]根據電磁學互易原理和垂直磁偶源產生的電磁場；求解出距形發(fā)射回線產生的頻率域電磁場；Nagendra Partap Singh等[5]利用超幾何分布函數，化簡了圓形發(fā)射回線電磁場表達式中貝塞爾函數的乘積，從而求解出頻率域電磁場。

作者從基本的電磁場理論出發(fā)，給出大回線發(fā)射源框內外不同位置的電磁響應表達式，利用一種新的算法計算電磁響應，并與傳統的線性數值濾波法計算結果進行了對比分析(對比所用的濾波權系數，采用Anderson 計算的801個系數)。

1 理論基礎

ε算法是一種遞歸方法，主要是基于以下三組關系式[7-8]：

每個區(qū)間上的積分，采用以下形式計算：

式中m是積分階數；w是與橫坐標x有關的權重向量。

通過對上述方程做簡單的變換，可以變換成與FHT算法一樣的形式：

2 回線源瞬變電磁響應理論分析

圖1為大回線源瞬變電磁法工作裝置圖。置于地表的發(fā)送線框激發(fā)的階躍波，會在地下形成渦流場，回線電流的通斷過程中，在回線內的磁通量會發(fā)生變化，通過分析接收線框接收到的磁通量變化信息，可以獲知地下介質的分布情況。

圖1 大定源回線瞬變電磁工作布置示意圖Fig．1 Working schematic of transient electromagnetic method in large loop source

在層狀介質情況下，由階躍電流激發(fā)，在發(fā)送線框內任意位置接收的頻率域電磁響應為[9]：

其中Hr(ω)、Eφ(ω)為徑向的磁場分量和電場分量；Hz(ω)為垂向磁場分量；a是回線等效半徑；I是發(fā)送電流；z0是表層阻抗；z(1)是總的波阻抗；h為發(fā)送線框距地面高度；r為發(fā)送線圈中心點與接收線框中心點之間的距離。波阻抗由以下公式計算[10]得到：

Hn為層狀介質第n層的厚度。

z(n)=znzj=-iωμ0/uj

z0=-iωμ0/u0=-iωμ0/λ

3 計算結果分析

計算都是基于均勻層狀介質模型進行的。ρ1=100 Ω·m，ρ2=1 000 Ω·m，a=100 m，H=40 m，發(fā)送電流為I=1 A，發(fā)送頻率為f= 25 Hz。

圖2和圖3描述了積分核函數的收斂情況。由圖2可知，對于h= 0 m的情況，零階和一階的第一類貝塞爾函數在λ增大的時候，振蕩很明顯，這是貝塞爾函數本身的性質。除去貝塞爾函數之后，核函數剩余部分的虛部是單調遞減的，而實部則是單調遞增的。在h= 5 m 時(圖3)，除去貝塞爾函數之后的核函數虛部仍然是單調遞減的，而實部則是隨λ增大而增大，達到一個最大值后開始隨λ增大而減小。

由Anderson 的研究可知，在利用快速漢克爾變換計算貝塞爾函數積分的時候，積分核函數(除去貝塞爾函數)必須是單調遞減的。所以對于h= 0 m的情況，用漢克爾變換不能得到穩(wěn)定的解(圖4)。

對于垂向磁場實部的計算：在h= 0 m時(圖4)，利用FHT算法得到的結果跳動比較厲害，而且在r>a時，計算結果誤差較大，而用QWE法得到的結果除了在r= 100 m附近出現一個極大值外，在框內、外的計算結果都比較穩(wěn)定；在h= 5 m時(圖5)，兩種方法的計算結果曲線趨勢是一樣的，不過用QWE算法得到的計算結果相對較穩(wěn)定。對于虛部的計算，兩者的計算結果吻合的都非常好。通過上面的分析可知，這是因為核函數的虛部是收斂，所以利用FHT算法計算能得出比較精確的結果。

此外，由于QWE計算積分時的求積點是不均勻分布的一些高斯積分點，計算的核函數的數量是隨r的變化而變化的，而數值濾波法由于濾波權系數數量是固定的(801個)，所以計算的核函數的數量是不隨r變化的。由圖6可以看出，對于大多數的r，用QWE算法所需計算的核函數數量都是少于用FHT算法的。

圖2 r =50 m，h =0 m 時，積分核函數(不含貝塞爾函數)及貝塞爾函數性質Fig．2 Characteristics of kernel functions (without the Bessel functions) and Bessel functions when r =50m, h=0 m

圖3 r =50 m，h =5 m 時，積分核函數(不含貝塞爾函數)及貝塞爾函數性質Fig．3 Characteristics of kernel functions (without the Bessel functions) and Bessel functions when r =50 m, h =5 m

圖4 h=0 m時，兩種方法在不同偏移距條件下實部和虛部響應對比圖Fig．4 Comparison of real and imaginary responses by each method in different offset, h=0 m

圖5 h = 5 m時，兩種方法在不同偏移距條件下實部和虛部響應對比圖Fig．5 Comparison of real and imaginary responses by each method in different offset, h=5 m

4 結論

作者對大回線源瞬變電磁法一維層狀正演算法進行了研究，將QWE方法與傳統的基于快速漢克爾變換的數字濾波法進行了比較：兩者對于虛部的計算吻合比較好；然而對于實部的計算，QWE算法計算的結果相對較穩(wěn)定，尤其是對于線框位于地表的情況，FHT算法對于不同偏移距的計算結果變化較大，而QWE算法的結果比較穩(wěn)定；此外QWE算法所需計算的核函數數量隨偏移距的不同是不斷變化的，在大多數偏移距條件下，所需計算的核函數數量均少于數字濾波法。

圖6 h =0 m時，兩種方法不同偏移距條件下核函數計算數量對比Fig．6 Comparison of the number of kernels evaluations required by each method in different offset

參考文獻：

[1] 薛國強, 秦克章, 黃樹峰,等.大回線源瞬變電磁技術在西藏山南地區(qū)探礦中的應用[J].地質與勘探,2011, 47(1):100-106.

[2] 徐貴東, 黃繼軍, 劉德軍,等.瞬變電磁在吉林省紅旗嶺銅鎳礦勘查中的應用[J].吉林地質, 2010, 29(4): 79-82.

[3] 翁愛華,李舟波,王雪秋.地表大回線源在任意層狀介質中產生磁場的計算[J].物探化探計算技術, 2000, 22(3):245-248.

[4] PODDARM A. Rectangular Loop Source of Current on a Two-layered Earth[J]. Geophys,Prosp, 1982, 30:101-114.

[5] NAGENDRA PARTAP SINGH, TORU MOGI. Electromagnetic Response of a Large Circular Loop Source on a Layered Earth: A New Computation Method[J]. Pure and Applied Geophysics,2005, 162:181.

[6] LONGMAN I M. Note on a method for computing infinite integrals of oscillatory functions: Mathematical Proceedings of the Cambridge Philosophical Society[J].Mathematical proceedings of the Cambrige Philosophical Society,1956,52(4): 764-768.

[7] SHANKS D. Nonlinear transformations of divergent and slowly convergent sequences[D]．Journal of Mathematical Physics, 1955, 34：1-42.

[8] KERRY KEY.Is the fast Hankel transform fast than quadrature[J]．Geophysics, 2012,77(3): 21-30.

[9] 考夫曼A.A, 凱勒.G.V.頻率域和時間域電磁測深[M].北京:地質出版社,1987.

[10] 牛之璉.時間域電磁法原理[M].長沙：中南大學出版社, 2007.

[11] ANDERSON W L. Fast Hankel-transfroms using related and lagged convolution[J]．ACM Transforms on Mathematical Software, 1982, 8(4)：344-368.

[12] GUPTASARMA D B. Singh. New digital linear filters for Hankel J(0)and J(1)transforms[J]. Geophysical Prospecting, 1997, 45:745-762.

[13] 方文藻,李予國,李貅.瞬變電磁測深法原理[M].西安:西北工業(yè)大學出版社,1993.

[14] 李貅.瞬變電磁測深的理論與應用[M].西安：陜西科學技術出版社, 2002.

[15] 薛國強,李貅,郭文波.大回線源瞬變電磁響應[J].石油地球物理勘探,2007,42(5):586-590.

[16] CHAVE A D. Numerical integration of related Hankel transforms by quadrature and continued fraction expansion[J].Geophysics, 1983, 48:1671-1686.

[17] 劉繼東, 方文藻. 用線性數字濾波法計算大回線源在地下形成的瞬變電磁場[J].物探化探計算技術, 1996, 18(4):231-237.

[18] 李建慧, 劉樹才, 李富,等.大定源瞬變電磁法矩形發(fā)射回線激發(fā)的電磁場[J]．物探化探計算技術, 2008, 30(2):154-157.