基于鉆孔TEM智能立體成像的快速掘進(jìn)超前探測(cè)方法

范 濤,張幼振，趙 睿,劉 磊,李博凡，郭建磊，李宇騰，田小超，蔣必辭

(1.中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710077; 2.陜西省重點(diǎn)科技創(chuàng)新團(tuán)隊(duì)(地球物理探測(cè)技術(shù)與裝備創(chuàng)新團(tuán)隊(duì))，陜西 西安 710077; 3.煤炭行業(yè)工程研究中心(物探技術(shù)與裝備)，陜西 西安 710077)

采掘失調(diào)是煤礦安全生產(chǎn)的重大隱患之一，而巷道快速掘進(jìn)是保障采掘平衡的重要手段[1]。煤礦井下巷道掘進(jìn)需要超前預(yù)報(bào)前方的水害隱患和地質(zhì)構(gòu)造。據(jù)統(tǒng)計(jì)，近年來(lái)我國(guó)煤礦重特大事故總體呈下降趨勢(shì)，但重特大水災(zāi)事故起數(shù)和死亡人數(shù)在事故占比中逐年增加，平均占比達(dá)18.05%，掘進(jìn)工作面是煤礦重特大水災(zāi)事故最易發(fā)生的突水地點(diǎn)，占比達(dá)51.16%[2]。

目前，巷道快速掘進(jìn)日進(jìn)尺可達(dá)50 m以上，而礦井地球物理超前探測(cè)的有效距離僅80～100 m，探測(cè)精度受井下復(fù)雜環(huán)境影響也較低，因此探掘接續(xù)緊張矛盾已成為制約巷道掘進(jìn)速度的最重要影響因素。

為了解決煤礦井下掘進(jìn)工作面前方的探測(cè)深度與探測(cè)精度的矛盾，利用定向鉆孔進(jìn)行瞬變電磁探測(cè)工作，可以在掘進(jìn)前開展遠(yuǎn)距離、高精度的隱伏水害超前預(yù)報(bào)，該方法的優(yōu)勢(shì)在于:① 發(fā)射回線在目標(biāo)體附近激發(fā)，能在防爆限制下最大限度激發(fā)煤層附近的目標(biāo)體;② 孔中接收裝置則既避開了巷道中的鐵磁性干擾，還能最大限度減少目標(biāo)體2次場(chǎng)因距離帶來(lái)的能量損耗;③ 通過(guò)對(duì)三分量數(shù)據(jù)的聯(lián)合處理可最大程度利用單鉆孔實(shí)現(xiàn)電阻率立體成像，提升瞬變電磁方法的解釋精度，超前準(zhǔn)確預(yù)報(bào)隱伏水體的位置和規(guī)模等空間分布信息;④ 一次性完成大于500 m超前探測(cè)，通過(guò)開鉆窩接續(xù)施工可以保障快速掘進(jìn)過(guò)程中不必停工等待短距離物探超前探測(cè)工作。

早在20世紀(jì)70年代，國(guó)外學(xué)者就已開展了利用鉆孔瞬變電磁三分量探測(cè)信號(hào)解釋孔旁地質(zhì)異常的研究，WOODS等[3]進(jìn)行了比例模型實(shí)驗(yàn)研究，總結(jié)出一套解釋板體模型不同參數(shù)的特征關(guān)系曲線，MACNAE等[4]闡述了導(dǎo)電背景中的井中響應(yīng)符號(hào)變化現(xiàn)象和特征，BUSELLI等[5]模擬了導(dǎo)電覆蓋層下多個(gè)目標(biāo)體的信號(hào)響應(yīng),KOZHEVNIKOV等[6]研究了鉆孔套管對(duì)孔中瞬變電磁響應(yīng)的影響。我國(guó)學(xué)者自20世紀(jì)80年代引入地-井瞬變電磁裝備后也開展了相關(guān)研究，胡平和石中英[7]開展了基于自由空間的球體和板體的地-井瞬變電磁響應(yīng)的理論計(jì)算，對(duì)國(guó)外已報(bào)導(dǎo)的結(jié)果進(jìn)行了補(bǔ)充，張杰[8]推導(dǎo)了矩形回線在空間任意點(diǎn)處產(chǎn)生的一次場(chǎng)表達(dá)式，提出三分量數(shù)據(jù)矢量交匯技術(shù)，楊毅等[9]提出基于導(dǎo)電薄板等效渦流的異常反演方法。孟慶鑫等[10]通過(guò)大地介質(zhì)影響下的正演模擬確定了圍巖背景場(chǎng)對(duì)于總響應(yīng)的影響結(jié)果，徐正玉等[11-13]采用時(shí)域有限差分法模擬研究了接觸帶埋深位置不同和接觸面兩側(cè)電阻率不同對(duì)信號(hào)的影響，楊海燕等[14]研究了覆蓋層影響下板狀體異常響應(yīng)規(guī)律，武軍杰等[15]定義了電性源地-井瞬變電磁全域視電阻率，陳衛(wèi)營(yíng)等[16]對(duì)電性源在地下激發(fā)的6個(gè)電磁場(chǎng)分量的擴(kuò)散、分布特性和探測(cè)能力進(jìn)行了分析研究。在隧/巷道內(nèi)工作的鉆孔瞬變電磁方法近幾年才被提出，相關(guān)研究資料較少，國(guó)外只有VELLA[17]曾將地-井瞬變電磁發(fā)射線圈移到金屬礦巷道中來(lái)探測(cè)含金塊狀黃鐵礦體，國(guó)內(nèi)王世睿[18]研究了隧道10 m以內(nèi)淺孔中的瞬變電磁響應(yīng)特征，提出利用移動(dòng)掘進(jìn)工作面上發(fā)射線框位置來(lái)定性判斷異常體方位的施工技術(shù)，孫懷鳳等[19]通過(guò)物理模擬試驗(yàn)證明了孔中瞬變電磁信號(hào)可用于判斷隧道掘進(jìn)工作面前方是否存在異常構(gòu)造，陳丁等[20]通過(guò)在全空間一維背景上增加三維異常體的積分方程數(shù)值模擬研究了煤礦巷道垂直孔中瞬變電磁特性，范濤[21-23]研究了鉆孔瞬變電磁的疊加超前探測(cè)方法、徑向探測(cè)數(shù)據(jù)的二維擬地震反演方法和短直鉆孔旁裂縫的偽立體成像方法。

綜合以上參考文獻(xiàn)可知，水平分量的形態(tài)組合和幅值差異對(duì)孔旁異常體位置敏感，結(jié)合垂直分量反演成像結(jié)果可對(duì)孔旁異常進(jìn)行立體解釋。但是，根據(jù)水平分量異常形態(tài)組合確定異常體所在象限需人工進(jìn)行識(shí)別和判斷，效率較低，尤其當(dāng)測(cè)點(diǎn)較多時(shí)，人工逐點(diǎn)逐測(cè)道識(shí)別異常曲線形態(tài)更是不現(xiàn)實(shí)的工作。因此，鉆孔瞬變電磁法當(dāng)前的立體解釋還處于定性水平，有必要引入機(jī)器學(xué)習(xí)中的聚類算法實(shí)現(xiàn)智能識(shí)別水平分量異常形態(tài)。

聚類算法在地震勘探領(lǐng)域應(yīng)用較多，刁桂苓等[24]利用系統(tǒng)聚類對(duì)海城地震序列中的24個(gè)震源機(jī)制解進(jìn)行了聚類分析，王偉濤和王寶善[25]基于層次聚類分析有效識(shí)別了汶川余震序列中的相似地震以及重復(fù)地震，張巖等[26]應(yīng)用結(jié)構(gòu)聚類字典學(xué)習(xí)有效壓制了地震數(shù)據(jù)隨機(jī)噪聲;在重磁資料的處理解釋中，張新兵等[27]提出了一種基于改進(jìn)K-means聚類分析的重磁局部異常自動(dòng)圈定方法，李斐等[28]基于聚類分析結(jié)果來(lái)優(yōu)化重力數(shù)據(jù)在不同區(qū)域的觀測(cè)密度，曹書錦等[29]引入自適應(yīng)模糊聚類算法實(shí)現(xiàn)了準(zhǔn)確確定多異常源;而在電磁數(shù)據(jù)處理解釋領(lǐng)域，在大地電磁方法中聚類分析應(yīng)用相對(duì)較多，楊生和楊彥峰[30]將其用于大地電磁曲線分類中，有效克服地質(zhì)推斷的多解性，李晉等[31]提出基于遞歸分析和聚類的信噪比辨識(shí)及分離算法，改善了低頻段的MT數(shù)據(jù)質(zhì)量，黃穎等[32]使用K-means聚類對(duì)MT三維反演結(jié)果中的地質(zhì)構(gòu)造進(jìn)行了識(shí)別和歸類。

參考以上資料，筆者通過(guò)機(jī)器學(xué)習(xí)中的聚類算法對(duì)大數(shù)據(jù)量水平分量異常響應(yīng)曲線形態(tài)進(jìn)行自動(dòng)分類，完成異?？臻g角度定位，并建立其與反演深度的映射關(guān)系，最終實(shí)現(xiàn)鉆孔瞬變電磁立體成像解釋。

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

1 孔旁異常體定位基本原理

1.1 不同象限異常體的水平分量特征

設(shè)計(jì)如圖1所示模型，發(fā)射線圈中心法線方向與鉆孔延伸方向(Z方向)一致，接收線圈中心法線方向與X,Y,Z正方向一致。規(guī)定X分量與Y分量正方向之間區(qū)域?yàn)榈?象限，順時(shí)依次定義為2,3,4象限，在Z=50 m處，分別放置16個(gè)規(guī)模為20 m×20 m×6 m的水平長(zhǎng)方異常體，異常體中心點(diǎn)組成的正方形邊長(zhǎng)為30 m，模型中煤層、巷道和異常體的電阻率分別為1 000,10 000,10 Ω·m。采用文獻(xiàn)[23]中的方法提取水平分量異常后，繪制16個(gè)長(zhǎng)方體的水平分量異常響應(yīng)多測(cè)道圖如圖2所示，圖中橫坐標(biāo)h為鉆孔深度。

圖2 模型的水平分量異常響應(yīng)Fig.2 Horizontal component response of the model

由圖2可知，以水平鉆孔為參考系，所有水平分量異常響應(yīng)形態(tài)均為“正弦曲線”或“反向正弦曲線”形態(tài)，且當(dāng)異常中心與鉆孔的連線與坐標(biāo)軸夾角為0°時(shí)，與該坐標(biāo)軸方向?qū)?yīng)的水平分量的響應(yīng)幅值達(dá)到最小。兩組水平分量形態(tài)組合與異常體所在象限之間的關(guān)系如圖3所示。

圖3 異常位于不同象限時(shí)水平分量響應(yīng)形態(tài)Fig.3 Horizontal component response curve in different quadrant

1.2 異常體XOY平面旋轉(zhuǎn)角計(jì)算方法

鉆孔瞬變電磁徑向探測(cè)時(shí)，異常體引起的二次場(chǎng)是矢量場(chǎng)，那么由水平渦流場(chǎng)的空間分布特征可知，在鉆孔中觀測(cè)到兩個(gè)水平分量Vx,Vy的矢量和Vxy，其方向一定是由鉆孔指向異常體的等效渦流中心上，那么只需要求出Vxy的方向，就知道異常體中心的具體方位。

如圖4所示，設(shè)Vxy與X軸夾角為θ，則

(1)

其中，Vx,Vy均為已知值，求反正弦即得到

(2)

圖4 異常偏轉(zhuǎn)角示意Fig.4 Schematic diagram of abnormal deflection angle window

最后根據(jù)異常體所在象限，可由如下公式求出對(duì)應(yīng)的XOY平面旋轉(zhuǎn)角α:① 異常體在第1象限:α=θ;② 異常體在第2象限:α=π-θ;③ 異常體在第3象限:α=π+θ;④ 異常體在第4象限:α=2π-θ。

根據(jù)以上內(nèi)容，可以看出，若將鉆孔旁空間中任意一點(diǎn)都視為一個(gè)小的地質(zhì)異常體，理論上可以對(duì)任意一點(diǎn)的電阻率進(jìn)行空間定位，實(shí)現(xiàn)孔旁電性信息立體成像。但顯然，水平分量異常場(chǎng)形態(tài)的確定工作量很大，依靠人工完成難以保證效率，必須引入人工智能算法實(shí)現(xiàn)曲線形態(tài)自動(dòng)分類。

2 孔旁異常體智能立體成像方法

2.1 水平分量異常曲線類型自動(dòng)識(shí)別方法

欲對(duì)水平分量異常場(chǎng)曲線進(jìn)行自動(dòng)分類，應(yīng)以每一個(gè)測(cè)點(diǎn)為插值窗口中心，采用Hermit插值求得所有窗口對(duì)應(yīng)的水平分量異常場(chǎng)。插值窗口大小的選擇可根據(jù)垂直分量中的主要異常區(qū)統(tǒng)計(jì)平均值確定。之后對(duì)提取出的異常數(shù)據(jù)進(jìn)行正規(guī)化，將曲線橫坐標(biāo)范圍統(tǒng)一，然后對(duì)區(qū)間正規(guī)化后的所有異常數(shù)據(jù)進(jìn)行特征提取，提出數(shù)據(jù)中極值對(duì)應(yīng)的正規(guī)化點(diǎn)號(hào)，最后以極大值點(diǎn)號(hào)為X軸，極小值點(diǎn)號(hào)為Y軸，形成特征值分布圖(圖5)。

圖5 特征值分布Fig.5 Distribution of eigenvalues

圖5為1組數(shù)值模擬數(shù)據(jù)和1組實(shí)測(cè)數(shù)據(jù)的水平分量異常場(chǎng)特征值分布，可以看出均具有明顯的二分類特性，與水平分量異常場(chǎng)的“正弦曲線”或“反向正弦曲線”形態(tài)存在顯著相關(guān)性。基于該特性，無(wú)需提前進(jìn)行標(biāo)簽樣本的監(jiān)督訓(xùn)練，可直接選用無(wú)監(jiān)督機(jī)器學(xué)習(xí)算法對(duì)數(shù)據(jù)進(jìn)行分類。

無(wú)監(jiān)督機(jī)器學(xué)習(xí)常常被應(yīng)用在數(shù)據(jù)挖掘領(lǐng)域，用于在大量無(wú)標(biāo)簽數(shù)據(jù)中發(fā)現(xiàn)規(guī)律。它的訓(xùn)練數(shù)據(jù)是無(wú)標(biāo)簽的，訓(xùn)練目標(biāo)是能對(duì)觀察值進(jìn)行分類或區(qū)分等。常用的無(wú)監(jiān)督學(xué)習(xí)算法主要有主成分分析方法、等距映射方法、局部線性嵌入方法、拉普拉斯特征映射方法、黑塞局部線性嵌入方法、局部切空間排列方法和最常用的聚類方法。

聚類算法是指將一堆沒(méi)有標(biāo)簽的數(shù)據(jù)自動(dòng)劃分成幾類的方法，這個(gè)方法要保證同一類的數(shù)據(jù)有相似的特征。筆者選擇K-means聚類算法，該算法是使用最大期望算法求解的高斯混合模型在正態(tài)分布的協(xié)方差為單位矩陣，且隱變量的后驗(yàn)分布為一組狄拉克δ函數(shù)時(shí)所得到的特例，它假設(shè)相同類別中數(shù)據(jù)之間的距離應(yīng)該都很近，即數(shù)據(jù)之間的相似度與它們之間的歐式距離成反比。

需要將n個(gè)水平分量純異常數(shù)據(jù){xk}聚為2類，令經(jīng)過(guò)聚類之后每個(gè)數(shù)據(jù)所屬的類別為{tk}，而這2個(gè)聚類的中心為{μl}，可定義如下的損失函數(shù)L:

(3)

式中，l為聚類數(shù)目。

實(shí)際計(jì)算時(shí)先隨機(jī)設(shè)置2個(gè)質(zhì)心把所有數(shù)據(jù)粗略分成2個(gè)初始類，計(jì)算所有數(shù)據(jù)與質(zhì)心的歐式距離，再根據(jù)平均值重新計(jì)算質(zhì)心和類別，對(duì)以上過(guò)程反復(fù)迭代，直至達(dá)到終止條件。終止條件可設(shè)置為簇中心點(diǎn)變化率ηk，即

(4)

將分類好的數(shù)據(jù)類別與“正弦曲線”或“反向正弦曲線”形態(tài)進(jìn)行對(duì)應(yīng)，根據(jù)圖3就可以準(zhǔn)確確定鉆孔瞬變電磁觀測(cè)數(shù)據(jù)中任意測(cè)點(diǎn)任意測(cè)道反映的電性信息所在的象限，再按照第1.2節(jié)所述的異常體XOY平面旋轉(zhuǎn)角計(jì)算方法，將每一測(cè)點(diǎn)每一測(cè)道對(duì)應(yīng)的電性信息視為異常體代入計(jì)算，即可獲取對(duì)應(yīng)的XOY平面旋轉(zhuǎn)角度。

2.2 垂直分量一維反演深度與XOY平面旋轉(zhuǎn)角的映射

得到每一測(cè)點(diǎn)每一測(cè)道的XOY平面旋轉(zhuǎn)角后，就需要求取相應(yīng)的反演電性信息。筆者使用的鉆孔瞬變電磁工作方法本質(zhì)上仍屬于中心回線裝置類型，其垂直分量實(shí)測(cè)數(shù)據(jù)曲線形態(tài)與礦井瞬變電磁探測(cè)數(shù)據(jù)曲線形態(tài)基本相同(圖6)，僅是因?yàn)榘l(fā)射線圈尺寸與匝數(shù)的原因而導(dǎo)致電感影響較大，因此數(shù)據(jù)處理方法可參考礦井瞬變電磁，采用文獻(xiàn)[33]中的預(yù)處理技術(shù)對(duì)電感影響進(jìn)行校正，對(duì)校正后的數(shù)據(jù)則可應(yīng)用Occam反演、虛擬波場(chǎng)反演等方法[34-35]進(jìn)行深度和電阻率反演。

圖6 礦井裝置與鉆孔裝置實(shí)測(cè)曲線對(duì)比Fig.6 Comparison of measured curves between mine device and borehole device

但因?yàn)榉囱莸玫降牡貙訉訑?shù)一般遠(yuǎn)小于觀測(cè)時(shí)間道數(shù)，因此由垂直分量反演得到的深度(鉆孔探測(cè)半徑)與采樣時(shí)間之間并不存在一一對(duì)應(yīng)關(guān)系，而計(jì)算得到的XOY平面旋轉(zhuǎn)角與采樣時(shí)間一一對(duì)應(yīng)，因此無(wú)法直接形成反演深度與XOY平面旋轉(zhuǎn)角的映射，限制了反演電阻率的立體化。

為解決這一問(wèn)題，首先引入瞬變電磁生產(chǎn)中經(jīng)常使用的基于趨膚深度原理的深度計(jì)算公式

(5)

式中，h為反演得到的深度(鉆孔探測(cè)半徑);C為深度系數(shù);ρ為反演得到的電阻率;t為采樣時(shí)間。

根據(jù)式(5)，計(jì)算大量層狀模型，并對(duì)比不同電阻率參數(shù)情況下的層狀模型分界面的計(jì)算深度與模型設(shè)定深度之間的差別，推導(dǎo)出深度系數(shù)C與電阻率ρ的關(guān)系如圖7所示。

圖7 深度系數(shù)與電阻率的關(guān)系Fig.7 Relationship between depth coefficient and resistivity

可以看出，在對(duì)數(shù)坐標(biāo)系下，深度系數(shù)與電阻率呈現(xiàn)線性關(guān)系，通過(guò)擬合可以得到該直線方程為

lgC=-0.500 093 65lgρ-0.001 374 49

(6)

由式(6)可逐層推導(dǎo)出層狀地層情況下反演深度與采樣時(shí)間之間的計(jì)算公式，建立2者之間的映射關(guān)系，進(jìn)而建立反演鉆孔探測(cè)半徑與XOY平面旋轉(zhuǎn)角的一一對(duì)應(yīng)關(guān)系，具體算法如下:

(2)根據(jù)lgC=-0.500 093 65lgρ-0.001 374 49計(jì)算第i層對(duì)應(yīng)的標(biāo)準(zhǔn)系數(shù)Ci;

(3)根據(jù)Ai=lgCi-lgC1,CRi=10Ai計(jì)算第i層對(duì)應(yīng)的相對(duì)系數(shù)CRi;

(4)除第1層外的第i層真深度與采樣時(shí)間可建立:

(7)

即在第1層外的第i層的采樣時(shí)間可以按如下公式計(jì)算:

(8)

此時(shí)由反演得到的深度(鉆孔探測(cè)半徑)和電阻率可反推出相應(yīng)的采樣時(shí)間，再通過(guò)插值可以獲得與實(shí)際采樣時(shí)間道對(duì)應(yīng)的反演鉆孔探測(cè)半徑和電阻率。

2.3 反演結(jié)果的立體成像算法

由鉆孔瞬變電磁探測(cè)的垂直分量可反演得到以鉆孔孔深為橫坐標(biāo)、以鉆孔探測(cè)半徑為縱坐標(biāo)的一維反演電阻率剖面成像圖，在已求得圖中每一測(cè)點(diǎn)每一測(cè)道電性信息對(duì)應(yīng)的XOY平面旋轉(zhuǎn)角信息時(shí)，假設(shè)鉆孔為直鉆孔，可通過(guò)三角函數(shù)關(guān)系將其投影在XY平面上，將剖面圖中每一測(cè)點(diǎn)的一維反演電性信息曲線轉(zhuǎn)換為二維電性信息平面，多個(gè)測(cè)點(diǎn)的二維成像結(jié)果組合即可實(shí)現(xiàn)鉆孔瞬變電磁反演電阻率立體成像。

單個(gè)測(cè)點(diǎn)電性信息一維坐標(biāo)擴(kuò)展至二維的計(jì)算公式為

(9)

式中，j為測(cè)點(diǎn)數(shù);p為測(cè)道數(shù);x為鉆孔瞬變電磁坐標(biāo)系中電性信息對(duì)應(yīng)的X方向坐標(biāo);y為鉆孔瞬變電磁坐標(biāo)系中電性信息對(duì)應(yīng)的Y方向坐標(biāo);r為鉆孔探測(cè)半徑;α為XOY平面旋轉(zhuǎn)角。

2.4 非直鉆孔立體成像成果空間坐標(biāo)校正方法

由于施工空間為定向長(zhǎng)鉆孔，其真實(shí)軌跡并非直線，有必要根據(jù)鉆孔軌跡對(duì)基于直線鉆孔坐標(biāo)系的電阻率立體成像的空間坐標(biāo)進(jìn)行校正，獲取與實(shí)際空間位置完全對(duì)應(yīng)的立體成像成果。

首先要計(jì)算鉆孔軌跡，實(shí)際的鉆孔軌跡是一條空間曲線，計(jì)算鉆孔軌跡坐標(biāo)時(shí)假設(shè)鉆孔軌跡由若干直線段組成，得到每個(gè)測(cè)段的坐標(biāo)增量，然后累加求得坐標(biāo)值進(jìn)行軌跡繪制。

以計(jì)算第2測(cè)量點(diǎn)坐標(biāo)為例說(shuō)明:測(cè)量起始位置處為原點(diǎn)(X1=0,Y1=0,Z1=0)，第2測(cè)量點(diǎn)各坐標(biāo)增量(ΔX1，ΔY1，ΔZ1)可表示為

(10)

γ=arccos[cosβ1cosβ2+sinβ1sinβ2cos(φ2-φ1)]

(11)

式中，ΔL為測(cè)點(diǎn)間距；β1為第1測(cè)量段傾斜角;β2為第2測(cè)量段傾斜角;φ1為第1測(cè)量段傾斜方位角;φ2為第2測(cè)量段傾斜方位角。

即第2個(gè)測(cè)量點(diǎn)處鉆孔軌跡坐標(biāo)為

X2=X1+ΔX1,Y2=Y1+ΔY1,Z2=Z1+ΔZ1

(12)

要進(jìn)一步實(shí)現(xiàn)三維探測(cè)成果的空間展示，還需要對(duì)探測(cè)成果坐標(biāo)系進(jìn)行三維轉(zhuǎn)換。以鉆孔軌跡為旋轉(zhuǎn)軸，以測(cè)點(diǎn)所在孔位的傾角和方位角為旋轉(zhuǎn)角進(jìn)行空間旋轉(zhuǎn)。新、原坐標(biāo)系相對(duì)位置如圖8所示，新坐標(biāo)系原點(diǎn)坐標(biāo)為(X0,Y0,Z0)，相對(duì)原坐標(biāo)系其單位坐標(biāo)矢量關(guān)系為

(13)

圖8 新、原坐標(biāo)系相對(duì)關(guān)系Fig.8 Relative relationship between the new and original coordinate systems

因此，要將三維坐標(biāo)從原XYZ坐標(biāo)系轉(zhuǎn)換成新X′Y′Z′坐標(biāo)系，可由以下兩步實(shí)現(xiàn):

步驟1:平移矩陣。

平移矩陣T可表示為

(14)

步驟2:構(gòu)造坐標(biāo)旋轉(zhuǎn)矩陣。

坐標(biāo)旋轉(zhuǎn)矩陣R可表示為

(15)

其坐標(biāo)變換可用矩陣表示為

(X′,Y′,Z′,1)=(X,Y,Z,1)·T·(-X0,

-Y0,-Z0)·R

(16)

對(duì)已得出的直線鉆孔坐標(biāo)系下的成像結(jié)果中每一點(diǎn)空間坐標(biāo)都進(jìn)行如上變換，即可實(shí)現(xiàn)與實(shí)際鉆孔軌跡對(duì)應(yīng)的電阻率立體成像處理。

3 模擬檢驗(yàn)

3.1 數(shù)值模擬

為驗(yàn)證基于水平分量異常特征聚類的鉆孔瞬變電磁立體成像方法探測(cè)效果，設(shè)計(jì)如圖9所示的三維模型，采用時(shí)域有限差分方法進(jìn)行了數(shù)值模擬。在鉆孔深度方向50 m處，第3象限偏轉(zhuǎn)30°放置1個(gè)規(guī)模為15 m×15 m×15 m的低阻異常體，異常體中心點(diǎn)距離鉆孔30 m，模型中煤層、巷道和異常體的電阻率分別為1 000,10 000,10 Ω·m。

圖9 數(shù)值模型示意Fig.9 Schematic diagram of the mode

對(duì)鉆孔深度20～80 m內(nèi)的測(cè)點(diǎn)數(shù)據(jù)進(jìn)行Occam反演成像可得到如圖10所示的沿鉆孔方向的電阻率剖面，由圖10可知在鉆孔深度Z=50 m、鉆孔徑向R=20～30 m位置有較為明顯的低阻異常響應(yīng)，但從圖10中無(wú)法反映異常的空間方位。

圖10 數(shù)值模型電阻率剖面Fig.10 Resistivity profile of numerical model

采用本文1.2和2.1節(jié)方法對(duì)模型水平分量進(jìn)行處理，再采用2.2和2.3節(jié)方法對(duì)反演電阻率進(jìn)行XOY平面旋轉(zhuǎn)角分配，可得到如圖11所示的單測(cè)點(diǎn)電阻率展開平面圖。由圖11可以看出，在孔深50 m的平面圖中XY平面第3象限有明顯的低阻異常響應(yīng)，其中心點(diǎn)與坐標(biāo)系原點(diǎn)(鉆孔)之間的距離為30 m，形狀規(guī)模與模型參數(shù)基本一致，而在孔深20 m和70 m的平面圖中則沒(méi)有明顯低阻異常顯示。

圖11 單測(cè)點(diǎn)電阻率展開平面Fig.11 Expanded plan view of resistivity of single measuring point

將每一測(cè)點(diǎn)的X,Y坐標(biāo)與鉆孔孔深Z坐標(biāo)組合，通過(guò)Voxler軟件進(jìn)行立體成像如圖12所示。圖中低阻異常體與模型設(shè)置參數(shù)一致，說(shuō)明基于水平分量異常特征聚類的方法對(duì)鉆孔徑向電阻率立體成像有效，準(zhǔn)確性較高。

3.2 物理模擬

為進(jìn)一步驗(yàn)證立體成像方法對(duì)鉆孔瞬變電磁實(shí)際數(shù)據(jù)的探測(cè)效果，在長(zhǎng)安大學(xué)地球物理專用的物理模擬實(shí)驗(yàn)水槽進(jìn)行了模擬試驗(yàn)，模型設(shè)置參考數(shù)值模型，如圖13所示，模型介質(zhì)水和銅板的電阻率分別為40,1.75×10-8Ω·m。采用多匝小線框激發(fā)、多匝小線框接收的施工方式取得數(shù)據(jù)，發(fā)射線圈邊長(zhǎng)為0.4 m，匝數(shù)為10匝，接收線圈面積約為0.9 m2，發(fā)射電流強(qiáng)度為1.5 A，施工布置如圖13所示，測(cè)線上有15個(gè)測(cè)點(diǎn)，測(cè)點(diǎn)間距為0.05 m。在模擬鉆孔深度方向0.35 m處，第3象限偏轉(zhuǎn)30°放置1個(gè)規(guī)模為0.2 m×0.2 m×0.002 m的銅板，銅板邊界距離鉆孔0.1 m。

對(duì)模擬鉆孔深度0～0.7 m的測(cè)點(diǎn)數(shù)據(jù)進(jìn)行Occam反演成像可以得到如圖14所示的沿模擬鉆孔方向的電阻率剖面圖，可以較為清晰的看到在鉆孔深度Z為0.35 m、鉆孔徑向R為0.1～0.3 m位置有較為明顯的低阻異常響應(yīng)，但從該成果圖中無(wú)法反映異常的空間方位。

采用本文1.2和2.1節(jié)方法對(duì)模型水平分量進(jìn)行處理，再采用2.2和2.3節(jié)方法對(duì)反演電阻率進(jìn)行XOY平面旋轉(zhuǎn)角分配，可得到如圖15所示的單測(cè)點(diǎn)電阻率展開平面圖。由圖15可看出，在孔深0.35 m的平面圖中XY平面第3象限有明顯的低阻異常響應(yīng)，其邊界與坐標(biāo)系原點(diǎn)(鉆孔)之間的距離為0.1 m，形狀規(guī)模與銅板參數(shù)基本一致，而在孔深0.15 m和0.55 m的平面圖中則沒(méi)有明顯低阻異常顯示。

將每一測(cè)點(diǎn)的X,Y坐標(biāo)與模擬鉆孔孔深Z坐標(biāo)組合，通過(guò)Voxler軟件進(jìn)行立體成像如圖16所示。

圖12 電阻率立體成像圖Fig.12 Stereo imaging diagram of resistivity

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

圖14 物理模型電阻率剖面Fig.14 Resistivity profile of physical model

圖中低阻異常體與銅板參數(shù)基本一致，說(shuō)明立體成像方法對(duì)鉆孔瞬變電磁實(shí)際數(shù)據(jù)同樣有效，準(zhǔn)確性較高。

圖15 單測(cè)點(diǎn)電阻率展開平面Fig.15 Expanded plan view of resistivity of single measuring point

圖16 電阻率立體成像Fig.16 Stereo imaging diagram of resistivity

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

山西某煤礦所處井田地質(zhì)構(gòu)造復(fù)雜，已開采工作面內(nèi)部陷落柱分布密集，為了保障03工作面回風(fēng)巷掘進(jìn)效率，在巷道開拓前采用定向鉆孔中瞬變電磁技術(shù)對(duì)掘進(jìn)巷道遭遇構(gòu)造情況進(jìn)行長(zhǎng)距離超前探查。

鉆場(chǎng)位于03工作面回風(fēng)聯(lián)絡(luò)巷，偏離03工作面回風(fēng)巷17.7 m,開孔位置上距03號(hào)煤底板4.5 m，下距02號(hào)煤頂板1.38 m，鉆孔深度468 m。根據(jù)已知地質(zhì)信息，兩層煤中間夾層主要是砂巖和泥巖。鉆孔瞬變電磁施工段為孔深12～450 m，測(cè)量點(diǎn)距為3 m，具體如圖17所示。

圖17 03工作面定向鉆孔開孔位置剖面Fig.17 Sectional view of the position of the directional drilling hole in 03 working face

對(duì)鉆孔內(nèi)的測(cè)點(diǎn)數(shù)據(jù)進(jìn)行Occam反演成像可以得到如圖18所示的沿鉆孔方向的電阻率剖面圖，可以較為清晰的看到在鉆孔深度Z為380～450 m、鉆孔徑向R為2～20 m位置有較為明顯的低阻異常響應(yīng)，但從該成果圖中無(wú)法反映異常的空間方位。

采用本文1.2和2.1節(jié)方法對(duì)實(shí)測(cè)數(shù)據(jù)水平分量進(jìn)行處理，再采用2.2,2.3和2.4節(jié)方法對(duì)反演電阻率進(jìn)行XOY平面旋轉(zhuǎn)角分配和空間坐標(biāo)校正，可得到如圖19所示的單測(cè)點(diǎn)電阻率展開平面圖。由圖19可以看出，在孔深432 m的平面圖中XY平面下方兩個(gè)象限均有明顯的低阻異常響應(yīng)，且在中部有一筍狀突起，而在孔深300 m的平面圖中則沒(méi)有明顯低阻異常顯示。

將每一測(cè)點(diǎn)的X,Y坐標(biāo)與鉆孔孔深Z坐標(biāo)組合，通過(guò)Voxler軟件進(jìn)行立體成像如圖20所示。圖中低阻異常體主要分布在鉆孔深度380～450 m，鉆孔下方2～20 m位置，在該區(qū)域，鉆孔自然伽馬測(cè)量值同樣跳動(dòng)較大，且鉆孔施工至433 m時(shí)返渣中出現(xiàn)黃色物質(zhì)，該物質(zhì)未在之前鉆探砂巖返渣中出現(xiàn)，因此綜合推測(cè)此處可能存在陷落柱構(gòu)造。

圖18 實(shí)測(cè)數(shù)據(jù)電阻率剖面Fig.18 Resistivity profile of measured data

圖19 單測(cè)點(diǎn)電阻率展開平面Fig.19 Expanded plan view of resistivity of single measuring point

圖20 電阻率立體成像Fig.20 Stereo imaging diagram of resistivity

根據(jù)本次物探工作結(jié)果，礦方立即對(duì)該疑似陷落柱做了工作預(yù)案，并開展快速掘進(jìn)工作?？梢?，基于定向鉆孔的瞬變電磁立體成像方法為礦方安全快速掘進(jìn)提供了遠(yuǎn)距離地質(zhì)超前預(yù)報(bào)方面的有力技術(shù)支撐。

5 結(jié) 論

(1)以鉆孔鉆進(jìn)方向?yàn)閆軸正方向，以孔口所在平面右向?yàn)閄軸正方向，下向?yàn)閅軸正方向，鉆孔瞬變電磁所有水平分量異常響應(yīng)形態(tài)均為“正弦曲線”或“反向正弦曲線”形態(tài)，且當(dāng)某一水平分量與該水平分量坐標(biāo)軸夾角為0°時(shí)，該水平分量的響應(yīng)幅值達(dá)到最小，因此，通過(guò)X,Y分量異常形態(tài)組合可判定異常體所在象限。

(2)由2組水平分量的幅值基于三角函數(shù)關(guān)系可計(jì)算得到異常體中心在異常所在象限內(nèi)的偏轉(zhuǎn)角度，結(jié)合異常象限可得出異常體中心XOY平面旋轉(zhuǎn)角。

(3)在對(duì)數(shù)坐標(biāo)系下，深度計(jì)算系數(shù)與電阻率呈線性關(guān)系，因此借由深度經(jīng)驗(yàn)公式可推導(dǎo)出深度-電阻率與采樣時(shí)間的映射關(guān)系，實(shí)現(xiàn)反演電阻率與測(cè)點(diǎn)測(cè)道的一一對(duì)應(yīng)。

(4)將測(cè)量數(shù)據(jù)每一測(cè)點(diǎn)每一測(cè)道對(duì)應(yīng)的電阻率視為一個(gè)獨(dú)立的異常體，采用K-means聚類算法對(duì)相應(yīng)的兩組水平分量異常曲線中的極值進(jìn)行二分類，可以自動(dòng)確定任意一點(diǎn)一道電阻率的分布象限，再由水平分量異常場(chǎng)幅值算出XOY平面旋轉(zhuǎn)角，即可結(jié)合垂直分量成像結(jié)果實(shí)現(xiàn)鉆孔徑向電阻率的立體成像。

(5)假設(shè)鉆孔軌跡由若干直線段組成，根據(jù)每個(gè)測(cè)段的坐標(biāo)增量累加可求得坐標(biāo)值得出定向鉆孔軌跡，再以該軌跡為旋轉(zhuǎn)軸，以異常位置角度為旋轉(zhuǎn)角進(jìn)行空間旋轉(zhuǎn)，可實(shí)現(xiàn)非直鉆孔徑向的電阻率立體成像。

(6)數(shù)值模擬和物理模擬對(duì)本文提出的方法做了充分的檢驗(yàn)，反演結(jié)果與模型吻合度較高，說(shuō)明該方法準(zhǔn)確、有效，提高了鉆孔瞬變電磁方法的解釋水平。

(7)井下巷道空間內(nèi)的探測(cè)實(shí)例通過(guò)立體成像遠(yuǎn)距離超前預(yù)報(bào)了一個(gè)威脅礦方安全生產(chǎn)、影響礦方巷道掘進(jìn)的陷落柱，并獲得了鉆孔自然伽馬測(cè)井資料和鉆探返渣的佐證，證明本文方法可以推廣應(yīng)用至煤礦實(shí)際生產(chǎn)中，為掘進(jìn)前工作面前方500 m距離內(nèi)隱伏水害遠(yuǎn)距離精準(zhǔn)超前探測(cè)提供技術(shù)支撐，進(jìn)一步通過(guò)接續(xù)鉆進(jìn)+探測(cè)可有效保障巷道快速掘進(jìn)。