李雯玉,鄭小霞
(上海電力大學自動化工程學院,上海 200090)
滾動軸承對各種機械設備而言是不可或缺的重要組件,在工作時不僅處于高速運行狀態(tài),而且受到周期性的部件沖擊,使其極易受到損害并產(chǎn)生故障[1-2]。一旦發(fā)生故障會影響機械設備的整體運行工作,使得安全保障埋下隱患且易造成經(jīng)濟成本損耗,因此滾動軸承狀態(tài)監(jiān)測的相關研究成為近年來的熱點[3-4]。由于旋轉(zhuǎn)機械傳動系統(tǒng)工作環(huán)境的復雜性和多樣性,滾動軸承故障信號往往伴隨著強背景信號與噪聲,并且滾動軸承的振動監(jiān)測信號往往是典型的非平穩(wěn)信號,比較難以準確辨別捕捉。因此,滾動軸承難以提取出有效的特征信息[5-6]。
特征指標的提取是退化狀態(tài)識別的關鍵步驟,所提取特征指標的優(yōu)劣影響著退化狀態(tài)識別結(jié)果的準確性[7]。許多先進的信號處理方法已被應用于滾動軸承的特征提取[8-10]。
李恒等[11]對軸承采集到的振動信號進行短時傅里葉變換(short time fourier transform,STFT)處理后從而提取到特征指標,但是STFT無法兼顧良好頻域分辨率和較高時間分辨率。GUO等[12]通過小波變換(wavelet transform,WT)得到軸承故障信號的時頻圖用于特征提取,但WT難以選擇小波基且信號相位信息不足。WANG等[13]提出了一種基于S變換(S-transform,ST)和卷積神經(jīng)網(wǎng)絡結(jié)合的滾動軸承故障診斷方法,從而解決了STFT存在固定時頻分辨率的難題;同時能維持信號的絕對相位信息;但ST的核心窗函數(shù)固定,使得ST在實際應用受到限制,廣義S變換通過引入窗寬調(diào)節(jié)因子使S變換具有更高的自適應性。
在滾動軸承的性能退化過程中,是無法直接觀察其狀態(tài),而是通過監(jiān)測到的觀測信息與隱含信息間的關系判斷。該過程特性與隱馬爾可夫模型(hidden markov model,HMM) 類似從而可用其來表述。WANG等[14]使用HMM模型對數(shù)據(jù)樣本進行分類用于故障檢測,證明了HMM的有效性。姜萬錄等[15]提出了結(jié)合支持向量聚類和連續(xù)HMM的軸承退化評估方法,說明HMM模型能有效地對軸承運行中的不同退化狀態(tài)做出診斷和預測。但傳統(tǒng)HMM中用于參數(shù)學習問題的Baum-Welch算法容易在訓練過程陷入局部最優(yōu)的情況。基于此,ZHENG等[16]將遺傳算法用于HMM的參數(shù)優(yōu)化,改進后的HMM對光伏逆變器故障的診斷速度更快、更準確,但遺傳算法存在早熟、局部尋優(yōu)能力差的問題。其次,HMM在實際應用中容易對數(shù)據(jù)進行過度擬合,通過采用變分貝葉斯(variational bayesian,VB)改進隱馬爾可夫模型解決該問題。
綜上所述,針對S變換和HMM模型存在的不足,本文提出一種基于廣義S變換和變分貝葉斯-隱馬爾可夫模型(VB-HMM)的滾動軸承性能評估方法。首先,利用廣義S變換對滾動軸承振動信號進行特征指標提取;其次,結(jié)合變分貝葉斯改進后的隱馬爾可夫模型評估滾動軸承性能退化過程;最后,利用軸承退化的實驗數(shù)據(jù)驗證所提出方法的可靠性和準確性。
連續(xù)S變換可基于小波變換或者短時傅里葉變換引出,對于一維信號x(t),其連續(xù)S變換S(τ,f)可定義為:
(1)
式中:S(τ,f)為信號的S變換形式,f為頻率,ω(t-τ,f)為高斯窗函數(shù),τ為對應窗函數(shù)的時間中心,其中:
(2)
S變換具備可逆性,其反變換式定義為:
(3)
由于S變換時頻分辨率受窗函數(shù)形式固定限制,無法靈活調(diào)節(jié)窗寬使其在實際應用中受到限制。因此引入調(diào)節(jié)系數(shù)k來對高斯窗進行改進,其表達式為:
(4)
改進后的廣義S變換表達式為:
(5)
通過改變調(diào)節(jié)系數(shù)k可以改變時頻分辨率,當k=1時,為標準S變換;當k>1時,此時時域分辨率與高斯窗的寬度成反比,窗寬越窄,在分析信號時時域分辨率更為優(yōu)秀;當k∈(0,1)時,在處理信號時的頻域分辨率與窗寬成正比。以靈活改變k的大小的方法來改善S變換使其更為廣泛地運用于實際應用。
為了驗證廣義S變換的效果,進行仿真信號分析:
(6)
圖1顯示了原始信號和當k取不同值時仿真信號的時頻譜。第2幅圖為k=1時的普通S變換;第3幅圖為k=0.4時,廣義S變換的頻域分辨率明顯優(yōu)于第1幅圖;第4幅圖為k=1.5時,此時高斯窗變窄,時域分辨率變高。
圖1 仿真信號的時頻譜
憑借快速傅里葉變換(FFT)使廣義S變換離散化以達到快速處理計算機上實際信號的目的,得到廣義S變換離散化公式,可表示為:
(7)
式中:N為x(t)的采樣點數(shù),T為x(t)的采樣間隔,假設時間序列{X[kT],k=0,1,2,…,N-1}是x(t)的離散時間序列,對被分析的信號X[KT]進行廣義S變換后得到的結(jié)果是一個復時頻矩陣。該矩陣被取模以達到簡化處理的目的,得到廣義S變換的模矩陣,此時矩陣的行列實際意義分別為離散頻率和采樣時間。
信息熵利用狀態(tài)的概率分布來估算系統(tǒng)的復雜度和隨機性[17]。假定若一個離散的隨機變量X=[X1,X2,…,Xm]以概率為P=P(xi)(i=1,2,…,m)的情況發(fā)生,且滿足:
(8)
則X的信息熵可表示為:
(9)
信息熵的熵值越大,表示該事件含有的信息量越少,發(fā)生的概率越低。
廣義S變換是一種具有多分辨率特性的振動信號時頻分析方法,信息熵可以衡量系統(tǒng)的復雜性,評價隨機信號在總體上的不確定性。綜上所述,將廣義S變換和信息熵相結(jié)合,提出結(jié)合兩者優(yōu)勢的GST-熵值指標作為評價滾動軸承的性能特征指標。
采用廣義S變換提取性能特征指標的基本過程如下:
(1)對采樣長度為N的時序信號X進行廣義S變換,得到一個相應矩陣,進一步取模處理,獲得一個N×N的時頻譜系數(shù)矩陣,記為Y。
(2)對矩陣Y的行、列分別按照式(9)進行信息熵處理,將所得N個特征值取平均,分別得到GST-頻率熵與GST-時間熵。
(3)對矩陣Y的時頻譜矩陣整體按照式(9)進行信息熵處理,得到GST-時頻熵。
HMM由兩個隨機過程組成,其中一個隨機過程序列中的狀態(tài)是不可觀測的,而是通過觀測值序列來估計。另一個則可直接通過可觀測量估計其隱藏的狀態(tài)[18]。
HMM包括3個概率矩陣和2個狀態(tài)集合,分別是隱含狀態(tài)S,可觀測狀態(tài)Y,隱含狀態(tài)轉(zhuǎn)移矩陣A,觀測狀態(tài)轉(zhuǎn)移概率矩陣B,模型初始時刻狀態(tài)概率分布π。
在滾動軸承性能評估過程中進行軸承狀態(tài)診斷時,無法通過觀察直接得出其實際狀態(tài),而是通過軸承工作時產(chǎn)生可監(jiān)測的的振動信號作為判斷其狀態(tài)的依據(jù)。該過程可以看作馬爾可夫模型的雙隨機過程。
變分貝葉斯(variational bayesian,VB)的產(chǎn)生是建立在傳統(tǒng)貝葉斯推斷與EM估計算法的基礎上,同時結(jié)合導入了變分近似理論。
模型統(tǒng)計中,用Y表示觀測數(shù)據(jù)樣本,θ、Z分別表示未知模型參數(shù)以及潛變量,從而得到后驗分布p(Z|Y)。在進行求解后驗概率分布時,憑借一種簡單分布q(Z)去代替后驗分布p(Z|Y),以達到消除后驗分布的復雜性對求解精確解的影響[19]。引入度量指標KL散度(kullback-leibler divergence)來描述q(Z)對p(Z|Y)的逼近程度。變分貝葉斯模型的對數(shù)似然函數(shù)可以等價表述為式(10):
(10)
定義L(q)和KL(q‖p)將式(10)的等價表達為式(13):
(11)
(12)
logp(Y)=L(q)+KL(q‖p)
(13)
式中:KL(q‖p)≥0,L(q)為變分自由能。當q(Z)與p(Z|Y)同分布時等號成立,此時KL散度達到最小值,即使變分自由能最大化,使q(Z)≈p(Z|Y),近似分布可認為等價于后驗分布。
傳統(tǒng)的隱馬爾可夫模型與貝葉斯算法存在以下兩個問題:
(1)隱馬爾可夫模型擬合時的參數(shù)量可能大于可行的數(shù)據(jù)量,模型內(nèi)存在著許多參數(shù)和觀測數(shù)據(jù)。在進行參數(shù)求解時,多數(shù)時候無法避免模型復雜度導致容易對數(shù)據(jù)過度擬合的情況,這是由于傳統(tǒng)的HMM通常選取極大似然估計的方法。
(2)貝葉斯算法無法在隱馬爾可夫模型中直接運用,通常而言,信號分析的有效性受參數(shù)邊緣化過程影響,由于該過程往往伴隨著積分計算產(chǎn)生相應的運算量。
為解決上述問題,可以采用變分貝葉斯代替極大似然估計來估計參量從而改進隱馬爾可夫模型。具體步驟為:
首先,式(14)和式(15)描述的似然函數(shù)分別對應隱馬爾可夫模型中的隱含狀態(tài)和觀測數(shù)據(jù)。
(14)
(15)
參數(shù)可在基于貝葉斯推斷的原理上將后驗密度表達為式(16):
(16)
此時變分貝葉斯的應用是為了進一步處理式(16)中無法計算的分母,在該步驟中為了獲取期望的估計q(θ,S)應當使下界L(q)取最大值,同時避免過擬合問題。邊緣似然為:
(17)
通過對式(17)兩邊取對數(shù)后再取后驗分布q(θ,S)的期望,得到:
logP(Y)=?q(S,θ)logP(Y,S,θ)dSdθ=
-?q(S,θ)logP(S,θ|Y)dSdθ
(18)
整理式(18),可得:
(19)
可將式(19)表達為式(11)~式(13)的形式。此時KL散度達到最小值,即下界L(q)最大,來實現(xiàn)logp(Y)的估計。
下界可表示為:
(20)
式中:S為狀態(tài),θ為參量。
變分貝葉斯改進的迭代過程中需要一個VBE步驟來推斷隱含變量,以及一個VBM步驟用來推斷模型參數(shù)的后驗分布。通過對式(20)求偏導函數(shù)獲取兩個步驟中運算所需要迭代的方程。
(1)VBM-步驟。
(21)
(22)
(23)
(2)VBE-步驟。
(24)
(25)
(26)
VB-HMM模型采用前向-后向算法,表達為:
(27)
(28)
若下界L(q)處于可以忽略不計其極其微小的變化時,此時滿足了VBEM算法的終止條件。
首先選取正常狀態(tài)下滾動軸承的數(shù)據(jù)作為訓練樣本,應用2.3節(jié)的變分貝葉斯-隱馬爾可夫模型建立性能評估模型,得到軸承正常階段的VB-HMM模型;然后將待測數(shù)據(jù)輸入已訓練好的VB-HMM模型,計算輸出其偏離值,實現(xiàn)對滾動軸承的性能退化分析。
具體的滾動軸承性能評估流程如圖2所示。
圖2 基于VB-HMM的性能評估流程
步驟1:根據(jù)圖2的特征指標提取流程,提取出正常狀態(tài)下的滾動軸承特征向量,訓練得到正常狀態(tài)的VB-HMM模型;
步驟2:對待測數(shù)據(jù)同樣經(jīng)歷步驟1的特征提取,得到樣本將其輸入到訓練好的模型中;
步驟3:通過模型計算輸出測試樣本的對數(shù)似然概率值,達到性能評估的目的。
根據(jù)故障位置的不同,滾動軸承的點蝕故障可以分為滾動體故障、內(nèi)圈故障和外圈故障。因此,假定滾動軸承承受不同沖擊強度A時故障數(shù)學模型可以表示為:
(29)
式中:Ai為調(diào)制幅值,τi為微小滑動,n(t)為高斯噪聲,fn為系統(tǒng)共振頻率,B為系統(tǒng)共振衰減系數(shù)。
根據(jù)美國辛辛那提大學(IMS)提供的軸承全壽命疲憊實驗數(shù)據(jù)進行實驗。該實驗裝置中含有4個施加了2700 kg徑向載荷的軸承,轉(zhuǎn)速保持在2000 r/min,軸承節(jié)圓直徑為7.15 cm,一個周期的數(shù)據(jù)長度為2048,滾動體的直徑為0.84 cm[20]。總共含有3組軸承的全壽命實驗數(shù)據(jù),本文使用數(shù)據(jù)集2的數(shù)據(jù),采樣頻率為20 480 Hz,軸承振動信號采樣間隔為10 min/次,共采集到984組樣本。沖擊強度取值為A0={1,5,10},根據(jù)故障類型參考表1分別取值沖擊幅值調(diào)制頻率fm和沖擊周期T。
表1 fm和T不同故障取值
當滾動軸承處于內(nèi)圈故障、外圈故障和滾動體故障時,對一個周期內(nèi)的仿真信號分別提取GST-時間熵、GST-頻率熵、GST-時頻熵特征指標。結(jié)果如圖3~圖5所示。處于3種不同的故障類型時,隨著沖擊強度加大,滾動軸承的故障程度逐漸加深,GST-時間熵、GST-頻率熵、GST-時頻熵-均逐漸減小。該變化驗證了GST熵值特征指標的有效性。
圖3 內(nèi)圈故障時不同GST熵值指標
圖4 外圈故障時不同GST熵值指標
圖5 滾動體故障時不同GST熵值指標
為了進一步驗證GST-時間熵、GST-頻率熵、GST-時頻熵3個特征指標在滾動軸承實際性能退化過程中的可行性,采用IMS提供的軸承全壽命疲憊實驗數(shù)據(jù)進行實驗驗證。波形指標是應用較為廣泛成熟的時域指標,表達式為:
(30)
同時將GST熵值指標與波形指標在滾動軸承全壽命周期中的變化情況進行對比。
圖6為辛辛那提軸承振動信號的時域波形和經(jīng)過廣義S變換后的時頻譜圖。
(a) 時域波形 (b) 二維時頻譜圖
(c) 三維時頻譜圖圖6 辛辛那提軸承振動信號GST時頻譜
圖6b和圖6c分別是軸承信號經(jīng)過廣義S變換后的二維時頻譜和三維時頻譜。圖6c中X軸為時間軸,Y軸為頻率軸,Z軸表示幅度。圖6c顯示了軸承振動信號經(jīng)廣義S變換后能量隨頻率和時間的變化,具有較高的時頻分辨率。
圖7分別是波形指標、GST時頻熵、GST時間熵和GST頻率熵在全壽命周期里隨時間變化的曲線。曲線變化過程共分為4個階段。第1階段曲線基本處于平穩(wěn)狀態(tài),該階段滾珠軸承處于正常工作的健康狀態(tài)。第2階段曲線開始上升或者下降,該階段滾動軸承處于早期微弱故障階段。第3階段曲線產(chǎn)生震蕩,但震蕩幅度較小,該階段滾動軸承故障程度隨著時間加劇。第4階段曲線開始劇烈震蕩,震蕩頻率較上一階段明顯增大,直到最后上升或者下降至新的極點,該階段滾動軸承經(jīng)過不斷磨損逐漸失效。整個曲線展示了滾動軸承從健康狀態(tài)到故障逐步加深直至失效的過程,曲線變化與滾動軸承實際變化較為一致,驗證了GST熵值指標的有效性。
圖7 辛辛那提全壽命周期不同特征指標曲線
表2為不同特征指標在滾動軸承不同階段的對比表,結(jié)合圖7分析可得以下結(jié)論:
表2 不同特征指標對比表
(1)波形指標相對于GST熵值指標,對于滾動軸承從早期故障到嚴重故障直至失效階段的變化并不明顯。
(2)GST-時間熵相對其他3種特征指標對滾動軸承早期微弱故障最為敏感,反應時間早,有助于初期的維護運行。
(3)GST-頻率熵對早期故障靈敏度相比于GST-時頻、GST-時間熵較低,但嚴重故障到失效階段反應時間有所提前,有助于在滾動軸承嚴重故障時及時檢測停機維修。
為了更好的描述滾動軸承的性能退化失效過程,便于劃分失效狀態(tài),方便后續(xù)數(shù)據(jù)處理更加便捷快速,通過式(31)對GST-熵值指標進行歸一化處理,獲取相應的健康指數(shù)(health index,HI)指標,將取值范圍映射到[0,1]之間。
(31)
經(jīng)過處理后的健康指數(shù)曲線波動范圍在0~1之間,初期健康指數(shù)曲線較為平穩(wěn)且數(shù)值較小,隨著滾動軸承不斷磨損故障加深,曲線開始上升并在達到一定之后開始波動。
為了進一步應用于實際應用中,將滾動軸承的健康狀態(tài)劃分為4個狀態(tài)。在健康指數(shù)曲線處于較為穩(wěn)定階段,實時健康等級評估為健康;曲線開始上升時,滾動軸承出現(xiàn)損耗進入早期微弱故障,實施健康等級評估為亞健康;當曲線出現(xiàn)波動后跌至新的最低點,此時滾動軸承性能嚴重退化,根據(jù)曲線震蕩幅度的大小可分為監(jiān)控運行與建議停機2個階段。
為了進一步描述滾動軸承性能退化過程,采用上文所說的變分貝葉斯-隱馬爾可夫狀態(tài)評估模型,選取上述數(shù)據(jù)集2中采集到的984組數(shù)據(jù)樣本進行實驗,根據(jù)圖2的VB-HMM性能評估流程,分別提取GST-時間熵健康指數(shù)、GST-頻率熵健康指數(shù)、GST-時頻熵健康指數(shù)和波形指標健康指數(shù)4個特征指標作為VB-HMM模型的輸入。將正常狀態(tài)下100組軸承數(shù)據(jù)用來訓練模型,整體數(shù)據(jù)為測試數(shù)據(jù)。
為了進一步驗證VB-HMM評估模型在性能退化中準確性,分別使用采用傳統(tǒng)HMM和VB-HMM兩種模型進行性能評估。
全壽命周期的HMM變換曲線如圖8a所示,在滾動軸承運行的1~696組期間,樣本的曲線處于平穩(wěn)狀態(tài),說明其運行屬于正常狀態(tài);在697~755組階段,曲線呈緩慢下降趨勢,但幅度較小,對數(shù)似然概率值逐漸減小,此時滾動軸承開始出現(xiàn)輕微磨損,性能開始劣化并保持該狀態(tài)繼續(xù)運行;在隨后的756~984組中,曲線開始處于較大幅度的波動,并在最后跌至最低點,說明該期間的滾動軸承磨損嚴重,出現(xiàn)較大故障,性能退化至失效。
(a) HMM (b) VB-HMM圖8 全壽命周期似然對數(shù)曲線
圖8b中,在1~475組階段中滾動軸承處于正常運行階段,曲線保持平穩(wěn)。在第476組樣本處出現(xiàn)了早期故障,曲線逐漸下降偏離正常狀態(tài)。退化曲線在684~852組過程中產(chǎn)生波動震蕩,此時滾動軸承開始故障惡化;從853組直至曲線結(jié)束,滾動軸承深度劣化直至損壞,過大的沖擊使曲線震蕩較上一階段更為頻繁,到最后完全失效對數(shù)似然概率值達到最小值。
基于兩種方法的對比結(jié)果如表3所示,結(jié)合圖8分析,HMM模型的對數(shù)似然概率值相差35個單位,VB-HMM的取值區(qū)間為120個單位,所以VB-HMM的變化趨勢更加清晰明顯。HMM模型在第696組樣本處開始出現(xiàn)早期故障,但曲線下降程度微弱,VB-HMM模型在第476組樣本處開始出現(xiàn)早期故障,說明VB-HMM模型在滾動軸承早期微弱故障階段更為敏感;同時,VB-HMM模型在第684本樣本后可分為兩個階段,684~852組為磨損程度逐漸加深的階段,853~984組為失效階段,而HMM模型在756~984組過程中曲線震蕩幅相差不大。由此可以看出:在滾動軸承從健康到劣化失效整個期間,VB-HMM模型比HMM模型更為清晰有效的表征退化過程。
表3 HMM和VB-HMM結(jié)果對比
且結(jié)合3.3節(jié)的健康指數(shù)分類實時評估,可得滾動軸承在1~475組階段為正常工作狀態(tài),處于健康狀態(tài)的時間為4760 min;476~683組滾動軸承為早期故障,處于亞健康狀態(tài)時間為2070 min;684組~852組滾動軸承故障加劇,處于監(jiān)控運行階段時間為1690 min;從853組直至結(jié)束,滾動軸承嚴重劣化已影響正常工作,處于建議停機的時間為1320 min。滾動軸承從監(jiān)控運行轉(zhuǎn)為建議停機狀態(tài)的轉(zhuǎn)折時刻可作為最佳預警的時間點,即平均提前50.1 h啟動預警。
本文VB-HMM模型對辛辛那提大學滾動軸承數(shù)據(jù)性能評估結(jié)果,與文獻[21]中:判斷1~539組為正常運行狀態(tài),處于540~699組時,滾動軸承為緩慢損壞階段,在700~984組階段,滾動軸承處于嚴重損壞直到失效階段。兩者在判斷正常運行狀態(tài)以及早期故障和故障加劇階段結(jié)論基本一致,但本文在判斷早期故障階段更為靈敏;并進一步區(qū)分了故障加劇和嚴重故障至失效階段,在實際檢測維修應用方面更為靈活。這驗證了本文變分貝葉斯-隱馬爾可夫模型的可行性。
計算ZA-2115軸承各部件的故障特征頻率,如表4所示。
表4 ZA-2115軸承故障特征頻率
為了驗證4.1節(jié)評估結(jié)果的正確性,對健康等級變換節(jié)點前后進行包絡譜分析,分別為第420組樣本、第476組樣本、第684組樣本和第853樣本,分析圖如圖9所示。
(a) 第420組 (b) 第476組
(c) 第684組 (d) 第853組圖9 健康等級變換節(jié)點數(shù)據(jù)包絡譜分析圖
由圖9可知,圖9a中,第420組數(shù)據(jù)樣本無故障特征頻率,滾動軸承為正常工作狀態(tài);圖9b中,第476組數(shù)據(jù)樣本包絡譜圖出現(xiàn)了236 Hz,與表4中外圈故障的頻率接近,可判斷該滾動軸承的故障類型是外圈故障,且滾動軸承處于早期故障階段;圖9c中,第684組數(shù)據(jù)樣本頻率幅值增大,且外圈故障1倍頻、2倍頻和3倍頻特征明顯,說明滾動軸承進一步劣化,為監(jiān)控運行階段;圖9d中,第853組數(shù)據(jù)樣本出現(xiàn)外圈故障1倍頻、2倍頻、3倍頻和4倍頻,此時故障嚴重直至失效,建議停機。
各個變化節(jié)點的包絡譜變化與4.1節(jié)VB-HMM模型實時評估結(jié)果是一致的,證明了該模型的有效性。
結(jié)合廣義S變換與信息熵對滾動軸承進行性能退化特征指標提取,并計算相應健康指數(shù),通過滾動軸承模擬信號仿真實驗證明了廣義S變換熵值指標的可行性;并通過辛辛那提大學的滾動軸承全壽命實驗數(shù)據(jù)進一步說明了廣義S變換熵值指標在整個性能退化過程的變化趨勢,該變化趨勢與滾動軸承從健康狀態(tài)到損壞狀態(tài)的動態(tài)變化趨勢相符合,驗證了其有效性。
最后,結(jié)合變分貝葉斯改進后的隱馬爾可夫模型,建立了滾動軸承的狀態(tài)評估模型,采用標準實驗數(shù)據(jù)對模型的有效性進行了驗證,其曲線變化良好的描述了軸承從正常運行到早期故障直至損壞的退化過程,對滾動軸承實時監(jiān)控和運行評估具有一定的實際應用價值。