基于小波包散布熵與Meanshift 概率密度估計的軸承故障識別方法研究

張雄，張逸軒，張明，萬書亭?，何玉靈，豆龍江

（1.河北省電力機械裝備健康維護與失效預防重點實驗室，河北 保定 071003；2.華北電力大學 機械工程系，河北 保定 071003）

滾動軸承是旋轉機械中最常見、故障率最高的零部件之一，它的運行狀態關系到整個機械設備的可靠性和安全性，因此，軸承故障診斷方法是近年來工程測試和信號處理領域的熱點[1-3].振動信號中含有大量的軸承周期性沖擊信息，在軸承故障診斷中有著廣泛的應用[3-6].

軸承故障診斷一般分為兩步.第一步是故障特征的提取過程.這一過程的核心是如何準確地抑制振動信號中的干擾信息，準確地提取故障特征元素.在這一過程中，通常采用小波變換（Wavelet Transform[7]、小波包變換（Wavelet Packet Transform）[8]、經驗小波變換（Empirical Wavelet Transform）[9]、經驗模態分解（Empirical Mode Decomposition）[10]、集成經驗模態分解（Ensemble Empirical Mode Decomposition）[11]、局部均值分解（Local Mean Decomposition）[12]等處理手段對信號進行濾波和增維處理，目的是提取能更有效反映軸承故障信息的模態分量.通過對分解后和濾波后的分量的動態特性進行統計學計算，構造出能夠反映軸承振動信號的特征矩陣.其中信息熵、排列熵、模糊熵等動力學指標常被用來反映信號的瞬態特征.陳法法等[13]提出一種基于信息熵與優化最小二乘支持向量機的軸承性能退化趨勢模糊粒子預測方法，用于提升軸承性能退化指標預測精度.Zhang 等[14]通過計算局部迭代分解濾波后固有模態分量的多尺度排列熵，構造歸一化特征向量，對不同工況條件下的軸承故障進行識別.鄭近德等[15]采用復合多尺度模糊熵和迭代拉普拉斯得分對變分模態分解升維后的信號進行敏感特征選擇，以支持向量機對不同故障類型進行劃分.第二步是利用機器學習方法將特征集作為訓練樣本和測試樣本進行模式識別.該部分的核心問題包括聚類、分類、回歸和降維.Li 等[16]對比分析了模糊C 均值（Fuzzy C-Means，FCM）、Gustafson -Kessel 算 法、FN -DBSCAN 和FCMFP 算法各自特點.Yu 等[17]利用Gath-Geva（GG）聚類對故障特征進行分類，得到各軸承狀態的聚類中心和隸屬度矩陣，進行模式識別.

構造能夠充分反映信號樣本屬性且具有良好類內聚集性的特征矩陣，并尋求具有自適應能力和邊界特征的樣本分類方法是模式識別領域的核心問題.本文提出了一種基于小波包散布熵與Meanshift 概率密度估計軸承故障特征矩陣構造方法，通過計算樣本小波包各子帶的散布熵值，構建特征矩陣；進而利用PCA 對特征矩陣進行可視化降維，選取貢獻度最高的兩個主成分；最后采用Meanshift 概率密度估計聚類中心位置.通過實驗數據分析，驗證了該方法能有效識別不同類型的故障和不同程度的故障.

1 基本理論

1.1 小波包散布熵

小波包變換能同時連續分解信號的高頻分量和低頻分量，并能自適應地確定不同頻段的分辨率，大大提高了信號的時頻局部分析能力，得到了廣泛的應用.小波包變換過程可用式（1）表示.

式中：xi，j表示第i 層的第j 子帶信號（其中，i 是分解層數，j 是對應層的信號數）；K 為序列長度；Ln和Gn分別是小波包的低通濾波器和高通濾波器.

為了解決樣本熵計算時間長、實時性差、排列熵不考慮平均振幅與振幅之差等問題，Rostaghi 等[18]提出了一種新的時間序列不規則性度量指標，稱為散布熵（Dispersion Entropy，DE）.與樣本熵和排列熵（Permutation Entropy，PE）類似，散布熵也是一種表征時間序列不規則性的方法.散布熵值越大，不規則度越高；散布熵值越小，不規則度越低.

對于長度為N 的時間序列x={xj，j=1，2，…，N}，散布熵的計算步驟如下：

1）通過正態分布函數用于將時間序列映射到y={yj，j=1，2，…，N}.

式中：μ 和σ2分別表示序列的期望和方差.

2）通過線性變換將y 映射到[1，2，…，c]范圍.

式中：c 為類別個數；R 為取整函數.

3）計算嵌入向量：

式中：m 和d 表示嵌入維數和時延.

4）計算散布模式πv0v1…vm-1（v=1，2，…，c），如果，則πv0v1…vm-1為對應散布模式.

5）計算散布模式πv0v1…vm-1的概率：

式中：Number（πv0v1…vm-1）表示zm，ci在πv0v1…vm-1中的映射個數.

6）類比香農熵定義，將原信號的散布熵定義為：

當所有散布模式具有相同的概率（如噪聲信號）時，散布熵取最大值lncm.相反，當只有一個p（πv0v1…vm-1）值不等于零時（如周期信號），則表示時間序列是完全規則或可預測的數據，散布熵取最小值.

1.2 Meanshift 概率密度估計

Meanshift 聚類算法是一種無參數的聚類算法，能夠在根據樣本點計算數據概率密度分布區間.該算法已成功應用于圖像平滑、圖像分割和運動目標跟蹤等領域.

設Rd為d 維空間，x={xi}（i=1，2，…，n）為離散數據集合.Meanshift 可以定義為：

式中：Sh（x）={y:（y-x）T（y-x）≤h2}為球體區域；h為半徑.

向量Mh（x）對數據的概率密度梯度具有指向性.由于不同距離的點具有不同的權重系數，引入核函數K（x），概率密度函數f（x）表示為：

核函數定義為：

式中：o 為正則化系數，用來保證∫k（x）dx=1.

通過求偏導得到概率密度函數f（x）極值點.

式中：g（x）=-k′（x），相應的核函數為G（x）=og，dg（||x||2）.公式前半部分是以G（x）為核函數的概率密度估計的概率密度估計，后半部分為Meanshift 所指向的最大概率密度梯度的方向，可以表示為

Meanshift 算法本質上是一種自適應遞增迭代搜索數據分布概率密度分布梯度峰值的運算.迭代次數為t，搜索窗口（空間）為r，給定任意初始點x.迭代過程可以表述如下：

1）初始化t，r，設定閾值σ；

2）計算第t 次迭代的概率密度梯度mh（xt）；

3）更新搜索空間r，xt+1=xt+mh（xt）；

4）重復步驟2 和步驟3，直至mh（xt）≤σ.

采用仿真數據對上述過程進行說明.給定一組以一定概率分布在二維空間中的數據點.設定Meanshift 算法參數為r=0.5，σ=1×10-4.迭代過程如圖1 所示，對所設定的高維球區域內中心位置到離散數據點的向量進行加權處理，合成迭代向量梯度方向（類似于力的合成），然后，更新搜索窗口位置.Meanshift 算法在不預先設定分類數的情況下，可以自適應地沿著概率密度梯度方向迭代，并最終找到聚類中心的位置.

圖1 Meanshift 的迭代過程Fig.1 The iterative process of Meanshift

1.3 故障特征表征及模式識別過程

本文提出的軸承故障診斷方法流程如圖2 所示，具體步驟如下.

圖2 故障診斷方法流程圖Fig.2 Flow chart of fault diagnosis method

1）構建特征矩陣.選取訓練樣本形成信號集x=（x1，x2，…，xm），對原始信號集中的各個元素進行小波包分解，計算每個小波包子帶的散布熵構建特征矩陣.

2）采用主成分分析法對特征矩陣進行降維.將特征矩陣投影到二維空間，選擇貢獻率最高的兩個主成分構造二維特征矩陣（選擇兩個主成分（Principal Component，PC）可以顯示為二維圖，三個PC 可以顯示為三維圖，本文數據特征樣本以二維平面圖的形式顯示，選擇貢獻率最高的前兩個PC 分量構造特征矩陣）.

3）建立了估計模型.設定Meanshift 參數（本文搜索半徑r 的取值原則為在保障聚類種數的前提下，選擇盡可能小的窗口半徑），對主成分空間坐標點進行概率密度估計，得到聚類類別和聚類中心.

4）對測試樣本進行估計.對測試樣本重復上述步驟1 和步驟2，得到主成分特征矩陣，并計算其與訓練樣本的聚類中心的歐式距離，得到相應的隸屬關系.

2 實測信號分析

為了驗證該方法對軸承不同故障類型和故障程度診斷的有效性，分別采用CWRU 實驗室開源數據和QPZZ-II 旋轉機械故障模擬實驗臺數據進行分析.

2.1 CWRU 滾動軸承實驗數據分析（不同故障程度）

故障源數據為驅動端SKF6205 軸承經電火花加工在內圈生成的四類故障程度樣本，故障尺寸分別為0.007 英寸，0.014 英寸，0.021 英寸和0.028 英寸（本文選用數據為美國凱斯西儲大學實驗臺數據，原數據說明中使用單位為英寸，故本文使用單位為英寸.轉換為國際單位后，四類樣本故障尺寸分別是0.017 78 cm，0.035 56 cm，0.053 34 cm 和0.071 12 cm）.電機轉速為1 750 r/min，采樣頻率為12 kHz，軸承實驗臺模型如圖3 所示.

圖3 CWRU 軸承故障模擬實驗臺Fig.3 The structure of the experimental platform in CWRU

對四類不同故障程度的振動信號數據各取一組樣本，其時域波形如圖4 所示.

圖4 不同故障程度振動信號的時域波形Fig.4 The time domain waveform of vibration signals with different fault degrees

驗證散布熵相較于排列熵的穩定性以及對于不同故障程度具有較好的區分度.對四類不同故障程度的振動信號劃分成不同數據長度構造數據節點，節點1 數據長度為512，節點2 數據長度為1 024（512×2），節點3 數據長度為2 048（512×4），節點4對應數據長度為3 072（512×6），以此類推.分別計算四類故障程度振動信號10 個節點數據的散布熵，結果如圖5 所示.可以看出，不同故障程度下散布熵隨數據點長度增長的走勢大體相近且變化平緩，四種故障程度在各節點具有較好的區分度.

圖5 不同故障程度振動信號各節點數據的散布熵Fig.5 DE of each node of vibration signals with different fault degrees

計算上述各節點的排列熵作為對比，結果如圖6所示，可以看出，不同故障程度下排列熵隨數據長度增長的走勢振蕩明顯，且存在交叉，說明數據長度的選擇在較大程度上影響類間區分度.

圖6 不同故障程度振動信號各節點數據的排列熵Fig.6 PE of each node of vibration signals with different fault degrees

對四類不同故障程度振動信號各取40 組分析樣本，其中20 組為訓練樣本，20 組為測試樣本，采用本文所提故障識別方法進行處理.首先利用小波包分解對訓練樣本數據進行升維處理，然后計算每個樣本小波包各子帶的散布熵值，構建特征矩陣，進而利用PCA 對特征矩陣進行可視化降維，選取貢獻度最高的兩個主成分，最后采用Meanshift 概率密度估計聚類中心位置，結果如圖7 所示.

圖7 不同故障程度訓練樣本散布熵分布及Meanshift 聚類中心估計結果Fig.7 The results of DE distribution and Meanshift clustering center of training samples with different fault degrees

對20 組測試樣本進行分析，采用同樣的方法計算小波包散布熵構造特征矩陣，并通過PCA 進行可視化降維，然后計算測試樣本點與上述聚類中心的歸一化歐氏距離，結果如圖8 所示.歸一化歐氏距離越小，說明樣本與該聚類中心的隸屬度越高，可以看出，測試樣本被較清晰的劃分到四類故障程度類別中.

圖8 不同故障程度測試樣本散布熵與各聚類中心的歐氏距離Fig.8 The normalized Euclidean distance between DE of test samples with different fault degrees and clustering centers

采用EEMD 排列熵構造特征矩陣進行對比分析，通過PCA 可視化降維和Meanshift 概率密度估計后的訓練樣本分布和聚類中心位置如圖9 所示.可以看出，數據分布的類間距較小，類內聚集性較差.測試樣本與各聚類中心的歸一化歐氏距離如圖10 所示，可以看出，測試樣本1 和測試樣本2 出現較為嚴重的混疊，難以明確其隸屬關系.

圖9 不同故障程度訓練樣本排列熵分布及Meanshift 聚類中心估計結果Fig.9 The results of PE distribution and Meanshift clustering center of training samples with different fault degrees

圖10 不同故障程度測試樣本排列熵與各聚類中心的歐氏距離Fig.10 The normalized Euclidean distance between PE of test samples with different fault degrees and clustering centers

2.2 QPZZ-II 旋轉機械故障模擬實驗臺數據分析（不同故障類型）

為進一步驗證所提方法的有效性，采用QPZZ-II軸承故障模擬實驗臺（電機功率0.55 kW，調速范圍75～1 450 r/min）進行數據分析，故障軸承型號6205E（利用線切割分別在內圈、外圈及滾動體植入故障），軸承座位置水平方向和垂直方向布置振動加速度傳感器（型號：東華1A116E，量程：50 g），測試系統采用DH5922N 型動態信號采集分析儀（16 通道/256 kHz），采樣頻率為12 800 Hz，實驗臺結構圖如圖11所示.對三類不同故障類型的振動信號數據各取一組樣本，其時域波形如圖12 所示.

圖11 QPZZ-II 實驗臺結構Fig.11 The structure of the experimental platform on QPZZ-II

圖12 不同故障類型振動信號的時域波形Fig.12 The time domain waveform of vibration signals with different fault types

對三類不同故障類型振動信號各取40 組分析樣本，其中20 組為訓練樣本，20 組為測試樣本，采用本文所提故障識別方法進行處理.對20 組訓練樣本構造小波包散布熵特征矩陣，利用PCA 進行可視化降維，并用Meanshift 概率密度估計聚類中心，結果如圖13 所示.對20 組測試樣本進行分析，計算測試樣本點與上述聚類中心的歸一化歐氏距離，結果如圖14 所示.可以看出，測試樣本被較清晰的劃分到三類故障程度類別中.采用EEMD 排列熵構造特征矩陣進行對比分析，訓練樣本分布和聚類中心位置如圖15 所示，測試樣本與各聚類中心的歸一化歐氏距離如圖16 所示，可以看出，測試樣本1 和測試樣本2 出現較為嚴重的混疊.

圖13 不同故障類型訓練樣本散布熵分布及Meanshift 聚類中心估計結果Fig.13 The results of DE distribution and Meanshift clustering center of training samples with different fault types

圖14 測試樣本散布熵與各聚類中心的歐氏距離Fig.14 The normalized Euclidean distance between test samples DE and clustering centers

圖15 不同故障類型訓練樣本排列熵分布及Meanshift 聚類中心估計結果Fig.15 The results of PE distribution and Meanshift clustering center of training samples with different fault types

圖16 測試樣本排列熵與各聚類中心的歐氏距離Fig.16 The normalized Euclidean distance between test samples PE and clustering centers

3 結論

本文針對軸承故障模式識別領域的兩類典型問題（不同故障類型和不同故障程度數據樣本識別）展開研究，提出一種基于小波包散布熵和Meanshift 概率密度估計的軸承故障模式識別方法，通過CWRU和QPZZ-II 實驗臺數據分析驗證了所構造的小波包散布熵特征矩陣能夠充分反映信號樣本屬性且具有較好類內聚集性，同時Meanshift 無參概率密度估計具有良好的聚類邊界和數據樣本模式識別能力.具體而言：

1）散布熵隨數據點長度增長的走勢相較于排列熵變化平緩，各節點具有較好的區分度，說明散布熵對截取的不同長度信號樣本具有更好的穩定性和適應性.

2）訓練樣本的小波包散布熵經PCA 降維后相較于同樣處理的EEMD 排列熵具有更穩定的聚類區域以及更大的類間距離.

3）Meanshift 無參概率密度估計能夠通過迭代準確識別樣本特征的聚類中心，通過計算測試樣本散布熵坐標與各聚類中心的歐氏距離可以實現對測試樣本隸屬關系的判別.