侯彥威，高小偉，李雄偉，郭建磊，姜 濤

(中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710077)

地空瞬變電磁法[1]結(jié)合了地面瞬變電磁法和航空瞬變電磁法各自的特點，相比地面瞬變電磁法能有效提高野外工作效率；地空瞬變電磁法采用地面激發(fā)、無人機(jī)搭載接收探頭空中接收的裝置形式，有效解決了地表情況復(fù)雜地區(qū)的瞬變電磁勘探施工問題，目前已逐步被用于環(huán)境地質(zhì)調(diào)查[2]、煤田采空區(qū)勘查[3，4]等方面。地空瞬變電磁法是Nabighian[5](1988年)基于水平電偶源的工作方式而提出；Smith等[6]早在2001年的時候在加拿大某礦區(qū)對地面、地空和航空瞬變電磁法的探測能力進(jìn)行了對比分析，證明了同等條件下地空瞬變電磁法同樣具有勘探深度大、數(shù)據(jù)信噪比高等優(yōu)點；Mogi等[7]研制了電性源作為發(fā)射源、直升機(jī)搭載磁場接收探頭的數(shù)據(jù)采集方式，即電性源地空瞬變電磁；吉林大學(xué)[8]研制的無人飛艇長導(dǎo)線源時域地-空電磁勘探系統(tǒng)采集的是磁場的時間導(dǎo)數(shù)，促進(jìn)了地空電磁勘探的進(jìn)步。在地空瞬變電磁正演及數(shù)據(jù)處理方面也有大量研究，如李賀和孫懷鳳[9，10]分別實現(xiàn)了三維有限元和三維時域有限差分正演；李肅義等[11]研究了小波分析壓制噪聲；李貅等[12-14]從解釋方面提出了全域視電阻率定義、逆合成孔徑成像、擬地震成像及聯(lián)合解釋等技術(shù)；張澎等[15]將一維自適應(yīng)正則化反演方法應(yīng)用于半航空時間域電磁數(shù)據(jù)反演解釋中，并實現(xiàn)了該反演算法的并行計算；王振榮等[4]利用地空時間域電磁系統(tǒng)在哈拉溝煤礦開展采空區(qū)勘查工作，采用視電阻率和視深度進(jìn)行繪圖。

綜上所述，隨著研究的深入，地空瞬變電磁逐步被推廣應(yīng)用，數(shù)據(jù)反演計算方法也逐步被探索應(yīng)用。本文采用正反演相結(jié)合，通過模擬計算的數(shù)據(jù)，分析認(rèn)識所建模型的地空瞬變電磁響應(yīng)特征，嘗試采用約束反演對數(shù)據(jù)進(jìn)行反演計算，驗證可行，并將其用于實際采空區(qū)探測數(shù)據(jù)處理中，解釋成果符合地質(zhì)電性規(guī)律，低阻異常區(qū)也獲得了驗證。

1 地空瞬變電磁工作原理

地空瞬變電磁工作原理與地面瞬變電磁近似，電性源地空瞬變電磁是在地面布置一定長度的導(dǎo)線源，采用大功率設(shè)備沿導(dǎo)線向地下發(fā)射階躍變化的電磁波信號，此信號即為激發(fā)場源(稱為“一次場”)，根據(jù)法拉第電磁感應(yīng)定律，當(dāng)發(fā)射電流被關(guān)斷時，大地或探測目標(biāo)體在一次場的作用下，其內(nèi)部會感生具有空間和時間特性的渦流，其大小與諸多因素有關(guān)，如目標(biāo)體的空間特征和電性特征、激發(fā)場的特征等，而且因為熱損耗的緣故會逐漸減弱直至消失。通過搭載在無人機(jī)上的接收設(shè)備接收這種渦流產(chǎn)生的電磁場(稱為“二次場”)，二次場具有強(qiáng)弱之分，也具有空間和時間特性，研究這些特性可分析地層或地下目標(biāo)體的電性，并據(jù)此完成探測。

地空瞬變電磁探測裝置的飛行高度可靈活調(diào)整，一次布設(shè)可在短時間內(nèi)完成大面積勘探。

2 正演理論與模擬驗證

在均勻、各向同性、有耗、非磁性、無源媒質(zhì)中，Maxwell方程組為：

(1)

式中，E表示電場強(qiáng)度；B表示磁感應(yīng)強(qiáng)度；H表示磁場強(qiáng)度；σ表示介質(zhì)電導(dǎo)率；ε表示介電常數(shù)；t表示時間。

瞬變電磁勘探計算中一般忽略位移電流，而為了獲得顯式的時間迭代格式，需在上方程組中的式(1b)中加入虛擬介電常數(shù)項，在有源區(qū)域加入電流密度項，式(1b)即為：

(2)

式中，γ表示具有介電常數(shù)的量綱；Js表示源電流密度。

地空瞬變電磁法的激發(fā)源一般采用電性源[1]，為便于采用時間域有限差分法進(jìn)行正演計算，常在位于直角坐標(biāo)系中晶胞網(wǎng)格的Y軸方向加載電性源。

在地空瞬變電磁正演過程中仍需考慮其在低頻情況下正演結(jié)果的正確性，因此需對磁場各分量采用低頻近似處理，將方程組(1)中的式(a)、(d)和式(2)在直角坐標(biāo)系下展開，得分項式為：

(3)

(4)

對式(3)和式(4)用差分代替微分，由于Euler前向差分對離散時間步的要求比較嚴(yán)格，故空間離散采用后向差分，均勻網(wǎng)格剖分中電/磁場的時間采樣恰好在兩相鄰磁/電場采樣時刻的中心，故時間離散采用中心差分，可以得到各個分量的迭代公式。其中，有源區(qū)域中的Ez分量迭代公式為：

(5)

采用第一類邊界條件，為保證電磁場傳播的因果關(guān)系，將需要計算的整個模型剖分得很大且在計算過程中時間域和空間域的穩(wěn)定性條件(式6)滿足時，在邊界處可將電場與磁場進(jìn)行強(qiáng)制性賦零。

(6)

在均勻半空間條件下構(gòu)建電性發(fā)射源與接收試驗?zāi)Ｐ?，基于三維時域有限差分算法進(jìn)行地空瞬變電磁正演計算，將計算結(jié)果與解析解進(jìn)行對比以驗證程序的可靠性。發(fā)射源與接收點相對位置如圖1所示，模型參數(shù)為：均勻半空間地層電阻率100Ω·m，電性發(fā)射源長度500m，接收點坐標(biāo)為(x=300m，y=-50m)，接收點高度50m，電流方向沿Y軸正向，電流1A，二次場采樣時間10ms，最小網(wǎng)格尺寸為10m且總剖分網(wǎng)格數(shù)為221×221×200。將三維時域有限差分模擬計算結(jié)果和一維解析解通過直接與誤差統(tǒng)計的方式進(jìn)行對比，繪制對比圖(圖2)，其中圖2(a)為直接對比圖，響應(yīng)曲線形態(tài)、特征基本一致；圖2(b)為誤差對比曲線圖，經(jīng)統(tǒng)計二者的最大誤差為2.176%，誤差小于1%的占比89.47%，可見三維時域有限差分模擬計算的結(jié)果滿足精度要求；同時由圖2(a)可見，地空瞬變電磁的感應(yīng)電動勢衰減曲線與地面瞬變電磁衰減曲線規(guī)律近似，因此可在地面瞬變電磁數(shù)據(jù)反演方法的基礎(chǔ)上對地空瞬變電磁數(shù)據(jù)進(jìn)行反演試算。

圖1 發(fā)射源與接收點相對位置俯視圖

圖2 瞬變電磁正演模擬與驗證對比

3 約束反演方法

地空瞬變電磁數(shù)據(jù)反演采用自適應(yīng)的約束反演方法，此方法在反演的過程中采用自主調(diào)整正則化因子的方案，反演矩陣條件數(shù)比較大時增大正則化因子，當(dāng)次反演擬合差大于上次反演的擬合差時減小正則化因子；為減小反演過程的病態(tài)性，反演的層厚按照等對數(shù)間隔進(jìn)行離散；同時對每一層電阻率進(jìn)行上下限約束，如果超過給定的上下限，那么就減小反演模型的步長，直到滿足閾值[16]。反演目標(biāo)函數(shù)及正則化因子為：

(7)

[(WJ1)TWJ1+μATA]Δmk=(WJ1)TWΔdk

(8)

(9)

4 地質(zhì)模型地空瞬變電磁響應(yīng)特征

為進(jìn)一步對比與驗證上述約束反演方法的適用性，特建立地質(zhì)模型分別如下：

1)均勻半空間模型：電阻率為100Ω·m。

2)與陜北侏羅系煤田縱向整體電性特征相似的K型地質(zhì)模型(圖3a)：第一層模擬第四系地層，電阻率為100Ω·m，厚度50m；第二層假設(shè)為侏羅系含煤地層，電阻率為1000Ω·m，厚度為50m；第三層電阻率為100Ω·m。

3)含低阻異常體的K型模型：在侏羅系含煤地層中部加入電阻率為10Ω·m的低阻異常體(假設(shè)為積水采空區(qū)，見圖3a)，低阻異常體體積為100m×100m×40m。

圖3 地空瞬變電磁K型模型與計算結(jié)果對比

模擬的電性發(fā)射源長度500m，接收點高度50m，接收點位于低阻異常體中心上方，電流1A，二次場采樣時間10ms，最小網(wǎng)格尺寸為10m且總剖分網(wǎng)格數(shù)為221×221×200。

對上述模型進(jìn)行三維時域有限差分模擬計算，計算結(jié)果如圖3(b)所示，圖中三條衰減曲線分別為均勻半空間模型感應(yīng)電動勢衰減曲線、K型地質(zhì)模型感應(yīng)電動勢衰減曲線和K型地質(zhì)模型中間層加入低阻異常體后的感應(yīng)電動勢衰減曲線。可見，三條衰減曲線隨時間均呈逐漸衰減趨勢，在早期略有交叉；中期K型地質(zhì)模型感應(yīng)電動勢值較低，加入低阻異常體后的感應(yīng)電動勢值有所增加但仍小于均勻半空間模型感應(yīng)電動勢值；晚期有重合趨勢，但數(shù)值仍有差異，K型模型模擬計算值始終小于均勻半空間模型值。以上衰減曲線的不同特征反映了不同地質(zhì)模型的電性分布規(guī)律，為反演取得好的結(jié)果奠定了基礎(chǔ)。

將圖3(b)中K型地質(zhì)模型及其加入低阻異常體后模型的感應(yīng)電動勢數(shù)據(jù)分別進(jìn)行約束反演計算，其結(jié)果如圖3(c)所示。可見，從上而下兩種模型的反演曲線在淺部均接近模型第一層電阻率值，隨深度增加，電阻率值也在逐漸增大，且逐步跨進(jìn)第二層高阻層，直至達(dá)到極值后又逐漸降低，約在-170m附近接近第三層電阻率值，整體由淺至深均呈“低阻-高阻-低阻”的變化趨勢，與K型地質(zhì)模型三層電性規(guī)律吻合較好。加入低阻異常體后的反演曲線所示電阻率值更低，可見該探測方法及約束反演對低阻異常體反應(yīng)明顯，以此通過橫向?qū)Ρ瓤蓪崿F(xiàn)積水采空區(qū)的探測。

5 積水采空區(qū)探測

5.1 測區(qū)地形與地空瞬變電磁工作布置

某煤礦地處陜北黃土高原，地形起伏較大、地表溝壑縱橫。探測區(qū)域內(nèi)海拔標(biāo)高最高點達(dá)1320m，最低點約為1160m，高差達(dá)160m(如圖4所示)，地面施工難度較大，所以采用地空瞬變電磁進(jìn)行勘探。

圖4 探測區(qū)地形

探測區(qū)域與工作布置如圖5所示，煤礦生產(chǎn)巷道由西南向東北進(jìn)入該探測區(qū)域(面積約0.7km2)，探測區(qū)域前方有小煤窯采空區(qū)分布，積水情況不明，采用地空瞬變電磁探測，設(shè)計線距40m，點距采用10m，沿設(shè)計線路迂回折返飛行，飛行高度為探測區(qū)域內(nèi)最高點以上50m；電性發(fā)射源位于測區(qū)東南向3km處，長1.5km。

圖5 探測工作布置

5.2 地層與地球物理特征

煤礦地層、相應(yīng)巖性與地球物理特征見表1。

表1 地層巖性及其地球物理特征

由表1可見，探測區(qū)域由淺至深地層電性整體呈“低阻—高阻—低阻”的變化特征。當(dāng)煤層被開采采空區(qū)未積水呈更高阻的特征；當(dāng)采空區(qū)有積水時，水體的浸潤使得周圍介質(zhì)電阻率明顯降低，其電阻率值將明顯降低，一般會小于30Ω·m，使得積水采空區(qū)與周圍介質(zhì)產(chǎn)生明顯的電性差異，上述導(dǎo)電性差異為地空瞬變電磁勘探提供了條件。

5.3 探測效果分析

按照本文前述反演方法對實測數(shù)據(jù)進(jìn)行計算，其結(jié)果可反映原生地層和被采動地層兩種情況，故本文舉例兩條典型斷面如圖6所示，圖6(a)為反映原生地層電性的3線反演電阻率斷面，圖6(b)為反映原生地層被采空區(qū)破壞且存在低阻異常區(qū)的7線反演電阻率斷面。斷面圖橫向為橫向水平距離，縱向為高程，高程1190m以淺為低阻層，3-3煤層附近為高阻層，深部1110m以下為低阻層，詳見斷面圖右側(cè)柱狀電阻率標(biāo)尺。

圖6 地空瞬變電磁探測反演電阻率斷面圖

由圖6可見，反演電阻率橫向上整體呈似層狀分布，淺部為低阻層，對應(yīng)新近系以淺地層；中部為高阻層，對應(yīng)侏羅系含煤地層，探測目標(biāo)層3-3煤層恰位于本層中部；深部為低阻層，對應(yīng)三疊系地層。此種由淺至深的電性變化規(guī)律與前述地球物理特征一致，表明反演結(jié)果反映地層實際。由圖6(a)可見含煤地層橫向上連續(xù)性好，電性分布均勻，表明地層為未受采動破壞的原生地層；圖6(b)中部高阻層橫向上的連續(xù)性在420～800m處被打破，電阻率值突然降低，推斷為低阻異常區(qū)，同時也表明此處地層應(yīng)受到了采動破壞且局部積水。

依據(jù)3-3煤層底板高程和全區(qū)反演電阻率數(shù)據(jù)，抽取本層電阻率數(shù)據(jù)并橫向取異常，如圖7所示，圖7中有顏色填充處為低阻區(qū)域。

圖7 地空瞬變電磁探測3-3煤層順層平面圖

由圖7可見，探測區(qū)內(nèi)大部分區(qū)域反演電阻率等值線較為稀疏，表明該目標(biāo)地層大部分區(qū)域電性值較為穩(wěn)定，可以此間接推測等值線稀疏區(qū)域目標(biāo)煤層未受采動破壞；同時也發(fā)現(xiàn)一大一小兩處低阻異常區(qū)，經(jīng)與后期提供的小煤窯老巷道對比，此兩處低阻異常區(qū)均與老巷道有對應(yīng)關(guān)系，推測其應(yīng)為小煤窯沿巷道向周圍開采煤層形成的采空區(qū)積水所致。從圖7中煤層底板等高線分布趨勢可見，大范圍的低阻異常區(qū)恰位于局部低洼區(qū)域，在有采空空間存在的情況下，利于水的積聚，所以此大范圍低阻異常區(qū)推測為采空積水區(qū)的依據(jù)充分，且經(jīng)鉆探(TC1號鉆孔)驗證，確為積水采空區(qū)。

6 結(jié) 論

1)嘗試使用地面瞬變電磁數(shù)據(jù)反演方法對地空瞬變電磁數(shù)據(jù)進(jìn)行反演計算，探測結(jié)果證明該約束反演算法可以有效的實現(xiàn)地空瞬變電磁數(shù)據(jù)的處理工作，為地空瞬變電磁數(shù)據(jù)處理提供一條新的思路。

2)利用地空瞬變電磁法實現(xiàn)煤礦積水采空區(qū)探測，將探測結(jié)果與采掘巷道資料及鉆孔資料對比，驗證了地空瞬變電磁法對于煤礦積水采空區(qū)探測的可行性與高效性，可為后期煤礦開采提供有效的參考資料。