基于半隨機(jī)濾波的線(xiàn)性退化設(shè)備剩余壽命預(yù)測(cè)

朱 旭，司小勝，胡昌華，莫仁鵬

(火箭軍工程大學(xué)，西安 710000)

0 引言

隨著工程實(shí)際中設(shè)備復(fù)雜度的不斷增高，系統(tǒng)因故障失效導(dǎo)致的損失也變得愈發(fā)不可承受，嚴(yán)重的設(shè)備故障不僅會(huì)造成經(jīng)濟(jì)損失，甚至?xí)a(chǎn)生人員傷亡。因此，追求更高精度的設(shè)備剩余壽命預(yù)測(cè)，成為提高設(shè)備可靠性和安全性、進(jìn)行合理的視情維修的關(guān)鍵。

現(xiàn)有的剩余壽命(RUL)預(yù)測(cè)研究，主要針對(duì)軍事裝備[1]、鋰電池[2]、大型民用機(jī)械(如船舶柴油機(jī))[3]等設(shè)備，由此發(fā)展出了基于機(jī)理模型的方法、基于機(jī)器學(xué)習(xí)的方法和統(tǒng)計(jì)數(shù)據(jù)驅(qū)動(dòng)的方法。其中，統(tǒng)計(jì)數(shù)據(jù)驅(qū)動(dòng)的設(shè)備剩余壽命預(yù)測(cè)方法中研究最多的是基于隨機(jī)退化過(guò)程建模的方法，這類(lèi)方法由于基于概率論和隨機(jī)過(guò)程理論，通過(guò)預(yù)測(cè)設(shè)備性能退化變量到達(dá)失效閾值的時(shí)間實(shí)現(xiàn)剩余壽命的預(yù)測(cè)，具有良好的數(shù)學(xué)性質(zhì)和統(tǒng)計(jì)特性，能夠較準(zhǔn)確地描述設(shè)備退化和壽命分布的數(shù)學(xué)特征，因而得到了廣泛的應(yīng)用。在最新的研究中，針對(duì)隨機(jī)退化設(shè)備剩余壽命預(yù)測(cè)問(wèn)題已經(jīng)拓展到兩階段[4]、非線(xiàn)性[5]、多重不確定性[6]和多源信息融合[7]等情形。然而，當(dāng)前的研究普遍未能解決失效閾值如何確定的問(wèn)題。盡管最近的研究中提出了隨機(jī)失效閾值[8]和多環(huán)境實(shí)驗(yàn)[9]等方法，但對(duì)失效閾值的確定仍然依賴(lài)退化實(shí)驗(yàn)，并且這些方法通常只考慮了歷史監(jiān)測(cè)信息，難以將失效數(shù)據(jù)有效地應(yīng)用于剩余壽命預(yù)測(cè)模型的建立。

在統(tǒng)計(jì)數(shù)據(jù)驅(qū)動(dòng)的剩余壽命預(yù)測(cè)方法之中，文獻(xiàn)[10]提出的半隨機(jī)濾波方法不需要失效閾值的相關(guān)信息，可以充分利用失效數(shù)據(jù)，并基于貝葉斯濾波理論將歷史監(jiān)測(cè)信息應(yīng)用于剩余壽命預(yù)測(cè)模型的建立，因此在應(yīng)用中具有明顯的優(yōu)勢(shì)。然而，該方法預(yù)測(cè)精度受制于半隨機(jī)濾波模型狀態(tài)變量的初始分布，即壽命分布?，F(xiàn)有的半隨機(jī)濾波相關(guān)研究中，普遍缺乏對(duì)壽命分布進(jìn)行估計(jì)的理論依據(jù)，沒(méi)有充分利用設(shè)備退化先驗(yàn)知識(shí)。在隨后對(duì)軸承等的研究[11]當(dāng)中，采用了數(shù)據(jù)擬合的方法對(duì)壽命分布進(jìn)行估計(jì)；隨后，文獻(xiàn)[12]在對(duì)風(fēng)電機(jī)組齒輪箱的剩余壽命預(yù)測(cè)中運(yùn)用核密度估計(jì)的方法提高了壽命分布的擬合精度，但該方法計(jì)算量較大且只能數(shù)值求解。

可以看出，基于隨機(jī)退化過(guò)程建模的方法與半隨機(jī)濾波方法互有優(yōu)勢(shì)，若能結(jié)合兩類(lèi)方法的優(yōu)勢(shì)，將有望提高對(duì)設(shè)備剩余壽命的預(yù)測(cè)精度。為此，本文提出了一種考慮設(shè)備退化先驗(yàn)知識(shí)的基于半隨機(jī)濾波的線(xiàn)性隨機(jī)退化設(shè)備剩余壽命預(yù)測(cè)方法，該方法依據(jù)線(xiàn)性隨機(jī)退化設(shè)備的退化失效特點(diǎn)，將描述線(xiàn)性隨機(jī)退化設(shè)備壽命的逆高斯分布作為半隨機(jī)濾波模型的初始分布，然后利用半隨機(jī)濾波技術(shù)預(yù)測(cè)設(shè)備的剩余壽命。利用歷史監(jiān)測(cè)數(shù)據(jù)和壽命數(shù)據(jù)通過(guò)極大似然方法估計(jì)確定模型參數(shù)。最后，基于疲勞裂紋增長(zhǎng)數(shù)據(jù)進(jìn)行了驗(yàn)證和對(duì)比，結(jié)果表明所提方法能夠有效提高剩余壽命預(yù)測(cè)精度。

1 符號(hào)定義與模型假設(shè)

為便于描述半隨機(jī)濾波方法的相關(guān)理論，將本文所使用的半隨機(jī)濾波方法預(yù)測(cè)模型的符號(hào)及其定義描述如下：tk為第k個(gè)監(jiān)測(cè)點(diǎn)的時(shí)間，稱(chēng)為tk時(shí)刻；xk為tk時(shí)刻的設(shè)備條件剩余壽命；yk為tk時(shí)刻獲得的設(shè)備監(jiān)測(cè)信息；Y1:k為從初始時(shí)刻到tk時(shí)刻的設(shè)備歷史監(jiān)測(cè)信息，Y1:k={y1,y2,…,yk}；p0(x0)為設(shè)備在t0時(shí)刻的剩余壽命概率密度函數(shù)，即設(shè)備壽命的概率密度函數(shù)；p(yk|xk)為tk時(shí)刻xk條件下yk的條件概率密度；p(xk|Y1:k)為tk時(shí)刻在歷史監(jiān)測(cè)信息Y1:k的條件下設(shè)備剩余壽命xk的條件概率密度。

根據(jù)工程實(shí)踐中的經(jīng)驗(yàn)和本文建模的需要，提出以下假設(shè)：1)對(duì)設(shè)備的監(jiān)測(cè)是在離散時(shí)間上進(jìn)行的；2)tk時(shí)刻的條件剩余壽命xk是一個(gè)隨機(jī)變量，以概率密度函數(shù)描述tk時(shí)刻的監(jiān)測(cè)信息yk與xk的關(guān)系，其僅與xk有關(guān)；3)相鄰兩次監(jiān)測(cè)期間，不會(huì)發(fā)生設(shè)備維修、替換等操作。

對(duì)于假設(shè)1)，對(duì)設(shè)備健康狀態(tài)的監(jiān)測(cè)是在離散時(shí)間上進(jìn)行的，該假設(shè)與監(jiān)測(cè)設(shè)備以及工程實(shí)際有關(guān)，也是工程實(shí)際中普遍的情況;對(duì)于假設(shè)2)，條件剩余壽命的隨機(jī)性和監(jiān)測(cè)信息與條件剩余壽命之間的隨機(jī)關(guān)系是符合先驗(yàn)知識(shí)的，只有存在隨機(jī)關(guān)系才能通過(guò)監(jiān)測(cè)數(shù)據(jù)預(yù)測(cè)設(shè)備的剩余壽命;對(duì)于假設(shè)3)，設(shè)備維修和備件替換會(huì)使得當(dāng)前時(shí)刻為止的歷史信息與替換維修后的設(shè)備無(wú)相關(guān)性，因此假定為無(wú)維修替換。綜上所述，本文的假設(shè)是合理有依據(jù)的，并且將在后續(xù)的模型建立中得到使用。

2 半隨機(jī)濾波模型的建立

在文獻(xiàn)[10]模型中，在tk-1到tk時(shí)刻的時(shí)間間隔內(nèi)，假如設(shè)備未失效，基于假設(shè)相鄰監(jiān)測(cè)期間未發(fā)生維修，則將相鄰時(shí)刻的剩余壽命關(guān)系描述為

(1)

由此可知，xk的存在性?xún)H依賴(lài)于xk-1≥tk-tk-1。將剩余壽命減少量直接描述為監(jiān)測(cè)時(shí)間間隔，并將剩余壽命視作狀態(tài)量，xk和xk-1之間所滿(mǎn)足的確定性關(guān)系即該方法被稱(chēng)為半隨機(jī)濾波的原因。

為估計(jì)tk時(shí)刻的剩余壽命xk，求得p(xk|Y1:k)是本文的核心目的?？蓪⒈磉_(dá)條件剩余壽命的概率密度函數(shù)p(xk|Y1:k)描述為

(2)

通過(guò)式(2)可見(jiàn)，為求得給定到當(dāng)前時(shí)刻tk的歷史信息Y1:k的條件下設(shè)備剩余壽命xk的概率密度函數(shù)，需要p(xk,yk|Y1:k-1)和p(yk|Y1:k-1)的函數(shù)?；诩僭O(shè)2)所提的條件，可得

p(xk,yk|Y1:k-1)=p(yk|xk,Y1:k-1)pk(xk|Y1:k-1)=p(yk|xk)pk(xk|Y1:k-1)

(3)

通過(guò)積分可推導(dǎo)出

(4)

由此，根據(jù)式(1)可得

(5)

將式(3)～(5)代入式(2)可得

(6)

即給定歷史信息條件下設(shè)備剩余壽命的概率密度函數(shù)式，式(1)～(6)構(gòu)成了基于半隨機(jī)濾波方法進(jìn)行設(shè)備剩余壽命預(yù)測(cè)的過(guò)程。

依據(jù)上述的半隨機(jī)濾波預(yù)測(cè)模型，通過(guò)初始時(shí)刻的壽命分布p0(x0)和所有的歷史監(jiān)測(cè)信息Y1:k，可迭代得到當(dāng)前時(shí)刻的條件剩余壽命分布的表達(dá)式。

3 模型與數(shù)據(jù)擬合

在基于半隨機(jī)濾波方法的設(shè)備剩余壽命預(yù)測(cè)中，需要壽命分布p0(x0)和后驗(yàn)概率p(yk|xk)的分布形式。布朗運(yùn)動(dòng)驅(qū)動(dòng)的線(xiàn)性漂移過(guò)程被稱(chēng)為維納(Wiener)過(guò)程，該過(guò)程因具有良好的數(shù)學(xué)性質(zhì)而被廣泛應(yīng)用于線(xiàn)性隨機(jī)退化設(shè)備的可靠性分析和剩余壽命預(yù)測(cè)領(lǐng)域[1]。工程中存在一類(lèi)線(xiàn)性退化設(shè)備，其退化過(guò)程可以用線(xiàn)性Wiener過(guò)程描述，那么它達(dá)到失效閾值的時(shí)間服從逆高斯分布，即壽命分布為逆高斯分布。經(jīng)過(guò)國(guó)內(nèi)外學(xué)者的廣泛研究，逆高斯分布描述的設(shè)備壽命分布已具有很廣泛的適用性。傳統(tǒng)的半隨機(jī)濾波方法中，受限于設(shè)備先驗(yàn)知識(shí)的缺乏，半隨機(jī)濾波方法預(yù)測(cè)精度又依賴(lài)于初始分布的準(zhǔn)確性，因此，選取了擬合性較好的威布爾分布以進(jìn)行具有一般性的數(shù)據(jù)壽命分布描述。然而，該方法對(duì)于這類(lèi)具有一定的退化過(guò)程先驗(yàn)知識(shí)的數(shù)據(jù)仍使用威布爾分布描述壽命分布，使得壽命分布的描述不夠準(zhǔn)確，進(jìn)而影響半隨機(jī)濾波方法的剩余壽命預(yù)測(cè)結(jié)果精度。

為解決該問(wèn)題，本文首先解決影響預(yù)測(cè)結(jié)果精度的壽命分布即初始分布的確定問(wèn)題。在Wiener過(guò)程中，設(shè)備的壽命分布被描述為服從參數(shù)為ω，λ和σ的逆高斯分布，然而，由于現(xiàn)實(shí)中存在設(shè)備失效閾值往往難以確定的問(wèn)題，在本文中將設(shè)備壽命分布描述為兩參數(shù)的逆高斯分布，即

(7)

式中：x0為初始時(shí)刻設(shè)備的剩余壽命，即設(shè)備壽命；λ和μ為描述分布的參數(shù)，兩參數(shù)的逆高斯分布描述形式規(guī)避了對(duì)失效閾值的依賴(lài)。

其次，半隨機(jī)濾波方法所需的后驗(yàn)概率密度p(yk|xk)是獲得條件剩余壽命分布的必要條件。由于威布爾分布具有良好的擬合特性，為不失一般性，在本文中將后驗(yàn)概率密度p(yk|xk)描述為兩參數(shù)的威布爾分布，即

p(yk|xk)=ρη(ρη)η-1exp{-(ρη)η}

(8)

在實(shí)際中，退化越嚴(yán)重，設(shè)備的剩余壽命越小。因此，設(shè)備的退化量和剩余壽命間存在反比關(guān)系。為描述這一反比關(guān)系，本文中選取ρ=1/(A+Be-Cxk)，使得監(jiān)測(cè)信息和剩余壽命之間滿(mǎn)足E(yk)∝(A+Be-Cxk)。

最后，根據(jù)以上模型設(shè)置與描述，線(xiàn)性隨機(jī)退化設(shè)備tk時(shí)刻的條件剩余壽命概率密度函數(shù)可以表示為

(9)

4 參數(shù)估計(jì)

為了實(shí)現(xiàn)剩余壽命預(yù)測(cè)方法，需要根據(jù)設(shè)備的歷史監(jiān)測(cè)數(shù)據(jù)對(duì)壽命分布p0(x0)的參數(shù)λ和μ以及后驗(yàn)概率p(yk|xk)中參數(shù)A，B，C和η進(jìn)行估計(jì)。為此，本文基于極大似然估計(jì)方法分別對(duì)這類(lèi)參數(shù)分布進(jìn)行估計(jì)。

4.1 p0(x0)的參數(shù)估計(jì)

根據(jù)設(shè)備歷史監(jiān)測(cè)數(shù)據(jù)，利用完整壽命數(shù)據(jù)tj f和截?cái)鄩勖鼣?shù)據(jù)tl可以依據(jù)p0(x0)構(gòu)建似然函數(shù)，本文由于采用逆高斯分布描述p0(x0)，因此似然函數(shù)可表示為

(10)

式中：參數(shù)n為完整壽命數(shù)據(jù)的組數(shù)；m為截?cái)鄩勖鼣?shù)據(jù)的組數(shù)；tj f為第j組數(shù)據(jù)的失效時(shí)間；tl為第i組數(shù)據(jù)的最后監(jiān)測(cè)時(shí)間。然后，通過(guò)使式(10)最大化來(lái)確定模型參數(shù)λ和μ。

4.2 p(yk|xk)的參數(shù)估計(jì)

對(duì)于后驗(yàn)概率密度p(yk|xk)中的參數(shù)，可以根據(jù)設(shè)備的歷史監(jiān)測(cè)信息數(shù)據(jù)構(gòu)建似然函數(shù)，即

(11)

式中:j為數(shù)據(jù)組數(shù);nj為第j組數(shù)據(jù)的監(jiān)測(cè)點(diǎn)個(gè)數(shù)；假如失效時(shí)間存在的話(huà)，則tj f表示第j組的失效時(shí)間。然后，通過(guò)最大化式(11)中的似然函數(shù)，可得到A,B,C和η的極大似然估計(jì)。

5 實(shí)例驗(yàn)證

為驗(yàn)證本文提出的方法，以文獻(xiàn)[13]提供的疲勞裂紋數(shù)據(jù)為依托進(jìn)行驗(yàn)證研究，該退化數(shù)據(jù)共有21組，在間隔相同的12個(gè)監(jiān)測(cè)時(shí)間獲取了監(jiān)測(cè)信息，監(jiān)測(cè)時(shí)間的單位為百萬(wàn)周期，退化數(shù)據(jù)的單位為m。選取其中一組數(shù)據(jù)用于檢驗(yàn)剩余壽命預(yù)測(cè)模型的準(zhǔn)確性，其余的數(shù)據(jù)用于第4章的參數(shù)估計(jì)。圖1所示為疲勞裂紋數(shù)據(jù)隨時(shí)間變化的曲線(xiàn)。

圖1 疲勞裂紋增長(zhǎng)數(shù)據(jù)Fig.1 Fatigue crack growth data

由圖1所示的數(shù)據(jù)隨時(shí)間變化的特點(diǎn)可以看出，其退化趨勢(shì)具有明顯的線(xiàn)性規(guī)律。因此，若采用線(xiàn)性Wiener過(guò)程建模疲勞裂紋數(shù)據(jù)，對(duì)應(yīng)的壽命分布即為逆高斯分布。因此，采用現(xiàn)有的半隨機(jī)濾波方法進(jìn)行剩余壽命預(yù)測(cè)時(shí)，文獻(xiàn)[10]模型方法使用威布爾分布擬合壽命分布時(shí)就會(huì)因?yàn)闆](méi)有充分利用先驗(yàn)知識(shí)而導(dǎo)致擬合結(jié)果不精確，進(jìn)而導(dǎo)致預(yù)測(cè)結(jié)果誤差較大的問(wèn)題。為克服這一問(wèn)題，本文將半隨機(jī)濾波模型中的初始?jí)勖植加媚娓咚狗植济枋?，使用疲勞裂紋增長(zhǎng)數(shù)據(jù)以驗(yàn)證模型的有效性，并與傳統(tǒng)的半隨機(jī)模型預(yù)測(cè)結(jié)果進(jìn)行比較。在本文后續(xù)的描述當(dāng)中，將文獻(xiàn)[10]模型方法描述為M1，將本文的改進(jìn)模型方法描述為M2。

5.1 參數(shù)估計(jì)

為了應(yīng)用半隨機(jī)模型，首先估計(jì)模型參數(shù)。依據(jù)先驗(yàn)知識(shí)，當(dāng)疲勞裂紋長(zhǎng)度增長(zhǎng)到0.040 64 m時(shí)就停止了后續(xù)監(jiān)測(cè)，認(rèn)為設(shè)備失效。因此，基于這一失效閾值，可以得到數(shù)據(jù)集中各個(gè)設(shè)備的失效時(shí)間，將其視為設(shè)備的壽命數(shù)據(jù)。具體地，將11組完整壽命數(shù)據(jù)和9組截?cái)鄶?shù)據(jù)代入式(10)，極大化后求得參數(shù)的估計(jì)值為λ=18.195 9，μ=0.116 3。再對(duì)后驗(yàn)概率的參數(shù)進(jìn)行估計(jì)，應(yīng)用式(11)，通過(guò)Matlab的多維搜索，即可得出后驗(yàn)概率p(yk|xk)中的參數(shù)估計(jì)值：A=-0.629 6，B=2.284 1，C=3.636 1，η=24.363 8。

5.2 方法比較

為了驗(yàn)證模型方法(M1)和本文方法(M2)在改進(jìn)后對(duì)壽命分布描述的準(zhǔn)確性，將參數(shù)估計(jì)得到的結(jié)果代回兩種模型方法中進(jìn)行對(duì)比。圖2所示為兩種模型方法下設(shè)備壽命分布概率密度函數(shù)曲線(xiàn)。由于半隨機(jī)濾波將剩余壽命減少量描述為監(jiān)測(cè)時(shí)間間隔，因此依據(jù)式(1)可知，檢驗(yàn)數(shù)據(jù)對(duì)應(yīng)的設(shè)備壽命為0.12個(gè)百萬(wàn)周期。

圖2 壽命分布概率密度函數(shù)Fig.2 PDF of life distribution

由圖2可以看出，模型方法(M1)的預(yù)測(cè)結(jié)果在對(duì)壽命分布的擬合上是不夠準(zhǔn)確的，比本文方法(M2)的預(yù)測(cè)結(jié)果要差一些，說(shuō)明了其精度不足，也體現(xiàn)出了傳統(tǒng)模型中壽命分布選擇不準(zhǔn)確問(wèn)題。

圖3所示為兩種模型方法預(yù)測(cè)的剩余壽命的條件概率密度p(xk|Y1:k)。

圖3 剩余壽命的條件概率密度函數(shù)Fig.3 PDF of the remaining useful life

由圖3可以看出，改進(jìn)模型方法的概率密度函數(shù)方差要小于傳統(tǒng)方法，說(shuō)明本文提出方法的預(yù)測(cè)不確定性要明顯好于傳統(tǒng)方法。

圖4所示為兩種模型方法在各個(gè)時(shí)刻的剩余壽命預(yù)測(cè)值的對(duì)比。

圖4 剩余壽命的均值Fig.4 Mean value of the remaining useful life

由圖4可知，本文改進(jìn)方法(M2)的預(yù)測(cè)結(jié)果優(yōu)于已有的模型方法(M1)，可以看出,改進(jìn)方法的預(yù)測(cè)結(jié)果與真實(shí)值之間更為接近。并且，隨著時(shí)間推移，新的監(jiān)測(cè)信息對(duì)預(yù)測(cè)結(jié)果修正以后，預(yù)測(cè)值愈發(fā)貼近真實(shí)值。

圖5所示為剩余壽命預(yù)測(cè)的均方誤差(MSE),本文采用MSE描述改進(jìn)方法和傳統(tǒng)方法預(yù)測(cè)剩余壽命的精度。

圖5 剩余壽命預(yù)測(cè)的MSEFig.5 MSE of the predicted remaining useful life

由圖5可知，本文提出的模型方法(M2)預(yù)測(cè)誤差要小于傳統(tǒng)模型方法(M1)。

6 結(jié)束語(yǔ)

本文針對(duì)工程中隨機(jī)退化方法中失效閾值難以確定且無(wú)法兼顧失效數(shù)據(jù)和退化數(shù)據(jù)、半隨機(jī)濾波方法中初始分布選擇缺乏理論依據(jù)等問(wèn)題，在文獻(xiàn)[10]預(yù)測(cè)模型的基礎(chǔ)上對(duì)其進(jìn)行改進(jìn)，將描述線(xiàn)性隨機(jī)退化設(shè)備壽命的逆高斯分布作為半隨機(jī)濾波模型的初始分布，然后利用半隨機(jī)濾波技術(shù)預(yù)測(cè)設(shè)備的剩余壽命，并給出了模型參數(shù)的極大似然估計(jì)方法，最后，基于疲勞裂紋增長(zhǎng)數(shù)據(jù)對(duì)本文提出的方法進(jìn)行了驗(yàn)證。結(jié)果表明：相比于文獻(xiàn)[10]的半隨機(jī)預(yù)測(cè)模型，本文所提方法能夠更有效地提高剩余壽命預(yù)測(cè)精度。