基于混合高斯輸出貝葉斯信念網(wǎng)絡(luò)模型的設(shè)備退化狀態(tài)識(shí)別與剩余使用壽命預(yù)測(cè)方法研究

張星輝， 康建設(shè)， 趙勁松,2， 肖 雷， 曹端超， 劉 浩

(1.軍械工程學(xué)院，石家莊 050003；2.軍事交通學(xué)院，天津 300161；3.重慶大學(xué) 重慶 400030)

基于混合高斯輸出貝葉斯信念網(wǎng)絡(luò)模型的設(shè)備退化狀態(tài)識(shí)別與剩余使用壽命預(yù)測(cè)方法研究

張星輝1， 康建設(shè)1， 趙勁松1,2， 肖 雷3， 曹端超1， 劉 浩1

(1.軍械工程學(xué)院，石家莊 050003；2.軍事交通學(xué)院，天津 300161；3.重慶大學(xué) 重慶 400030)

提出了基于混合高斯輸出貝葉斯信念網(wǎng)絡(luò)模型的設(shè)備退化狀態(tài)識(shí)別與剩余使用壽命預(yù)測(cè)新方法，將變量消元和期望最大化算法相結(jié)合對(duì)模型進(jìn)行推理，應(yīng)用聚類評(píng)價(jià)指標(biāo)對(duì)狀態(tài)數(shù)進(jìn)行優(yōu)化，通過(guò)計(jì)算待識(shí)別特征向量的概率值來(lái)確定設(shè)備當(dāng)前的退化狀態(tài)，在退化狀態(tài)識(shí)別的基礎(chǔ)上，提出了剩余使用壽命預(yù)測(cè)方法。最后，分別應(yīng)用50組軸承全壽命仿真數(shù)據(jù)和3組軸承全壽命實(shí)驗(yàn)數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證。結(jié)果表明，該模型可有效地識(shí)別設(shè)備的退化狀態(tài)并對(duì)剩余使用壽命進(jìn)行預(yù)測(cè)。

混合高斯輸出貝葉斯信念網(wǎng)絡(luò)模型；退化狀態(tài)識(shí)別；剩余使用壽命預(yù)測(cè)；軸承

基于狀態(tài)的維修(Condition Based Maintenance，CBM)是指由安裝在設(shè)備外部的傳感器獲取設(shè)備運(yùn)行時(shí)的狀態(tài)信息后，通過(guò)對(duì)狀態(tài)信息進(jìn)行分析，判斷設(shè)備當(dāng)前所處的退化狀態(tài)并對(duì)設(shè)備的剩余使用壽命(Remaining Useful Life，RUL)進(jìn)行預(yù)測(cè)，在此基礎(chǔ)上對(duì)設(shè)備的維修活動(dòng)進(jìn)行規(guī)劃。目前，RUL預(yù)測(cè)方法主要有兩種：基于物理模型的方法和基于數(shù)據(jù)驅(qū)動(dòng)的方法。Bechhoefer等[1]應(yīng)用Paris法則構(gòu)建了裂紋與振動(dòng)監(jiān)測(cè)指標(biāo)之間的關(guān)系并對(duì)RUL進(jìn)行預(yù)測(cè)。Jin等[2]研究了衛(wèi)星動(dòng)量輪的故障機(jī)理并建立了RUL預(yù)測(cè)模型。對(duì)于復(fù)雜設(shè)備的部件而言，由于其故障機(jī)理非常復(fù)雜，很難建立有效的物理模型對(duì)RUL進(jìn)行預(yù)測(cè)。而基于數(shù)據(jù)驅(qū)動(dòng)的方法則可以直接利用狀態(tài)監(jiān)測(cè)數(shù)據(jù)對(duì)RUL進(jìn)行預(yù)測(cè)。其主要方法有：人工神經(jīng)網(wǎng)絡(luò)(Artificial Neural Network，ANN)[3-4]、隱馬爾可夫族模型(Hidden Markov Models，HMM)[5-11]、支持向量機(jī)(Support Vector Machine，SVM)[12-15]、比例風(fēng)險(xiǎn)模型(Proportional Hazards Model，PHM)[16-17]和狀態(tài)空間模型(State Space Model，SSM)[18-20]等。

如果按照退化狀態(tài)進(jìn)行分類，基于數(shù)據(jù)驅(qū)動(dòng)的方法可以分為兩類：連續(xù)退化狀態(tài)模型和離散退化狀態(tài)模型。對(duì)于連續(xù)退化狀態(tài)模型，其需要用最新的退化特征對(duì)模型參數(shù)不斷的更新。如果退化特征非平穩(wěn)且有較大的噪聲，那么模型的性能將會(huì)下降并導(dǎo)致虛警率升高。而離散退化狀態(tài)模型則能克服這些不足，對(duì)退化特征的平穩(wěn)性要求較低。與已有的離散退化狀態(tài)模型(HMM等)相比，貝葉斯信念網(wǎng)絡(luò)(Bayesian Belief Networks，BBN)模型[21-24]結(jié)構(gòu)更易于擴(kuò)展，且有較強(qiáng)的退化狀態(tài)識(shí)別能力。孟光磊等[25]構(gòu)建了混合貝葉斯網(wǎng)絡(luò)對(duì)空域目標(biāo)的威脅進(jìn)行評(píng)估。王利民等[26]提出了具有抗噪音能力的增量式混合貝葉斯網(wǎng)絡(luò)。王雙成等[27-28]對(duì)混合貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的學(xué)習(xí)方法進(jìn)行了研究。Tobon-Mejia等[11]應(yīng)用混合高斯輸出HMM(Mixture of Gaussian Hidden Markov Model, MoG-HMM)對(duì)設(shè)備的退化過(guò)程建模并對(duì)其RUL進(jìn)行預(yù)測(cè)。借鑒混合高斯輸出HMM的模型結(jié)構(gòu)，可以構(gòu)建與之有相同功能且結(jié)構(gòu)更為簡(jiǎn)化的混合高斯輸出貝葉斯信念網(wǎng)絡(luò)(Mixture of Gaussian Bayesian Belief Networks，MoG-BBN)。且實(shí)驗(yàn)分析證明MoG-BBN退化狀態(tài)識(shí)別效果要優(yōu)于MoG-HMM，因此該研究采用MoG-BBN對(duì)設(shè)備退化過(guò)程進(jìn)行建模。

退化狀態(tài)識(shí)別需要提取退化特征向量，該特征向量可由時(shí)域特征、頻域特征和時(shí)頻域特征中的一種或者多種共同組成。本文應(yīng)用小波包分解提取退化特征向量。目前，如何合理劃分設(shè)備的退化狀態(tài)數(shù)還需要進(jìn)行更為深入的研究，它是退化狀態(tài)識(shí)別和RUL預(yù)測(cè)乃至維修決策的前提。為此，本文應(yīng)用聚類指標(biāo)確定退化狀態(tài)數(shù)。然后應(yīng)用已知退化狀態(tài)序列的設(shè)備歷史數(shù)據(jù)對(duì)模型進(jìn)行訓(xùn)練。模型訓(xùn)練好后即可對(duì)新設(shè)備的退化狀態(tài)進(jìn)行識(shí)別并結(jié)合預(yù)測(cè)模型對(duì)RUL進(jìn)行預(yù)測(cè)。最后，仿真數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了方法的有效性。

1 MoG-BBN基本原理

圖1 貝葉斯信念網(wǎng)絡(luò)Fig.1 Bayesian Belief Network

BBN是一種概率推理網(wǎng)絡(luò)，它以圖形節(jié)點(diǎn)表示隨機(jī)變量，節(jié)點(diǎn)之間的有向箭頭表示隨機(jī)變量之間的因果關(guān)系。BBN也是一種在不確定條件下強(qiáng)有力的知識(shí)表達(dá)和推理方法。BBN中節(jié)點(diǎn)所代表隨機(jī)變量的取值可以是連續(xù)的，也可以是離散的。連續(xù)變量可以服從任意分布，離散變量的取值可以是兩個(gè)或多個(gè)。本文構(gòu)建的MoG-BBN模型，其網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示：

MoG-BBN模型參數(shù)描述如下：

(1)X表示隱狀態(tài)，代表設(shè)備的退化狀態(tài)，取值為1,2,…,a；

(2)M表示混合數(shù)，代表組成設(shè)備某個(gè)退化狀態(tài)的高斯函數(shù)的個(gè)數(shù)，取值范圍為1,2,…,b；

在圖1所示的貝葉斯網(wǎng)絡(luò)中，X是根節(jié)點(diǎn)(M和Y的父節(jié)點(diǎn))，M是中間節(jié)點(diǎn)(X的子節(jié)點(diǎn)，Y的父節(jié)點(diǎn))，Y是葉節(jié)點(diǎn)(X和M的子節(jié)點(diǎn))。

該模型和HMM、SVM類似，都是具有模式識(shí)別功能的模型。利用已知類別的數(shù)據(jù)({yi,xi}，i=1,2,…,n)對(duì)模型進(jìn)行訓(xùn)練，再利用該模型對(duì)未知類別的數(shù)據(jù)(yj，j=1,2,…,n)進(jìn)行分類。在退化狀態(tài)識(shí)別中，X表示設(shè)備的退化狀態(tài)，Y表示經(jīng)過(guò)信號(hào)處理得到的反映某時(shí)刻設(shè)備退化的特征向量，y表示該隨機(jī)變量的值。識(shí)別過(guò)程即求取P(X|Y)的過(guò)程，其可以通過(guò)VE-EM算法實(shí)現(xiàn)。

2 VE-EM算法

2.1 變量消元(Variable Elimination，VE)算法

根據(jù)貝葉斯信念網(wǎng)鏈規(guī)則得：

(1)

根據(jù)條件概率公式得：

(2)

式(2)中P(X,Y=y)可表示為

(3)

式(2)中P(Y=y)可表示為

(4)

結(jié)合式(1)、(3)、(4)，式(2)可進(jìn)一步表示為

(5)

2.2 期望最大化(ExpectationMaximization，EM)算法

假設(shè)訓(xùn)練數(shù)據(jù)集為{y(1), …,y(n)}，n表示樣本數(shù)。由式(5)可以最終求得數(shù)據(jù)集所對(duì)應(yīng)的退化狀態(tài)序列，因此必須求得P(X)，P(M|X)和P(Y|X,M)，分別對(duì)應(yīng)四個(gè)參數(shù)：

(1)π：隱狀態(tài)(設(shè)備退化狀態(tài))概率分布，πi=P(X=i)，i∈{1, 2, …,a}；

(2)C：隱狀態(tài)混合系數(shù)，Cij=P(M=j|X=i)，i∈{1, 2, …,a}，j∈{1, 2, …,b}；

(3)μ：隱狀態(tài)產(chǎn)生高斯分布的均值；

(4) ∑：隱狀態(tài)產(chǎn)生高斯分布的方差。

這四個(gè)參數(shù)可以表示為：λ=(π,C,μ, ∑)。

(6)

s.t.

(7)

其中，為引入q(X)的未知分布，并給出它的約束條件。根據(jù)條件概率公式和貝葉斯公式，(7)式可表示為[29]：

(8)

其中的兩項(xiàng)分別為

(9)

(10)

根據(jù)Jensen不等式，(10)式可表示為

(11)

Q(λt,λt-1)+const

(12)

其中：

(13)

(14)

E步：

對(duì)所有的(i,j,l)

(15)

i∈{1, 2, …,a}，j∈{1, 2, …,b}，l∈{1, 2, …,n}。

根據(jù)貝葉斯公式，式(15)可表示為

(16)

M步：更新參數(shù)

(17)

(18)

(19)

(20)

3 基于MoG-BBN的退化識(shí)別與RUL預(yù)測(cè)

設(shè)備從正常到失效是一個(gè)逐漸退化的過(guò)程，即存在“潛在故障—功能故障”間隔。該間隔也稱為P-F間隔，如圖2所示?！癙”點(diǎn)是潛在故障點(diǎn)，是能夠發(fā)現(xiàn)設(shè)備缺陷的點(diǎn)，它通常會(huì)以更快的速度退化到功能故障點(diǎn)“F”點(diǎn)。

圖2 P-F間隔期Fig.2 P-F interval

設(shè)備的全壽命過(guò)程可分為若干個(gè)退化狀態(tài)，有些退化狀態(tài)處于潛在故障點(diǎn)“P”之前，退化指標(biāo)相比于設(shè)備剛開(kāi)始運(yùn)行時(shí)并沒(méi)有明顯的變化。因此，RUL預(yù)測(cè)不考慮“P”點(diǎn)之前的階段。對(duì)于CBM而言，為確保將設(shè)備的停機(jī)時(shí)間控制到最短，提高設(shè)備的可用度，最大化生產(chǎn)效益，需要RUL能夠滿足維修活動(dòng)規(guī)劃(備件訂購(gòu)、運(yùn)輸、維修隊(duì)實(shí)施更換等)。也就是說(shuō)，只要RUL能夠滿足這些要求，就說(shuō)明該預(yù)測(cè)方法是有效的。

利用MoG-BBN模型對(duì)設(shè)備退化狀態(tài)進(jìn)行識(shí)別和RUL預(yù)測(cè)包括：確定最優(yōu)退化狀態(tài)數(shù)、狀態(tài)識(shí)別和RUL預(yù)測(cè)。

3.1 退化狀態(tài)數(shù)優(yōu)化

設(shè)備退化狀態(tài)數(shù)的確定是利用MoG-BBN模型進(jìn)行退化狀態(tài)識(shí)別和RUL預(yù)測(cè)的前提。多數(shù)研究都是靠經(jīng)驗(yàn)來(lái)確定退化狀態(tài)數(shù)。董明根據(jù)油液的污染度將泵軸承分為四個(gè)退化狀態(tài)[7-8]。文獻(xiàn)[9]在應(yīng)用HMM預(yù)測(cè)時(shí)，假設(shè)設(shè)備的退化狀態(tài)數(shù)為兩個(gè)(好和壞)。退化狀態(tài)數(shù)的確定大致可分為三種方法：一是由專家根據(jù)經(jīng)驗(yàn)確定；二是利用交叉驗(yàn)證思想以分類器的分類錯(cuò)誤率最低來(lái)確定；三是應(yīng)用聚類評(píng)價(jià)指標(biāo)確定。專家經(jīng)驗(yàn)是通過(guò)大量實(shí)驗(yàn)積累的，對(duì)于高可靠、價(jià)格昂貴的部件并不適用。利用交叉驗(yàn)證的思想確定最優(yōu)退化狀態(tài)數(shù)需要對(duì)不同的退化狀態(tài)數(shù)都訓(xùn)練分類器并檢驗(yàn)其分類錯(cuò)誤率，計(jì)算時(shí)間長(zhǎng)，效率低[5]。而應(yīng)用聚類評(píng)價(jià)指標(biāo)確定最優(yōu)退化狀態(tài)數(shù)則簡(jiǎn)便易行，計(jì)算效率高。因此，本文采用第三種方法。

首先，假設(shè)聚類數(shù)(退化狀態(tài)數(shù))為a，a∈{2, 3…,amax}?？扇我膺x擇一種模糊聚類方法將提取的特征矩陣分為a類，計(jì)算各聚類指標(biāo)的值，每個(gè)聚類數(shù)值都有一組評(píng)價(jià)指標(biāo)值與其對(duì)應(yīng)，由此可以確定最優(yōu)的聚類數(shù)。聚類評(píng)價(jià)指標(biāo)分別為

分割系數(shù)[30](Partition Coefficient，PC)，其定義式為

(21)

分割指標(biāo)[31](Partition Index，SC)，其定義式為

(22)

分離指標(biāo)[31](Separation Index，SI)，其定義式為

(23)

Xie和Beni指標(biāo)[32](Xie and Beni’s Index，XB)：其定義式為

(24)

以上表達(dá)式中，uil表示第l個(gè)樣本屬于第i類的概率，y(l)表示第l個(gè)樣本，vi表示第i類的聚類中心，vk表示第k類的聚類中心，1≤i,k≤a，l表示樣本數(shù)目。對(duì)于最優(yōu)的聚類數(shù)目而言，PC取值越大越好，SC、SI和XB取值越小越好。

3.2 基于MoG-BBN的退化狀態(tài)識(shí)別

圖3 全壽命過(guò)程退化狀態(tài)識(shí)別Fig.3 Degradation state recognition of full life cycle

3.3 基于MoG-BBN的RUL預(yù)測(cè)

設(shè)備全壽命退化過(guò)程狀態(tài)識(shí)別結(jié)果表示為設(shè)備處于每個(gè)狀態(tài)的概率，如圖4所示。它表示當(dāng)設(shè)備最開(kāi)始運(yùn)行時(shí)，其處于正常狀態(tài)(第一退化狀態(tài))的概率為1，然后其概率逐漸下降，到達(dá)一定時(shí)刻時(shí)，第二退化狀態(tài)的概率逐漸上升，依次往復(fù)，當(dāng)最后一個(gè)退化狀態(tài)的概率為1時(shí)，設(shè)備失效。這些求取的退化狀態(tài)概率值是RUL預(yù)測(cè)的基礎(chǔ)。

圖4 全壽命過(guò)程退化狀態(tài)識(shí)別結(jié)果示意Fig.4 Degradation state identification results of full life cycle

t時(shí)刻設(shè)備的RUL可以通過(guò)式(25)計(jì)算得出：

(25)

其中：

(26)

式中，Dui表示由歷史數(shù)據(jù)得出的設(shè)備在第i個(gè)退化狀態(tài)的駐留時(shí)間，RULi表示由歷史數(shù)據(jù)得出的設(shè)備處于第i個(gè)退化狀態(tài)時(shí)的RUL。k為當(dāng)前狀態(tài)駐留時(shí)間系數(shù)，用于調(diào)整RUL預(yù)測(cè)的準(zhǔn)確度。

假設(shè)設(shè)備退化狀態(tài)的駐留時(shí)間服從正態(tài)分布，由歷史數(shù)據(jù)可以統(tǒng)計(jì)得出各個(gè)退化狀態(tài)駐留時(shí)間的均值和方差。最后，根據(jù)2Sigma準(zhǔn)則得出RUL預(yù)測(cè)值的范圍。待新設(shè)備從正常運(yùn)行到失效后，用所得退化狀態(tài)識(shí)別和RUL預(yù)測(cè)結(jié)果對(duì)歷史數(shù)據(jù)進(jìn)行更新，重新計(jì)算各退化狀態(tài)駐留時(shí)間的均值和方差。

基于MoG-BBN的退化狀態(tài)識(shí)別和RUL預(yù)測(cè)過(guò)程可用圖5表示。

圖5 基于MoG-BBN的設(shè)備退化狀態(tài)識(shí)別與RUL預(yù)測(cè)框架Fig.5 Framework of degradation states identification and RUL prediction based on MoG-BBN

4 仿真數(shù)據(jù)分析

4.1 數(shù)據(jù)仿真

根據(jù)文獻(xiàn)[33-34]中的方法對(duì)軸承外圈故障進(jìn)行仿真，所仿真軸承的幾何參數(shù)為：節(jié)徑39.039 8 mm、滾動(dòng)體9個(gè)、滾動(dòng)體直徑15.001 2 mm、接觸角0°。通過(guò)軸承故障特征頻率計(jì)算公式可求得軸承的頻率倍數(shù)如表1所示。

表1 軸承各部件故障頻率倍數(shù)

仿真軸承的轉(zhuǎn)速為1 800 r/min，采樣頻率為12 kHz，采樣持續(xù)時(shí)間10 s。

圖6所示為仿真信號(hào)的時(shí)域表示，圖7為該信號(hào)的包絡(luò)譜。

由于指數(shù)函數(shù)經(jīng)常用來(lái)描述機(jī)械產(chǎn)品的退化規(guī)律，可見(jiàn)文獻(xiàn)[35-37]。因此，為使得信號(hào)的峰值不斷增長(zhǎng)直至超過(guò)預(yù)先設(shè)定的閾值，可在仿真時(shí)，將軸承外圈故障信號(hào)乘以一個(gè)指數(shù)函數(shù)。假設(shè)軸承的壽命服從威布爾分布，指數(shù)函數(shù)可由式(27)確定。

eTi/βi=閾值

(27)

Ti表示軸承壽命，βi表示與軸承壽命相對(duì)應(yīng)的參數(shù)，當(dāng)閾值和軸承壽命已知時(shí)，可求得βi，從而進(jìn)行軸承全壽命數(shù)據(jù)的仿真。取威布爾分布的尺度參數(shù)和形狀參數(shù)為3 500和5，設(shè)定峰值20 m/s2為閾值。仿真后，得到50組軸承全壽命數(shù)據(jù)。圖8所示為仿真的前5組軸承全壽命數(shù)據(jù)的峰值。

圖6 軸承外圈故障仿真信號(hào)時(shí)域表示Fig．6Simulatedsignalofbearingwithanoutracedefect圖7 軸承外圈故障仿真信號(hào)包絡(luò)譜Fig．7Frequencyspectrumafterenvelopedemodulation圖8 前5組軸承全壽命數(shù)據(jù)的峰值Fig．8Peakvaluesof5simulatedbearingtest?to?failuredata

4.2 仿真數(shù)據(jù)退化狀態(tài)識(shí)別

(1) 狀態(tài)數(shù)優(yōu)化

用‘db4’小波對(duì)軸承仿真數(shù)據(jù)進(jìn)行小波包3層分解，提取各頻帶能量作為特征向量，從而得到每組數(shù)據(jù)全壽命過(guò)程特征矩陣F。對(duì)軸承1(第1組數(shù)據(jù))的特征矩陣F進(jìn)行聚類時(shí)，取聚類數(shù)2～7。分別計(jì)算PC、SC、SI、XB的值，結(jié)果如表2所示。從表中可以看出，

聚類數(shù)目為5時(shí)，PC值最大，SC、SI和XB值最小。因此將軸承的全壽命過(guò)程分為5個(gè)退化狀態(tài)。應(yīng)用K均值聚類算法將軸承仿真數(shù)據(jù)的特征矩陣聚為5類，獲得類別序號(hào)及聚類中心向量。

為了驗(yàn)證提出的聚類指標(biāo)確定退化狀態(tài)數(shù)方法的優(yōu)越性，與交叉驗(yàn)證方法進(jìn)行對(duì)比分析。按照文獻(xiàn)[5]中的步驟，可以用交叉驗(yàn)證方法求得軸承1每種退化狀態(tài)數(shù)相對(duì)應(yīng)的分類錯(cuò)誤率，如圖9所示。

表2 四種聚類評(píng)價(jià)指標(biāo)不同狀態(tài)數(shù)目劃分取值

圖9 軸承退化狀態(tài)數(shù)優(yōu)化Fig.9 Degradation state number optimization

從圖9可以看出，退化狀態(tài)數(shù)為5時(shí)，分類錯(cuò)誤率最低。和聚類評(píng)價(jià)指標(biāo)確定的退化狀態(tài)數(shù)一致。但交叉驗(yàn)證方法確定最優(yōu)狀態(tài)數(shù)所需的時(shí)間為1 113.54 s，而聚類評(píng)價(jià)指標(biāo)確定最優(yōu)狀態(tài)數(shù)僅需13.75 s，大大縮短了計(jì)算時(shí)間。對(duì)于待優(yōu)化退化狀態(tài)數(shù)范圍的選取，一般從2開(kāi)始，上限則需要根據(jù)具體設(shè)備的退化過(guò)程確定。對(duì)于本文中的軸承仿真數(shù)據(jù)，從其峰值觀測(cè)，其不可能有太多的退化狀態(tài)，因此將上限定為7。

(2) MoG-BBN模型訓(xùn)練和退化狀態(tài)識(shí)別

首先，用前40組軸承仿真數(shù)據(jù)中每個(gè)退化狀態(tài)所對(duì)應(yīng)的特征矩陣訓(xùn)練MoG-BBN模型，應(yīng)用訓(xùn)練好的模型對(duì)第41組軸承仿真數(shù)據(jù)進(jìn)行退化狀態(tài)識(shí)別后，該組數(shù)據(jù)就加入歷史數(shù)據(jù)，重新對(duì)MoG-BBN模型進(jìn)行訓(xùn)練。依次類推，對(duì)剩余的9組軸承仿真數(shù)據(jù)進(jìn)行狀態(tài)識(shí)別。圖10所示為對(duì)第41組數(shù)據(jù)的狀態(tài)識(shí)別結(jié)果。

圖10所示的結(jié)果與圖4相似，軸承開(kāi)始運(yùn)行時(shí)，處于正常狀態(tài)(第一退化狀態(tài))的概率為1， 然后逐漸下降，與此同時(shí)，下一退化狀態(tài)的概率逐漸上升，依次往復(fù)，直至軸承失效。

圖10 第41組仿真數(shù)據(jù)的退化狀態(tài)識(shí)別結(jié)果Fig.10 State recognition results of 41th bearing simulation data

4.3 仿真數(shù)據(jù)RUL預(yù)測(cè)

結(jié)合第41組仿真數(shù)據(jù)退化狀態(tài)識(shí)別的結(jié)果，可認(rèn)為“P”點(diǎn)即是退化狀態(tài)2起始時(shí)刻。因此，對(duì)第41-50組仿真組據(jù)的RUL預(yù)測(cè)都是從退化狀態(tài)2起始時(shí)刻開(kāi)始的。每組數(shù)據(jù)預(yù)測(cè)完成后，根據(jù)狀態(tài)識(shí)別結(jié)果，對(duì)退化狀態(tài)駐留時(shí)間的均值和方差重新估計(jì)，用于下一組數(shù)據(jù)的預(yù)測(cè)。以2Sigma作為RUL預(yù)測(cè)的范圍。預(yù)測(cè)結(jié)果如圖11所示。

圖11 第41～50組軸承全壽命仿真數(shù)據(jù)RUL預(yù)測(cè)結(jié)果Fig.11 RUL prediction results of the 41st～50th bearing test-to-failure simulation data

從圖11中可以看出，對(duì)于多數(shù)數(shù)據(jù)而言，其RUL預(yù)測(cè)結(jié)果與實(shí)際RUL都較為接近。由于第48組軸承的壽命與平均軸承壽命偏差太大，因此其每個(gè)退化狀態(tài)的駐留時(shí)間與平均駐留時(shí)間相差也較大。導(dǎo)致預(yù)測(cè)得到的RUL比實(shí)際RUL值小很多，但當(dāng)軸承將近失效時(shí)，預(yù)測(cè)的RUL與實(shí)際RUL值較為接近。像第48組軸承壽命與平均壽命偏差較大的情況屬于少數(shù)。因此，對(duì)于絕大多數(shù)軸承而言，該預(yù)測(cè)方法還是非常有效的。另外在制定維修決策時(shí)，充分考慮退化狀態(tài)識(shí)別結(jié)果可以避免RUL預(yù)測(cè)不準(zhǔn)確帶來(lái)的損失。

5 實(shí)驗(yàn)數(shù)據(jù)分析

5.1 預(yù)測(cè)模型實(shí)驗(yàn)驗(yàn)證

軸承全壽命數(shù)據(jù)來(lái)自IMS中心[38]。實(shí)驗(yàn)臺(tái)如圖12所示，1個(gè)軸上安裝了4套R(shí)exnord(萊克斯諾)ZA-2115雙列滾子軸承，每列滾子數(shù)量為16，滾子組節(jié)圓直徑為75.501 mm，滾子直徑為8.407 4 mm，接觸角為15.17°。軸的轉(zhuǎn)速為2 000 r/min，采樣頻率為20 kHz。每個(gè)軸承座都安裝2個(gè)PCB加速度傳感器，分別從水

平和垂直方向采集數(shù)據(jù)。振動(dòng)信號(hào)由NI公司DAQCard-6062E數(shù)據(jù)采集卡采集，每次采樣長(zhǎng)度為20 480個(gè)點(diǎn)。軸承2和軸承3額外增加27 210 N載荷。

圖12 軸承實(shí)驗(yàn)臺(tái)Fig.12 Bearing Test Rig

共進(jìn)行了3組實(shí)驗(yàn)，獲得了12個(gè)軸承的壽命數(shù)據(jù)?？梢詫⒌谝唤M實(shí)驗(yàn)的軸承編號(hào)為：軸承1、軸承2、軸承3和軸承4；第二組實(shí)驗(yàn)的軸承編號(hào)為：軸承5、軸承6、軸承7和軸承8；第三組實(shí)驗(yàn)的軸承編號(hào)為：軸承9、軸承10、軸承11和軸承12。12組軸承壽命過(guò)程的峰值如圖13所示。

圖13 12組軸承壽命過(guò)程的峰值Fig．13Peakvalueoftwelvebearing’slifeprocess圖14 軸承4RUL預(yù)測(cè)Fig．14RULpredictionresultsofbearing4圖15 軸承4均方根值Fig．15RootMeanSquarevalueofbearing4

假如設(shè)定軸承失效閾值為1.7 m/s2，從圖13中可以看出，只有軸承3、軸承4和軸承5的峰值超過(guò)了失效閾值。用軸承3和軸承5的數(shù)據(jù)對(duì)MoG-BBN訓(xùn)練，用軸承4的數(shù)據(jù)對(duì)模型進(jìn)行驗(yàn)證。為了提取有效的軸承退化特征，首先對(duì)數(shù)據(jù)進(jìn)行包絡(luò)分析，然后用‘db4’小波對(duì)包絡(luò)信號(hào)進(jìn)行小波包3層分解，以第三層分解系數(shù)的能量作為退化特征向量。應(yīng)用聚類評(píng)價(jià)指標(biāo)確定最優(yōu)退化狀態(tài)數(shù)為3。按照?qǐng)D5所示的流程可以對(duì)軸承4的RUL進(jìn)行預(yù)測(cè)。結(jié)果如圖14所示。

對(duì)比圖14和15可知，應(yīng)用MoG-BBN模型對(duì)軸承4的RUL預(yù)測(cè)是非常有效的。從時(shí)間點(diǎn)1 200/10 min開(kāi)始，軸承4開(kāi)始加速退化，其RUL預(yù)測(cè)值與真實(shí)值也開(kāi)始接近。

5.2 MoG-BBN與MoG-HMM模型對(duì)比分析

為了對(duì)MoG-BBN和MoG-HMM的退化識(shí)別能力進(jìn)行對(duì)比分析，對(duì)凱斯西楚大學(xué)公開(kāi)軸承實(shí)驗(yàn)數(shù)據(jù)[39]進(jìn)行分析。選擇驅(qū)動(dòng)端同工況下不同程度的內(nèi)圈退化數(shù)據(jù)對(duì)模型進(jìn)行對(duì)比，數(shù)據(jù)的采樣頻率為12 kHz。軸承型號(hào)為6205-2RS JEM SKF，損傷直徑分別為0.177 8、0.355 6和0.533 4 mm，分別對(duì)應(yīng)軸承的三種退化狀態(tài)。工況為1 730 r/min，3馬力。將每個(gè)退化狀態(tài)的信號(hào)分為24段(每段0.42 s采樣)，應(yīng)用‘db4’小波對(duì)信號(hào)進(jìn)行小波包3層分解，以第三層分解系數(shù)的能量作為退化特征向量。每個(gè)退化狀態(tài)用4組數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，其它20組數(shù)據(jù)作為測(cè)試數(shù)據(jù)。對(duì)所有的數(shù)據(jù)的數(shù)據(jù)進(jìn)行歸一化。MoG-HMM應(yīng)用PMTK[40]工具箱實(shí)現(xiàn)。MoG-HMM和MoG-BBN的初始參數(shù)都隨機(jī)選擇。混合數(shù)設(shè)為2。退化識(shí)別結(jié)果如圖16和17所示。P0.177 8、P0.355 6和P0.533 4分別代表三種退化狀態(tài)的概率，概率最高的即為識(shí)別結(jié)果。

圖16 MoG-HMM識(shí)別結(jié)果Fig.16 Degradation state identification using MoG-HMM

圖17 MoG-BBN識(shí)別結(jié)果Fig.17 Degradation state identification using MoG-BBN

對(duì)比分析圖14和圖15可知，MoG-BBN退化識(shí)別效果要優(yōu)于MoG-HMM。對(duì)其它實(shí)驗(yàn)數(shù)據(jù)的分析得到同樣的結(jié)果。通常，MoG-HMM需要更多的退化特征才能達(dá)到較好的效果。而當(dāng)提取的退化特征有限時(shí)，MoG-HMM的識(shí)別效果就會(huì)較差。與此相的是MoG-BBN即使是在退化特征有限的條件下，仍然能夠達(dá)到較好的退化識(shí)別效果。因此，充分說(shuō)明了應(yīng)用MoG-BBN進(jìn)行退化狀態(tài)識(shí)別和RUL預(yù)測(cè)的有效性。由于篇幅關(guān)系不能和所有的方法一一比較。

6 結(jié) 論

(1) 軸承全壽命仿真數(shù)據(jù)和全壽命實(shí)驗(yàn)數(shù)據(jù)較好地驗(yàn)證了MoG-BBN模型退化狀態(tài)識(shí)別和結(jié)合預(yù)測(cè)模型進(jìn)行RUL預(yù)測(cè)的有效性。在具體實(shí)施過(guò)程中，當(dāng)前狀態(tài)駐留時(shí)間系數(shù)k的選取對(duì)RUL預(yù)測(cè)結(jié)果有較大影響。在第41-50組軸承全壽命仿真數(shù)據(jù)RUL預(yù)測(cè)中，k的取值分別為0.3、0.9、0.4、0.3、0.3、0.3、0.5、0、0.4、0.5；軸承全壽命實(shí)驗(yàn)數(shù)據(jù)RUL預(yù)測(cè)中，k的取值為0.7。在下一步研究中，如何選取最優(yōu)的系數(shù)k以提高平均預(yù)測(cè)準(zhǔn)確率是主要研究?jī)?nèi)容。對(duì)于混合數(shù)M的選取，目前還沒(méi)有有效的方法。在下一步的研究中，將以退化識(shí)別結(jié)果為目標(biāo)對(duì)混合數(shù)進(jìn)行優(yōu)化。在工程實(shí)際中，同一種型號(hào)的軸承其失效模式不盡相同，因此其退化過(guò)程也不相同，論文中的退化識(shí)別和預(yù)測(cè)方法需要結(jié)合實(shí)際做進(jìn)一步的改進(jìn)。

(2) 基于狀態(tài)的維修決策優(yōu)化主要是在RUL預(yù)測(cè)的基礎(chǔ)上對(duì)備件訂購(gòu)、運(yùn)輸、庫(kù)存及維修隊(duì)派遣等一系列活動(dòng)進(jìn)行優(yōu)化，旨在降低設(shè)備停機(jī)時(shí)間，提高設(shè)備的可用度，降低維修費(fèi)用，最大化生產(chǎn)效益。因此，只要在RUL預(yù)測(cè)值與實(shí)際值相近的時(shí)刻，還能留有足夠的時(shí)間用于維修決策優(yōu)化活動(dòng)就說(shuō)明該方法是有實(shí)際意義的。從軸承實(shí)驗(yàn)數(shù)據(jù)研究結(jié)果可知，預(yù)測(cè)值與實(shí)際值相近的時(shí)間約占全壽命的28%。由文獻(xiàn)[41]可知，風(fēng)電機(jī)組的平均故障間隔時(shí)間約為10 000～15 000小時(shí)。那么，按照MoG-BBN預(yù)測(cè)的RUL有效時(shí)間可以達(dá)到2 800～4 200小時(shí)，該段時(shí)間內(nèi)能夠保障企業(yè)有充足時(shí)間對(duì)各種維修活動(dòng)進(jìn)行安排。案例研究發(fā)現(xiàn)退化特征的提取對(duì)于RUL預(yù)測(cè)的有效時(shí)間(P-F間隔)有重要影響，如果提取的特征能夠更早的具有明顯的增長(zhǎng)趨勢(shì)(也就是更早地發(fā)現(xiàn)“P”點(diǎn))，那么，RUL的有效預(yù)測(cè)區(qū)間也能更大，就有更加充足的時(shí)間對(duì)維修活動(dòng)進(jìn)行安排。因此，提取更為有效的退化特征是下一步的主要研究?jī)?nèi)容。

本文提出了基于MoG-BBN模型的設(shè)備退化狀態(tài)識(shí)別與RUL預(yù)測(cè)方法，建立了基于VE-EM的模型推理算法。50組軸承全壽命仿真數(shù)據(jù)和3組軸承全壽命實(shí)驗(yàn)數(shù)據(jù)較好地驗(yàn)證了論文提出方法的有效性，為設(shè)備的健康管理提供了科學(xué)依據(jù)，對(duì)于提高設(shè)備的可用度，最大化生產(chǎn)效益具有重要的現(xiàn)實(shí)意義。

[ 1 ] Bechhoefer E, Bernhard A, He D. Use of paris law for prediction of component remaining life[C]. IEEE Aerospace Conference, 2008, 1-9.

[ 2 ] Jin G, Matthews D, Fan Y, et al. Physics of failure-based degradation modeling and lifetime prediction of the momentum wheel in a dynamic covariate environment[J]. Engineering Failure Analysis, 2013, 28(3):222-240.

[ 3 ] Gebraeel N, Lawley M, Liu R, et al. Residual life predictions from vibration-based degradation signals: a neural network approach[J]. IEEE Transactions on Industrial Electronics, 2004, 51(3):694-700.

[ 4 ] Tian Z G, Wong L, Safaei N. A neural network approach for remaining useful life prediction utilizing both failure and suspension histories[J]. Mechanical Systems and Signal Processing, 2010, 24(5):1542-1555.

[ 5 ] 滕紅智,趙建民,賈希勝,等. 基于CHMM的齒輪箱狀態(tài)識(shí)別研究[J]. 振動(dòng)與沖擊, 2012, 31(5):92-96. TENG Hong-zhi,ZHAO Jian-min, JIA Xi-sheng, et al. Gearbox state recognition based on Continuous Hidden Markov Model[J]. Journal of Vibration and Shock, 2012, 31(5):92-96.

[ 6 ] 胡海峰,安茂春,秦國(guó)軍,等. 基于隱半Markov模型的故障診斷和故障預(yù)測(cè)方法研究[J]. 兵工學(xué)報(bào),2009,30(1):69-75. HU Hai-feng, AN Mao-chun, QIN Guo-jun，et al. Study on fault diagnosis and prognosis methods based on hidden semi-markov model[J]. Acta Armamentarii, 2009,30(1):69-75.[ 7 ] Dong M, He D. A segmental Hidden Semi-Markov Model (HSMM)-based diagnostics and prognostics framework and methodology[J]. Mechanical Systems and Signal Processing, 2007,21(5):2248-2266.

[ 8 ] Peng Y,Dong M. A prognosis method using age-dependent hidden semi-Markov model for equipment health prediction[J]. Mechanical Systems and Signal Processing, 2011, 25(1):237-252.

[ 9 ] Zhou Z J, Hu C H, Xu D L,et al. A model for real-time failure prognosis based on hidden Markov model and belief rule base[J]. European Journal of Operational Research, 2010, 207(1):269-283.

[10] Baruah P,Chinnam R B. HMMs for diagnostics and prognostics in machining processes[J]. International Journal of Production Research, 2005, 43(6):1275-1293.

[11] Tobon-Mejia D A, Medjaher K, Zerhouni N, et al. A data-driven failure prognostics method based on mixture of gaussians hidden Markov models[J]. IEEE Transactions on Reliability, 2012, 61(2):491-503.

[12] Kim H E, Tan A C C, Mathew J, et al. Bearing fault prognosis based on health state probability estimation[J]. Expert Systems with Applications, 2012,39(5):5200-5213.

[13] Caesarendra W, Widodo A, Thom P H, et al. Combined probability approach and indirect data-driven method for bearing degradation prognostics[J]. IEEE Transactions on Reliability, 2011, 60(1):14-20.

[14] Kang J S, Zhang X H, Zhao J M, et al. Gearbox fault diagnosis based on CHMM and SVM[C]. Proceedings of the 2012 International Conference on Quality, Reliability, Risk, Maintenance, and Safety Engineering, 703-708.

[15] Widodo A, Yang B S. Machine health prognostics using survival probability and support vector machine[J]. Expert Systems with Applications, 2011, 38(7):8430-8437.

[16] Carr M J,Wang W. A case comparision of a proportional hazards model and a stochastic filter for condition-based maintenance applications using oil-based condition monitoring information[J]. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, 2008, 222(1):47-55.

[17] Samrout M, Chatelet E, Kouta R, et al. Optimization of maintenance policy using the proportional hazard model[J]. Reliability Engineering and System Safety, 2009, 94(1):44-52.

[18] Sun J Z, Zuo H F, Wang W B, et al. Application of a state space modeling technique to system prognostics based on a health index for condition-based maintenance[J]. Mechanical Systems and Signal Processing, 2012, 28(4):585-596.

[19] Wang W B. A prognosis model for wear prediction based oil-based monitoring[J]. Journal of the Operational Research Society, 2007, 58(7):887-893.

[20] Wang W B. Overview of a semi-stochastic filtering approach for residual life estimation with applications in condition based maintenance[J]. Proceedings of the Institution of Mechanical Engineers, Part O: Journal of Risk and Reliability, 2011, 225(2), 185-197.

[21] Xu B G. Intelligent fault inference for rotating flexible rotors using Bayesian belief network[J]. Expert Systems with Applications, 2012,39(1):816-822.

[22] Dey S, Stori J A. A Bayesian network approach to root cause diagnosis of process variations[J]. International Journal of Machine Tools and Manufacture, 2005,45(1):75-91.

[23] Verron S, Li J, Tiplica T. Fault detection and isolation of faults in a multivariate process with Bayesian network[J]. Journal of Process Control, 2010,20(8):902-911.

[24] Sahin F, Yavuz M C, Arnavut Z, et al. Fault diagnosis for airplane engines using Bayesian networks and distributed particle swarm optimization[J]. Parallel Computing, 2007,33(2):124-143.

[25] 孟光磊,龔光紅.基于混合貝葉斯網(wǎng)的空域目標(biāo)威脅評(píng)估方法[J].系統(tǒng)工程與電子技術(shù),2010,32(11):2398-2401. MENG Guang-lei, GONG Guang-hong. Threat assessment of aerial targets based on hybrid Bayesian network[J]. Systems Engineering and Electronics, 2010,32(11):2398-2401.

[26] 王利民,苑森淼. 具有抗噪音能力的增量式混合貝葉斯網(wǎng)絡(luò)[J]. 儀器儀表學(xué)報(bào),2005,26(3):221-225. WANG Li-min, YUAN Sen-miao. Incremental hybrid bayesian network with noise resistance ability[J]. Chinese Journal of Scientific Instrument, 2005,26(3):221-225.

[27] 王雙成,李小琳,侯彩虹.用于因果分析的混合貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)學(xué)習(xí)[J].智能系統(tǒng)學(xué)報(bào),2007,2(6):82-89. WANG Shuang-cheng, LI Xiao-lin, HOU Cai-hong. Learning in a hybrid Bayesian network structure for causal analysis[J]. CAAIT Transactions on Intelligent Systems, 2007,2(6):82-89.

[28] 王雙成.混合貝葉斯網(wǎng)絡(luò)隱藏變量學(xué)習(xí)研究[J]. 計(jì)算機(jī)學(xué)報(bào),2005,28(6):1564-1569. WANG Shuang-cheng. Research on learning the hidden variables of hybrid Bayesian network[J]. Chinese Journal of Computers, 2005,28(6):1564-1569.

[29] Bishop C M. Pattern recognition and machine learning[M]. Springer, 2006:435-455.

[30] Bezdek J C. Pattern recognition with fuzzy objective function algorithms[M]. Plenum Press, 1981.

[31] Bensaid A M, Hall L O, Bezdek J C, et al. Validity-guided (re) clustering with applications to image segmentation[J]. IEEE Transactions on Fuzzy Systems, 1996, 4(2):112-123.

[32] Xie X L, Beni G. A validity measure for fuzzy clustering[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1991, 13(8):841-847.

[33] McFadden P D, Smith J D. Model for the vibration produced by a single point defect in rolling element bearing[J]. Journal of Sound and Vibration, 1984, 96(1):69-82.

[34] Wang Y F, Kootsookos P J. Modelling of low shaft speed bearing faults for condition monitoring[J]. Mechanical System and Signal Processing, 1998, 12(3): 415-426.

[35] Shao Y, Nezu K. Prognosis of remaining bearing life using neural networks[C]. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering, 2000: 217-230.

[36] Gebraeel N, Lawley M, Liu R, et al. Residual life predictions from vibration-based degradation signals: A neural network approach[J]. IEEE Transactions on Industrial Electronics, 2004, 51(3):694-700.

[37] Heng A S. Intelligent prognostics of machinery health utilising suspended condition monitoring data[D]. Thesis of Queenland University of Technology, Australia,2009.

[38] NASA Prognostic Data Repository: Bearing Data Set NSF I/UCRC Center for Intelligent Maintenance Systems[DB/OL]. 2010, Available at http://ti.arc.nasa.gov/tech/dash/pcoe/prognostic-data-repository.

[39] Bearing Data Center (BDC), Case Western Reserve University Bearing Data[DB/OL]. Available at http://csegroups.case.edu/bearingdatacenter/pages/welcome-case-western-reserve-university-bearing-data-center-website,2013.

[40] PMTK: Avaiable at http://code.google.com/p/pmtk3/ (2013).

[41] Tavner P J, Xiang J, Spinato F. Reliability analysis of wind turbines[J]. Wind Energy, 2007, 10(1):1-18.

Equipment degradation state identification and residual life prediction based on MoG-BBN

ZHANG Xing-hui1, KANG Jian-she1, ZHAO Jin-song1,2, XIAO Lei3, CAO Duan-chao1, LIU Hao1

(1.Ordnance Engineering College, Shijiazhuang 050003, China; 2.Military Transportation College, Tianjin 300161, China; 3. Chongqing University, Chongqing 400030, China)

A new approach for equipment health state identification and residual life prediction based on mixture of Gaussian Bayesian belief network (MoG-BBN) was presented. The inference algorithm was established based on variable elimination (VE) and expectation maximization (EM). State number was optimized based on cluster validity indexes. The equipment degradation state was determined through calculating the probability of eigenvectors to be identified. Then, the residual life prediction method was presented based on identifying the degradation state. Finally, 50 bearings whole life simulation data and 3 bearings whole life test data were used to demonstrate the proposed methods. The results showed that the proposed method can be used to identify degradation states of an equipment and predict its residual life effectively.

mixture of Gaussian Bayesian belief network; health state identification; residual life prediction; bearing

國(guó)家自然科學(xué)基金(51035008);軍隊(duì)“十二五”武器裝備預(yù)先研究項(xiàng)目(51327020101，51327020304)

2013-02-05 修改稿收到日期：2013-05-17

張星輝 男，博士生，講師，1984年6月生

TH165.3

A

10.13465/j.cnki.jvs.2014.08.030