K-means＋＋算法優化及其在地震前兆分析中的應用研究

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

(防災科技學院 應急管理學院，河北 燕郊065201)

0 引 言

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

但是，由于K-means++依舊存在難以確定合適的K值問題，可能導致聚類結果會產生局部最優解。基于此，本文針對K值難于確定的問題，采用3種方法聯合確定最佳的聚類中心個數，從而改進Kmeans++算法，并將優化算法應用于地震地磁數據聚類分析中。

1 K-means＋＋聚類算法

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

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

K-means++聚類算法實現步驟如下:

（1）從樣本集X中隨機選擇一個初始聚類中心點，加入到聚類中心集C中；

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

依次計算每個樣本點的概率，如式（2）所示:

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

（3）以各個簇內樣本點的中心點更新聚類中心，得到K個新聚類中心；

（4）計算各樣本點到K個中心點的距離，按照最短距離歸屬原則歸類，形成K個新類；

（5）重復步驟（3）、（4）直到聚類中心點不變或者達到迭代次數結束計算，并記錄最終的K個聚類中心點和簇內數據。

2 K-means＋＋算法改進

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

2.1 拐肘法

拐肘法的核心思想，是通過計算各個K值時的畸變程度，畫出畸變程度與K值的關系圖，找出拐點。畸變程度就是每個簇的質點與簇內樣本點的均方差。一個簇畸變程度的高低，表明簇內成員的松散程度高低。隨著K值的增大，到達拐點后，畸變程度會得到很大改善，下降幅度趨于平緩，此時這個拐點就可以考慮為聚類性能較好的點，可以確定其對應K值作為聚類個數，如式（3）所示:

其中，SSE是畸變程度；ci代表第i個簇；p是ci中的樣本點；mi是ci的質心，即ci中所有樣本的均值。

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

2.2 輪廓系數法

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

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

所有點的輪廓系數求平均，得到最終的平均輪廓系數。對于現有的分類數，求取輪廓系數的最大值Sk，與之對應的K值就是最佳聚類數［9］。

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

2.3 CH指標法

CH指標由分離度與緊密度的比值得到［10］。假設數據集被劃分為K個類，CH指標計算如式（5）所示:

其中，N為數據集總樣本數；K為類別數；tr表示對矩陣求跡，即矩陣主對角線元素之和；B（K）為類別之間的協方差矩陣；W（K）為類別內部數據的協方差矩陣。

CH越大代表類內越密集，類間越分散，聚類結果性能越好，對應的K值也是最優的聚類數。利用CH指標法獲得的K值記作K3。

2.4 聯合法確定K值

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

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

依據四舍五入原則確定最終的K值，這樣求取的K值既兼顧了簇內集中簇間分散的要求，又兼顧了距離最小原則，達到最優聚類結果。

3 地震地磁前兆數據分析

聚類分析方法在地震監測、地震裂縫分布預測［11］等地震數據分析領域應用較多，能夠從大量的地震數據中獲取規律性的知識，避免了不同學者的主觀性問題。地震地磁前兆數據，是地震臺站通過地磁儀長期觀測到的時間序列數據，單位為奧斯特。根據地震學家的研究結果表明，在地震前1個月內，震中一定范圍內會出現地磁異常變化。因此可以通過長期的觀測并通過數據挖掘方法搜索離群點，從而發現地震發生前期地磁變化，為地震預測帶來可能。

3.1 數據來源及預處理

本次實驗數據選取2008年1～6月份四川省成都市的地震地磁前兆數據。由于地磁前兆數據存在數據缺失、非標準化等問題，需要進行插值、整理、規范化處理等預處理操作，以保證數據的可用性。缺失數據采用均值法補齊，規范化采用z-score標準化方法進行完成，并去掉量綱。經過處理的數據符合標準正態分布，按照時間排列構成標準化的數據集。

3.2 計算結果與分析

首先以拐肘法、輪廓系數法和CH指標法分別計算K1、K2和K3值，如圖1所示。

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

根據聚類性能標準評估原理，從圖1（a）的拐肘法得到K1值為4，從圖1（b）的輪廓系數法得到K2值為5，從圖1（c）的CH指標法得到K3值為5。

其次求取平均值作為最優K值，K＝（4+5+5）／3＝4.7≈5。因此K＝5是最佳的聚類類別個數，此時可以達到最好的聚類效果，計算結果如圖2所示。

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

由圖2可見，存在5個聚類中心。即圖中大圓點，檢測到2個異常點，即圖中五角星，正好與2次地震對應。詳盡信息見表1。

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

作為對比，同樣以K＝5作為聚類數，利用傳統K-means算法對成都地震地磁前兆數據進行聚類分析，結果如圖3所示。

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

從圖3可知，傳統K-means聚類方法沒有找到異常離群點。實際上，在實驗過程中，K-means聚類方法因初始聚類中心選擇的隨機性問題，該方法難以發現離群點，并且每次計算結果差異較大。這也說明，同樣的K值下，K-means++聚類方法比Kmeans方法效果更好。

4 結束語

文中在介紹了K-means++算法基礎上，通過拐肘法、輪廓系數法和CH指標法聯合求取K值，優化了K值的確定方法，從而很好地解決了K-means++聚類算法中K值難以確定的問題。利用改進后的K-means++聚類算法對2008年1～6月份四川省成都市的地震地磁前兆數據進行聚類分析，發現了2個離群點，正好與發生的2次地震相對應，效果明顯優于傳統K-means結果。得到結論如下:

優化后的K-means++聚類算法優于傳統Kmeans算法，K-means++聚類算法可以較好地解決地震地磁的聚類問題。通過離群點檢索可以發現可能發生的地震，對地震監測預報工作具有重要的現實意義。