范 彬，胡 雷，胡蔦慶

(國防科技大學裝備綜合保障技術重點實驗室，湖南長沙410073)

科學技術的加速發(fā)展使得現(xiàn)代復雜工程系統(tǒng)的結構越來越復雜，同時也使得保障其長期安全運行的難度越來越高，一旦發(fā)生事故都將造成巨大的損失和危害。為避免事故發(fā)生，就需要對設備剩余使用壽命進行預測，從而針對性地制定維護方案，提高系統(tǒng)的可靠性和安全性，降低維修成本。因此，在許多重要的領域，對設備的故障和剩余使用壽命進行預測已經受到越來越高的重視［1］。

近年來，國內外學者開發(fā)了很多可用于設備剩余使用壽命預測的方法［2－10］，如神經網絡［5］、自回歸滑動平均(Auto-Regressive and Moving Average，ARMA)時間序列模型［6］、隱半馬爾可夫模型［7］、相關向量機［8］等基于數(shù)據驅動的方法，以及直接對具體零部件解析建模或仿真進行預測的基于物理模型方法［9，10］等。

其中，粒子濾波方法是一種基于數(shù)據與模型混合的預測方法，相比于傳統(tǒng)的數(shù)據驅動預測方法，其預測結果中的不確定性信息對于決策者十分重要。粒子濾波基于大數(shù)定理，利用大量的“粒子”樣本來近似狀態(tài)變量的真實后驗概率密度函數(shù)。盡管算法中的概率分布只是真實分布的一種近似，但因其非參數(shù)化的特點，得以擺脫解決非線性濾波問題時隨機量必須滿足高斯分布的制約，能表達比高斯模型更廣泛的分布，也對變量參數(shù)的非線性特性有更強的建模能力。粒子濾波已在目標跟蹤、航空航天、機器人定位、故障檢測等領域得到成功應用，近年來也開始被學者應用于故障與壽命預測［11－13］。

盡管如此，粒子濾波目前仍不夠成熟，除了算法本身的重要性函數(shù)選擇、粒子退化和樣本枯竭等問題［14］，在其應用于壽命預測時，還存在預測性能過度依賴于預測模型、對模型參數(shù)的初始分布過于敏感等問題［14，15］。因為對于不同的目標系統(tǒng)，選擇的退化特征不同，需要針對性地選擇合適的預測模型。而預測模型的復雜程度和準確程度大致成正比關系，使粒子濾波方法在剩余使用壽命預測上的應用及數(shù)據觀測受到了明顯的限制。除此之外，觀測數(shù)據通常還受外界噪聲等眾多因素影響，增加了建模難度。

為解決緩變型的退化失效預測問題，提出一種基于退化速率跟蹤粒子濾波(Degradation Rate Tracking-based Particle Filter，DRT-PF)的剩余使用壽命預測框架，該預測框架采用的預測模型只與預測特征的退化速率相關。由歷史數(shù)據預測特征計算得到的退化速率信息，既可以用于模型的迭代更新，同時也可以輔助外推預測。由于預測模型不需要根據具體系統(tǒng)和特征進行調整，因此該框架對所有緩變型退化系統(tǒng)的預測都適用。

1 基本的粒子濾波預測方法

粒子濾波是一種基于序貫蒙特卡洛方法的非線性濾波方法，其基本思想是:首先依據系統(tǒng)狀態(tài)變量的經驗條件分布，在狀態(tài)空間產生一組隨機“粒子”，然后根據測量值不斷調整粒子的權重和位置，通過調整后的粒子信息修正最初的經驗條件分布。其實質是用粒子及其權重組成的離散隨機測度近似相關的概率分布，并且根據算法遞推更新離散隨機測度。當樣本容量很大時，這種蒙特卡洛描述就近似于狀態(tài)變量真實的后驗概率密度函數(shù)。

1.1 最優(yōu)貝葉斯濾波

假設動態(tài)系統(tǒng)可以由數(shù)學模型描述，系統(tǒng)狀態(tài)方程和測量方程分別為

式中，xk為k時刻的狀態(tài)值，yk為觀測值，f(·)為系統(tǒng)狀態(tài)xk－1的非線性函數(shù)，{nk－1，k∈N}是平穩(wěn)噪聲序列，h(·)為系統(tǒng)狀態(tài)xk的非線性函數(shù)，{ωk，k∈N}是平穩(wěn)噪聲序列。

貝葉斯估計的遞推過程分為預測和更新兩步。

第一步，預測。假設在k－1時刻，狀態(tài)的后驗概率分布是已知的，則對于一階馬爾科夫過程，由根據系統(tǒng)的狀態(tài)轉移概率，推導出狀態(tài)的先驗概率

第二步，更新。即:

1.2 基本粒子濾波預測方法的實現(xiàn)

利用粒子濾波對某系統(tǒng)狀態(tài)進行預測，所采用的策略是首先指定預測模型，然后根據已有的觀測值對模型參數(shù)進行粒子濾波，通過迭代更新，使參數(shù)盡可能逼近實際值，最后將得到的模型參數(shù)估計值代入預測模型進行外推，實現(xiàn)對系統(tǒng)狀態(tài)的預測。具體的實現(xiàn)步驟如下:

1)提取退化特征時間序列。根據歷史觀測數(shù)據，提取適用于預測的退化特征x。

2)獲取先驗模型知識。通過對歷史數(shù)據退化特征的回歸擬合或根據經驗確定預測模型以及模型參數(shù)的初始分布。

3)粒子樣本初始化。由先驗概率p(x0)隨機產生粒子群，令初始粒子權值

4)采樣。相當于預測過程，即用狀態(tài)方程f預測k+1時刻的狀態(tài)

5)粒子狀態(tài)更新。根據當前時刻的觀測值，計算粒子樣本與觀測值的似然概率密度函數(shù)，以此更新粒子權值并對其進行歸一化:

可得k時刻狀態(tài)x的最小均方估計為

6)重采樣。根據粒子權值大小進行重采樣，得到新的粒子樣本集

7)判斷是否有新的觀測值。在k+1時刻，如果有新的觀測值，則轉到第4步;否則，繼續(xù)第8步。

8)將迭代更新完成的模型帶入狀態(tài)方程，對狀態(tài)x進行外推預測。

9)預測結果統(tǒng)計。根據預測模型和更新得到的模型參數(shù)估計值，對所有粒子進行外推預測，至指定時刻或達到指定閾值。對N個粒子在指定時刻對應的預測值xi或達到指定閾值的時間進行統(tǒng)計，得到x在指定時刻的概率密度分布或x達到閾值時間的概率密度分布。

1.3 存在的問題

盡管基本的粒子濾波預測方法已經在某些領域得到了應用，但在其實現(xiàn)過程中，還存在一些不足之處:

1)退化特征的“質量”直接影響預測結果的準確程度。所謂的“質量”，指的是用于預測的退化特征反映目標系統(tǒng)退化趨勢的能力。多數(shù)情況下，通過直接測量或初步提取得到的退化特征會受噪聲等多種因素影響，只能反映系統(tǒng)的大致退化趨勢。而將這樣的退化特征用于粒子濾波預測，其觀測噪聲對模型的迭代更新會產生較大影響，使模型參數(shù)容易產生發(fā)散。

2)預測模型不易選擇。在基本的粒子濾波預測方法中，預測模型的選擇是必不可少的一個步驟，而且預測結果的準確程度很大程度取決于預測模型的準確與否。預測模型的細微改變，就可能導致預測結果的明顯差異［15］，這嚴重影響了預測方法的整體穩(wěn)定性。

3)預測模型參數(shù)較多，多維狀態(tài)變量的粒子濾波難度加大。預測模型的選擇，通常是根據對歷史數(shù)據的擬合誤差大小來判斷的。簡單的預測模型常常不能很準確地刻畫歷史數(shù)據的退化趨勢，而準確的預測模型又比較復雜，模型參數(shù)多，會加大方法的實現(xiàn)難度，同時還會使算法的運算量成倍增長。

4)預測結果對模型參數(shù)的初始分布設置過于敏感。由于通常的預測模型參數(shù)在迭代更新階段是沒有實際觀測值的，因此這些參數(shù)的初始分布設置直接決定了參數(shù)是否能夠收斂及其收斂速度。盡管在文獻［14］中，利用Dempster-Shafer理論成功給出了合適的模型參數(shù)初始分布，但實際上這種方法并不能保證始終得到合適的結果。

2 基于DRT－PF的預測框架

模型的選擇及其參數(shù)的設置在傳統(tǒng)的粒子濾波預測方法中占據了過于重要的地位，不僅占用了大量計算時間，而且在辨識過程中存在的較大不確定性也會影響預測結果的準確性。為了解決這個問題，從預測模型的角度，提出了一種DRT-PF預測方法。該方法為了提高預測方法的通用性，采用了一種根據退化速率來迭代更新退化狀態(tài)的通用模型。對于不同的目標系統(tǒng)，其差異性體現(xiàn)為退化速率序列的不同，而預測模型是一致的。

2.1 通用預測模型

采用通用的預測模型式(6)代替之前需要進行實現(xiàn)辨識的預測模型，簡化預測流程。

式中:yk為觀測值，對應于原始的退化特征;xk為狀態(tài)值，對應于預處理后的退化特征;ak為退化狀態(tài)的退化速率，即后一時刻狀態(tài)與前一時刻狀態(tài)的比值xk+1/xk，也是狀態(tài)方程中唯一的模型參數(shù);ωk為觀測誤差。通用的預測模型將使粒子狀態(tài)更加簡單，只包括狀態(tài)x和退化速率a，大大降低了計算量。另外，退化速率a在迭代更新階段可以計算得到實測值，所以也能夠對a進行重采樣，因此能夠明顯降低模型參數(shù)的初始分布設置對于預測結果的影響。

2.2 原始特征的預處理

在實際預測中，通常只關心退化特征反映的系統(tǒng)退化趨勢，而并不關心退化特征所包含的其他細節(jié)信息。而由于噪聲等外界因素的影響，對于原始退化特征，退化速率的取值范圍難以確定。所以要盡可能地去除噪聲等外界因素的影響，這樣得到的特征才能更直觀地反映目標系統(tǒng)的退化趨勢，同時有利于確定退化速率取值范圍。另外，設備的退化過程是單調不可逆的。因此對初步提取的退化特征依次進行平滑和單調化的預處理，能夠更適用于粒子濾波預測。圖1為預處理前后的退化特征。

圖1 預處理前后的退化特征對比示意圖Fig.1 Comparison between raw degradation features and the pretreated features

2.3 歷史數(shù)據的先驗信息獲取

在以上預測模型的基礎上，從歷史數(shù)據中可獲取兩方面的先驗信息。一是預測模型中各參數(shù)的初始分布信息，與粒子樣本的初始化相關;二是退化速率隨狀態(tài)x變化的分布情況，主要應用于外推預測階段。先驗信息的獲取流程如圖2所示。

圖2 先驗信息獲取流程圖Fig.2 The flowchart for obtaining prior information

1)模型參數(shù)的初始分布。在進行預測之前，需要根據模型參數(shù)的初始分布進行粒子初始化，所以首先應確定預測模型參數(shù)的初始分布信息。粒子初始化使用一致分布，根據式(6)，需要進行初始化的參數(shù)包括:退化狀態(tài)x，觀測噪聲標準差σ(ωk的標準差)，退化速率a。其中x與a在迭代過程中會進行重采樣，所以其初始分布只需根據歷史數(shù)據的初始觀測值給出略大的區(qū)間即可。而σ則可以根據各組歷史數(shù)據的原始信號與預處理后的信號殘差進行統(tǒng)計得到。

2)不同階段的退化速率分布圖。a是預測模型中唯一的參數(shù)，體現(xiàn)系統(tǒng)的退化速度，更直接決定了對系統(tǒng)剩余壽命的預測結果，因此需要在預測過程中準確地估計退化速率的演變規(guī)律。一方面，因為設備退化是對時間單調變化的，所以退化特征經過平滑和單調化預處理后，退化速率a的取值上邊界或下邊界很容易確定(根據退化特征不同，可得a≤1或a≥1)，同時可以根據歷史數(shù)據中a的極值情況確定其另一個取值邊界。另一方面，由歷史數(shù)據統(tǒng)計信息可以給出退化速率a隨系統(tǒng)狀態(tài)x變化的先驗分布p(a|x)，作為預測階段粒子狀態(tài)的指導信息。首先將歷史數(shù)據的退化速率對應到退化特征坐標軸上，再利用插值方法計算退化特征軸上的退化速率趨勢曲線。最后對所有趨勢曲線進行統(tǒng)計，在每個退化特征點上，對其相應的退化速率進行正態(tài)分布擬合，可以得到退化速率隨特征變化的分布信息。

2.4 監(jiān)測數(shù)據的預測

在獲取了所需的先驗信息之后，就可以對已有的監(jiān)測數(shù)據進行預測。首先根據模型參數(shù)的初始分布信息，對N個粒子樣本初始化，然后依次進行迭代估計和外推預測兩個步驟，最后對所有粒子樣本的狀態(tài)進行統(tǒng)計，得到最終的預測結果。預測流程圖如圖3所示。

1)迭代估計階段。首先對由已知的監(jiān)測數(shù)據提取的退化特征進行平滑和單調化的預處理。然后利用基于通用預測模型的粒子濾波方法對模型參數(shù)進行迭代估計。在迭代估計過程中，將對退化狀態(tài)值x和退化速率a分別進行重采樣。

2)外推預測階段。利用已知的觀測數(shù)據完成對模型參數(shù)的迭代估計之后，再結合預測模型以及之前得到的退化速率分布圖進行外推預測。設當前退化狀態(tài)預測值為

那么，根據p(a|x)，隨機生成當前狀態(tài)下的先驗退化速率aik(i=1，2，…，N)。然后再根據預測模型，可以得到下一步的退化狀態(tài)值。

依次迭代預測，當退化狀態(tài)預測值大于或小于失效閾值時，預測結束。對所有粒子的當前狀態(tài)進行統(tǒng)計，就能得到失效時間預測的概率密度分布以及預測結果的置信區(qū)間等信息。

圖3 監(jiān)測數(shù)據預測流程圖Fig.3 The flowchart of prediction scheme for monitoring data

3 實驗驗證與討論

3.1 預測結果分析

圖4所示為一個加速壽命實驗的滾動軸承的剩余使用壽命預測結果。在圖4(a)中，前500個數(shù)據點(數(shù)據點之間間隔10 s)被用來迭代更新預測模型，再用更新之后的模型結合退化速率分布圖，對軸承500個數(shù)據點之后的狀態(tài)進行預測。實際的失效時間位于預測結果的90%置信區(qū)間以內，表示預測結果有效。但因為用于迭代更新的數(shù)據點相對較少，所以預測結果的不確定度大。在圖4(b)中可以看到，隨著已知數(shù)據點的增加，預測結果逐漸逼近實際值，置信區(qū)間的寬度隨已知數(shù)據點的增加而變窄，表示預測結果的準確性在升高，而不確定性在逐漸減低。

圖4 軸承剩余使用壽命預測結果Fig.4 Prediction results of bearing remaining useful life

圖5 電池剩余使用壽命預測結果Fig.5 Prediction results of battery remaining useful life

圖5所示為NASA鋰離子電池的循環(huán)壽命預測結果。鋰離子電池循環(huán)壽命的退化特征為電池容量，呈逐步下降的趨勢，與軸承的退化特征剛好相反，但對本文提出的預測方法沒有影響。因為受到松弛效應的影響，鋰離子電池的原始退化特征存在明顯的非單調下降的部分，但采用文獻［16］中的方法消除松弛效應之后，即可按照文中提出的預測框架進行處理。從圖5(a)中可以看出，該方法能夠對鋰離子電池的循環(huán)壽命進行有效估計。而在圖5(b)中，可以看到預測結果收斂速度更快，準確性更高。在圖示范圍內，預測中值與實際值的最大誤差僅為11個循環(huán)。這是因為相對于圖4，鋰離子電池的預測特征與系統(tǒng)退化過程相關性更強，而且特征受其他因素影響更小，所以相應的預測效果也更好。

3.2 討論

如前所述，通過對預測特征的退化速率進行跟蹤，能夠使粒子濾波對目標系統(tǒng)的狀態(tài)估計迅速收斂，在預測階段也能比較準確地逼近真實的退化趨勢。預測結果的準確度隨已知數(shù)據長度逐漸增減，不確定度隨之下降。而且，對于無論是代表劣化程度(逐漸下降)還是代表健康程度(逐漸升高)的退化特征，該預測框架都同樣適用，表現(xiàn)出了其良好的通用性，并在一定程度上簡化了粒子濾波預測方法。

但是，對比兩種不同系統(tǒng)的預測結果，可以發(fā)現(xiàn)預測特征仍然是對算法影響很大的一個因素，對于機械設備的預測尤為重要。由振動信號提取的預測特征往往受外界噪聲等影響嚴重，如果預測特征和系統(tǒng)的退化趨勢相關程度較低，那么算法的預測效果也會相對下降，所以如何提取與系統(tǒng)高度相關的預測特征還值得深入研究。

4 結論

本文提出了一種基于DRT-PF的預測框架。該方法對原始預測特征進行平滑和單調化處理，再利用歷史數(shù)據中的退化速率信息輔助粒子濾波的迭代更新與外推預測。由于采用了通用的退化預測模型，因此該方法理論上對于任何緩變型退化系統(tǒng)的預測都適用。當擁有與系統(tǒng)退化趨勢高度相關的預測特征和足夠的歷史數(shù)據時，該方法能夠有效預測目標系統(tǒng)的剩余使用壽命。

References)

［1］Kruzic J J.Predicting fatigue failures［J］.Science，2009，325(5937):156－158.

［2］Bagul Y G，Zeid I，Kamarthi S V.Overview of remaining useful life methodologies［C］//Proceedings of the ASME 2008 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference，American Society of Mechanical Engineers，New York，2008:1391－1400.

［3］王小林，程志軍，郭波.基于維納過程金屬化膜電容器的剩余壽命預測［J］.國防科學技術大學學報，2011，33(4):146－151.WANG Xiaolin，CHENG Zhijun，GUO Bo.Residual life forecasting of metallized film capacitor based on wiener process［J］.Journal of National University of Defense Technology，2011，33(4):146－151.(in Chinese)

［4］鄧士杰，張英波，康海英，等.隨機濾波模型在變速箱剩余壽命預測中的應用研究［J］.軍械工程學院學報，2013(1):29－32.DENG Shijie，ZHANG Yingbo，KANG Haiying，et al.Research of stochastic filtering model for gear box residual useful life prediction［J］.Journal of Ordnance Engineering College，2013(1):29－32.(in Chinese)

［5］Shao Y，Nezu K.Prognosis of remaining bearing life using neural networks［J］.Proceedings of the Institution of Mechanical Engineers，Part I:Journal of Systems and Control Engineering，2000，214(3):217－230.

［6］Wang T Y，Yu J B，Lee J，et al.A similarity-based prognostics approach for remaining useful life estimation of engineered systems［C］//Proceedings of the International Conference on Prognostics and Health Management，IEEE，2008:1－6.

［7］Dong M，He D.A segmental hidden semi-Markov model-based diagnostics and prognostics framework and methodology［J］.Mechanical Systems and Signal Processing，2007，21(5):2248－2266.

［8］Di Maio F，Tsui K L，Zio E.Combining relevance vector machines and exponential regression for bearing residual life estimation［J］.Mechanical Systems and Signal Processing，2012，31:405－427.

［9］Li C J，Lee H.Gear fatigue crack prognosis using embedded model，gear dynamic model and fracture mechanics［J］.Mechanical Systems and Signal Processing，2005，19(4):836－849.

［10］Feng Z P，Zuo M J，Chu F L.Application of regularization dimension to gear damage assessment［J］.Mechanical Systems and Signal Processing，2010，24(4):1081－1098.

［11］Orchard M E.A particle filtering-based framework for on-line fault diagnosis and failure prognosis［D］.Atlanta:Georgia Institute of Technology，2007.

［12］Zio E，Peloni G.Particle filtering prognostic estimation of the remaining useful life of nonlinear components［J］.Reliability Engineering and System Safety，2011，96(3):403－409.

［13］Jin G，Matthews D E，Zhou Z B.A bayesian framework for on-line degradation assessment and residual life prediction of secondary batteries in spacecraft［J］.Reliability Engineering and System Safety，2013，113:7－20.

［14］He W，Williard N，Osterman M，et al.Prognostics of lithium-ion based on dempster-shafer theory and the bayesian monte carlo method［J］.Journal of Power Sources，2011，196(23):10134－10321.

［15］Wang D，Miao Q，Pecht M.Prognostics of lithium-ion batteries based on relevance vectors and a conditional threeparameter capacity degradation model［J］.Journal of Power Sources，2013，239(1):253－264.

［16］Tang S J，Yu C Q，Wang X，et al.Remaining useful life prediction of lithium-ion batteries based on the wiener process with measurement error［J］.Energies，2014，7(2):520－547.