多尺度排列熵及其在滾動軸承故障診斷中的應(yīng)用

鄭近德 程軍圣 楊 宇

湖南大學(xué)汽車車身先進設(shè)計制造國家重點實驗室，長沙，410082

0 引言

時頻分析方法由于能夠提供振動信號時域和頻域局部信息而在故障診斷領(lǐng)域得到了廣泛應(yīng)用，是目前故障診斷的重要手段［1－2］。該方法通過對信號進行小波分析［3］或經(jīng)驗?zāi)B(tài)分解［4－5］，將非平穩(wěn)信號分解為若干個簡單的平穩(wěn)信號之和，然后對每個分量進行處理，提取時頻域信息，進而得到原始信號的完整時頻信息。然而，由于機械運轉(zhuǎn)過程中的摩擦、振動以及負載等因素，機械系統(tǒng)振動信號往往表現(xiàn)出非線性行為，采用時頻分析的方法，將信號分解為平穩(wěn)信號，難免有一定的局限性。非線性分析的方法可以不經(jīng)過對原始信號分解，而能夠直接提取隱藏在機械系統(tǒng)振動信號中其他方法無法提取的故障信息［6］。

近年來，越來越多的非線性分析方法被用于機械故障診斷，如關(guān)聯(lián)維數(shù)、近似熵、樣本熵、多尺度熵等。徐玉秀等［7］研究了旋轉(zhuǎn)機械的分形特征及故障診斷，Yan等［8］將近似熵應(yīng)用于機械系統(tǒng)健康狀態(tài)監(jiān)測，Zhang等［9］將多尺度熵應(yīng)用于滾動軸承的故障診斷。然而，分形維數(shù)的計算依賴數(shù)據(jù)長度，且較耗時，不適合在線監(jiān)測；近似熵相對一致性較差［9］；多尺度熵［10－11］是基于樣本熵而定義的，計算較耗時，且受時間序列的非平穩(wěn)性和異常值的影響。

排列熵［12］是一種新的隨機性和動力學(xué)突變的檢測方法，它具有計算簡單，抗噪能力強，且得到較穩(wěn)定的系統(tǒng)特征值所需時間序列短，適合在線監(jiān)測等優(yōu)點，在肌電信號處理［13］、心率信號處理［14］、氣溫復(fù)雜度［15］等方面都取得了良好的效果。Yan等［6］將其應(yīng)用于旋轉(zhuǎn)機械振動信號的特征提取，并將其與近似熵和Lempel－Ziv復(fù)雜度進行了對比，結(jié)果表明，排列熵能夠有效地檢測和放大振動信號的動態(tài)變化，并且能夠表征滾動軸承在不同狀態(tài)下的工況特征。然而，與傳統(tǒng)的基于單一尺度分析的非線性參數(shù)類似，排列熵只是檢測時間序列在單一尺度上的隨機性和動力學(xué)突變。Aziz等［16］提出了多尺度排列熵（multiscale permutation entropy，MPE）的概念，用于衡量時間序列在不同尺度下的復(fù)雜性和隨機性，并通過分析生理信號，將其與多尺度熵進行了對比，結(jié)果表明，相對于多尺度熵，MPE更具有魯棒性。由于機械系統(tǒng)比較復(fù)雜，振動信號不僅在單一尺度上包含有重要信息，而且在其他尺度上也包含有重要信息，因此，對振動信號進行多尺度分析是一種有效的方法。

本文將MPE引入到機械故障診斷領(lǐng)域，應(yīng)用于滾動軸承的故障特征的提取。由于正常滾動軸承的振動信號是隨機振動信號，而當滾動軸承發(fā)生故障時，振動信號隨機性和動力學(xué)行為會發(fā)生變化。因此，本文考慮用MPE來衡量振動信號的隨機性變化和動力學(xué)突變，并將MPE值作為特征參數(shù)，提取滾動軸承的故障特征。在此基礎(chǔ)上，結(jié)合支持向量機（support vector mechine，SVM）［17－18］作為模式識別分類，提出一種基于MPE和SVM的滾動軸承故障診斷方法，并將其應(yīng)用于滾動軸承實驗數(shù)據(jù)的分析，結(jié)果表明，新提出的方法能夠有效地診斷滾動軸承的故障類型。

1 排列熵及多尺度排列熵原理和算法

1.1 排列熵算法

排列熵的原理在于不考慮數(shù)據(jù)具體值，而是基于相鄰數(shù)據(jù)的對比。下面說明其計算方法。

考慮長度為N 的時間序列｛x（i），i＝1，2，…，N｝，對其進行相空間重構(gòu)，得到如下的時間序列：

式中，m為嵌入維數(shù)；λ為時延。

將X（i）的m個數(shù)據(jù)按照升序重新排列，即

如果存在x（i＋ （ji1－1）λ）＝x（i＋ （ji2－1）λ），此時按j值的大小來進行排序，即當jk1＜jk2，有x（i＋（ji1－1）λ）≤x（i＋（ji2－1）λ），所以，任意一個數(shù)據(jù)X（i）都可以得到一組符號序列：

其中，g ＝1，2，…，k，k ≤ m！，m 個不同的符號｛j1，j2，…，jm｝共有m！種不同的排列，對應(yīng)地，共有m！種不同的符號序列，S（g）是m！種符號序列中的一種。計算每一種符號序列出現(xiàn)的概率＝1，此時，時間序列｛x（i），i＝1，2，…，N｝的排列熵就可以按照Shannon熵的形式定義為

注意到，當Pg＝1／m！時，Hp（m）達到最大值ln（m?。虼?，可以通過ln（m！）將排列熵Hp（m）進行標準化處理，即

顯然，Hp的取值范圍是0≤Hp≤1。Hp值的大小表示時間序列的復(fù)雜和隨機程度。Hp越大，說明時間序列越隨機，反之，則說明時間序列越規(guī)則。Hp值的變化反映和放大了時間序列的局部細微變化。

1.2 排列熵參數(shù)的選取

在排列熵的計算中，有3個參數(shù)值需要考慮和設(shè)定，即時間序列長度N、嵌入維數(shù)m和時延λ。Bandt等［12］建議嵌入維數(shù)m取3～7，因為如果m等于1或2，此時重構(gòu)的序列中包含太少的狀態(tài)，算法失去意義和有效性，不能檢測時間序列的動力學(xué)突變。但是，如果m取值過大，也不合適，因為相空間的重構(gòu)將會均勻化時間序列，此時不僅計算比較耗時，而且也無法反映序列的細微變化［8，15］。

為研究數(shù)據(jù)長度N對排列熵值的影響，以長度分別為128、256、512、1024和2048的高斯白噪聲信號為例，求得對應(yīng)排列熵值，分別記為PE1～PE5，如圖1所示，它們在不同嵌入維數(shù)下差值如表1所示。

圖1 不同長度的高斯白噪聲的排列熵

表1 不同長度高斯白噪聲信號排列熵在不同嵌入維數(shù)下的差值

由圖1和表1可以發(fā)現(xiàn)，以嵌入維數(shù)m＝6為例，數(shù)據(jù)長度分別為1024和512時，熵值相差0.0659，而數(shù)據(jù)長度分別為2048和1024時，則熵值僅相差0.0309，因此，此時選數(shù)據(jù)長度為1024較合適。而對m＝5而言，數(shù)據(jù)長度分別為1024和256的信號的熵值僅相差0.0502，此時，數(shù)據(jù)長度為256已經(jīng)可以估計合理的排列熵值。一般，嵌入維數(shù)較小時，數(shù)據(jù)長度則要求越小。

時延λ對時間序列的計算影響較小，以長度為512的高斯白噪聲信號為例，在不同λ下的排列熵值隨嵌入維數(shù)的變化關(guān)系如圖2所示，由圖可以看出，時延對信號熵值的影響較小，因此，本文取λ＝1。

圖2 高斯白噪聲信號在不同時延下的排列熵

1.3 多尺度排列熵定義

多尺度排列熵定義為不同尺度下的排列熵，計算方法如下：

（1）考慮時間序列｛x（i），i＝1，2，…，N｝，對其進行粗粒化處理，得到粗?；蛄械谋磉_式為

其中，［N／τ］表示對N／τ取整；τ為尺度因子，τ＝1，2，…。顯然τ＝1，粗?；蛄屑礊樵夹蛄?；τ＞1時，原始序列被粗粒化為長度為［N／τ］的粗粒序列。

（2）計算每個粗粒序列的排列熵，并畫成尺度因子的函數(shù)，上述過程即稱為多尺度排列熵分析。

尺度因子的最大值一般取大于10即可，但要保證粗?；蛄虚L度［N／τ］不影響熵值的計算。

為了選取合適的計算MPE的嵌入維數(shù)，仍以高斯白噪聲為例，數(shù)據(jù)長度為2048，尺度因子最大值為12，λ＝1，在m 分別為4、5、6和7時，求得它們的MPE，相對耗時分別為0.1880s、0.6710s、3.8290s和27.6710s，將MPE畫成尺度因子的函數(shù)，如圖3所示。

由圖3可以看出，若m 取值太小，則PE值隨尺度因子的增大而減小，但m越大，計算越耗時，因此，本文選取m＝6。此外，由圖3可以看出，高斯白噪聲的MPE隨著尺度因子的增大而單調(diào)遞減，這說明白噪聲只在最小尺度上包含有主要信息。

圖3 高斯白噪聲在不同嵌入維數(shù)下的MPE

2 基于MFE的滾動軸承故障診斷方法及應(yīng)用

在上述理論的基礎(chǔ)上，本文提出基于MPE和SVM的滾動軸承故障診斷方法。首先，從滾動軸承的原始振動信號從提取MPE；其次，依據(jù)MPE提取合適的故障特征向量；第三，采用SVM進行故障分類，從而實現(xiàn)滾動軸承故障類別的診斷。

為了說明本文方法的有效性，本文將該方法應(yīng)用于實驗數(shù)據(jù)分析。本文實驗數(shù)據(jù)采用Case Western Reserve University（CWRU）軸承數(shù)據(jù)中心提供的滾動軸承試驗數(shù)據(jù)［9］。測試軸承為6205－2RS JEM SKF深溝球軸承，電機功率約為2206.4963W，轉(zhuǎn)速為1730r／min，采用電火花加工技術(shù)在軸承上布置單點故障，故障直徑為0.5334mm，深度為0.2794mm。在此情況下采集到正常（normal，簡稱 NORM）、內(nèi)圈故障（inner race fault，IRF）、外圈故障（outer race fault，ORF）和滾動體故障（rolling element fault，REF）4種狀態(tài)的振動信號，各30組數(shù)據(jù)，數(shù)據(jù)長度為2048，采樣頻率為12kHz，4種狀態(tài)軸承的振動信號時域波形如圖4所示。

圖4 正常和不同故障軸承振動信號的時域波形

從圖4不易發(fā)現(xiàn)正常和故障軸承振動信號的明顯區(qū)別，尤其是正常和滾動體故障，以及內(nèi)圈故障和外圈故障。因此，本文首先對振動信號進行MPE分析，取嵌入維數(shù)m＝6，時延λ＝1，最大尺度因子為12，4種狀態(tài)滾動軸承的多尺度排列熵畫成尺度因子的函數(shù)，如圖5所示。

圖5 正常和故障滾動軸承振動信號的多尺度排列熵

在尺度因子等于1時，即為原始振動信號的排列熵，由于熵值比較接近，無法明顯地區(qū)別三種故障和正常軸承的類型，因此有必要對振動信號進行多尺度分析。為此，以多尺度排列熵值為特征參數(shù)，同時建立基于SVM的分類器，進行訓(xùn)練和測試。如果采用全部的12個特征值進行訓(xùn)練，會造成信息的冗余，且訓(xùn)練比較耗時，也需要較多的訓(xùn)練樣本，且由圖5也可以看出，前幾個尺度的熵值表征了振動信號的主要信息，因此，采用前4個尺度的排列熵值作為特征向量，即T＝（PE1，PE2，PE3，PE4）。因此，本文的方法如下：

首先，提取特征參數(shù)，即對振動信號進行MPE分析，提取特征參數(shù)T。正常、滾動體故障、內(nèi)圈故障和外圈故障4種狀態(tài)，每種狀態(tài)取30個樣本，故每種狀態(tài)可得到30個表征故障特征的特征向量，共得到120個特征向量。

其次，訓(xùn)練分類器。由于有3種故障狀態(tài)和正常狀態(tài)，因此，需建立3個SVM，其中SVM1為正常對3種故障分類器，SVM2為內(nèi)圈故障對滾動體和外圈故障分類器，SVM3為滾動體故障和外圈故障分類器。每種狀態(tài)隨機抽取10個樣本進行訓(xùn)練，并將每組30個樣本用來測試。經(jīng)過訓(xùn)練，SVM1和SVM3采用徑向基核函數(shù)，SVM2采用多項式核函數(shù)。基于SVM的多故障分類器如圖6所示。

圖6 多類故障支持向量機分類器示意圖

最后，測試分類器。對已訓(xùn)練的SVM1、SVM2和SVM3，用全部樣本進行測試，詳細測試樣本輸出結(jié)果如表2所示。

表2 測試樣本SVM分類器的輸出結(jié)果

由表2可以看出，本文提出的方法有很好的效果，在全部樣本用來測試中，只有一組外圈故障的樣本被錯分為滾動體故障，其他都得到了正確的分類，正確識別率為99.17%。為了比較，下面建立以 BP 神 經(jīng) 網(wǎng) 絡(luò) 為 基 礎(chǔ) 的 多 分 類 器［9，19－20］，BP分類器除輸入層外，第一層隱含層有8個節(jié)點，第二層輸出層有4個節(jié)點。為表述方便，標記正常為1類，內(nèi)圈故障為2類，滾動體為3類，外圈故障為4類。BP分類器的訓(xùn)練和測試樣本與支持向量機相同，其分類結(jié)果如圖7所示。

圖7 BP神經(jīng)網(wǎng)絡(luò)分類器輸出結(jié)果

圖7中，訓(xùn)練樣本作測試時，全部分類正確，而測試樣本作測試時有3組分類錯誤，準確率為97.5%。這說明支持向量機的分類效果要優(yōu)于BP神經(jīng)網(wǎng)絡(luò)。而且，在訓(xùn)練時間上，支持向量機也比BP神經(jīng)網(wǎng)絡(luò)短得多。

為了說明進行多尺度分析的必要性，下面選取尺度因子等于1時（即原始信號）的排列熵值作為特征參數(shù)，分別通過SVM和BP神經(jīng)網(wǎng)絡(luò)分類器進行訓(xùn)練和測試。其中，SVM分類器中，有6組樣本分類錯誤，如表3所示；而BP神經(jīng)網(wǎng)絡(luò)分類器也有6組樣本分類錯誤，如圖8所示。

表3 測試樣本SVM分類器的輸出結(jié)果

圖8 BP神經(jīng)網(wǎng)絡(luò)分類器輸出結(jié)果

從上文可以看出，原始信號單一尺度的排列熵作為特征參數(shù)，分類效果不理想，這說明了單一尺度上原始振動信號的排列熵并不能反映故障的本質(zhì)，而多個尺度上的排列熵值則能夠更好地實現(xiàn)分類。另外易發(fā)現(xiàn)，特征值的選取對分類結(jié)果的影響尤為關(guān)鍵。本文選取特征值為前4個尺度上的特征值，主要基于以下原因考慮：如果特征值過少，不能完全反映故障的特征信息，而特征值過多會造成信息冗余，且需要增加訓(xùn)練樣本和訓(xùn)練時間，因此，本文選取了前4個尺度上的特征值。文獻［9］中選擇多尺度熵值的統(tǒng)計量時，將最大值、最小值、代數(shù)平均、幾何平均和標準差作為特征向量，但統(tǒng)計量方法忽略了特征值之間的內(nèi)在關(guān)系，因此，采用前4個尺度因子的排列熵值作為特征參數(shù)。

3 結(jié)束語

機械系統(tǒng)發(fā)生故障時，振動信號會在不同尺度上表現(xiàn)出不同程度的隨機性和動力學(xué)突變，基于此，本文提出了一種新的基于多尺度排列熵和支持向量機的滾動軸承故障診斷方法，并將支持向量機與BP神經(jīng)網(wǎng)絡(luò)分類效果進行了對比。結(jié)果表明，支持向量機在訓(xùn)練時間和準確率方面，都優(yōu)于BP神經(jīng)網(wǎng)絡(luò)。此外，本文還將特征向量包含多個尺度上熵值與特征向量僅包含原始信號單一尺度上排列熵值進行了對比，結(jié)果表明，原始信號單一尺度上排列熵值的不能全部反映故障的本質(zhì)，而多尺度的排列熵值則有很好的診斷效果。本文提出的方法為故障診斷提供了一種新的思路和手段。

［1］于德介，程軍圣，楊宇.機械故障診斷的 Hilbert－Huang變換方法［M］.北京：科學(xué)出版社，2007.

［2］何正嘉，陳進，王太勇，等.機械故障診斷理論及應(yīng)用［M］.北京：高等教育出版社，2010.

［3］程軍圣，于德介，鄧乾旺，等.連續(xù)小波變換在滾動軸承故障診斷中的應(yīng)用［J］.中國機械工程，2003，14（23）：2037－2040.Cheng Junsheng，Yu Dejie，Deng Qianwang，et al.Rolling Bearing Fault Diagnosis Using Continuous Wavelet Transform［J］.China Mechanical Engineering，2003，14（23）：2037－2040.

［4］Huang N E，Wu Z.A Review on Hilbert－Huang Transform：Method and Its Applications to Geophysical Studies［J］.Advances in Adaptive Data Analysis，2009，1：1－23.

［5］Yu Dejie，Cheng Junsheng，Yang Yu.Application of EMD Method and Hilbert Spectrum to the Fault Diagnosis of Roller Bearings［J］.Mechanical Systems and Signal Processing，2005，19：259－270.

［6］Yan Ruqiang，Liu Yongbin，Gao R X.Permutation Entropy：A Nonlinear Statistical Measure for Status Characterization of Rotary Machines［J］.Mechanical Systems and Signal Processing，2012，29：474－484.

［7］徐玉秀，鐘建軍，聞邦椿.旋轉(zhuǎn)機械動態(tài)特性的分形特征及故障診斷［J］.機械工程學(xué)報，2005，41（12）：186－189.Xu Yuxiu，Zhong Jianjun，Wen Bangchun.Fractal Fault Diagnosis and Classification to Modal Characteristic of Rotor System［J］.Journal of Chinese Mechanical Engineering，2005，41（12）：186－189.

［8］Yan Ruqiang，Gao R X.Approximate Entropy as a Diagnostic Tool for Machine Health Monitoring［J］.Mech.Syst.Signal Process，2007，21：824－839.

［9］Zhang Long，Xiong Guoliang，Liu Hesheng.Bearing Fault Diagnosis Using Multi－scale Entropy and Adaptive Neuro－fuzzy Inference［J］.Expert Systems with Applications，2010，37：6077－6085.

［10］Richman J S，Moorman J R.Physiological Timeseries Analysis Using Approximate Entropy and Sample Entropy［J］.American Journal of Physiology－Heart and Circulatory Physiology，2000，278：2039－2049.

［11］Costa M，Goldberger A L，Peng C K.Multiscale Entropy Analysis of Physiologic Time Series［J］.Physical Review Letters，The American Physiological Society，2002：068102（1－4）.

［12］Bandt C，Pompe B.Permutation Entropy：a Natural Complexity Measure for Time Series［J］.Physical Review Letters，The American Physiological Society，2002：174102（1－4）.

［13］袁明，羅志增.基于排列組合熵的表面肌電信號特征分析［J］.杭州電子科技大學(xué)學(xué)報，2012，32（1）：64－67.Yuan Ming，Luo Zhizeng.Feature Analysis of SEMG Based on Permutation Entropy［J］.Journalof Hangzhou Dianzi University，2012，32（1）：64－67.

［14］馬千里，卞春華.改進排列熵方法及其在心率變異復(fù)雜度分析中的應(yīng)用［J］.中國組織工程研究與臨床康復(fù)，2010，52（14）：9781－9785.Ma Qianli，Bian Chunhua.Application of Modified Permutation Entropy in Heart Rate Variability Analysis［J］.Journal of Clinical Rehabilitative Tissue Engineering Research，2010，52（14）：9781－9785.

［15］侯威，封國林，董文杰，等.利用排列熵檢測近40年華北地區(qū)氣溫突變的研究［J］.物理學(xué)報，2006，55（55）：2663－2668.Hou Wei，F(xiàn)eng Guolin，Dong Wenjie，et al.A Technique for Distinguishing Dynamical Species in the Temperature Time Series of North China［J］.Acta Physica Sinica，2006，55（55）：2663－2668.

［16］Aziz W，Arif M.Multiscale Permutation Entropy of Physiological Time Series［C］／／Proceeding of IEEE International Multi－topic Conference，INMIC，2005.

［17］李岳，陶利民，溫熙森.用于滾動軸承故障檢測與分類的支持向量機方法［J］.中國機械工程，2005，16（6）：498－501.Li Yue，Tao Limin，Wen Xisen.Support Vector Machines Based Approach for Ball Bearing Fault Detection and Classification［J］.China Mechanical Engineering，2005，16（6）：498－501.

［18］Yang Yu，Yu Dejie，Cheng Junsheng.A Fault Diagnosis Approach for Roller Bearing Based on IMF Envelope Spectrum and SVM［J］.Measurement，2007，40：943－950.

［19］李永強，劉杰，侯祥林，等.人工神經(jīng)網(wǎng)絡(luò)的混合算法及其工程應(yīng)用［J］.機械工程學(xué)報，2004，40（1）：127－130.Li Yongqiang，Liu Jie，Hou Xianglin，et al.Mixed Method of Artificial Neural Network and Its Application on Fault Diagnosis for Rotational Machine［J］.Journal of Chinese Mechanical Engineering，2004，40（1）：127－130.

［20］飛思科技產(chǎn)品研發(fā)中心.神經(jīng)網(wǎng)絡(luò)理論與 MATLAB7實現(xiàn)［M］.北京：電子工業(yè)出版社，2005.