基于性能退化數(shù)據(jù)融合的航空發(fā)動(dòng)機(jī)剩余壽命預(yù)測(cè)方法*

劉君強(qiáng) 謝吉偉

(南京航空航天大學(xué)民航學(xué)院 南京 211106)

0 引 言

航空發(fā)動(dòng)機(jī)的健康狀態(tài)大多數(shù)會(huì)反映在性能參數(shù)上，當(dāng)其性能參數(shù)超過(guò)特定閾值時(shí)，發(fā)動(dòng)機(jī)就會(huì)下發(fā)或者進(jìn)行返廠維修.所以，在航空發(fā)動(dòng)機(jī)的使用過(guò)程中如何利用性能參數(shù)對(duì)其進(jìn)行可靠性評(píng)估與剩余壽命預(yù)測(cè)成為亟待解決的問(wèn)題.對(duì)此，國(guó)內(nèi)外學(xué)者做了相關(guān)的研究.王燕萍等[1]在多種先驗(yàn)分布下，基于MCMC方法，利用Bayes方法將失效、截尾時(shí)間數(shù)據(jù)進(jìn)行融合并對(duì)發(fā)動(dòng)機(jī)可靠性進(jìn)行評(píng)估；孫見(jiàn)忠等[2]提出了一種基于動(dòng)態(tài)線性模型的航空發(fā)動(dòng)機(jī)在翼壽命預(yù)測(cè)方法；王華偉等[3]基于Gamma過(guò)程建立了航空發(fā)動(dòng)機(jī)性能可靠性評(píng)估模型；朱磊等[4]提出了基于Wiener過(guò)程的航空發(fā)動(dòng)機(jī)性能可靠性評(píng)估與剩余壽命預(yù)測(cè)方法.由于航空發(fā)動(dòng)機(jī)在復(fù)雜環(huán)境下工作，其性能參數(shù)變化具有較強(qiáng)的隨機(jī)性，采用隨機(jī)過(guò)程模型對(duì)航空發(fā)動(dòng)機(jī)的性能參數(shù)進(jìn)行建模是最合適的.在隨機(jī)過(guò)程模型中，Wiener過(guò)程模型由于其自身的計(jì)算優(yōu)勢(shì)及良好特性，是目前基于性能退化可靠性分析中應(yīng)用最為廣泛的模型.

雖然，Wiener過(guò)程模型在可靠性建模中有著廣泛的應(yīng)用，但仍然存在如下幾方面的問(wèn)題：①當(dāng)性能退化數(shù)據(jù)較少時(shí)，預(yù)測(cè)精度無(wú)法保證；②模型計(jì)算結(jié)果只能反映總體特征，無(wú)法反映個(gè)體的差異性.

針對(duì)第一個(gè)問(wèn)題，劉琦等[5]利用Bayes方法融合火箭發(fā)動(dòng)機(jī)實(shí)時(shí)運(yùn)行數(shù)據(jù)和地面試驗(yàn)數(shù)據(jù)對(duì)其進(jìn)行了可靠性分析；彭寶華等[6-7]利用Bayes方法融合某產(chǎn)品的歷史失效時(shí)間數(shù)據(jù)與實(shí)時(shí)性能退化數(shù)據(jù)對(duì)其進(jìn)行壽命預(yù)測(cè)；王小林等[8-9]提出了一種融合歷史性能退化數(shù)據(jù)、歷史壽命數(shù)據(jù)與實(shí)時(shí)性能退化數(shù)據(jù)的可靠性評(píng)估方法；Yan等[10-11]提出了利用經(jīng)驗(yàn)Bayes融合性能退化數(shù)據(jù)與壽命數(shù)據(jù)的剩余壽命預(yù)測(cè)方法.但由于Bayes統(tǒng)計(jì)分析涉及到大量的模擬分析、數(shù)值積分等問(wèn)題，通常需要借助MCMC方法進(jìn)行計(jì)算.

針對(duì)第二個(gè)問(wèn)題，Wang等[12]將Wiener過(guò)程的漂移參數(shù)和擴(kuò)散參數(shù)都看作隨機(jī)變量，假設(shè)參數(shù)先驗(yàn)分布為正態(tài)-逆伽馬分布，并利用期望最大化(EM)算法確定參數(shù)先驗(yàn)分布中的未知超參數(shù)；Peng等[13]在對(duì)金屬化膜電容器進(jìn)行剩余壽命預(yù)測(cè)研究時(shí)，也同樣采用了隨機(jī)化參數(shù)的Wiener過(guò)程進(jìn)行建模，并利用EM算法的對(duì)參數(shù)先驗(yàn)分布進(jìn)行估計(jì)，Gebraeel等[14]則采用先驗(yàn)矩的方法，對(duì)隨機(jī)參數(shù)的退化模型參數(shù)先驗(yàn)分布進(jìn)行求解.

結(jié)合上述研究的優(yōu)缺點(diǎn)，本文首先采用隨機(jī)參數(shù)的Wiener過(guò)程對(duì)航空發(fā)動(dòng)機(jī)性能退化進(jìn)行建模，選取參數(shù)的先驗(yàn)分布為Jeffrey無(wú)信息先驗(yàn)分布，利用Bayes方法融合航空發(fā)動(dòng)機(jī)的現(xiàn)場(chǎng)性能退化數(shù)據(jù)與截尾時(shí)間數(shù)據(jù)，當(dāng)性能退化數(shù)據(jù)較少時(shí)，借助bootstrap自助法得到性能退化數(shù)據(jù)的自助樣本，然后根據(jù)MCMC方法對(duì)模型參數(shù)的后驗(yàn)分布進(jìn)行抽樣估計(jì)，并使用Kolmogorov-Smirnov檢驗(yàn)確定模型參數(shù)的最優(yōu)擬合后驗(yàn)分布.最后，由模型參數(shù)估計(jì)值推導(dǎo)出航空發(fā)動(dòng)機(jī)的剩余壽命分布，并對(duì)其剩余壽命進(jìn)行預(yù)測(cè).

1 航空發(fā)動(dòng)剩余壽命預(yù)測(cè)模型

1.1 Wiener過(guò)程

如果連續(xù)時(shí)間隨機(jī)過(guò)程{X(t),t>0}服從漂移系數(shù)為μ和擴(kuò)散系數(shù)為σ的Wiener過(guò)程，那么X(t)滿足以下條件：

1)X(0)=0并且X(t)在t=0處連續(xù).

2) 對(duì)于兩個(gè)不相交的時(shí)間區(qū)間[t1,t2]和[t3,t4]，且滿足t1

3) 時(shí)刻t到t+Δt之間的增量X(t+Δt)-X(t)服從均值為μΔt和方差為σ2Δt的正態(tài)分布.

1.2 剩余壽命分布

假設(shè)航空發(fā)動(dòng)機(jī)的性能退化趨勢(shì)服從Wiener過(guò)程且失效閾值D已知.根據(jù)失效物理分析，當(dāng)航空發(fā)動(dòng)機(jī)性能退化量X(t)首次超過(guò)預(yù)定閾值D時(shí)，航空發(fā)動(dòng)機(jī)會(huì)下發(fā)進(jìn)行維修.設(shè)T表示航空發(fā)動(dòng)機(jī)的在翼壽命，則在翼壽命可以表示為T=inf{x:X(t)≥D;t>0}.由性能退化導(dǎo)致的在翼壽命T服從逆高斯分布，其概率密度函數(shù)形式與累積概率分布函數(shù)為

(1)

(2)

式中：Φ(x)為標(biāo)準(zhǔn)正態(tài)累積分布函數(shù).

如果某臺(tái)航空發(fā)動(dòng)機(jī)運(yùn)行到τ時(shí)刻累積性能退化量仍未超過(guò)預(yù)定閾值D，記其當(dāng)前性能退化量為xτ(xτ≤D)，根據(jù)Wiener過(guò)程的齊次馬爾可夫性質(zhì)，航空發(fā)動(dòng)機(jī)的剩余在翼壽命Tτ為

Tτ=inf{t|X(t)≥D-xτ,t>0}

(3)

由式(3)可知，剩余在翼壽命Tτ也服從逆高斯分布，其概率密度函數(shù)形式為

(4)

當(dāng)已知性能退化模型參數(shù)后，單臺(tái)航空發(fā)動(dòng)機(jī)的平均剩余在翼壽命為期望值為

(5)

2 性能退化建模與參數(shù)估計(jì)

2.1 性能退化數(shù)據(jù)預(yù)處理

對(duì)于一臺(tái)新投入使用的航空發(fā)動(dòng)機(jī)，制造商可以通過(guò)比較實(shí)際EGTM測(cè)量值與設(shè)計(jì)的EGTM閾值判斷航空發(fā)動(dòng)機(jī)的使用情況.隨著航空發(fā)動(dòng)機(jī)性能逐漸衰退，EGTM逐漸減少，通過(guò)EGTM可以衡量航空發(fā)動(dòng)機(jī)的性能衰退程度，當(dāng)EGTM到達(dá)規(guī)定閾值時(shí)，航空發(fā)動(dòng)機(jī)就會(huì)下發(fā)或者修.由于這種方法簡(jiǎn)單，容易操作，在航空公司中得到了廣泛的應(yīng)用.設(shè)yt為EGTM在t時(shí)刻的測(cè)量值，xt為由性能衰退導(dǎo)致的EGTM退化值，其中：

yt=xt+vt

(6)

式中：vt為觀測(cè)噪聲，包含由模型誤差、傳感器噪聲以及其他影響因素造成的誤差.

由于實(shí)際測(cè)量得到的EGTM數(shù)據(jù)含有較高的噪聲，直接使用EGTM測(cè)量數(shù)據(jù)進(jìn)行建模會(huì)引進(jìn)過(guò)多的不確定性使得剩余壽命預(yù)測(cè)誤差過(guò)大.因此，本文采用移動(dòng)平均平滑方法通過(guò)yt對(duì)進(jìn)行xt估計(jì).其中，移動(dòng)平均平滑法為

(7)

2.2 基于性能退化數(shù)據(jù)建模

假設(shè)獲得了某臺(tái)航空發(fā)動(dòng)機(jī)的性能退化樣本{xt}，設(shè)初始時(shí)刻的性能退化量為x0=0，分別在時(shí)刻t1,t2,…,tm得到其在當(dāng)前時(shí)刻的性能退化量為x1,x2,…,xm，記性能退化增量為Δxi(1

Δxi～N(μΔti,σ2Δti)

L(x,x*;μ,σ2)=

(8)

2.3 基于壽命數(shù)據(jù)建模

L(T;μ,σ)=

(9)

2.4 基于Bayes的參數(shù)估計(jì)

由于考慮到航空發(fā)動(dòng)機(jī)個(gè)體的差異性，認(rèn)為發(fā)動(dòng)機(jī)性能退化模型中的漂移系數(shù)和擴(kuò)散系數(shù)是隨機(jī)變量，且服從同一分布類型.因此，文中采用Bayes方法對(duì)漂移系數(shù)和擴(kuò)散系數(shù)的后驗(yàn)分布進(jìn)行計(jì)算.在已有的研究基礎(chǔ)中，大多數(shù)都采用正態(tài)-逆伽馬分布作為參數(shù)的先驗(yàn)分布，但事前假定參數(shù)先驗(yàn)分布的類型會(huì)影響最終剩余壽命預(yù)測(cè)結(jié)果.為避免主觀經(jīng)驗(yàn)造成的誤差，文中根據(jù)Jeffrey無(wú)信息先驗(yàn)原則，利用參數(shù)Fisher信息量的平方根作為參數(shù)的無(wú)信息先驗(yàn)分布.

π(μ)=[I(μ)]1/2∝1

由于參數(shù)(μ,σ2)無(wú)信息先驗(yàn)相互獨(dú)立，可以得到參數(shù)聯(lián)合無(wú)信息先驗(yàn)分布為

(10)

根據(jù)Bayes公式可得，參數(shù)(μ,σ2)的聯(lián)合后驗(yàn)分布為

π(μ,σ2|x,x*,T)∝π(μ,σ2)·L(x,x*,T;μ,σ2)∝

(11)

考慮一種特殊情況：當(dāng)無(wú)壽命數(shù)據(jù)時(shí)，即不存在失效時(shí)間數(shù)據(jù)與截尾時(shí)間數(shù)據(jù)時(shí)，漂移系數(shù)和擴(kuò)散系數(shù)的聯(lián)合后驗(yàn)分布為

π(μ,σ2|x,x*)∝π(μ,σ2)·L(x,x*;μ,σ)∝

(12)

式中：

由式(12)可知，后驗(yàn)分布與正態(tài)-逆伽馬分布的核形式相同，所以參數(shù)后驗(yàn)分布服從正態(tài)-逆伽馬分布，即

(13)

根據(jù)文獻(xiàn)[10]，漂移系數(shù)的邊緣后驗(yàn)分布為非中心t分布，擴(kuò)散系數(shù)的邊緣后驗(yàn)分布為逆伽馬分布，根據(jù)平方損失函數(shù)最小原則，參數(shù)Bayes估計(jì)值分別為

(14)

(15)

當(dāng)樣本數(shù)據(jù)中存在截尾時(shí)間數(shù)據(jù)時(shí)，由于參數(shù)后驗(yàn)分布形式較為復(fù)雜，難以給出參數(shù)后驗(yàn)分布的具體表達(dá)形式.故采用Gibbs抽樣方法產(chǎn)生后驗(yàn)分布的隨機(jī)樣本，從而對(duì)后驗(yàn)分布進(jìn)行統(tǒng)計(jì)推斷.Gibbs抽樣方法是從某個(gè)初始點(diǎn)出發(fā)通過(guò)滿條件分布的循環(huán)抽樣產(chǎn)生馬氏鏈.其中，參數(shù)μ的滿條件分布為

(16)

參數(shù)σ2的滿條件分布為

(17)

由于滿條件分布不容直接抽取，在抽樣過(guò)程中借助Metropolis-Hasting算法進(jìn)行抽樣，本文采用正態(tài)分布作為抽樣建議分布.

在滿足一定的迭代次數(shù)后，根據(jù)由Gibbs抽樣方法所產(chǎn)生的Monte Carlo樣本，對(duì)退化模型參數(shù)的后驗(yàn)分布進(jìn)行統(tǒng)計(jì)推斷.假設(shè)模型參數(shù)的邊緣后驗(yàn)分布形式已知，根據(jù)平方損失函數(shù)最小原則，參數(shù)Bayes估計(jì)值分別為

(18)

(19)

3 實(shí)例計(jì)算

選用某航空公司的CFM56-3型航空發(fā)動(dòng)機(jī)進(jìn)行實(shí)例計(jì)算.由于CFM56-3型發(fā)動(dòng)機(jī)屬于視情維修使用，沒(méi)有固定的壽命限制，該型發(fā)動(dòng)機(jī)第一次返廠維修時(shí)間都由用戶自行決定.在實(shí)際運(yùn)行過(guò)程中，航空公司一般都會(huì)在發(fā)動(dòng)機(jī)未超過(guò)預(yù)定閾值時(shí)便送廠維修，在實(shí)際情況下難以獲得該型號(hào)發(fā)動(dòng)機(jī)的真實(shí)失效時(shí)間數(shù)據(jù)[15]，只能夠獲得其對(duì)應(yīng)的截尾時(shí)間數(shù)據(jù).文中選取排氣溫度裕度(EGTM)作為該型發(fā)動(dòng)機(jī)的性能退化特征量，已知某臺(tái)航空發(fā)動(dòng)機(jī)的實(shí)際EGTM數(shù)據(jù)，見(jiàn)圖1.

圖1 EGTM測(cè)量數(shù)據(jù)與平滑數(shù)據(jù)

首先利用對(duì)EGTM數(shù)據(jù)進(jìn)行平滑處理，然后對(duì)過(guò)平滑后的EGTM數(shù)據(jù)進(jìn)行正態(tài)分布擬合優(yōu)度檢驗(yàn)，檢驗(yàn)結(jié)果P值為0.80，表明性能退化數(shù)據(jù)近似服從正態(tài)分布，因此，可以使用Wiener過(guò)程對(duì)該型發(fā)動(dòng)機(jī)性能退化進(jìn)行建模.

選取其中第2 500至第5 000個(gè)值作為模型參數(shù)的抽樣樣本，使用Kolmogorov-Smirnov檢驗(yàn)對(duì)參數(shù)可能的分布類型進(jìn)行擬合優(yōu)度檢驗(yàn)，最終確定漂移系數(shù)最優(yōu)擬合于正態(tài)分布，擴(kuò)散參數(shù)最優(yōu)擬合于對(duì)數(shù)正態(tài)分布.最終通過(guò)極大似然估計(jì)可以確定模型參數(shù)的后驗(yàn)分布形式，表1為在綜合考慮性能退化數(shù)據(jù)，截尾時(shí)間數(shù)據(jù)的情況下的參數(shù)點(diǎn)估計(jì)值及其對(duì)應(yīng)的95%置信區(qū)間：

表1 參數(shù)點(diǎn)估計(jì)值與95%置信區(qū)間

當(dāng)退化模型參數(shù)已知后，可以確定性能退化數(shù)據(jù)樣本的95%置信上下限，見(jiàn)圖2.

圖2 性能退化數(shù)據(jù)的95%置信區(qū)間

由圖2可知，MCMC方法得到的性能退化數(shù)據(jù)95%置信區(qū)間長(zhǎng)度要比利用極大似然估計(jì)方法得到的95%置信區(qū)間長(zhǎng)度小，這說(shuō)明了利用MCMC方法綜合考慮性能退化數(shù)據(jù)與失效時(shí)間數(shù)據(jù)所建立的性能退化模型可信度較高.

下面進(jìn)一步分析在只考慮性能退化數(shù)據(jù)、截尾時(shí)間數(shù)據(jù)和融合兩類數(shù)據(jù)這3種不同數(shù)據(jù)源情況下對(duì)剩余壽命估計(jì)的影響.航空發(fā)動(dòng)機(jī)的剩余壽命概率密度曲線與可靠性曲線見(jiàn)圖3.

圖3 概率密度和可靠性曲線

由圖3可知，只考慮性能退化數(shù)據(jù)得到的剩余壽命估計(jì)結(jié)果較小，只考慮失效時(shí)間數(shù)據(jù)得到的剩余壽命估計(jì)結(jié)果較大，而融合兩類數(shù)據(jù)得到的剩余壽命估計(jì)結(jié)果介于單獨(dú)利用一類數(shù)據(jù)得出的結(jié)果之間.為了更加直觀地反應(yīng)不同數(shù)據(jù)源對(duì)剩余壽命估計(jì)結(jié)果的影響，表2為利用不同數(shù)據(jù)源得到的平均剩余壽命預(yù)測(cè)值及剩余壽命90%置信下限.

D為僅考慮性能退化數(shù)據(jù)，T為僅考慮壽命數(shù)據(jù)，TD為融合性能退化數(shù)據(jù)與壽命數(shù)據(jù).由表2可知，融合性能退化數(shù)據(jù)與壽命數(shù)據(jù)的方法能夠較準(zhǔn)確地預(yù)測(cè)航空發(fā)動(dòng)機(jī)個(gè)體的剩余壽命.另外，由該方法得到的剩余壽命90%置信下限能夠覆蓋大部分時(shí)間截尾數(shù)據(jù)，因此該方法可以在航空公司實(shí)際運(yùn)行過(guò)程中使用.

表2 性能比較表

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

針對(duì)航空發(fā)動(dòng)機(jī)在實(shí)際運(yùn)行過(guò)程中無(wú)完全失效時(shí)間數(shù)據(jù)，建立了能夠融合性能退化數(shù)據(jù)、時(shí)間截尾數(shù)據(jù)的航空發(fā)動(dòng)機(jī)剩余壽命預(yù)測(cè)模型.利用MCMC方法對(duì)模型參數(shù)后驗(yàn)分布進(jìn)行估計(jì)，并比較了不同數(shù)據(jù)源對(duì)剩余壽命預(yù)測(cè)結(jié)果的影響.最終實(shí)例計(jì)算表明，所提出的綜合考慮性能退化數(shù)據(jù)與截尾時(shí)間數(shù)據(jù)的剩余壽命預(yù)測(cè)方法精度較高，能夠?yàn)楹娇展绢A(yù)測(cè)實(shí)際下發(fā)時(shí)間和制定維修計(jì)劃提供依據(jù).