時間域航空電磁系統(tǒng)瞬變?nèi)珪r響應正演模擬

殷長春，黃 威，賁 放

吉林大學地球探測科學與技術學院，長春 130026

1 引 言

航空電磁法（Airborne Electromagnetics，簡稱AEM），是用來快速地質填圖、普查良導電金屬礦、地下水及環(huán)境和工程監(jiān)測等領域的航空物探方法［1］.該方法優(yōu)勢在于速度快、成本低、通行性好，尤其是在地形地質條件復雜地區(qū)，比如沙漠地區(qū)、山區(qū)、森林、沼澤，以及大量植被覆蓋地區(qū)，可以達到一般物探所不能達到的效果［2］.航空電磁法根據(jù)發(fā)射波形的不同分為時間域和頻率域系統(tǒng).

時間域直升機吊艙系統(tǒng)是以直升機為載體，在飛行過程中通過安裝在吊艙上的發(fā)射線圈，向地下發(fā)送一次脈沖磁場，即為一次場；在其激發(fā)下，地下的地質體中激發(fā)出感應渦流，將產(chǎn)生隨時間變化的感應二次場.通過采集地下返回的電磁信號，可達到探測地下地質體的目的［3］.該方法由于直升機飛行速度慢，數(shù)據(jù)采樣密集，具有很高的橫向分辨率，可以較好地探測地下地質構造.由于時間域系統(tǒng)進行純異常觀測，且發(fā)射功率很大，可探測較大深度的勘探目標.

國外的航空系統(tǒng)目前已較發(fā)達，加拿大Fugro航空地球物理公司的HeliGEOTEM系統(tǒng)，Geotech公司的VTEM系統(tǒng)，Aeroquest公司的AeroTEM系統(tǒng)，丹麥Aarhus大學的SkyTEM系統(tǒng)等已廣泛應用于地質填圖、地下水和礦產(chǎn)資源調(diào)查［4］.有關航空電磁系統(tǒng)響應的理論和應用研究國外也起步較早.上世紀50年代釋［7］；1998年，Morrison等將早年有關地面瞬變電磁法的研究成果應用到航空電磁存在局部不均勻體的異常響應特征研究［8］；同年，Barongo通過將頻率域的計算結果轉換到時間域研究了航空瞬變電磁法的一維正演問題［9］；2003年，Balch和Rudd等對三角波激勵下AeroTEM系統(tǒng)的on－time數(shù)據(jù)進行反演解釋［10］；2008年，Yin等對半正弦和梯形波激勵均勻半空模型的on－time和off－time電磁響應特征進行模擬和分析［11］.

我國航空電磁勘探開發(fā)研究始于20世紀70年代.1979年航空物探總隊從國外引進雙水獺翼尖硬架三頻電磁系統(tǒng)，目前已退役.后來，國內(nèi)自行研制的Y－11飛機翼尖硬架雙頻航電系統(tǒng).2006年，吉林大學開始自主研制時間域直升機吊艙系統(tǒng)，項目組研制出JLATEM吊艙式時間域直升機電磁勘查系統(tǒng)樣機，目前仍處于試驗階段.國內(nèi)引進的系統(tǒng)目前僅有國土資源部航遙中心從加拿大Aeroquest公司引進的AeroTEM－III型時間域航空電磁系統(tǒng).國內(nèi)很多學者很早就開始研究航空電磁法，1980年樸化榮等發(fā)表了均勻大地上空的時間域電磁響應，并解釋了不同飛行高度觀測的時間域航空電磁資料［12］；隨后，更多學者投入到航空電磁研究領域.2003年羅延鐘等基于電磁勘探理論推導出層狀大地時間域航空電磁法偶極－偶極裝置的正演計算公式和算法［13］；2005年阮百堯等對均勻水平大地上頻率域垂直磁偶源電磁場數(shù)字濾波解法進行研究［14］；2010年吳小平等利用δ單脈沖函數(shù)作為激勵源的時間域電磁響應以實現(xiàn)電磁測深［15］；2011年王世隆等對同心補償式直升機時間域航空電磁法吊艙校準裝置進行了研究［16］.

本文推廣了Yin［11］等人的算法理論，不再僅僅計算均勻半空間電磁響應而是計算任意層狀介質的電磁響應，并將數(shù)值轉換方法推廣到三維電磁響應的正演模擬.我們首先利用快速漢克爾變換數(shù)值濾波算法計算頻率域電磁響應，再通過傅里葉變換將頻率域電磁響應轉換到時間域，得到航空電磁系統(tǒng)的階躍響應.最后，通過和實際發(fā)射波形的褶積獲得對于任意發(fā)射波形的正演模擬算法.我們以目前世界上最先進的航空電磁直升機吊艙系統(tǒng)HeliGEOTEM為例來說明算法的有效性.

2 航空直升機吊艙系統(tǒng)的正演理論

2.1 垂直磁偶極子的電磁響應公式

在目前航空電磁系統(tǒng)中，垂直磁偶極子是最常用的發(fā)射源裝置.在這種裝置下可以測量各磁場分量.這里主要研究同線x方向和垂直z方向上的磁場分量的正演計算.我們假設直角坐標系的原點位于垂直磁偶極正下方的地表，z坐標向下為正，磁偶極子發(fā)射源的高度為h.參考圖1，在測點（x，0，z）處沿著x方向的磁場水平分量為［17］

而磁場垂直分量HZ為

其中z±＝h±z，R＋＝，m＝nIS是發(fā)射偶極矩，而I，n，S分別表示發(fā)射電流、發(fā)射線圈匝數(shù)和面積；T（z-）是漢克爾積分，即

上式的積分稱為漢克爾積分，J0（kr）和J1（kr）分別為零階和一階貝塞爾函數(shù)，核函數(shù)η（k，ξ）的循環(huán)算法參見 Yin和 Hodges［17］.本文應用了 Yin和 Hodges［17］給出的100點漢克爾數(shù)值濾波系數(shù)計算 （3）～（4）式中的漢克爾積分.

2.2 頻率域與時間域的轉換關系

如果頻率域中場值已知，就可通過傅里葉變換將頻率域航空電磁響應轉換到時間域

其中，B（ω）和BI（t）分別為航空電磁系統(tǒng)的頻率域響應和對應的脈沖響應.本節(jié)為簡化起見，我們省略場表達式中的空間坐標.當發(fā)射電流采用負階躍電流i（t）＝u（t）時，系統(tǒng)觀測到的響應稱為階躍響應.由于階躍函數(shù)的傅里葉變換為－1／iω，我們可得到頻率域中諧變電磁場分量B（ω）和時間域中負階躍電流激勵下瞬變電磁響應BS（t）的對應關系［18］為

上式經(jīng)過簡單變換變成

我們利用Yin等［11］給出的160點漢克爾濾波系數(shù)計算（7）式中的半整數(shù)階漢克爾積分.

2.3 時間域航空電磁響應

時間域航空電磁響應傳統(tǒng)的算法基于傅里葉變換，因此有關系式［19］

上式中B0（ω）是系統(tǒng)的單位脈沖響應，而I（ω）為發(fā)射電流的傅氏譜.利用卷積定理以及階躍響應和脈沖響應之間的關系，得到任意發(fā)射電流時間域電磁響應的褶積公式：

其中，I（t）是發(fā)射電流，而BI（t）和BS（t）分別是航空電磁系統(tǒng)的脈沖響應和階躍響應.這里值得注意的是任何系統(tǒng)階躍響應和脈沖響應均存在積分／微分的關系.方程（10）和（11）將 Yin等［17］給出的計算公式做了進一步延伸.

公式（10）和（11）中的褶積可以通過高斯積分來實現(xiàn)［20］.因此，可得到

其中xi和wi分別為高斯抽樣點坐標及對應的加權系數(shù).

3 一維層狀介質脈沖響應和階躍響應特征分析

3.1 模型參數(shù)介紹

假設垂直磁偶極發(fā)射源位于直角坐標系原點的正上方（取Z坐標垂直向下，如圖1），測點坐標為（x，0，z）.h＝30m 為發(fā)射線圈離地面的高度，z＝50m為接收線圈離地面的高度，r＝10m為發(fā)射線圈和接收線圈之間的水平距離.航空電磁系統(tǒng)的發(fā)射偶極矩為615000Am2.圖1給出電磁系統(tǒng)模型示意圖.我們分別討論表1中均勻半空間、兩層和三層介質的情況.

表1 模型參數(shù)Table 1 Model parameters

3.2 脈沖函數(shù)和階躍函數(shù)的電磁響應

圖2給出航空電磁系統(tǒng)的脈沖響應和階躍響應.由圖可以看出，無論是同線分量還是垂直分量，航空電磁系統(tǒng)的脈沖響應在時間接近0時出現(xiàn)奇異性［圖2a，b］，而階躍響應（脈沖響應的積分）則保持平穩(wěn)狀態(tài)并隨時間減小趨近于零［圖2c，d］.因此，利用系統(tǒng)脈沖響應的褶積運算不穩(wěn)定，而利用系統(tǒng)階躍響應的褶積運算是穩(wěn)定的.這更說明，盡管在公式（10）、（11）中的褶積在數(shù)學上是等價的，但由于脈沖響應在t＝0時刻出現(xiàn)奇異性，公式（10）、（11）中只有最后的褶積方法是穩(wěn)定可行的.鑒于此，本文利用（10）和（11）中的最后兩式計算任意發(fā)射波形的航空電磁全時響應.

圖1 模型示意圖Fig.1 Sketch map of AEM system and the layered earth model

4 一維層狀介質任意發(fā)射波形的電磁響應

4.1 發(fā)射波形參數(shù)和精度檢驗

下面以半正弦波為例討論時間域航空電磁系統(tǒng)全時響應的正演模擬問題，參考Fugro的HeliGEOTEM系統(tǒng).如圖3，假設半正弦波的基頻為f＝30Hz，脈沖寬度為4ms，發(fā)射偶極矩為615000Am2.發(fā)射和接收機幾何參數(shù)與圖2所述相同.由于目前尚未見到國內(nèi)外有關時間域航空電磁系統(tǒng)全時響應計算結果，我們針對Yin等［11］給出的半空間模型，利用本文提出的層狀模型計算結果與其進行比較.如圖4所示，吻合程度令人滿意.

圖2 層狀介質的脈沖響應和階躍響應（a）和（b）是x和z方向磁場的脈沖響應，而（c）和（d）是x和z方向磁場的階躍響應.Fig.2 Impulse and step responses of AEM system for a layered medium，where（a）and（b）display the impulse system responses，while（c）and（d）display the step response.

圖3 航空電磁半正弦發(fā)射波形Fig.3 Airborne EM half-sine transmitting waveform

4.2 一維層狀介質的電磁響應

圖5和6分別給出航空電磁系統(tǒng)全時響應B和dB／dt（地電模型參數(shù)見表1）.為比較起見，同時給出利用系統(tǒng)階躍響應和脈沖響應進行褶積的計算結果.

從圖可以看出，無論是同線還是垂直分量，利用系統(tǒng)階躍響應和發(fā)射電流褶積計算的磁場B和磁感應dB／dt之間保持很好的積分／導數(shù)關系，而利用脈沖響應計算的B和dB／dt不存在該積分／導數(shù)關系.另外，由脈沖響應計算的電磁響應與電阻率基本無關，響應曲線基本重合，不能反映地下電性分布特征.這說明算法存在問題.相反，由階躍響應計算的電磁場響應很好地反映地下電性特征的變化.從圖中可以比較明顯地看出不同大地導電率對應不同的晚期衰減速率.大地越良導，電磁信號衰減越慢，反之大地越高阻，電磁信號衰減越快.

比較圖5a和圖5c、圖6a和圖6c可以看出，所有模型中B和dB／dt電磁響應x方向的幅值比z方向的幅值均小一個數(shù)量級，說明z方向的電磁響應具有較高垂向分辨率和勘探深度.進一步比較圖5a和圖5c可以看出，地下電性越良導，供電和斷電（全時域）電磁響應信號B越強，供、斷電期間的峰峰值越大且衰減越慢.因而時間域航空電磁系統(tǒng)可有效應用于地下良導體目標的勘探.比較圖6a和圖6c可以看出，地下電性越良導，電磁響應信號dB／dt越強，供、斷電期間的峰峰值越大且衰減越慢.然而，和磁場響應B相比，磁感應dB／dt曲線形態(tài)更復雜，異常分辨能力明顯減弱.

圖4 本文計算結果和Yin［11］等的計算結果比較Fig.4 Modeling results of this paper compared with those from Yin et al.［11］

圖5 航空電磁層狀介質全時域電磁響應B計算結果其中（a）和（c）是利用系統(tǒng)階躍響應進行褶積，而（b）和（d）是利用系統(tǒng)脈沖響應進行褶積.Fig.5 AEM full-time Bresponses over a layered medium for a half-sine transmitting wave，where （a）and （c）are magnetic field responses calculated by convolving the transmitting current with step response，while （b）and （d）are calculated by convolving with impulse response.

5 三維地質體的電磁響應

下面討論利用積分方程方法計算三維不均勻體的時間域航空電磁響應.本文采用的頻率域正演模擬是基于Raiche［21］提出的算法.它是利用并矢格林函數(shù)理論建立二次感應場和一次場的關系，求解異常體內(nèi)感應電流分布，并通過對異常體內(nèi)感應電流的體積分計算航空電磁系統(tǒng)的頻率響應.通過漢克爾變換，將系統(tǒng)的頻率響應轉換到時間域，并通過和發(fā)射電流波形或其導數(shù)的褶積計算出航空電磁系統(tǒng)的全時域瞬變電磁響應.

5.1 三維模型及參數(shù)介紹

假設垂直磁偶極發(fā)射源位于三維異常體的正上方.航空電磁系統(tǒng)參數(shù)見圖1，地下三維目標體的埋深為d，目標體的圍巖電阻率為100Ωm.我們分別計算了兩種情況：同一深度不同目標體電阻率ρ＝1、10、100Ωm，及同一電阻率不同目標體深度d＝10、30、50m.假設異常體的體積為50m×50m×200m，正演模擬時將其剖分為5m×5m×20m的單元.

5.2 三維地質體的電磁響應

圖8給出三維地質體在同一電阻率ρ＝1Ωm不同深度的電磁響應，而圖9是三維地質體位于同一深度d＝30m，其電阻率不同時的電磁響應.我們利用系統(tǒng)的階躍響應與發(fā)射波形進行褶積計算出x和z方向磁場B和磁感應dB／dt響應.從圖8a、8b和圖9a、9b可以看出早期異常較大，說明供電期間的航空電磁信號易于分辨地下良導體；比較兩圖可以看出z方向幅值比x方向幅值至少高出一個數(shù)量級，這與一維層狀介質的情況類似.比較圖8c、8d可以看出不同埋深良導體無論在磁感應dB／dt還是在磁場響應B上均得到很好地反映.無論對于同線分量還是垂直分量，B和dB／dt之間均有很好的積分／導數(shù)一致性關系.

圖6 航空電磁層狀介質全時域電磁響應dB／dt計算結果（a）和（c）是利用系統(tǒng)階躍響應進行褶積，而（b）和（d）是利用系統(tǒng)脈沖響應進行褶積.Fig.6 AEM full-time dB／dt responses over a layered medium for a half-sine transmitting wave，where（a）and （c）are dB／dt calculated by convolving the transmitting current with step response，while（b）and（d）are calculated by convolving with impulse response.

6 結 論

本文研究了時間域航空電磁系統(tǒng)瞬變?nèi)珪r域電磁響應的正演模擬問題，對于系統(tǒng)供、斷電時產(chǎn)生的B和dB／dt的電磁響應給出了一個穩(wěn)定的算法.較之于傳統(tǒng)的算法，本方法計算更精確、穩(wěn)定，并很好地保持了磁場B和磁感應dB／dt的一致性關系.文中給出的算法不僅適用于一維地質體，也適用于三維復雜地質模型電磁響應的正演模擬.

從理論模擬結果可見，航空瞬變電磁響應的垂直分量比同線分量信號強的多，因而具有較大勘探深度和較高分辨率.三維異常體的供電（on－time）電磁響應強，并且對異常體的導電性非常敏感，因此全時域航空電磁響應觀測和解釋易于識別地下良導體.

必須指出的是：本文提出的算法是利用系統(tǒng)階躍響應和發(fā)射電流的導數(shù)進行褶積計算航空電磁響應，非常適合理論模型的正反演計算；對于實際發(fā)射波形，首先需要對其進行擬合，從而得到其導數(shù)才可計算出航空電磁響應.

圖7 三維地質體模型Fig.7 Three－dimensional geological model

圖8 同一電阻率不同深度三維地質體的航空電磁響應Fig.8 AEM responses for a three－dimensional geological body with the same resistivity but different depths

現(xiàn)在探討一種新的思路：通過利用航空電磁系統(tǒng)接收機在飛行到很高高度時獲得的參考信號來代替發(fā)射電流波形，并和系統(tǒng)階躍響應進行褶積來計算航空電磁系統(tǒng)全時域電磁響應.這將是未來的研究方向.

（References）

［1］強建科，羅延鐘，湯井田等.航空瞬變電磁法的全時域視電阻率計算方法.地球物理學進展，2010，25（5）：1657－1661.Qiang J K，Luo Y Z，Tang J T，et al.The algorithm of alltime calculation of apparent resistivity for Airborne Transient Electromagnetic（ATEM）survey.Progress in Geophysics（in Chinese），2010，25（5）：1657－1661.

圖9 三維地質體在同一深度但具有不同電阻率的航空電磁響應Fig.9 AEM responses for a three－dimensional geological body with the same depth but different resistivities

［2］雷棟，胡祥云，張素芳.航空電磁法的發(fā)展現(xiàn)狀.地質找礦論叢，2006，21（1）：40－53.Lei D，Hu X Y，Zhang S F.Development status quo of airborne electromagnetic.Contributions to Geology and Mineral Resources Research （in Chinese），2006，21（1）：40－53.

［3］牛之璉.時間域電磁法原理.長沙：中南大學出版社，2007.Niu Z L.Principle of Time－domain Electromagnetic Method（in Chinese）.Changsha：Central South University Publishing House，2007.

［4］Smith R.Airborne electromagnetic methods：applications to minerals，water and hydrocarbon exploration.CSEG Distinguished Lecture Tour，2010.

［5］Dobrin M B.Introduction to Geophysical Prospecting.3rd ed.New York：McGraw－Hill，1976.

［6］Nerson P H and Morris D B.Theoretical response of a timedomain airborne，electromagnetic system.Geophysics，1969，34（5），729－738.

［7］Annan A P，Smith R S，Lemieux J，et al.Resistive－limit time－domain AEM apparent conductivity.Geophysics，1996，61（1）：93－99.

［8］Morrison H F，Becker A，Hoversten G M.Physics of airborne EM systems.Exploration Geophysics，1998，29（2）：97－102.

［9］Barongo J O.Selection of an appropriate model for the interpretation of time－domain airborne electromagnetic data for geological mapping.Exploration Geophysics，1998，29（2）：107－110.

［10］Balch S J，Boyko W P，Paterson N R.The AeroTEM airborne electromagnetic system.The leading Edge，2003，22（6）：562－566.

［11］Yin C，Smith R S，Hodges G，Anne P.Modeling results of on－and off－time B and dB／dt for time－domain airborne EM systems.Extended Abstract，70th Annual EAGE Conference and Exhibition，Rome，2008，1－4.

［12］樸化榮，沙樹琴，王延良.均勻大地上空的時間域電磁響應.地球物理學報，1980，23（2）：207－218.Piao H R， Sha S Q， Wang Y L. Time domain electromagnetic response above the surface of a homogeneous earth.Chinese J.Geophys.（in Chinese），1980，23（2）：207－218.

［13］羅延鐘，張勝業(yè)，王衛(wèi)平.時間域航空電磁法－維正演研究.地球物理學報，2003，46（5）：719－724.Luo Y Z，Zhang S Y，Wang W P.A research on one－dimension forward for aerial electromagnetic method in time domain.Chinese J.Geophys.（in Chinese），2003，46（5）：719－724.

［14］阮百堯.均勻水平大地上頻率域垂直磁偶源電磁場數(shù)值濾波解法.桂林工學院學報，2005，25（1）：14－18.Run B Y.Digital filter method of evaluating electromagnetic field from a vertical magnetic dipole above the homogeneous earth.Journal of Guilin University of Technology （in Chinese），2005，25（1）：14－18.

［15］吳小平，柳建新.一種廣義時間域電磁測深法——單脈沖電磁測深方法.地球物理學進展，2010，25（6）：2137－2143.Wu X P，Liu J X.A generalized time－domain electromagnetic sounding——theδfunction electromagnetic sounding.Progress in Geophysics（in Chinese），2010，25（6）：2137－2143.

［16］王世隆，王言章，隨陽軼等.同心補償式直升機時間域航空電磁法吊艙校準裝置研究.地球物理學報，2011，54（9）：2397－2406.Wang S R，Wang Y Z，Sui Y Y，et al.A bird calibration device of helicopter－borne TEM with concentric bucking loop.Chinese J.Geophys.（in Chinese），2011，54（9）：2397－2406.

［17］Yin C，Hodges G.Simulated annealing for airborne EM inversion.Geophysics，2007，72（4）：F189－F196.

［18］Yin C.Electromagnetic induction in a layered conductor with arbitrary anisotropy ［Ph.D.Thesis］Germany：Technical University of Braunschweig，1999.

［19］樸華榮.電磁測深法原理.北京：地質出版社，1990.Piao H R.The theory of Electromagnetic Sounding （in Chinese）.Beijing：Geological Publishing House，1990.

［20］何光渝，周宇斌，張國鳳.FORTRAN 77算法手冊.北京：科學出版社，1993.He G Yu，Zhou Y B，Zhang G F.Numerical Recipes in Fortran 77（in Chinese）.Beijing：Science Press，1993.

［21］Raiche A. 3DEM modeling using integral equation algorithm，AMIRA Project Report P223E，2001.