多元隨機(jī)過(guò)程的剩余壽命預(yù)測(cè)

馮媛媛，楊玉桃，朱向洪

(南昌大學(xué)理學(xué)院，江西 南昌 330031)

現(xiàn)代工業(yè)產(chǎn)品和設(shè)備可靠性分析及剩余壽命預(yù)測(cè)工作已經(jīng)貫穿系統(tǒng)設(shè)計(jì)、生產(chǎn)、使用和維護(hù)各個(gè)階段，對(duì)于系統(tǒng)正常使用和健康管理具有重要意義。如今工程系統(tǒng)可靠性通常較高，應(yīng)用傳統(tǒng)基于故障數(shù)據(jù)的方法存在一定的困難；而在系統(tǒng)使用過(guò)程中，通常會(huì)伴隨著一個(gè)或者多個(gè)性能指標(biāo)的退化，退化數(shù)據(jù)的獲取相對(duì)簡(jiǎn)單且成本較低，因此基于退化數(shù)據(jù)的方法研究越來(lái)越多。

現(xiàn)有的研究工作大多數(shù)考慮一維退化數(shù)據(jù)。目前常用的方法包括基于退化軌跡、隨機(jī)過(guò)程和非參數(shù)模型。王美男等[1]利用粒子濾波方法對(duì)設(shè)備退化模型的參數(shù)進(jìn)行更新,通過(guò)建立與傳統(tǒng)方法不同的粒子濾波量測(cè)方程,使依賴于更新后的退化模型參數(shù)的預(yù)測(cè)曲線符合已知退化趨勢(shì),得到更加準(zhǔn)確的壽命預(yù)測(cè)結(jié)果。除了退化軌跡模型外，有學(xué)者嘗試用非參數(shù)模型的方法來(lái)預(yù)測(cè)系統(tǒng)剩余壽命?？凳貜?qiáng)等[2]提出一種基于無(wú)監(jiān)督深度模型遷移的滾動(dòng)軸承剩余使用壽命預(yù)測(cè)方法。該方法首先對(duì)滾動(dòng)軸承數(shù)據(jù)提取均方根特征，引入時(shí)間序列分割算法對(duì)振動(dòng)信號(hào)進(jìn)行狀態(tài)信息標(biāo)記，并在神經(jīng)網(wǎng)絡(luò)訓(xùn)練得到狀態(tài)識(shí)別模型的基礎(chǔ)上，利用狀態(tài)概率估計(jì)法結(jié)合狀態(tài)識(shí)別模型建立滾動(dòng)軸承壽命預(yù)測(cè)模型。與退化軌跡模型和非參數(shù)模型相比，基于隨機(jī)過(guò)程的剩余壽命受到了很多學(xué)者的青睞。一個(gè)重要的原因就是隨機(jī)過(guò)程模型能夠很好地描述退化過(guò)程隨時(shí)間的不確定性。而在這個(gè)方面，目前的研究也大部分只考慮了一維退化情況。李月鋅等[3]對(duì)電池容量退化服從非線性維納過(guò)程建立狀態(tài)空間模型,并認(rèn)為參數(shù)是服從共軛分布的隨機(jī)變量,增加了模型不確定性使之更加符合鋰離子電池容量的退化過(guò)程。利用自助法和粒子濾波對(duì)每一時(shí)刻的參數(shù)及退化狀態(tài)進(jìn)行估計(jì)和更新,進(jìn)而預(yù)測(cè)電池的剩余壽命。黃亮等[4]考慮了發(fā)動(dòng)機(jī)性能退化非線性的特點(diǎn),并采用多階段Wiener過(guò)程建立發(fā)動(dòng)機(jī)的性能退化模型。根據(jù)發(fā)動(dòng)機(jī)的歷史性能監(jiān)測(cè)數(shù)據(jù),利用極大似然估計(jì)和貝葉斯方法對(duì)參數(shù)進(jìn)行評(píng)估更新，進(jìn)而得到個(gè)體發(fā)動(dòng)機(jī)剩余壽命的實(shí)時(shí)預(yù)測(cè)值。

以上研究雖從各個(gè)方面考慮復(fù)雜系統(tǒng)的實(shí)際剩余壽命情況，但其僅僅只關(guān)注一維性能指標(biāo)對(duì)于系統(tǒng)退化過(guò)程的影響，而這種潛在假設(shè)在某些情況下不現(xiàn)實(shí)。如為了不同的目的(即照明和顏色)，一個(gè)發(fā)光二極管系統(tǒng)可能會(huì)經(jīng)歷不止一個(gè)降解過(guò)程[5]；在銣電燈工作時(shí)，其性能退化指標(biāo)應(yīng)綜合考慮銣的消耗和燈的強(qiáng)度降低[6]。在上述情況下，若僅關(guān)注一元性能指標(biāo)對(duì)系統(tǒng)的影響，將對(duì)進(jìn)一步的系統(tǒng)剩余壽命預(yù)測(cè)造成巨大的偏差。從而為了更加全面地反映系統(tǒng)的退化狀態(tài)并準(zhǔn)確地預(yù)測(cè)其剩余壽命，往往需要同時(shí)利用多個(gè)性能指標(biāo)。經(jīng)以上現(xiàn)象啟發(fā)，本文將構(gòu)造一種新權(quán)重融合方法，綜合考慮多個(gè)非平穩(wěn)伽馬過(guò)程，并采用馬爾科夫蒙特卡洛方法在線更新模型參數(shù)，得出系統(tǒng)剩余壽命概率密度函數(shù)。另外通過(guò)裂紋實(shí)驗(yàn)驗(yàn)證所提模型的有效性并得出部分結(jié)論。

1 模型假設(shè)

基于系統(tǒng)性能退化數(shù)據(jù)，構(gòu)建合理的退化過(guò)程模型是進(jìn)行系統(tǒng)剩余壽命預(yù)測(cè)的基礎(chǔ)。在描述性能指標(biāo)退化過(guò)程中，由于隨機(jī)過(guò)程描述變量在時(shí)間軸上不確定性的優(yōu)點(diǎn)和Gamma過(guò)程連續(xù)非減的良好性質(zhì)，所以本文選用非平穩(wěn)伽馬過(guò)程描述性能指標(biāo)退化過(guò)程。

現(xiàn)假定存在m個(gè)系統(tǒng)且收集了每個(gè)系統(tǒng)中l(wèi)個(gè)性能指標(biāo)在正常工作時(shí)間ni內(nèi)的退化數(shù)據(jù)，其中Xik(tj) 表示第i系統(tǒng)第k個(gè)性能指標(biāo)在tj時(shí)刻的退化量。

假設(shè)1第i個(gè)系統(tǒng)在退化過(guò)程中，第k個(gè)性能指標(biāo)Xik服從非平穩(wěn)伽馬過(guò)程，既而其性能指標(biāo)隨著時(shí)間的變化，其增量可表示為：

Xik(tj+c)-Xik(tj)～Ga(ηk(tj+c;θk),βk)

其中，ηk(t;θk)是關(guān)于時(shí)間的非減函數(shù)的形狀參數(shù)，且Δηk(tj+c;θk)=ηk(tj+c;θk)-ηk(tj;θk)；βk是服從伽馬過(guò)程Ga(δk,γk)的尺度參數(shù)，用來(lái)描述系統(tǒng)中指標(biāo)間的差異。

假設(shè)2若第i個(gè)系統(tǒng)中，任意一個(gè)性能指標(biāo)Xik超過(guò)了其對(duì)應(yīng)的失效閾值Dik，則認(rèn)定該系統(tǒng)失效。

假設(shè)3規(guī)定第i個(gè)系統(tǒng)的失效閾值Di由其各個(gè)性能指標(biāo)所對(duì)應(yīng)的失效閾值Dik綜合表示，即Di=wi1Di1+…+wilDil。其中Dik、wik(k=1,2,…,l)分別表示第i個(gè)系統(tǒng)中第k個(gè)性能指標(biāo)對(duì)應(yīng)的失效閾值和權(quán)重。

假設(shè)4在采集系統(tǒng)的各個(gè)性能指標(biāo)時(shí)，同個(gè)系統(tǒng)的數(shù)據(jù)采集時(shí)間必須相同，而對(duì)于不同系統(tǒng)而言，數(shù)據(jù)采集時(shí)間則不限定。

2 一元性能指標(biāo)剩余壽命

在工程數(shù)據(jù)中,通過(guò)多個(gè)性能指標(biāo)退化數(shù)據(jù)共同描述一個(gè)復(fù)雜系統(tǒng)的實(shí)際退化過(guò)程,并利用首達(dá)的概念對(duì)系統(tǒng)進(jìn)行剩余壽命預(yù)測(cè)。因此，有必要在進(jìn)行系統(tǒng)剩余壽命預(yù)測(cè)之前對(duì)一元性能指標(biāo)的剩余壽命進(jìn)行研究。

假設(shè)th+c表示Xk(th+c)退化過(guò)程首次達(dá)到失效閾值Dk的時(shí)間，另在本節(jié)中規(guī)定

xk=Xk(th+c)-Xk(th)～Ga(Δηk(th+c;θk),βk),從而第k個(gè)性能指標(biāo)的剩余壽命可重新定義為

根據(jù)伽馬過(guò)程的性質(zhì)可知，產(chǎn)品退化過(guò)程中其各個(gè)性能指標(biāo)到達(dá)失效閾值的首達(dá)時(shí)間th+c的概率分布函數(shù)：

(1)

為了數(shù)值求解方便，Padget等[7]提出了利用疲勞損傷分布(BS)近似th+c的CDF?，F(xiàn)定義Zk(t)=ηk(t,θk)且Yk(Zk)=Xk(t)，從而Yk(Zk)服平穩(wěn)伽馬過(guò)程且其失效時(shí)間分布可近似為：

(2)

(3)

在給定X1:h的情況下，概率密度函數(shù)中尺度參數(shù)βk的更新過(guò)程是f(βk|X1:h)。利用全概率公式可得：

3 權(quán)重

系統(tǒng)的實(shí)際退化情況是綜合考慮了各個(gè)性能指標(biāo)的結(jié)果，所以各個(gè)性能指標(biāo)間的關(guān)系研究就顯得十分必要。而不同性能指標(biāo)在描述系統(tǒng)實(shí)際退化情況中的占比不同，就需要權(quán)重來(lái)衡量每個(gè)性能指標(biāo)對(duì)系統(tǒng)的重要程度。本文在兩兩一致性指標(biāo)概念基礎(chǔ)上提出了一種新權(quán)重計(jì)算方法(即第i個(gè)系統(tǒng)的第p個(gè)性能指標(biāo)的權(quán)重wip)，其構(gòu)造過(guò)程如下。

1) 隸屬度函數(shù)。

(4)

2) 構(gòu)造權(quán)重wip。

本小節(jié)將在兩兩一致性指標(biāo)的基礎(chǔ)上構(gòu)建權(quán)重wip。首先利用隸屬度函數(shù)計(jì)算得出每個(gè)性能指標(biāo)對(duì)應(yīng)退化量，并利用文獻(xiàn)[8]中計(jì)算方法求出第i個(gè)系統(tǒng)中第p維性能指標(biāo)分別與其他維度的性能指標(biāo)間的一致性指標(biāo)。其次分別將一致性指標(biāo)除以其對(duì)應(yīng)維度性能指標(biāo)退化數(shù)據(jù)的平均值，可得第p維性能指標(biāo)和其他維性能指標(biāo)的一致性程度。最后用1減去第p維的累加一致性程度,即兩兩不一致性指標(biāo)作為權(quán)重wip。其構(gòu)造方程如下：

其中，l表示性能指標(biāo)的總維數(shù)，Sμ(μip,μiq)=sup min{μip(xip(tj)),μiq(xiq(tj))}是第i個(gè)系統(tǒng)中第p個(gè)和第q個(gè)性能變量的一致性指標(biāo)[8]。之后對(duì)rik進(jìn)行歸一化處理:

(5)

4 產(chǎn)品剩余壽命

在假設(shè)3的基礎(chǔ)上，利用權(quán)重計(jì)算方法融合多維性能指標(biāo)可得出目標(biāo)產(chǎn)品退化過(guò)程到達(dá)失效閾值D的首達(dá)時(shí)間th+ζ分布函數(shù)：

(6)

從式(6)卷積計(jì)算可以看出，融合多個(gè)性能指標(biāo)作為系統(tǒng)總體退化指標(biāo)再計(jì)算產(chǎn)品剩余壽命，由于公式的復(fù)雜性，計(jì)算的難度和復(fù)雜程度將非常大。進(jìn)而本文將簡(jiǎn)化產(chǎn)品剩余壽命計(jì)算：綜合式(3)所得的各個(gè)性能指標(biāo)剩余壽命密度函數(shù)作為系統(tǒng)的剩余壽命密度函數(shù)，即：

(7)

在給定X1:h的情況下，概率密度函數(shù)中的尺度參數(shù)βk的更新過(guò)程是f(βk|X1:h)。利用全概率公式可得：

5 參數(shù)評(píng)估

由于已有的性能指標(biāo)數(shù)據(jù)能對(duì)參數(shù)產(chǎn)生影響,所以在確定產(chǎn)品剩余壽命之前需利用已有數(shù)據(jù)更新參數(shù)?，F(xiàn)假設(shè)目標(biāo)系統(tǒng)th時(shí)仍可正常工作且可測(cè)得在th前的所有性能指標(biāo),令Xk,1:h=(xk(t1)…xk(th))'。表示目標(biāo)系統(tǒng)的第k個(gè)性能指標(biāo)在th之前的退化量。在給定目標(biāo)系統(tǒng)X1:h=(x1,1:h,x2,1:h,…,xl,1:h)'的情況下，根據(jù)貝葉斯理論可得在th時(shí)的參數(shù)后分布。

(8)

其中ΔXk(tj)=Xk(tj)-Xk(tj-1)，同樣地，在測(cè)得目標(biāo)系統(tǒng)在th+1時(shí)的性能退化數(shù)據(jù)后，由式(8)可以實(shí)現(xiàn)實(shí)時(shí)的更新尺度參數(shù)(β1,…,βl)。

本節(jié)將采用兩階段方式對(duì)參數(shù)進(jìn)行估計(jì)，步驟如下：

第1階段估計(jì)參數(shù)Θ1=(θ1,…,θl),Θ2=(β11,…,βm1),Θ3=(β12,…,βm2),…,Θl+1=(β1l,…,βml)。

對(duì)于第i個(gè)系統(tǒng)，(θ1,βi1,…,θl,βil)完全似然函數(shù)可表示為：

(9)

其中ψ(.)表示雙伽馬函數(shù)，為Γ(.)算法的導(dǎo)函數(shù)。

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

以裂紋增長(zhǎng)數(shù)據(jù)為例，對(duì)上述模型和方法進(jìn)行有效性驗(yàn)證。假設(shè)系統(tǒng)存在兩個(gè)裂紋增長(zhǎng)點(diǎn)，為了便于分析，將裂紋長(zhǎng)度作為系統(tǒng)性能指標(biāo)(P1和P2)。退化數(shù)據(jù)來(lái)源于文獻(xiàn)[9]，其中包括21個(gè)裂紋增長(zhǎng)實(shí)驗(yàn)數(shù)據(jù)序列，測(cè)量間隔為0.01百萬(wàn)轉(zhuǎn)，測(cè)量終止時(shí)間為0.09百萬(wàn)轉(zhuǎn)。

為便于驗(yàn)證分析，選取20個(gè)數(shù)據(jù)序列作為實(shí)驗(yàn)系統(tǒng)退化數(shù)據(jù)，并將其平均分為兩組，分別作為兩個(gè)性能指標(biāo)的退化數(shù)據(jù)。此外，兩個(gè)性能指標(biāo)的初始退化量為2.286 cm，相應(yīng)的故障閾值設(shè)為(P1=4.064,P2=3.302)。其具體的退化情況如圖1所示。

時(shí)間/百分次

表1 參數(shù)估計(jì)

為了驗(yàn)證模型的有效性，模擬了基于以上的參數(shù)估計(jì)的目標(biāo)系統(tǒng)的一元預(yù)測(cè)和二元預(yù)測(cè)，從而得出所提多元預(yù)測(cè)方法優(yōu)于一元預(yù)測(cè)的結(jié)論。現(xiàn)已知目標(biāo)系統(tǒng)的性能指標(biāo)退化情況如圖2。

在得到目標(biāo)系統(tǒng)的退化數(shù)據(jù)后，根據(jù)上文所提出的權(quán)重計(jì)算方法分別算出P1和P2得出表2。

表2 權(quán)重及調(diào)節(jié)參數(shù)

從圖2中可知當(dāng)P1超過(guò)對(duì)應(yīng)的失效閾值1.6時(shí)，對(duì)應(yīng)的時(shí)間為0.097 1百萬(wàn)轉(zhuǎn)(時(shí)間間隔為0.01百萬(wàn)轉(zhuǎn))，在本文中用t1:9={0.01,0.02,…,0.09}來(lái)驗(yàn)證模型的正確性。

本文利用最大誤差Δ=(|xreal-xpre|)評(píng)判標(biāo)準(zhǔn)，用以比較兩種預(yù)測(cè)方式的優(yōu)劣性(見(jiàn)圖3)。由目標(biāo)系統(tǒng)數(shù)據(jù)計(jì)算可得一元預(yù)測(cè)的最大誤差分別為P1=0.140 1、P2=0.262 6。二元綜合預(yù)測(cè)的最大誤差為0.061 5。從而可以得出本文所提出的綜合多元性能指標(biāo)的方法更具優(yōu)勢(shì)。

時(shí)間/百分次

(a) β1在0.08百萬(wàn)轉(zhuǎn)時(shí)的采樣情況

應(yīng)用多元變量模型建立其目標(biāo)系統(tǒng)的剩余壽命方程，并得出了系統(tǒng)在9個(gè)時(shí)間點(diǎn)的RL估計(jì)值，如圖5所示，其中虛線表示實(shí)際系統(tǒng)剩余壽命，實(shí)線表示模型評(píng)估的剩余壽命。

圖5 剩余壽命的概率密度函數(shù)

從圖(6)可知隨著時(shí)間的推移，剩余壽命的相對(duì)誤差越小(即使用的退化數(shù)據(jù)越多，相對(duì)誤差就越小)。因此為了獲得準(zhǔn)確的估計(jì)，有必要在獲得新的退化數(shù)據(jù)后重新估計(jì)結(jié)果。

7 結(jié)論

1) 相較于一元退化變量進(jìn)行剩余壽命預(yù)測(cè)，多變量預(yù)測(cè)方法能夠更為有效地估計(jì)系統(tǒng)剩余壽命；

2) 能夠準(zhǔn)確地在線預(yù)測(cè)系統(tǒng)的剩余壽命，并且隨著退化數(shù)據(jù)累積，系統(tǒng)剩余壽命估計(jì)更為精確，降低維修成本提高維修效率。