劉德彪，李夕兵，李 響，尚雪義

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

0 引言

微震監(jiān)測技術(shù)在深部開采的礦山中，得到了越來越廣泛的應(yīng)用[1-4]。微震事件數(shù)目龐大，其區(qū)域分布特征人工劃分具有較大的主觀性。聚類分析借助聚類算法實現(xiàn)微震事件的劃分，可降低人工劃分的主觀性，發(fā)現(xiàn)潛在的微震集群，從而有效的分析微震事件分布特征和活動規(guī)律。

目前，國外學者對地震事件聚類分析進行了較多研究，且在國內(nèi)也越來越受到關(guān)注。Zaliapin等[5]采用K-means對人為和地質(zhì)構(gòu)造引起的地震事件進行了區(qū)分；Weatherill和Burton[6]采用K-means證明了震源模型可由地震目錄信息計算得到；Morales等[7]提出1種基于自適應(yīng)馬氏距離的K-means聚類方法，該方法可用于球形簇和橢球形簇地震聚類；Ramdani等[8]將地震深度屬性引入到K-means聚類分析，改善了地震事件動態(tài)變化過程的圖像分辨率；吳愛祥等[9]用最短距離聚類法分析了礦山微震活動的時空分布，區(qū)分了礦山的微震聚集區(qū)；Wang等[10]采用模糊C均值聚類得出了微震事件的活動特征與三維波速之間的關(guān)系；Shang等[11]提出采用S-KL指標選擇最佳聚類數(shù)，并解釋了聚類簇與地質(zhì)構(gòu)造的關(guān)系；劉棟等[12-13]采用時空共享近鄰聚類算法(STSNN)和具有噪聲的基于密度的聚類方法(DBSCAN)分析了巖體的活動性。

由上述分析可知，在地震聚類分析中K-means運用最為廣泛。K-means聚類算法具有較高的計算效率和較強的靈活性，適用于球形簇結(jié)構(gòu)的數(shù)據(jù)集，但數(shù)據(jù)集中的異常事件會對算法的結(jié)果產(chǎn)生較大的影響。初始聚類中心的選擇也容易影響聚類結(jié)果，使得聚類算法很難達到全局最優(yōu)值。Wang等[14]提出了1種改進的K-means聚類算法，先用基于局部離群因子(Local Outlier Factor，LOF)的方法檢測異常并移除事件，然后用所有數(shù)據(jù)的均值作為第1個初始聚類中心，再依次計算其他的初始聚類中心，該方法提高了聚類的準確性。

本文借鑒Wang等[14]提出的算法，提出了1種新的基于LOF的K-means聚類算法(LOF-K-means)，并用該方法對礦山微震事件的分布特征進行分析。首先，采用LOF檢測離群微震事件和選取初始聚類中心，再利用Krzanowski-Lai(KL)指標確定最佳聚類分組數(shù)。采用簇內(nèi)誤差平方和(within-cluster Sum of Squared Errors, SSE)比較本文方法、文獻[14] K-means聚類方法和傳統(tǒng)K-means聚類的聚類效果。最后，采用本文方法對用沙壩礦微震事件進行聚類分析，根據(jù)聚類簇的分布特征對礦山的微震活動性作出評價。

1 局部離群因子算法

局部離群因子算法[15]是1種基于密度來進行異常事件檢測的算法。局部離群因子的大小反映對象xi對于局部中心的偏離程度，局部離群因子的值越大，說明對象xi偏離局部中心的程度越多，所在位置局部密度越?。痪植侩x群因子的值越小，說明對象xi偏離局部中心的程度越少，所在位置局部密度越大，越接近局部中心。在本文中，用dist(xi,p)表示對象xi與對象p之間的距離，局部離群因子LOFk(xi)的具體定義如下：

定義1：對象p的第k距離k-dist(p)在數(shù)據(jù)集X中，若滿足：

1)存在至多k-1 個對象xi′∈X(xi′≠xi)且dist(xi′,p)

2)存在至少k個對象xi′∈X(xi′≠xi)且dist(xi′,p)≤dist(xi,p)，則對象p的第k距離記為k-dist(p)=dist(xi,p)。圖1說明了k=6時，對象p1的第k距離k-dist(p1)。

圖1 k=6時的第k距離領(lǐng)域與可達距離Fig.1 k-distance neighbourhood and reachability distance, when k=6

定義2：對象p的第k距離領(lǐng)域Nk(p)

Nk(p)指對象p的第k距離內(nèi)所有對象組成的集合(包括第k距離)。|Nk(p)|是指對象p的第k距離領(lǐng)域內(nèi)所有對象的個數(shù)，且|Nk(p)|≥k。圖1對象p1的|Nk(p1)|為7，對象p2的|Nk(p2)|為6。

定義3：對象xi相對于對象p的可達距離

令k為正整數(shù)，對象xi相對于對象p的可達距離計算如下：

reach-distk(xi,p)=max(k-dist(p),dist(xi,p))

(1)

圖1中，對象xi相對于對象p1的可達距離為reach-distk(xi,p1)=k-dist(p1)，對象xi相對于p2的可達距離為reach-distk(xi,p2)=dist(xi,p2)。

定義4：對象xi的局部離群因子

對象xi的局部離群因子LOFk(xi)的定義如下：

(2)

式中：lrdk(xi)是指對象xi的局部可達密度(Local Reachability Density)，定義如下：

(3)

由公式(2)～(3)可知，對象xi的LOF值越小，說明其局部可達密度越大，越接近局部中心；對象xi的LOF值越大，說明其局部可達密度越小，越接近局部邊緣，為異常事件的可能性越大。

2 聚類算法的優(yōu)化

2.1.1 異常事件檢測

數(shù)據(jù)集中的異常事件容易影響K-means聚類結(jié)果，因此聚類前需檢測和剔除異常事件。首先，利用公式(2)計算數(shù)據(jù)集X中每個對象xi的LOF值，再將所有對象xi的LOF值按升序排列，并進行歸一化處理；然后，將所有對象xi的歸一化值進行升序排序，并計算出其拐點值，將此值作為異常事件和正常事件的臨界值[16]。如果對象xi的歸一化值大于拐點值，則剔除對象xi；反之，則留下對象xi。高斯分布的歸一化公式如下[17]：

(4)

2.1.2 初始聚類中心的選擇

本文選取初始聚類中心的思路為：選擇去除異常事件數(shù)據(jù)集中，LOF值最小的對象xi作為聚類算法的第1個初始聚類中心；然后，將與第1個初始聚類中心距離較遠且全局密度較大的對象作為第2個初始聚類中心；接著，將與前2個初始聚類中心的均值相距較遠且全局密度較大的對象作為下1個初始聚類中心，直到計算得到與聚類分組數(shù)相同的聚類中心數(shù)為止。具體如下：

2)計算第2個初始聚類中心，直到第K個初始聚類中心。

(5)

2.1.3 LOF-K-means聚類算法

LOF-K-means聚類算法實現(xiàn)過程如下：

Input：數(shù)據(jù)集X={x1,x2,…,xn}，聚類分組數(shù)K，第k領(lǐng)域值k。

Step2：用2.1.2節(jié)的方法選取初始聚類中心。

LOF-K-means聚類算法采用拐點值判別異常事件，較人為判別更加客觀。同時，用LOF值最小的對象，即所在區(qū)域密度最大的對象作為第1個聚類中心，可以適應(yīng)更廣泛類型的數(shù)據(jù)集，避免出現(xiàn)數(shù)據(jù)集類型或者大小的改變對聚類算法產(chǎn)生的影響。

2.1.4 評價指標

1)聚類效果評價

采用函數(shù)簇內(nèi)誤差平方和(within-cluster Sum of Squared Errors, SSE)評價聚類效果，SSE值越小說明各類間分隔越明顯，聚類結(jié)果越好。

(6)

2)聚類數(shù)選取指標

采用KL指標確定聚類數(shù)，Krzanowski和Lai[18]通過計算2個連續(xù)不同分組數(shù)的聚類結(jié)果的簇內(nèi)協(xié)方差矩陣的跡來確定最佳聚類分組數(shù)。KL指數(shù)值越高，其對應(yīng)分組數(shù)的聚類結(jié)果越好。對于聚類分組數(shù)K≥2的數(shù)據(jù)集X={x1,x2,…,xn}，xi∈Rd，KL指標的定義如下：

(7)

3 模擬測試

為測試本文方法的優(yōu)越性，選取傳統(tǒng)K-means聚類和文獻[14] K-means聚類算法作為對比。不同聚類方法的聚類效果可能受數(shù)據(jù)集的大小影響，本文選取數(shù)據(jù)集包含100，250，500，750和1 000個對象進行討論。同一數(shù)據(jù)集分別運用上述3種聚類算法進行聚類計算，使用的聚類參數(shù)為x, y軸坐標，第k領(lǐng)域設(shè)為20。首先計算上述3種聚類方法在同一數(shù)據(jù)集、同一個分組數(shù)下的SSE值，SSE值由小至大分別記為①，②和③，根據(jù)表1進行比較得到評分，見表1。例如：LOF-K-means聚類、傳統(tǒng)K-means聚類和文獻[14] K-means聚類算法的SSE值分別為100，200和150時，那么該SSE值由小至大排列后，對應(yīng)表1中的評價工況4，且LOF-K-means聚類、傳統(tǒng)K-means聚類和文獻[14] K-means聚類算法評分分別為2，0和1。再將每1種聚類算法在2～10個聚類分組下得到的評分相加，得到總評分作為該聚類算法的綜合SSE評價指標，其值越大則說明該聚類算法越好。為減少個別聚類結(jié)果對不同聚類方法的影響，將每種數(shù)量規(guī)模的數(shù)據(jù)集分別隨機生成100次進行聚類計算，得到這3種聚類方法的綜合SSE評分。圖2為數(shù)據(jù)的聚類過程?？芍跏季垲愔行牡倪x取與數(shù)據(jù)集的分布有關(guān)，且與數(shù)據(jù)集的局部密度緊密聯(lián)系。

表1 不同SSE工況下評分Table 1 Scoring for different SSE conditions

注：①～③分別對應(yīng)SSE值從小到大排列的3種聚類方法。

圖2 數(shù)據(jù)集的聚類過程Fig.2 Clustering process of dataset

模擬測試的結(jié)果如圖3所示，當數(shù)據(jù)集事件數(shù)為100，250和1 000時，LOF-K-means聚類綜合SSE評分的中位數(shù)最大；當事件數(shù)為500和750時，其中位數(shù)與文獻[14] K-means聚類的相同，上四位數(shù)與文獻[14] K-means聚類的相同或較大，兩者均比傳統(tǒng)K-means聚類的評分大。

圖3 不同數(shù)據(jù)集大小的聚類結(jié)果Fig.3 Clustering results for different sizes of dataset

總的來說，對于不同事件數(shù)大小的數(shù)據(jù)集，LOF-K-means聚類算法的聚類效果最好，文獻[14] 的聚類方法次之，傳統(tǒng)K-means聚類方法效果最差。

4 礦山微震監(jiān)測應(yīng)用

4.1 用沙壩礦微震系統(tǒng)

本文實例所用的數(shù)據(jù)來源于用沙壩礦微震監(jiān)測系統(tǒng)，其位于貴州省中部，東經(jīng)106.81°，北緯27.08°。礦體呈穩(wěn)定的層狀，礦體厚度穩(wěn)定，沿走向和傾向連續(xù)性較好，傾角為10°～55°，磷礦的年產(chǎn)量超過200萬t且已探明的礦石儲量為4 000多萬t。2013年開始，在用沙壩礦區(qū)建立了礦山IMS微震監(jiān)測系統(tǒng)，主要用來監(jiān)測礦區(qū)內(nèi)的微震活動。IMS微震監(jiān)測系統(tǒng)由28個傳感器組成 ，如圖4(a)中三角形區(qū)域，其中單向傳感器有26個，三向傳感器有2個；主要分布在920中段，1080中段和1120中段。圖4(a)說明礦山微震事件主要分布在巷道區(qū)域，圖4(b)為礦山生產(chǎn)活動區(qū)域和斷層分布，可知微震事件分布與生產(chǎn)活動區(qū)域和斷層分布具有較好的吻合性。一般認為較大震級微震主要受斷層影響，本文選取了較大震級的事件(M≥0)展開分析[11]，嘗試將微震事件與斷層活動聯(lián)系起來。

圖4 用沙壩礦區(qū)的微震系統(tǒng)與斷層分布Fig.4 Microseismic system and fault distribution in Yongshaba mine

4.2 礦山微震事件的聚類分析

本文使用的數(shù)據(jù)集為2014年1-6月測得的1 649個矩震級大于等于0級的微震事件，如圖4(a)中圓形所示。采用LOF-K-means聚類算法，使用的聚類參數(shù)為微震事件的x,y,z軸坐標，所取的第k領(lǐng)域為k=90。聚類后的KL指數(shù)如圖5所示，可知較好的聚類分組數(shù)為2，5和7。

圖5 LOF-K-means聚類的微震事件KL指數(shù)Fig.5 KL index of microseismic events by LOF-K-means clustering

圖6給出了聚類分組數(shù)K為2，5和7的微震事件聚類結(jié)果，參照文獻[11]，選取的一些重要聚類參數(shù)值見表2～4。從圖6(a)可知，K=2時，礦區(qū)的微震事件由斷層F310c和斷層F350劃分為左右兩簇。分析可知，C1簇主要受主斷層F310，F(xiàn)313和斷層F316，F(xiàn)331的影響；C2簇主要受主斷層F309，F(xiàn)350和F302的影響。結(jié)合表2，可以解釋這2個聚類簇主要依據(jù)區(qū)域斷層結(jié)構(gòu)間的作用，而引起微震事件的影響程度進行劃分。C1簇和C2簇的mEs /Ep值大于10，且C2簇的(Es/Ep)0.5為8.07比C1簇的大，說明C2簇受斷層的影響更大。

表2 K=2時，不同聚類簇的微震參數(shù)Table 2 Microseismic parameter of different clusters，when K=2

表3 K=5時，不同聚類簇的微震參數(shù)Table 3 Microseismic parameter of different clusters，when K=5

表4 K=7時，不同聚類簇的微震參數(shù)Table 4 Microseismic parameter of different clusters，when K=7

注：N指簇內(nèi)微震事件數(shù)；Mmax指簇內(nèi)最大的震級；N1指簇內(nèi)震級≥1.0的微震事件數(shù)；N1.5指簇內(nèi)震級≥1.5的微震事件數(shù)；mEs/Ep為簇內(nèi)S波與P波能量比的均值；(Es/Ep)0.5為簇內(nèi)S波與P波能量比的中位數(shù)；Mmean為簇內(nèi)平均震級。

圖6 用沙壩礦微震事件LOF-K-means的聚類結(jié)果Fig.6 LOF-K-means clustering results of microseismic events in Yongshaba mine

從圖6(b)可知，K=5時，C1簇受斷層I的影響，C2簇和C3簇沿著主斷層F310a，F(xiàn)316和F313劃分；C4簇受主斷層F350和F302影響；C5簇受主斷層F309和斷層XVII，XVI的影響。結(jié)合表3，C1簇的事件數(shù)最少，但N1.5有5個；mEs/Ep大于10，(Es/Ep)0.5與其他簇相比處于較大值且Mmean最大，說明C1簇主要受斷層滑移的剪切作用影響，推測主斷層F310或F331可能延伸到C1簇區(qū)域。C2簇的事件數(shù)最多，mEs/Ep大于10，N1有18個，N1.5有2個，說明其受到斷層滑移的影響；但Mmean最小且(Es/Ep)0.5較小，說明其也受到礦山生產(chǎn)活動的影響，且影響作用比斷層滑移的大。C3簇與C2簇類似，受到礦山生產(chǎn)活動和斷層滑移共同影響。C4簇與C5簇所在區(qū)域較難描述，但C5簇mEs/Ep大于10且(Es/Ep)0.5接近10，說明其基本受斷層滑移的剪切作用影響。

從圖6(c)可知，K=7時，C1簇與K=5時的C1簇基本相同，C2簇主要受主斷層F310a和F331的影響，C3簇與C4簇沿著主斷層F313和F316劃分；C5簇受主斷層F350的影響；C6簇受主斷層F302的影響且沿著F309與C7簇劃分。C2 簇與C3簇的各項微震參數(shù)基本相同，且與K=5時的C2簇相似，說明主要受到各斷層滑移和礦山生產(chǎn)活動共同影響。C4簇與C3簇類似，但Mmax和Mmean較大，說明其受到各斷層滑移和礦山生產(chǎn)活動共同影響，相對受斷層滑移的影響程度較大。C5簇的事件數(shù)較少但N1.5有2個，且mEs/Ep和(Es/Ep)0.5都大于3不到10，說明該區(qū)域由斷層滑移和礦山生產(chǎn)活動共同作用；且從微震事件數(shù)量上看，礦山生產(chǎn)的影響較小，斷層滑移的影響較大。C6簇可作為由礦山生產(chǎn)造成微震事件典型的集群，事件數(shù)較多且Mmean和Mmax最小，說明其由礦山生產(chǎn)活動影響，雖然其mEs/Ep和(Es/Ep)0.5均大于3，但綜合分析仍可認為C6簇主要處于礦山生產(chǎn)活躍區(qū)域。C7簇mEs/Ep和(Es/Ep)0.5均大于10，N1.5有2個，Mmean為0.31，說明主要受斷層滑移剪切作用的影響?？偟膩碚f，K=7時各簇微震事件的分布效果較好，可以更好的解釋微震活動性特征。

K=7時聚類簇事件的時鐘矢量圖如圖7所示。由圖7可知，C1，C5和C6簇微震事件初期發(fā)生時間集中在12∶00方向，后逐漸轉(zhuǎn)向13∶00～14∶00；C2，C3，C4和C7簇微震事件發(fā)生時間主要在13∶00～14∶00。7個聚類簇微震事件的軌跡曲線都超出圓，說明礦區(qū)微震事件受生產(chǎn)活動影響。用沙壩礦區(qū)的生產(chǎn)爆破時間集中在11∶00～13∶00，則爆破期間和爆破后的1～3 h內(nèi)是微震事件頻發(fā)期。

圖7 K=7時聚類簇事件的時鐘矢量Fig.7 Clock vectors of clustering events ,when K=7

5 結(jié)論

1)針對K-means聚類易受異常事件和初始聚類中心影響的問題，引入了LOF算法進行異常事件的檢驗和初始聚類中心的選擇，提高了聚類結(jié)果的有效性。

2)建立了聚類算法的綜合SSE評價指標，通過計算模擬比較了在不同數(shù)據(jù)集大小下，LOF-K-means聚類算法、文獻[14] K-means聚類算法和傳統(tǒng)K-means聚類算法的綜合SSE評價指標，得到LOF-K-means聚類算法最優(yōu)。

3)將LOF-K-means聚類算法用于分析用沙壩礦微震事件分布特征，得出最佳分組數(shù)為7。其中C1，C7簇主要受斷層滑移的影響；C6簇主要受礦山生產(chǎn)活動的影響，為礦山微震活動性分析提供了一種有效的方法。