粒子濾波參數估計方法在齒輪箱剩余壽命預測中的應用研究

孫 磊，賈云獻，蔡麗影，張星輝

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

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

目前，很多學者通過建立復雜的系統退化模型來預測設備的剩余壽命，在這些模型的基礎上制定維修策略［2-5］。例如，Vlok等［6］將狀態監測值作為協變量，用于比例風險模型(PHM)預測;Kopnov，Pulkkinen等［7，8］假定完全監測情況下，建立隨機退化過程模型來制定最優維修策略。在實際應用中，因為設備的退化狀態很難被直接監測到，對設備的剩余壽命預測往往比較困難，而且監測值也常常被噪聲所“污染”。因此，需要建立基于噪聲量測序列的系統退化模型，然后通過設備量測值對其內在狀態值進行估計，進而計算設備剩余壽命的預測值。

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

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

1 理論基礎

1.1 問題的提出

基于狀態空間模型的系統剩余壽命預測基本原理，是根據系統已知的先驗信息得到的未知狀態的先驗分布，結合觀測數據，來估計系統的后驗分布及其剩余壽命，如圖1所示。

圖1 基于狀態空間模型的設備RUL預測Fig.1 Concept of dynamic model-based prognostics

一般地，用下面的狀態方程來描述系統模型及其狀態向量的演化過程:

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

系統的量測方程可由連續觀測值tk的觀測序列{zk，k∈N}來描述:

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

1.2 傳統的比例風險剩余壽命預測模型

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

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

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

相應的可靠度函數為

則平均剩余壽命為

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

1.3 貝葉斯遞歸和粒子濾波

根據1.1節，我們想要得到的結果是后驗概率分布p(xk|z0，k)，在貝葉斯框架下，已知k-1時刻的概率分布p(xk-1|z0，k-1)，根據系統狀態方程(1)來預測系統狀態xk先驗概率分布，由Chapman-Kolmogorov方程得:

k時刻收集到新的觀測值zk，根據貝葉斯準則更新系統狀態的后驗分布，從而得到當前系統狀態xk的后驗概率分布:

其中的常數項:

遞推公式(7)和公式(8)構成了貝葉斯遞歸求解的基礎。然而，除了極少數情況，包括線性高斯狀態空間模型(卡爾曼濾波)和有限狀態空間隱馬爾科夫鏈(Wohnam濾波)，用解析方法很難求解上述分布，因為計算需要大量運算和高維積分運算。

這就需要其它有效求解方法——蒙特卡羅采樣。接下來，本文只分析蒙特卡羅方法的基本步驟，詳細推導參考文獻［14-16］。

一般來說，已知量測向量z0;k，系統整個狀態序列x0;k的后驗分布可以寫成下述形式:

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

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

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

其估計式為:

其中，

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

其中，

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

其中，推導過程中再次用到了系統狀態方程(1)符合一階馬爾科夫過程這個假設，在系統狀態序列已知的情況下，量測方程(2)是條件獨立的。則重要性函數可以選取如下:

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

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

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

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

歸一化權重

重采樣

如果重采樣，則根據權重選擇N個粒子序列

否則，

End

其重要性函數的權重為:

在每一步迭代過程都要進行重采樣wik-1=1/Ns，經過歸一化和更新后的權重等于量測量的似然值，即:

2 基于粒子濾波的剩余壽命預測方法

機械設備的振動信號可以作為描述其退化狀態的指標，為了準確預測設備的退化趨勢及其剩余壽命，本文提出了基于設備振動信號和粒子濾波參數估計的預測算法，見圖2所示。具體包括以下四步:

步驟一:(數據采集)首先通過振動傳感器采集機械設備的原始信號，包括其從正常狀態到發生故障的全壽命數據。

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

步驟三:(預測和濾波)通過逐步采樣，應用粒子濾波算法，計算粒子的重要性權重，然后重采樣技術，實現設備的壽命預測。

步驟四:(效果評估)為了確定所提方法的有效性和準確性，首先對預測結果是否落于其95%置信區間進行分析，然后提出預測效果評估指標和相應計算方法，進一步對預測的結果進行效果評估。

為了對上述預測算法的效果進行綜合定量的評估，本文選用均方根誤差(RMSE)、平均絕對誤差(MAE)和方差絕對誤差(VAE)作為絕對誤差指標;用平均相對誤差(MARE)和方差相對誤差(VRE)作為相對誤差指標。它們的計算公式如下:

其中，xt表示狀態真實值表示狀態預測值，上述指標值越小，表示預測的精度越高。

3 應用案例分析:齒輪箱全壽命實驗及其剩余壽命預測

3.1 實驗設置

齒輪箱由于其結構緊湊、傳動精確等優點，廣泛應用于設備傳動系統中，由于其工作環境惡劣等原因，齒輪箱容易出現故障，因此本文選擇齒輪箱作為實驗和研究對象，對其進行剩余壽命預測，討論本文所提方法的實際應用問題。

實驗設備和實驗所用齒輪箱的結構如圖3所示，傳感器測點布置如圖4所示。動力源為電磁調速電機，型號YCT180-4A，磁粉制動器為齒輪箱提供載荷，型號FZ200.K/F型風冷磁粉制動器。實驗工況，扭矩為額定負載的2～2.5倍，目的是為了減少實驗時間。實驗所用齒輪箱主要參數見表1。

圖2 粒子濾波預測算法流程圖Fig.2 Flowchart of particle filter algorithm for prognosis

圖3 實驗臺和齒輪箱結構示意圖Fig.3 Sketch map of the test-bed and gearbox structure

圖4 齒輪箱振動傳感器安置位置Fig.4 Accelerometers mounting

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

實驗齒輪箱為本實驗的核心部分(圖4)，它是一種二級斜齒輪箱，額定傳輸功率為0.75 kW;在箱體的不同部位上共裝有4個振動加速度傳感器，能夠通過數據采集系統同時采集4路不同部位的振動信號。通過實驗發現，本次被測齒輪箱工作約450小時后，主要故障形式是齒輪齒面的嚴重磨損，實驗結果如圖5所示。

圖5 齒輪實驗前后對Fig.5 Comparison of gear before and after experiment

3.2 實驗數據采集和特征提取

實驗過程中對齒輪箱振動信號進行采集，采樣頻率為20 k，采集長度為40 k，每小時采樣1次，累計采樣448 h。

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

表2 不同檢測時間輪齒磨損量Tab.2 The wear of gear teeth at different inspections

圖6 齒輪箱振動信號特定頻帶能量Fig.6 Special band energy of vibration

3.3 基于PHM的齒輪箱剩余壽命預測

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

圖7 異常數據識別結果Fig.7 The results of abnormal data detection

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

通過計算，可得齒輪箱在不同時刻預測的平均剩余壽命，表3為不同檢測時間的剩余壽命預測結果。

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

3.4 基于PF的齒輪箱剩余壽命預測

基于PF的齒輪箱剩余壽命預測模型與傳統PHM不同，其訓練樣本綜合采用特定頻帶能量E和磨損量wl。利用PF對齒輪箱輪齒磨損過程進行建模。由于Gamma過程具有平穩、獨立增量等退化建模所需的所有屬性，被認為是描述產品退化過程的首選方法［18］，磨損量演化的狀態空間模型如下:

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

圖8 參數收斂過程Fig.8 Convergence of the parameters

由于參數求解過程中，濾波粒子通過先驗分布隨機生成，導致每次參數求解結果具有一定的隨機性，因此進行1 000次重復參數求解后可得各參數的分布圖及參數均值和標準差，具體結果見表4。

表4 參數評估結果Tab.4 The results of parameter estimation

獲得模型參數后，可對不同時刻的磨損情況進行評估，利用狀態方程可實現對未來裂紋概率密度的預測，預測公式為:

其中，p(xk|)為xk在粒子條件下的先驗概率論，為粒子所對應的權值，粒子集{，由粒子濾波算法獲得。則通過式(29)，可實現對齒輪箱磨損量的預測，預測結果如圖9所示。

圖9 磨損量預測結果Fig.9 The prediction results of wear

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

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

采用同樣的方法，得到不同時刻的振動檢測序列后，確定相應的模型參數，表5為不同時刻測模型參數。

表5 不同檢測時刻的模型參數Tab.5 Estimated parameters at different inspection time

獲得不同時刻模型參數后，可計算得到剩余壽命的概率密度函數和平均剩余壽命，如圖11所示。表6為不同時刻的剩余壽命預測結果及95%的置信區間。根據公式(23)～(27)對預測的結果進行效果評估(見表7)。

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

表6 剩余壽命預測值及置信區間Tab.6 Mean value and 95%confidence interval of RUL

表7 剩余壽命預測效果評估結果Tab.7 Performance evaluation results of RUL prediction

3.5 結果分析

表3和表6分別給出了兩種剩余壽命預測方法的預測結果。從表6可知，齒輪箱實際剩余壽命均落于置信度為0.95的置信區間之內，因此基于PF的剩余壽命預測結果具有一定的合理性。表7給出了兩種方法RUL預測效果評估結果，從上述5個指標對比可以看出，基于PF的剩余壽命預測無論是在預測的準確度還是精確度方面均較PHM有明顯優勢。這主要是由于基于PF的剩余壽命預測時綜合考慮了磨損量和振動特征值，使模型能夠更為準確的反映齒輪箱的磨損過程;而PHM單純采用振動特征進行剩余壽命預測，由于振動特征并不能完全準確地反映齒輪磨損狀態，因此預測結果較差。但PHM較基于PF的剩余壽命預測來說，模型結構比較簡單，計算復雜度較低，模型求解過程更為容易。

4 結論

本文針對非線性系統剩余壽命預測問題，提出了基于粒子濾波算法的剩余壽命預測框架，并根據系統的狀態監測數據，結合粒子濾波算法，提出了剩余壽命預測方法，針對預測效果評估問題，進行了準確性和精確性驗證。最后通過齒輪箱全壽命實驗進行了實例分析，結果表明，與傳統的PHM預測方法進行對比分析，基于粒子濾波模型的剩余壽命預測方法更加有效和準確。

本文的研究結論為復雜噪聲條件下非高斯非線性系統的剩余壽命預測提供了一種可行的解決方法，具有一定的普適性。為增強該方法實際應用的實時性，結合Rao-blackwellized進行算法簡化研究和實時預測研究;以及設備多故障模式條件下的故障追蹤和預測問題，將做為進一步研究的重點方向。

［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.