蘇倸玉，李 忠，朱 婷，張 偉

(防災(zāi)科技學(xué)院 應(yīng)急管理學(xué)院，河北 燕郊065201)

0 引 言

K-means++算法是2007年由David Arthur和Sergei Vassilvitskii提出的，是K-means算法的改進(jìn)版本。傳統(tǒng)K-means算法的聚類效果以及運(yùn)行時(shí)間在很大程度上受到初始聚類中心選擇的影響，Kmeans++算法對(duì)此進(jìn)行了改進(jìn)。雖然K-means++算法在初始化簇中心時(shí)，增加了計(jì)算量，但在整個(gè)聚類過程中，能夠顯著地提升計(jì)算效率和改善聚類結(jié)果誤差［1］，被廣泛應(yīng)用于分類［2］、聚類［3-4］等領(lǐng)域。

但是，由于K-means++依舊存在難以確定合適的K值問題，可能導(dǎo)致聚類結(jié)果會(huì)產(chǎn)生局部最優(yōu)解?；诖?，本文針對(duì)K值難于確定的問題，采用3種方法聯(lián)合確定最佳的聚類中心個(gè)數(shù)，從而改進(jìn)Kmeans++算法，并將優(yōu)化算法應(yīng)用于地震地磁數(shù)據(jù)聚類分析中。

1 K-means＋＋聚類算法

經(jīng)典的K-means算法是無監(jiān)督的聚類方法［5］，具有簡(jiǎn)單、高效、易于實(shí)現(xiàn)等特點(diǎn)。局部搜索能力較強(qiáng)，且對(duì)于大數(shù)據(jù)樣本空間的聚類有較高的效率［6］，在數(shù)據(jù)聚類分析領(lǐng)域應(yīng)用廣泛［7］。但是，Kmeans也存在2個(gè)致命缺點(diǎn):一是K值多以使用者的經(jīng)驗(yàn)來確定，因此存在很大的主觀性，在不同應(yīng)用場(chǎng)景使用時(shí)造成很大困擾；二是初始值采用隨機(jī)選取方式，可能造成算法收斂速度較慢、計(jì)算效率較低的問題。

基于上述問題，有學(xué)者提出了K-means++算法。在初始值選取時(shí)，采用樣本點(diǎn)與中心點(diǎn)的距離作為概率值，距離越遠(yuǎn)則被選中的概率越大，距離越近則概率越小，從而有效解決了上述第二個(gè)問題。

K-means++聚類算法實(shí)現(xiàn)步驟如下:

（1）從樣本集X中隨機(jī)選擇一個(gè)初始聚類中心點(diǎn)，加入到聚類中心集C中；

（2）遍歷樣本集，計(jì)算每一個(gè)樣本點(diǎn)xi與已存在聚類中心點(diǎn)ck的距離dis（xi，ck），選取最短距離D（xi），如式（1）所示:

依次計(jì)算每個(gè)樣本點(diǎn)的概率，如式（2）所示:

其中，P為一個(gè)樣本點(diǎn)xi被選中為聚類中心的概率；D（xi）代表一個(gè)樣本點(diǎn)xi到已有聚類中心的最短距離；代表所有樣本點(diǎn)到每個(gè)已存在聚類中心的最小距離的平方和。選擇概率最大的樣本點(diǎn)作為下一個(gè)聚類中心點(diǎn)，直至找到K個(gè)初始聚類中心；

（3）以各個(gè)簇內(nèi)樣本點(diǎn)的中心點(diǎn)更新聚類中心，得到K個(gè)新聚類中心；

（4）計(jì)算各樣本點(diǎn)到K個(gè)中心點(diǎn)的距離，按照最短距離歸屬原則歸類，形成K個(gè)新類；

（5）重復(fù)步驟（3）、（4）直到聚類中心點(diǎn)不變或者達(dá)到迭代次數(shù)結(jié)束計(jì)算，并記錄最終的K個(gè)聚類中心點(diǎn)和簇內(nèi)數(shù)據(jù)。

2 K-means＋＋算法改進(jìn)

從上述介紹可見，K-means++算法解決了Kmeans算法的初始中心選擇問題，但對(duì)K的選擇依然沒有改變思路。為解決這個(gè)問題，本文采用拐肘法（Elbow Method）、輪 廓 系 數(shù) 法（Silhouette coefficient）和Calinski-Harabaz指標(biāo)法（CH）聯(lián)合確定K值。

2.1 拐肘法

拐肘法的核心思想，是通過計(jì)算各個(gè)K值時(shí)的畸變程度，畫出畸變程度與K值的關(guān)系圖，找出拐點(diǎn)?；兂潭染褪敲總€(gè)簇的質(zhì)點(diǎn)與簇內(nèi)樣本點(diǎn)的均方差。一個(gè)簇畸變程度的高低，表明簇內(nèi)成員的松散程度高低。隨著K值的增大，到達(dá)拐點(diǎn)后，畸變程度會(huì)得到很大改善，下降幅度趨于平緩，此時(shí)這個(gè)拐點(diǎn)就可以考慮為聚類性能較好的點(diǎn)，可以確定其對(duì)應(yīng)K值作為聚類個(gè)數(shù)，如式（3）所示:

其中，SSE是畸變程度；ci代表第i個(gè)簇；p是ci中的樣本點(diǎn)；mi是ci的質(zhì)心，即ci中所有樣本的均值。

利用拐肘法獲得的K值記作K1。

2.2 輪廓系數(shù)法

輪廓系數(shù)（Silhouette Coefficient）是由Kaufman等人提出的一種用來評(píng)價(jià)算法聚類質(zhì)量的有效性指標(biāo)。該指標(biāo)結(jié)合了凝聚度和分離度，不僅用以評(píng)價(jià)聚類質(zhì)量，還可用來獲取最佳聚類數(shù)［8］。假設(shè)樣本集X包含n個(gè)樣本點(diǎn)x1，x2，....，xn，將樣本集分為K個(gè)簇，每個(gè)樣本點(diǎn)的輪廓系數(shù)如式（4）所示:

其中，Si為第i樣本點(diǎn)的輪廓系數(shù)；ai是樣本點(diǎn)i到其所屬簇中所有其它點(diǎn)的平均距離；bi為樣本點(diǎn)i到其它簇中所有點(diǎn)的最小距離。Si介于［-1，1］，接近-1則說明樣本點(diǎn)i更應(yīng)該分類到另外的簇，接近0則說明樣本i在2個(gè)簇的邊界附近，越趨近于1則聚類效果越優(yōu)。

所有點(diǎn)的輪廓系數(shù)求平均，得到最終的平均輪廓系數(shù)。對(duì)于現(xiàn)有的分類數(shù)，求取輪廓系數(shù)的最大值Sk，與之對(duì)應(yīng)的K值就是最佳聚類數(shù)［9］。

利用輪廓系數(shù)法獲得的K值記作K2。

2.3 CH指標(biāo)法

CH指標(biāo)由分離度與緊密度的比值得到［10］。假設(shè)數(shù)據(jù)集被劃分為K個(gè)類，CH指標(biāo)計(jì)算如式（5）所示:

其中，N為數(shù)據(jù)集總樣本數(shù)；K為類別數(shù)；tr表示對(duì)矩陣求跡，即矩陣主對(duì)角線元素之和；B（K）為類別之間的協(xié)方差矩陣；W（K）為類別內(nèi)部數(shù)據(jù)的協(xié)方差矩陣。

CH越大代表類內(nèi)越密集，類間越分散，聚類結(jié)果性能越好，對(duì)應(yīng)的K值也是最優(yōu)的聚類數(shù)。利用CH指標(biāo)法獲得的K值記作K3。

2.4 聯(lián)合法確定K值

從上述3種求K值的算式和求取過程可以看出:拐肘法實(shí)際上計(jì)算聚類的最小距離，理論上距離越小越好，但聚類個(gè)數(shù)太多就失去意義了。因此，K值是在下降率突然變緩時(shí)刻的值，認(rèn)為此時(shí)聚類效果最佳。輪廓系數(shù)法強(qiáng)調(diào)的是簇內(nèi)部凝聚度最大、簇間分離度最大時(shí)聚類效果最佳，符合客觀實(shí)際，其值域在［-1，1］，越靠近1說明對(duì)應(yīng)的K值越佳。CH指標(biāo)法實(shí)質(zhì)上要求類別內(nèi)部數(shù)據(jù)的協(xié)方差越小越好，類別之間的協(xié)方差越大越好，類似于輪廓系數(shù)法，而且指標(biāo)值越大，對(duì)應(yīng)K值的聚類效果越佳。

鑒于3種方法從不同的角度求取K值，因此本文將3種方法獲得的K值綜合求取平均值，即:

依據(jù)四舍五入原則確定最終的K值，這樣求取的K值既兼顧了簇內(nèi)集中簇間分散的要求，又兼顧了距離最小原則，達(dá)到最優(yōu)聚類結(jié)果。

3 地震地磁前兆數(shù)據(jù)分析

聚類分析方法在地震監(jiān)測(cè)、地震裂縫分布預(yù)測(cè)［11］等地震數(shù)據(jù)分析領(lǐng)域應(yīng)用較多，能夠從大量的地震數(shù)據(jù)中獲取規(guī)律性的知識(shí)，避免了不同學(xué)者的主觀性問題。地震地磁前兆數(shù)據(jù)，是地震臺(tái)站通過地磁儀長(zhǎng)期觀測(cè)到的時(shí)間序列數(shù)據(jù)，單位為奧斯特。根據(jù)地震學(xué)家的研究結(jié)果表明，在地震前1個(gè)月內(nèi)，震中一定范圍內(nèi)會(huì)出現(xiàn)地磁異常變化。因此可以通過長(zhǎng)期的觀測(cè)并通過數(shù)據(jù)挖掘方法搜索離群點(diǎn)，從而發(fā)現(xiàn)地震發(fā)生前期地磁變化，為地震預(yù)測(cè)帶來可能。

3.1 數(shù)據(jù)來源及預(yù)處理

本次實(shí)驗(yàn)數(shù)據(jù)選取2008年1～6月份四川省成都市的地震地磁前兆數(shù)據(jù)。由于地磁前兆數(shù)據(jù)存在數(shù)據(jù)缺失、非標(biāo)準(zhǔn)化等問題，需要進(jìn)行插值、整理、規(guī)范化處理等預(yù)處理操作，以保證數(shù)據(jù)的可用性。缺失數(shù)據(jù)采用均值法補(bǔ)齊，規(guī)范化采用z-score標(biāo)準(zhǔn)化方法進(jìn)行完成，并去掉量綱。經(jīng)過處理的數(shù)據(jù)符合標(biāo)準(zhǔn)正態(tài)分布，按照時(shí)間排列構(gòu)成標(biāo)準(zhǔn)化的數(shù)據(jù)集。

3.2 計(jì)算結(jié)果與分析

首先以拐肘法、輪廓系數(shù)法和CH指標(biāo)法分別計(jì)算K1、K2和K3值，如圖1所示。

圖1 不同算法求得K值折線圖Fig.1 Line graph of K-value by different algorithms

根據(jù)聚類性能標(biāo)準(zhǔn)評(píng)估原理，從圖1（a）的拐肘法得到K1值為4，從圖1（b）的輪廓系數(shù)法得到K2值為5，從圖1（c）的CH指標(biāo)法得到K3值為5。

其次求取平均值作為最優(yōu)K值，K＝（4+5+5）／3＝4.7≈5。因此K＝5是最佳的聚類類別個(gè)數(shù)，此時(shí)可以達(dá)到最好的聚類效果，計(jì)算結(jié)果如圖2所示。

圖2 K-means＋＋的聚類結(jié)果Fig.2 Clustering results of K-means＋＋algorithm

由圖2可見，存在5個(gè)聚類中心。即圖中大圓點(diǎn)，檢測(cè)到2個(gè)異常點(diǎn)，即圖中五角星，正好與2次地震對(duì)應(yīng)。詳盡信息見表1。

表1 兩次地震情況Tab.1 Details of two earthquakes

作為對(duì)比，同樣以K＝5作為聚類數(shù)，利用傳統(tǒng)K-means算法對(duì)成都地震地磁前兆數(shù)據(jù)進(jìn)行聚類分析，結(jié)果如圖3所示。

圖3 K-means的聚類結(jié)果Fig.3 Clustering results of K-means algorithm

從圖3可知，傳統(tǒng)K-means聚類方法沒有找到異常離群點(diǎn)。實(shí)際上，在實(shí)驗(yàn)過程中，K-means聚類方法因初始聚類中心選擇的隨機(jī)性問題，該方法難以發(fā)現(xiàn)離群點(diǎn)，并且每次計(jì)算結(jié)果差異較大。這也說明，同樣的K值下，K-means++聚類方法比Kmeans方法效果更好。

4 結(jié)束語

文中在介紹了K-means++算法基礎(chǔ)上，通過拐肘法、輪廓系數(shù)法和CH指標(biāo)法聯(lián)合求取K值，優(yōu)化了K值的確定方法，從而很好地解決了K-means++聚類算法中K值難以確定的問題。利用改進(jìn)后的K-means++聚類算法對(duì)2008年1～6月份四川省成都市的地震地磁前兆數(shù)據(jù)進(jìn)行聚類分析，發(fā)現(xiàn)了2個(gè)離群點(diǎn)，正好與發(fā)生的2次地震相對(duì)應(yīng)，效果明顯優(yōu)于傳統(tǒng)K-means結(jié)果。得到結(jié)論如下:

優(yōu)化后的K-means++聚類算法優(yōu)于傳統(tǒng)Kmeans算法，K-means++聚類算法可以較好地解決地震地磁的聚類問題。通過離群點(diǎn)檢索可以發(fā)現(xiàn)可能發(fā)生的地震，對(duì)地震監(jiān)測(cè)預(yù)報(bào)工作具有重要的現(xiàn)實(shí)意義。