范 濤,李鴻泰,劉 磊，趙 睿,李博凡,郭建磊，李宇騰

(1. 中煤科工集團西安研究院有限公司，陜西 西安 710077; 2. 長安大學(xué) 地球物理場多參數(shù)綜合模擬實驗室(中國地球物理學(xué)會重點實驗室)，陜西 西安 710054； 3. 四川省地質(zhì)礦產(chǎn)勘查開發(fā)局物探隊，四川 成都 610072)

0 引 言

隨著中國煤炭開發(fā)逐漸向深部發(fā)展，復(fù)雜地質(zhì)構(gòu)造和隱蔽致災(zāi)因素對煤礦安全生產(chǎn)影響加重，尤其是在煤礦巷道掘進過程中，面臨復(fù)雜未知的地質(zhì)條件，對地質(zhì)保障技術(shù)的先進性與精準(zhǔn)性提出更高的要求。2018年頒布的《煤礦防治水細(xì)則》中就對地球物理超前探測技術(shù)提出了新的要求，其中明確提出探測點20 m范圍內(nèi)不得有積水和金屬物體，然而目前主流的超前探測技術(shù)均不符合該項規(guī)定，急需研發(fā)新的超前探測手段[1]。

研究表明，利用掘進工作面的鉆孔進行瞬變電磁探測工作，可以在掘進前開展遠距離、高精度的隱伏水害超前預(yù)報。鉆孔瞬變電磁法基本施工裝置如圖1所示。圖1中，發(fā)射線圈和接收線圈共同集成在孔內(nèi)探管中，通過孔口控制一次場發(fā)射與關(guān)斷，激發(fā)探管附近的低阻異常體產(chǎn)生二次場，再接收二次場的三分量信號，最后解釋得到鉆孔周圍的電性分布。該方法的優(yōu)勢在于：①發(fā)射線圈在目標(biāo)體附近激發(fā)，能在防爆限制下最大限度激發(fā)煤層附近的目標(biāo)體；②孔中接收線圈既避開了巷道中的鐵磁性干擾，還能最大限度減少目標(biāo)體二次場因距離帶來的能量損耗；③通過對水平分量的處理定位孔旁地質(zhì)異常體所在方位，可利用單鉆孔探測資料對異常體進行立體解釋。

圖1 鉆孔瞬變電磁法基本施工裝置示意圖Fig.1 Schematic View of Basic Construction Device forBorehole Transient Electromagnetic Method

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

從前人研究成果可以看出，水平分量形態(tài)組合和幅值差異對孔旁地質(zhì)異常體位置敏感，結(jié)合垂直分量視電阻率成像結(jié)果可對異常體進行立體解釋。但是，根據(jù)水平分量異常形態(tài)組合確定異常體所在象限需人工進行識別和判斷，效率較低，尤其當(dāng)測點較多時，人工逐點逐道識別異常曲線形態(tài)更是不現(xiàn)實的工作。因此，鉆孔瞬變電磁法當(dāng)前的立體解釋還處于定性水平，水平分量異常曲線形態(tài)自動識別是阻礙鉆孔瞬變電磁法實現(xiàn)立體成像的關(guān)鍵因素，本文通過機器學(xué)習(xí)中的K-means聚類算法實現(xiàn)對大數(shù)據(jù)水平分量異常曲線形態(tài)的自動分類。

1 方法原理

1.1 垂直分量數(shù)據(jù)處理解釋方法

本文使用的鉆孔瞬變電磁法本質(zhì)上仍屬于中心回線裝置類型，其垂直(z)分量實測數(shù)據(jù)曲線形態(tài)與礦井瞬變電磁探測數(shù)據(jù)曲線形態(tài)基本相同(圖2)，僅是發(fā)射線圈尺寸與匝數(shù)不同導(dǎo)致電感影響較大，因此，數(shù)據(jù)處理方法可參考礦井瞬變電磁，采用范濤等提出的預(yù)處理技術(shù)[29]對電感影響進行校正，對校正后的數(shù)據(jù)可采用晚期視電阻率計算和層厚累加法進行處理解釋。

B為磁感應(yīng)強度圖2 礦井與鉆孔瞬變電磁探測數(shù)據(jù)曲線對比Fig.2 Comparison of Transient Electromagnetic DetectionData Curves Between Mine and Borehole

晚期視電阻率計算公式為

(1)

式中：ρ為晚期視電阻率；μ0為真空磁導(dǎo)率，取值為4π×10-7H·m-1；ST為發(fā)射線圈面積；SR為接收線圈面積；t為測道時間，V(t)/I為歸一化感應(yīng)電動勢，I為電流；C為視電阻率修正系數(shù)。

深度計算公式為

(2)

(3)

式中：Hi為第i道深度；K為深度修正系數(shù)；hi為第i道雙程傳播深度；ρi為第i道晚期視電阻率；ti為第i道采樣時間。

1.2 水平分量與異常體所在象限的規(guī)律

設(shè)計如圖3所示模型，規(guī)定水平x與y分量正方向之間區(qū)域為第一象限，逆時針依次定義為第二、三、四象限，在鉆孔深度方向50 m處，依次放置16個規(guī)模為20 m×20 m×6 m的板狀異常體，異常體中心點距離鉆孔15 m，模型其他參數(shù)見表1。繪制16個模型的水平分量異常響應(yīng)多測道圖(圖4～11)。

表1 模型參數(shù)

圖(a)中數(shù)字為異常體編號圖3 模型示意圖Fig.3 Schematic Views of Model

圖4～11分別為4個象限水平分量異常響應(yīng)。從圖4～11可以看出，以水平鉆孔為參考系，所有水平分量異常響應(yīng)均呈“S”或“反S”形態(tài)，且當(dāng)某一水平分量與該水平分量坐標(biāo)軸夾角為0°時，該水平分量的響應(yīng)幅值達到最小。兩組水平分量形態(tài)組合與異常體所在象限之間的關(guān)系如圖12所示。

1.3 異常體工具面角計算方法

鉆孔瞬變電磁法徑向探測時，異常體引起的二次場是矢量場。根據(jù)水平渦流場的空間分布特征，假設(shè)異常體為板狀體，在鉆孔中觀測到兩個磁場水平分量Vx、Vy的矢量和Vxy，其方向一定處于由鉆孔指向異常體的等效渦流中心上，那么只需要求出Vxy的方向，就知道異常體中心的具體方位[6]。

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

(4)

其中，Vx、Vy均為已知值，求反正弦即得到θ的表達式為

(5)

最后，根據(jù)異常體所在象限，可分別求出對應(yīng)的工具面角α：①異常體在第一象限，α=θ；②異常體在第二象限，α=π-θ；③異常體在第三象限，α=π+θ；④異常體在第四象限，α=2π-θ。

1.4 水平分量異常曲線類型自動識別方法

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

圖14是一組數(shù)值模擬數(shù)據(jù)的水平分量異常場特征值分布。從圖14可以看出其具有明顯的二分類特性，與水平分量異常場的“S”或“反S”形態(tài)存在顯著相關(guān)性。基于該特性，無需提前進行標(biāo)簽樣本的監(jiān)督訓(xùn)練，可直接選用無監(jiān)督機器學(xué)習(xí)算法對數(shù)據(jù)進行分類。

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

圖4 第一象限水平x分量異常響應(yīng)Fig.4 Horizontal x Component Anomaly Responses ofthe First Quadrant

圖5 第一象限水平y(tǒng)分量異常響應(yīng)Fig.5 Horizontal y Component Anomaly Responses ofthe First Quadrant

圖6 第二象限水平x分量異常響應(yīng)Fig.6 Horizontal x Component Anomaly Responses ofthe Second Quadrant

圖7 第二象限水平y(tǒng)分量異常響應(yīng)Fig.7 Horizontal y Component Anomaly Responses ofthe Second Quadrant

圖8 第三象限水平x分量異常響應(yīng)Fig.8 Horizontal x Component Anomaly Responses ofthe Third Quadrant

圖9 第三象限水平y(tǒng)分量異常響應(yīng)Fig.9 Horizontal y Component Anomaly Responses ofthe Third Quadrant

圖10 第四象限水平x分量異常響應(yīng)Fig.10 Horizontal x Component Anomaly Responses ofthe Fourth Quadrant

圖11 第四象限水平y(tǒng)分量異常響應(yīng)Fig.11 Horizontal y Component Anomaly Responses ofthe Fourth Quadrant

圖12 異常體位于不同象限時水平分量響應(yīng)形態(tài)Fig.12 Horizontal Component Response Forms ofAnomalous Body in Different Quadrants

圖13 異常偏轉(zhuǎn)角示意圖Fig.13 Schematic View of Anomalous Deflection Angle

圖14 特征值分布Fig.14 Distribution of Eigenvalues

聚類算法是指將一堆沒有標(biāo)簽的數(shù)據(jù)自動劃分成幾類的方法，這個方法要保證同一類數(shù)據(jù)有相似的特征。考慮到水平分量異常場特征值本身就有明顯的二分類特性，無需采用復(fù)雜的聚類算法，因此，本文選擇應(yīng)用廣泛、速度快、魯棒性強的K-means聚類算法。K-means聚類算法使用最大期望算法求解的高斯混合模型在正態(tài)分布的協(xié)方差為單位矩陣，且是在隱變量的后驗分布為一組狄拉克δ函數(shù)時所得到的特例，其假設(shè)相同類別中數(shù)據(jù)之間的距離應(yīng)該都很近，即數(shù)據(jù)之間的相似度與它們之間的歐式距離成反比。

本文需要將N個水平分量純異常數(shù)據(jù){xn}聚為2類，令經(jīng)過聚類之后每個數(shù)據(jù)所屬的類別為{tn}，而這2個聚類的中心為{μm}，則損失函數(shù)(L)計算公式為

(6)

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

(7)

將分類好的數(shù)據(jù)類別分別同“S”與“反S”形態(tài)進行對應(yīng)，根據(jù)圖12就可以準(zhǔn)確確定鉆孔瞬變電磁法觀測數(shù)據(jù)中任意測點任意測道反映的視電阻率信息所在的象限，再按照異常體工具面角計算方法，將每一測點每一測道對應(yīng)的視電阻率信息視為異常體帶入計算，即可獲取對應(yīng)的工具面角。

1.5 視電阻率立體成像算法

應(yīng)用1.1節(jié)方法可由鉆孔瞬變電磁法探測的垂直分量計算并繪制以鉆孔深度為橫坐標(biāo)、以鉆孔徑向探測距離為縱坐標(biāo)的視電阻率剖面成像圖。在已知圖中每一測點每一測道視電阻率對應(yīng)的工具面角信息時，假設(shè)鉆孔為直鉆孔，可通過三角函數(shù)關(guān)系將其投影在xOy平面上，將剖面圖中每一測點的一維視電阻率信息曲線轉(zhuǎn)換為二維視電阻率信息平面，多個測點的二維成像結(jié)果組合即可實現(xiàn)鉆孔瞬變電磁視電阻率立體成像。

單個測點視電阻率信息一維坐標(biāo)擴展至二維坐標(biāo)的計算公式為

(8)

(9)

式中：j是測點數(shù)；i是測道數(shù)；x是鉆孔瞬變電磁坐標(biāo)系中視電阻率信息對應(yīng)的x方向坐標(biāo)；y是鉆孔瞬變電磁坐標(biāo)系中視電阻率信息對應(yīng)的y方向坐標(biāo)；r是鉆孔探測半徑。

2 數(shù)值模擬

為驗證基于水平分量異常場特征聚類的鉆孔瞬變電磁立體成像方法探測效果，設(shè)計如圖15所示的三維模型，使用孫懷鳳等提出的三維時域有限差分正演程序[30]進行數(shù)值模擬。在鉆孔深度方向40 m、第一象限45°處放置1個規(guī)模為5 m×5 m×5 m

圖15 數(shù)值模型示意圖Fig.15 Schematic Views of Numerical Model

的低阻異常體，異常體中心點距離鉆孔15 m，模型其他參數(shù)如表1所示。

對26～54 m的測點垂直分量數(shù)據(jù)采用晚期視電阻率計算和層厚累加法進行視電阻率成像，視電阻率修正系數(shù)取1.58，深度修正系數(shù)取1.00，可以得到如圖16所示的沿鉆孔方向的視電阻率剖面。從圖16可以較為清晰地看到，在鉆孔深度40 m、鉆孔徑向15 m位置有較為明顯的低阻異常響應(yīng)，但從該成果圖中無法反映異常的空間方位。

圖16 視電阻率剖面Fig.16 Profile of Apparent Resistivity

采用1.3、1.4節(jié)方法對模型水平x、y分量數(shù)據(jù)進行處理，可得到如圖17所示的單測點視電阻率展開平面圖。從圖17可以看出，在孔深40 m的平面圖中，xOy平面第一象限有明顯的低阻異常響應(yīng)，與坐標(biāo)系原點O(鉆孔)之間的距離為15 m，形狀規(guī)模與模型參數(shù)基本一致，而在孔深30 m和50 m的平面圖中則沒有明顯低阻異常顯示。

圖17 單測點視電阻率展開平面圖Fig.17 Expanded Plan Views for Apparent Resistivity ofSingle Measuring Point

將每一測點的x、y坐標(biāo)與鉆孔孔深z坐標(biāo)組合，通過Voxler軟件進行立體成像如圖18所示。圖18中，由藍綠色表示的低阻異常體與模型設(shè)置參數(shù)一致，說明基于水平分量異常場的聚類算法對鉆孔徑向視電阻率立體成像有效，準(zhǔn)確性較高。

圖18 視電阻率立體成像Fig.18 Stereo Imaging of Apparent Resistivity

3 結(jié) 語

(1)以水平鉆孔鉆進方向為z軸正方向，以孔口所在平面右向為x軸正方向，下向為y軸正方向，鉆孔瞬變電磁法所有水平分量異常響應(yīng)均呈“S”或“反S”形態(tài)，且當(dāng)某一水平分量與該水平分量坐標(biāo)軸夾角為0°時，該水平分量的響應(yīng)幅值達到最小，因此，通過水平x、y分量異常形態(tài)組合可判定異常體所在象限。

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

(3)將測量數(shù)據(jù)每一測點每一測道對應(yīng)的視電阻率視為一個獨立的異常體，采用K-means聚類算法對相應(yīng)的兩組水平分量異常曲線中的特征值進行二分類，可以自動確定任意測點任意測道視電阻率的分布象限，再由水平分量異常場幅值計算出工具面角，即可結(jié)合垂直分量成像結(jié)果實現(xiàn)鉆孔徑向視電阻率的立體成像。

(4)采用三維時域有限差分?jǐn)?shù)值模擬數(shù)據(jù)對基于水平分量異常場的聚類算法進行了檢驗，成像結(jié)果與模型吻合度較高，說明該方法準(zhǔn)確、有效地提高了鉆孔瞬變電磁法的解釋水平，可以推廣應(yīng)用至煤礦實際生產(chǎn)中，為掘進工作面隱伏水害精準(zhǔn)超前探測提供技術(shù)保障。

謹(jǐn)以此文慶祝母校七十周年華誕，祝愿母校光輝歷程更輝煌，人才輩出代代強，桃李滿天揚四海，碩果累累振中華！2001年9月到2008年6月，我在母校度過了七年的求學(xué)時光。剛?cè)雽W(xué)的時候，羅馬廣場還是食堂，八號學(xué)生公寓還沒竣工，運動操場全是泥地，軍訓(xùn)時經(jīng)常一陣風(fēng)刮過就卷起一片沙塵；之后短短幾年時間，母校就迎來了高速發(fā)展，渭水校區(qū)基本建成，交通科技大廈、地學(xué)科技大廈也都拔地而起，呈現(xiàn)一片欣欣向榮的景象?；叵肫鹪?jīng)那個“土土”的校園和無憂無慮、一起學(xué)習(xí)進步的同學(xué)，總讓人想穿越回去再次感受那份善良與純真、真誠與熱情！