基于二元混合隨機(jī)過程的軸承剩余壽命預(yù)測 ①

金曉航 李建華 郭遠(yuǎn)晶 賈 虹*

(*特種裝備制造與先進(jìn)加工技術(shù)教育部重點實驗室(浙江工業(yè)大學(xué)) 杭州 310023)(**浙江工業(yè)大學(xué)機(jī)械工程學(xué)院 杭州 310023)(***寧?？h浙工大科學(xué)技術(shù)研究院 寧海 315600)(****浙江工業(yè)大學(xué)之江學(xué)院 紹興 312030)

0 引 言

裝備(如盾構(gòu)機(jī)、風(fēng)力發(fā)電機(jī)、水泵水輪發(fā)電機(jī)組等)在當(dāng)今社會建設(shè)和經(jīng)濟(jì)發(fā)展中發(fā)揮了重要作用，裝備一旦發(fā)生故障會造成重大事故和巨額損失[1-5]。為避免災(zāi)難事故的發(fā)生，同時降低運(yùn)維成本，需要對裝備健康狀態(tài)進(jìn)行有效的監(jiān)測、評估并準(zhǔn)確地預(yù)測其剩余壽命(remaining useful life， RUL)。目前用于剩余壽命預(yù)測的方法有貝葉斯濾波算法、神經(jīng)網(wǎng)絡(luò)算法、隨機(jī)過程模型等[6-7]?；陔S機(jī)過程的裝備剩余壽命預(yù)測方法能夠較好地用于分析裝備在工作過程中受到外圍環(huán)境、運(yùn)行條件、制造工藝等不同因素造成的不確定性影響，且具有良好的計算分析特性，成為近年來研究的熱點[8-9]。

為了更加全面地反映裝備的健康狀態(tài)和準(zhǔn)確地預(yù)測其剩余壽命，通常需要同時利用多個性能指標(biāo)。在利用多個性能指標(biāo)共同表征裝備的健康狀態(tài)時，需考慮不同性能指標(biāo)間的相關(guān)性問題[10]。如Bai等人[11]在對應(yīng)力-強(qiáng)度模型的可靠性分析中利用Gumbel Copula函數(shù)描述應(yīng)力和強(qiáng)度變量之間的相關(guān)關(guān)系，并通過工程案例進(jìn)行分析，證明了所建立模型的實用性。金曉航等人[12]利用Copula函數(shù)分析兩個性能指標(biāo)間的相關(guān)關(guān)系，提出了一種基于二元維納過程的軸承剩余壽命預(yù)測方法。上述的工作都是基于多個性能指標(biāo)服從同一種隨機(jī)過程的假設(shè)下開展的。但在一些情況下，所選取的性能指標(biāo)性質(zhì)差異較大時，不同指標(biāo)無法同時服從同一種隨機(jī)過程。

基于上述分析，考慮一種將維納過程和伽馬過程融合分析的二元混合隨機(jī)過程，用于預(yù)測裝備存在兩個性質(zhì)不同的性能指標(biāo)時的剩余壽命。首先基于維納過程和伽馬過程對不同性能指標(biāo)的退化過程分別進(jìn)行建模分析，然后利用赤池信息準(zhǔn)則(Akaike information criterion，AIC)選取合適的Copula函數(shù)分析不同性能指標(biāo)間的相關(guān)性，采用分步極大似然估計法在線更新模型參數(shù)，實現(xiàn)裝備的剩余壽命預(yù)測。

1 基于二元混合隨機(jī)過程的退化建模

目前應(yīng)用較為廣泛的隨機(jī)過程模型有維納過程和伽馬過程等[9]。維納過程適合分析退化過程連續(xù)且呈現(xiàn)隨機(jī)波動的性能指標(biāo)；伽馬過程適合分析退化過程嚴(yán)格單調(diào)遞增(或遞減)的性能指標(biāo)[13-14]。因此，需要結(jié)合裝備性能指標(biāo)變化的特點選擇合適的隨機(jī)過程建立相應(yīng)的退化模型。本文針對裝備兩個不同特點的性能指標(biāo)，分別利用維納過程和伽馬過程進(jìn)行建模分析。

1.1 基于維納過程的退化模型

若裝備退化過程的性能指標(biāo)表現(xiàn)出非嚴(yán)格單調(diào)遞增(或遞減)的過程，可以利用維納過程進(jìn)行建模分析。模型的表達(dá)式如式(1)所示[12]。

X(t)=X(0)+μt+σB(t)

(1)

式中，X(t)表示為t時刻的性能指標(biāo)值，μ為漂移系數(shù)，σ為擴(kuò)散系數(shù)，B(t)為標(biāo)準(zhǔn)維納過程。

利用維納過程分析的退化過程{X(t),t≥0}具有如下性質(zhì)。

(1)X(0)=0；

(2) {X(t),t≥0}具有平穩(wěn)獨(dú)立增量；

(3) 任意時刻間的增量服從正態(tài)分布，即：

X(t+Δt)-X(t)～N(μΔt,σ2Δt)

根據(jù)維納過程的定義，增量ΔX(t)服從正態(tài)分布，令Δt=1，其概率密度函數(shù)可表示為

(2)

1.2 基于伽馬過程的退化模型

若裝備退化過程的性能指標(biāo)表現(xiàn)為嚴(yán)格單調(diào)遞增(或遞減)的情況時，可以利用伽馬過程進(jìn)行建模分析[15]。設(shè){G(t),t≥0}是一個形狀參數(shù)為α，尺度參數(shù)為β的伽馬過程，利用伽馬過程分析的退化過程應(yīng)滿足以下性質(zhì)。

(1)G(0)=0；

(2) {G(t),t≥0}具有平穩(wěn)獨(dú)立增量；

(3) 任意時刻間的增量服從伽馬分布，即：

G(t+Δt)-G(t)～Ga(αΔt,β)

根據(jù)伽馬過程的定義，增量ΔG(t)服從伽馬分布Ga(αΔt,β)，令Δt=1，其概率密度函數(shù)為

(3)

1.3 相關(guān)性分析

利用多個性能指標(biāo)同時表征裝備的健康狀態(tài)時，需要考慮不同性能指標(biāo)之間的相關(guān)特性。Copula函數(shù)能夠有效地構(gòu)建出各個性能指標(biāo)的聯(lián)合分布函數(shù)和邊緣分布函數(shù)之間的關(guān)系，可被用于分析不同性能指標(biāo)之間的相關(guān)性。因此，引入Copula函數(shù)分析裝備多個性能指標(biāo)之間的相關(guān)關(guān)系，如下式所示。

F(x1,x2,…,xn)=C(u1,u2,…,un;θ)

(4)

式中，F(xiàn)(x1,x2,…,xn)為多個變量的聯(lián)合分布函數(shù)；C(u1,u2,…,un;θ)和θ是分別為Copula函數(shù)和Copula函數(shù)中的參數(shù)；u1,u2,…,un為各個變量的邊緣分布函數(shù)。

當(dāng)各變量的邊緣分布函數(shù)已知時，它們的聯(lián)合分布函數(shù)可根據(jù)Sklar定理[15]獲得。用Copula函數(shù)描述兩個退化性能指標(biāo)之間的相關(guān)特性時，聯(lián)合分布函數(shù)可表示為

F(x1,x2)=C(F(x1),F(x2);θ)

(5)

考慮到常用的Copula函數(shù)有Gaussian Copula、Frank Copula、Gumbel Copula和Clayton Copula等[16]，因此需要結(jié)合裝備實際退化過程選擇合適的Copula函數(shù)。鑒于赤池信息準(zhǔn)則(AIC)是一種廣泛被用于評價模型分析效果優(yōu)劣的方法[17]，文章中利用AIC信息準(zhǔn)則選擇合適Copula函數(shù)：

AIC=-2log(A)+2m

(6)

式中，A為模型對應(yīng)的極大似然函數(shù)，m為模型中參數(shù)個數(shù)，AIC的值越小說明所選擇的模型越合適。

2 剩余壽命預(yù)測

由于同時利用了裝備的兩個性能指標(biāo)(性能指標(biāo)1服從維納過程，性能指標(biāo)2服從伽馬過程)進(jìn)行剩余壽命預(yù)測，所以定義當(dāng){X1(t),t≥0}和{X2(t),t≥0}中的任意一個性能指標(biāo)超過其所對應(yīng)的失效閾值ω1、ω2時，即認(rèn)為裝備發(fā)生故障[18]。因此，當(dāng)利用兩個性能指標(biāo)預(yù)測分析裝備的剩余壽命時，剩余壽命T定義為

T=inf{t：X1(t)≥ω1或X2(t)≥ω2}

(7)

2.1 剩余壽命分布

裝備退化時，若其性能指標(biāo)變化的情況服從維納過程，則裝備的剩余壽命服從逆高斯分布[19]。確定維納過程模型中的參數(shù)μ、σ后，即可得到剩余壽命TW的概率密度函數(shù)：

fTW(t|ωW,μ,σ)=

(8)

若裝備的性能指標(biāo)變化服從伽馬過程時(其性能指標(biāo)初始值為0，失效閾值為ωG)，根據(jù)伽馬過程的性質(zhì)可知，裝備退化過程中其性能指標(biāo)到達(dá)失效閾值的首達(dá)時間TG為

FTG(t)=P{TG≤t}=P{G(t)≥ωG}

(9)

通過式(9)，TG的概率分布函數(shù)可以表示為

(10)

式中，Γ(a,b)為不完全伽馬函數(shù)，可表示為

(11)

裝備剩余壽命TG的概率密度函數(shù)為

(12)

由于通過式(12)直接計算fTG(t|ωG,α,β)比較困難，可以通過離散化的方式獲得裝備剩余壽命概率密度函數(shù)的近似表示。取tn=n/α,n=0,1,…，并令qn=P{tn

=αqn

(13)

基于不同性能指標(biāo)獲得裝備剩余壽命概率密度函數(shù)后，利用Copula函數(shù)分析性能指標(biāo)間的相關(guān)特性，即可獲得裝備剩余壽命T的聯(lián)合概率函數(shù)：

fT(t|ωW,ωG,μ,σ,α,β,θ)=c(FTW(t),FTG(t);θ)

·fTW(t|ωW,μ,σ)

·fTG(t|ωG,α,β)

(14)

式中c(·)為C(·)的密度函數(shù)。

2.2 參數(shù)估計

由式(14)可知剩余壽命預(yù)測模型中包含5個未知參數(shù)(μ,σ,α,β,θ)，可通過分步極大似然估計法求解各個參數(shù)：(1)更新不同性能指標(biāo)退化模型中的參數(shù)(μ,σ,α,β)；(2)更新Copula函數(shù)中的參數(shù)θ，具體過程如下。

第1步。由維納過程的性質(zhì)知道，服從維納過程的性能指標(biāo)增量ΔXi～N(μΔti,σ2Δti)，可得到模型參數(shù)(μ,σ)的似然函數(shù)為

(15)

基于式(15)，分別對漂移系數(shù)μ和擴(kuò)散系數(shù)σ求偏導(dǎo)，可得μ、σ的極大似然估計值分別為

(16)

(17)

由伽馬過程的性質(zhì)知道ΔGi～Ga(αΔti,β)，可獲得模型參數(shù)(α,β)的似然函數(shù)為

L(α,β)=

(18)

由極大似然估計法，令

(19)

(20)

3 仿真算例和實例分析

3.1 仿真算例分析

通過仿真算例來驗證分析本文所提出的基于二元混合隨機(jī)過程的剩余壽命預(yù)測方法的可行性。分別利用維納過程和伽馬過程的性質(zhì)(維納過程的增量服從正態(tài)分布，伽馬過程的增量服從伽馬分布)，在Matlab軟件中生成兩個性能指標(biāo)數(shù)據(jù)，失效閾值分別定義為(ω1,ω2)=(55.18, 31.89)，如圖1所示。兩個性能指標(biāo)的增量分別滿足正態(tài)分布N(0.15, 0.3)和伽馬分布Ga(0.2, 0.4)，如圖2所示。

圖1 兩模擬性能指標(biāo)退化趨勢

圖2 兩個模擬退化增量

利用AIC信息準(zhǔn)則選擇Copula函數(shù)，結(jié)果如表1所示。從表1中可以看到，針對該組數(shù)據(jù)，F(xiàn)rank Copula函數(shù)所對應(yīng)的AIC值最小，因此選擇Frank Copula函數(shù)分析性能指標(biāo)間的相關(guān)特性。利用式(8)和(13)分別計算基于不同性能指標(biāo)的裝備剩余壽命的邊緣概率密度函數(shù)。

表1 4種Copula函數(shù)的AIC值

根據(jù)Sklar定理構(gòu)建出剩余壽命的聯(lián)合概率密度函數(shù)，利用兩步極大似然估計法更新模型參數(shù)。首先，基于兩個性能指標(biāo)各個時刻的測量值，結(jié)合極大似然估計法在線更新退化模型中的參數(shù)(μ,σ,α,β)，進(jìn)而確定兩個不同性能指標(biāo)剩余壽命的邊緣分布函數(shù)。然后，將兩個邊緣分布函數(shù)作為Copula函數(shù)的輸入，再次利用極大似然估計法更新Copula函數(shù)中的參數(shù)θ，獲得各個時刻參數(shù)的估計結(jié)果，如圖3所示。

圖3 各個監(jiān)測時刻點的參數(shù)估計值

獲得模型參數(shù)的估計值后，結(jié)合式(14)計算得到裝備的剩余壽命概率密度函數(shù)，如圖4所示。觀察不同時刻的概率密度函數(shù)可以看出，隨著觀測數(shù)據(jù)的增多，剩余壽命預(yù)測結(jié)果的分布范圍逐漸縮小，說明其預(yù)測結(jié)果的精確度越來越高。將各個時刻概率密度函數(shù)最大值所對應(yīng)的時間作為剩余壽命的預(yù)測值，可得到整個過程的剩余壽命預(yù)測的結(jié)果，如圖5所示，從圖中可以看出剩余壽命的預(yù)測值逐漸逼近真實值。

圖4 不同時刻點剩余壽命的概率密度函數(shù)

利用均方根誤差(root mean square error, RMSE)對不同方法預(yù)測結(jié)果的誤差進(jìn)行分析，結(jié)果如表2所示。RMSE值越小說明預(yù)測效果越好[20]。從表2中可以看出，相較于基于一元隨機(jī)過程的剩余壽命預(yù)測方法，基于二元混合隨機(jī)過程的剩余壽命預(yù)測方法的預(yù)測效果更好。

圖5 剩余壽命預(yù)測結(jié)果

表2 不同預(yù)測方法預(yù)測結(jié)果的RMSE

3.2 實例分析

利用PRONOSTIA實驗平臺上的軸承全壽命周期的實驗數(shù)據(jù)對基于二元混合隨機(jī)過程的剩余壽命預(yù)測方法做進(jìn)一步的驗證分析[21]。實驗過程中，利用Dytran3035B的加速度傳感器采集振動信號，通過NI數(shù)據(jù)采集系統(tǒng)每隔10 s，以25.6 kHz的采樣頻率，采集并存儲一個時長為0.1 s的振動數(shù)據(jù)；利用Pt100的溫度傳感器連續(xù)采集溫度信號，采樣頻率為10 Hz。利用軸承振動信號的有效值和溫度信號的平均值作為表征其健康狀態(tài)的性能指標(biāo)。圖6(a)、(b)是軸承性能退化階段，兩個性能指標(biāo)歸一化后的變化過程，分別利用維納和伽馬過程分析振動和溫度的數(shù)據(jù)。對軸承兩個性能指標(biāo)建模分析之前，需要判斷它們是否符合隨機(jī)過程的性質(zhì)，通過Jarque-Bera驗證分析，振動信號有效值的增量服從正態(tài)分布N(0.0045, 0.0483)，如圖7(a)所示。由于受環(huán)境溫度以及實驗平臺其他部件傳熱等因素的影響，軸承溫度信號的實際增量并不一定符合伽馬分布，因此需要對溫度信號進(jìn)行預(yù)處理分析(即對時序的溫度信號進(jìn)行包絡(luò)分析，且令后一時刻的溫度值大于等于前一時刻的值)。結(jié)合伽馬過程的性質(zhì)，預(yù)處理后的溫度性能指標(biāo)如圖6(c)所示，經(jīng)檢驗，該信號的增量服從Ga(1.2645, 0.0015)的伽馬分布，如圖7(b)所示。圖7的結(jié)果表明針對軸承所構(gòu)建的兩個性能指標(biāo)可以分別利用維納和伽馬過程進(jìn)行建模分析。

圖6 軸承兩性能指標(biāo)退化趨勢

圖7 軸承兩個性能指標(biāo)增量的直方圖

通過AIC信息準(zhǔn)則選擇Copula函數(shù)分析不同性能指標(biāo)間的相關(guān)特性，如表1所示。從表1中可知Frank Copula函數(shù)所對應(yīng)的AIC值為-3 856，相較于其他3個Copula函數(shù)，F(xiàn)rank Copula函數(shù)的AIC值最小，因此選擇Frank Copula函數(shù)分析軸承兩個性能指標(biāo)間的相關(guān)特性。進(jìn)而構(gòu)建軸承剩余壽命的聯(lián)合概率密度函數(shù)，利用分步極大似然估計法更新模型中的參數(shù)，在軸承性能退化階段的第100、120、140、160、180時刻的模型參數(shù)估計值，如表3所示。

表3 不同監(jiān)測時刻點所對應(yīng)的參數(shù)估計值

獲得各個時刻模型參數(shù)的估計值后，結(jié)合式(14)即可獲得對應(yīng)時刻的剩余壽命概率密度函數(shù)，如圖8所示。圖中實線表示的是概率密度函數(shù)最大值對應(yīng)的時刻，即為所預(yù)測的剩余壽命；虛線表示的是當(dāng)前時刻剩余壽命的真實值，實線與虛線之間的距離代表預(yù)測結(jié)果的誤差。圖8中(a)到(d)依次為軸承退化階段第100、130、150、170時刻的剩余壽命概率密度函數(shù)圖，從圖中可以看出，隨著測量時刻點的增加，實線與虛線之間的距離不斷減小，說明預(yù)測誤差逐漸縮小、預(yù)測精度不斷提高。

圖8 不同時刻剩余壽命概率密度函數(shù)

將不同時刻剩余壽命概率密度函數(shù)的最大值對應(yīng)的時間提取出來，獲得如圖9所示的軸承剩余壽命預(yù)測結(jié)果。從圖中可以看出，相較于單獨(dú)利用維納過程或伽馬過程的剩余壽命預(yù)測方法，基于二元混合隨機(jī)過程的剩余壽命方法的預(yù)測效果更好。

為了客觀地衡量不同剩余壽命方法的預(yù)測效果，利用均方根誤差和誤差區(qū)間分布對預(yù)測結(jié)果的精度分別進(jìn)行了比較分析。不同方法預(yù)測結(jié)果的RMSE值如表2所示，從表中可以看到基于二元混合隨機(jī)過程的軸承剩余壽命預(yù)測方法的RMSE值為

圖9 剩余壽命預(yù)測結(jié)果

3.6535，小于其他兩種方法的RMSE值。同時將所有預(yù)測結(jié)果的誤差值從大到小劃分為[-300, 30]、[-100, 20]、[-50, 10]和[-10, 5]4個不同的區(qū)間，結(jié)果如表4所示。從表中可以看出基于二元混合隨機(jī)過程的剩余壽命預(yù)測方法在各個區(qū)間中所占的比例分別為100%、82.61%、54.35%、27.17%，比其他剩余壽命預(yù)測方法在各個不同誤差區(qū)間中所占的比例更高。綜合上述分析，說明基于二元混合隨機(jī)過程的軸承剩余壽命預(yù)測方法的預(yù)測精度較高。

表4 不同剩余壽命預(yù)測方法的預(yù)測結(jié)果誤差區(qū)間對比分析

4 結(jié) 論

針對裝備存在多個性能指標(biāo)且指標(biāo)隨時間變化的不同特點，提出利用維納和伽馬過程相融合的二元混合隨機(jī)過程用于預(yù)測裝備的剩余壽命。

利用AIC信息準(zhǔn)則選擇合適的Copula函數(shù)分析不同性能指標(biāo)間的相關(guān)特性，通過分步參數(shù)估計方法更新模型參數(shù)并預(yù)測裝備未來時刻的健康狀態(tài)，進(jìn)而實現(xiàn)剩余壽命的預(yù)測。

分別利用仿真和軸承實驗數(shù)據(jù)對所提方法進(jìn)行驗證分析，結(jié)果表明基于二元混合隨機(jī)過程的剩余壽命預(yù)測方法的預(yù)測效果比單獨(dú)利用維納過程或伽馬過程進(jìn)行建模分析的剩余壽命預(yù)測方法更好。