基于小波包和聚類算法的滾動(dòng)軸承故障檢測研究

楊 健，張永平

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

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

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

1 基本理論

1.1 小波包分解

小波包分解是一種改進(jìn)小波變換缺點(diǎn)的算法。小波變換在低頻段具有較好的頻率分辨率，但在高頻段的頻率分辨率卻很差，所以信號(hào)的細(xì)節(jié)在小波變換時(shí)并不能很好地表示出來。小波包分解的基本思路是把信號(hào)分解為多個(gè)層次，即對(duì)各個(gè)頻段進(jìn)行分解包括分解高頻部分的頻段，然后根據(jù)分解后的特征選取對(duì)應(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ù)集合中隨機(jī)選擇k個(gè)數(shù)據(jù)作為k個(gè)聚類的初始中心點(diǎn)；然后按照距離最近的原則將數(shù)據(jù)集中其余數(shù)據(jù)分配給最近的一個(gè)聚類，直到數(shù)據(jù)集是空集為止；最后，通過各個(gè)聚類的樣本均值對(duì)聚類中心進(jìn)行更新，直到聚類中心點(diǎn)沒有變化或達(dá)到終止條件后再輸出最終的聚類中心和k個(gè)聚類劃分。

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

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

1.3 I-Canopy-Kmeans聚類算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

總能量E為：

如以每個(gè)子頻段的能量百分比表示該信號(hào)的特征向量，則信號(hào)的特征向量T可以表示為：

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

2.2 特征向量歸一化

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

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

3 實(shí)驗(yàn)過程及結(jié)果分析

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

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

表1 特征提取的部分?jǐn)?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)下信號(hào)的能量分布是有差異的，其中正常狀態(tài)下的滾軸能量主要分布在第1、2 頻段，故障狀態(tài)下的能量主要分布在第4、5 頻段，且滾動(dòng)體及外圈故障的能量分布在第4、5 頻段更加集中，滾軸內(nèi)圈故障的能量分布則稍微分散于1、2、3、6 頻段；在滾動(dòng)體及外圈故障中，外圈故障在第6個(gè)頻段的能量分布比前3個(gè)頻段更多，而滾動(dòng)體故障在前3 個(gè)頻段的能量分布比第6個(gè)頻段更多。

3.2 實(shí)驗(yàn)過程及結(jié)果分析

基于表1 的特征提取數(shù)據(jù)，采用傳統(tǒng)的Kmeans 算法和I-Canopy-Kmeans 算法分別計(jì)算兩種算法下各自經(jīng)歷1 000 次實(shí)驗(yàn)的滾軸平均故障檢測準(zhǔn)確率，結(jié)果如表2 所示。由表2 可知，采用傳統(tǒng)Kmeans 算法得到的聚類效果評(píng)價(jià)指標(biāo)中標(biāo)準(zhǔn)化互信息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 檢測算法的各項(xiàng)指標(biāo)得到明顯提高，其中，IAR提高最多，從0.553 8提高到0.775 4，提高了40.01%。

表2 算法評(píng)價(jià)指標(biāo)Table 2 Algorithm evaluation metrics

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

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

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

4 結(jié)束語

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