基于多尺度排列熵和線性局部切空間排列的機(jī)械故障診斷特征提取

趙建崗， 寧 靜,2， 寧云志， 陳春俊， 李艷萍

(1. 西南交通大學(xué) 機(jī)械工程學(xué)院， 成都 610031； 2. 軌道交通運(yùn)維技術(shù)與裝備四川省重點(diǎn)實(shí)驗(yàn)室， 成都 610031)

過去的三十年中，基于信號處理方法的特征提取算法在機(jī)械領(lǐng)域被廣泛研究[1]，在文獻(xiàn)[2]中，用近似熵(approximate entropy,ApEn)監(jiān)測滾動軸承的健康狀況。然而ApEn嚴(yán)重依賴于數(shù)據(jù)長度，并且其估計(jì)值總是低于處理短數(shù)據(jù)集時的預(yù)期值。為了克服這種不足，提出了樣本熵(sample entropy,SampEn)[3]，由于SampEn值對數(shù)據(jù)長度不敏感并且不受信號中的噪聲影響，因此它受到了很多關(guān)注。然而ApEn和SampEn僅在單一尺度上估計(jì)復(fù)雜性，當(dāng)他們應(yīng)用于多時間尺度時效果極差。為克服這個弱點(diǎn)，在文獻(xiàn)[4]中提出了多尺度熵(multiscale entropy, MSE)。但MSE對故障診斷復(fù)雜性的估計(jì)不準(zhǔn)確。另外用它處理長時間序列時特別耗時[5]。

近年來，引入排列熵(permutation entropy,PE)[6]來計(jì)算機(jī)械系統(tǒng)的復(fù)雜性[7-8]。但是與AnEn和SampEn一樣，PE只能通過單一尺度分析時間序列，即可能會丟失嵌入其他尺度的許多有用的故障信息。但是某些地區(qū)機(jī)車和機(jī)床的操作環(huán)境也非?？量毯蛷?fù)雜[9-10]，監(jiān)控系統(tǒng)收集的大量信號通常是非靜態(tài)和非線性信號。由于機(jī)械組件之間的相互作用和耦合效應(yīng)，這些獲得的信號包含多種自然振動模式，這導(dǎo)致單尺度特征提取方法不能理想地表達(dá)這些信號[11]。為克服這一缺點(diǎn)，提出了多尺度排列熵(multiscale permutation entropy,MPE)[12]來估算不同尺度上時間序列的復(fù)雜性。在文獻(xiàn)[13]中MPE用于提取滾動軸承的故障特征，并且驗(yàn)證了MPE的性能優(yōu)于PE。

通常在使用MPE進(jìn)行特征提取之后，特征矩陣的維度將非常高。提高故障識別的準(zhǔn)確性，需要降低高維特征矩陣的維數(shù)以獲得低維主要特征，即降維[14]。主成分分析(principal component analysis,PCA)[15]是經(jīng)典的降維算法，但其僅對線性和高斯數(shù)據(jù)有效，并且在處理非線性或非平穩(wěn)信號時經(jīng)常不能達(dá)到理想的效果。最近流形學(xué)習(xí)被用來提取嵌入在非線性高維數(shù)據(jù)中的低維內(nèi)在結(jié)構(gòu)，而且它已被用于機(jī)械故障診斷[16-17]。其中以線性局部切線空間排列(linear local tangent space alignment,LLTSA)[18]性能較好，本文采用LLTSA進(jìn)一步提取嵌入在由MPE計(jì)算的高維特征矩陣中的低維主要特征。

丁吉等[19]用同步壓縮小波變換對發(fā)電機(jī)的連桿軸承進(jìn)行時頻分析，為發(fā)電機(jī)相關(guān)性能和工況的保障提供了重要的估計(jì)方法。孫自強(qiáng)等[20]用小波分形方法計(jì)算了風(fēng)機(jī)組軸承時域波形的關(guān)鍵維數(shù)，可以確定軸承的故障形式。向丹[21]等利用LLTSA對特征進(jìn)行降維，并通過SVM對故障類別進(jìn)行辨識。蘇祖強(qiáng)等[22]通過LLTSA與加權(quán)k最近鄰分類器提高了故障分類的精度。熊慶等[23]研究了多重分形趨勢波動分析特征參數(shù)的敏感性和穩(wěn)定性，并成功運(yùn)用于滾動軸承的定量故障診斷。李永健等[24]通過改進(jìn)的多尺度排列熵對列車的軸箱軸承進(jìn)行了診斷研究，取得了較好的效果。

本文將MPE算法在多尺度信息處理上的優(yōu)越性和LLTSA在非線性數(shù)據(jù)強(qiáng)適應(yīng)性相結(jié)合應(yīng)用于軸承的故障診斷，用兩個軸承箱數(shù)據(jù)對MPE-LLTSA方法進(jìn)行了測試，測試結(jié)果表明了該方法在機(jī)械模式分類和故障識別領(lǐng)域的應(yīng)用潛力。

1 理論背景

1.1 MPE的簡介

1.1.1 MPE算法

多尺度排列熵(MPE)被定義為不同尺度時間序列的一組PE值。MPE算法包括兩個步驟：① 應(yīng)用粗粒度過程從原始時間序列獲得多尺度時間序列;② 計(jì)算每個粗粒度時間序列的PE值。這兩個步驟可以簡要總結(jié)如下[25]。

(1)

(2) 在MPE分析中，每個粗粒度時間序列的PE值由式(1)～(5)計(jì)算并且比例因子s的函數(shù)可以表示為

(2)

1.1.2 MPE的參數(shù)選擇

在MPE算法中，需要設(shè)置四個參數(shù)，嵌入維度m，時間延遲τ，時間序列長度N和比例因子s。其中嵌入維度m決定輸入的狀態(tài)數(shù)量m!，并且PE值高度依賴于m的選擇。如果m太小該方法將無法工作，因?yàn)椴煌瑺顟B(tài)太少。但是如果m太大，由于條件N≥5m!，可能會對隨后的多尺寸研究造成相當(dāng)大的不便[26]。Bandt等指出尺寸m應(yīng)滿足3≤m≤7。為了檢測信號的動態(tài)變化，通常通過信息損失和計(jì)算復(fù)雜性之間的權(quán)衡關(guān)系來選擇嵌入維數(shù)m。因此本文中m設(shè)為4。時間延遲τ對結(jié)果幾乎沒有影響，所以我們在本文中將τ設(shè)置為1。時間序列N的長度對PE值的估計(jì)有很大影響，N越大計(jì)算時間越久。而太小的N不能滿足條件N≥5m!?？紤]到這些約束，我們在第3章的模擬試驗(yàn)中將數(shù)據(jù)長度N設(shè)置為400，將比例因子s設(shè)置為7。在第4章的研究中，N設(shè)置為1 200，s設(shè)置為12。

1.2 線性局部切線空間排列

線性局部切線空間排列(LLTSA)[27]是一種經(jīng)典的流形學(xué)習(xí)算法，它可以提取高維數(shù)據(jù)集的固有幾何信息，即挖掘隱藏在高維觀察空間中嵌入的低維流形，從而不會丟失數(shù)據(jù)中的重要信息[28-29]。

假設(shè)一個高維數(shù)據(jù)集X={x∈(R)D,i=1,2,…,N}，它是從屬于(R)d的較低維特征空間的基礎(chǔ)流形中采樣的，其中D和d(d

Y=ATXHN

(3)

式中:HN=I-eeT/N表示中心矩陣;I表示單位矩陣;e表示所有元素等于1的列向量。 LLTSA的基本假設(shè)是每個樣本附近的局部結(jié)構(gòu)信息可以用其局部切線空間表示。因此為構(gòu)建從高維輸入空間到低維特征空間的投影矩陣，可以在全局低維特征空間中重新排列所有樣本的局部切線空間。

在LLTSA算法中，根據(jù)樣本之間的歐幾里德距離構(gòu)造采樣點(diǎn)xi(i=1,2,…,N)的鄰域xi=[xi1,xi2,…,xik]，其中k表示最近樣本的數(shù)量。然后用局部變換矩陣Qi來將鄰域xi映射到局部低維切線空間。因此樣本點(diǎn)xi附近的局部結(jié)構(gòu)可以近似表示為

(4)

獲得局部結(jié)構(gòu)后，為得到X的全局低維表示Y，將所有樣本的局部切向空間重新排列在全局低維特征空間中，選擇矩陣為S=[S1,S2,…,SN]其中Si(i=1,2,…,N)是0-1的選擇向量，因此Yi=YSi，其中Yi=[yi1,yi2,…,yik]是Xj的全局低維表示。此步驟的目標(biāo)函數(shù)可以轉(zhuǎn)換為以下最小問題

(5)

(6)

式中,B=SWWTST。用ATXHNBHNXTA=Id確定唯一的Y。最后式(4)的解轉(zhuǎn)化為廣義特征值問題為

ATXHNBHNXTα=λXHNXTα

(7)

可以由對應(yīng)于式(7)的特征向量λ1≤λ2≤…≤λd的d個特征向量α1,α2,…,αd得到轉(zhuǎn)換矩陣A，即A=α1,α2,…,αd。

2 基于MPE-LLTSA的特征提取方法

故障診斷中，對于經(jīng)常出現(xiàn)在機(jī)械領(lǐng)域中的非線性信號，特征選擇是一個難題。EEMD因?yàn)槠渥赃m應(yīng)性而被廣泛用于處理這種非平穩(wěn)信號。流形學(xué)習(xí)在非線性信號中具有更強(qiáng)的降維能力[30]，因此本文引入了流形學(xué)習(xí)算法LLTSA。本文將該特征提取方法縮寫為EEMD-LLTSA, EEMD-LLTSA的技術(shù)路線及其識別方法如圖1(a)所示。

(a) EEMD-LLTSA

然而在復(fù)雜的機(jī)械設(shè)備中，監(jiān)測系統(tǒng)收集的大量非線性信號通常包含多種自然振動模式，這導(dǎo)致單尺度特征提取通常無法表征這些信號。因此引入多尺度方法MPE來改善傳統(tǒng)基于單尺度的非線性動態(tài)參數(shù)的性能。首先通過使用MPE計(jì)算信號x(t)以獲得一組熵值{HMN1,HMN2,…,HMNs},其中s是比例因子。然后構(gòu)建具有高維度的樣本矩陣V，用LLSTA挖掘高維樣本特征矩陣V中的固有結(jié)構(gòu)以適應(yīng)低維特征矩陣。最后用LSSVM建立診斷模型進(jìn)行狀態(tài)分類。本文將該特征提取方法縮寫為MPE-LLTSA，MPE-LLTSA特征提取方法及其條件識別的技術(shù)路線如圖1(b)所示。詳細(xì)步驟描述如下：

步驟1參數(shù)選擇：時間序列x(i),k=1,2,…,N，相應(yīng)參數(shù)有嵌入維數(shù)m，時間延遲τ，時間序列N的長度和比例因子s。

步驟2MPE算法計(jì)算信號的多尺度熵值為

(8)

步驟3對所有樣本進(jìn)行步驟1和2，構(gòu)建高維樣本特征矩陣V為

(9)

式中,p代表樣本數(shù)量。

步驟4通過流形學(xué)習(xí)降低維數(shù)：使用LLTSA降低樣本特征矩陣V的維數(shù)以獲得低維特征向量，如下面的公式所示

(10)

式中:fn表示第n維特征;d是嵌入維度。

步驟5訓(xùn)練和識別：將獲得的低維流形特征分為訓(xùn)練樣本和測試樣本，并使用LSSVM實(shí)現(xiàn)狀態(tài)識別。

3 試驗(yàn)探究

3.1 模擬試驗(yàn)

3.1.1 信號模擬

非線性信號通常具有調(diào)幅[31]的特性，這經(jīng)常出現(xiàn)在機(jī)械領(lǐng)域。根據(jù)振動信號構(gòu)造兩個調(diào)幅信號來驗(yàn)證MPE-LLTSA特征提取方法的性能。

w(t)=am·[1+mcos(2πfmt)]·sin(2πfct)

(11)

式中，am=1，m=8,fm=0.3 Hz，fm=0.3 Hz。對于第一信號w1(t)，fc,1=3，對于第二信號w2(t)，fc,2=3.1 Hz，采樣頻率為20 Hz。

下面用隨機(jī)生成的混合矩陣模擬三個復(fù)雜信號。

(12)

通過上述矩陣M對w1(t)和w2(t)執(zhí)行矩陣運(yùn)算

(13)

添加三個隨機(jī)高斯噪聲n1(t),n2(t) 和n3(t)以分別生成新信號a(t),b(t)和c(t)，如圖2所示。

(a)

(14)

3.1.2 特征評估

如圖2所示，箭頭處的那三個模擬信號的幅度是不同的，但差異不是那么明顯。根據(jù)它們的波形無法準(zhǔn)確區(qū)分這三個信號。為了通過頻率特征區(qū)分三個信號，將FFT方法引入上述三個模擬信號，分別生成其頻域圖，如圖3所示?？梢钥闯?，通過FFT獲得的這三個頻域圖之間的差異并不明顯。通過FFT，我們同樣無法準(zhǔn)確地區(qū)分這些有調(diào)制特性的非平穩(wěn)信號。最后引入Hilbert-Huang變換(HHT)來分析這三類信號的時頻特性，譜圖的結(jié)果如圖4所示。它表明這些簡單的時頻分析方法無法準(zhǔn)確識別出三個信號。

(a)

(a)

上述時頻分析方法無法準(zhǔn)確識別這三種信號。為了解決這個問題，引入了MPE-LLTSA特征提取方法。為每種類型的信號選擇60個樣本，總共180個。每個樣本的數(shù)據(jù)長度N設(shè)置為400，時間長度為20秒。如1.1節(jié)所示，在MPE算法中，分別設(shè)置三個參數(shù)，嵌入維數(shù)m，時間延遲τ和比例因子s。在本次研究中，我們設(shè)定m=

4,τ=1和s=7。如1.2節(jié)所示，在算法LLTSA中，根據(jù)文獻(xiàn)[32]，兩個參數(shù)鄰居數(shù)k和內(nèi)在維度d分別設(shè)置為7和3。a(t),b(t)和c(t) 的每個樣本的MPE結(jié)果顯示在圖5中。通過MPE計(jì)算所有樣本以重建180×12高維特征矩陣。然后通過使用LLTSA方法獲得180×3的低維特征矩陣。降維后的結(jié)果如圖6(a)所示?？梢钥闯鲞@三個模擬信號可以通過MPE-LLTSA算法完全分離，得到的不同信號的三維特征高度集中，聚類效果極佳。為了證明LLTSA方法的優(yōu)越性，還引入了PCA方法。聚類效果如圖6(b)所示。結(jié)果表明，MPE-PCA特征提取方法 還可以準(zhǔn)確識別三個模擬信號，聚類效果也非常好。它表明MPE特性可以有效地克服信號中的噪聲影響，適用于非線性信號處理。

圖5 三個模擬信號的MPE特征

(a) MPE-LLTSA

為了證明MPE-LLTSA優(yōu)于EEMD-LLTSA方法，引入了機(jī)械中常用的能量和SampEn功能進(jìn)行仿真分析。根據(jù)文獻(xiàn)[32]在EEMD算法中集合數(shù)和信噪比(SNR)兩個參數(shù)分別設(shè)置為100和0.3。首先所有樣本都通過EEMD進(jìn)行分解，α(t),b(t)和c(t)的每個樣本的結(jié)果如圖7所示。由于EEMD是主成分分析方法，因此信號的主要信息包含在前幾個IMF中。此外所有樣本獲得的IMF數(shù)量須大于6。因此為了確保EEMD-LLTSA方法的后續(xù)處理，提取前六個IMF的能量和SampEn特征以分別構(gòu)造180×6高維特征矩陣。然后通過LLTSA算法減少特征矩陣以生成180×3的低維特征矩陣(參見圖8)。從圖7和圖8可知，EEMD-energy-LLTSA特征提取方法可以分離三類信號，但聚類效果明顯差于MPE-LLTSA和MPE-PCA。EEMD- SampEn-LLTSA方法甚至無法分離這三個模擬信號，因?yàn)樵谔卣魈崛∵^程中，單尺度特征提取方法僅考慮信號的一個單尺度信息。在某種程度上，它忽略了其他維度中包含的故障信息。這種單尺度特征提取方法經(jīng)常無法在具有嚴(yán)重故障信息耦合的機(jī)械設(shè)備實(shí)現(xiàn)期望的結(jié)果。

圖7 模擬信號α(t),b(t)和c(t)的EEMD結(jié)果

(a) EEMD-energy-LLTSA

3.1.3 聚類評價(jià)

為了準(zhǔn)確地評估聚類能力，引入了三個評價(jià)指標(biāo)，即類內(nèi)散射Sw，類間散射Sb和比率Sw/Sb以描述聚類效應(yīng)。對于特征向量[f1，f2,…,fd]，d是嵌入維度，參數(shù)Sw和Sb的定義為

(15)

(16)

通過計(jì)算類內(nèi)離散度Sw，類間離散度Sb和比率Sw/Sb，可以獲得3.1.2節(jié)中實(shí)現(xiàn)的特征的準(zhǔn)確評估，如表1所示。結(jié)果表明，MPE-LLTSA方法具有良好的分類和聚類能力，優(yōu)于MPE-PCA和EEMD-LLTSA方法。

表1 對三種信號的不同特征進(jìn)行聚類評估

3.2 機(jī)械故障診斷試驗(yàn)結(jié)果

從滾動軸承系統(tǒng)測量的信號非常復(fù)雜。為了進(jìn)一步驗(yàn)證MPE-LLTSA方法的有效性，在這一節(jié)中將使用Case Western Reserve University提供的滾動軸承的試驗(yàn)數(shù)據(jù)[33]。軸承試驗(yàn)系統(tǒng)及其示意圖如圖9所示。它由2馬力的馬達(dá)(左)，扭矩傳感器(中心)，測功機(jī)(右)和控制電子設(shè)備(未示出)組成。試驗(yàn)軸承是6205-2RS JEM SKF深溝球軸承，局部單點(diǎn)損壞分別使用電火花加工機(jī)設(shè)置在外滾道，內(nèi)滾道和滾動體的位置。本文中的振動信號取自安裝在驅(qū)動端的加速度計(jì)，采樣頻率為12 kHz。

圖9 滾動軸承試驗(yàn)系統(tǒng)草圖

3.2.1 四類情況

這種情況下，分析包括健康，內(nèi)在缺陷(IRD)，滾動元件缺陷(RED)和外在缺陷(ORD)在內(nèi)的四類情況。本節(jié)中的振動信號從損壞0.178 cm的軸承中收集。數(shù)據(jù)長度N設(shè)置為1 200。訓(xùn)練和測試樣本的數(shù)量在表2中給出。引入MPE-LLTSA來識別這四個信號。如1.1.3節(jié)所述，在MPE方法中我們設(shè)置m=4，τ=1和s=12。四個信號的波形如圖10(a)及其MPE結(jié)果如圖10(b)所示。對所有樣品通過MPE構(gòu)建高維樣本特征矩陣之后用LLTSA降維得到低維特征。在LLTSA方法中，我們設(shè)置k=7和d=3。計(jì)算結(jié)果如圖11(a)所示，表明MPE-LLTSA可以完全分離這四個信號。為了驗(yàn)證我們提出的MPE-LLTSA方法的優(yōu)越性，還計(jì)計(jì)算了MPE-PCA和EEMD-LLTSA結(jié)果，如圖11(b)，圖11(c)和圖11(d)所示。為得到精確的評估，還計(jì)算了類內(nèi)離散度Sw，類間離散度Sb和Sw/Sb，如表3中所示。最后引入LSSVM方法識別四種狀態(tài)，正確的識別率如表3所示。從圖11和表3可知MPE-LLTSA方法可以區(qū)分這四種狀態(tài)，且其聚類效果優(yōu)于MPE-PCA和EEMD-LLTSA。

表2 不同檢測的試驗(yàn)數(shù)據(jù)的詳細(xì)描述

(a)

(a) MPE-LLTSA

表3 不同檢測下不同特征的聚類評價(jià)

3.2.2 五種狀態(tài)

在表4所示的情況下，分析了滾動體的五種不同的滾動軸承故障(健康的，0.178 cm，0.356 cm，0.533 cm和1.02 cm)，詳細(xì)信息顯示在表4中。數(shù)據(jù)長度N選擇為1 200，訓(xùn)練和測試樣本的數(shù)量分別為50和50。在MPE方法中，我們設(shè)置m=4，τ=1和s=12。在LLTSA方法中，我們設(shè)置k=7和d=3。五個信號的波形圖如圖12(a)及其MPE結(jié)果如圖12(b)所示。MPE-LLTSA的計(jì)算結(jié)果如圖13(a)所示。為了做比較，在圖13(b)，圖13(c)和圖13(d)中列出了MPE-PCA和EEMD-LLTSA的結(jié)果。為了評估更精確，還在表5中計(jì)算了三個參數(shù)Sw，Sb和Sw/Sb。最后引入LSSVM方法識別五種狀態(tài)，并在表5中列出正確的識別率。從圖13和表5，我們可以看到MPE-LLTSA可以完全分離出五類故障，分類結(jié)果和聚類效果優(yōu)于MPE-PCA和EEMD-LLTSA。測試結(jié)果表明，MPE-LLTSA方法可用于軸承缺陷嚴(yán)重性評估。

表4 五個檢測嚴(yán)重程度的試驗(yàn)數(shù)據(jù)的詳細(xì)描述

(a)

(a) MPE-LLTSA

3.3 討 論

我們提出的MPE-LLTSA方法的優(yōu)點(diǎn)主要是由于以下兩個原因：

(1) 能量和SampEn等單尺度特征提取方法只能分析單個尺度的時間序列，導(dǎo)致隱藏在其他尺度中的有用信息丟失。然而非線性或非平穩(wěn)信號通常包含多種自然振蕩模式，并且故障信息通常是復(fù)雜的并且總是隱藏在多尺度時間序列中。MPE方法可以估算不同尺度的時間序列的復(fù)雜度，可以有效地克服信號中的噪聲影響。

(2) 在模式識別領(lǐng)域，由EEMD，MPE等計(jì)算出的特征矩陣的維數(shù)通常很高，難以直接識別這些高維特征。為了提高識別精度，有必要進(jìn)一步提取嵌入在高維特征空間中的低維主要特征。LLTSA具有良好的減小高維非線性數(shù)據(jù)維度的能力。

4 結(jié) 論

本文提出了一種用于機(jī)械故障診斷的MPE-LLTSA特征提取方法。由于機(jī)械設(shè)備結(jié)構(gòu)復(fù)雜，操作環(huán)境惡劣，監(jiān)測系統(tǒng)采集的大量非線性信號通常包含多種自然振動信號，單尺度特征提取方法往往無法表征這些信號。為了解決這個問題，引入MPE來提取非線性信號的多尺度特征并構(gòu)建高維特征矩陣。隨后用LLTSA挖掘嵌入的內(nèi)在結(jié)構(gòu)并進(jìn)行低維特征提取。最后用LSSVM進(jìn)行機(jī)械故障診斷。同時還引入了一些特征提取方法，包括MPE-PCA和EEMD-LLTSA進(jìn)行比較。仿真和實(shí)際案例的測試結(jié)果表明，我們提出的MPE-LLTSA方法具有良好的分類和聚類能力，與其他方法的比較結(jié)果證明了我們提出的MPE-LLTSA方法在機(jī)械故障診斷領(lǐng)域的價(jià)值和應(yīng)用潛力。