基于聚類算法的鉆孔瞬變電磁視電阻率立體成像方法

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

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

0 引 言

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

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

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

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

從前人研究成果可以看出，水平分量形態(tài)組合和幅值差異對孔旁地質(zhì)異常體位置敏感，結合垂直分量視電阻率成像結果可對異常體進行立體解釋。但是，根據(jù)水平分量異常形態(tài)組合確定異常體所在象限需人工進行識別和判斷，效率較低，尤其當測點較多時，人工逐點逐道識別異常曲線形態(tài)更是不現(xiàn)實的工作。因此，鉆孔瞬變電磁法當前的立體解釋還處于定性水平，水平分量異常曲線形態(tài)自動識別是阻礙鉆孔瞬變電磁法實現(xiàn)立體成像的關鍵因素，本文通過機器學習中的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ù)不同導致電感影響較大，因此，數(shù)據(jù)處理方法可參考礦井瞬變電磁，采用范濤等提出的預處理技術[29]對電感影響進行校正，對校正后的數(shù)據(jù)可采用晚期視電阻率計算和層厚累加法進行處理解釋。

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

晚期視電阻率計算公式為

(1)

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

深度計算公式為

(2)

(3)

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

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

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

表1 模型參數(shù)

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

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

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

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

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

(4)

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

(5)

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

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

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

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

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

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

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

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

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

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

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

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

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

圖12 異常體位于不同象限時水平分量響應形態(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

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

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

(6)

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

(7)

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

1.5 視電阻率立體成像算法

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

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

(8)

(9)

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

2 數(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位置有較為明顯的低阻異常響應，但從該成果圖中無法反映異常的空間方位。

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

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

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

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

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

3 結 語

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

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

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

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

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