楊 健，張永平

（鹽城工學(xué)院 信息工程學(xué)院，江蘇 鹽城 224051）

隨著現(xiàn)代工業(yè)技術(shù)的飛速發(fā)展，機器設(shè)備日益龐大復(fù)雜，機器故障特別是其中的軸承故障也日益增多，因此實時監(jiān)測軸承運行狀況并進行故障診斷變得越來越重要。基于無監(jiān)督學(xué)習(xí)的故障診斷方法是近年來故障診斷領(lǐng)域的一個研究熱點，Kmeans 算法是其中比較經(jīng)典的算法之一，由于該算法簡單、實用以及能快速收斂等特點，被運用于諸多領(lǐng)域。Kmeans 算法通常會隨機選擇k個點作為初始聚類中心，以此來確定初始劃分，如果初始聚類中心的選取不合理，很容易陷入局部最優(yōu)解［1］。因此，很多學(xué)者都在Kmeans 算法上作了一些改進，并取得了一系列成果。鄧海等［2］將密度法與最大最小原則相結(jié)合，對Kmeans的初始聚類中心點選擇進行了優(yōu)化，使得算法的準確率有所提高，但時間復(fù)雜度也比較高，耗費時間較長；趙慶［3］采用Canopy 算法粗聚類數(shù)據(jù)，克服了傳統(tǒng)Kmeans中心點隨機選取的盲目性，提升了模型的精確度，但Canopy 算法的初始閾值是人工設(shè)定的，導(dǎo)致聚類結(jié)果不太穩(wěn)定；劉紀偉等［4］用密度思想優(yōu)化了Kmeans初始中心點的選取，同時引入聚類有效性判別函數(shù)確定值，提高了算法的準確度，但也增加了算法的運行時間，執(zhí)行效率較低；李曉瑜等［5］用MapReduce 分布式框架并行化實現(xiàn)改進的Canopy-Kmeans 算法，具有很好的精確度和可擴展性，但人工設(shè)定Canopy 算法的初始閾值問題依然存在；李琪等［6］提出一種基于密度峰值的M-Canopy-Kmeans 算法，為Canopy-Kmeans算法初始中心點的隨機選取、算法受噪聲點影響等問題提供了解決思路；陳勝發(fā)等［7］提出基于密度加權(quán)的Canopy 改進K-medoids 算法，改善了算法的精準度；王海燕等［8］提出使用Canopy+算法用于實現(xiàn)對T1，T2的改進。

Canopy 算法的全局尋優(yōu)性具有很好的聚類效果，可以彌補傳統(tǒng)Kmeans算法隨機選擇初始聚類中心的缺陷。因此，本文在傳統(tǒng)Kmeans算法的基礎(chǔ)上，提出一種結(jié)合小波包分解和優(yōu)化Canopy-Kmeans 算法的滾軸故障檢測方法，并在西儲大學(xué)滾軸數(shù)據(jù)集上進行測試，結(jié)果表明該方法可有效提高滾軸故障檢測的準確性和穩(wěn)定性。

1 基本理論

1.1 小波包分解

小波包分解是一種改進小波變換缺點的算法。小波變換在低頻段具有較好的頻率分辨率，但在高頻段的頻率分辨率卻很差，所以信號的細節(jié)在小波變換時并不能很好地表示出來。小波包分解的基本思路是把信號分解為多個層次，即對各個頻段進行分解包括分解高頻部分的頻段，然后根據(jù)分解后的特征選取對應(yīng)的頻段。圖1為小波包三層分解結(jié)構(gòu)圖。

圖1 小波包三層分解結(jié)構(gòu)圖Fig. 1 Structure diagram of wavelet packet three-layer decomposition

1.2 Canopy-Kmeans算法

Kmeans 算法基本思路是：在給定的數(shù)據(jù)集合中隨機選擇k個數(shù)據(jù)作為k個聚類的初始中心點；然后按照距離最近的原則將數(shù)據(jù)集中其余數(shù)據(jù)分配給最近的一個聚類，直到數(shù)據(jù)集是空集為止；最后，通過各個聚類的樣本均值對聚類中心進行更新，直到聚類中心點沒有變化或達到終止條件后再輸出最終的聚類中心和k個聚類劃分。

Canopy-Kmeans 算法是在Kmeans 算法的基礎(chǔ)上加入Canopy 算法的一種優(yōu)化方法，流程如圖2 所示。在圖2 的Canopy-Kmeans 算法中，Canopy算法首先對數(shù)據(jù)集進行“粗”聚類，然后將得到的n個Canopy 子集的中心點作為Kmeans 算法的初始中心點，再用Kmeans 算法進行“細”聚類，最終得到聚類結(jié)果。

Canopy-Kmeans 算法雖然克服了Kmeans 算法中初始中心點隨機選擇和人工設(shè)定的問題，但也存在如下問題：Canopy 算法的初始閾值需要經(jīng)驗法或交叉驗證法等進行確定，使得聚類結(jié)果的穩(wěn)定性下降；同時，由于時間復(fù)雜度高，導(dǎo)致串行運算速度慢。

1.3 I-Canopy-Kmeans聚類算法

I-Canopy-Kmeans 算法是在Canopy-Kmeans算法的基礎(chǔ)上，對初始聚類中心的選取和閾值獲取方式方面進行了優(yōu)化，主要表現(xiàn)為：在初始聚類中心的選取方面，I-Canopy-Kmeans 算法是按照“最遠最近原則”進行的，即在得到n個Canopy時，任意兩個Canopy 中心點之間的距離都要盡可能地遠，且第n個Canopy 中心點應(yīng)該是其他數(shù)據(jù)點與前面n-1 個中心點最遠距離中最小的一個；在閾值獲取方式方面，主要是使用歐氏距離求出所有數(shù)據(jù)點的均值點，再計算均值點到所有數(shù)據(jù)點的距離，并用L1表示最遠距離、L2表示最近距離，然后將賦值給閾值T1、賦值給閾值T2。具體地說，I-Canopy-Kmeans 算法在對初始聚類中心的選取和閾值獲取方式方面進行優(yōu)化的主要流程如下：

步驟1 根據(jù)給定的數(shù)據(jù)集D，計算D中所有數(shù)據(jù)間的平均距離dM。

步驟2 將平均距離dM在數(shù)據(jù)集D中對應(yīng)的點稱為均值點，計算均值點和所有數(shù)據(jù)之間的歐氏距離dij，并將最遠的距離dij，max記為L1、最近的距離dij，min記為L2，同時將賦值給閾值T1、賦值給閾值T2。

步驟3 選取D中距離均值點最近的點c1，作為第1 個Canopy 的聚類中心，并將c1添加到中心點集合C中，即C=｛c1｝。

步驟4 計算D中所有數(shù)據(jù)與c1之間的歐氏距離dij。如果dij＜T2，就把它從數(shù)據(jù)集D中刪除；如果dij＜T1，就把它分配到c1所在的Canopy 中，同時從數(shù)據(jù)集D中刪除；如果dij≥T1，則比較各個dij大小，并將最大者對應(yīng)的點作為第2 個Canopy 聚類中心c2，同時將該點添加到中心點集合C，使C=｛c1，c2｝。

步驟5 計算數(shù)據(jù)集D中所有數(shù)據(jù)（點c1除外）與c2之間的歐氏距離dij。如果dij＜T2，就把它從數(shù)據(jù)集D中刪除；如果dij＜T1，就把它分配到c2所在的Canopy中，同時從數(shù)據(jù)集D中刪除。

步驟6 在數(shù)據(jù)集D中尋找距離c1最遠的點d1和距離c2最遠的點d2，從中選擇距離更近的點作為第3 個Canopy 聚類中心c3，并將該點添加到中心點集合C，此時C=｛c1，c2，c3｝。

步驟7 計算數(shù)據(jù)集D中所有數(shù)據(jù)（點c1、c2除外）與c3點之間歐氏距離dij。如果dij＜T2，就把它從數(shù)據(jù)集D中刪除；如果dij＜T1，就把它分配到c3所在的Canopy中，同時從數(shù)據(jù)集D中刪除。

步驟8 在數(shù)據(jù)集D中尋找距離c1最遠的點d1、距離c2最遠的點d2、距離c3最遠的點d3，從中選擇距離更近的點作為第4 個Canopy 聚類中心c4，并將該點添加到中心點集合C，此時C=｛c1，c2，c3，c4｝。

步驟9 計算數(shù)據(jù)集D中所有數(shù)據(jù)（點c1、c2、c3除外）與c4點之間歐氏距離dij。如果dij＜T2，就把它從數(shù)據(jù)集D中刪除；如果dij＜T1，就把它分配到c4所在的Canopy中，同時從數(shù)據(jù)集D中刪除。

步驟10 按照步驟8、步驟9 繼續(xù)在數(shù)據(jù)集D中尋找其他聚類中心點，直到數(shù)據(jù)集D中只剩下各個中心點為止。

步驟11 將生成的Canopy 數(shù)量賦值給K，Canopy 的中心點作為聚類初始中心點進行Kmeans聚類。

上述優(yōu)化流程中計算數(shù)據(jù)集中數(shù)據(jù)之間的歐氏距離如式（1）所示，計算數(shù)據(jù)集中所有數(shù)據(jù)間的平均距離如式（2）所示。

式中：dM為數(shù)據(jù)集X中所有樣本元素間的平均距離，其中X＝｛x1，x2，…，xn｝是一個包含n個樣本對象的數(shù)據(jù)集，每個樣本對象包含d維特征屬性；dij為數(shù)據(jù)集X中數(shù)據(jù)對象xi和xj之間的歐氏距離；xip（i＝1，2，…，n；p=1，2，…，d）表示第i個數(shù)據(jù)對象的第p維屬性，xjp（j＝1，2，…，n；p=1，2，…，d）表示第j個數(shù)據(jù)對象的第p維屬性。

2 I-Canopy-Kmeans 聚類算法的故障檢測方法

I-Canopy-Kmeans 聚類算法的滾軸故障檢測流程如圖3 所示，主要包括基于小波包分解的數(shù)據(jù)預(yù)處理、特征向量歸一化和I-Canopy-Kmeans算法的故障檢測模型訓(xùn)練3個方面。在具體的滾軸故障檢測中，先用三層小波包分解滾軸數(shù)據(jù)，求出滾動軸承特征向量；將特征向量歸一化以獲得新的樣本集，并將樣本集劃分為訓(xùn)練集和測試集；使用I-Canopy-Kmeans 算法對訓(xùn)練集中的數(shù)據(jù)進行訓(xùn)練，建立滾軸故障檢測模型后，再將測試集中的數(shù)據(jù)引入到滾軸故障檢測模型中，對模型進行有效性檢驗。

圖3 滾軸故障檢測流程Fig. 3 Flow of roller fault detection

2.1 基于小波包分解的數(shù)據(jù)預(yù)處理

一般來說，滾軸通常有4種運行狀態(tài)，分別為正常狀態(tài)、外圈故障、內(nèi)圈故障和滾動體故障。由于滾軸振動信號不同頻段的能量分布代表滾軸不同的運行狀態(tài)［9-10］，本文采用三層小波包分解對滾軸振動信號進行分析，得到8個小波包的分解頻段，再通過計算得到各頻段的能量及能量所占比例，進而根據(jù)能量分布構(gòu)建滾軸不同運行狀態(tài)的特征向量。

按照帕塞瓦爾（Parseval）定律，信號時間內(nèi)的總能量等于其頻率范圍內(nèi)的總能量。小波包分解僅改變了信號的形態(tài)，將信號的高頻成分與低頻成分分開，分解前后的總能量保持不變。若對信號x（t）進行小波包分解，得到n層頻帶信號，則各子頻帶信號的能量計算如式（3）所示。

式中：Enm表示第n層第m個頻段的能量，J；Sa，b表示二維系數(shù)矩陣，其中a是尺度參數(shù)，b是位移參數(shù)；xn，m(k)表示第n層第m個頻段中第k個小波包的系數(shù)值；x（t）為原始信號，t為時間變量，s。

總能量E為：

如以每個子頻段的能量百分比表示該信號的特征向量，則信號的特征向量T可以表示為：

通過小波包分解，獲得了多個特征向量，該特征向量包括8 個頻段的能量，這些特征向量將取代原有的龐大振動數(shù)據(jù)集，作為新的數(shù)據(jù)集來使用。

2.2 特征向量歸一化

為提高檢測模型的訓(xùn)練速度，需要將小波包分解后得到的特征向量T按式（6）進行歸一化處理。

式中：xm′(m= 1，2，…，8)是特征向量中第m個頻段歸一化后的結(jié)果；xm是特征向量中第m個頻段的能量值；xmax、xmin分別是特征向量中的最大值與最小值。

3 實驗過程及結(jié)果分析

3.1 振動數(shù)據(jù)特征提取

基于美國凱斯西儲大學(xué)提供的滾動軸承數(shù)據(jù)集［11］，采用I-Canopy-Kmeans 聚類算法的滾軸故障檢測方法，對滾軸振動數(shù)據(jù)使用三層小波包分解并進行特征提取，得到不同狀態(tài)下不同頻段的能量百分比，如圖4 所示。其中特征提取的部分數(shù)據(jù)如表1所示。

表1 特征提取的部分數(shù)據(jù)Table 1 Partial data of feature extraction

圖4 滾軸不同狀態(tài)下的能量譜Fig. 4 Energy spectrum of the rollers in different states

由圖4 可以看出，不同狀態(tài)下信號的能量分布是有差異的，其中正常狀態(tài)下的滾軸能量主要分布在第1、2 頻段，故障狀態(tài)下的能量主要分布在第4、5 頻段，且滾動體及外圈故障的能量分布在第4、5 頻段更加集中，滾軸內(nèi)圈故障的能量分布則稍微分散于1、2、3、6 頻段；在滾動體及外圈故障中，外圈故障在第6個頻段的能量分布比前3個頻段更多，而滾動體故障在前3 個頻段的能量分布比第6個頻段更多。

3.2 實驗過程及結(jié)果分析

基于表1 的特征提取數(shù)據(jù)，采用傳統(tǒng)的Kmeans 算法和I-Canopy-Kmeans 算法分別計算兩種算法下各自經(jīng)歷1 000 次實驗的滾軸平均故障檢測準確率，結(jié)果如表2 所示。由表2 可知，采用傳統(tǒng)Kmeans 算法得到的聚類效果評價指標中標準化互信息INM為0.709 3、調(diào)整互信息IAM為0.705 7、蘭德系數(shù)IR為0.849 1、調(diào)整蘭德系數(shù)IAR為0.553 8，采用I-Canopy-Kmeans 得到的INM為0.782 8、IAM為0.781 3、IR為0.9093、IAR為0.775 4。顯然，采用I-Canopy-Kmeans 檢測算法的各項指標得到明顯提高，其中，IAR提高最多，從0.553 8提高到0.775 4，提高了40.01%。

表2 算法評價指標Table 2 Algorithm evaluation metrics

經(jīng)過比較分析，發(fā)現(xiàn)兩種算法下得到的聚類中心對應(yīng)的特征向量均有變化，分別如表3、表4所示。表3、表4 說明傳統(tǒng)的Kmeans 算法可能得到了局部最優(yōu)解，I-Canopy-Kmeans 算法則更可能得到了全局最優(yōu)解。

表3 使用Kmeans得到的聚類中心對應(yīng)的特征向量Table 3 Feature vectors corresponding to the clustering centers obtained using Kmeans

表4 使用I-Canopy-Kmeans得到的聚類中心對應(yīng)的特征向量Table 4 Eigenvectors corresponding to clustering centers obtained using I-Canopy-Kmeans

4 結(jié)束語

I-Canopy-Kmeans 算法在Canopy-Kmeans 算法的基礎(chǔ)上，改進了初始聚類中心的選取問題，優(yōu)化了閾值的選取方式。實驗結(jié)果表明，與傳統(tǒng)的Kmeans 算法相比，I-Canopy-Kmeans 算法的各項評價指標均有提高，其中指標IAR提高最多，提高了40.01%。