基于TVAR-HMM 的滾動(dòng)軸承故障診斷

王國(guó)鋒，李玉波，秦旭達(dá)，喻 秀，李啟銘

(天津大學(xué)機(jī)械工程學(xué)院，天津 300072)

機(jī)械故障診斷在現(xiàn)代工業(yè)生產(chǎn)中起著越來(lái)越重要的作用．對(duì)反映機(jī)械設(shè)備狀況的振動(dòng)信號(hào)進(jìn)行信號(hào)分析及特征提取，然后利用數(shù)學(xué)模型進(jìn)行故障識(shí)別是目前故障診斷的主要方法．

滾動(dòng)軸承是機(jī)械設(shè)備中故障常發(fā)部件，多變的工況條件以及設(shè)備自身的非線性，使反映軸承故障的動(dòng)態(tài)信號(hào)表現(xiàn)出較強(qiáng)的非平穩(wěn)特性．目前信號(hào)主要的時(shí)頻分析方法有短時(shí)傅里葉變換、Wigner-Ville 分布、小波尺度分析等．這些方法的時(shí)頻分辨率都要受到不確定性原理的限制，因此無(wú)法同時(shí)在時(shí)域和頻域均得到很好的分辨率[1]．

故障識(shí)別模型以神經(jīng)網(wǎng)絡(luò)模型應(yīng)用最廣，但利用神經(jīng)網(wǎng)絡(luò)模型進(jìn)行故障識(shí)別需要大量的訓(xùn)練樣本，在有限的訓(xùn)練樣本情況下，泛化能力差，且優(yōu)化過(guò)程存在陷入局部極值的可能[2]；另外，大部分神經(jīng)網(wǎng)絡(luò)模型處理的是靜態(tài)模式分類問(wèn)題，而旋轉(zhuǎn)機(jī)械的故障行為是一個(gè)動(dòng)態(tài)變化的過(guò)程．

為此，筆者提出了應(yīng)用時(shí)變自回歸(time-varying autoregressive ，TVAR) 和 隱 馬 爾 科 夫 模 型 (hidden Markov model，HMM)的滾動(dòng)軸承故障診斷方法．時(shí)變自回歸模型是一種基于參數(shù)模型的時(shí)頻分析方法．它是在 AR(autoregressive)建模的基礎(chǔ)上發(fā)展起來(lái)的，TVAR 方法得出的時(shí)頻譜圖具有分辨率高、無(wú)交叉干擾項(xiàng)以及計(jì)算速度快等優(yōu)點(diǎn)，非常適合處理非平穩(wěn)信號(hào)，已在故障診斷領(lǐng)域有所應(yīng)用[3-5]．隱馬爾科夫模型是信號(hào)動(dòng)態(tài)時(shí)間序列統(tǒng)計(jì)模型，能夠?qū)σ粋€(gè)時(shí)間跨度上的信息進(jìn)行統(tǒng)計(jì)建模和分類，具有嚴(yán)謹(jǐn)?shù)臄?shù)據(jù)結(jié)構(gòu)和可靠的計(jì)算性能．它能通過(guò)較少的樣本訓(xùn)練出可靠的模型，并按模式匹配原理，尋找與未知信號(hào)最相似的模式作為識(shí)別結(jié)果，現(xiàn)已成為語(yǔ)音識(shí)別的主流技術(shù)[6]．傳統(tǒng)故障診斷方法一直停留在靜態(tài)分析基礎(chǔ)上，忽略了故障變化發(fā)展的動(dòng)態(tài)信息，而 HMM 用于故障診斷，可對(duì)故障變化的動(dòng)態(tài)信息進(jìn)行統(tǒng)計(jì)分析，能夠盡早地發(fā)現(xiàn)故障發(fā)展的跡象．與神經(jīng)網(wǎng)絡(luò)相比，保留了更多訓(xùn)練數(shù)據(jù)的統(tǒng)計(jì)信息，具有更高的識(shí)別率和魯棒性[7]．

筆者結(jié)合了兩種模型的優(yōu)點(diǎn)，研究了非平穩(wěn)信號(hào)的時(shí)變參數(shù)建模以及隱馬爾科夫建模故障識(shí)別過(guò)程，并通過(guò)實(shí)驗(yàn)驗(yàn)證該方法在滾動(dòng)軸承故障診斷中的可靠性．

1 時(shí)變自回歸模型

{x(t),t=1,2,…,N }是一非平穩(wěn)時(shí)間序列，TVAR是利用先前值 x (t ? i)來(lái)估計(jì)現(xiàn)在值 x(t)，可寫(xiě)成

式中：x (t)為信號(hào)的采集值；ai(t)(i=1,2,…,p)為時(shí)變參數(shù)；p 為模型階數(shù)；e(t)為模型殘差，服從零均值正態(tài)分布．

將時(shí)變參數(shù){ai(t),i=1,2,…,p}作為一組基時(shí)間函數(shù) gj(t )的線性組合，即

式中：m 為擴(kuò)展維數(shù)；aij為組合的時(shí)不變權(quán)值．

將式(2)代入式(1)，得

式中

這樣模型時(shí)變參數(shù)的估計(jì)問(wèn)題就變成常參數(shù) aij的估計(jì)，從而把一個(gè)線性非平穩(wěn)問(wèn)題轉(zhuǎn)化為一個(gè)線性時(shí)不變問(wèn)題．

文獻(xiàn)[8]闡述了利用最小二乘法估計(jì)模型參數(shù) A的過(guò)程，在估計(jì)出參數(shù) A 之后就可以根據(jù)式(2)計(jì)算出時(shí)變參數(shù) ai(t)(i =1,2,…,p)，繼而求解模型方差，即

瞬時(shí)頻率能夠通過(guò) TVAR 參數(shù) ai(t )提?。酶鱾€(gè)時(shí)刻的參數(shù)依次計(jì)算出每個(gè)時(shí)刻的參數(shù)譜圖，疊加在一起就可以構(gòu)造一個(gè)三維時(shí)頻譜圖，因此TVAR 譜函數(shù)定義為

2 隱馬爾科夫模型

HMM 是一個(gè)雙內(nèi)嵌式隨機(jī)過(guò)程，即 HMM 由兩個(gè)隨機(jī)過(guò)程組成：一是隱狀態(tài)轉(zhuǎn)移序列，它對(duì)應(yīng)一個(gè)單純的 Markov 過(guò)程；二是與隱狀態(tài)有關(guān)的觀測(cè)序列．其中隱狀態(tài)轉(zhuǎn)移序列是不可觀測(cè)的，只能通過(guò)另一個(gè)隨機(jī)過(guò)程的輸出觀測(cè)序列進(jìn)行推斷，所以稱之為隱馬爾科夫模型[9]．

HMM 可以記為 λ={Nλ，Mλ，πλ，Aλ，Bλ}或簡(jiǎn)寫(xiě)為λ={πλ，Aλ，Bλ}．其中 Nλ為模型中 Markov 鏈的狀態(tài)數(shù)目；Mλ為每個(gè)狀態(tài)對(duì)應(yīng)的可能的觀測(cè)值數(shù)目；πλ為初始概率分布矢量；Aλ為狀態(tài)轉(zhuǎn)移概率矩陣；Bλ為觀測(cè)值概率矩陣，對(duì)于連續(xù)的 HMM，Bλ是一組觀測(cè)值概率函數(shù)．

圖1為1 個(gè)三狀態(tài)的離散HMM 的例子，由圖可以看出，1 個(gè)HMM 可以由1 組狀態(tài)轉(zhuǎn)移概率以及與每一個(gè)狀態(tài)有關(guān)的觀測(cè)分布概率來(lái)定義．

圖1 三狀態(tài)HMM示意Fig.1 Diagram of three-state HMM

一個(gè)標(biāo)準(zhǔn)的HMM 模型需要解決概率計(jì)算、模型訓(xùn)練和隱狀態(tài)估計(jì)3 個(gè)基本問(wèn)題．

(1)概率計(jì)算．給定模型參數(shù) λ={πλ，Aλ，Bλ}、觀測(cè)序列 O={o1，o2，…，ot}，計(jì)算觀測(cè)序列 O 在給定模型下的概率，即計(jì)算似然 P(O︱λ)或者對(duì)數(shù)似然ln,P(O︱λ)，它表示參數(shù) λ 擬合數(shù)據(jù) O 的準(zhǔn)確程度．可通過(guò)前向-后向算法對(duì)HMM 進(jìn)行遍歷實(shí)現(xiàn)．

(2)模型訓(xùn)練，即模型參數(shù)估計(jì)問(wèn)題．也就是說(shuō)，對(duì)于初始模型和給定用于訓(xùn)練的觀測(cè)序列 O={o1，o2，…，ot}，如何調(diào)整模型 λ={πλ，Aλ，Bλ}的參數(shù)，使其能夠最好地?cái)M合觀測(cè)數(shù)據(jù) ? argmaxλ=f (O?)λ，這可以通過(guò)EM 算法或Baum-Welch 算法實(shí)現(xiàn)．

(3)隱狀態(tài)估計(jì)．給定模型參數(shù) λ={πλ，Aλ，Bλ}和觀測(cè)序列 O={o1，o2，…，ot}，基于某種最優(yōu)準(zhǔn)則估計(jì)最可能的隱狀態(tài)序列 Q={q1，q2，…，qt}，即估計(jì)產(chǎn)生觀測(cè)序列的最可能經(jīng)過(guò)路徑，這可以通過(guò) Viterbi 算法估計(jì)實(shí)現(xiàn)．

3 TVAR-HMM故障診斷步驟

圖2所示為利用 TVAR-HMM 進(jìn)行故障診斷的主要步驟．

圖2 TVAR-HMM故障診斷流程Fig.2 Flow chart of fault diagnosis based on TVAR-HMM

(1)對(duì)不同故障類型的振動(dòng)信號(hào)進(jìn)行 TVAR 建模，得到時(shí)頻譜圖．建模過(guò)程中，需要選定基函數(shù)gj(t)、模型階數(shù) p 和擴(kuò)展維數(shù)m ，文獻(xiàn)[10]詳細(xì)介紹了它們的選擇方法．

(2)基于能量法提取故障特征．對(duì)信號(hào) TVAR 建模后，故障信息完全體現(xiàn)在時(shí)頻譜圖上．用能量法提取故障特征可以更加有效地提取時(shí)頻圖上的故障信息．文獻(xiàn)[11]已研究了這種方法的有效性．其主要過(guò)程是：將時(shí)頻空間劃分為若干區(qū)域，分別計(jì)算其時(shí)頻局部能量，將時(shí)頻局部能量作為元素，組成特征向量．

(3)用特征向量訓(xùn)練 HMM 模型．特征向量為HMM 模型的觀測(cè)序列，故障類型為隱狀態(tài)，1 種故障類型訓(xùn)練1 個(gè)HMM 模型．

(4)利用訓(xùn)練好的模型對(duì)測(cè)試樣本進(jìn)行識(shí)別，即對(duì)于未知的觀測(cè)序列，分別計(jì)算在各個(gè)模型下的輸出概率，然后進(jìn)行比較，選取輸出概率最大的模型，即代表識(shí)別出的故障類型．

4 實(shí)驗(yàn)研究

為了驗(yàn)證本文方法在實(shí)際故障診斷中的應(yīng)用情況，筆者進(jìn)行了實(shí)驗(yàn)分析．實(shí)驗(yàn)在一個(gè)故障模擬平臺(tái)上進(jìn)行，包括驅(qū)動(dòng)電機(jī)、軸承座、扭矩加載器和采集設(shè)備，見(jiàn)圖 3．通過(guò)更換正常軸承和故障軸承進(jìn)行故障診斷實(shí)驗(yàn)．實(shí)驗(yàn)測(cè)試的軸承是深溝球軸承，型號(hào)為6205-2RS，品牌為SKF，其內(nèi)圈直徑為 25,mm，外圈直徑為52,mm，滾子直徑為7.94,mm．軸承故障用電火花加工而成，故障直徑分別為 0.177,8,mm、0.355,6,mm、0.533,4,mm,3 個(gè)等級(jí)，深度都為 0.279,4,mm．電機(jī)轉(zhuǎn)速為 1,797,r/min．振動(dòng)信號(hào)由加速度傳感器拾取，分別采集安裝無(wú)故障軸承、外圈故障軸承、內(nèi)圈故障軸承、滾子故障軸承時(shí)的振動(dòng)信號(hào)，各測(cè)得數(shù)據(jù) 50 組，每組2,048 點(diǎn)，采樣頻率為12,000,Hz．

圖3 實(shí)驗(yàn)裝置Fig.3 Experimental device

4.1 TVAR建模與特征提取

對(duì)無(wú)故障、內(nèi)圈故障、外圈故障、滾子故障4 種狀態(tài)信號(hào) TVAR 建模時(shí)，參照文獻(xiàn)[10]中的論述，基函數(shù)選擇傅里葉函數(shù)，采用AIC準(zhǔn)則確定模型階數(shù)p．在確定模型階數(shù)前必須先確定擴(kuò)展維數(shù) m，為了保證信號(hào)時(shí)頻譜的高分辨率和計(jì)算速度，一般模型的擴(kuò)展維數(shù)取為 8～12．本次實(shí)驗(yàn)選取擴(kuò)展維數(shù) m=10，依據(jù)AIC 準(zhǔn)則確定軸承無(wú)故障、內(nèi)圈故障、外圈故障、滾子故障狀態(tài)下模型階數(shù) p 分別為 11、15、15、12．

4 種類型的振動(dòng)信號(hào)TVAR 建模后，可以得到相應(yīng)的時(shí)頻譜，如圖 4(a)～(d)所示．由時(shí)頻譜中可以看出，不同故障類型軸承的振動(dòng)信號(hào)的頻率成分是不同的，并且軸承存在故障時(shí)，時(shí)頻譜中存在時(shí)變性與沖擊性．圖 4(e)為軸承外圈故障信號(hào)時(shí)域波形，圖4(f)為軸承外圈故障信號(hào)經(jīng)傅里葉變換后得到的頻譜．通過(guò)圖4(d)～(f)的比較可以明顯看出圖4(d)與4(f)中信號(hào)的頻率成分保持一致，圖 4(d)與 4(e)中頻率成分發(fā)生沖擊的時(shí)間及時(shí)間間隔基本一致．這說(shuō)明經(jīng) TVAR 建模得到的時(shí)頻譜具有較高的時(shí)頻分辨率．

基于能量法對(duì)時(shí)頻圖進(jìn)行特征提取時(shí)，考慮到速度與精度雙重因素，經(jīng)過(guò)多次實(shí)驗(yàn)分析，選擇將時(shí)頻譜分為16 部分(過(guò)多會(huì)影響HMM 訓(xùn)練速度；過(guò)少會(huì)影響測(cè)試的準(zhǔn)確度)，求每一部分能量的均值，將均值序列歸一化作為特征向量．通過(guò)特征提取，每種狀態(tài)得到 50 組 16 維的特征向量，4 種狀態(tài)共得到 4×50組特征向量．

圖4 4種狀態(tài)的時(shí)頻譜Fig.4 Time-frequency spectra of four states

4.2 HMM建模與故障識(shí)別

建模時(shí)，采用 4 個(gè)隱狀態(tài)代表無(wú)故障軸承、內(nèi)圈故障軸承、外圈故障軸承和滾子故障軸承，分別記作λ1、λ2、λ3、λ4．模型的觀測(cè)序列是 TVAR 建模后得到的 16 維特征向量．初始概率分布矢量、初始狀態(tài)轉(zhuǎn)移概率矩陣和初始觀測(cè)值概率矩陣均由 Matlab 中rand 函數(shù)隨機(jī)取得，然后歸一化．

選用120 組特征向量(每種狀態(tài)30 組)作為訓(xùn)練樣本，訓(xùn)練前需使用Lloyd 算法對(duì)特征向量進(jìn)行標(biāo)量量化．訓(xùn)練算法為Baum-Welch 算法． HMM 的訓(xùn)練過(guò)程中，隨著迭代次數(shù)的增加，最大對(duì)數(shù)似然估計(jì)值也在不斷增加，直到達(dá)到收斂誤差為止．訓(xùn)練結(jié)束后，得到4 個(gè)隱狀態(tài)對(duì)應(yīng)的4 個(gè)HMM 識(shí)別模型．圖5 給出了4 種狀態(tài)的HMM 訓(xùn)練曲線，所有狀態(tài)在迭代次數(shù)20 都達(dá)到了收斂，收斂速度快．

圖5 HMM訓(xùn)練曲線Fig.5 HMM training curves

在訓(xùn)練好各狀態(tài) HMM 以后，將剩余的 80 組特征向量(每種狀態(tài)20 組)作為測(cè)試樣本輸入到各狀態(tài)HMM 模型中(測(cè)試前需使用 Lloyd 算法對(duì)特征向量進(jìn)行標(biāo)量量化編碼)，在各個(gè)模型下能夠各自輸出 1個(gè)對(duì)數(shù)似然概率值．對(duì)數(shù)似然概率值反映了特征向量與各 HMM 的相似程度，對(duì)數(shù)似然概率值越大，特征向量越接近該狀態(tài) HMM，特征向量屬于使輸出對(duì)數(shù)似然概率值最大的模型所對(duì)應(yīng)的故障類型．識(shí)別算法選用 Viterbi 算法．表 1 給出了部分測(cè)試樣本的測(cè)試結(jié)果，通過(guò)得到的對(duì)數(shù)似然概率值，可以準(zhǔn)確識(shí)別出軸承的故障類型；對(duì)于同一種故障類型，故障程度越高，得到的對(duì)數(shù)似然概率值越大．

為了比較基于HMM 和基于BP 網(wǎng)絡(luò)、RBF 網(wǎng)絡(luò)故障診斷模型的識(shí)別效果，將本實(shí)驗(yàn)的數(shù)據(jù)分別用BP 網(wǎng)絡(luò)、RBF 網(wǎng)絡(luò)故障診斷模型來(lái)識(shí)別．網(wǎng)絡(luò)的層數(shù)設(shè)置為3 層，輸入層神經(jīng)元個(gè)數(shù)為16 個(gè)，輸出層神經(jīng)元個(gè)數(shù)為 4 個(gè)．采用與基于 HMM 故障診斷模型同樣的訓(xùn)練數(shù)據(jù)和測(cè)試數(shù)據(jù)分別訓(xùn)練和測(cè)試 BP 網(wǎng)絡(luò)和RBF 網(wǎng)絡(luò)，得出診斷結(jié)果．

圖6為本次實(shí)驗(yàn)中 3 種模型的故障識(shí)別正確率比較．由圖6 可見(jiàn)，基于 HMM 故障診斷模型的有效識(shí)別率明顯好于基于 BP 網(wǎng)絡(luò)或 RBF 網(wǎng)絡(luò)故障診斷模型．實(shí)驗(yàn)結(jié)果表明，基于 HMM 故障診斷模型可以成功地診斷旋轉(zhuǎn)機(jī)械的軸承故障，識(shí)別率較高，穩(wěn)定性好．

表1 部分測(cè)試結(jié)果Tab.1 Some testing results

圖6 HMM、BP、RBF識(shí)別有效性比較Fig.6 Diagnosis accuracies of HMM、BP、RBF

5 結(jié) 論

(1)TVAR 模型非常適合非平穩(wěn)信號(hào)的處理，時(shí)變參數(shù)體現(xiàn)了信號(hào)頻率的時(shí)變性和沖擊性．TVAR 建模分析得到的信號(hào)時(shí)頻譜的分辨率高，且無(wú)交叉項(xiàng)的干擾，計(jì)算速度快．

(2)HMM 模型可以成功地對(duì)旋轉(zhuǎn)機(jī)械的軸承故障進(jìn)行診斷，識(shí)別率高．由于其為動(dòng)態(tài)時(shí)間序列統(tǒng)計(jì)模型，非常適合識(shí)別信息量大、非平穩(wěn)性強(qiáng)、特征重復(fù)再現(xiàn)性不佳的故障信號(hào)，在旋轉(zhuǎn)機(jī)械故障診斷中有很好的應(yīng)用前景．

［1］李 強(qiáng). 非平穩(wěn)信號(hào)特征提取新方法和實(shí)用診斷技術(shù)的研究[D]. 天津：天津大學(xué)機(jī)械工程學(xué)院，2003.Li Qiang. Research on New Method of Non-Stationary Signal Feature Extraction and Practical Diagnostic Techniques[D]. Tianjin：School of Mechanical Engineering，Tianjin University，2003(in Chinese).

［2］王奉濤. 非平穩(wěn)信號(hào)故障特征提取與智能診斷方法的研究與應(yīng)用[D]. 大連：大連理工大學(xué)機(jī)械工程學(xué)院，2003.Wang Fengtao. Research on Non-Stationary Signal Feature Extraction and Intelligent Diagnostic Method[D].Dalian：School of Mechanical Engineering，Dalian University of Technology，2003(in Chinese).

［3］陳 慧，張 龍，李 嶷. 基于 TVAR 的非平穩(wěn)工況轉(zhuǎn)子故障診斷技術(shù)研究[J]. 華東交通大學(xué)學(xué)報(bào)，2006，23(5)：116-118.Chen Hui，Zhang Long，Li Yi. TVAR based fault diagnosis for rotation machine under non-stationary condition[J]. Journal of East China Jiaotong University，2006，23(5)：116-118(in Chinese).

［4］Wang Guofeng，Luo Zhigao. Fault identification and classification of rolling element bearing based on timevarying autoregressive spectrum [J]. Mechanical Systems and Signal Progressing，2008，22(4)：934-947.

［5］張 龍，熊?chē)?guó)良，柳和生，等. 基于時(shí)變自回歸模型與支持向量機(jī)的旋轉(zhuǎn)機(jī)械故障診斷方法[J]. 中國(guó)電機(jī)工程學(xué)報(bào)，2007，27(9)：99-103.Zhang Long，Xiong Guoliang，Liu Hesheng，et al. TVAR and SVM based fault diagnosis method for rotation machine [J]. Proceedings of the CSEE，2007，27(9)：99-103(in Chinese).

［6］Abramovich Y I，Spencer N K，Turley M D E. Timevarying autoregressive(TVAR)models for multiple radar observations [J]. IEEE Transactions on Signal Processing，2007，55(4)：1298-1311.

［7］Tai A H，Ching Wai-Ki，Chan L Y. Detection of machine failure：Hidden Markov model approach[J]. Computers& Industrial Engineering，2009，57(2)：608-619.

［8］Hall M G，Oppenheim A V，Willsky A S. Time varying parametric modeling of speech [J]. Signal Processing，1983，5(3)：267-285.

［9］Ocak H，Loparo K A. A new bearing fault detection and diagnosis scheme based on hidden Markov modeling of vibration signals [C]//IEEE International Conference on Acoustics，Speech，and Signal Processing. Salt Lake City，USA，2001，5：3141-3144.

［10］Zhan Y M，Jardine A K S. Adaptive autoregressive modeling of non-stationary vibration signals under distinct gear states(Part 1)：Modeling [J]. Journal of Sound and Vibration，2005，286：429-450.

［11］劉 偉. 基于能量的時(shí)頻特征提取研究[J]. 紅外與激光工程，2004，33(3)：296-299.Liu Wei. Research on the time-frequency feature extraction based on energy[J]. Infrared and Laser Project，2004，33(3)：296-299(in Chinese).