張雄，張逸軒，張明，萬(wàn)書亭?，何玉靈，豆龍江

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

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

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

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

1 基本理論

1.1 小波包散布熵

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

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

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

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

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

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

2）通過(guò)線性變換將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為對(duì)應(yīng)散布模式.

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

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

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

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

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

Meanshift 聚類算法是一種無(wú)參數(shù)的聚類算法，能夠在根據(jù)樣本點(diǎn)計(jì)算數(shù)據(jù)概率密度分布區(qū)間.該算法已成功應(yīng)用于圖像平滑、圖像分割和運(yùn)動(dòng)目標(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）對(duì)數(shù)據(jù)的概率密度梯度具有指向性.由于不同距離的點(diǎn)具有不同的權(quán)重系數(shù)，引入核函數(shù)K（x），概率密度函數(shù)f（x）表示為：

核函數(shù)定義為：

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

通過(guò)求偏導(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.迭代過(guò)程可以表述如下：

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

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

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

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

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

圖1 Meanshift 的迭代過(guò)程Fig.1 The iterative process of Meanshift

1.3 故障特征表征及模式識(shí)別過(guò)程

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

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

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

2）采用主成分分析法對(duì)特征矩陣進(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ù)的前提下，選擇盡可能小的窗口半徑），對(duì)主成分空間坐標(biāo)點(diǎn)進(jìn)行概率密度估計(jì)，得到聚類類別和聚類中心.

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

2 實(shí)測(cè)信號(hào)分析

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

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

故障源數(shù)據(jù)為驅(qū)動(dòng)端SKF6205 軸承經(jīng)電火花加工在內(nèi)圈生成的四類故障程度樣本，故障尺寸分別為0.007 英寸，0.014 英寸，0.021 英寸和0.028 英寸（本文選用數(shù)據(jù)為美國(guó)凱斯西儲(chǔ)大學(xué)實(shí)驗(yàn)臺(tái)數(shù)據(jù)，原數(shù)據(jù)說(shuō)明中使用單位為英寸，故本文使用單位為英寸.轉(zhuǎn)換為國(guó)際單位后，四類樣本故障尺寸分別是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)臺(tái)模型如圖3 所示.

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

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

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

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

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

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

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

對(duì)四類不同故障程度振動(dòng)信號(hào)各取40 組分析樣本，其中20 組為訓(xùn)練樣本，20 組為測(cè)試樣本，采用本文所提故障識(shí)別方法進(jìn)行處理.首先利用小波包分解對(duì)訓(xùn)練樣本數(shù)據(jù)進(jìn)行升維處理，然后計(jì)算每個(gè)樣本小波包各子帶的散布熵值，構(gòu)建特征矩陣，進(jìn)而利用PCA 對(duì)特征矩陣進(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

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

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

采用EEMD 排列熵構(gòu)造特征矩陣進(jìn)行對(duì)比分析，通過(guò)PCA 可視化降維和Meanshift 概率密度估計(jì)后的訓(xùn)練樣本分布和聚類中心位置如圖9 所示.可以看出，數(shù)據(jù)分布的類間距較小，類內(nèi)聚集性較差.測(cè)試樣本與各聚類中心的歸一化歐氏距離如圖10 所示，可以看出，測(cè)試樣本1 和測(cè)試樣本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 不同故障程度測(cè)試樣本排列熵與各聚類中心的歐氏距離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)臺(tái)數(shù)據(jù)分析（不同故障類型）

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

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

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

對(duì)三類不同故障類型振動(dòng)信號(hào)各取40 組分析樣本，其中20 組為訓(xùn)練樣本，20 組為測(cè)試樣本，采用本文所提故障識(shí)別方法進(jìn)行處理.對(duì)20 組訓(xùn)練樣本構(gòu)造小波包散布熵特征矩陣，利用PCA 進(jìn)行可視化降維，并用Meanshift 概率密度估計(jì)聚類中心，結(jié)果如圖13 所示.對(duì)20 組測(cè)試樣本進(jìn)行分析，計(jì)算測(cè)試樣本點(diǎn)與上述聚類中心的歸一化歐氏距離，結(jié)果如圖14 所示.可以看出，測(cè)試樣本被較清晰的劃分到三類故障程度類別中.采用EEMD 排列熵構(gòu)造特征矩陣進(jìn)行對(duì)比分析，訓(xùn)練樣本分布和聚類中心位置如圖15 所示，測(cè)試樣本與各聚類中心的歸一化歐氏距離如圖16 所示，可以看出，測(cè)試樣本1 和測(cè)試樣本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 測(cè)試樣本散布熵與各聚類中心的歐氏距離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 測(cè)試樣本排列熵與各聚類中心的歐氏距離Fig.16 The normalized Euclidean distance between test samples PE and clustering centers

3 結(jié)論

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

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

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

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