礦井巷道-鉆孔瞬變電磁二維擬地震反演方法及應(yīng)用

范 濤

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

積水采空區(qū)對(duì)于煤礦井下采掘生產(chǎn)危害極大，必須進(jìn)行治理。山西由于歷史原因存在大量使用巷道穿采技術(shù)的小煤窯，遺留下大量巷道尺度的采空區(qū)，這些采空區(qū)往往與含水層連通，極易引發(fā)掘進(jìn)工作面水害事故[1-3]。在我國(guó)煤礦井下，采空區(qū)井下超前探測(cè)目前以鉆探為主，物探為輔。而采空巷道規(guī)模遠(yuǎn)小于一般采空區(qū)，正常布設(shè)的超前探放水鉆孔密度很難對(duì)其進(jìn)行覆蓋，常規(guī)的物探超前探測(cè)技術(shù)也難以對(duì)這一尺度的異常體進(jìn)行準(zhǔn)確定位，更勿論對(duì)其形態(tài)、規(guī)模等參數(shù)做出精確解釋[4-7]?？紤]到鉆孔孔內(nèi)空間距離地質(zhì)異常體很近，我可以利用它遠(yuǎn)離巷道、靠近地質(zhì)異常這一特性，在巷道孔口處發(fā)射瞬變電磁信號(hào)，在鉆孔孔內(nèi)逐點(diǎn)接收，探測(cè)鉆孔徑向可能存在的小規(guī)模地質(zhì)異常，并對(duì)其做出較為準(zhǔn)確的地質(zhì)解釋?zhuān)Ｕ厦旱V的安全生產(chǎn)。

利用鉆孔進(jìn)行瞬變電磁工作，以地-井裝置研究最多，該裝置屬于二維工作裝置，隨收發(fā)距改變數(shù)據(jù)特征差異很大，對(duì)資料的處理解釋帶來(lái)了很大困難，所以熱點(diǎn)研究方向主要集中在數(shù)據(jù)特征的總結(jié):① 孔中信號(hào)的響應(yīng)特征。通過(guò)物理模擬和多種數(shù)值模擬手段總結(jié)全空間條件下簡(jiǎn)單異常體(球體或板體)不同參數(shù)情況下的特征關(guān)系曲線(xiàn)及符號(hào)變化規(guī)律[8-10]。② 導(dǎo)電介質(zhì)對(duì)信號(hào)的影響規(guī)律。通過(guò)理論推導(dǎo)和數(shù)值模擬研究圍巖背景場(chǎng)或覆蓋層對(duì)孔中瞬變電磁信號(hào)的影響，以及導(dǎo)電覆蓋層下多個(gè)目標(biāo)體的信號(hào)響應(yīng)[11-14]。與正演工作相比，地-井瞬變電磁處理解釋方面的研究工作開(kāi)展很少，很難滿(mǎn)足實(shí)際需求，具體可分為兩類(lèi):① 從觀測(cè)數(shù)據(jù)直接解釋。根據(jù)正演研究得出的數(shù)據(jù)特征和規(guī)律，提出一些異常體定性或半定量解釋方法，如通過(guò)多測(cè)道曲線(xiàn)變化判斷異常體的大致方位和利用三分量信號(hào)矢量交匯定位異常體中心位置[15-16]。② 近似反演方法。通過(guò)一定的前提假設(shè)對(duì)地質(zhì)情況進(jìn)行簡(jiǎn)化，提出一些反演方法，如需要假定介質(zhì)關(guān)于鉆孔柱對(duì)稱(chēng)且源為磁偶極源的局部非線(xiàn)性近似快速反演和基于導(dǎo)電薄板等效渦流的異常反演方法[17-18]。在隧道或巷道中開(kāi)展孔中瞬變電磁工作的相關(guān)資料更少，最早是國(guó)外VELLA L(1997)將地-井裝置移至金屬礦巷道中探測(cè)黃鐵礦[19]，之后國(guó)內(nèi)王世睿(2016)研究了隧道10 m深鉆孔中的瞬變電磁信號(hào)特性，提出根據(jù)發(fā)射線(xiàn)框改變產(chǎn)生的異常變化可定性解釋隧道掘進(jìn)工作面前方異常體方位[20]。孫懷鳳等通過(guò)物理模擬試驗(yàn)證明了孔中瞬變電磁信號(hào)可用于判斷隧道掘進(jìn)工作面前方是否存在異常構(gòu)造[21]。筆者(2017)研究了煤礦井下孔中動(dòng)源定接收裝置的超前探測(cè)技術(shù)，利用鉆孔提高了常規(guī)井下瞬變電磁超前探測(cè)精度[22]。陳丁(2018)通過(guò)在全空間一維背景上增加三維異常體的積分方程數(shù)值模擬研究了煤礦巷道垂直孔中瞬變電磁特性[23]。

總體來(lái)看，盡管理論上巷道-鉆孔/地-井瞬變電磁接收裝置更加靠近異常，能通過(guò)地層壓制電磁干擾，應(yīng)該比礦井、地面方法有更優(yōu)的異常反映和解釋精度，但由于在于定源動(dòng)接收情況下，孔中瞬變電磁其實(shí)是一組二維裝置，尚缺乏好的反演成像技術(shù)，實(shí)際資料處理精度反而遠(yuǎn)落后于常規(guī)方法，僅能作為補(bǔ)充手段。所以必須另辟蹊徑，避開(kāi)收發(fā)距帶給擴(kuò)散場(chǎng)的復(fù)雜影響，采用數(shù)學(xué)辦法將其轉(zhuǎn)換為虛擬波場(chǎng)，針對(duì)特征簡(jiǎn)單的波場(chǎng)記錄進(jìn)行高效精準(zhǔn)的擬地震反演。

瞬變電磁擬地震處理的基礎(chǔ)是LEE等(1989)建立的電磁擴(kuò)散方程與虛擬波動(dòng)方程之間的數(shù)學(xué)關(guān)系[24]，研究熱點(diǎn)主要集中在:① 波場(chǎng)反變換。目前實(shí)現(xiàn)方法有全時(shí)段預(yù)條件正則化算法和分段奇異值分解法[25-28]。正則化法能基本還原波場(chǎng)特征，抗干擾能力強(qiáng)，但展寬效應(yīng)和受正則化約束類(lèi)型影響均較大;奇異值分解結(jié)果展寬效應(yīng)小，但抗噪能力弱，結(jié)果穩(wěn)定性差。② 波場(chǎng)分辨率增強(qiáng)。該研究主要為提高波場(chǎng)反變換效果，包括不同層位波幅的能量均衡、消除展寬效應(yīng)的反褶積和壓制噪音的相干疊加等方法[29-30]。此舉雖能提高分辨率，但也必然損失部分有用信息。③ 虛擬波場(chǎng)速度分析、反演及成像。筆者之前通過(guò)等效導(dǎo)電平面法和聲波全波形反演方法實(shí)現(xiàn)了中心回線(xiàn)裝置數(shù)據(jù)的初始速度建模、速度分析和擬電阻率反演[31-32]。成像方面國(guó)外用電磁波旅行時(shí)實(shí)現(xiàn)了虛擬波場(chǎng)層析成像[33]，國(guó)內(nèi)則是對(duì)中心回線(xiàn)裝置數(shù)據(jù)進(jìn)行了Kirchhoff積分偏移成像[34-35]。

筆者將對(duì)井下巷道-鉆孔定源動(dòng)接收裝置數(shù)據(jù)進(jìn)行特征分析，通過(guò)基于相關(guān)疊加的滑動(dòng)時(shí)窗波場(chǎng)反變換算法將二維瞬變電磁數(shù)據(jù)轉(zhuǎn)換為虛擬波場(chǎng)數(shù)據(jù)，并研究其動(dòng)校正技術(shù)，最終利用筆者以前完成的全波形反演方法實(shí)現(xiàn)巷道-鉆孔數(shù)據(jù)的擬地震反演成像，精細(xì)解釋鉆孔徑向存在的小規(guī)模地質(zhì)異常，為井下準(zhǔn)確探測(cè)積水采空巷道提供技術(shù)支撐。

1 方法原理

1.1 工作裝置

筆者提出的工作裝置如圖1所示，在孔口巷道處布設(shè)一多匝小線(xiàn)框發(fā)射源，線(xiàn)框邊長(zhǎng)不大于2 m，接收使用磁探頭。將發(fā)射線(xiàn)圈固定放置在孔口，法線(xiàn)指向鉆孔延伸方向，之后沿鉆孔按一定點(diǎn)距(如2 m)向孔中推送接收探頭并逐點(diǎn)進(jìn)行二次場(chǎng)測(cè)量，直至孔底。對(duì)于單個(gè)鉆孔可觀測(cè)三分量數(shù)據(jù)，通過(guò)垂直(Z)分量數(shù)據(jù)進(jìn)行反演成像，水平(X,Y)分量對(duì)異常體中心進(jìn)行定位;對(duì)于多個(gè)鉆孔，可以通過(guò)對(duì)不同鉆孔的垂直分量反演成像結(jié)果進(jìn)行綜合分析，聯(lián)合解釋地質(zhì)異常體位置。

圖1 施工裝置示意Fig.1 Schematic diagram of construction device

1.2 孔中數(shù)據(jù)特征

設(shè)置如圖2所示模型，可研究巷道-鉆孔瞬變電磁裝置的數(shù)據(jù)特征，以指導(dǎo)資料處理方法研究。模型將一個(gè)2 m×2 m的發(fā)射線(xiàn)圈布置在掘進(jìn)工作面孔口處，發(fā)射電流為1 A，掘進(jìn)巷道規(guī)格為20 m×4 m×4 m，鉆孔長(zhǎng)度為80 m，在鉆孔前方60 m處偏上5 m位置放置一個(gè)低阻異常體(積水巷道)，鉆孔從異常體下方經(jīng)過(guò)，異常體規(guī)格為4 m×20 m×4 m，分別在鉆孔孔深0,20,40,60,80 m處設(shè)置5個(gè)接收點(diǎn)。設(shè)置巷道電阻率為10 000 Ω·m，煤層電阻率為1 000 Ω·m，異常體電阻率為10 Ω·m。

圖2 模型示意Fig.2 Schematic diagram of the model

圖3 Z分量數(shù)據(jù)對(duì)比Fig.3 Comparison of Z component data

對(duì)圖2模型用時(shí)域有限差分三維正演算法模擬了含異常體和不含異常體兩種情況，剖分網(wǎng)格邊長(zhǎng)為0.5 m。5個(gè)接收點(diǎn)不同情況下得到的Z分量數(shù)據(jù)對(duì)比情況如圖3所示。從圖3(a)可明顯看出，隨著入孔深度的增加，感應(yīng)電動(dòng)勢(shì)曲線(xiàn)由一條單調(diào)曲線(xiàn)逐漸變?yōu)橐粭l單峰曲線(xiàn)，孔深越深，對(duì)應(yīng)測(cè)點(diǎn)的峰值越低，峰值對(duì)應(yīng)的時(shí)間越晚。從圖3(b)可看出，由于異常體體積較小，僅距離異常體最近的60 m孔深處數(shù)據(jù)曲線(xiàn)上對(duì)其有明顯的拱起反映，其余孔深數(shù)據(jù)曲線(xiàn)上肉眼無(wú)法識(shí)別出異常響應(yīng)特征，說(shuō)明距離異常體越近，反映越強(qiáng)，這意味著采用孔內(nèi)測(cè)量的方式確實(shí)有利于提高異常分辨能力。

但因測(cè)量曲線(xiàn)已不再符合單調(diào)衰減特征，所以很難對(duì)該裝置的視電阻率給出類(lèi)似中心回線(xiàn)裝置晚期定義形式的簡(jiǎn)單公式，其精細(xì)的分辨能力就無(wú)法通過(guò)成像進(jìn)行高精度呈現(xiàn)，必須考慮新的反演成像辦法。

1.3 波場(chǎng)反變換

擴(kuò)散方程到波動(dòng)方程轉(zhuǎn)換[27]的公式為

(1)

從式(1)可以看出，擴(kuò)散場(chǎng)f(x,y,z,t)和波動(dòng)場(chǎng)u(x,y,z,τ)中包含x,y,z三個(gè)空間坐標(biāo)元素，說(shuō)明該公式對(duì)空間任意位置的兩種場(chǎng)量變換均適用，可將該公式用于非中心回線(xiàn)裝置的其他瞬變電磁工作裝置。

式(1)離散表達(dá)式為

(2)

式中，f(x,y,z,ti)為瞬變電磁場(chǎng)量;u(x,y,z,τj)為f(x,y,z,ti)所對(duì)應(yīng)的虛擬波場(chǎng)量;τ為與瞬變電磁場(chǎng)的時(shí)間t相對(duì)應(yīng)的虛擬波場(chǎng)時(shí)間;hj為積分系數(shù);i為瞬變電磁場(chǎng)數(shù)據(jù)時(shí)間測(cè)道數(shù);j為虛擬波場(chǎng)數(shù)據(jù)時(shí)間測(cè)道數(shù)。

(3)

是式(2)的核函數(shù)。

式(1)的穩(wěn)定性很差，數(shù)值化后的方程(2)屬于極度病態(tài)的類(lèi)型，其病態(tài)程度與方程組階數(shù)是高度相關(guān)的，因此，必須大幅降低波場(chǎng)反變換式的線(xiàn)性方程組階數(shù)，有效改善方程的不適定性。

最初提出的處理方法被稱(chēng)為分時(shí)段波場(chǎng)反變換算法，該方法認(rèn)為計(jì)算過(guò)程中對(duì)方程組階數(shù)影響較大的參數(shù)主要有兩個(gè)，一是積分系數(shù)的數(shù)量，二是對(duì)數(shù)值積分離散時(shí)使用的采樣點(diǎn)的數(shù)量，只要將這2個(gè)量降低，就能把方程組階數(shù)降至較穩(wěn)定的水平。這種辦法提出根據(jù)采樣時(shí)間將數(shù)據(jù)分段，對(duì)每一段數(shù)據(jù)來(lái)說(shuō)，其積分系數(shù)和采樣點(diǎn)個(gè)數(shù)均較少，可以得到較為穩(wěn)定的計(jì)算結(jié)果[28]。

圖4就是在均勻介質(zhì)條件下采用分時(shí)段波場(chǎng)反變換計(jì)算的結(jié)果，介質(zhì)電阻率為100 Ω·m，圖中縱坐標(biāo)A表示歸一化的波場(chǎng)振幅，由圖4可以看出:由于模型不存在地電界面，7個(gè)時(shí)段中的黑色直線(xiàn)即為理論上應(yīng)得到的波場(chǎng)數(shù)據(jù)，點(diǎn)狀線(xiàn)是按照分時(shí)段波場(chǎng)反變換方法計(jì)算得到的結(jié)果，顯然2者存在誤差，且前后時(shí)段的時(shí)間重疊部分的反變換結(jié)果也存在差異，7個(gè)時(shí)段波場(chǎng)數(shù)據(jù)的拼接也存在問(wèn)題(圖5)。

圖4 均勻介質(zhì)分時(shí)段反變換的波場(chǎng)圖像Fig.4 Results of wavefield transformation for timephased in Homogeneous medium

圖5 7段波場(chǎng)反變換結(jié)果疊加Fig.5 Superposition graph of 7 periods results

可以看出，分時(shí)段波場(chǎng)反變換算法只是權(quán)宜之計(jì)，要想更為精確得一次性完成瞬變電磁數(shù)據(jù)的波場(chǎng)反變換計(jì)算，必須使用一種能顯著降低病態(tài)方程條件數(shù)的算法。對(duì)矩陣條件數(shù)進(jìn)行預(yù)處理是一個(gè)可行的辦法，使用超松弛預(yù)條件子能有效將矩陣條件數(shù)降至原來(lái)的平方，從而保證解方程(2)得到可靠且穩(wěn)定的計(jì)算結(jié)果[26]。

該預(yù)條件全時(shí)段波場(chǎng)反變換算法效果如圖6所示。圖6(a)為在均勻介質(zhì)條件下采用全時(shí)段波場(chǎng)反變換計(jì)算的結(jié)果，介質(zhì)電阻率為100 Ω·m，與圖4,5相比計(jì)算結(jié)果精度更高，并且可以一次將整個(gè)時(shí)段波場(chǎng)結(jié)果計(jì)算出來(lái)，避免了各個(gè)時(shí)段重疊部分的處理。圖6(b)為一個(gè)H型地電模型采用全時(shí)段波場(chǎng)反變換得到的波形，模型參數(shù)見(jiàn)表1。由圖6(b)可以看出:全時(shí)段波場(chǎng)反變換結(jié)果能反映模型設(shè)置的2個(gè)地電層位，但波形較寬，尤其是第2個(gè)波形，展寬現(xiàn)象較為嚴(yán)重，導(dǎo)致前后2個(gè)波形的分界點(diǎn)很難確定。

圖6 全時(shí)段波場(chǎng)反變換的波場(chǎng)圖像Fig.6 Results of wavefield transformation for full time

表1 H型模型參數(shù)Table 1 Parameter of H model

綜合來(lái)看，使用分時(shí)段算法和全時(shí)段算法均可以實(shí)現(xiàn)瞬變電磁數(shù)據(jù)的波場(chǎng)反變換計(jì)算，但2者均有明顯的不足之處，分時(shí)段算法不夠穩(wěn)定，結(jié)果的精度不高，最大的問(wèn)題是各個(gè)時(shí)段反變換得到的波場(chǎng)值在相互重疊的部分并不一致，甚至差異較大，僅通過(guò)平均算法壓制的效果不佳，造成重疊段出現(xiàn)虛假子波;全時(shí)段算法穩(wěn)定，計(jì)算精度相對(duì)較高，沒(méi)有虛假子波的問(wèn)題，但單個(gè)波形寬度較大，在虛擬波場(chǎng)時(shí)間上所包含的波數(shù)會(huì)較少，反映的地電層位界面不清晰，分辨率并不能實(shí)現(xiàn)顯著提高解釋精度的目的。

因此，筆者考慮結(jié)合以上2種算法的優(yōu)勢(shì)，既能通過(guò)時(shí)間分段壓縮波形，提高界面分辨能力，又能保留全時(shí)段算法的穩(wěn)定和精度，還能保證段與段重疊部分的真實(shí)可靠。

具體實(shí)施如圖7所示，首先對(duì)分時(shí)段的波場(chǎng)反變換做一個(gè)改進(jìn)，不再將時(shí)間分為固定段，而是設(shè)置一個(gè)固定長(zhǎng)度的時(shí)間窗口，讓它從第1個(gè)時(shí)間處按一定的時(shí)間步長(zhǎng)向最后1個(gè)時(shí)間處滑動(dòng)，每次滑動(dòng)都對(duì)窗口內(nèi)的瞬變電磁數(shù)據(jù)進(jìn)行波場(chǎng)反變換，同時(shí)也要使用1.2節(jié)所述方法計(jì)算整個(gè)時(shí)段的波場(chǎng)反變換結(jié)果，之后將每個(gè)時(shí)間窗口計(jì)算的波場(chǎng)結(jié)果與全時(shí)段計(jì)算的波場(chǎng)結(jié)果做相關(guān)，根據(jù)相關(guān)性強(qiáng)弱來(lái)決定各個(gè)時(shí)間窗口波場(chǎng)結(jié)果是否與全時(shí)段波場(chǎng)結(jié)果進(jìn)行疊加處理，所有與全時(shí)段結(jié)果相關(guān)性強(qiáng)的時(shí)窗結(jié)果與全時(shí)段結(jié)果的疊加結(jié)果即為最終的波場(chǎng)反變換結(jié)果。

圖7 時(shí)間窗口滑動(dòng)示意Fig.7 Schematic diagram of sliding time window

相關(guān)系數(shù)可以定義為

(4)

式中，w1，w2為兩個(gè)獨(dú)立時(shí)窗反變換得到的虛擬波動(dòng)場(chǎng)結(jié)果，k為它們之間的相關(guān)系數(shù)。顯然，k的值越大，則w1，w2的形態(tài)也就越接近。

這種基于相關(guān)算法的滑動(dòng)時(shí)窗波場(chǎng)反變換算法結(jié)果的好壞，主要取決于時(shí)間窗口的大小。如果時(shí)間窗口設(shè)置太大，即對(duì)整個(gè)時(shí)域分段太少，必然導(dǎo)致每次計(jì)算的不適定性增強(qiáng)，結(jié)果精度降低，同時(shí)各個(gè)時(shí)窗的反變換后波形的相關(guān)性會(huì)很大，起不到壓制偽波形的目的;如果時(shí)間窗口設(shè)置太小，參與計(jì)算的采樣點(diǎn)個(gè)數(shù)太少將導(dǎo)致對(duì)反變換方程的約束不足，由于方程奇異產(chǎn)生的數(shù)據(jù)跳變表現(xiàn)明顯，反而掩蓋了由地下電性變化引起的真實(shí)波場(chǎng)信息。

顯然，使用滑動(dòng)時(shí)窗算法時(shí)，必然有一個(gè)最優(yōu)的時(shí)間窗口存在，它既能使反變換得到的波場(chǎng)最大限度的反映真實(shí)電性變化，又能完美壓制不需要的偽波形。經(jīng)過(guò)大量的正演模擬檢驗(yàn)，得到了最優(yōu)時(shí)間窗口的2個(gè)選取標(biāo)準(zhǔn):① 一個(gè)瞬變電磁時(shí)窗轉(zhuǎn)化得到的虛擬波動(dòng)場(chǎng)時(shí)窗中至少要囊括一個(gè)整波，這使得該窗口包含的采樣點(diǎn)個(gè)數(shù)能夠?qū)Ψ醋儞Q方程產(chǎn)生足夠的約束能力;② 在已經(jīng)符合標(biāo)準(zhǔn)1的情況下，滑動(dòng)時(shí)間窗口要盡可能縮小，這使得后期能夠參與相關(guān)運(yùn)算的窗口數(shù)量足夠多，有效壓制偽波形的出現(xiàn)。

確定了相關(guān)系數(shù)和最優(yōu)時(shí)間窗口，很容易可以實(shí)現(xiàn)滑動(dòng)時(shí)窗波場(chǎng)反變換算法，最終可得到精度很高的波場(chǎng)反變換結(jié)果。

圖8為對(duì)應(yīng)表1參數(shù)的H型地電模型采用滑動(dòng)時(shí)窗波場(chǎng)反變換計(jì)算的結(jié)果，圖8中黑色實(shí)線(xiàn)是理論上應(yīng)得到的波場(chǎng)，點(diǎn)狀線(xiàn)是滑動(dòng)時(shí)窗波場(chǎng)反變換算法計(jì)算得到的波場(chǎng)，可以看出:波場(chǎng)反變換結(jié)果的幅值、位置、波形寬度均與理論波形吻合度很高，僅在波形兩側(cè)有少量震蕩，相比全時(shí)段波場(chǎng)反變換結(jié)果受展寬效應(yīng)影響較小。

圖8 滑動(dòng)時(shí)窗波場(chǎng)反變換的波場(chǎng)圖像Fig.8 Results of Wavefield transformation for sliding time window

顯然基于相關(guān)算法的滑動(dòng)時(shí)窗波場(chǎng)反變換算法綜合了分時(shí)段和全時(shí)段波場(chǎng)反變換算法各自?xún)?yōu)勢(shì)，是通過(guò)瞬變電磁信號(hào)求取虛擬波動(dòng)場(chǎng)信號(hào)的最優(yōu)算法。

1.4 二維“一發(fā)多收”數(shù)據(jù)的動(dòng)校正

常規(guī)瞬變電磁采用中心回線(xiàn)裝置施工，經(jīng)過(guò)1.3節(jié)的工作，其數(shù)據(jù)可轉(zhuǎn)化為類(lèi)似于地震自激自收的一維虛擬波場(chǎng)數(shù)據(jù)，已有較好的反演成像方法。巷道-鉆孔瞬變電磁雖然單次測(cè)量是“一發(fā)一收”，但因?yàn)榘l(fā)射裝置位置、能量均未發(fā)生變化，故整體上類(lèi)似地震勘探中 “一發(fā)多收”施工裝置，本質(zhì)屬于二維裝置，需要對(duì)其轉(zhuǎn)化得到的二維虛擬波場(chǎng)數(shù)據(jù)進(jìn)行動(dòng)校正后方可反演。

二維瞬變電磁“一發(fā)多收”施工與二維地震勘探工作方式相似，只是將地震炮點(diǎn)的炸藥激發(fā)改為用線(xiàn)圈或接地電極激發(fā)，將地震的檢波器改為三分量瞬變電磁接收線(xiàn)圈或電極(圖9)。

圖9 二維瞬變電磁工作裝置Fig.9 2D TEM working device

為研究二維瞬變電磁“一發(fā)多收”施工裝置數(shù)據(jù)的波場(chǎng)特征，設(shè)計(jì)一個(gè)D型模型采用1.2節(jié)中使用過(guò)的時(shí)域有限差分進(jìn)行三維數(shù)值模擬，層厚與電阻率等模型參數(shù)如圖10(a)所示，發(fā)射框邊長(zhǎng)定為100 m，供電電流1 A，發(fā)射框中心位于1號(hào)接收點(diǎn)左側(cè)400 m，接收點(diǎn)距5 m。

圖10 D型模型和起伏界面示意Fig.10 Schematic diagram of the D model and undulating formation

對(duì)模擬的垂直分量感應(yīng)電動(dòng)勢(shì)數(shù)據(jù)做波場(chǎng)反變換，可得如圖11(a)所示的波形，對(duì)比一般地震單炮數(shù)據(jù)可以看出:瞬變電磁虛擬波場(chǎng)的“單炮”記錄與地震反射波有明顯不同，由于波形記錄沿測(cè)點(diǎn)方向斜率很小，并未表現(xiàn)出明顯的雙曲線(xiàn)特征，而更符合地震直達(dá)波的直線(xiàn)特征，而隨測(cè)點(diǎn)距離的增大，虛擬波場(chǎng)數(shù)據(jù)的能量也依次減弱，這與地震單炮記錄一致。

圖11 D型模型二維波場(chǎng)反變換結(jié)果Fig.11 Wavefield transformation result of D model and undulating formation

為檢驗(yàn)圖11(a)中的波場(chǎng)時(shí)距曲線(xiàn)特征是否與直觀感受到的一致，對(duì)反射波峰值點(diǎn)分別進(jìn)行了直線(xiàn)擬合和雙曲線(xiàn)擬合，結(jié)果如圖12所示，圖中紅色實(shí)線(xiàn)為直線(xiàn)的結(jié)果，藍(lán)色虛線(xiàn)為雙曲線(xiàn)的結(jié)果，兩條曲線(xiàn)幾乎重合，差異很小。從擬合度上看，直線(xiàn)擬合度為0.976 8，雙曲線(xiàn)擬合度為0.958 4，2者均超過(guò)0.9，都可以作為虛擬波場(chǎng)的時(shí)距曲線(xiàn)特征，但考慮到直線(xiàn)擬合度更高，且計(jì)算簡(jiǎn)單，后續(xù)研究將以虛擬波場(chǎng)的直線(xiàn)型時(shí)距曲線(xiàn)作為基準(zhǔn)。

圖12 D型模型虛擬波場(chǎng)時(shí)距曲線(xiàn)擬合結(jié)果Fig.12 Fitting results of time curve for virtual wavefield of D model

考慮到射線(xiàn)理論，地震單炮記錄還有一個(gè)重要特性，即反射點(diǎn)對(duì)應(yīng)的是炮點(diǎn)與檢波點(diǎn)的中點(diǎn)，但瞬變電磁波2次場(chǎng)在二維條件下是否也遵循射線(xiàn)理論，仍需要數(shù)值模擬結(jié)果的檢驗(yàn)。設(shè)計(jì)一個(gè)D型起伏地層模型采用1.2節(jié)中使用過(guò)的時(shí)域有限差分進(jìn)行三維數(shù)值模擬，層厚與電阻率等模型參數(shù)如圖10(b)所示，在水平位置100～200 m的層面有一梯形凸起，梯形凸起下邊長(zhǎng)100 m，上邊長(zhǎng)75 m，高度30 m，發(fā)射框邊長(zhǎng)定為100 m，供電電流1 A，發(fā)射框中心位于1號(hào)接收點(diǎn)左側(cè)400 m，接收點(diǎn)距10 m。

對(duì)模擬的垂直分量感應(yīng)電動(dòng)勢(shì)數(shù)據(jù)做波場(chǎng)反變換，可得如圖11(b)所示的波形，瞬變電磁虛擬波場(chǎng)的“單炮”記錄仍然與圖11(a)類(lèi)似，并未表現(xiàn)出明顯的雙曲線(xiàn)特征，而表現(xiàn)為直線(xiàn)特征;隨測(cè)點(diǎn)距離的增大，虛擬波場(chǎng)數(shù)據(jù)的能量也依次減弱。值得注意的是:在11～20號(hào)測(cè)點(diǎn)之間，波形記錄上表現(xiàn)出層位的起伏特征，該位置與圖10(b)所示模型設(shè)計(jì)在水平位置100～200 m的凸起位置吻合，證明波形記錄點(diǎn)與接收點(diǎn)(檢波點(diǎn))一致，這也與常規(guī)地震的單炮記錄不同。

因此，根據(jù)數(shù)值模擬結(jié)果可認(rèn)為瞬變電磁二維擬地震單炮數(shù)據(jù)波形記錄具有兩個(gè)基本特征:① 波形記錄點(diǎn)即接收點(diǎn)，即接收點(diǎn)處的虛擬波場(chǎng)記錄反映的就是該接收點(diǎn)正下方的地電信息;② 隨著收發(fā)距的增大，同一層位的反射波記錄時(shí)間也逐步變大，時(shí)距曲線(xiàn)特征與直線(xiàn)特征相符。按照這兩個(gè)特征，再參考地震動(dòng)校正的計(jì)算方法，可以推導(dǎo)瞬變電磁二維擬地震單炮數(shù)據(jù)的動(dòng)校正計(jì)算方法如下。

假設(shè)收發(fā)距為x的接收點(diǎn)波形記錄時(shí)間為t(x)，收發(fā)距為0的接收點(diǎn)波形記錄時(shí)間為t(0)，虛擬波場(chǎng)波速為v，可以得到

t(x)=t(0)+x/v

(5)

動(dòng)校正量可寫(xiě)為

Δt=t(x)-t(0)=x/v

(6)

虛擬波場(chǎng)波速v可由下一節(jié)中的計(jì)算方法獲得。

得到Δt后，對(duì)于收發(fā)距為x的接收點(diǎn)來(lái)說(shuō)，從記錄時(shí)間中減去該值即可完成動(dòng)校正。

對(duì)于如圖10(a)所示D型模型，按照上方確定的動(dòng)校正算法進(jìn)行校正，可以由圖11(a)得到圖13(a)動(dòng)校正后的波場(chǎng)結(jié)果圖，可以看到不同接收點(diǎn)反映層位信息的正峰波形記錄時(shí)間不再傾斜，呈水平層狀特征，與設(shè)計(jì)的數(shù)值模型一致，說(shuō)明動(dòng)校正方法有效。

圖13 D型模型和起伏地層模型動(dòng)校正結(jié)果Fig.13 Dynamic correction result of D model and undulating formation

對(duì)于如圖10(b)所示D型起伏地層模型，同樣按照上方確定的動(dòng)校正算法進(jìn)行校正，可以由圖11(b)得到圖13(b)的動(dòng)校正后的波場(chǎng)結(jié)果圖，可以看到不同接收點(diǎn)的反映層位信息的正峰波形記錄時(shí)間基本不再表現(xiàn)出整體傾斜的特征，除11～12號(hào)測(cè)點(diǎn)之間的起伏特征外，其他測(cè)點(diǎn)大體呈水平層狀特征。而11～20號(hào)測(cè)點(diǎn)位置與模型設(shè)計(jì)在水平位置100～200 m的凸起位置吻合，說(shuō)明動(dòng)校正方法有效。

1.5 初始速度模型

根據(jù)文獻(xiàn)[31]，可以從等效導(dǎo)電平面法出發(fā)，建立虛擬波場(chǎng)速度與電阻率的轉(zhuǎn)換關(guān)系。

等效導(dǎo)電平面法是根據(jù)視縱向電導(dǎo)曲線(xiàn)的特征值直觀地劃分地層的一種近似方法，該方法可以形象地理解為:隨著時(shí)間t的增減，等效導(dǎo)電平面以一定的速度上下“浮動(dòng)”?？烧J(rèn)為此速度即為虛擬波場(chǎng)的傳播速度，可由如下公式計(jì)算

(7)

式中,σi為電導(dǎo)。

1.6 全波形反演

經(jīng)過(guò)1.4節(jié)的動(dòng)校正，二維“一發(fā)多收”瞬變電磁數(shù)據(jù)已轉(zhuǎn)化為類(lèi)似于地震自激自收的一維波形數(shù)據(jù)，可參考文獻(xiàn)[32]對(duì)其進(jìn)行全波形反演，實(shí)現(xiàn)擬地震反演成像。反演迭代公式為

(8)

式中，vp為縱波速度;E為誤差泛函;α為迭代步長(zhǎng);k為迭代次數(shù)，運(yùn)算時(shí)以目標(biāo)函數(shù)梯度的反向作為迭代方向，在第k次的模型上把誤差函數(shù)展開(kāi)就可以得到迭代步長(zhǎng)。

2 模擬檢驗(yàn)

2.1 數(shù)值模擬

為驗(yàn)證巷道-鉆孔瞬變電磁二維工作方法的探測(cè)效果，設(shè)計(jì)如圖14(a)所示的三維數(shù)值模型。在迎頭前方有一超前鉆孔，采空區(qū)異常體中心位于鉆孔孔深37 m處左側(cè)30 m，模型參數(shù)見(jiàn)表2。發(fā)射線(xiàn)圈邊長(zhǎng)為2 m，發(fā)射電流強(qiáng)度為1 A，測(cè)線(xiàn)布置如圖14(b)所示，鉆孔中測(cè)線(xiàn)上有73個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)間距1 m。

對(duì)正演數(shù)據(jù)進(jìn)行二維波場(chǎng)反變換后，可得到圖15(a)所示的波形圖，可以清晰看出虛擬波場(chǎng)對(duì)采空區(qū)異常體前后兩個(gè)地電界面的響應(yīng)，只是該波形沿測(cè)點(diǎn)方向有時(shí)間上的后移。經(jīng)過(guò)動(dòng)校正后可以得到圖15(b)所示的波形，圖中采空區(qū)左右界面的兩層波形已基本校平，與模型設(shè)置基本吻合，可以進(jìn)行下一步全波形反演工作。

圖14 模型示意Fig.14 Schematic diagram of the mode

表2 模型參數(shù)Table 2 Model parameter table

圖15 三維數(shù)值模型數(shù)據(jù)虛擬波場(chǎng)波形Fig.15 Waveform of virtual wavefield for 3D numerical model

對(duì)動(dòng)校正后波形進(jìn)行反演，可得如圖16所示的結(jié)果?？梢钥闯鲈诳咨?5～50 m內(nèi)，鉆孔徑向15～40 m內(nèi)有一明顯低阻異常，異常體形狀接近正方形，邊長(zhǎng)約為25 m，這一低阻異常的厚度、位置均與模型設(shè)置吻合較好，塊狀體的特征反映非常明顯，解釋精度較高。

圖16 三維數(shù)值模型數(shù)據(jù)反演結(jié)果Fig.16 Inversion results of 3D numerical model

2.2 物理模擬

為進(jìn)一步驗(yàn)證方法對(duì)巷道-鉆孔二維施工實(shí)際數(shù)據(jù)的探測(cè)效果，在長(zhǎng)安大學(xué)地球物理專(zhuān)用的物理模擬實(shí)驗(yàn)水槽進(jìn)行了模擬試驗(yàn)，模型設(shè)置如圖17所示。采用多匝小線(xiàn)框激發(fā)、多匝小線(xiàn)框接收的施工方式取得數(shù)據(jù)，發(fā)射線(xiàn)圈邊長(zhǎng)為0.4 m，匝數(shù)為10匝，接收線(xiàn)圈面積約為0.9 m2，發(fā)射電流強(qiáng)度為1.5 A，施工布置如圖17所示，測(cè)線(xiàn)上有14個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)間距為0.05 m。異常體用一銅板代替，銅板規(guī)格為0.1 m×0.1 m×0.002 m，銅板布設(shè)在圖17中坐標(biāo)系相對(duì)發(fā)射、接收線(xiàn)圈中心軸線(xiàn)位置的第1象限內(nèi)，垂直測(cè)線(xiàn)放置，8號(hào)測(cè)點(diǎn)徑向約0.2 m位置。

圖17 模型施工示意Fig.17 Schematic diagram of the model

對(duì)圖17中的觀測(cè)數(shù)據(jù)進(jìn)行二維波場(chǎng)反變換后，可得到圖18(a)所示的波形圖，可以清晰看到8號(hào)測(cè)點(diǎn)處有一組反映低阻體前后兩個(gè)界面的強(qiáng)幅值波形，較其他測(cè)點(diǎn)波幅差距很大，顯然是銅板異常的響應(yīng)，只是該波形對(duì)應(yīng)的時(shí)間較晚，直接反演得到的深度會(huì)比模型布置的參數(shù)偏深。經(jīng)過(guò)動(dòng)校正后可以得到圖18(b)所示的波形，圖中主要異常處的波形時(shí)間經(jīng)過(guò)了校準(zhǔn)，可以進(jìn)行下一步全波形反演工作。

圖18 水槽物理模型數(shù)據(jù)虛擬波場(chǎng)波形Fig.18 Waveform of virtual wavefield for flume physical model

對(duì)動(dòng)校正后波形進(jìn)行反演，可得如圖19所示的結(jié)果?？梢钥闯鲈跍y(cè)線(xiàn)延伸方向長(zhǎng)度0.35 m處得到一明顯的豎條狀低阻異常，沿徑向長(zhǎng)度約為0.1 m，與銅板寬度一致，距離發(fā)射、接收線(xiàn)圈中心軸線(xiàn)的距離約為0.2 m，與銅板放置位置吻合較好。

圖19 水槽物理模型數(shù)據(jù)反演結(jié)果Fig.19 Inversion results of flume physical model

水槽模型相對(duì)比較理想，為接近實(shí)際工況，選擇山西大同某煤礦井下開(kāi)展了1∶1物理模擬試驗(yàn)，施工布置如圖20所示，鉆孔為裸孔，無(wú)套管，使用不導(dǎo)電碳纖維桿人工將探頭送入鉆孔。采用多匝小線(xiàn)框激發(fā)、孔中探頭接收的施工方式取得數(shù)據(jù)，發(fā)射線(xiàn)圈邊長(zhǎng)為2 m，匝數(shù)為40匝，接收探頭等效面積約為1 000 m2，發(fā)射電流強(qiáng)度為10 A，施工鉆孔深72 m，孔中從3 m開(kāi)始有24個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)間距為3 m。認(rèn)為鉆孔左側(cè)掘進(jìn)工作面前方放置的掘進(jìn)機(jī)和膠帶為異常體，掘進(jìn)機(jī)在鉆場(chǎng)前方10～15 m范圍，距離鉆孔徑向距離約5 m。

圖20 施工示意Fig.20 Schematic diagram

對(duì)圖20中的觀測(cè)數(shù)據(jù)進(jìn)行二維波場(chǎng)反變換后，可得到圖21(a)所示的波形圖，可以清晰看到前5個(gè)測(cè)點(diǎn)上有明顯的異常體波形響應(yīng)，與掘進(jìn)機(jī)、皮帶的位置對(duì)應(yīng)較好，且波形記錄清晰表現(xiàn)出“一發(fā)多收”裝置的時(shí)距曲線(xiàn)特征。經(jīng)過(guò)動(dòng)校正后可以得到圖21(b)所示的波形，圖中對(duì)應(yīng)掘進(jìn)機(jī)和膠帶的異常波形已基本校平，與現(xiàn)場(chǎng)實(shí)際情況基本吻合，可以進(jìn)行下一步全波形反演工作。

圖21 井下物理模型數(shù)據(jù)虛擬波場(chǎng)波形Fig.21 Waveform of virtual wavefield for mine physical model

對(duì)動(dòng)校正后波形進(jìn)行反演，可得如圖22所示的結(jié)果?？梢钥闯鲈跍y(cè)線(xiàn)徑向5 m處、延伸方向長(zhǎng)度3～16 m范圍內(nèi)有明顯條帶狀低阻異常，其中10～16 m范圍內(nèi)異常幅值最強(qiáng)，與掘進(jìn)機(jī)位置基本一致，后方弱異常則與膠帶位置吻合較好。

圖22 井下物理模型數(shù)據(jù)反演結(jié)果Fig.22 Inversion results of mine physical model

3 探測(cè)實(shí)例

山西陽(yáng)泉某煤礦15103工作面掘進(jìn)階段，其回風(fēng)巷道掘進(jìn)前方為一已關(guān)閉煤礦。經(jīng)調(diào)查該礦存在越界開(kāi)采行為，且采空區(qū)存在大量積水。由于越界開(kāi)采具體范圍無(wú)法摸清，采空積水成為該回風(fēng)巷道掘進(jìn)過(guò)程中的重要安全隱患。15103工作面回風(fēng)巷道設(shè)計(jì)長(zhǎng)度1 020 m，已掘進(jìn)750 m，根據(jù)調(diào)查推測(cè)再向前掘進(jìn)60 m極有可能揭露鄰礦越界開(kāi)采的積水采空巷道(圖23)。

圖23 15103工區(qū)平面示意Fig.23 Schematic diagram of 15103 work area

為避免揭露采空巷道，提前對(duì)工作面布局進(jìn)行調(diào)整，利用掘進(jìn)工作面已有鉆孔設(shè)計(jì)了巷孔瞬變電磁探測(cè)工作，具體工作設(shè)計(jì)如圖24所示，施工時(shí)每個(gè)鉆孔發(fā)射線(xiàn)框中心法線(xiàn)方向均指向鉆孔鉆進(jìn)方向。施工選擇2號(hào)、3號(hào)和7號(hào)鉆孔施工，其中2號(hào)、3號(hào)鉆孔為順煤層水平孔，7號(hào)鉆孔為斜向上頂板孔，鉆孔均為裸孔，無(wú)套管，使用不導(dǎo)電碳纖維桿人工將探頭送入鉆孔?？字惺┕c(diǎn)距均為3 m，其中2號(hào)孔孔深75 m，實(shí)際測(cè)量72 m，3號(hào)孔孔深75 m，實(shí)際測(cè)量75 m，7號(hào)孔孔深65 m，實(shí)際測(cè)量63 m。

圖24 施工示意Fig.24 Schematic diagram

因篇幅關(guān)系，僅列出3號(hào)鉆孔觀測(cè)數(shù)據(jù)的二維波場(chǎng)反變換結(jié)果和動(dòng)校正結(jié)果，如圖25所示。由圖25(a)中可以清晰看出自10號(hào)測(cè)點(diǎn)開(kāi)始出現(xiàn)束狀波形，符合采空巷道特征，但波形隨收發(fā)距增大出現(xiàn)明顯時(shí)間延遲現(xiàn)象，需要?jiǎng)有Ｕ幚?。?jīng)過(guò)動(dòng)校正后得到圖25(b)所示的波形，時(shí)間延遲現(xiàn)象依然較為明顯，結(jié)合施工巷道、鉆孔以及推測(cè)采空巷道的空間位置分析，此時(shí)的時(shí)間差異應(yīng)與實(shí)際采空巷道位置相關(guān)，可以進(jìn)行下一步全波形反演工作。

圖25 3號(hào)鉆孔數(shù)據(jù)虛擬波場(chǎng)波形Fig.25 Waveform of virtual wavefield for No.3 borehole

對(duì)3個(gè)施工鉆孔動(dòng)校正后波形進(jìn)行反演，可得如圖26所示的結(jié)果?？梢钥闯?個(gè)圖上均有明顯的條帶狀低阻異常，符合采空巷道特征。該異常均出現(xiàn)在3個(gè)鉆孔孔深30 m之后，異常走向與2號(hào)孔基本平行。異常幅值與該異常和鉆孔的距離存在顯著相關(guān)性，圖26(a)中的異常幅值最強(qiáng)，圖26(b),(c)中的異常幅值較弱。

圖26 3個(gè)鉆孔實(shí)測(cè)數(shù)據(jù)反演結(jié)果Fig.26 Inversion results of measured data from 3 boreholes

綜合3個(gè)鉆孔的反演結(jié)果，通過(guò)交匯定位可基本確定低阻異常體位置，得到如圖27所示的積水采空巷道平面解釋結(jié)果圖。相比調(diào)查結(jié)果，積水采空巷道位置更偏北，角度也向北偏轉(zhuǎn)更大，距離掘進(jìn)工作面距離比預(yù)計(jì)的60 m小一半，僅有30 m。

圖27 15103工區(qū)探測(cè)結(jié)果平面示意Fig.27 Schematic diagram of detection results for 15103 work area

本次物探工作結(jié)果提交后，礦方邀請(qǐng)專(zhuān)家召開(kāi)了會(huì)議討論，認(rèn)為物探成果可靠，立即采取措施，停止掘進(jìn)，做好防治水準(zhǔn)備工作，并重新設(shè)計(jì)了開(kāi)切眼和工作面?？梢?jiàn)，巷道-鉆孔瞬變電磁二維擬地震反演方法為礦方安全掘進(jìn)提供了技術(shù)支撐，為礦方下一步工作的決策提供了可靠依據(jù)。

4 結(jié) 論

(1)在距離異常體較近位置觀測(cè)可獲得更強(qiáng)的異常體響應(yīng)，巷道-鉆孔瞬變電磁工作裝置有利于提高低阻異常體的探測(cè)精度，但測(cè)量曲線(xiàn)已不再符合單調(diào)衰減特征，其視電阻率定義非常困難。

(2)傳統(tǒng)的分時(shí)段波場(chǎng)反變換和全時(shí)段波場(chǎng)反變換算法均能將瞬變電磁數(shù)據(jù)轉(zhuǎn)換為虛擬波場(chǎng)數(shù)據(jù)，且各有長(zhǎng)處和不足，基于相關(guān)疊加的滑動(dòng)時(shí)窗波場(chǎng)反變換算法結(jié)合了2者的優(yōu)勢(shì)，又規(guī)避了各自的缺陷，波場(chǎng)轉(zhuǎn)換效果更優(yōu)。

(3)“一發(fā)多收”二維裝置下瞬變電磁虛擬波場(chǎng)數(shù)據(jù)特征與地震單炮記錄特征不同，不是雙曲線(xiàn)形態(tài)，而是直線(xiàn)形態(tài)，且記錄點(diǎn)與接收點(diǎn)相同，據(jù)此得出的動(dòng)校正算法，可將“一發(fā)多收”二維裝置虛擬波場(chǎng)數(shù)據(jù)轉(zhuǎn)化為“自激自收”一維數(shù)據(jù)，之后可采用較為成熟的全波形反演方法對(duì)其實(shí)現(xiàn)高精度的擬地震反演成像。

(4)數(shù)值模擬和物理模擬對(duì)本文提出的方法做了充分的檢驗(yàn)，反演結(jié)果與模型吻合度較高，說(shuō)明巷道-鉆孔瞬變電磁擬地震反演方法能夠突出鉆孔周?chē)惓Ｌ卣鳎瑪U(kuò)大鉆孔有效控制范圍。

(5)井下巷道空間內(nèi)的探測(cè)實(shí)例通過(guò)多個(gè)鉆孔的綜合探測(cè)，精細(xì)解釋了一條威脅礦方安全生產(chǎn)的積水采空巷道，通過(guò)與前期調(diào)查結(jié)果的對(duì)比，證明巷道-鉆孔瞬變電磁擬地震反演方法能提高采空巷道這類(lèi)小目標(biāo)體的探測(cè)準(zhǔn)確率，對(duì)煤礦安全生產(chǎn)決策具有重要意義。