軌道交通采空區探測中的航空電磁全域視電阻率成像算法研究

賈亞娟　鄭建波　周紅芳

摘 要：目前采空區探測的方法主要有地震勘探、地面瞬變電磁法，高密度電法等，然而在地形地貌復雜的地區開展傳統地面勘探時，工作難度大、效率低、成本高，與上述方法相比，時間域航空電磁法可以快速地、大面積地進行采空區的勘查，工作效率大大提高，特別是在地形條件復雜的山區等，更能凸顯其優勢。

關鍵詞：采空區 航空電磁 視電阻率 成像

1 引言

近年來，我國的軌道交通有了日新月異的發展，在幾條貫穿全國的主干線基礎上已經建成了四通八達的交通網。為了使人們的出行更加快速、方便，開展城際高速鐵路和城際軌道交通建設成為當前的工作重點。然而部分城市在過去很長一段時間對地下煤炭和礦產資源的開采留下了大量采空區，這對城際軌道交通建設埋下了巨大安全隱患。對采空區進行精細探測成為保障軌道交通安全的關鍵。然而全域視電阻率成像方法并未應用于時間域航空電磁數據成像當中，與地面和井中電磁法相比，航空電磁法的發射波形十分復雜，不再是簡單的階躍電流波形，而當前的全域視電阻率成像方法并未考慮實際發射電流波形的影響，這將導致時間域航空電磁數據成像精度受到很大影響。為此，本文開展任意發射波形時間域航空電磁數據全域視電阻率成像算法研究，對提高時間域航空的電磁數據解釋水平具有積極的促進作用。

2 時間域航空電磁響應模擬

針對當前軌道交通采空區勘探過程中航空電磁觀測數據成像技術的不足，開展任意發射波形時間域航空電磁快速成像算法研究，利用時間域褶積算法模擬任意電流發射波形時間域航空電磁響應，基于反函數定理開發考慮發射波形的全域視電阻率成像方法和數據處理軟件，最終實現時間域航空電磁數據快速、高精度成像，進而推動軌道交通采空區航空電磁勘探技術的發展。

2.1 時間域航空電磁階躍響應快速數值模擬

從頻率域麥克斯韋方程組出發推導一維介質模型頻率域航空電磁響應表達形式，利用漢克爾變換算法獲得頻率航空電磁響應。采用基于正弦變換的頻時轉換技術將頻率域航空電磁響應轉換到時間域，完成階躍發射波形時間域航空電磁響應快速數值模擬。在階躍響應模擬過程中考慮航空電磁觀測系統的裝置形式和飛行高度的影響。

2.2 基于時間域褶積算法的任意電流發射波形時間域航空電磁響應模擬

時間域航空電磁系統實際發射波形十分復雜，而且發射波形對觀測響應影響嚴重，要獲得準確的航空電磁成像結果必須考慮發射波形的影響，為此本項目采用時間域階躍響應與發射電流時間導數褶積的方法模擬任意發射波形時間域航空電磁響應。開發任意發射波形時間域航空電磁響應模擬算法是實現考慮發射波形的全域視電阻率成像算法的基礎。

2.3 基于反函數定理的全域視電阻率成像方法

通過分析時間域航空電磁響應隨地下介質電性的變化規律，可知地下介質電阻率與時間域航空電磁觀測響應間存在一一對應關系。基于這一單調特性，利用反函數定理和泰勒展開公式將積分形式轉化為級數形式，進而建立一種基于迭代格式的全域視電阻率成像方法。由于在航空電磁響應數值模擬過程中考慮了航空電磁系統的裝置形式、飛行高度和發射電流波形，因此利用其定義的全域視電阻率不再受到裝置形式、飛行高度和發射電流波形的影響，從而可以實現快速、高精度成像技術。開發考慮電流發射波形的全域視電阻率定義方法是對時間域航空電磁數據進行高精度成像的關鍵。

3 技術方案及算法研究

3.1 階躍波形時間域航空電磁數值模擬方法

從頻率域麥克斯韋方程組出發

其中E和H分別為電場和磁場，ω是角頻率，σ、μ、ε分別是地下介質電導率、磁導率和介電常數。通過引入謝坤諾夫勢和格林函數可以獲得一維頻率域航空電磁響應：

其中m為發射源磁矩，h為系統飛行高度，z是接收點坐標，rTE是地層對電磁波的反射系數，與地層厚度和電阻率有關，u0是關于頻率的函數，J0是0階貝塞爾函數。利用漢克爾變換技術可以將（3）式寫成數值累加形式，通過合理設置漢克爾變換系數即可完成頻率域航空電磁響應數值模擬。

利用頻率域響應實部和虛部的奇偶性，由傅里葉逆變換可以將頻率域航空電磁響應轉換到間域：

其中Im[ ]表示取虛部。利用正弦函數和貝塞爾函數之間的關系可以得到時間域航空電磁階躍響應：

J1/2為0.5階貝塞爾函數。在次利用漢克爾變換即可完成時間域航空電磁階躍響應計算。

3.2 任意電流發射波形時間域航空電磁響應模擬

對于任意波形發射電流的時間域電磁響應一般采用輸入電流與系統響應進行褶積求得。然而Smith和Lee（2004）通過對航空電磁法的研究發現脈沖響應是不穩定的，當時間t趨近于0時，半空間的脈沖響應趨近于t-1/2，所以脈沖響應在t=0處存在一個奇點，這對于瞬變電磁響應B（t）及其對時間的導數dB/dt在斷電瞬間的穩定性有很大影響。為此本項目采用階躍響應與電流波形時間導數褶積的方法實現任意電流發射波形時間域航空電磁響應模擬。利用脈沖響應與階躍響應間的微積分關系可得任意電流發射波形時間域航空電磁響應：

上式中是發射電流，是時間域階躍響應。可以利用高斯積分方法完成褶積計算。

3.3 基于反函數定理的全域視電阻率成像方法

在實際工作過程中，探測結束后觀測裝置和飛行高度已經成為已知量，此時航空電磁響應可以簡寫為。通過分析時間域航空電磁響應隨地下介質電性的變化關系可知，由于響應是關于電阻率ρ的單調函數，根據反函數定理知，必然存在一個電阻率值ρ唯一地對應著一個值，為此，采用泰勒展開的辦法，將磁感應強度的積分表達式展成級數形式，取其線性主部，建立起迭代關系的視電阻率定義式。在初始電阻率的鄰域內對進行泰勒展開：

取其線性主部建立迭代關系的視電阻率定義式完成視電阻率計算：

其中。通過設定迭代終止條件：

其中是實測航空電磁數據，是電阻率為的半空間模型航空電磁理論響應，ε是迭代終止誤差。

基于反函數定理和泰勒展開公式建立一種基于迭代格式的全域視電阻率成像方法，在此過程中考慮實際航空電磁觀測系統工作過程中的裝置形式和飛行高度，可以實現時間域航空電磁數據快速成像。在時間域航空電磁全域視電阻率成像過程中，考慮實際電流發射波形的影響，通過時間域褶積算法模擬實際電流發射波形時間域航空電磁響應，實現任意發射波形航空電磁響應視電阻率成像，有效提高了航空電磁成像精度。

4 結論

本文研究的航空電磁全域視電阻率成像算法不受裝置形式、飛行高度的影響，同時考慮了實際發射電流波形，可以實現對時間域航空電磁數據快速、高精度成像。該成果不但可以為城際鐵路建設過程中的采空區探測提供更加準確的數據解釋結果，同樣可以應用于礦產勘查、環境地質調查等諸多領域，因此具有比較廣泛的應用和推廣前景。

陜西省教育廳科研計劃項目資助（項目編號：21JK0746）。

參考文獻：

[1]Huang H， Rudd J. 2008. Conductivity-depth imaging of helicopter-borne TEM data based on a pseudolayer half-space model[J].Geophysics，73（3）：F115-F120.

[2]曾友強，謝國勝，智慶全等. 2019. 小回線源TEM法全域視電阻率成像在某金礦水害中的應用[J]. 地球物理學進展，34（2）：0826-0831.

[3]雷棟，胡祥云，張素芳.2006.航空電磁法的發展現狀[J].地質找礦論叢，21（1）： 40-44.

[4]武軍杰，李貅，智慶全等.2017. 電性源地-井瞬變電磁全域視電阻率定義[J].. 地球物理學報， 60（4）：1595-1605.

[5]殷長春，張博，劉云鶴，等. 2015. 航空電磁勘查技術發展現狀及展望[J]. 地球物理學報，58（8）：2637-2653.

[6]張瑩瑩，李貅，姚偉華，智慶全，李佳.多輻射場源地空瞬變電磁法多分量全域視電阻率定義[J].地球物理學報，2015，58（8）：2745-2758.

[7]趙越，李貅，王祎鵬. 2015. 大回線源瞬變電磁全域視電阻率定義[J].地球物理學進展，30（4）：1856-1863.