范 彬，胡 雷，胡蔦慶

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

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

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

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

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

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

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

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

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

假設(shè)動態(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)噪聲序列。

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

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

第二步，更新。即:

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

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

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

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

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

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

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

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

6)重采樣。根據(jù)粒子權(quán)值大小進行重采樣，得到新的粒子樣本集

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

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

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

1.3 存在的問題

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

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

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

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

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

2 基于DRT－PF的預測框架

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

2.1 通用預測模型

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

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

2.2 原始特征的預處理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實驗驗證與討論

3.1 預測結(jié)果分析

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

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

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

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

3.2 討論

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

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

4 結(jié)論

本文提出了一種基于DRT-PF的預測框架。該方法對原始預測特征進行平滑和單調(diào)化處理，再利用歷史數(shù)據(jù)中的退化速率信息輔助粒子濾波的迭代更新與外推預測。由于采用了通用的退化預測模型，因此該方法理論上對于任何緩變型退化系統(tǒng)的預測都適用。當擁有與系統(tǒng)退化趨勢高度相關(guān)的預測特征和足夠的歷史數(shù)據(jù)時，該方法能夠有效預測目標系統(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］.國防科學技術(shù)大學學報，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］鄧士杰，張英波，康海英，等.隨機濾波模型在變速箱剩余壽命預測中的應(yīng)用研究［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.