鄭小霞,劉 靜,魏彥彬,蔣海生
(上海電力大學(xué) 自動(dòng)化工程學(xué)院,上海 200082)
滾動(dòng)軸承在旋轉(zhuǎn)機(jī)械上應(yīng)用比較廣泛,它的工作狀態(tài)影響著運(yùn)行機(jī)械的工作精度[1]。滾動(dòng)軸承是發(fā)電機(jī)傳動(dòng)系統(tǒng)中最易受損的主要部件之一,若發(fā)生故障會對整個(gè)機(jī)組的運(yùn)行狀態(tài)產(chǎn)生很大影響,造成很大的經(jīng)濟(jì)損失[2]。因此,為保證發(fā)電機(jī)組的可靠性和安全性,降低維修成本,并提高其運(yùn)行效率,對滾動(dòng)軸承進(jìn)行健康狀態(tài)監(jiān)測具有重要的意義[3]。
滾動(dòng)軸承的振動(dòng)監(jiān)測信號是非平穩(wěn)的復(fù)雜信號,其頻率成分復(fù)雜,且根據(jù)實(shí)際狀況的影響,該信號會不斷發(fā)生變化,諧波、干擾與故障特征頻率之間也容易出現(xiàn)重疊現(xiàn)象。因此,采用傳統(tǒng)的方法難以提取其故障信息[4]。
小波包分解是一種典型的復(fù)雜非平穩(wěn)信號處理方法,能對信號進(jìn)行精確的處理,被廣泛用于各類設(shè)備的狀態(tài)監(jiān)測中[5,6]。文獻(xiàn)[7]通過小波包分解提取信號特征,進(jìn)行了發(fā)電機(jī)組齒輪箱的故障診斷。文獻(xiàn)[8]為了解決運(yùn)行中監(jiān)測信號噪聲的影響,基于小波包分解,提出了一種評價(jià)不同頻帶故障的方法。
互補(bǔ)聚合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)是由YEH J R[9]提出的一種常用的信號特征提取方法,能對信號進(jìn)行很好的分解,并減少重構(gòu)誤差的產(chǎn)生。文獻(xiàn)[10]通過CEEMD方法,將風(fēng)速劃分為一組固有模態(tài)函數(shù),并結(jié)合極限學(xué)習(xí)機(jī)(ELM)實(shí)現(xiàn)了對風(fēng)速的預(yù)測,但是該方法不能完全消除模態(tài)混疊。
對于提取的故障特征信號,需要選擇合適的模型來進(jìn)行設(shè)備的健康狀態(tài)評價(jià)。隱馬爾可夫模型(HMM)是由LeonardE.Baum提出的一種描述隨機(jī)過程的概率模型,該模型對于復(fù)雜的時(shí)間序列有很好的建模能力,被廣泛應(yīng)用于設(shè)備的壽命預(yù)測和狀態(tài)評價(jià)[11]。文獻(xiàn)[12]通過隱馬爾可夫模型進(jìn)行了齒輪箱的功能退化檢測,推導(dǎo)出了其性能的下降規(guī)律。但在實(shí)際應(yīng)用中,隱馬爾可夫模型容易出現(xiàn)過度擬合問題;而變分貝葉斯(VB)優(yōu)化的隱馬爾可夫(HMM)模型對發(fā)電機(jī)組軸承部件則有著更好的適用性。
綜上所述,筆者提出一種基于小波包-互補(bǔ)聚合經(jīng)驗(yàn)?zāi)B(tài)分解和變分貝葉斯-隱馬爾可夫模型的滾動(dòng)軸承健康狀態(tài)評價(jià)方法,并利用軸承退化的實(shí)驗(yàn)數(shù)據(jù)對該模型進(jìn)行驗(yàn)證。
互補(bǔ)聚合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)適合非線性信號分析,可自適應(yīng)地將復(fù)雜信號分解成由高頻到低頻的特征模態(tài)函數(shù)(IMF)[13]。CEEMD算法原理與EEMD算法一致,只是分解過程中加入了成對的白噪聲,能夠減少由白噪聲引起的誤差,更準(zhǔn)確地對原始信號進(jìn)行重構(gòu)[14]。
由于原始信號中的頻率成分比較復(fù)雜,CEEMD很難得到單分量IMF。針對這一問題,可選擇小波包分解(WPD)對原始振動(dòng)信號進(jìn)行預(yù)處理,以降低振動(dòng)信號的復(fù)雜性,抑制CEEMD分解過程中出現(xiàn)模態(tài)混疊問題。WPD不僅能對低頻部分進(jìn)行分解,對高頻信號也有具很好的分解能力,保證其良好的分辨率。
應(yīng)用小波包分解時(shí),首先要根據(jù)分解的數(shù)據(jù)類型進(jìn)行小波基函數(shù)的選取。筆者主要依據(jù)降噪后信號的完整度和準(zhǔn)確度選擇小波基函數(shù)。另外,小波包的分解層數(shù)的選取也很重要,分解層數(shù)過小,不能體現(xiàn)信號的細(xì)節(jié)信息;分解層數(shù)過大,計(jì)算量會增加。所以根據(jù)研究和計(jì)算,筆者采用db3小波基對原始信號進(jìn)行3層分解,并計(jì)算第3層各節(jié)點(diǎn)信號的能量占比,取總能量高于80%,且能量較大的信號進(jìn)行重構(gòu),由此可以消除大部分的干擾信號,為下一步的CEEMD分解奠定基礎(chǔ)。
采用改進(jìn)的CEEMD方法,首先要對原始信號進(jìn)行小波包降噪,以降低信號的頻率分量,再利用CEEMD將降噪后的信號分解為IMF分量,根據(jù)相關(guān)系數(shù)選擇有效分量進(jìn)行重構(gòu),最后提取出故障的特征信號。
HMM模型是一種用來描述含有隱藏未知參數(shù)的馬爾可夫過程的統(tǒng)計(jì)模型,其狀態(tài)序列為隱藏的馬爾可夫鏈隨機(jī)生成的不可觀測序列;每個(gè)狀態(tài)生成一個(gè)觀測,再由此產(chǎn)生的觀測的隨機(jī)序列,稱為觀測序列[15]。
假設(shè)觀測序列Y={y1,y2,…yr},可以表示為二變量過程(S,Y)={(st,yt),t=1,2,…,T}。隱含過程{st}是一個(gè)馬爾可夫鏈。
隱馬爾可夫模型可用θ=(π,A,B)來描述,分別定義為:
(1)狀態(tài)轉(zhuǎn)移概率矩陣:A={aij,i,j∈S},aij=P(st=j|st-1=i);
(2)發(fā)射概率矩陣:B={bim,i,m∈S},bjm=P(yt=m|st=i);
(3)初始狀態(tài)分布:π={πi,i∈S},πi=P(s1=i)
假設(shè)A、B和π之間相互獨(dú)立,先驗(yàn)概率為:
P(θ)=P(π)P(A)P(B)
(1)
由于狄利克雷分布與似然項(xiàng)共軛,所以選擇參數(shù)π、A、B的先驗(yàn)概率為Dirichlet分布:
(2)
(3)
(4)
式中:uA,uB,uπ—狄利克雷分布參數(shù);M—觀測數(shù)目;I—狀態(tài)數(shù)。
在貝葉斯算法中,對參數(shù)邊緣化過程中需要進(jìn)行積分運(yùn)算,這個(gè)計(jì)算會影響信號處理的可行性,所以此處采用變分貝葉斯來解決這個(gè)問題。
logp(Y)=logp(Y,Z)-logp(Z|Y)
(5)
logp(Y)=log{p(Y,Z)/q(Z)}-
log{(p(Z│Y)/q(Z)}
(6)
兩邊積分為:
(8)
后驗(yàn)分布p(Z|Y)的估計(jì)和模型證據(jù)p(Y)可表示為:
logp(Y)=L(q)+KL(q‖p)
(9)
其中,定義:
(10)
(11)
式中:L(q)—下界;KL(q‖p)—q和真實(shí)后驗(yàn)p之間的KL。
使KL最小,相當(dāng)于使L(q)最大,最大化過程即是通過迭代潛在變量和參量的后驗(yàn)概率。
隱馬爾可夫模型求解參數(shù)時(shí)采用的是極大似然估計(jì)(maximum likelihood estimate,MLE),沒有考慮模型的復(fù)雜性。而對于一個(gè)隱馬爾可夫模型而言,其復(fù)雜性與模型中的狀態(tài)轉(zhuǎn)移矩陣A的連通度、隱藏狀態(tài)數(shù)量t和發(fā)射矩陣B中,每個(gè)隱藏狀態(tài)對應(yīng)觀察狀態(tài)的概率分布有關(guān)。所以,在實(shí)際應(yīng)用中,如果有很多觀察數(shù)據(jù),擬合參數(shù)的數(shù)量可能會超過可用的數(shù)據(jù)量,造成結(jié)果的不準(zhǔn)確。
為解決這個(gè)問題,筆者提出了一種變分貝葉斯算法改進(jìn)的隱馬爾可夫模型,采用變分貝葉斯(VB)代替極大似然估計(jì)(MLE),來求解HMM中的模型參數(shù),對發(fā)電機(jī)組復(fù)雜信號的建模具有更好的效果。
變分貝葉斯(VB)是一種在貝葉斯估計(jì)和機(jī)器學(xué)習(xí)領(lǐng)域中用于近似復(fù)雜積分的技術(shù)[16]。在貝葉斯算法中,在參數(shù)邊緣化時(shí)往往需要積分運(yùn)算,而在進(jìn)行實(shí)際信號處理時(shí),該運(yùn)算會影響信號處理的可行性。采用變分貝葉斯則可解決了一問題。
變分貝葉斯(VB)本質(zhì)上是一種近似算法,即對難以求解的后驗(yàn)分布進(jìn)行近似求解。在變分推斷中,經(jīng)常用一個(gè)簡單分布q(Z:φ)來近似復(fù)雜分布p(Z|Y,θold),從而得到一個(gè)局部最優(yōu),且有確定解的近似后驗(yàn)分布。
變分貝葉斯改進(jìn)隱馬爾可夫的具體步驟為:
在上文的HMM中,隱含狀態(tài)和觀測數(shù)據(jù)的似然函數(shù)為:
(12)
(13)
由貝葉斯公式得出參數(shù)的后驗(yàn)密度為:
(14)
其中,分母中的積分是無法處理的,因此,此處引進(jìn)變分貝葉斯作為一種處理方法。
邊緣似然表示為:
(15)
兩邊對其取對數(shù),再取關(guān)于后驗(yàn)分布q(S,θ)的期望,有:
(16)
式中:q(S,θ)—潛在變量和參數(shù)的變分后驗(yàn)。
把它們的組合表示為Z(S,θ)。logP(Y)的估計(jì)可以通過使KL散度最小來實(shí)現(xiàn),也等價(jià)為使下界最大。
下界可表示為:
(17)
(18)
對上式關(guān)于參量θ狀態(tài)S的后驗(yàn)求偏導(dǎo),可得到VBM-step和VBE-step中的迭代方程:
(1)VBM-step:
(19)
(20)
(21)
(2)VBE-step:
(22)
(23)
(24)
此處仍應(yīng)用前向-后向算法,有:
(25)
(26)
設(shè)前向概率值為DI,其可表示為:
(27)
發(fā)電機(jī)組運(yùn)行環(huán)境特殊,且難以取得完備的運(yùn)行數(shù)據(jù),故筆者將采集的健康狀態(tài)信號作為模型訓(xùn)練樣本,建立變分貝葉斯-隱馬爾可夫(VB-HMM)健康狀態(tài)模型;再將監(jiān)測到的當(dāng)前數(shù)據(jù)輸入到模型中,得到當(dāng)前狀態(tài)與健康狀態(tài)的接近度,實(shí)現(xiàn)對軸承的健康狀態(tài)評價(jià)。
在所建立的VB-HMM模型中,每一組數(shù)據(jù)對應(yīng)一個(gè)對數(shù)似然概率值,來表示健康狀態(tài)評價(jià)結(jié)果。
建立發(fā)電機(jī)組滾動(dòng)軸承健康狀態(tài)評價(jià)模型的具體步驟如下:
(1)采用小波包分解對軸承振動(dòng)信號進(jìn)行預(yù)處理,提取能量較大且總能量高于80%的節(jié)點(diǎn)信號;
(2)將步驟(1)提取的節(jié)點(diǎn)信號進(jìn)行CEEMD分解,取相關(guān)系數(shù)大于0.9的IMF分量作為特征信號;
(3)用健康數(shù)據(jù)訓(xùn)練健康狀態(tài)模型;
(4)將樣本數(shù)據(jù)輸入到健康狀態(tài)模型中,得到對數(shù)似然概率值,對訓(xùn)練結(jié)果進(jìn)行分析。
此處筆者利用美國辛辛那提大學(xué)智能維護(hù)系統(tǒng)中心(IMS Center)提供的軸承全壽命疲勞試驗(yàn)數(shù)據(jù)進(jìn)行軸承的退化過程驗(yàn)證。
該實(shí)驗(yàn)裝置包括4個(gè)雙列滾子軸承,筆者在每個(gè)軸承外座處安裝一個(gè)采樣頻率為20 kHz的PCB353B33加速度傳感器,并對軸承施加266 689 N的徑向力。
實(shí)驗(yàn)中軸承傳感器的布置示意圖如圖1所示。
圖1 傳感器布置示意圖
筆者將軸承的部分故障原始信號進(jìn)行小波包分解;采用db3小波基,將信號分解為3層,此時(shí)可保證故障信號完整,且可最大限度地消除干擾信號。
進(jìn)行小波包分解后,第3層的8個(gè)頻帶節(jié)點(diǎn)的能量分布如圖2所示。
圖2 軸承故障振動(dòng)信號頻帶能量分布
由圖2可看出:軸承故障振動(dòng)信號的能量主要集中在節(jié)點(diǎn)3和節(jié)點(diǎn)7,占到總能量的80%以上,因此,此處對節(jié)點(diǎn)3和節(jié)點(diǎn)7進(jìn)行重構(gòu)。
原始信號和重構(gòu)后的信號頻域波形如圖3所示。
圖3 原始信號和小波包分解后的頻域波形
對軸承的原始振動(dòng)信號進(jìn)行小波包降噪預(yù)處理后,筆者根據(jù)能量分布提取能量較高的節(jié)點(diǎn)3和節(jié)點(diǎn)7,再進(jìn)行CEEMD分解,以消除其他干擾信號和虛假分量,提取出故障的特征信號。
節(jié)點(diǎn)3和節(jié)點(diǎn)7分解后,IMF分量的時(shí)域波形如圖4所示。
圖4 IMF分量時(shí)域波形
兩節(jié)點(diǎn)重構(gòu)后的IMF分量與小波分解后的原信號的相關(guān)系數(shù)如表1所示。
表1 IMF分量與小波分解后的原信號的相關(guān)系數(shù)
提取相關(guān)系數(shù)高于0.9的分量作為故障特征信號,其余IMF信號與原信號相似性較低,作為干擾信號舍去。
筆者選擇辛辛那提大學(xué)軸承疲勞數(shù)據(jù)中的SET-2中的849組數(shù)據(jù)作為實(shí)驗(yàn)數(shù)據(jù),建立軸承的全壽命周期曲線。每組數(shù)據(jù)有20 480個(gè)數(shù)據(jù)。為方便后面的計(jì)算,把每組數(shù)據(jù)分為20段,每段1 024個(gè)數(shù)據(jù);選擇其中的200組健康狀態(tài)數(shù)據(jù)訓(xùn)練健康狀態(tài)的模型參數(shù),剩余的649組作為測試數(shù)據(jù);將所有數(shù)據(jù)進(jìn)行小波包分解,并根據(jù)上文中的分析提取3、7節(jié)點(diǎn)的數(shù)據(jù)進(jìn)行CEEMD分解,提取相關(guān)系數(shù)高于0.9的IMF分量作為特征向量,進(jìn)行健康模型訓(xùn)練。
為比較隱馬爾可夫模型(HMM)和變分貝葉斯-隱馬爾可夫模型(VB-HMM)的區(qū)別,此處采用相同的數(shù)據(jù)分別對兩種模型進(jìn)行訓(xùn)練和測試,得到的結(jié)果如圖5所示。
圖5 似然對數(shù)曲線
圖5中,為方便進(jìn)行比較,兩個(gè)曲線圖取相同的縱坐標(biāo)。從縱坐標(biāo)來看,圖5(a)曲線的極值相差70個(gè)單位的DI值,圖5(b)中極值相差150個(gè)單位的DI值,所以圖5(b)中的曲線變化趨勢更加明顯。
由圖5可明顯看出:200組后,圖5(b)曲線開始明顯下降,但圖5(a)曲線沒有表現(xiàn)出下降趨勢;330組后,圖5(b)曲線開始有大幅度振蕩,圖5(a)中只有較小幅度振蕩。由此可見,變分貝葉斯-隱馬爾可夫(VB-HMM)模型更適合軸承健康模型的建立。
再用圖5(b)對實(shí)驗(yàn)臺軸承的健康狀態(tài)進(jìn)行分析。圖5(b)中,進(jìn)行軸承健康狀態(tài)檢測的有649組數(shù)據(jù)。在實(shí)驗(yàn)開始的1組~200組中,曲線很平穩(wěn),且對數(shù)似然概率值(DI)較大,說明這200組數(shù)據(jù)是健康數(shù)據(jù),軸承處于健康狀態(tài);隨著時(shí)間推移,201組~330組的曲線開始緩慢下降,表明軸承性能逐漸偏離健康狀態(tài),開始劣化;在實(shí)驗(yàn)后期,331組~649組數(shù)據(jù)的曲線開始有較大波動(dòng)且下降嚴(yán)重,并在最后下降到最低點(diǎn),表明軸承性能處于嚴(yán)重受損狀態(tài)。
由此可以看出:整個(gè)曲線可以描述軸承從健康狀態(tài)到逐步劣化,再到失效的整個(gè)過程。所以,本文的健康狀態(tài)評價(jià)模型能夠準(zhǔn)確地評價(jià)軸承的健康狀態(tài)。
為了進(jìn)一步說明采用變分貝葉斯-隱馬爾可夫模型對軸承數(shù)據(jù)進(jìn)行準(zhǔn)確評價(jià)的有效性,筆者采用滾動(dòng)軸承加速實(shí)驗(yàn)平臺的數(shù)據(jù)進(jìn)行驗(yàn)證。采用測量裝置采集滾動(dòng)軸承整個(gè)壽命周期的振動(dòng)信號。其中,加速度計(jì)型號為DYTRAN3035B,振動(dòng)信號的采樣頻率為25.6 kHz;每次采樣間隔為10 s,采樣個(gè)數(shù)為2 560個(gè);滾動(dòng)軸承的轉(zhuǎn)速為1 800 r/min,載荷為4 000 N。
筆者將軸承振動(dòng)信號的部分故障數(shù)據(jù)進(jìn)行特征提取。首先對原始振動(dòng)數(shù)據(jù)進(jìn)行小波包分解預(yù)處理,分解后的能量分布如圖6所示。
圖6 軸承故障振動(dòng)信號頻帶能量分布
圖6中,該軸承信號的能量集中在節(jié)點(diǎn)6,所以筆者提取節(jié)點(diǎn)6進(jìn)行信號重構(gòu),并進(jìn)行CEEMD分解。
各IMF分量與小波分解后的信號的相關(guān)系數(shù)如表2所示。
表2 IMF分量與小波分解后的原信號的相關(guān)系數(shù)
筆者取相關(guān)系數(shù)高于0.9的IMF分量,作為特征信號進(jìn)行模型訓(xùn)練;選取1 860組平臺數(shù)據(jù)建立變分貝葉斯-隱馬爾可夫模型。其中,取600組健康數(shù)據(jù)進(jìn)行模型訓(xùn)練,取1 260組數(shù)據(jù)進(jìn)行測試。
筆者將每組數(shù)據(jù)分成10段,每段256個(gè)點(diǎn),得到的軸承的對數(shù)似然曲線如圖7所示。
圖7 似然對數(shù)曲線
從圖7的實(shí)驗(yàn)臺軸承全壽命周期曲線可知:
(1)前400組數(shù)據(jù)的DI值較大,表明軸承處于健康狀態(tài);(2)第400組開始,數(shù)據(jù)的DI值開始降低,表明此時(shí)軸承開始劣化并保持劣化狀態(tài)繼續(xù)運(yùn)行;(3)第900組數(shù)據(jù),DI值進(jìn)一步降低,表明軸承狀態(tài)持續(xù)劣化;(4)從1 200組數(shù)據(jù)到運(yùn)行結(jié)束,軸承的DI值持續(xù)降低并出現(xiàn)最低點(diǎn),表明軸承加速劣化直至失效。
由此可見,本文中的實(shí)驗(yàn)方法對實(shí)測平臺數(shù)據(jù)仍然適用。
結(jié)合小波包-CEEMD特征提取和VB-HMM的優(yōu)點(diǎn),筆者提出了一種基于小波包-互補(bǔ)聚合經(jīng)驗(yàn)?zāi)B(tài)分解和變分貝葉斯-隱馬爾可夫模型的滾動(dòng)軸承健康狀態(tài)評價(jià)方法,建立了發(fā)電機(jī)組滾動(dòng)軸承部件的健康狀態(tài)評價(jià)模型;采用標(biāo)準(zhǔn)軸承數(shù)據(jù)和實(shí)驗(yàn)臺實(shí)測數(shù)據(jù)進(jìn)行了驗(yàn)證,很好地反映了軸承部件的性能退化過程,為制定運(yùn)行維護(hù)策略提供了參考。
研究結(jié)果表明:該模型輸出的對數(shù)似然概率值可作為指標(biāo)進(jìn)行狀態(tài)評價(jià);對數(shù)似然概率曲線的變化均能反映軸承的健康、劣化和失效的運(yùn)行狀態(tài),該實(shí)驗(yàn)方法具有很好的適用性。