范 彬,胡 雷,胡蔦慶
(國防科技大學(xué)裝備綜合保障技術(shù)重點實驗室,湖南長沙410073)
科學(xué)技術(shù)的加速發(fā)展使得現(xiàn)代復(fù)雜工程系統(tǒng)的結(jié)構(gòu)越來越復(fù)雜,同時也使得保障其長期安全運行的難度越來越高,一旦發(fā)生事故都將造成巨大的損失和危害。為避免事故發(fā)生,就需要對設(shè)備剩余使用壽命進(jìn)行預(yù)測,從而針對性地制定維護(hù)方案,提高系統(tǒng)的可靠性和安全性,降低維修成本。因此,在許多重要的領(lǐng)域,對設(shè)備的故障和剩余使用壽命進(jìn)行預(yù)測已經(jīng)受到越來越高的重視[1]。
近年來,國內(nèi)外學(xué)者開發(fā)了很多可用于設(shè)備剩余使用壽命預(yù)測的方法[2-10],如神經(jīng)網(wǎng)絡(luò)[5]、自回歸滑動平均(Auto-Regressive and Moving Average,ARMA)時間序列模型[6]、隱半馬爾可夫模型[7]、相關(guān)向量機(jī)[8]等基于數(shù)據(jù)驅(qū)動的方法,以及直接對具體零部件解析建?;蚍抡孢M(jìn)行預(yù)測的基于物理模型方法[9,10]等。
其中,粒子濾波方法是一種基于數(shù)據(jù)與模型混合的預(yù)測方法,相比于傳統(tǒng)的數(shù)據(jù)驅(qū)動預(yù)測方法,其預(yù)測結(jié)果中的不確定性信息對于決策者十分重要。粒子濾波基于大數(shù)定理,利用大量的“粒子”樣本來近似狀態(tài)變量的真實后驗概率密度函數(shù)。盡管算法中的概率分布只是真實分布的一種近似,但因其非參數(shù)化的特點,得以擺脫解決非線性濾波問題時隨機(jī)量必須滿足高斯分布的制約,能表達(dá)比高斯模型更廣泛的分布,也對變量參數(shù)的非線性特性有更強(qiáng)的建模能力。粒子濾波已在目標(biāo)跟蹤、航空航天、機(jī)器人定位、故障檢測等領(lǐng)域得到成功應(yīng)用,近年來也開始被學(xué)者應(yīng)用于故障與壽命預(yù)測[11-13]。
盡管如此,粒子濾波目前仍不夠成熟,除了算法本身的重要性函數(shù)選擇、粒子退化和樣本枯竭等問題[14],在其應(yīng)用于壽命預(yù)測時,還存在預(yù)測性能過度依賴于預(yù)測模型、對模型參數(shù)的初始分布過于敏感等問題[14,15]。因為對于不同的目標(biāo)系統(tǒng),選擇的退化特征不同,需要針對性地選擇合適的預(yù)測模型。而預(yù)測模型的復(fù)雜程度和準(zhǔn)確程度大致成正比關(guān)系,使粒子濾波方法在剩余使用壽命預(yù)測上的應(yīng)用及數(shù)據(jù)觀測受到了明顯的限制。除此之外,觀測數(shù)據(jù)通常還受外界噪聲等眾多因素影響,增加了建模難度。
為解決緩變型的退化失效預(yù)測問題,提出一種基于退化速率跟蹤粒子濾波(Degradation Rate Tracking-based Particle Filter,DRT-PF)的剩余使用壽命預(yù)測框架,該預(yù)測框架采用的預(yù)測模型只與預(yù)測特征的退化速率相關(guān)。由歷史數(shù)據(jù)預(yù)測特征計算得到的退化速率信息,既可以用于模型的迭代更新,同時也可以輔助外推預(yù)測。由于預(yù)測模型不需要根據(jù)具體系統(tǒng)和特征進(jìn)行調(diào)整,因此該框架對所有緩變型退化系統(tǒng)的預(yù)測都適用。
粒子濾波是一種基于序貫蒙特卡洛方法的非線性濾波方法,其基本思想是:首先依據(jù)系統(tǒng)狀態(tài)變量的經(jīng)驗條件分布,在狀態(tài)空間產(chǎn)生一組隨機(jī)“粒子”,然后根據(jù)測量值不斷調(diào)整粒子的權(quán)重和位置,通過調(diào)整后的粒子信息修正最初的經(jīng)驗條件分布。其實質(zhì)是用粒子及其權(quán)重組成的離散隨機(jī)測度近似相關(guān)的概率分布,并且根據(jù)算法遞推更新離散隨機(jī)測度。當(dāng)樣本容量很大時,這種蒙特卡洛描述就近似于狀態(tài)變量真實的后驗概率密度函數(shù)。
假設(shè)動態(tài)系統(tǒng)可以由數(shù)學(xué)模型描述,系統(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)噪聲序列。
貝葉斯估計的遞推過程分為預(yù)測和更新兩步。
第一步,預(yù)測。假設(shè)在k-1時刻,狀態(tài)的后驗概率分布是已知的,則對于一階馬爾科夫過程,由根據(jù)系統(tǒng)的狀態(tài)轉(zhuǎn)移概率,推導(dǎo)出狀態(tài)的先驗概率
第二步,更新。即:
利用粒子濾波對某系統(tǒng)狀態(tài)進(jìn)行預(yù)測,所采用的策略是首先指定預(yù)測模型,然后根據(jù)已有的觀測值對模型參數(shù)進(jìn)行粒子濾波,通過迭代更新,使參數(shù)盡可能逼近實際值,最后將得到的模型參數(shù)估計值代入預(yù)測模型進(jìn)行外推,實現(xiàn)對系統(tǒng)狀態(tài)的預(yù)測。具體的實現(xiàn)步驟如下:
1)提取退化特征時間序列。根據(jù)歷史觀測數(shù)據(jù),提取適用于預(yù)測的退化特征x。
2)獲取先驗?zāi)P椭R。通過對歷史數(shù)據(jù)退化特征的回歸擬合或根據(jù)經(jīng)驗確定預(yù)測模型以及模型參數(shù)的初始分布。
3)粒子樣本初始化。由先驗概率p(x0)隨機(jī)產(chǎn)生粒子群,令初始粒子權(quán)值
4)采樣。相當(dāng)于預(yù)測過程,即用狀態(tài)方程f預(yù)測k+1時刻的狀態(tài)
5)粒子狀態(tài)更新。根據(jù)當(dāng)前時刻的觀測值,計算粒子樣本與觀測值的似然概率密度函數(shù),以此更新粒子權(quán)值并對其進(jìn)行歸一化:
可得k時刻狀態(tài)x的最小均方估計為
6)重采樣。根據(jù)粒子權(quán)值大小進(jìn)行重采樣,得到新的粒子樣本集
7)判斷是否有新的觀測值。在k+1時刻,如果有新的觀測值,則轉(zhuǎn)到第4步;否則,繼續(xù)第8步。
8)將迭代更新完成的模型帶入狀態(tài)方程,對狀態(tài)x進(jìn)行外推預(yù)測。
9)預(yù)測結(jié)果統(tǒng)計。根據(jù)預(yù)測模型和更新得到的模型參數(shù)估計值,對所有粒子進(jìn)行外推預(yù)測,至指定時刻或達(dá)到指定閾值。對N個粒子在指定時刻對應(yīng)的預(yù)測值xi或達(dá)到指定閾值的時間進(jìn)行統(tǒng)計,得到x在指定時刻的概率密度分布或x達(dá)到閾值時間的概率密度分布。
盡管基本的粒子濾波預(yù)測方法已經(jīng)在某些領(lǐng)域得到了應(yīng)用,但在其實現(xiàn)過程中,還存在一些不足之處:
1)退化特征的“質(zhì)量”直接影響預(yù)測結(jié)果的準(zhǔn)確程度。所謂的“質(zhì)量”,指的是用于預(yù)測的退化特征反映目標(biāo)系統(tǒng)退化趨勢的能力。多數(shù)情況下,通過直接測量或初步提取得到的退化特征會受噪聲等多種因素影響,只能反映系統(tǒng)的大致退化趨勢。而將這樣的退化特征用于粒子濾波預(yù)測,其觀測噪聲對模型的迭代更新會產(chǎn)生較大影響,使模型參數(shù)容易產(chǎn)生發(fā)散。
2)預(yù)測模型不易選擇。在基本的粒子濾波預(yù)測方法中,預(yù)測模型的選擇是必不可少的一個步驟,而且預(yù)測結(jié)果的準(zhǔn)確程度很大程度取決于預(yù)測模型的準(zhǔn)確與否。預(yù)測模型的細(xì)微改變,就可能導(dǎo)致預(yù)測結(jié)果的明顯差異[15],這嚴(yán)重影響了預(yù)測方法的整體穩(wěn)定性。
3)預(yù)測模型參數(shù)較多,多維狀態(tài)變量的粒子濾波難度加大。預(yù)測模型的選擇,通常是根據(jù)對歷史數(shù)據(jù)的擬合誤差大小來判斷的。簡單的預(yù)測模型常常不能很準(zhǔn)確地刻畫歷史數(shù)據(jù)的退化趨勢,而準(zhǔn)確的預(yù)測模型又比較復(fù)雜,模型參數(shù)多,會加大方法的實現(xiàn)難度,同時還會使算法的運算量成倍增長。
4)預(yù)測結(jié)果對模型參數(shù)的初始分布設(shè)置過于敏感。由于通常的預(yù)測模型參數(shù)在迭代更新階段是沒有實際觀測值的,因此這些參數(shù)的初始分布設(shè)置直接決定了參數(shù)是否能夠收斂及其收斂速度。盡管在文獻(xiàn)[14]中,利用Dempster-Shafer理論成功給出了合適的模型參數(shù)初始分布,但實際上這種方法并不能保證始終得到合適的結(jié)果。
模型的選擇及其參數(shù)的設(shè)置在傳統(tǒng)的粒子濾波預(yù)測方法中占據(jù)了過于重要的地位,不僅占用了大量計算時間,而且在辨識過程中存在的較大不確定性也會影響預(yù)測結(jié)果的準(zhǔn)確性。為了解決這個問題,從預(yù)測模型的角度,提出了一種DRT-PF預(yù)測方法。該方法為了提高預(yù)測方法的通用性,采用了一種根據(jù)退化速率來迭代更新退化狀態(tài)的通用模型。對于不同的目標(biāo)系統(tǒng),其差異性體現(xiàn)為退化速率序列的不同,而預(yù)測模型是一致的。
采用通用的預(yù)測模型式(6)代替之前需要進(jìn)行實現(xiàn)辨識的預(yù)測模型,簡化預(yù)測流程。
式中:yk為觀測值,對應(yīng)于原始的退化特征;xk為狀態(tài)值,對應(yīng)于預(yù)處理后的退化特征;ak為退化狀態(tài)的退化速率,即后一時刻狀態(tài)與前一時刻狀態(tài)的比值xk+1/xk,也是狀態(tài)方程中唯一的模型參數(shù);ωk為觀測誤差。通用的預(yù)測模型將使粒子狀態(tài)更加簡單,只包括狀態(tài)x和退化速率a,大大降低了計算量。另外,退化速率a在迭代更新階段可以計算得到實測值,所以也能夠?qū)進(jìn)行重采樣,因此能夠明顯降低模型參數(shù)的初始分布設(shè)置對于預(yù)測結(jié)果的影響。
在實際預(yù)測中,通常只關(guān)心退化特征反映的系統(tǒng)退化趨勢,而并不關(guān)心退化特征所包含的其他細(xì)節(jié)信息。而由于噪聲等外界因素的影響,對于原始退化特征,退化速率的取值范圍難以確定。所以要盡可能地去除噪聲等外界因素的影響,這樣得到的特征才能更直觀地反映目標(biāo)系統(tǒng)的退化趨勢,同時有利于確定退化速率取值范圍。另外,設(shè)備的退化過程是單調(diào)不可逆的。因此對初步提取的退化特征依次進(jìn)行平滑和單調(diào)化的預(yù)處理,能夠更適用于粒子濾波預(yù)測。圖1為預(yù)處理前后的退化特征。
圖1 預(yù)處理前后的退化特征對比示意圖Fig.1 Comparison between raw degradation features and the pretreated features
在以上預(yù)測模型的基礎(chǔ)上,從歷史數(shù)據(jù)中可獲取兩方面的先驗信息。一是預(yù)測模型中各參數(shù)的初始分布信息,與粒子樣本的初始化相關(guān);二是退化速率隨狀態(tài)x變化的分布情況,主要應(yīng)用于外推預(yù)測階段。先驗信息的獲取流程如圖2所示。
圖2 先驗信息獲取流程圖Fig.2 The flowchart for obtaining prior information
1)模型參數(shù)的初始分布。在進(jìn)行預(yù)測之前,需要根據(jù)模型參數(shù)的初始分布進(jìn)行粒子初始化,所以首先應(yīng)確定預(yù)測模型參數(shù)的初始分布信息。粒子初始化使用一致分布,根據(jù)式(6),需要進(jìn)行初始化的參數(shù)包括:退化狀態(tài)x,觀測噪聲標(biāo)準(zhǔn)差σ(ωk的標(biāo)準(zhǔn)差),退化速率a。其中x與a在迭代過程中會進(jìn)行重采樣,所以其初始分布只需根據(jù)歷史數(shù)據(jù)的初始觀測值給出略大的區(qū)間即可。而σ則可以根據(jù)各組歷史數(shù)據(jù)的原始信號與預(yù)處理后的信號殘差進(jìn)行統(tǒng)計得到。
2)不同階段的退化速率分布圖。a是預(yù)測模型中唯一的參數(shù),體現(xiàn)系統(tǒng)的退化速度,更直接決定了對系統(tǒng)剩余壽命的預(yù)測結(jié)果,因此需要在預(yù)測過程中準(zhǔn)確地估計退化速率的演變規(guī)律。一方面,因為設(shè)備退化是對時間單調(diào)變化的,所以退化特征經(jīng)過平滑和單調(diào)化預(yù)處理后,退化速率a的取值上邊界或下邊界很容易確定(根據(jù)退化特征不同,可得a≤1或a≥1),同時可以根據(jù)歷史數(shù)據(jù)中a的極值情況確定其另一個取值邊界。另一方面,由歷史數(shù)據(jù)統(tǒng)計信息可以給出退化速率a隨系統(tǒng)狀態(tài)x變化的先驗分布p(a|x),作為預(yù)測階段粒子狀態(tài)的指導(dǎo)信息。首先將歷史數(shù)據(jù)的退化速率對應(yīng)到退化特征坐標(biāo)軸上,再利用插值方法計算退化特征軸上的退化速率趨勢曲線。最后對所有趨勢曲線進(jìn)行統(tǒng)計,在每個退化特征點上,對其相應(yīng)的退化速率進(jìn)行正態(tài)分布擬合,可以得到退化速率隨特征變化的分布信息。
在獲取了所需的先驗信息之后,就可以對已有的監(jiān)測數(shù)據(jù)進(jìn)行預(yù)測。首先根據(jù)模型參數(shù)的初始分布信息,對N個粒子樣本初始化,然后依次進(jìn)行迭代估計和外推預(yù)測兩個步驟,最后對所有粒子樣本的狀態(tài)進(jìn)行統(tǒng)計,得到最終的預(yù)測結(jié)果。預(yù)測流程圖如圖3所示。
1)迭代估計階段。首先對由已知的監(jiān)測數(shù)據(jù)提取的退化特征進(jìn)行平滑和單調(diào)化的預(yù)處理。然后利用基于通用預(yù)測模型的粒子濾波方法對模型參數(shù)進(jìn)行迭代估計。在迭代估計過程中,將對退化狀態(tài)值x和退化速率a分別進(jìn)行重采樣。
2)外推預(yù)測階段。利用已知的觀測數(shù)據(jù)完成對模型參數(shù)的迭代估計之后,再結(jié)合預(yù)測模型以及之前得到的退化速率分布圖進(jìn)行外推預(yù)測。設(shè)當(dāng)前退化狀態(tài)預(yù)測值為
那么,根據(jù)p(a|x),隨機(jī)生成當(dāng)前狀態(tài)下的先驗退化速率aik(i=1,2,…,N)。然后再根據(jù)預(yù)測模型,可以得到下一步的退化狀態(tài)值。
依次迭代預(yù)測,當(dāng)退化狀態(tài)預(yù)測值大于或小于失效閾值時,預(yù)測結(jié)束。對所有粒子的當(dāng)前狀態(tài)進(jìn)行統(tǒng)計,就能得到失效時間預(yù)測的概率密度分布以及預(yù)測結(jié)果的置信區(qū)間等信息。
圖3 監(jiān)測數(shù)據(jù)預(yù)測流程圖Fig.3 The flowchart of prediction scheme for monitoring data
圖4所示為一個加速壽命實驗的滾動軸承的剩余使用壽命預(yù)測結(jié)果。在圖4(a)中,前500個數(shù)據(jù)點(數(shù)據(jù)點之間間隔10 s)被用來迭代更新預(yù)測模型,再用更新之后的模型結(jié)合退化速率分布圖,對軸承500個數(shù)據(jù)點之后的狀態(tài)進(jìn)行預(yù)測。實際的失效時間位于預(yù)測結(jié)果的90%置信區(qū)間以內(nèi),表示預(yù)測結(jié)果有效。但因為用于迭代更新的數(shù)據(jù)點相對較少,所以預(yù)測結(jié)果的不確定度大。在圖4(b)中可以看到,隨著已知數(shù)據(jù)點的增加,預(yù)測結(jié)果逐漸逼近實際值,置信區(qū)間的寬度隨已知數(shù)據(jù)點的增加而變窄,表示預(yù)測結(jié)果的準(zhǔn)確性在升高,而不確定性在逐漸減低。
圖4 軸承剩余使用壽命預(yù)測結(jié)果Fig.4 Prediction results of bearing remaining useful life
圖5 電池剩余使用壽命預(yù)測結(jié)果Fig.5 Prediction results of battery remaining useful life
圖5所示為NASA鋰離子電池的循環(huán)壽命預(yù)測結(jié)果。鋰離子電池循環(huán)壽命的退化特征為電池容量,呈逐步下降的趨勢,與軸承的退化特征剛好相反,但對本文提出的預(yù)測方法沒有影響。因為受到松弛效應(yīng)的影響,鋰離子電池的原始退化特征存在明顯的非單調(diào)下降的部分,但采用文獻(xiàn)[16]中的方法消除松弛效應(yīng)之后,即可按照文中提出的預(yù)測框架進(jìn)行處理。從圖5(a)中可以看出,該方法能夠?qū)︿囯x子電池的循環(huán)壽命進(jìn)行有效估計。而在圖5(b)中,可以看到預(yù)測結(jié)果收斂速度更快,準(zhǔn)確性更高。在圖示范圍內(nèi),預(yù)測中值與實際值的最大誤差僅為11個循環(huán)。這是因為相對于圖4,鋰離子電池的預(yù)測特征與系統(tǒng)退化過程相關(guān)性更強(qiáng),而且特征受其他因素影響更小,所以相應(yīng)的預(yù)測效果也更好。
如前所述,通過對預(yù)測特征的退化速率進(jìn)行跟蹤,能夠使粒子濾波對目標(biāo)系統(tǒng)的狀態(tài)估計迅速收斂,在預(yù)測階段也能比較準(zhǔn)確地逼近真實的退化趨勢。預(yù)測結(jié)果的準(zhǔn)確度隨已知數(shù)據(jù)長度逐漸增減,不確定度隨之下降。而且,對于無論是代表劣化程度(逐漸下降)還是代表健康程度(逐漸升高)的退化特征,該預(yù)測框架都同樣適用,表現(xiàn)出了其良好的通用性,并在一定程度上簡化了粒子濾波預(yù)測方法。
但是,對比兩種不同系統(tǒng)的預(yù)測結(jié)果,可以發(fā)現(xiàn)預(yù)測特征仍然是對算法影響很大的一個因素,對于機(jī)械設(shè)備的預(yù)測尤為重要。由振動信號提取的預(yù)測特征往往受外界噪聲等影響嚴(yán)重,如果預(yù)測特征和系統(tǒng)的退化趨勢相關(guān)程度較低,那么算法的預(yù)測效果也會相對下降,所以如何提取與系統(tǒng)高度相關(guān)的預(yù)測特征還值得深入研究。
本文提出了一種基于DRT-PF的預(yù)測框架。該方法對原始預(yù)測特征進(jìn)行平滑和單調(diào)化處理,再利用歷史數(shù)據(jù)中的退化速率信息輔助粒子濾波的迭代更新與外推預(yù)測。由于采用了通用的退化預(yù)測模型,因此該方法理論上對于任何緩變型退化系統(tǒng)的預(yù)測都適用。當(dāng)擁有與系統(tǒng)退化趨勢高度相關(guān)的預(yù)測特征和足夠的歷史數(shù)據(jù)時,該方法能夠有效預(yù)測目標(biāo)系統(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]王小林,程志軍,郭波.基于維納過程金屬化膜電容器的剩余壽命預(yù)測[J].國防科學(xué)技術(shù)大學(xué)學(xué)報,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ī)濾波模型在變速箱剩余壽命預(yù)測中的應(yīng)用研究[J].軍械工程學(xué)院學(xué)報,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.