考慮隨機沖擊影響的自適應(yīng)Wiener過程剩余壽命預(yù)測方法

董青,鄭建飛,*,胡昌華,余銅輝,牟含笑

1. 火箭軍工程大學(xué) 導(dǎo)彈工程學(xué)院,西安 710025 2. 火箭軍駐西安地區(qū)第一軍事代表室,西安 710100

隨著高新技術(shù)的日益發(fā)展,現(xiàn)代工業(yè)設(shè)備正朝著高度自動化、集成化發(fā)展。由于壽命周期中受到振動沖擊、工況切換、機械磨損、化學(xué)侵蝕、負(fù)載變化以及能量消耗等內(nèi)外部因素影響的綜合作用,此類設(shè)備的性能與健康狀態(tài)不可避免的出現(xiàn)退化現(xiàn)象,最終導(dǎo)致設(shè)備失效,甚至造成難以挽回的生命財產(chǎn)及環(huán)境損失。如果能在設(shè)備失效之前預(yù)測其壽命或剩余壽命,進(jìn)而采用相應(yīng)的維護(hù)策略,則可有效避免災(zāi)難性事故,從而達(dá)到降低成本和保護(hù)環(huán)境的目的。作為預(yù)測與健康管理(Prognostics and Health Management, PHM)的核心技術(shù),剩余壽命預(yù)測(Remaining Useful Life, RUL)是對設(shè)備進(jìn)行科學(xué)維修決策的前提和基礎(chǔ)。因此,準(zhǔn)確預(yù)測設(shè)備剩余壽命來降低突發(fā)失效風(fēng)險具有重要的理論意義和工程應(yīng)用價值。

經(jīng)過幾十年的發(fā)展,剩余壽命預(yù)測技術(shù)取得了豐碩的理論成果并得到廣泛應(yīng)用。Pecht在關(guān)于PHM的專著中,將剩余壽命預(yù)測方法主要分為基于機制建模和基于數(shù)據(jù)驅(qū)動建模的方法2類。柴天佑指出復(fù)雜的工業(yè)過程往往具有多變量、強耦合、強非線性、動態(tài)工況變化、難以用數(shù)學(xué)模型描述等復(fù)雜特性,數(shù)據(jù)驅(qū)動的方法主要包括人工智能方法和統(tǒng)計數(shù)據(jù)驅(qū)動的方法。Si等將數(shù)據(jù)驅(qū)動的方法分為基于直接監(jiān)測數(shù)據(jù)的方法和基于間接監(jiān)測數(shù)據(jù)的方法。其中直接監(jiān)測數(shù)據(jù)的方法又可分為基于隨機系數(shù)回歸模型的方法、基于Wiener過程的方法、基于Gamma過程的方法、基于馬爾科夫鏈的方法以及逆高斯過程的方法。相對于其他隨機過程,Wiener過程是一類具有高斯獨立分布增量的非單調(diào)退化過程。憑借其良好的數(shù)學(xué)特性,已被廣泛應(yīng)用于設(shè)備可靠性分析和壽命預(yù)測領(lǐng)域。

在工程實際中,設(shè)備除去正常工作外,還會受到單粒子效應(yīng)、溫度振動沖擊、腐蝕等外部環(huán)境因素的影響,而這些額外因素可以看作廣義的隨機沖擊,對設(shè)備退化過程造成一定影響。張延靜等針對受到隨機沖擊影響的單部件系統(tǒng),考慮沖擊對退化量的影響,建立競爭失效可靠性模型,并運用到維修決策中。孫富強等針對存在沖擊韌性的系統(tǒng),考慮沖擊對退化量和退化率2方面的影響,建立非線性Wiener過程的競爭失效模型進(jìn)行可靠性評估。王浩偉等針對導(dǎo)彈武器系統(tǒng),假定突發(fā)失效概率與整彈性能退化相關(guān),采用Gamma過程和Weibull分布分別作為退化失效模型和突發(fā)失效模型進(jìn)行可靠性評估和RUL預(yù)測。王華偉等利用航空發(fā)動機失效信息,分別對性能退化失效和突發(fā)失效建立RUL預(yù)測模型。通過建立混合Weibull可靠性模型,量化性能退化失效對突發(fā)失效的影響,實現(xiàn)突發(fā)失效退化情況下的航空發(fā)動機剩余壽命預(yù)測。白燦等提出一種考慮隨機沖擊影響的非線性退化設(shè)備RUL預(yù)測方法,在首達(dá)時間意義下給出壽命和剩余壽命的近似解表達(dá)式,極大縮短了計算時間,為后續(xù)開展設(shè)備健康維護(hù)奠定了良好基礎(chǔ)。

基于上述分析,考慮沖擊與退化之間的相互影響已經(jīng)開展許多工作,但大多應(yīng)用于設(shè)備可靠性評估領(lǐng)域,而在壽命預(yù)測領(lǐng)域仍需要進(jìn)一步深入拓展。在現(xiàn)有基于Wiener過程進(jìn)行退化建模和剩余壽命預(yù)測方法研究中,對于設(shè)備壽命周期中現(xiàn)實存在的同時具有隨機沖擊影響、監(jiān)測間隔不均勻、與先驗數(shù)據(jù)測量頻率不一致的情況尚未考慮。同時,對于模型漂移系數(shù)的更新,僅局限于實時監(jiān)測數(shù)據(jù)下的更新,而未考慮未來RUL預(yù)測中該模型自適應(yīng)漂移的可變性。

綜上,本文考慮隨機沖擊對退化過程的影響,提出一種基于自適應(yīng)Wiener過程的非線性退化建模和剩余壽命預(yù)測方法。首先,通過自適應(yīng)Wiener過程描述自身退化過程；其次,通過正態(tài)分布刻畫單次沖擊對設(shè)備退化水平的影響,并將隨機沖擊融入到自身退化建模中,在首達(dá)時間意義下,推導(dǎo)出設(shè)備隨機沖擊影響下剩余壽命分布解析式；再次,通過構(gòu)建狀態(tài)空間模型,基于Kalman濾波框架和期望最大化算法,利用退化數(shù)據(jù)實現(xiàn)參數(shù)和RUL的自適應(yīng)更新；最后,通過數(shù)值仿真、慣導(dǎo)系統(tǒng)陀螺儀以及鋰電池實例,驗證所提方法進(jìn)行剩余壽命預(yù)測的精確性和實用性。

1 考慮隨機沖擊影響的自適應(yīng)Wiener過程退化建模

1.1 自適應(yīng)Wiener過程退化模型

令()表示設(shè)備在運行中的退化過程,現(xiàn)有常用的Wiener過程模型為

()=++()

(1)

式中：、分別為退化過程的漂移系數(shù)和擴(kuò)散系數(shù)；()為標(biāo)準(zhǔn)Brownian運動；為退化初值,通常假設(shè)=0。

由于上述Wiener過程模型存在測量間隔不均勻、測量頻率不一致、忽略自適應(yīng)漂移可變性3點不足,當(dāng)不考慮隨機沖擊影響時,對于設(shè)備退化過程{(),≥0},采用自適應(yīng)Wiener過程來描述

(2)

式中：()為一個符合Wiener過程且隨時間變化的漂移系數(shù);為初始漂移率,>0;為自適應(yīng)漂移項的擴(kuò)散系數(shù);()為設(shè)備監(jiān)測的退化量;為退化過程的擴(kuò)散系數(shù);()為1個獨立于()的標(biāo)準(zhǔn)Brownian運動;為時間的1個參數(shù)；(,)為1個隨時間增長的非線性函數(shù),為了方便表示,后續(xù)可簡寫為()。當(dāng)=0,=1時,設(shè)備呈線性退化,此時退化模型變?yōu)槿缡?1)所示的Wiener過程。

1.2 隨機沖擊模型

現(xiàn)有考慮退化與沖擊相互影響的研究中,大多采用泊松過程來描述系統(tǒng)受到的隨機沖擊過程。相比于其他隨機過程,泊松過程主要有以下3點優(yōu)勢：① 泊松過程是一種點過程,用來表征隨機沖擊這種單事件效應(yīng)現(xiàn)象是合理的；② 泊松過程具有無記憶屬性,即沖擊是隨機發(fā)生的；③ 泊 松密度()可為任意形式,能較好地描述隨機沖擊的出現(xiàn)頻次。

(3)

式中：()為隨機沖擊引起的退化水平變化量之和；為單次沖擊引起的退化水平突變量。

上述隨機沖擊模型已經(jīng)被應(yīng)用于考慮隨機沖擊影響的退化建模中,而本文基于此模型,與非線性自適應(yīng)Wiener過程相結(jié)合,建立退化模型,具體方法在1.3節(jié)展開介紹。

1.3 考慮隨機沖擊影響的退化模型

將隨機沖擊模型融入設(shè)備退化模型中,如圖1 所示。圖1中:橫坐標(biāo)為退化時間,縱坐標(biāo)為退化量(),為失效閾值,(=1,2,3,…)為第次狀態(tài)監(jiān)測時間,為設(shè)備剩余壽命分布。當(dāng)()達(dá)到預(yù)先設(shè)定的失效閾值時,即認(rèn)定設(shè)備失效。

圖1 考慮隨機沖擊影響時設(shè)備退化過程Fig.1 Degradation process of equipment with random shook

用隨機過程{(),>0}表示設(shè)備退化過程,考慮隨機沖擊影響,設(shè)備實際退化過程分為2部分：設(shè)備正常退化和隨機沖擊對退化水平的影響。此時,設(shè)備退化模型為

(4)

2 隨機沖擊影響下設(shè)備剩余壽命預(yù)測

2.1 設(shè)備剩余壽命預(yù)測

在隨機退化建?？蚣芟?根據(jù)首達(dá)時間定義設(shè)備的壽命為

=inf{：()≥|(0)<)

(5)

式中：()為設(shè)備退化狀態(tài)；(0)為退化初值;為設(shè)定的失效閾值。

由于Brownian運動的隨機性,壽命為服從逆高斯分布的隨機變量,根據(jù)文獻(xiàn)[15]自適應(yīng)Wiener過程壽命分布的概率密度函數(shù)(Probability Density Function, PDF)為

(6)

若～(,),、都為常數(shù),為正實數(shù),則以下期望公式成立：

(7)

根據(jù)引理1,可通過全概率公式,得到考慮隨機沖擊影響下自適應(yīng)Wiener過程的壽命分布PDF為

(8)

根據(jù)隨機過程首達(dá)時間的概念,當(dāng){(),≥0}首次達(dá)到失效閾值時,即認(rèn)為設(shè)備失效。因此,基于監(jiān)測數(shù)據(jù)0:={,,…,},將設(shè)備在時刻的剩余壽命定義為

=inf{：(+)≥|()<)

(9)

由式(8)可得,考慮隨機沖擊影響時自適應(yīng)Wiener過程的剩余壽命分布PDF為

(10)

2.2 設(shè)備剩余壽命預(yù)測在線更新

為實現(xiàn)同時考慮時變不確定性、個體差異性、隨機沖擊影響等因素時設(shè)備RUL預(yù)測的在線更新,需構(gòu)建退化過程的離散模型。此外,當(dāng)受到隨機沖擊時,設(shè)備的退化速率可能會受到不同程度的影響,并且同樣的沖擊對不同退化速率設(shè)備的影響也不盡相同。因此,在自適應(yīng)Wiener過程中,表征設(shè)備退化速率的漂移系數(shù)()和隨機沖擊影響的退化量可能會相互影響,且該影響隨著設(shè)備運行的時間不斷變化。鑒于此,本節(jié)將式(4)中的退化模型改寫為

(11)

(12)

2) 狀態(tài)估計

(13)

3) 協(xié)方差更新

(14)

若～(,)的元正態(tài)分布,且、為已知的常實數(shù),、為已知的維常值向量,則

(15)

根據(jù)引理2,通過全概率公式可求得同時考慮隨機沖擊和個體差異性時,首達(dá)時間意義下設(shè)備RUL分布的PDF為

(16)

基于泊松分布的概率密度函數(shù),產(chǎn)生個隨機變量表示隨機沖擊次數(shù),結(jié)果記為(),(),…,()，有

(17)

計算步驟3得到RUL分布的均值,令其作為隨機沖擊影響下設(shè)備剩余壽命的PDF

(18)

3 參數(shù)估計

(19)

(20)

E步：

(21)

M步：

(22)

根據(jù)狀態(tài)空間模型,基于貝葉斯定理和條件概率公式,未知參數(shù)的對數(shù)聯(lián)合似然函數(shù)為

(23)

基于構(gòu)建的狀態(tài)空間模型,以上條件期望值可通過RTS(Rauch-Tung-Striebel)平滑濾波計算。具體步驟如下：

1) 平滑濾波

(24)

2) 協(xié)方差初始值

(25)

3) 后向迭代

(26)

(27)

交替執(zhí)行E步和M步,直至估計的結(jié)果滿足預(yù)設(shè)的算法終止條件,即最大迭代次數(shù)或參數(shù)收斂。當(dāng)有新的狀態(tài)監(jiān)測數(shù)據(jù)時,將之前估計得到的參數(shù)作為EM算法的初值,并基于新監(jiān)測數(shù)據(jù)執(zhí)行EM算法,以獲得相關(guān)參數(shù)的最新估計值。

4 實例驗證

4.1 數(shù)值仿真

由表1可見,相比于白燦等所提方法、不考慮隨機沖擊影響的方法,本文所提方法的參數(shù)估計結(jié)果更為準(zhǔn)確,說明了本文所提方法的合理性和有效性,為后續(xù)準(zhǔn)確預(yù)測剩余壽命奠定基礎(chǔ)。

表1 模型參數(shù)估計值Table 1 Estimated values of model parameters

此外,設(shè)備退化過程首次超過預(yù)先設(shè)定的閾值即為首達(dá)時間,經(jīng)過大量模擬仿真可得到壽命分布直方圖如圖2所示,并與白燦等所提方法得到的壽命分布PDF、本文所提方法得到的壽命分布PDF進(jìn)行比較。為了后續(xù)便于表示,令本文所提模型記為M0,白燦等提出的隨機沖擊退化模型記為M1,不考慮隨機沖擊影響的退化模型記為M2。

圖2 不同模型壽命分布PDFFig.2 PDF of life distribution with different models

圖3～圖5的結(jié)果反映了隨機沖擊對壽命分布的影響,即： ① 增大,加速設(shè)備退化,設(shè)備壽命隨之不斷下降,且壽命預(yù)測的不確定性增高； ② 增大,壽命分布的均值幾乎不變,但壽命預(yù)測的不確定性增高； ③ 減小,設(shè)備壽命下降且預(yù)測不確定性也隨之下降。

圖3 不同均值uI下壽命分布PDFFig.3 PDF of life distribution with different mean values

圖4 不同方差下壽命分布PDFFig.4 PDF of life distribution with different variances

圖5 不同泊松密度μ下壽命分布PDFFig.5 PDF of life distribution with different Poisson densities

4.2 陀螺儀實例分析

陀螺儀是導(dǎo)彈等運載體慣性導(dǎo)航系統(tǒng)的核心設(shè)備,其性能的優(yōu)劣直接決定了導(dǎo)航精度和任務(wù)的成敗,同時具有結(jié)構(gòu)復(fù)雜、精度要求高、價格昂貴等特點。通常,陀螺儀在安裝到導(dǎo)彈后,會隨彈經(jīng)歷長時間儲存、測試、拆裝、運輸?shù)炔煌h(huán)境過程,易受到隨機振動沖擊等影響,對表征陀螺儀精度性能的陀螺儀誤差模型系數(shù)的退化過程產(chǎn)生加速退化影響,最終影響陀螺儀服役時間。因此,為了準(zhǔn)確預(yù)測陀螺儀剩余壽命,為后續(xù)維護(hù)管理提供科學(xué)依據(jù),必須考慮隨機沖擊對陀螺儀退化的影響,某機械轉(zhuǎn)子陀螺儀在壽命周期中漂移系數(shù)退化數(shù)據(jù)如圖6所示,其中的數(shù)據(jù)突變是由運輸振動沖擊導(dǎo)致的退化加速。

在工程實際中,該陀螺儀漂移系數(shù)超過預(yù)先設(shè)定閾值=0.37 (°)/h,認(rèn)為其不再滿足使用需求,即認(rèn)為發(fā)生失效故障。為了驗證所提模型RUL預(yù)測的有效性,用前67個數(shù)據(jù)估計模型參數(shù),并用剩余數(shù)據(jù)對陀螺儀的RUL更新。利用上一節(jié)提出的EM算法估計未知參數(shù)。

圖6 陀螺儀退化數(shù)據(jù)Fig.6 Gyro degradation data

圖7 參數(shù)更新過程Fig.7 Update process of parameter

圖8 參數(shù)P0|0更新過程Fig.8 Update process of parameter P0|0

圖9 陀螺儀漂移數(shù)據(jù)擬合效果Fig.9 Gyro drift data fitting effect

在此,給出3種模型各自的RUL預(yù)測結(jié)果,具體如圖10、圖11所示。結(jié)果表明,相比于模型M2,模型M0考慮了設(shè)備運行過程中隨機沖擊對退化量和退化率的影響,更加符合實際退化情況,因此能顯著提高預(yù)測準(zhǔn)確度；與模型M1相比,模型M0基于非線性自適應(yīng)Wiener過程進(jìn)行退化建模,能克服測量間隔不均勻、測量頻率不一致的影響,同時該模型將漂移系數(shù)視為Wiener過程,在未來RUL預(yù)測時漂移參數(shù)自適應(yīng)變化。此外,MO模型考慮到同一批次設(shè)備不同個體間退化率的差異性的影響,并將隨機沖擊對退化率的影響融入狀態(tài)空間模型中,實現(xiàn)RUL在線更新。因此,本文所提模型M0剩余壽命分布的PDF窄而尖銳,預(yù)測結(jié)果的不確定性較小,能獲得更好的預(yù)測結(jié)果。

圖10 不同模型RUL預(yù)測結(jié)果Fig.10 RUL prediction results of different models

圖11 不同監(jiān)測點RULFig.11 RUL prediction at different monitoring points

為了對不同模型預(yù)測結(jié)果進(jìn)行量化比較,進(jìn)而驗證所提模型M0的有效性。定義設(shè)備的均方根誤差(Root Mean Squared Error, RMSE)為

(28)

圖12表明,隨著設(shè)備不斷退化,預(yù)測結(jié)果的RMSE也隨之減小。所提模型M0對應(yīng)RMSE最小,進(jìn)一步證明了M0較M1、M2具有更好的模型擬合性,能有效降低預(yù)測結(jié)果的不確定性。

圖12 不同模型RUL預(yù)測的均方根誤差Fig.12 Root mean square error of RUL prediction with different models

為了定量分析預(yù)測結(jié)果,在此引入絕對誤差(Absolute Error, AE)和評分函數(shù)(Scoring Function)進(jìn)行性能評價。評分函數(shù)的計算公式為

(29)

由表2可知,3種模型RUL預(yù)測結(jié)果的絕對誤差在不斷減小,且模型M0對應(yīng)的絕對誤差最小,即預(yù)測結(jié)果更準(zhǔn)確。評分函數(shù)主要用于描述預(yù)測值與真實值高估與低估的關(guān)系。一般評分函數(shù)分值越小,預(yù)測結(jié)果越準(zhǔn)確。由式(29)可以計算得到,3種模型對應(yīng)結(jié)果分別為0.124 1、4.816 9、25.097 0,模型M0對應(yīng)的評分函數(shù)分值最小。接下來,采用美國國家航空航天局(NASA)公開的鋰電池數(shù)據(jù)集進(jìn)一步驗證本文所提方法的有效性和通用性。

表2 不同監(jiān)測點絕對誤差Table 2 Absolute error at different monitoring points

4.3 鋰電池實例分析

鋰電池的容量退化數(shù)據(jù)隨著充放電循環(huán)而不斷降低,在退化過程中,由于充放電循環(huán)的中止,電池在沒有電子壓力的情況下,其容量數(shù)據(jù)可能會有所恢復(fù)。而這種容量的恢復(fù)可認(rèn)為一種負(fù)向的隨機沖擊施加到退化過程中。本文從公開的鋰電池數(shù)據(jù)#5、#6、#7、#18中,選擇#5電池進(jìn)行驗證,鋰電池退化數(shù)據(jù)如圖13所示。在工程實際中,一般認(rèn)為電池容量衰減30%將失效。

圖13 鋰電池退化數(shù)據(jù)Fig.13 Lithium battery degradation data

通過初始電池容量值減去每次充放電循環(huán)后的電池容量,得到符合Wiener過程的鋰電池真實退化數(shù)據(jù)和預(yù)測退化數(shù)據(jù),具體如圖14所示。圖14結(jié)果表明,本文方法預(yù)測數(shù)據(jù)與真實數(shù)據(jù)擬合效果較好,能較好地反映鋰電池退化情況。

圖14 鋰電池退化數(shù)據(jù)擬合效果Fig.14 Lithium battery degradation data fitting effect

根據(jù)參數(shù)估計結(jié)果代入模型RUL預(yù)測解析式中,可得到模型M0和模型M1的RUL預(yù)測結(jié)果,具體如圖15、圖16所示。

由圖15、圖16可見,與模型M1相比,本文所提模型M0的預(yù)測效果更為準(zhǔn)確。尤其在退化初期,模型M0能動態(tài)調(diào)整漂移參數(shù),預(yù)測精度較高。此外,模型M0同時考慮時變不確定性、同一批次不同設(shè)備間的個體差異性以及隨機沖擊影響,考慮情況更加全面,預(yù)測結(jié)果也更為準(zhǔn)確。

圖15 不同模型RUL預(yù)測結(jié)果Fig.15 RUL prediction results of different models

圖16 不同監(jiān)測點RULFig.16 RUL prediction at different monitoring points

同樣地,為了定性和定量分析本文所提方法的有效性,給模型M0、M1的RUL預(yù)測結(jié)果的均方根誤差和絕對誤差,分別如圖17、表3所示。

圖17 不同模型RUL預(yù)測的均方根誤差Fig.17 Root mean square error of RUL prediction with different models

表3 不同監(jiān)測點絕對誤差Table 3 Absolute error at different monitoring points

圖17表明,隨著設(shè)備運行時間的推進(jìn),模型M0和M1預(yù)測結(jié)果的RMSE不斷減小。在模型M0和M1中,所提模型M0對應(yīng)RMSE較小,表明模型M0能有效降低預(yù)測結(jié)果的不確定性。

由表3可見,模型M0的絕對誤差始終小于模型M1,其預(yù)測結(jié)果更為準(zhǔn)確。綜上,本文所提模型M0預(yù)測準(zhǔn)確度高、不確定性小,并且具有一定的適用性,可為后續(xù)設(shè)備的維護(hù)決策提供理論依據(jù)。

5 結(jié) 論

針對考慮隨機沖擊影響的隨機退化設(shè)備,本文提出一種基于非線性自適應(yīng)Wiener過程的剩余壽命預(yù)測方法。

1) 利用非線性自適應(yīng)Wiener過程和正態(tài)分布分別描述設(shè)備正常退化和隨機沖擊對設(shè)備退化量的影響,建立了能有效融合隨機沖擊影響的自適應(yīng)Wiener過程隨機退化模型。

2) 在首達(dá)時間意義下推導(dǎo)出剩余壽命解析式,構(gòu)建狀態(tài)空間模型實現(xiàn)RUL在線更新。通過仿真實驗以及慣性導(dǎo)航系統(tǒng)陀螺儀、鋰電池數(shù)據(jù)驗證所提方法的有效性。

3) 本文所提方法基于非線性自適應(yīng)Wiener過程,能夠克服退化設(shè)備測量間隔分布不均勻、監(jiān)測數(shù)據(jù)的測量頻率與歷史數(shù)據(jù)頻率不一致的情況,并且考慮將來退化過程中自適應(yīng)漂移的可變性以及隨機沖擊對設(shè)備退化的影響,可為存在隨機沖擊影響的設(shè)備剩余壽命預(yù)測提供理論依據(jù),且具有一定的潛在應(yīng)用價值。