范 濤

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

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

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

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

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

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

1 方法原理

1.1 工作裝置

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

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

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

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

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

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

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

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

1.3 波場反變換

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

(1)

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

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

(2)

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

(3)

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

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

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

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

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

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

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

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

圖6 全時段波場反變換的波場圖像Fig.6 Results of wavefield transformation for full time

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

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

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

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

圖7 時間窗口滑動示意Fig.7 Schematic diagram of sliding time window

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

(4)

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

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

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

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

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

圖8 滑動時窗波場反變換的波場圖像Fig.8 Results of Wavefield transformation for sliding time window

顯然基于相關(guān)算法的滑動時窗波場反變換算法綜合了分時段和全時段波場反變換算法各自優(yōu)勢，是通過瞬變電磁信號求取虛擬波動場信號的最優(yōu)算法。

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

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

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

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

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

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

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

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

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

圖12 D型模型虛擬波場時距曲線擬合結(jié)果Fig.12 Fitting results of time curve for virtual wavefield of D model

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

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

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

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

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

(5)

動校正量可寫為

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

(6)

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

得到Δt后，對于收發(fā)距為x的接收點(diǎn)來說，從記錄時間中減去該值即可完成動校正。

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

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

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

1.5 初始速度模型

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

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

(7)

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

1.6 全波形反演

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

(8)

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

2 模擬檢驗(yàn)

2.1 數(shù)值模擬

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

對正演數(shù)據(jù)進(jìn)行二維波場反變換后，可得到圖15(a)所示的波形圖，可以清晰看出虛擬波場對采空區(qū)異常體前后兩個地電界面的響應(yīng)，只是該波形沿測點(diǎn)方向有時間上的后移。經(jī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ù)虛擬波場波形Fig.15 Waveform of virtual wavefield for 3D numerical model

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

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

2.2 物理模擬

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

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

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

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

對動校正后波形進(jìn)行反演，可得如圖19所示的結(jié)果。可以看出在測線延伸方向長度0.35 m處得到一明顯的豎條狀低阻異常，沿徑向長度約為0.1 m，與銅板寬度一致，距離發(fā)射、接收線圈中心軸線的距離約為0.2 m，與銅板放置位置吻合較好。

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

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

圖20 施工示意Fig.20 Schematic diagram

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

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

對動校正后波形進(jìn)行反演，可得如圖22所示的結(jié)果?？梢钥闯鲈跍y線徑向5 m處、延伸方向長度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 探測實(shí)例

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

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

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

圖24 施工示意Fig.24 Schematic diagram

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

(4)數(shù)值模擬和物理模擬對本文提出的方法做了充分的檢驗(yàn)，反演結(jié)果與模型吻合度較高，說明巷道-鉆孔瞬變電磁擬地震反演方法能夠突出鉆孔周圍異常特征，擴(kuò)大鉆孔有效控制范圍。

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