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

賈亞娟　鄭建波　周紅芳

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

關(guān)鍵詞：采空區(qū) 航空電磁 視電阻率 成像

1 引言

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

2 時間域航空電磁響應(yīng)模擬

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

2.1 時間域航空電磁階躍響應(yīng)快速數(shù)值模擬

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

2.2 基于時間域褶積算法的任意電流發(fā)射波形時間域航空電磁響應(yīng)模擬

時間域航空電磁系統(tǒng)實際發(fā)射波形十分復(fù)雜，而且發(fā)射波形對觀測響應(yīng)影響嚴(yán)重，要獲得準(zhǔn)確的航空電磁成像結(jié)果必須考慮發(fā)射波形的影響，為此本項目采用時間域階躍響應(yīng)與發(fā)射電流時間導(dǎo)數(shù)褶積的方法模擬任意發(fā)射波形時間域航空電磁響應(yīng)。開發(fā)任意發(fā)射波形時間域航空電磁響應(yīng)模擬算法是實現(xiàn)考慮發(fā)射波形的全域視電阻率成像算法的基礎(chǔ)。

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

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

3 技術(shù)方案及算法研究

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

從頻率域麥克斯韋方程組出發(fā)

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

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

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

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

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

3.2 任意電流發(fā)射波形時間域航空電磁響應(yīng)模擬

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

上式中是發(fā)射電流，是時間域階躍響應(yīng)?？梢岳酶咚狗e分方法完成褶積計算。

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

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

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

其中。通過設(shè)定迭代終止條件：

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

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

4 結(jié)論

本文研究的航空電磁全域視電阻率成像算法不受裝置形式、飛行高度的影響，同時考慮了實際發(fā)射電流波形，可以實現(xiàn)對時間域航空電磁數(shù)據(jù)快速、高精度成像。該成果不但可以為城際鐵路建設(shè)過程中的采空區(qū)探測提供更加準(zhǔn)確的數(shù)據(jù)解釋結(jié)果，同樣可以應(yīng)用于礦產(chǎn)勘查、環(huán)境地質(zhì)調(diào)查等諸多領(lǐng)域，因此具有比較廣泛的應(yīng)用和推廣前景。

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

參考文獻(xiàn)：

[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]曾友強(qiáng)，謝國勝，智慶全等. 2019. 小回線源TEM法全域視電阻率成像在某金礦水害中的應(yīng)用[J]. 地球物理學(xué)進(jìn)展，34（2）：0826-0831.

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

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

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

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

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