基于STSNN聚類算法的用沙壩礦微震事件活動特征研究*

劉　棟，李夕兵，劉志祥，董隴軍，周勇勇，陳光輝

(中南大學 資源與安全工程學院，湖南 長沙 410083)

0　引言

采礦活動向深部發展，破壞了深部的原巖應力，由此引起的地質災害狀況頻發，微震監測技術是目前監測災害的有效手段[1-2]。近幾年，微震監測技術已經在臺站布網優化[3-4]、P波[5-6]、S波拾取[7]、微震定位[8-9]及微震爆破事件識別[10]等方面有了較大的進展，并取得了一系列重要的研究成果。為了掌握巖體破壞的規律，對微震事件的聚集特征及活動規律研究顯得尤為重要。在此方面，國內外學者做了許多創新性的研究，例如吳愛祥[11]等利用最短距離聚類法對某礦山微地震活動的時空分布進行了分析，識別了井下微震活動聚集區域；謝和平[12]等采用分形理論研究微震活動的時空分布，發現分形維數在一定程度上反映了微震活動的變化狀態；Martin[13]利用空間聚類方法對微震事件進行聚類分析并解釋了微震活動規律。縱觀這些研究仍存在一定的不足：在劃定微震活動的聚集區域時只考慮微震事件的空間分布特征，沒有將時間特征考慮在內。對此，本文采用時空共享近鄰聚類(STSNN)算法，將微震事件的時間空間特征相結合并展開聚類研究。2012年Liu等[14]將共享近鄰概念引入時空聚類，提出STSNN聚類算法，該算法采用實體的有效近鄰數目來定義其時空密度，更適于發現不同密度的時空簇，是時空聚類算法中較先進的聚類算法。

開磷用沙壩礦存在嚴重的潛在地壓問題，于2013年安裝使用IMS微震監測系統來監測井下災害。本文以用沙壩礦為工程背景，以用沙壩礦IMS微震監測系統2014年收集到微震數據作為研究對象，通過P波、S波拾取后得到5 425個有效微震作為本文的數據集。采用STSNN算法并選擇合適的參數對其進行聚類分析，計算其聚類有效性，并對最大類簇的微震事件活動規律進行探究，根據其變化規律發現產生較大事件的預警點。通過對微震事件活動規律的研究可為大事件的產生提供有效的預判信息，為保證礦山安全生產發揮重要的作用。

1　時空共享近鄰聚類(STSNN算法)

1.1　基本定義

為了闡述時空近鄰共享聚類算法，首先介紹幾個重要的概念：

定義1時空事件：一個時空事件STi是具有空間坐標和確定的時間標志，表示為STi(Li,ti)。

定義2時空近鄰：時空事件STi的時空近鄰定義為以ΔD作為底面半徑、2ΔT作為高的圓柱體，ΔT為時間窗口。

定義3時空密度：時空事件STi的時空密度定義為STi時空直接可達的事件數量。

定義4時空可達：假如時空事件數據集ST1,ST2,…STi,…STn，STn是一個時空核心，如果STi-1是從STi時空直接可達，那么ST1也可視為從STn時空可達，表示為STn→ST1。

定義5噪聲：假設時空事件STi不是一個時空核心，在精確時空近鄰內不存在時空核心，那么可將STi視為噪聲。

定義6時空類簇：對于一個時空數據集STD，時空類簇STC是滿足下列2個條件的非空集。

條件1：?STi∈STD,STj∈STD:STi→STj

條件2：?STi∈STD,STj∈STD:STi～STj

1.2　STSNN算法及參數選擇

STSNN算法需要用戶設置4個參數：k,kT,MinPts和ΔT，k是時空近鄰的事件數量，kT是判別兩個事件是否時空直接可達的閾值，MinPts用來識別時空核心，ΔT是定義時空近鄰的時間窗口。基于上述定義，STSNN算法的執行可分為以下5個步驟：

Step1：對每個時空事件確定時空近鄰k。

Step2：對每個時空事件搜索其時空共享近鄰。

Step3：計算每個時空事件的時空密度以及進一步探測時空核心。

Step4：擴大時空類簇。選擇一個時空核心并使用遞歸算法將所有時空可達和時空相連的事件增加到時空類簇。

Step5：識別噪聲事件。

STSNN參數kT和MinPts通常設置為0.5k。另外2個參數，k和ΔT在本次研究中，ΔT初選值設定為[5,6]，k的初選值[5,6,7…17,18,19]，共有30種組合，為了方便后面表達，每一種組合的類別符合的如表1所示。

表1　30種組合的類別符號

2　聚類結果及有效性評價

分別利用STSNN算法對上述k和ΔT的30種組合進行聚類運算。當k和ΔT分別選擇不同的值時，可以得到不同的聚類類簇，如圖1所示，當ΔT一定k增大時，得到的類簇數量越來越少，而最大類簇事件數沒有呈現規律性變化。圖2表示不同類別的噪聲率，除ID1(k=5，ΔT=5)和ID16(k=5，ΔT=6)外，其他類別的噪聲率相對較低。

圖1　不同類別的簇數及最大類簇事件數Fig.1　The number of clusters of different classes and the maximum number of cluster events

圖2　不同類別的噪聲率Fig.2　Noise rates of different classes

聚類有效性評價是為了確定對數據劃分多少個類簇是最合適的。一個理想的聚類結果應該滿足2方面要求[15]:(1)緊密性，即簇內實體應盡可能相似；(2)分離性，即不同簇間實體差異盡可能大。一個好的聚類結果應該同時滿足緊密性與分離性。本文選擇RMSSDT指數[16]和PBM指數[17]作為有效性評價指標。RMSSDT指數旨在度量實體間的離散程度，其值越小，則表示簇的均質性越好；PBM指數采用融合多個度量指標的策略來進行聚類結果的有效性度量，其值越大，表示聚類結果越可靠，如圖3所示。

注：折線代表PBM；柱狀圖代表RMSSDT圖3　聚類有效性評價結果Fig.3　Clustering effectiveness evaluation results

利用RMSSDT指數和PBM指數對每一類別的有效性進行評價，如圖3所示單純從PBM評價指標(折線)來看，ID17(k=6，ΔT=6)、ID19(k=8，ΔT=6)和ID21(k=10，ΔT=6)聚類效果較好，優勢較明顯，進而從RMSSDT指數(柱狀圖)來比較，ID17(k=6，ΔT=6)聚類結果最好，并且其噪聲率不高，因此本文后續研究將利用ID17(k=6，ΔT=6)的聚類結果。該聚類共得到98個類簇，最大類簇的事件數為544，噪聲率為24.8%。最大的5個類簇事件數分別為544，246，234，161和100。圖4表示最大類簇在礦山中的主要聚集區域，從圖中可以發現，該類簇主要集中在用沙壩礦的斷層附近。

圖4　最大類簇微震事件的主要聚集區域Fig.4　The events aggregation area of the biggest cluster group

3　微震事件活動特征

3.1　微震事件24 h分布

圖5　微震事件24 h分布(m表示震級)Fig.5　The distribution of micro-seismic in 24 hours

圖5表示該區域聚類事件的24 h分布圖，微震事件主要發生在每天的11時到15時，大于0級的事件也主要發生在此時間段。用沙壩礦爆破時間主要集中在每天11時到14時。因此該區域的微震事件大部分分布在斷層區域，受采礦活動誘發而產生。在非爆破時間，該區域的微震事件在時間分布上較為均勻，但偶爾也有震級大于-0.5級事件發生。

3.2　微震事件活動率與視體積

微震事件活動率反映了巖體內部微裂紋的擴展變化趨勢，表現了微裂紋的產生和發展速度，微震事件的突變特征體現了巖體開挖導致圍巖破壞的前兆信息。微震事件的視體積表征地震非彈性變形的體積與它釋放的能量和體變勢相關，累計視體積隨時間變化曲線的斜率通常被認為是表示巖體應變速率的重要指標，曲線的斜率變化反映應變率的微小變化，巖體失穩變形的前兆是經歷一段時間的加速變形過程。通過觀察圖6，可以按微震活動率分為3個階段，AB階段微震事件活動率降低，說明此階段可能有大事件發生，能量得到釋放；BC段微震活動率變化緩慢，說明此階段巖體處于穩定狀態，開始生成新的裂隙并呈穩定發展階段；CD階段為巖體內裂隙發展迅速，微震活動率增加，因此可能發生較大事件。累積視體積曲線在B點(5月3日附近)和C(5月10日附近)點明顯上升，表明B，C點附近可能有破壞發生，導致巖體失穩。在CD段開始急劇上升，表明此階段大量的微裂隙張開占據巖石空間導致巖石體積膨脹，視體積曲線突然增大的點可以作為發生大事件的前兆信息，因此在D(5月16日)點之后，很可能產生較大的破壞性事件。實際上在5月19日現場有一次震級為0.7級的微震事件發生。因此可以認為微震事件率急劇下降并且累積視體積曲線忽然上升是巖石發生破壞的預兆現象。

圖6　微震活動率與累積視體積事件序列曲線Fig.6　Time series curve of cumulative apparent volume and event activity

3.3　施密特數與勁度系數

施密特數(Schmidt Number)是運動學黏度與擴散率之比，能測量震動巖石流變的時空復雜程度，也是反映巖石不穩定的參數，施密特數越低，穩定性越差。勁度系數(Stiffness Modulus)表征系統阻止地震變形和應力增加的能力。如圖7所示施密特數曲線與勁度系數曲線保持相同的變化趨勢，分別在5月2日、5月9和5月16日附近急劇上升然后迅速下降，降低程度幾乎相同，施密特數處于較低水平時該區域的巖體穩定性較差，產生破壞的可能性比較大。在這3個時間，從勁度系數曲線可以看出，巖體阻止系統變形的能力也瞬間下降，進而也表明該區域在這3個時間較容易發生危險事故。而實際現場監測中，5月10日期間，發生了2次震級超過0.5級的微震事件，與5月9日附近施密特數勁度系數先升后降的預測基本對應。因此可以把施密特數和勁度系數先升后降的點作為大事件將要發生的預警信息。

圖7　施密特數與勁度系數事件序列曲線Fig.7　Time series curve of Schmidt Number and Stiffness Modulus

4　結論

1)利用時空共享近鄰(STSNN)聚類算法對用沙壩礦微震事件進行聚類分析，可以發現微震事件的時空聚集類簇，并對其聚類有效性進行分析，確定聚類參數選擇的合理性。對微震事件聚類分析來確定所研究的區域，比人為劃定的事件聚集區域更客觀有效。

2)對聚類最大類簇的微震事件活動特征進行分析，通過24 h分布圖可以發現用沙壩礦微震事件主要發生在每天的11時到15時，大于0級的時間也主要發生在此時間段，微震活動率急劇下降并且累積視體積曲線忽然上升、施密特數和勁度系數先升后降的點作為巖體失穩發生破壞性事件的預警點。

[1]李夕兵. 巖石動力學基礎及應用[M]. 北京，科學出版社，2014.

[2]王春來,吳愛祥,劉曉暉.深井開采巖爆災害微震監測預警及控制技術[M].北京,冶金工業出版社,2013.

[3]周勇勇,李夕兵,劉志祥,等. 高精度三維地下微震臺網模糊優選方法[J]. 中國安全生產科學技術,2016, 12(7):82-86.

ZHOU Yongyong, LI Xibing, LIU Zhixiang, et al.A fuzzy optimum approach for high precision three-dimensional underground micro-seismic Network [J]. Journal of Safety Science and Technology, 2016, 12(7):82-86.

[4]張楚旋,李夕兵,董隴軍,等. 微震監測傳感器布設方案評價模型及應用[J]. 東北大學學報(自然科學版),2016(4):594-598,608.

ZHANG Chuxuan, LI Xibing, DONG Longjun, et al. Evaluation model of micro-seismic monitoring sensor layout scheme and its application[J].Journal of Northeastern University (Natural Science), 2016(4):594-598,608.

[5]LI X, SHANG X, WANG Z, et al. Identifying P-phase arrivals with noise: An improved Kurtosis method based on DWT and STA/LTA[J]. Journal of Applied Geophysics, 2016, 133: 50-61.

[6]LI X, SHANG X, Morales-Esteban A, et al. Identifying P phase arrival of weak events: the Akaike Information Criterion picking application based on the Empirical Mode Decomposition [J]. Computers & Geosciences, 2016.

[7]張楚旋,李夕兵,董隴軍,等. 三函數四指標礦震信號S波到時拾取方法及應用[J]. 巖石力學與工程學報,2015(8):1650-1659.

ZHANG Chuxuan, LI Xibing, DONG Longjun, et al. A S-wave phase picking method with four indicators of three functions for micro-seismic signal in mines [J]. Chinese Journal of Rock Mechanics and Engineering, 2015(8):1650-1659.

[8]Linqi HUANG, Xibing LI, Longjun DONG, et al. Relocation method of micro-seismic source in deep mines [J]. Transactions of Nonferrous Metals Society of China, 2016(11):2988-2996.

[9]張飛, 劉德峰, 張衡,等. 基于IMS微震監測系統的微震事件定位精度分析[J]. 中國安全生產科學技術, 2013, 9(6):21-26.

ZHANG Fei, LIU Defeng, ZHANG Heng, et al. Location accuracy analysis of the seismic events based on the IMS micro-seismic monitoring system [J]. Journal of Safety Science and Technology, 2013, 9(6):21-26.

[10]GUOyan ZHAO, Ju MA, Long-jun DONG, et al. Classification of mine blasts and micro-seismic events using starting-up features in seismograms[J]. Transactions of Nonferrous Metals Society of China, 2015(10): 3410-3420.

[11]吳愛祥,武力聰,劉曉輝,等. 礦山微地震活動時空分布[J]. 北京科技大學學報,2012,34(6):609-613.

WU Aixiang，WU Licong，LIU Xiaohui，et al. Space-time distribution of micro-seismic activities in mines[J].Journal of University of Science and Technology Beijing，2012,34(6):609-613.

[12]謝和平,W. G. Pariseau. 巖爆的分形特征和機理 [J]. 巖石力學與工程學報,1993,12(1):28-37.

XIE Heping，Pariseau W G. Fractal character and mechanism of rockbursts [J]. Chinese Journal of Rock Mechanics and Engineering,1993,12(1): 28-37.

[13]Hudyma MR. Analysis and interpretation of clusters of seismic events [D]. University of Western Australia, Perth, Australia, 2008.

[14]Liu Q., Deng M., Bi J., Yang, W. A Novel Method for Discovering Spatio-temporal Clusters of Different Sizes, Shapes and Densities in the Presence of Noise [J].International Journal of Digital Earth，2012, 7(2): 138-257.

[15]Berry M, Linoff G. Data Mining Techniques: For Marking, Sales and Customer Support[M]. New York: Wiley, 1997.

[16]Halkidi M, Batistakis Y, Vazirgiannis M. Clustering algorithm and validity measures[A]. Proceedings of 13th International Conference on Scientific and Statistics Database Management[C]. 2001, 3-22.

[17]Pakhira M K, Bandyopadbyay A, Maulik U. Validity index for crisp and fuzzy clustering[J]. Pattern Recognition, 2004, 37(3):487-501.