地空瞬變電磁在積水采空區(qū)探測(cè)中的應(yīng)用

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

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

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

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

1 地空瞬變電磁工作原理

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

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

2 正演理論與模擬驗(yàn)證

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

(1)

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

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

(2)

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

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

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

(3)

(4)

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

(5)

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

(6)

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

圖1 發(fā)射源與接收點(diǎn)相對(duì)位置俯視圖

圖2 瞬變電磁正演模擬與驗(yàn)證對(duì)比

3 約束反演方法

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

(7)

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

(8)

(9)

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

為進(jìn)一步對(duì)比與驗(yà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型模型與計(jì)算結(jié)果對(duì)比

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

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

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

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

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

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

圖4 探測(cè)區(qū)地形

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

圖5 探測(cè)工作布置

5.2 地層與地球物理特征

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

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

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

5.3 探測(cè)效果分析

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

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

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

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

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

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

6 結(jié) 論

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

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