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

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

（1.河北省電力機(jī)械裝備健康維護(hù)與失效預(yù)防重點(diǎn)實(shí)驗(yàn)室，河北 保定 071003；2.華北電力大學(xué) 機(jī)械工程系，河北 保定 071003）

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

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

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

1 基本理論

1.1 小波包散布熵

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

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

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

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

1）通過正態(tài)分布函數(shù)用于將時(shí)間序列映射到y(tǒng)={yj，j=1，2，…，N}.

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

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

式中：c 為類別個(gè)數(shù)；R 為取整函數(shù).

3）計(jì)算嵌入向量：

式中：m 和d 表示嵌入維數(shù)和時(shí)延.

4）計(jì)算散布模式πv0v1…vm-1（v=1，2，…，c），如果，則πv0v1…vm-1為對應(yīng)散布模式.

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

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

6）類比香農(nóng)熵定義，將原信號的散布熵定義為：

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

1.2 Meanshift 概率密度估計(jì)

Meanshift 聚類算法是一種無參數(shù)的聚類算法，能夠在根據(jù)樣本點(diǎn)計(jì)算數(shù)據(jù)概率密度分布區(qū)間.該算法已成功應(yīng)用于圖像平滑、圖像分割和運(yùn)動目標(biāo)跟蹤等領(lǐng)域.

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

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

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

核函數(shù)定義為：

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

通過求偏導(dǎo)得到概率密度函數(shù)f（x）極值點(diǎn).

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

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

1）初始化t，r，設(shè)定閾值σ；

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

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

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

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

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

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

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

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

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

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

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

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

2 實(shí)測信號分析

為了驗(yàn)證該方法對軸承不同故障類型和故障程度診斷的有效性，分別采用CWRU 實(shí)驗(yàn)室開源數(shù)據(jù)和QPZZ-II 旋轉(zhuǎn)機(jī)械故障模擬實(shí)驗(yàn)臺數(shù)據(jù)進(jìn)行分析.

2.1 CWRU 滾動軸承實(shí)驗(yàn)數(shù)據(jù)分析（不同故障程度）

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

圖3 CWRU 軸承故障模擬實(shí)驗(yàn)臺Fig.3 The structure of the experimental platform in CWRU

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

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

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

圖5 不同故障程度振動信號各節(jié)點(diǎn)數(shù)據(jù)的散布熵Fig.5 DE of each node of vibration signals with different fault degrees

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

圖6 不同故障程度振動信號各節(jié)點(diǎn)數(shù)據(jù)的排列熵Fig.6 PE of each node of vibration signals with different fault degrees

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

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

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

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

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

圖9 不同故障程度訓(xùn)練樣本排列熵分布及Meanshift 聚類中心估計(jì)結(jié)果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 旋轉(zhuǎn)機(jī)械故障模擬實(shí)驗(yàn)臺數(shù)據(jù)分析（不同故障類型）

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

圖11 QPZZ-II 實(shí)驗(yàn)臺結(jié)構(gòu)Fig.11 The structure of the experimental platform on QPZZ-II

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

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

圖13 不同故障類型訓(xùn)練樣本散布熵分布及Meanshift 聚類中心估計(jì)結(jié)果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 不同故障類型訓(xùn)練樣本排列熵分布及Meanshift 聚類中心估計(jì)結(jié)果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 結(jié)論

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

1）散布熵隨數(shù)據(jù)點(diǎn)長度增長的走勢相較于排列熵變化平緩，各節(jié)點(diǎn)具有較好的區(qū)分度，說明散布熵對截取的不同長度信號樣本具有更好的穩(wěn)定性和適應(yīng)性.

2）訓(xùn)練樣本的小波包散布熵經(jīng)PCA 降維后相較于同樣處理的EEMD 排列熵具有更穩(wěn)定的聚類區(qū)域以及更大的類間距離.

3）Meanshift 無參概率密度估計(jì)能夠通過迭代準(zhǔn)確識別樣本特征的聚類中心，通過計(jì)算測試樣本散布熵坐標(biāo)與各聚類中心的歐氏距離可以實(shí)現(xiàn)對測試樣本隸屬關(guān)系的判別.