基于聚類的煤礦井下鉆孔瞬變電磁異常響應(yīng)邊界成像方法

范 濤，李 萍，張幼振，趙 睿，房 哲，樊依林，劉 磊，王 程，李宇騰

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

近10 年煤礦事故總量逐年下降，但事故下降幅度趨緩。較大事故總起數(shù)、死亡人數(shù)總體呈下降趨勢，但從2013 年開始曲線斜率明顯趨平，2016 年后呈鋸齒形下降，2017 和2019 年等一些年份反彈，2020 年再度下降。重特大事故呈鋸齒形下降，同時波動幅度較大，2018 年后出現(xiàn)反彈，2020 年重特大事故總量與2019 年持平。事故中水害較大以上事故上升幅度大，2020 年發(fā)生較大以上水害事故3 起、死亡21 人，同比增加1 起、多死亡12 人，分別上升50%和133.3%，其中，較大水害事故2 起、死亡8 人，同比起數(shù)持平、少死亡1 人，重大事故1 起、死亡13 人，同比增加1 起、多死亡13 人，占全國煤礦較大以上事故起數(shù)和死亡人數(shù)的37.5%，經(jīng)濟(jì)損失和社會影響非常嚴(yán)重[1]。而掘進(jìn)工作面則是煤礦水害事故發(fā)生最多的地點[2]。因此，掘進(jìn)工作面前方隱伏水害超前探測是亟待解決的技術(shù)難題。

水害隱患超前探查代表性的方法主要是瞬變電磁法、直流電法等電磁類方法[3-5]，但是，煤礦井下探測裝備功率受煤礦本安防爆限制，且井下環(huán)境的電磁干擾較大等因素的存在，致使礦井電磁類方法的探測距離短、精度偏低、多解性強(qiáng)[6-8]。為了解決煤礦井下掘進(jìn)工作面前方的探測深度與探測精度的矛盾，許多學(xué)者逐步開始利用井下掘進(jìn)工作面的鉆孔進(jìn)行瞬變電磁探測工作，該方法可以在掘進(jìn)前開展遠(yuǎn)距離、高精度的隱伏水害超前預(yù)報[9-12]。

鉆孔瞬變電磁的數(shù)據(jù)處理一般采用電阻率反演成像方法，考慮到反演存在多解性，因此，反演擬合過程中為了提高精度，往往會選擇較為光滑的模型擬合實測數(shù)據(jù)，使得最終反演結(jié)果是一個電阻率連續(xù)光滑漸變的成像模型，在地質(zhì)體邊界處的對比程度和變化情況較為模糊，很難清晰地反映地質(zhì)異常體響應(yīng)與背景值的差異，對異常體規(guī)模、形態(tài)的解釋工作常常需要經(jīng)驗豐富的專家人為干預(yù)，也不利于生產(chǎn)中準(zhǔn)確指導(dǎo)物探工作之后的鉆探和掘進(jìn)工作。鑒于此，有必要研究一種提高成像結(jié)果中電阻率值聚合度，進(jìn)而突出電性邊界的成像方法。

顯然，對電阻率的分區(qū)聚合屬于典型的分類問題，可以選用無監(jiān)督機(jī)器學(xué)習(xí)方法進(jìn)行解決。數(shù)據(jù)挖掘領(lǐng)域經(jīng)常使用無監(jiān)督機(jī)器學(xué)習(xí)算法，主要是采用它來發(fā)現(xiàn)大量無標(biāo)簽數(shù)據(jù)的分布規(guī)律，實現(xiàn)對數(shù)據(jù)的區(qū)分或分類。經(jīng)常使用的無監(jiān)督機(jī)器學(xué)習(xí)方法主要有局部線性嵌入算法、主成分分析、局部切空間排列算法、拉普拉斯特征映射、等距映射和應(yīng)用最多最廣的聚類算法[13]。

聚類算法是指基于一定的優(yōu)化標(biāo)準(zhǔn)將一堆無標(biāo)簽數(shù)據(jù)對象自動劃分成若干類別的方法，這個方法要求同一類別的數(shù)據(jù)具有相似的特征，不同類別的數(shù)據(jù)具有不同的特征[14]。

近年來，聚類方法已經(jīng)廣泛被應(yīng)用于地球物理領(lǐng)域的數(shù)據(jù)處理和特征挖掘中。尤其在地震勘探領(lǐng)域的應(yīng)用較多，王偉濤等[15]通過層次聚類分析對汶川大地震的余震序列中的近似地震和重復(fù)地震進(jìn)行了有效辨識；張巖等[16]采用結(jié)構(gòu)聚類字典學(xué)習(xí)方法進(jìn)行了地震數(shù)據(jù)隨機(jī)噪聲壓制方面的研究；S.Scitovski[17]采用基于密度的聚類方法對地震記錄進(jìn)行了不同地震類型的劃分。在重磁資料的處理解釋中，張新兵等[18]基于改進(jìn)的K-Means 聚類分析方法實現(xiàn)了重磁局部異常的自動圈定；李斐等[19]在優(yōu)化不同區(qū)域的重力觀測密度方面應(yīng)用了聚類分析方法；曹書錦等[20]開展了自適應(yīng)模糊聚類對多異常源的精準(zhǔn)確定工作。在電磁數(shù)據(jù)處理解釋領(lǐng)域，楊生等[21]在大地電磁曲線類型劃分中使用了聚類方法，對削弱地質(zhì)推斷的多解性有較好的作用；SongYuchen 等[22]提出一種應(yīng)用自適應(yīng)密度聚類分類電測深曲線類型的方法；M.Audebert 等[23]結(jié)合了K-Means 聚類和電阻率CT 成像，提出一種多次反演解釋垃圾場滲流區(qū)的有效方法；李晉等[24]提出一種聯(lián)合聚類及遞歸的信噪比分離方法，對MT 低頻部分的數(shù)據(jù)質(zhì)量有所改善。

本文參考以上資料，考慮不同的地質(zhì)體存在電阻率差異，基于統(tǒng)計思想提出將具有相近電阻率值的成像元素劃分至一個類別，進(jìn)而實現(xiàn)快速識別地質(zhì)異常響應(yīng)邊界的方法。本文討論了選用聚類方法的原則，介紹了聚類的方法原理，研究了確定最佳聚類條件的方法，最后使用三維數(shù)值模擬和井下實測數(shù)據(jù)檢驗了方法的有效性和實用性。

1 方法原理

1.1 聚類方法的選取

對電阻率進(jìn)行聚類本質(zhì)上是對一個一維數(shù)據(jù)體進(jìn)行聚類，相當(dāng)于對電阻率進(jìn)行自適應(yīng)層級劃分。絕大多數(shù)聚類算法都是針對二維及以上的數(shù)據(jù)，針對一維數(shù)據(jù)的聚類算法非常少，主要有K-Means 方法和Jenks Natural Breaks 方法。

K-Means 算法，也被稱為K-平均或K-均值算法，是一種廣泛使用的聚類算法。K-Means 算法以距離作為數(shù)據(jù)對象間相似性度量的標(biāo)準(zhǔn)，即數(shù)據(jù)對象間的距離越小，則它們的相似性越高，則它們越有可能在同一個類簇。之所以被稱為K-Means 是因為它可以發(fā)現(xiàn)K個不同的簇，且每個簇的中心采用簇中所含值的均值計算而成。

Jenks Natural Breaks 算法，也就是自然斷點分類，分類的原則就是將方差接近的放在一起，分成若干類。通過計算每類的方差和方差和，用方差和的大小來比較分類的好壞，值最小的就是最優(yōu)的分類結(jié)果(但并不唯一)。

K-Means 和Jenks Natural Breaks 在處理一維數(shù)據(jù)時完全等價。它們的目標(biāo)函數(shù)一樣，但是算法的步驟不完全相同。K-Means 是先設(shè)定好K個初始隨機(jī)點。而Jenks Natural Breaks 則是用遍歷的方法，一個點一個點地移動，直到達(dá)到最小值。顯然，n個數(shù)分成k類，Jenks Natural Breaks 就要從n-1 個數(shù)中找k-1 個組合。當(dāng)數(shù)據(jù)量較大時，如果分類又多，計算量會顯著增加，因此，一般針對一維數(shù)組的聚類，更多選用計算量小的K-Means 算法。

1.2 K-Means 聚類算法

K-Means 算法的思想很簡單，對于給定的樣本集，按照樣本之間的距離大小，將樣本集劃分為K個簇。讓簇內(nèi)的點盡量緊密地連在一起，而讓簇間的距離盡量大。

假設(shè)數(shù)據(jù)集為：

式中：n為數(shù)據(jù)個數(shù)；xi為除質(zhì)心外的其他點。

劃分的簇為：

式中：K為預(yù)先設(shè)置的分類數(shù)。

此時，最小化平方誤差E即為：

式中：μi是簇Ci的均值向量，也稱為質(zhì)心，表達(dá)式為：

計算開始時，先采用隨機(jī)方式選定K個質(zhì)心對所有數(shù)據(jù)進(jìn)行初始分類，共分為K個初始類別，之后對所有數(shù)據(jù)計算它們與質(zhì)心之間的歐氏距離，再依據(jù)這些距離的平均值更新質(zhì)心和劃分類別，最后反復(fù)迭代該更新過程，直到滿足迭代的停止條件為止。

迭代的停止條件一般是質(zhì)心變化率滿足下式：

1.3 聚類中關(guān)鍵條件的確定方法

K-Means 聚類算法需要提前給出的條件主要有3 個：初始質(zhì)心、聚類時的距離計算規(guī)則和聚類數(shù)目。

由于K-Means 聚類易陷入局部極小，不同的初始質(zhì)心可能會導(dǎo)致不同的結(jié)果，本文根據(jù)聚類數(shù)目K的不同，選擇距離盡可能遠(yuǎn)的K個點為初始質(zhì)心，具體做法為：隨機(jī)選擇一個點作為第一個初始簇質(zhì)心，然后選擇距離該點最遠(yuǎn)的那個點作為第二個初始簇質(zhì)心，然后再選擇距離前2 個點的最近距離最大的點作為第3 個初始簇的質(zhì)心，以此類推，直至選出K個初始類質(zhì)心。

考核類簇之間的相似性程度主要依靠各類簇之間的距離，較為常見的距離計算方法有歐氏距離、曼哈頓距離、夾角余弦距離和相關(guān)距離等。一般來說，歐氏距離最為簡單直觀，也更能反映數(shù)據(jù)在數(shù)值特征上的差別，因此，本文選擇使用歐氏距離。

聚類數(shù)目是K-Means 聚類算法中最重要的參數(shù)，因為不準(zhǔn)確的聚類數(shù)目會明顯導(dǎo)致分類效果變差，但是目前并沒有依托數(shù)學(xué)原理的完美評價標(biāo)準(zhǔn)，本文采用基于組內(nèi)平方誤差和(Sum of Squared Error，SSE)的肘部法則來確定最佳聚類數(shù)目。

組內(nèi)平方誤差和(SSE)表示一個類簇內(nèi)各點與該類質(zhì)心的平方誤差之和，可由下式計算。

ESS越小則說明各個類簇越收斂，但顯然不是越小越好。因為考慮一種極端情況：將所有的樣本點均視作單獨類簇，此時ESS為0，而并未達(dá)到分類的目的。因此，需要在聚類數(shù)目和ESS之間尋找一個平衡點。

肘部法則就是這樣一種方法。指定一個j值，即可能的最大聚類數(shù)目。然后將聚類數(shù)目從1 開始一直遞增到j(luò)，計算出j個ESS。隨著聚類數(shù)目增多，每一個類簇中數(shù)據(jù)點數(shù)量越來越少，距離越來越近，因此，ESS值必然隨著聚類數(shù)目增多而減少。但當(dāng)ESS減少得很緩慢時，可以認(rèn)為進(jìn)一步增大聚類數(shù)目分類效果也并不能增強(qiáng)，這個“肘點(拐點)”就是最佳聚類數(shù)目。

1.4 異常響應(yīng)邊界成像步驟

應(yīng)用K-Means 聚類方法實現(xiàn)對異常響應(yīng)邊界成像的具體步驟如下：

(1) 對實測數(shù)據(jù)的Z分量數(shù)據(jù)進(jìn)行視電阻率計算或反演成像，獲得(視)電阻率參數(shù)。

(2) 計算(視)電阻率參數(shù)不同聚類數(shù)目情況下的組內(nèi)平方誤差和(ESS)，應(yīng)用肘部法則尋找最佳聚類數(shù)目。

(3) 按照最佳聚類數(shù)目對(視)電阻率參數(shù)進(jìn)行聚類。

(4) 將同一類中的(視)電阻率值全部更新為該類的質(zhì)心值，并與原空間坐標(biāo)重新對應(yīng)。

(5) 對新的聚類后電性參數(shù)數(shù)據(jù)文件進(jìn)行成像。

2 數(shù)值模擬效果

為驗證本文方法的成像效果，設(shè)計如圖1 所示的三維模型，采用時域有限差分方法進(jìn)行了數(shù)值模擬。在鉆孔深度方向40 m 處，第一象限45°放置1 個規(guī)模為10 m×10 m×10 m 的低阻異常體，異常體中心點距離鉆孔15 m，采樣時間范圍為1×10-6～6.136×10-4s，模型其他參數(shù)見表1。

圖1 數(shù)值模型Fig.1 Schematic diagram of the numerical model

表1 模型參數(shù)Table 1 Model parameters

對26～54 m 范圍內(nèi)的測點Z分量數(shù)據(jù)采用晚期視電阻率計算和層厚累加法進(jìn)行視電阻率成像，可以得到如圖2 所示的沿鉆孔方向的視電阻率剖面圖，圖中橫坐標(biāo)z為鉆孔深度，縱坐標(biāo)r為徑向探測距離?？梢暂^為清晰地看到在鉆孔深度40 m、鉆孔徑向15 m 位置有較為明顯的低阻異常響應(yīng)(藍(lán)、綠色區(qū)域)。

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

如果將藍(lán)色區(qū)域解釋為異常，則異常響應(yīng)明顯小于實際異常體，二者面積比值約為0.624；如果將綠色區(qū)域解釋為異常，則異常響應(yīng)明顯大于實際異常體，二者面積比值約為2.458；從圖2 中無法清晰反映異常響應(yīng)的邊界，需要根據(jù)經(jīng)驗人為解釋確定。

采用本文1.3 節(jié)肘部法則對模型視電阻率數(shù)據(jù)進(jìn)行分析，可得到如圖3 所示的ESS曲線圖。由圖3 可知，在聚類數(shù)目為2 時曲線出現(xiàn)肘點，因此，選擇K=2，采用1.2 節(jié)中的K-Means 算法對模型視電阻率數(shù)據(jù)進(jìn)行聚類處理，可得到圖4 所示的異常響應(yīng)邊界成像結(jié)果。由圖4 可知，在鉆孔深度40 m、鉆孔徑向15 m 位置有明顯的孤立低阻異常響應(yīng)(墨藍(lán)色區(qū)域)，異常響應(yīng)邊界非常清晰，受晚期視電阻率計算和圖像插值算法的影響，異常響應(yīng)形狀表現(xiàn)為直徑約10 m 的近似圓形，其與實際異常體(紅色虛線)的面積比約為0.985，整體參數(shù)與模型設(shè)置參數(shù)吻合較好。

圖3 數(shù)值模型ESS 曲線Fig.3 ESS curve of the numerical model

圖4 數(shù)值模型異常響應(yīng)邊界成像結(jié)果Fig.4 Imaging result of anomaly response boundary of the numerical model

3 應(yīng)用實例

為了驗證本文異常響應(yīng)邊界成像方法對實際生產(chǎn)中煤礦積水采空區(qū)的實際解釋能力，在陜北某礦開展了鉆孔瞬變電磁探測試驗。該礦井屬于侏羅紀(jì)煤田，之前由于大量使用以掘代采的方式采煤，產(chǎn)生較多的小型采空區(qū)，且資料保留不完全，位置不明。礦區(qū)煤層上部有砂巖裂隙水，部分區(qū)段還有第四系松散層潛水，因此，這些采空區(qū)大部分為積水采空區(qū)，對該煤礦的安全生產(chǎn)造成了重要影響。

鉆孔瞬變電磁探測測點范圍為孔深60～128 m。對測點數(shù)據(jù)進(jìn)行電阻率反演成像可以得到如圖5 所示的沿鉆孔方向的電阻率剖面圖。由圖5 可知，在鉆孔深度100 m、鉆孔徑向20 m 位置有較為明顯的低阻異常響應(yīng)(藍(lán)、綠色區(qū)域)，與數(shù)值模擬結(jié)果類似，從該成果圖中無法清晰反映異常響應(yīng)的邊界，需要根據(jù)經(jīng)驗人為解釋確定。

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

同樣，采用本文1.3 節(jié)肘部法則對電阻率數(shù)據(jù)進(jìn)行分析，可得到如圖6 所示的ESS曲線圖，由圖6 可知，在聚類數(shù)目為3 時曲線出現(xiàn)肘點，因此，選擇K=3。采用1.2 節(jié)中的K-Means 算法對電阻率數(shù)據(jù)進(jìn)行聚類處理，可得到如圖7 所示的異常響應(yīng)邊界成像結(jié)果。由圖7 可知，在鉆孔深度100 m、鉆孔徑向20 m 位置有明顯的孤立低阻異常響應(yīng)(藍(lán)色區(qū)域)，異常響應(yīng)邊界非常清晰，異常響應(yīng)形狀表現(xiàn)為邊長約20 m 的三角形。礦方后期設(shè)計了上仰鉆孔對該異常進(jìn)行勘查驗證，在進(jìn)尺98 m 附近出現(xiàn)掉鉆和出水現(xiàn)象，最終推斷該異常為上組煤遺留的小煤窯采空區(qū)。

圖6 實測數(shù)據(jù)ESS 曲線Fig.6 ESS curve of the measured data

圖7 實測數(shù)據(jù)異常響應(yīng)邊界成像結(jié)果Fig.7 Imaging result of anomaly response boundary of the measured data

驗證鉆孔開孔角度為上仰8°，此時鉆孔到達(dá)異常響應(yīng)邊界的距離為96.89 m，與實際揭露采空區(qū)的距離相差1.11 m，誤差為1.13%。

4 結(jié) 論

a.鉆孔瞬變電磁反演結(jié)果一般為連續(xù)光滑漸變的電阻率成像模型，對背景值和異常地質(zhì)體分界面的對比度表現(xiàn)得較為模糊，難以準(zhǔn)確解釋異常響應(yīng)的邊界。采用無監(jiān)督機(jī)器學(xué)習(xí)中的聚類方法可以通過將相近電阻率值進(jìn)行聚合分類實現(xiàn)異常響應(yīng)邊界快速成像。

b.針對電阻率這一個特征值，從計算量考慮選擇K-Means 聚類算法，按照電阻率樣本之間的距離大小，對樣本進(jìn)行分類。聚類計算中，基于最遠(yuǎn)距離原則確定初始質(zhì)心，距離計算方法選用歐氏距離方法，聚類數(shù)目應(yīng)用基于組內(nèi)平方誤差和(ESS)的肘部法則進(jìn)行確定。

c.數(shù)值模擬數(shù)據(jù)的處理效果說明這種基于聚類的成像方法有效突出了異常響應(yīng)邊界，清晰顯示了異常響應(yīng)形狀；煤礦井下應(yīng)用實例解釋的積水采空區(qū)經(jīng)過實際鉆探驗證，說明了該方法在現(xiàn)實探查中的準(zhǔn)確性和有效性，提高了巷道超前探測技術(shù)的精度和可靠度，為超前探測工程技術(shù)難題提供了重要的技術(shù)保障。

d.本文方法本質(zhì)上是一種統(tǒng)計方法，并不僅僅適用于瞬變電磁法，也不僅僅適用于電阻率成像結(jié)果，對于以熱力圖形式表現(xiàn)的成像數(shù)據(jù)，如音頻電透視及無線電波透視成像結(jié)果、直流電法探測成果，也均可應(yīng)用。