粒子濾波參數(shù)估計(jì)方法在齒輪箱剩余壽命預(yù)測(cè)中的應(yīng)用研究

孫 磊，賈云獻(xiàn)，蔡麗影，張星輝

(1.軍械工程學(xué)院裝備指揮與管理系，石家莊 050003;2.石家莊軍械技術(shù)研究所，石家莊 0500031)

設(shè)備剩余壽命預(yù)測(cè)是實(shí)現(xiàn)基于狀態(tài)維修(CBM)的關(guān)鍵技術(shù)之一，準(zhǔn)確的剩余壽命預(yù)測(cè)對(duì)維修決策的優(yōu)化起著重要的指導(dǎo)作用，因而受到了廣泛地關(guān)注和研究。特別地，當(dāng)設(shè)備的未來(lái)狀態(tài)能夠準(zhǔn)確預(yù)測(cè)后，對(duì)確定其維修檢測(cè)間隔期和開(kāi)展以可靠性為中心的維修(RCM)將大為方便［1］。而且現(xiàn)代維修策略的制定越來(lái)越依賴(lài)于設(shè)備當(dāng)前狀態(tài)而不是按照維修計(jì)劃按部就班，這就需要準(zhǔn)確預(yù)測(cè)出設(shè)備未來(lái)退化的發(fā)展趨勢(shì)。因此，對(duì)設(shè)備進(jìn)行早期故障識(shí)別并預(yù)測(cè)其剩余壽命(RUL)顯得尤為重要。

目前，很多學(xué)者通過(guò)建立復(fù)雜的系統(tǒng)退化模型來(lái)預(yù)測(cè)設(shè)備的剩余壽命，在這些模型的基礎(chǔ)上制定維修策略［2-5］。例如，Vlok等［6］將狀態(tài)監(jiān)測(cè)值作為協(xié)變量，用于比例風(fēng)險(xiǎn)模型(PHM)預(yù)測(cè);Kopnov，Pulkkinen等［7，8］假定完全監(jiān)測(cè)情況下，建立隨機(jī)退化過(guò)程模型來(lái)制定最優(yōu)維修策略。在實(shí)際應(yīng)用中，因?yàn)樵O(shè)備的退化狀態(tài)很難被直接監(jiān)測(cè)到，對(duì)設(shè)備的剩余壽命預(yù)測(cè)往往比較困難，而且監(jiān)測(cè)值也常常被噪聲所“污染”。因此，需要建立基于噪聲量測(cè)序列的系統(tǒng)退化模型，然后通過(guò)設(shè)備量測(cè)值對(duì)其內(nèi)在狀態(tài)值進(jìn)行估計(jì)，進(jìn)而計(jì)算設(shè)備剩余壽命的預(yù)測(cè)值。

貝葉斯理論非常適合解決這類(lèi)狀態(tài)估計(jì)問(wèn)題，它可以將過(guò)程監(jiān)測(cè)數(shù)據(jù)(一般是序貫觀測(cè)值)作為先驗(yàn)信息，對(duì)設(shè)備未來(lái)的狀態(tài)值進(jìn)行遞推估計(jì)。典型的貝葉斯估計(jì)方法是卡爾曼濾波(KF)，它可以求得高斯噪音下，線性狀態(tài)空間模型的最優(yōu)解［9］。然而在實(shí)際應(yīng)用中，絕大多數(shù)動(dòng)態(tài)系統(tǒng)往往是非線性的，量測(cè)噪聲也往往是非高斯的。為解決上述問(wèn)題，各種改進(jìn)算法和新算法相繼提出，例如，通過(guò)擴(kuò)展卡爾曼濾波(EKF)和無(wú)跡卡爾曼濾波(UKF)解決非線性高斯系統(tǒng)的狀態(tài)估計(jì)問(wèn)題［10］。近年來(lái)，隨著計(jì)算技術(shù)和硬件存儲(chǔ)技術(shù)的快速發(fā)展，序貫蒙特卡羅(SMC)方法，也叫粒子濾波(PF)算法，被越來(lái)越多地用于估計(jì)非線性非高斯條件下的系統(tǒng)狀態(tài)。文獻(xiàn)［11］將它應(yīng)用于馬爾科夫跳變線性系統(tǒng)，產(chǎn)生了很好的效果。

近年來(lái)，基于PF理論的故障診斷研究已經(jīng)得到了越來(lái)越多的關(guān)注，但是基于PF的剩余壽命預(yù)測(cè)問(wèn)題研究的還相對(duì)較少。為此本文結(jié)合粒子濾波對(duì)非線性非高斯系統(tǒng)的處理能力，提出了一套基于粒子濾波理論的設(shè)備剩余壽命預(yù)測(cè)框架。用非線性狀態(tài)空間模型(SSM)來(lái)表示系統(tǒng)的故障演化過(guò)程，通過(guò)粒子濾波算法預(yù)測(cè)出系統(tǒng)狀態(tài)的概率密度函數(shù)(PDF)。本文所提出的方法不僅能夠提供設(shè)備剩余壽命的預(yù)測(cè)結(jié)果，還能夠計(jì)算出結(jié)果的期望值和置信區(qū)間。另一方面，PF算法的研究仍處于發(fā)展階段，在實(shí)際系統(tǒng)壽命預(yù)測(cè)中的應(yīng)用還不成熟，絕大部分都是用于仿真領(lǐng)域，即人工時(shí)間序列模型產(chǎn)生輸入數(shù)據(jù)［12］。針對(duì)上述問(wèn)題，我們對(duì)PF算法在工程實(shí)際中的預(yù)測(cè)應(yīng)用開(kāi)展了實(shí)驗(yàn)研究。

1 理論基礎(chǔ)

1.1 問(wèn)題的提出

基于狀態(tài)空間模型的系統(tǒng)剩余壽命預(yù)測(cè)基本原理，是根據(jù)系統(tǒng)已知的先驗(yàn)信息得到的未知狀態(tài)的先驗(yàn)分布，結(jié)合觀測(cè)數(shù)據(jù)，來(lái)估計(jì)系統(tǒng)的后驗(yàn)分布及其剩余壽命，如圖1所示。

圖1 基于狀態(tài)空間模型的設(shè)備RUL預(yù)測(cè)Fig.1 Concept of dynamic model-based prognostics

一般地，用下面的狀態(tài)方程來(lái)描述系統(tǒng)模型及其狀態(tài)向量的演化過(guò)程:

其中，fk∶Rnx×Rnw→Rnx是非線性函數(shù)，時(shí)間步長(zhǎng)為tk=kΔt，{wk，k∈N}是已知分布形式的狀態(tài)噪聲向量序列。狀態(tài)向量序列{xk，k∈N}服從一階馬爾科夫過(guò)程。假設(shè)系統(tǒng)初始狀態(tài)p(x0)已知，狀態(tài)轉(zhuǎn)移概率p(xk|xk-1)由狀態(tài)方程(1)和分布類(lèi)型已知的噪聲向量wk定義。

系統(tǒng)的量測(cè)方程可由連續(xù)觀測(cè)值tk的觀測(cè)序列{zk，k∈N}來(lái)描述:

其中，hk∶Rnx×Rnv→Rnx是給定的非線性函數(shù)，{vk，k∈N}是已知分布形式的量測(cè)噪聲向量序列。并且，量測(cè)向量序列{zk，k∈N}狀態(tài)獨(dú)立于狀態(tài)過(guò)程{xk，k∈N}。也就是說(shuō)已知狀態(tài)xk時(shí)刻t，量測(cè)值z(mì)k的概率值不依賴(lài)于前一時(shí)刻的量測(cè)序列z0∶k-1=(z0，…，zk-1)。p(zk|xk，z0∶k-1)=p(zk|xk)可由量測(cè)方程(2)和噪聲向量vk確定。

1.2 傳統(tǒng)的比例風(fēng)險(xiǎn)剩余壽命預(yù)測(cè)模型

比例風(fēng)險(xiǎn)模型(PHM)作為一種剩余壽命預(yù)測(cè)方法，首先由Cox在1972年提出，并很快成為一種統(tǒng)計(jì)數(shù)據(jù)分析工具［13］。該模型綜合考慮設(shè)備運(yùn)行的壽命信息和各種設(shè)備狀態(tài)信息，從而有效地將狀態(tài)信息用于設(shè)備可靠性分析及剩余壽命預(yù)測(cè)，其基本形式如下:

其中，λ0(t)和g(X)都可能含未知參數(shù)，λ0(t)可以理解為g(X)=1下的標(biāo)準(zhǔn)故障率函數(shù)。由于威布爾分布可作為許多類(lèi)型設(shè)備(如真空管，滾珠軸承和電器的絕緣材料等)的壽命分布模型，因此下面主要研究Weibull-PHM剩余壽命預(yù)測(cè)模型。

Weibull-PHM的故障率函數(shù)定義如下［13］:

相應(yīng)的可靠度函數(shù)為

則平均剩余壽命為

其中，X為設(shè)備狀態(tài)向量，X=(x1，x2，…，xp)';λ(t|X)為設(shè)備在狀態(tài)X的故障率;α、δ為威布爾分布中的尺度參數(shù)和形狀參數(shù);X是p維狀態(tài)向量，反映設(shè)備的狀態(tài)信息;β為與X相對(duì)應(yīng)的變量系數(shù)，由上述公式可以看出，只要求得參數(shù)α、δ和β的值，即可確定Weibull-PHM的具體形式，模型參數(shù)用極大似然估計(jì)求解即可。

1.3 貝葉斯遞歸和粒子濾波

根據(jù)1.1節(jié)，我們想要得到的結(jié)果是后驗(yàn)概率分布p(xk|z0，k)，在貝葉斯框架下，已知k-1時(shí)刻的概率分布p(xk-1|z0，k-1)，根據(jù)系統(tǒng)狀態(tài)方程(1)來(lái)預(yù)測(cè)系統(tǒng)狀態(tài)xk先驗(yàn)概率分布，由Chapman-Kolmogorov方程得:

k時(shí)刻收集到新的觀測(cè)值z(mì)k，根據(jù)貝葉斯準(zhǔn)則更新系統(tǒng)狀態(tài)的后驗(yàn)分布，從而得到當(dāng)前系統(tǒng)狀態(tài)xk的后驗(yàn)概率分布:

其中的常數(shù)項(xiàng):

遞推公式(7)和公式(8)構(gòu)成了貝葉斯遞歸求解的基礎(chǔ)。然而，除了極少數(shù)情況，包括線性高斯?fàn)顟B(tài)空間模型(卡爾曼濾波)和有限狀態(tài)空間隱馬爾科夫鏈(Wohnam濾波)，用解析方法很難求解上述分布，因?yàn)橛?jì)算需要大量運(yùn)算和高維積分運(yùn)算。

這就需要其它有效求解方法——蒙特卡羅采樣。接下來(lái)，本文只分析蒙特卡羅方法的基本步驟，詳細(xì)推導(dǎo)參考文獻(xiàn)［14-16］。

一般來(lái)說(shuō)，已知量測(cè)向量z0;k，系統(tǒng)整個(gè)狀態(tài)序列x0;k的后驗(yàn)分布可以寫(xiě)成下述形式:

這里，δ(·)是狄拉克函數(shù)。假設(shè)系統(tǒng)真實(shí)后驗(yàn)概率p(x0∶k|z0∶k)已知，且能被采樣，則公式(10)可以由下式估計(jì):

其中，，i=1，2，…，Ns是從p(x0∶k|z0∶k)采樣得到的獨(dú)立隨機(jī)樣本集。

實(shí)際上，由于p(x0∶k|z0∶k)可能是多變量、非標(biāo)準(zhǔn)的，通常很難寫(xiě)成解析分布函數(shù)的組合形式，抽樣過(guò)程比較困難。為了克服這個(gè)問(wèn)題，可以借助重要性函數(shù)抽樣。所謂重要性函數(shù)就是指概率分布與p(x0∶k|z0∶k)相同，概率密度分布 π(x0∶k|z0∶k)已知，而且容易從中抽樣的分布函數(shù)，則公式(10)可以變換為:

其估計(jì)式為:

其中，

是系統(tǒng)狀態(tài)序列，i=1，2，…，Ns的權(quán)重，它可以通過(guò) π(|z0∶k)采樣得到;p(z0∶k|)是觀測(cè)序列的似然值。實(shí)際計(jì)算中，該權(quán)重很難求解，它需要知道p(z0∶k)= ∫p(z0∶k|x0∶k)p(x0∶k)dx0∶k，但上式很難表示為閉合形式進(jìn)行解析求解。為了解決這個(gè)問(wèn)題，后驗(yàn)分布p(x0∶k|z0∶k)可以由下式計(jì)算得到:

其中，

這樣，根據(jù)系統(tǒng)前k-1時(shí)刻狀態(tài)的分布p(x0∶k-1|z0∶k-1)來(lái)估計(jì)k時(shí)刻狀態(tài)的分布p(x0∶k|z0∶k)。從而，再根據(jù)貝葉斯濾波理論，由公式(8)遞推得到p(x0∶k|z0∶k):

其中，推導(dǎo)過(guò)程中再次用到了系統(tǒng)狀態(tài)方程(1)符合一階馬爾科夫過(guò)程這個(gè)假設(shè)，在系統(tǒng)狀態(tài)序列已知的情況下，量測(cè)方程(2)是條件獨(dú)立的。則重要性函數(shù)可以選取如下:

然后，可以得到非歸一化的和權(quán)值遞推公式如下:

重要性函數(shù)的選取是一個(gè)非常關(guān)鍵的問(wèn)題，選取原則之一就是使得重要性權(quán)值的方差最小。文獻(xiàn)［16］指出，選 擇 π (xk|x0∶k-1，z0∶k)=p(xk|x0∶k-1，z0∶k)可以滿足重要性權(quán)值的方程最小原則，但是這種選擇方法在實(shí)際中往往比較難以實(shí)現(xiàn)。從實(shí)際應(yīng)用角度來(lái)講，多數(shù)文獻(xiàn)普遍選擇采用 π(xk|x0∶k-1，z0∶k)=p(xk|xk-1)，這種方法盡管不是最優(yōu)方法，但是較容易實(shí)現(xiàn)。

然而，重采樣算法仍然存在嚴(yán)重的局限:在迭代很少幾步之后，重要性權(quán)值有可能集中到少數(shù)粒子上。這使得大量的更新運(yùn)算對(duì)最后的估計(jì)幾乎不起作用。通常，在迭代幾步之后，除了一個(gè)粒子權(quán)重外，其余粒子的權(quán)重都幾乎趨向于零。

為了避免這種退化，可以采用 bootstrap采樣技術(shù)［16］。對(duì)離散估計(jì)p(x0∶k|z0∶k)進(jìn)行重采樣Ns次，得到一系列粒子{xi*k，i=1，2，…，Ns}。本文采用了序貫重要性重采樣算法(SIRs)，具體過(guò)程如下:

歸一化權(quán)重

重采樣

如果重采樣，則根據(jù)權(quán)重選擇N個(gè)粒子序列

否則，

End

其重要性函數(shù)的權(quán)重為:

在每一步迭代過(guò)程都要進(jìn)行重采樣wik-1=1/Ns，經(jīng)過(guò)歸一化和更新后的權(quán)重等于量測(cè)量的似然值，即:

2 基于粒子濾波的剩余壽命預(yù)測(cè)方法

機(jī)械設(shè)備的振動(dòng)信號(hào)可以作為描述其退化狀態(tài)的指標(biāo)，為了準(zhǔn)確預(yù)測(cè)設(shè)備的退化趨勢(shì)及其剩余壽命，本文提出了基于設(shè)備振動(dòng)信號(hào)和粒子濾波參數(shù)估計(jì)的預(yù)測(cè)算法，見(jiàn)圖2所示。具體包括以下四步:

步驟一:(數(shù)據(jù)采集)首先通過(guò)振動(dòng)傳感器采集機(jī)械設(shè)備的原始信號(hào)，包括其從正常狀態(tài)到發(fā)生故障的全壽命數(shù)據(jù)。

步驟二:(特征提取)提取反映設(shè)備退化趨勢(shì)的特征指標(biāo)，本文選用的是特定頻帶能量和機(jī)械設(shè)備的磨損量。

步驟三:(預(yù)測(cè)和濾波)通過(guò)逐步采樣，應(yīng)用粒子濾波算法，計(jì)算粒子的重要性權(quán)重，然后重采樣技術(shù)，實(shí)現(xiàn)設(shè)備的壽命預(yù)測(cè)。

步驟四:(效果評(píng)估)為了確定所提方法的有效性和準(zhǔn)確性，首先對(duì)預(yù)測(cè)結(jié)果是否落于其95%置信區(qū)間進(jìn)行分析，然后提出預(yù)測(cè)效果評(píng)估指標(biāo)和相應(yīng)計(jì)算方法，進(jìn)一步對(duì)預(yù)測(cè)的結(jié)果進(jìn)行效果評(píng)估。

為了對(duì)上述預(yù)測(cè)算法的效果進(jìn)行綜合定量的評(píng)估，本文選用均方根誤差(RMSE)、平均絕對(duì)誤差(MAE)和方差絕對(duì)誤差(VAE)作為絕對(duì)誤差指標(biāo);用平均相對(duì)誤差(MARE)和方差相對(duì)誤差(VRE)作為相對(duì)誤差指標(biāo)。它們的計(jì)算公式如下:

其中，xt表示狀態(tài)真實(shí)值表示狀態(tài)預(yù)測(cè)值，上述指標(biāo)值越小，表示預(yù)測(cè)的精度越高。

3 應(yīng)用案例分析:齒輪箱全壽命實(shí)驗(yàn)及其剩余壽命預(yù)測(cè)

3.1 實(shí)驗(yàn)設(shè)置

齒輪箱由于其結(jié)構(gòu)緊湊、傳動(dòng)精確等優(yōu)點(diǎn)，廣泛應(yīng)用于設(shè)備傳動(dòng)系統(tǒng)中，由于其工作環(huán)境惡劣等原因，齒輪箱容易出現(xiàn)故障，因此本文選擇齒輪箱作為實(shí)驗(yàn)和研究對(duì)象，對(duì)其進(jìn)行剩余壽命預(yù)測(cè)，討論本文所提方法的實(shí)際應(yīng)用問(wèn)題。

實(shí)驗(yàn)設(shè)備和實(shí)驗(yàn)所用齒輪箱的結(jié)構(gòu)如圖3所示，傳感器測(cè)點(diǎn)布置如圖4所示。動(dòng)力源為電磁調(diào)速電機(jī)，型號(hào)YCT180-4A，磁粉制動(dòng)器為齒輪箱提供載荷，型號(hào)FZ200.K/F型風(fēng)冷磁粉制動(dòng)器。實(shí)驗(yàn)工況，扭矩為額定負(fù)載的2～2.5倍，目的是為了減少實(shí)驗(yàn)時(shí)間。實(shí)驗(yàn)所用齒輪箱主要參數(shù)見(jiàn)表1。

圖2 粒子濾波預(yù)測(cè)算法流程圖Fig.2 Flowchart of particle filter algorithm for prognosis

圖3 實(shí)驗(yàn)臺(tái)和齒輪箱結(jié)構(gòu)示意圖Fig.3 Sketch map of the test-bed and gearbox structure

圖4 齒輪箱振動(dòng)傳感器安置位置Fig.4 Accelerometers mounting

表1 齒輪箱主要參數(shù)Tab.1 Prime parameters of the gearbox

實(shí)驗(yàn)齒輪箱為本實(shí)驗(yàn)的核心部分(圖4)，它是一種二級(jí)斜齒輪箱，額定傳輸功率為0.75 kW;在箱體的不同部位上共裝有4個(gè)振動(dòng)加速度傳感器，能夠通過(guò)數(shù)據(jù)采集系統(tǒng)同時(shí)采集4路不同部位的振動(dòng)信號(hào)。通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，本次被測(cè)齒輪箱工作約450小時(shí)后，主要故障形式是齒輪齒面的嚴(yán)重磨損，實(shí)驗(yàn)結(jié)果如圖5所示。

圖5 齒輪實(shí)驗(yàn)前后對(duì)Fig.5 Comparison of gear before and after experiment

3.2 實(shí)驗(yàn)數(shù)據(jù)采集和特征提取

實(shí)驗(yàn)過(guò)程中對(duì)齒輪箱振動(dòng)信號(hào)進(jìn)行采集，采樣頻率為20 k，采集長(zhǎng)度為40 k，每小時(shí)采樣1次，累計(jì)采樣448 h。

對(duì)采集到的數(shù)據(jù)，利用小波包頻帶能量提取方法［17］，獲取不同時(shí)刻振動(dòng)信號(hào)特定頻帶能量(3～4.5 kHz，通過(guò)觀測(cè)，該頻帶的振動(dòng)信號(hào)能量隨工作時(shí)間逐漸增大，趨勢(shì)較為明顯，反映齒輪箱健康狀態(tài)的指標(biāo)較為理想){Ei;i=1，2，…，448}作為觀測(cè)特征值，如圖6所示。同時(shí)對(duì)輪齒磨損量進(jìn)行不完全數(shù)據(jù)采集，即僅對(duì)部分運(yùn)行時(shí)刻的磨損狀態(tài)進(jìn)行檢測(cè)，所獲得的磨損數(shù)據(jù)如表2。

表2 不同檢測(cè)時(shí)間輪齒磨損量Tab.2 The wear of gear teeth at different inspections

圖6 齒輪箱振動(dòng)信號(hào)特定頻帶能量Fig.6 Special band energy of vibration

3.3 基于PHM的齒輪箱剩余壽命預(yù)測(cè)

為了克服模型在小樣本且存在異常數(shù)據(jù)條件下預(yù)測(cè)精度較差的缺陷，需要選取一定的方法識(shí)別狀態(tài)信息樣本中的異常數(shù)據(jù)，以提高PHM的預(yù)測(cè)精度。本文采用支持限量機(jī)(SVM)算法，對(duì)狀態(tài)信息樣本中的異常數(shù)據(jù)進(jìn)行識(shí)別及分離。

圖7 異常數(shù)據(jù)識(shí)別結(jié)果Fig.7 The results of abnormal data detection

利用SVM對(duì)樣本進(jìn)行處理及平滑后，進(jìn)一步利用PHM進(jìn)行齒輪箱剩余壽命預(yù)測(cè)。通過(guò)極大似然算法可獲得最有參數(shù)為:α =438.9，δ=4.216，β =2.567 ×10-5，則齒輪箱的故障率函數(shù)為

通過(guò)計(jì)算，可得齒輪箱在不同時(shí)刻預(yù)測(cè)的平均剩余壽命，表3為不同檢測(cè)時(shí)間的剩余壽命預(yù)測(cè)結(jié)果。

表3 齒輪箱剩余壽命Tab.3 The remaining useful life of the gearbox

3.4 基于PF的齒輪箱剩余壽命預(yù)測(cè)

基于PF的齒輪箱剩余壽命預(yù)測(cè)模型與傳統(tǒng)PHM不同，其訓(xùn)練樣本綜合采用特定頻帶能量E和磨損量wl。利用PF對(duì)齒輪箱輪齒磨損過(guò)程進(jìn)行建模。由于Gamma過(guò)程具有平穩(wěn)、獨(dú)立增量等退化建模所需的所有屬性，被認(rèn)為是描述產(chǎn)品退化過(guò)程的首選方法［18］，磨損量演化的狀態(tài)空間模型如下:

取448 h的磨損量0.3 mm為故障閾值WLf，ε～(0，σ)為觀測(cè)誤差。為便于計(jì)算，取a(0)的初始值為1，分別可得初始參數(shù)參數(shù)ξ(0)=2 105、c(0)=48 744和σ(0)=968.7。令Wu表示所獲得的第u個(gè)磨損量觀測(cè)值，相應(yīng)的監(jiān)測(cè)時(shí)間為tu，通過(guò)PF方法可計(jì)算模型參數(shù)，其收斂過(guò)程如圖8所示。

圖8 參數(shù)收斂過(guò)程Fig.8 Convergence of the parameters

由于參數(shù)求解過(guò)程中，濾波粒子通過(guò)先驗(yàn)分布隨機(jī)生成，導(dǎo)致每次參數(shù)求解結(jié)果具有一定的隨機(jī)性，因此進(jìn)行1 000次重復(fù)參數(shù)求解后可得各參數(shù)的分布圖及參數(shù)均值和標(biāo)準(zhǔn)差，具體結(jié)果見(jiàn)表4。

表4 參數(shù)評(píng)估結(jié)果Tab.4 The results of parameter estimation

獲得模型參數(shù)后，可對(duì)不同時(shí)刻的磨損情況進(jìn)行評(píng)估，利用狀態(tài)方程可實(shí)現(xiàn)對(duì)未來(lái)裂紋概率密度的預(yù)測(cè)，預(yù)測(cè)公式為:

其中，p(xk|)為xk在粒子條件下的先驗(yàn)概率論，為粒子所對(duì)應(yīng)的權(quán)值，粒子集{，由粒子濾波算法獲得。則通過(guò)式(29)，可實(shí)現(xiàn)對(duì)齒輪箱磨損量的預(yù)測(cè)，預(yù)測(cè)結(jié)果如圖9所示。

圖9 磨損量預(yù)測(cè)結(jié)果Fig.9 The prediction results of wear

通過(guò)計(jì)算，分別得到400 h時(shí)的齒輪箱剩余壽命累積分布函數(shù)、概率密度函數(shù)和剩余壽命(如圖10，其中實(shí)際剩余壽命為48 h，預(yù)測(cè)值為60.7 h)。

圖10 400 h時(shí)的剩余壽命累積分布及概率密度函數(shù)Fig.10 The CDF and PDF of RUL at 400 h

采用同樣的方法，得到不同時(shí)刻的振動(dòng)檢測(cè)序列后，確定相應(yīng)的模型參數(shù)，表5為不同時(shí)刻測(cè)模型參數(shù)。

表5 不同檢測(cè)時(shí)刻的模型參數(shù)Tab.5 Estimated parameters at different inspection time

獲得不同時(shí)刻模型參數(shù)后，可計(jì)算得到剩余壽命的概率密度函數(shù)和平均剩余壽命，如圖11所示。表6為不同時(shí)刻的剩余壽命預(yù)測(cè)結(jié)果及95%的置信區(qū)間。根據(jù)公式(23)～(27)對(duì)預(yù)測(cè)的結(jié)果進(jìn)行效果評(píng)估(見(jiàn)表7)。

圖11 不同檢測(cè)時(shí)刻剩余壽命概率密度函數(shù)和均值Fig.11 Probability distribution function and mean of RUL at different inspection time

表6 剩余壽命預(yù)測(cè)值及置信區(qū)間Tab.6 Mean value and 95%confidence interval of RUL

表7 剩余壽命預(yù)測(cè)效果評(píng)估結(jié)果Tab.7 Performance evaluation results of RUL prediction

3.5 結(jié)果分析

表3和表6分別給出了兩種剩余壽命預(yù)測(cè)方法的預(yù)測(cè)結(jié)果。從表6可知，齒輪箱實(shí)際剩余壽命均落于置信度為0.95的置信區(qū)間之內(nèi)，因此基于PF的剩余壽命預(yù)測(cè)結(jié)果具有一定的合理性。表7給出了兩種方法RUL預(yù)測(cè)效果評(píng)估結(jié)果，從上述5個(gè)指標(biāo)對(duì)比可以看出，基于PF的剩余壽命預(yù)測(cè)無(wú)論是在預(yù)測(cè)的準(zhǔn)確度還是精確度方面均較PHM有明顯優(yōu)勢(shì)。這主要是由于基于PF的剩余壽命預(yù)測(cè)時(shí)綜合考慮了磨損量和振動(dòng)特征值，使模型能夠更為準(zhǔn)確的反映齒輪箱的磨損過(guò)程;而PHM單純采用振動(dòng)特征進(jìn)行剩余壽命預(yù)測(cè)，由于振動(dòng)特征并不能完全準(zhǔn)確地反映齒輪磨損狀態(tài)，因此預(yù)測(cè)結(jié)果較差。但PHM較基于PF的剩余壽命預(yù)測(cè)來(lái)說(shuō)，模型結(jié)構(gòu)比較簡(jiǎn)單，計(jì)算復(fù)雜度較低，模型求解過(guò)程更為容易。

4 結(jié)論

本文針對(duì)非線性系統(tǒng)剩余壽命預(yù)測(cè)問(wèn)題，提出了基于粒子濾波算法的剩余壽命預(yù)測(cè)框架，并根據(jù)系統(tǒng)的狀態(tài)監(jiān)測(cè)數(shù)據(jù)，結(jié)合粒子濾波算法，提出了剩余壽命預(yù)測(cè)方法，針對(duì)預(yù)測(cè)效果評(píng)估問(wèn)題，進(jìn)行了準(zhǔn)確性和精確性驗(yàn)證。最后通過(guò)齒輪箱全壽命實(shí)驗(yàn)進(jìn)行了實(shí)例分析，結(jié)果表明，與傳統(tǒng)的PHM預(yù)測(cè)方法進(jìn)行對(duì)比分析，基于粒子濾波模型的剩余壽命預(yù)測(cè)方法更加有效和準(zhǔn)確。

本文的研究結(jié)論為復(fù)雜噪聲條件下非高斯非線性系統(tǒng)的剩余壽命預(yù)測(cè)提供了一種可行的解決方法，具有一定的普適性。為增強(qiáng)該方法實(shí)際應(yīng)用的實(shí)時(shí)性，結(jié)合Rao-blackwellized進(jìn)行算法簡(jiǎn)化研究和實(shí)時(shí)預(yù)測(cè)研究;以及設(shè)備多故障模式條件下的故障追蹤和預(yù)測(cè)問(wèn)題，將做為進(jìn)一步研究的重點(diǎn)方向。

［1］Marseguerra M， ZioE， PodofilliniL. Condition-based maintenance optimization by means of genetic algorithms and Monte Carlo simulation［J］.Reliability Engineering and System Safety，2002，77:151-65.

［2］ Samanta P K，Vesely W E，Hsu F，et al.Degradation modeling with application to aging and maintenance effectivenessevaluations［R］. NUREG/CR-5612. US Nuclear Regulatory Commission，1991.

［3］Kopnov V A.Optimal degradation process control by twolevel policies［J］.Reliability Engineering and System Safety，1999，66:1-11.

［4］Lam C， Yeh R. Optimal maintenance policies for deteriorating systems under various maintenance strategies［J］.IEEE Trans Reliability，1994，43:423-430.

［5］Hontelez J A M，Burger H H，Wijnmalen D J D.Optimum condition-based maintenance policies for deteriorating systems with partial information［J］.Reliability Engineering and System Safety，1996，51:267-274.

［6］ Vlok P J，Coetzee J L，Banjevic D，et al.An application of vibration monitoring in proportional hazards models for optimal component replacement decisions［J］.Journal of Operational Research Sociaty，2002，53:193-202.

［7］Kopnov V A.Optimal degradation process control by twolevel policies［J］.Reliability Engineering and System Safety，1999，66:1-11.

［8］Pulkkinen U，Uryasev S.Optimal operational strategies for an inspected component［C］.In:Petersen K E，Rasmussen B，editors.Safety and Reliability'92， Proceedings of the European Safety and Reliability Conference'92.London，1992:896-907.

［9］Anderson B D，Moore J B.Optimal filtering［M］.Englewood Cliffs(NJ):Prentice Hall，1979.

［10］ Myotyri E，Pulkkinen U，Simola K.Application of stochastic filtering for lifetime prediction［J］.Reliability Engineering and System Safety，2006，91:200-208.

［11］ Zio E.Parameter identification in degradation modeling by reversible-jump Markov Chain Monte Carlo［J］. IEEE Tranctions on reliability，2009，58(1):123-131.

［12］ Guo D，Wang X，Chen R.New sequential Monte Carlo methods for nonlinear dynamic systems［J］.Statistics and Computing，2005，15:135-147.

［13］ Newby M.Perspective on Weibull proportional hazards models［J］.IEEE Transactions on Reliability，1994，43(2):217-223.

［14］ Orchard M E，Tang L，Goebel K，et al.A novel RSPF approach to prediction of high-risk，low-probability failure event［C］.Annual Conference of the Prognostics and Health Management Society，2009:1-8.

［15］ Cadini F，Zio E.A Monte Carlo method for the model-based estimation of nuclear reactor dynamics［J］.Ann.Nuclear Energy 2007，34(10):773-781.

［16］Doucet A，de Freitas N，Gordon N.Sequential Monte Carlo in practice［M］.New York:Springer-Verlag，2000.

［17］ Hu Q，He Z J.Fault diagnosis of rotating machinery based on improved wavelet package transform and SVMs ensemble［J］.Mechanical Systemsand SignalProcessing， 2007，21:688-705.

［18］ Van Noortwijk J M.A survey of the application of gamma processes in maintenance［J］.Reliability Engineering and System Safety，2009，94:2-21.

［19］ Boser B E，Guyon I M，Vapnik V N.A training algorithm for optimal margin classifiers［C］.The 5th AnnualACM Workshop on Computation Learning Theory，ACM Press，1992:144-152.