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

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

0 引 言

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

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

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

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

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

1.1 基于維納過(guò)程的退化模型

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

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

(1)

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

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

(1)X(0)=0；

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

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

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

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

(2)

1.2 基于伽馬過(guò)程的退化模型

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

(1)G(0)=0；

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

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

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

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

(3)

1.3 相關(guān)性分析

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

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

(4)

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

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

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

(5)

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

AIC=-2log(A)+2m

(6)

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

2 剩余壽命預(yù)測(cè)

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

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

(7)

2.1 剩余壽命分布

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

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

(8)

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

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

(9)

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

(10)

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

(11)

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

(12)

由于通過(guò)式(12)直接計(jì)算fTG(t|ωG,α,β)比較困難，可以通過(guò)離散化的方式獲得裝備剩余壽命概率密度函數(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ù)估計(jì)

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

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

(15)

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

(16)

(17)

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

L(α,β)=

(18)

由極大似然估計(jì)法，令

(19)

(20)

3 仿真算例和實(shí)例分析

3.1 仿真算例分析

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

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

圖2 兩個(gè)模擬退化增量

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

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

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

圖3 各個(gè)監(jiān)測(cè)時(shí)刻點(diǎn)的參數(shù)估計(jì)值

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

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

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

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

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

3.2 實(shí)例分析

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

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

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

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

表3 不同監(jiān)測(cè)時(shí)刻點(diǎn)所對(duì)應(yīng)的參數(shù)估計(jì)值

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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