基于K-means法的微震事件活動(dòng)聚類分析

周立云 馮興隆 林潔

摘? 要：以自然崩落法開(kāi)采為背景，通過(guò)微震監(jiān)測(cè)系統(tǒng)對(duì)開(kāi)采區(qū)域的微震事件進(jìn)行監(jiān)測(cè)。利用K-means法對(duì)微震事件活動(dòng)進(jìn)行聚類分析，結(jié)合現(xiàn)場(chǎng)實(shí)際情況，得出3 736 m拉底水平和3 720 m出礦水平是微震事件最為集中的區(qū)域，采礦活動(dòng)是該區(qū)域微震事件產(chǎn)生的主要原因。在不同分類數(shù)條件下，微震事件的集中區(qū)域是一致的，具體應(yīng)根據(jù)現(xiàn)場(chǎng)實(shí)際選擇合適的分類數(shù)。結(jié)果表明，基于K-means法的微震事件聚類分析能夠較清晰的反映出微震活動(dòng)規(guī)律，可提前預(yù)判大事件發(fā)生范圍，能為礦山的安全生產(chǎn)提供依據(jù)，具一定借鑒意義。

關(guān)鍵詞：自然崩落法;微震監(jiān)測(cè);聚類分析;K-means法

礦山開(kāi)采過(guò)程中，自然崩落法是唯一可與露天開(kāi)采相媲美的采礦方法，它成本低、高效率，尤其針對(duì)一些礦體厚大、品位較低和礦巖完整性差的礦山[1]。崩落過(guò)程中，崩落面和崩落下來(lái)的礦石堆之間留有較大空隙，如果放礦速度過(guò)快，上部礦巖突然發(fā)生大范圍崩落，極有可能產(chǎn)生空氣沖擊波，對(duì)坑內(nèi)人員、設(shè)備產(chǎn)生巨大影響。澳大利亞北帕克斯礦曾發(fā)生過(guò)該類地災(zāi)[2]，因未準(zhǔn)確監(jiān)測(cè)崩落頂板的位置，在底部持續(xù)大量出礦使崩落面和存隆面之間距離過(guò)大，造成大范圍的頂板突然垮塌和嚴(yán)重的沖擊波災(zāi)害。自然崩落法開(kāi)采過(guò)程中，礦山井下絕大多數(shù)的人員、設(shè)備均在底部結(jié)構(gòu)內(nèi)進(jìn)行生產(chǎn)活動(dòng)，底部結(jié)構(gòu)穩(wěn)定性也是決定礦山生產(chǎn)成敗的關(guān)鍵因素[3]。因此，必須通過(guò)相關(guān)監(jiān)測(cè)手段，監(jiān)測(cè)崩落的發(fā)生和發(fā)展，監(jiān)測(cè)出礦水平設(shè)施和其他基礎(chǔ)設(shè)施的穩(wěn)定性，保障安全、降低風(fēng)險(xiǎn)。

微震監(jiān)測(cè)作為一種三維可視化空間監(jiān)測(cè)技術(shù)，對(duì)微震事件發(fā)生的時(shí)間、空間、強(qiáng)度信息反映了巖體微破裂的發(fā)展。因此，許多研究者試圖找到微震活動(dòng)與開(kāi)采擾動(dòng)之間的關(guān)系，并為微震預(yù)測(cè)提供定量分析的基礎(chǔ)。相較于其它方法，通過(guò)觀察微震事件發(fā)展的時(shí)空分布規(guī)律是最有效直觀的方式之一。除通過(guò)直觀觀察微震活動(dòng)的時(shí)空分布探尋巖石損傷的規(guī)律外，許多學(xué)者還在微震事件的聚集度和活動(dòng)特征方面做了更多的研究。吳愛(ài)祥等利用最短距離聚類法對(duì)某礦山微地震活動(dòng)的時(shí)空分布進(jìn)行分析[4]，識(shí)別了井下微震活動(dòng)聚集區(qū)域;Martin利用空間聚類方法對(duì)微震事件進(jìn)行聚類分析并解釋了微震活動(dòng)規(guī)律[5];程愛(ài)平等利用未確知聚類優(yōu)化法[6]，以樣本均值為聚類中心判別權(quán)重，建立動(dòng)態(tài)預(yù)測(cè)模型;劉棟、李夕兵等提出時(shí)空共享近鄰聚類算法（STSNN），對(duì)沙壩礦微震事件進(jìn)行了聚類分析[7]，揭示了微震事件活動(dòng)規(guī)律的研究可為大事件的產(chǎn)生提供有效的預(yù)判信息，為保證礦山安全生產(chǎn)發(fā)揮重要作用。

1? K-means算法

聚類分析是通過(guò)定量的方法，按照事物的相似性屬性，可以是類與類之間的距離、相關(guān)系數(shù)等[8]。K-means算法采用距離作為相似性指標(biāo)，這符合以微震事件點(diǎn)的坐標(biāo)位置劃分疏密關(guān)系的聚類分析為目的[9]。首先，設(shè)立分類目標(biāo)，確定將整個(gè)數(shù)據(jù)集分成K類，每一類有一個(gè)聚集中心點(diǎn)，聚集中心點(diǎn)由聚類密集程度劃分，一般采用歐式距離作為相似度指標(biāo)，聚類的目標(biāo)就是最小化，即平方和達(dá)到最小[10]：

式中：J為距離平方和;k為分類數(shù); n為n個(gè)數(shù)據(jù)點(diǎn); xi為每個(gè)點(diǎn); uk為聚類中心。

K-means算法中涉及了最小二乘法和拉格朗日原理，以聚類中心的反復(fù)迭代為目的，找到密集度高的分類區(qū)域。迭代步驟：① 系統(tǒng)通過(guò)空間點(diǎn)的位置坐標(biāo)分析疏密程度，按照目標(biāo)分類，選取k個(gè)聚類中心作為初始中心。②按照歐式距離的計(jì)算標(biāo)準(zhǔn)計(jì)算空間位置，反復(fù)迭代，找出聚類中心附近最近的點(diǎn)，予以劃分。③ 更新聚類初始中心，對(duì)每個(gè)分類中的對(duì)象進(jìn)行平均值的驗(yàn)算;分別計(jì)算作為聚類中心，計(jì)算目標(biāo)函數(shù)的值。④ 判斷聚類中心和目標(biāo)函數(shù)值是否發(fā)生改變，若不變，則輸出結(jié)果，若改變，則返回步驟②。給定N個(gè)三維坐標(biāo)點(diǎn)，指定分為k類，從N個(gè)三維坐標(biāo)點(diǎn)覆蓋的三維體積中挑出k個(gè)點(diǎn)，稱為種子元素，然后遍及所有元素點(diǎn)，不斷按照歐式距離計(jì)算與種子元素相近的點(diǎn)，更新聚類中心，計(jì)算聚類對(duì)象平均值，通過(guò)反復(fù)迭代的過(guò)程，找出最適宜的分類[11]。

2? K-means法微震活動(dòng)聚類分析

探尋微震監(jiān)測(cè)的時(shí)空分布規(guī)律是預(yù)測(cè)預(yù)警中最為有效的方法，通常微震活動(dòng)呈一種較分散的狀態(tài)，人為劃分區(qū)域有一定主觀性。同時(shí)，對(duì)于微震活動(dòng)的時(shí)間響應(yīng)也沒(méi)有很好的體現(xiàn)。本文運(yùn)用聚類分析來(lái)定量分析微震活動(dòng)時(shí)空分布，并初步劃出聚類區(qū)域。對(duì)微震活動(dòng)密集的聚類區(qū)域進(jìn)行時(shí)間序列分析，據(jù)時(shí)空分布特征，對(duì)微震事件所在區(qū)域進(jìn)行穩(wěn)定性評(píng)價(jià)[12]。

對(duì)首采區(qū)監(jiān)測(cè)范圍內(nèi)微震監(jiān)測(cè)點(diǎn)的三維坐標(biāo)進(jìn)行聚類分析。微震事件的發(fā)生具有離散性，人為劃分區(qū)域有一定局限性。此次分析使用三維K-means聚類分析，找出聚類中心，將點(diǎn)坐標(biāo)之間的親疏距離相近的分為一類定量的劃分區(qū)域[13]。結(jié)合開(kāi)采區(qū)的工作狀況，可對(duì)危險(xiǎn)區(qū)域進(jìn)行劃分，再對(duì)劃分的重點(diǎn)區(qū)域進(jìn)行時(shí)空演化規(guī)律研究[14]。

通過(guò)收集某銅礦2017年9月至2018年3月間發(fā)生的689次微震事件的監(jiān)測(cè)數(shù)據(jù)做聚類分析。聚類前將689次微震事件發(fā)生的時(shí)間、空間點(diǎn)三維坐標(biāo)、震級(jí)等參數(shù)統(tǒng)計(jì)，畫出空間位置散點(diǎn)圖，其中圓球即代表一個(gè)微震事件（圖1）[15]。

圖1中的圓球即代表微震事件，y為礦體走向，x為礦體傾向，從北至南分為N4～S9的13條穿脈，3 720m出礦水平和3 736 m拉底水平的穿脈之間的聚礦槽和掘進(jìn)工作面有相應(yīng)的爆破作業(yè)，可以看出3 720m出礦水平和3 736 m拉底水平附近微震事件較多，有向四周及頂板向上擴(kuò)展的趨勢(shì)，頂板最高的位置已到達(dá)3 910 m。

運(yùn)用Matlab軟件相關(guān)程序按照K-means方法對(duì)其空間分布進(jìn)行聚類分析，按照分為3類和5類的2種情況得出以下結(jié)果（圖2，圖3），其中黑色方點(diǎn)為每一類的一個(gè)聚類中心點(diǎn)，圖4和圖5是為便于分析的切面圖，同時(shí)分別對(duì)2種分類結(jié)果進(jìn)行統(tǒng)計(jì)（表1，2）。對(duì)各聚類中微震事件的時(shí)間分布進(jìn)行分析（圖6，7）。

以上是通過(guò)K-means劃分為3類的情況，并分別從水平和豎直方向進(jìn)行切面分析，借助水平圖和側(cè)視圖，判斷聚類中心點(diǎn)的位置，從而劃分危險(xiǎn)區(qū)域，探討聚類中心產(chǎn)生原因，分析微震活動(dòng)性和生產(chǎn)活動(dòng)的關(guān)系。從圖2和圖4中可以看出，微震活動(dòng)主要分布在沿礦體走向3.102 8×106 m至3.103 3×106 m的位置，傾向1.597 1×107 m至1.597 6×107 m的位置，垂直高度為3 550 m至3 910 m的范圍內(nèi)，與礦山開(kāi)采活動(dòng)的范圍較為一致，說(shuō)明3 736 m拉底出礦水平和3 720 m出礦水平聚礦槽爆破是產(chǎn)生微震事件聚集的主要原因。

圖6中，橫坐標(biāo)為時(shí)間列，第一點(diǎn)表示的是2017年8月9日，最后一個(gè)點(diǎn)是2018年3月9日，縱坐標(biāo)表示了聚類1、聚類2、聚類3在時(shí)間序列中的發(fā)展變化[16]。通過(guò)各聚類中的微震事件時(shí)間分布可看出，微震事件的聚集主要集中在2017年10月、11月、2018年3月，由出礦量統(tǒng)計(jì)可知，微震事件與大量出礦有關(guān)。

3 720 m水平出礦量要遠(yuǎn)大于3 736 m水平，2017年9月出礦量較大，2018年3月3 720 m水平出礦量達(dá)到歷史新高（圖8），大量出礦就要重點(diǎn)監(jiān)測(cè)底部結(jié)構(gòu)穩(wěn)定性、頂板崩落等情況，發(fā)現(xiàn)問(wèn)題需及時(shí)調(diào)整放礦計(jì)劃。同時(shí)出礦量代表了采礦活動(dòng)性，出礦量的統(tǒng)計(jì)有利于結(jié)合微震監(jiān)測(cè)數(shù)據(jù)分析采礦活動(dòng)與微震事件活動(dòng)的關(guān)系。

通過(guò)聚類圖導(dǎo)出聚類中心三維點(diǎn)坐標(biāo)可看出（表1），聚類1有269個(gè)事件，聚類2有217個(gè)事件，聚類3有211個(gè)事件。各聚類之間界限分明，距離聚類中心較遠(yuǎn)的點(diǎn)可忽略不計(jì)，帶有偶然性因素，可能是誤差或不穩(wěn)定影響因素。找到3個(gè)聚類中心點(diǎn)，分別是聚類中心1、聚類中心2、聚類中心3，聚類中心點(diǎn)圍繞在3 736 m水平和3 720 m水平上下，是受到拉底出礦水平聚礦槽爆破等影響。而z軸豎直方向可以看到3 751 m、3 687 m、3 672 m三個(gè)聚類點(diǎn)，位于3 736 m拉底水平上較近的聚類1，位于3 720 m出礦水平下的聚類2和聚類3。除了爆破和出礦的因素，還有應(yīng)力重分布的影響，有向四周擴(kuò)散的趨勢(shì)。以下為k=5的聚類結(jié)果。

以上是k=5的聚類結(jié)果，k=5類的結(jié)果與k=3類的結(jié)果相似，在3 720 m水平和3 736 m水平附近的微震事件活躍度較高，微震活動(dòng)主要分布在沿礦體走向3.102 8×106 m至3.103 3×106 m的位置，傾向1.597 1×107 m至1.59 76×107 m的位置，垂直高度為3 550 m到3 910 m的范圍內(nèi)，與開(kāi)采活動(dòng)范圍依然保持一致。

圖7表示了聚類1、聚類2、聚類3、聚類4、聚類5微震事件的時(shí)間分布，微震事件主要集中在2017年11月、2018年3月間，這與大量出礦有關(guān)系，具體聚類結(jié)果統(tǒng)計(jì)如表2所示。

將聚類3類分為更細(xì)的情況，統(tǒng)計(jì)分類數(shù)為5的情況，可以看出聚類1為164個(gè)事件，聚類2為145個(gè)事件，聚類3為117個(gè)事件，聚類4為119個(gè)事件，聚類5為152個(gè)事件，各聚類中心統(tǒng)計(jì)可看出（表2），聚類中心點(diǎn)在3 736 m拉底水平和3 720 m出礦水平上下，這多受首采區(qū)監(jiān)測(cè)區(qū)開(kāi)采活動(dòng)影響，3 736 m拉底水平上有3個(gè)聚類區(qū)域，372 0m出礦水平下有兩個(gè)聚類區(qū)域。

3? 結(jié)論

本文通過(guò)某銅礦自然崩落法現(xiàn)場(chǎng)開(kāi)采情況，借助微震監(jiān)測(cè)系統(tǒng)對(duì)礦山開(kāi)采過(guò)程進(jìn)行監(jiān)測(cè)，并結(jié)合聚類分析手段對(duì)開(kāi)采前后擾動(dòng)區(qū)域進(jìn)行分析，得出以下結(jié)論：

（1） 運(yùn)用K-means算法對(duì)首采區(qū)一定時(shí)間內(nèi)的所有微震事件的點(diǎn)進(jìn)行了聚類分析。探尋空間內(nèi)微震事件點(diǎn)的聚集度，用聚類中心點(diǎn)進(jìn)行描述，得出了聚類中心點(diǎn)圍繞在3 736 m拉底水平和3 720 m出礦水平上下。說(shuō)明采礦活動(dòng)是引起微震事件的主要因素，同時(shí)周圍的擴(kuò)散可以看出，采礦擾動(dòng)導(dǎo)致的應(yīng)力重分布的范圍，對(duì)該重點(diǎn)區(qū)域應(yīng)加強(qiáng)監(jiān)測(cè)。

（2） 通過(guò)聚類分析，分為3類和分為5類得出的結(jié)果較為一致，說(shuō)明聚類分析中并不是分類越細(xì)致，結(jié)果越可靠。聚類分析必須與實(shí)際相結(jié)合，同時(shí)兩種情況下的分類聚集區(qū)和開(kāi)采活動(dòng)范圍基本對(duì)應(yīng)，能夠?qū)ΦV山安全生產(chǎn)提供指導(dǎo)性意見(jiàn)。

（3） 微震事件時(shí)間分布分析可預(yù)測(cè)巖體穩(wěn)定性，一段時(shí)間內(nèi)的微震事件的突然聚集預(yù)示巖體失穩(wěn)的可能性，分析得出幾個(gè)時(shí)期的活躍期在11月和3月，和實(shí)際大量出礦時(shí)間吻合，此時(shí)應(yīng)該結(jié)合生產(chǎn)實(shí)際，分析微震事件聚集原因，為安全生產(chǎn)提供指導(dǎo)。

參考文獻(xiàn)

[1]? ? 袁海平，曹平.我國(guó)自然崩落法發(fā)展現(xiàn)狀與應(yīng)用展望[J].金屬礦山，2004（08）：25-28.

[2]? ? 沈南山，顧曉春，尹升華.國(guó)內(nèi)外自然崩落采礦法技術(shù)現(xiàn)狀[J].采礦技術(shù)，2009，9（4）：1-4+15.

[3]? ? Korenromp M H A . Comparative study of the S-shaped and Hoek-Brown failure criterion by finite-difference modelling， Northparkes， NSW， Australia：[J]. Delft University of Technology .2013.

[4]? ? 吳愛(ài)祥，武力聰，劉曉輝，等.礦山微地震活動(dòng)時(shí)空分布[J].北京科? 技大學(xué)學(xué)報(bào)，2012，34（6）：609-613.

[5]? ?Hudyma MR.Analysis and interpretation of clusters of seismicevents[D].University of Western Australia，2008.

[6]? ? 程愛(ài)平，高永濤，梁興旺，等.基于未確知聚類法的底板采動(dòng)破壞深度動(dòng)態(tài)預(yù)測(cè)[J].采礦與安全工程學(xué)報(bào)，2014，31（5）：739-744.

[7]? ? 劉棟，李夕兵，劉志祥，等.基于STSNN聚類算法的用沙壩礦微震事件活動(dòng)特征研究[J].中國(guó)安全生產(chǎn)科學(xué)技術(shù)，2017，13（2）：74-78.

[8]? ? 魯婧婧，張晉昕，袁向東，等.歐氏距離的加權(quán)處理對(duì)K-means法聚類效果的改進(jìn)[J].中國(guó)醫(yī)院統(tǒng)計(jì)，2008（1）：9-12.

[9]? ? 郭倩，張景飛，郭攀，等.基于K-means的煤巖孔裂隙微觀結(jié)構(gòu)SEM表征[A]. 中國(guó)煤炭學(xué)會(huì)礦山建設(shè)與巖土工程專業(yè)委員會(huì).礦山建設(shè)與巖土工程技術(shù)新進(jìn)展——2017年全國(guó)礦山建設(shè)學(xué)術(shù)年會(huì)論文集[C].中國(guó)煤炭學(xué)會(huì)礦山建設(shè)與巖土工程專業(yè)委員會(huì)：中國(guó)煤炭學(xué)會(huì)，2017：6.

[10]? 劉歡，閻鎮(zhèn).基于改進(jìn)的K-means法的高分辨率遙感影像道路提取[J].計(jì)算機(jī)與現(xiàn)代化，2017（11）：55-61.

[11]? JOUBERT P J. Microseismic monitoring of hydraulic fractures in block cave mines[J]. Transactions of the Institutions of Mining and Metallurgy，Section A：Mining Technology，2010，119（3）：193-197.

[12]? REYES-MONTES J M，PETTITT W S，PIERCE M E，et al. Microseismic validation of jointed rock models in cave mining[C]. 44th US Rock Mechanics Symposium，2010，63-74.

[13]? 陳通，王雙明，王悅，等.長(zhǎng)壁工作面采動(dòng)引起的微震活動(dòng)分布規(guī)律研究[J].采礦與安全工程學(xué)報(bào)，2018，35（4）：795-800.

[14]? 程關(guān)文，王悅，馬天輝，等.煤礦頂板巖體微震分布規(guī)律研究及其在頂板分帶中的應(yīng)用——以董家河煤礦微震監(jiān)測(cè)為例[J].巖石力學(xué)與工程學(xué)報(bào)，2017，36（S2）：4036-4046.

[15]? 李彪，徐奴文，戴峰，等.烏東德水電站地下廠房開(kāi)挖過(guò)程微震監(jiān)測(cè)與圍巖大變形預(yù)警研究[J].巖石力學(xué)與工程學(xué)報(bào)，2017，36（S2）：4102-4112.

[16]? 徐奴文，戴峰，李彪，等.猴子巖水電站地下廠房開(kāi)挖過(guò)程微震特征與穩(wěn)定性評(píng)價(jià)[J].巖石力學(xué)與工程學(xué)報(bào)，2016，35（S1）：3175-3186.