航空瞬變電磁正演中頻時轉換方法分析及系數的選取

杜興忠，朱海東

（中國電建集團貴陽勘測設計研究院有限公司，貴陽 550081）

航空瞬變電磁正演中頻時轉換方法分析及系數的選取

杜興忠，朱海東

（中國電建集團貴陽勘測設計研究院有限公司，貴陽 550081）

在瞬變電磁法正演計算中，常常需要頻時轉換來獲得垂直磁場隨時間的變化率。目前已經提出了多種頻時轉換方法，其中G－S變換法只涉及實數計算，具有較快的速度，因此已經成為航空瞬變電磁正演模擬的候選方法。由于航空瞬變電磁數據巨大，時間數較多，因此在保證精度的條件下怎樣提高計算速度成為當前亟待解決的問題。鑒于大量頻點的正弦變換和余弦變換具有較高的精度，這里使用三層和四層典型地電模型，通過比較余弦變換、正弦變換和不同點數G－S變換的正演計算結果來確定最佳的G－S變換系數。實驗結果比較和分析認為，采用8個系數的G－S變換法不僅有較快的計算速度，而且也有較高的計算精度。

航空瞬變電磁法；G－S變換；正弦變換；余弦變換

0 引言

在瞬變電磁法中，常常需要通過頻時轉換來實現時間域的電磁響應計算。通常的頻時轉換法有G －S變換法［1－2］、余弦變換多項式近似法［3］、余弦變換法［4］、正弦變換法［4－5］、數值線性濾波法［6］、滯后褶積法［7］、延遲譜法［8］、混合算法［9－10］等。盡管在上述方法的基礎上，許多研究者也提出了相應的改進方法，羅延鐘［11］基于拉氏變換的延遲定理，建立了一種新的快速Gaver－Stehfest拉氏逆變換；阮百堯［12］提出了改進的Guptasarma算法并用于瞬變電磁的正演計算中；為了提高計算效率，蔣淑芬［13］等提出了一種基于復積分的正余弦變換算法；徐振平等［14］給出了改進的余弦變換方法，但是都未能在計算效率上有顯著改善。昌彥君［15］等通過比較余弦變換多項式近似法、G－S變換法和Guptasarma濾波算法，認為Guptasarma濾波算法計算速度快，但是晚期誤差范圍大，適用的時間范圍最小；余弦變換法精度高，適用的時間范圍最寬，但是耗時最多；G －S變換法的計算精度和耗時介于兩者之間［15］。此外，孫業發［5］等通過比較正弦變換、余弦變換和余弦變換的折線逼近法，最后認為正弦變換在瞬變電磁響應晚期的轉換精度明顯高于其他兩種方法。

由于航空瞬變電磁法數據龐大，涉及該方法的正演計算耗時較大，為此作者以大量頻點（300個）的正弦變換法和余弦變換法為基準，通過與不同系數的G－S變換法的正演計算結果進行比較，確定出G－S變換法的最佳系數，以保證航空瞬變電磁正演計算在擁有較高計算精度的條件下，具有較高的計算速度。

1 航空瞬變電磁法正演計算

航空瞬變電磁法是以飛機為運載工具，發射脈沖電流，在發射機關斷電流時，觀測二次磁場隨時間的衰減特性，進而推斷地下不同介質的分布，從而解決各類地質問題的方法［16］。

圖1 航空瞬變電磁正演模型Fig.1 The forward model of airborne transient electromagnetic

對于圖1所示的航空瞬變電磁正演模型，中心回線裝置的航空瞬變電磁響應計算公式表示為［17］式（1）：

其中：J1表示一階第一類貝塞爾函數；a表示發射線圈半徑；λ是積分變量；RTE表示反射系數，具體為：

并且

式中：ω是角頻率；kn是第n層的波數；σn是第n層的電導率；通常取大地的導磁率μn和自由空間的導磁率μ0相等，在準靜態情況下，通常取u0＝λ。

根據上述公式可以實現航空瞬變電磁法一維正演模擬，其計算步驟如下：

1）設置系統參數和地電模型參數。

2）計算反射系數。

3）漢克爾變換計算。

4）頻時轉換計算。

5）計算電磁響應。

在上述計算步驟，影響電磁響應精度和計算效率的關鍵是漢克爾變換和頻時轉換。

漢克爾變換數值濾波計算方法有很多，這里采用Guptasarma和Singh［18］所提出的數字濾波方法。

漢克爾積分變換的一般形式可以表示為式（6）：

其中：K（λ）表示待變換的函數；Ji表示i階第一類貝塞爾函數；r是一個已給出的定值；λ表示被積函數。數字濾波器有一組濾波系數Wi和兩個常數（位移a和采樣間隔s）。位移a決定了輸入函數采樣的起始點，間隔s決定抽樣時的時間間隔。文獻［18］給出了漢克爾變換的計算公式：

一般常用的漢克爾變換濾波系數有安德森［19］引入的283個系數和801個系數，盡管其精度較高，但是計算速度較慢，所以這里主要采用文獻［18］中給出的濾波系數，即零階有61或120個系數，一階有47或120個系數。

2 頻時轉換

2.1 正弦、余弦變換

根據航空瞬變電磁響應計算公式（1）可以求得Hz。令激發波為階躍波：

對式（8）進行傅氏變換，可得式（9）。

將式（8）代入式（9）再經過傅氏反變換，得到式（10）與式（11）。

根據正余弦變換數值濾波算法［4］，式（10）和式（11）可以分別化為：

2.2 G－S變換

Gaver－Stehfest逆拉氏變換（簡稱G－S變換）是純實數運算，只需對較少的拉氏變換量s值做計算，因而計算速度較快。根據文獻［1］，可以獲得用G－S變換作逆拉氏變換的算法為：對給定的時間t，時間域的響應值f（t），由拉氏變換域中的變量（其中m＝1，2，…，n）的拉氏域響應值F（sm）算出：

表1 G－S變換系數Tab.1 G－S transformation coefficients

式中：n是決定于計算機位數的正偶整數；Km是G －S變換系數，且可由公式（15）計算得出。

式中：m＝1，2，…，n；求和下限i1是的整數部分。

在航空瞬變電磁響應計算中，為了在保證計算精度的條件下獲得較高的計算效率，作者考查濾波系數個數分別取n＝6、8、10、12、14、16時所對應的響應值。根據式（15）求得不同濾波系數個數所對應的G－S變換系數（表1）。

對于G－S變換系數個數n的取值問題，目前存在多種不同的觀點，Wooden［20］等認為n取8或者10比較合適，當n取得更大時，會使得變換變得不穩定，數據結果會發生振蕩；羅延鐘［21］等認為n取16合適；昌彥君和張桂青等［15］則認為n取12合適。作者結合航空瞬變電磁數據特性，通過正演模擬計算確定G－S變換系數個數n的取值，同時和使用正、余弦變換的響應計算結果進行對比分析，驗證n的取值的正確性。

3 電磁響應計算結果及分析

3.1 G－S變換系數對應的電磁響應對比實驗

為了比較不同數量的G－S變換系數對應的電磁響應，設置系統參數如下：發射電流為階躍電流，其歸一化值為1A，發射線圈半徑為1m；接收線圈有效面積歸一化為1m2；飛行高度為30m。設置三層K型地電模型的參數為：第一層電阻率為10Ω· m，層厚為50m，第二層電阻率為100Ω·m，層厚為50m，第三層的電阻率為50Ω·m。在上述系統參數和地電模型的參數設置下，使用表1中不同的G－S變換系數，通過正演計算，獲得的電磁響應如圖2所示。

從圖2中可以看出，不同系數的G－S變換所對應的電磁響應差異主要在晚期，說明G－S變換系數對航空瞬變電磁階躍響應的晚期影響較大，而早期響應幾乎不受影響。此外，6點G－S變換和16點G－S變換的電磁響應在晚期發生振蕩，可以認為6點G－S變換和16點G－S變換不適用于航空瞬變電磁響應計算，實際上，我們在實驗中也發現，小于6點或大于16點的G－S變換系數在電磁響應晚期都發生振蕩，因此對于航空瞬變電磁正演模擬，只余8點、10點、12點和14點G－S變換可以考慮作為候選頻時轉換方法。同時從圖2中也發現，10點、12點、14點G－S變換的階躍響應值幾乎一致，而8點G－S變換的階躍響應值和這3種G －S變換在晚期存在一定差異。

圖2 不同系數個數的G－S變換所對應的階躍響應Fig.2 The step response values corresponding to the G－S transform of different coefficients

對于航空瞬變電磁響應計算，我們在8點、10點、12點和14點G－S變換進行選擇，分別將8點、10點、12點、14點G－S變換和余弦變換的K型地電模型階躍響應值進行對比（圖3）。從圖3中可以看出，在早期，8～14點G－S變換和余弦變換所分別對應的航空瞬變電磁響應的形態和趨勢基本一致，只是幅值上稍有差異；但是在晚期，10～14點G －S變換和余弦變換所對應的電磁響應則有很大差異，而8點G－S變換和余弦變換的電磁響應比較接近，而當頻點足夠多時，余弦變換的晚期響應值更準確和穩定并且具有更寬的時間范圍［15］，這里假設8點G－S變換更能適用于航空瞬變電磁法。

3.2 G－S變換和余弦變換實驗對比

為了驗證8點G－S變換選取的正確性和適用性，這里考慮三層地電模型A型、H型、Q型和四層地電模型AA型、AK型、KH型、KQ型、QQ型、QH型、HA型和HK型，并對比8點G－S變換、10點G－S變換和余弦變換得到的階躍響應值。系統參數的設置同上，地電模型的設置為：A型地電模型的電阻率從上到下依次為1 0Ω·m、5 0Ω·m、1 0 0 Ω·m；H型地電模型的電阻率從上到下依次為50 Ω·m、10Ω·m、100Ω·m；K型地電模型的電阻率從上到下依次為10Ω·m、100Ω·m、50Ω·m；Q型地電模型的電阻率從上到下依次為100 Ω·m、50Ω·m、10Ω·m；AA型地電模型的電阻率從上到下依次為10Ω·m、50Ω·m、100Ω·m、200Ω·m；AK型地電模型的電阻率從上到下依次為10Ω·m、50Ω·m、100Ω·m、50Ω·m；KH型地電模型的電阻率從上到下依次為10Ω·m、50 Ω·m、10Ω·m、50Ω·m；KQ型地電模型的電阻率從上到下依次為10Ω·m、100Ω·m、50Ω·m、10Ω·m；QQ型地電模型的電阻率從上到下依次為200Ω·m、100Ω·m、50Ω·m、10Ω·m；QH型地電模型的電阻率從上到下依次為100Ω·m、50Ω·m、10Ω·m、50Ω·m；HA型地電模型的電阻率從上到下依次為50Ω·m、10Ω·m、50 Ω·m、100Ω·m；HK型地電模型的電阻率從上到下依次為50Ω·m、10Ω·m、50Ω·m、10Ω·m。每種地電模型各層之間的層厚都為50m。在典型的三層地電模型和四層地電模型中所獲得的電磁響應值分別如圖4和圖5所示。

圖3 K型地電模型中G－S變換和余弦變換的階躍響應值Fig.3 The step response values of G－S transform and cosine transform in K Geoelectric model

圖4 8、10點G－S變換和余弦變換在三層地電模型中的階躍響應值Fig.4 The step response values of 8，10points G－S transform and cosine transform in three－layer model

圖5 8、10點G－S變換和余弦變換在四層地電模型中的階躍響應值Fig.5 The step response values of 8，10points G－S transform and cosine transform in four－layer model

圖6 8、10點G－S變換與余弦變換和正弦變換在K型和H型上的階躍響應Fig.6 The step response values of 8，10points G－S transform，cosine transform and sine transform in K－type and H－type geoelectric model respectively （a）K型；（b）H型

在圖4和圖5中，藍色實線為余弦變換得到的階躍響應曲線，綠色實線為8點G－S變換得到的階躍響應曲線，紅色虛線為10點G－S變換得到的階躍響應曲線。從A型、K型、AA型、AK型地電模型圖中可以看出，在晚期的時候綠色實線較紅色虛線明顯更靠近藍色實線，而在其他幾種地電模型圖中綠色實線和紅色虛線幾乎重合，這說明8點和10點G－S變換在一些模型中，其階躍響應基本一致；而在另一些模型中，8點G－S變換所得到的階躍響應在晚期更接近于余弦變換所得到的階躍響應，這就證明了8點G－S變換適用于航空瞬變電磁法，因為其不僅擁有較高的計算精度，同時能獲得更快的計算速度。

3.3 G－S變換、余弦和正弦變換實驗對比

孫業發［5］等認為正弦變換在瞬變電磁響應晚期的轉換精度明顯高于余弦變換和余弦變換的折線逼近法。為了考查8點G－S變換的轉換精度，這里以三層K型和H型地電模型為例，對8點G－S變換、10點G－S變換、余弦變換和正弦變換的電磁響應進行比較。在地電模型的參數設置中，K型模型的電阻率由上至下分別為：10Ω·m、100Ω·m、50 Ω·m，H型模型的電阻率由上至下分別為：100 Ω·m、10Ω·m、100Ω·m，兩種地點模型各層之間的層厚均為50m。在K型和H型地電模型上，8 點G－S變換、10點G－S變換、余弦變換和正弦變換的電磁響應如圖6所示。

從圖6中可以看出，正弦變換得到的響應結果和余弦變換、G－S變換得到的響應結果在早期走勢基本一致，而在晚期它們之間存在一定差異。可以看到在晚期的時候8點G－S變換（綠色實線）比10點G－S變換（粉紅色虛線）更靠近正弦變換（紅色實線）或者余弦變換（藍色實線），驗證了8點GS變換具有較高的計算精度。

從運算時間上來說，在同等計算機平臺上（CPU：i3，內存：3G；操作系統：Windows 7；語言：Matlab 2010b），G－S變換耗時為0.063s，而正弦變換和余弦變換耗時分別為4.768s和4.766s。這說明正、余弦變換盡管有較寬的時間范圍，但是其計算量太大，而G－S變換具有較高的計算效率，這是由于G－S變換系數較少，且不涉及復數運算，而正、余弦變換需要進行復數計算。

綜合上述的計算精度和計算時間分析，可以看到，在計算量較大的航空瞬變電磁法中，8點G－S變換比較適用，因為其不僅具有較高的計算精度，而且還有很高的計算效率。

4 結論

在典型的三層和四層地電模型中，對8點、10點、12點、14點、16點G－S變換以及正弦變換和余弦變換的階躍響應進行對比分析，獲得了如下結論：

1）當需求較高的轉換精度時，正弦變換和余弦變換所需的頻點數較多，因此導致計算量增大，是G －S變換耗時的幾十倍。

2）在不同點數的G－S變換中，8點G－S變換具有較高的轉換精度，僅次于高頻點的正弦變換和余弦變換。

3）在計算效率上，和其他的G－S變換以及正、余弦變換相比，8點G－S變換具有最高的計算效率。

4）由于8點G－S變換在正演計算中不僅具有較快的計算速度，而且也具有較高的計算精度，所以在數據量巨大的航空瞬變電磁正演計算中，應該選擇其作為頻時轉換方案。

［1］ KNIGHT J H，RAICHE A P.Transient electromagnetic calculations using the Gaver－Stehfest inverse Laplace transform method［J］.GEOPHYSICS，1982，47（1）：47－50.

［2］ 樸華榮.電磁測深法原理［M］.北京：地質出版社，1990.

PIAO H R.Principle of electromagnetic sounding method［M］.Beijing：Geological Publishing，1990.（In Chinese）

［3］ 考夫曼A A，凱勒G V.頻率域和時間域電磁測深［M］.地質出版社，1987.

KAUFMAN A A，KELLER G V.Frequency and transient soundings［M］.Geology Press，1987.（In Chinese）

［4］ 王華軍.正余弦變換的數值濾波算法［J］.工程地球物理學報，2004，1（4）：329－335.

WANG HJ.Digital filter algorithm of the sine and cosine transform［J］.Chinese Journal of Engineering Geophysics，2004，1（4）：329－335.（In Chinese）

［5］ 孫業發，李桐林，范翠松，等.中心回線瞬變電磁響應計算的頻－時域轉換方法研究［J］.地球物理學進展，2014，29（3）：1243－1247.

SUN Y F，LI T L，FAN C S，et，al.Research on frequency－time domain transform method in centralloop transient electromagnetic field response computation［J］.Progress in Geophysics，2014，29（3）：1243 －1247.（In Chinese）

［6］ GUPTASARMA D.Computation of the time－domain response of a polarizable ground［J］.GEOPHYSICS，1982，47（11）：1574－1576.

［7］ ANDERSON W L.Fast Hankel transforms using related and lagged convolutions［J］.ACM TRANS.ON MATH.SOFTWARE，1982（8）：344－368.

［8］ NEWMAN G A，HOHMANN G W，ADNERSON W L.Transient electromagnetic response of a three－dimensional body in a layered earth［J］.GEOPHYSICS，1986，51（8）：1608－1626.

［9］ GOLDMAN M M.The integral－finite－difference method for calculating transient electromagnetic fields in a horizontally stratified medium［J］.GEOPHYSICAL PROSPECTING，1983，31：664－686.

［10］李汝傳，王書民，孫鴻雁.用混合算法計算電偶源瞬變電磁響應［J］.物探與化探，2002，26（3）：215－217.

LI R CH，WANG SH M，SUN H Y.The application of mixed algorithm to the Calculation of transient electromagnetic response of couple source.Geophysical ＆Geochemical Exploration，2002，26（3）：215－217.（In Chinese）

［11］羅延鐘，昌彥君.G－S變換的快速算法［J］.地球物理學報，2000，43（5）：684－690.

LUO Y ZH，CHANG Y J.A rapid algorithm for GS transform.Ch inese Journal of Geophysics，2000，43 （5）：684－690.（In Chinese）

［12］阮百堯.Guptasarma算法在瞬變電磁正演計算中的應用［J］.桂林工學院學報，1996，16（2）：167－170.RUAN BY.Applied of Guptasarma algorithm in TEM ＇s forward［J］.Journal of Guilin Institute of Technology，1996，16（2）：167－170.（In Chinese）

［13］蔣淑芬，向淑晃.一種正余弦變換的高效算法［J］.工程地球物理學報，2007，4（5）：512－515.

JIANG SH F，XIANG SH H.An Eifficient Method of Calculating the Sine and Cosine Transform.Chinese Journal of Engineering Geophysics，2007，4（5）：512 －515.（In Chinese）

［14］徐振平，胡文寶.基于快速余弦變換的低頻電磁場數值濾波法正演［J］.石油天然氣學報（江漢石油學院學報），2011，33（1）：63－67.

XU ZH P，HU W B.Low frequency electromagnetic field numerical filtering method based on fast cosine transform［J］.Journal of Oil and Gas Technology，2011，33（1）：63－67.（In Chinese）

［15］昌彥君，張桂青.電磁場從頻率域轉換到時間域的幾種算法比較［J］.物探化探計算技術，1995，17（3）：25 －29.

CHANG Y J，ZHANG G Q.Comparison among three transformation algorithms of electromagnetic field from frequency domain to time domain.Computing Techniques for Geophysical and Geochemical Exploration，1995，17（3）：25－29.（In Chinese）

［16］NABIGHIAN M N.Electromagnetic Methods in Applied Geophysics：Volume 2，Application，Parts A and B［M］.Tulsa：Society of Exploration Geophysics，1991.

［17］NABIGHIAN M N.Electromagnetic Methods in Applied Geophysics：Volume 1，Theory［M］.Tulsa：Society of Exploration Geophysics.1988.

［18］GUPTASARMA D，SINGH B.New digital linear filter for Hankel J0and J1transforms［J］.GEOPHYSICAL PROSPECTING，1997，45（5）：745－762.

［19］ANDERSON W L.Improved digital filters for evaluating Fourier and Hankel transform integrals［R］.U.S.G.S.REPORT，GD－75－012，PAGES，223，1975.

［20］WOODEN B，AZARI M，SOLIMAN M.Well test analysis benefits from new method of Laplace space inversion［J］.OIL ＆GAS JOURNAL，1992，90（29）：108－110.

［21］羅延鐘，張勝業，王衛平.時間域航空電磁一維正演研究［J］.地球物理學報，2003，46（5）：719－724.

LUO Y ZH，ZHANG SH Y，WANG W P.A research on one－dimension forward for aerial electromagnetic method in time domain［J］.Chinese Journal of Geophysics，2003，46（5）：719－724.（In Chinese）

Analysis and coefficients selection of frequency－time conversion methods for airborne transient electromagnetic forward modeling

DU Xing－zhong，ZHU Hai－dong
（Guiyang Hydropower investigation design and Research Institute，CHECC，Guiyang 550081，China）

The change rates of vertical magnetic field with time are usually obtained by the frequency－time conversion method in transient electromagnetic modeling.At present，despite many kinds of frequency－time conversion methods have been proposed，the G－S transform with faster computation speed has become an important candidate for airborne transient electromagnetic method.Due to the huge data and the long time of airborne electromagnetic method，how to improve the computation speed of modeling under the certain accuracy becomes a key problem to be solved.Due to the higher frequency－time conversion precision of sine and cosine transform with the larger number of frequency coefficients，this paper may determine the best coefficients of G－S transform by the comparison of the forward modeling results of sine transform，cosine transform and G－S transform with different frequency based on the three－layer and four－layer model.By comparison and analysis of the experimental results，this paper shows that the G－S transform with 8coefficients not only has faster speed，but also a high degree of accuracy.

airborne transient electromagnetic；G－S transform；sine transform；cosine transform

P 631.3

A

10.3969／j.issn.1001－1749.2015.06.01

1001－1749（2015）06－0671－09

2015－08－30改回日期：2015－09－25

貴州省科學技術基金（［2013］2297）

杜興忠（1973－），男，博士，高級工程師，主要從事工程物探技術研究和應用，E－mail：xingzhongdu＠163.com。