基于LOF的K-means聚類方法及其在微震監(jiān)測(cè)中的應(yīng)用*

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

(中南大學(xué) 資源與安全工程學(xué)院，湖南 長(zhǎng)沙 410083)

0 引言

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

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

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

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

1 局部離群因子算法

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

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

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

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

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

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

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

定義3：對(duì)象xi相對(duì)于對(duì)象p的可達(dá)距離

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

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

(1)

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

定義4：對(duì)象xi的局部離群因子

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

(2)

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

(3)

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

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

2.1.1 異常事件檢測(cè)

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

(4)

2.1.2 初始聚類中心的選擇

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

2)計(jì)算第2個(gè)初始聚類中心，直到第K個(gè)初始聚類中心。

(5)

2.1.3 LOF-K-means聚類算法

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

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

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

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

2.1.4 評(píng)價(jià)指標(biāo)

1)聚類效果評(píng)價(jià)

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

(6)

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

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

(7)

3 模擬測(cè)試

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

表1 不同SSE工況下評(píng)分Table 1 Scoring for different SSE conditions

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

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

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

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

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

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

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

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

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

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

本文使用的數(shù)據(jù)集為2014年1-6月測(cè)得的1 649個(gè)矩震級(jí)大于等于0級(jí)的微震事件，如圖4(a)中圓形所示。采用LOF-K-means聚類算法，使用的聚類參數(shù)為微震事件的x,y,z軸坐標(biāo)，所取的第k領(lǐng)域?yàn)閗=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é)果，參照文獻(xiàn)[11]，選取的一些重要聚類參數(shù)值見表2～4。從圖6(a)可知，K=2時(shí)，礦區(qū)的微震事件由斷層F310c和斷層F350劃分為左右兩簇。分析可知，C1簇主要受主斷層F310，F(xiàn)313和斷層F316，F(xiàn)331的影響；C2簇主要受主斷層F309，F(xiàn)350和F302的影響。結(jié)合表2，可以解釋這2個(gè)聚類簇主要依據(jù)區(qū)域斷層結(jié)構(gòu)間的作用，而引起微震事件的影響程度進(jìn)行劃分。C1簇和C2簇的mEs /Ep值大于10，且C2簇的(Es/Ep)0.5為8.07比C1簇的大，說明C2簇受斷層的影響更大。

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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