王璽，胡昌華,*，任子強(qiáng)，熊薇

1. 火箭軍工程大學(xué) 導(dǎo)彈工程學(xué)院，西安 710025 2. 國(guó)防大學(xué) 聯(lián)合勤務(wù)學(xué)院，北京 100089

航空發(fā)動(dòng)機(jī)作為飛機(jī)的“心臟”，不僅為飛機(jī)提供動(dòng)力，同時(shí)也為航空事業(yè)的發(fā)展注入強(qiáng)大動(dòng)力。然而，航空發(fā)動(dòng)機(jī)由于長(zhǎng)期運(yùn)行在復(fù)雜多變的惡劣環(huán)境下，其性能會(huì)隨之產(chǎn)生衰減，并在累積到某種量級(jí)時(shí)，導(dǎo)致發(fā)動(dòng)機(jī)失效[1]。目前，航空發(fā)動(dòng)機(jī)在航空領(lǐng)域發(fā)揮著重要的作用，一旦發(fā)生由失效引起的事故，所造成的后果將難以估計(jì)。因此，通過(guò)監(jiān)測(cè)的航空發(fā)動(dòng)機(jī)性能衰減數(shù)據(jù)，建立相應(yīng)的性能衰減模型，進(jìn)而實(shí)現(xiàn)航空發(fā)動(dòng)機(jī)的剩余壽命(Remanining Useful Life,RUL)預(yù)測(cè)是近年來(lái)熱門的一項(xiàng)關(guān)鍵技術(shù)[2]。然而，航空發(fā)動(dòng)機(jī)的壽命有多種定義，如首翻壽命、大修壽命、總壽命等，同時(shí)又受多種因素約束，如性能衰減、關(guān)鍵件壽命等。鑒于此，本文主要研究的是在首翻壽命的定義下，基于性能衰減建模的航空發(fā)動(dòng)機(jī)剩余壽命預(yù)測(cè)。

目前，基于性能衰減建模的剩余壽命預(yù)測(cè)方法已經(jīng)廣泛地應(yīng)用于航空發(fā)動(dòng)機(jī)[3-5]、導(dǎo)彈[6-8]、電機(jī)[9]等。例如，黃亮等[3]基于多階段Wiener過(guò)程建立了航空發(fā)動(dòng)機(jī)的性能衰減模型，通過(guò)融合發(fā)動(dòng)機(jī)的歷史性能衰減數(shù)據(jù)和實(shí)時(shí)監(jiān)測(cè)數(shù)據(jù)，實(shí)現(xiàn)對(duì)發(fā)動(dòng)機(jī)的剩余壽命預(yù)測(cè)。王浩偉等[6]針對(duì)導(dǎo)彈多個(gè)部件間的競(jìng)爭(zhēng)退化，綜合考慮了退化失效和突發(fā)失效2種競(jìng)爭(zhēng)模式，提高了對(duì)導(dǎo)彈剩余壽命預(yù)測(cè)的精度。袁慶洋等[9]針對(duì)電機(jī)的退化特性，建立了多段Wiener過(guò)程模型，實(shí)現(xiàn)了電機(jī)的剩余壽命預(yù)測(cè)。

然而，對(duì)于航空發(fā)動(dòng)機(jī)這類復(fù)雜的機(jī)械系統(tǒng)，在性能衰減過(guò)程中普遍存在著三源不確定性，即時(shí)變不確定性[10]、個(gè)體差異不確定性[11]、測(cè)量不確定性[12]。時(shí)變不確定性是指性能衰減過(guò)程隨著時(shí)間變化的固有不確定性；個(gè)體差異不確定性是指同類型設(shè)備之間的性能衰減差異，即同類型設(shè)備的不同個(gè)體在性能衰減過(guò)程中雖然具有相似的衰減路徑，但實(shí)際的衰減率是不相同的；測(cè)量不確定性是指測(cè)量的性能衰減數(shù)據(jù)與實(shí)際的性能衰減狀態(tài)之間的差異。在實(shí)際中，由于受噪聲、干擾、儀器精度等影響，想要對(duì)性能衰減狀態(tài)的精確測(cè)量幾乎是不可能的，因此測(cè)量中難免會(huì)引入測(cè)量誤差，造成測(cè)量的不確定性。

目前，有關(guān)同時(shí)考慮三源不確定性的研究較少。為此，文獻(xiàn)[13]基于線性的Wiener過(guò)程建立了同時(shí)考慮三源不確定性的性能衰減模型，并使得預(yù)測(cè)的剩余壽命分布可以隨著性能衰減數(shù)據(jù)在線更新。然而，對(duì)于航空發(fā)動(dòng)機(jī)而言，在性能衰減過(guò)程中不僅具有三源不確定性，還具有強(qiáng)非線性。鑒于此，文獻(xiàn)[14]在非線性下考慮了三源不確定性，并且證明了同時(shí)考慮三源不確定性可以有效地提高剩余壽命預(yù)測(cè)的精度。文獻(xiàn)[15]在此基礎(chǔ)上，利用粒子濾波實(shí)現(xiàn)了對(duì)鋰電池的剩余壽命預(yù)測(cè)。但是，以上有關(guān)考慮三源不確定性的工作中存在一個(gè)共性問(wèn)題，即假定有多組同類型設(shè)備的歷史性能衰減數(shù)據(jù)來(lái)離線地辨識(shí)模型初始參數(shù)，并且這些參數(shù)一旦確定，就不再隨著新的性能衰減數(shù)據(jù)進(jìn)行更新。因此，剩余壽命預(yù)測(cè)的精度受限于數(shù)據(jù)量，并且預(yù)測(cè)的結(jié)果在很大程度上依賴于選擇的模型初始參數(shù)。然而，針對(duì)新研發(fā)的航空發(fā)動(dòng)機(jī)，往往是缺乏此類相關(guān)的歷史性能衰減數(shù)據(jù)和先驗(yàn)信息，因此，以上方法具有一定的局限性，且目前相關(guān)的研究較少。

對(duì)于目前航空發(fā)動(dòng)機(jī)剩余壽命預(yù)測(cè)存在的問(wèn)題，本文建立了同時(shí)考慮三源不確定性和非線性的性能衰減模型，并在首達(dá)時(shí)間下推導(dǎo)出剩余壽命的分布。然后，針對(duì)新研發(fā)的航空發(fā)動(dòng)機(jī)，提出了一種基于Kalman濾波和條件期望最大化(Expectation Conditional Maximization,ECM)算法的參數(shù)估計(jì)方法，以克服缺乏歷史性能衰減數(shù)據(jù)和先驗(yàn)信息的問(wèn)題，并當(dāng)新的監(jiān)測(cè)數(shù)據(jù)可用時(shí)，能夠自適應(yīng)地估計(jì)模型參數(shù)，進(jìn)而更新剩余壽命的分布，實(shí)現(xiàn)了對(duì)新研發(fā)航空發(fā)動(dòng)機(jī)剩余壽命的在線預(yù)測(cè)。最后，通過(guò)實(shí)例驗(yàn)證了本文方法可以有效地提高預(yù)測(cè)的精度并降低預(yù)測(cè)的不確定性，且更加適用于新研發(fā)的航空發(fā)動(dòng)機(jī)。

1 模型描述

針對(duì)航空發(fā)動(dòng)機(jī)在性能衰減過(guò)程中存在的非線性和不確定性，本文采用一類一般性的非線性Wiener過(guò)程來(lái)建模航空發(fā)動(dòng)機(jī)的性能衰減過(guò)程。具體地，性能衰減過(guò)程{X(t),t≥0}由標(biāo)準(zhǔn)Brown運(yùn)動(dòng)B(t)驅(qū)動(dòng)，可以表示為

(1)

Y(t)=X(t)+ε

(2)

在以上框架下，本文的主要目標(biāo)是針對(duì)缺乏歷史數(shù)據(jù)和先驗(yàn)信息的新研發(fā)航空發(fā)動(dòng)機(jī)，利用實(shí)時(shí)監(jiān)測(cè)的性能衰減數(shù)據(jù)來(lái)自適應(yīng)地預(yù)測(cè)航空發(fā)動(dòng)機(jī)的剩余壽命，并且在獲得新的性能衰減數(shù)據(jù)后實(shí)現(xiàn)對(duì)剩余壽命的在線更新。具體地，首先假定獲得的性能衰減數(shù)據(jù)是在離散時(shí)間點(diǎn)0=t0

(3)

2 剩余壽命自適應(yīng)預(yù)測(cè)

考慮到表征設(shè)備健康水平的潛在性能衰減狀態(tài){X(t),t≥0}一旦首次達(dá)到設(shè)定的失效閾值ω(一般由工業(yè)標(biāo)準(zhǔn)設(shè)定)時(shí)，則認(rèn)為設(shè)備失效。因此，基于首達(dá)時(shí)間的概念[1-5]，設(shè)備在tn時(shí)刻的剩余壽命Ln定義為

Ln=inf{ln>0:X(ln+tn)≥ω}

(4)

相應(yīng)地，剩余壽命Ln的概率密度函數(shù)(PDF)為fLn|Y0:n(ln|Y0:n)。因此，預(yù)測(cè)航空發(fā)動(dòng)機(jī)剩余壽命的關(guān)鍵是基于監(jiān)測(cè)的性能衰減數(shù)據(jù)Y0:n={y0,y1,…,yn}，求解fLn|Y0:n(ln|Y0:n)，并在獲得新的性能衰減數(shù)據(jù)后，實(shí)現(xiàn)剩余壽命分布的在線更新。

為了實(shí)現(xiàn)剩余壽命的自適應(yīng)預(yù)測(cè)，首先將式(3)進(jìn)一步改寫為

(5)

(6)

式中：δn～N(0,φn)，且有

定義基于監(jiān)測(cè)數(shù)據(jù)Y0:n估計(jì)的zn的期望和協(xié)方差分別為

(7)

同樣地，定義一步前向的期望和協(xié)方差為

(8)

狀態(tài)估計(jì)

(9)

協(xié)方差更新

Pn|n=Pn|n-1-KnCPn|n-1

(10)

濾波初值

(11)

根據(jù)建立的狀態(tài)空間模型式(6)，基于監(jiān)測(cè)的性能衰減數(shù)據(jù)Y0:n，計(jì)算剩余壽命的PDF為

fLn|Y0:n(ln|Y0:n)?

(12)

3 基于ECM算法的模型參數(shù)估計(jì)

(13)

根據(jù)極大似然估計(jì)(MLE)理論，性能衰減數(shù)據(jù)Y0:n關(guān)于模型參數(shù)向量Θ的對(duì)數(shù)似然函數(shù)為

(14)

(15)

然而，由于將潛在性能衰減狀態(tài)xn和隨機(jī)參數(shù)λn視作“隱含”狀態(tài)，因此無(wú)法直接最大化式(15)。為解決這一問(wèn)題，自然會(huì)采用期望最大化(EM)算法進(jìn)行參數(shù)估計(jì)，通過(guò)最大化聯(lián)合似然函數(shù)p(zn,Y0:n|Θ)去逼近模型參數(shù)的MLE。EM算法的估計(jì)過(guò)程可以分為以下2步。

E步：

(16)

M步：

(17)

首先，為了計(jì)算對(duì)數(shù)聯(lián)合似然函數(shù)，可以通過(guò)條件概率公式得到

2lnp(zn,Y0:n|Θ)=

2[lnp(Y0:n|zn,Θ)+lnp(zn|Θ)]=

(18)

然后，通過(guò)E步對(duì)式(18)求條件期望，有

(19)

式中：

(20)

顯然，為了推導(dǎo)出E步結(jié)果，需要先計(jì)算式(20)。這里，本文采用Rauch-Tung-Striebel(RTS)平滑濾波算法[19]計(jì)算Ω、Ξ、Φ：

平滑濾波

(21)

協(xié)方差初值

Mn|n=(I2×2-KnC)AnPn-1|n-1

(22)

后向迭代

(23)

式中：Mi|n=cov(zi,zi-1|Y0:n)?；赗TS算法，條件期望Ω、Ξ、Φ可以計(jì)算為

(24)

(25)

(26)

(27)

(28)

CM2：為了基于性能衰減數(shù)據(jù)Y0:n估計(jì)參數(shù)θ，首先將式(25)和式(28)一起代入到式(19)中，便可以通過(guò)最大化來(lái)估計(jì)參數(shù)即

(29)

由此可見，與傳統(tǒng)的EM算法相比，本文的ECM算法通過(guò)CM1和CM2兩步，可以將原本需要優(yōu)化的3個(gè)參數(shù)簡(jiǎn)化為一個(gè)參數(shù)θ，降低了計(jì)算復(fù)雜度，提高了在線性能。此外，本文提出的參數(shù)估計(jì)方法可以在監(jiān)測(cè)到一個(gè)新的性能衰減數(shù)據(jù)時(shí)，實(shí)現(xiàn)對(duì)模型參數(shù)的自適應(yīng)估計(jì)和在線更新，克服了新研發(fā)航空發(fā)動(dòng)機(jī)歷史數(shù)據(jù)和先驗(yàn)信息不足的限制，這對(duì)于工程實(shí)際具有重要意義。

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

根據(jù)相關(guān)文獻(xiàn)可知，航空發(fā)動(dòng)機(jī)的性能衰減和故障通常表現(xiàn)在排氣溫度裕度(EGTM)的衰減[22-23]。因此，EGTM是表征航空發(fā)動(dòng)機(jī)健康狀態(tài)的一個(gè)重要指標(biāo)。EGTM的定義為航空發(fā)動(dòng)機(jī)出廠前測(cè)定的排氣溫度(EGT)限制值與全推力(或全功率)起飛時(shí)實(shí)際的EGT之間的差值。隨著航空發(fā)動(dòng)機(jī)工作時(shí)間的增加，EGT會(huì)逐漸升高，反映到EGTM就會(huì)逐漸衰退，因此準(zhǔn)確預(yù)測(cè)EGTM對(duì)提高航空發(fā)動(dòng)機(jī)的可靠性至關(guān)重要。

實(shí)際中，由于外部環(huán)境和內(nèi)部因素的隨機(jī)性，航空發(fā)動(dòng)機(jī)在某一時(shí)刻的性能衰減量也是隨機(jī)的，而且同類型航空發(fā)動(dòng)機(jī)EGTM的衰減軌跡雖然具有相似性，但實(shí)際上不同個(gè)體之間仍然存在著差異性。因此，監(jiān)測(cè)的EGTM數(shù)據(jù)不僅會(huì)受到時(shí)變不確定性的影響，還會(huì)受到個(gè)體差異不確定性的影響。此外，EGT數(shù)據(jù)來(lái)源于航空發(fā)動(dòng)機(jī)上安裝的傳感器，經(jīng)過(guò)空地?cái)?shù)據(jù)鏈、報(bào)文解析程序等步驟進(jìn)入到性能監(jiān)測(cè)系統(tǒng)，在這整個(gè)過(guò)程中不可避免地會(huì)引入測(cè)量誤差，造成測(cè)量的不確定性。根據(jù)EGTM的定義可知，EGTM是基于EGT計(jì)算的，因此EGTM數(shù)據(jù)中必然存在測(cè)量誤差。通過(guò)上述分析發(fā)現(xiàn)，監(jiān)測(cè)的EGTM數(shù)據(jù)中存在著三源不確定性共存的現(xiàn)象，因此同時(shí)考慮三源不確定性是必要的，也是更加符合實(shí)際需求的。下面以某型號(hào)航空發(fā)動(dòng)機(jī)在6年內(nèi)監(jiān)測(cè)的EGTM數(shù)據(jù)來(lái)驗(yàn)證本文方法的有效性，收集的性能衰減數(shù)據(jù)如圖1所示。根據(jù)工程經(jīng)驗(yàn)，設(shè)定EGTM的失效閾值ω=0 ℃，并且考慮以下冪函數(shù)形式的性能衰減模型：

圖1 航空發(fā)動(dòng)機(jī)的性能衰減路徑

X(t)=x0+λ·tθ+σBB(t)

(30)

這里，根據(jù)EGTM數(shù)據(jù)的特點(diǎn)，選擇式(1)中非線性函數(shù)的形式為ρ(t;θ)=θ·tθ-1。到目前為止，此類模型已經(jīng)在航空發(fā)動(dòng)機(jī)、軸承、電池等復(fù)雜系統(tǒng)的性能衰減建模中得到了廣泛的應(yīng)用[3,12,17]。

如圖1所示，監(jiān)測(cè)間隔為100次循環(huán)，在監(jiān)測(cè)到第4 000次循環(huán)時(shí)，EGTM首次達(dá)到失效閾值ω，期間共監(jiān)測(cè)了41次。顯然，隨著循環(huán)次數(shù)的增加，航空發(fā)動(dòng)機(jī)的性能衰減路徑整體呈下降趨勢(shì)。然而，本文所提出的剩余壽命預(yù)測(cè)方法是基于性能衰減路徑呈上升趨勢(shì)的情況下推導(dǎo)出來(lái)的，因此需要對(duì)原始性能衰減數(shù)據(jù)進(jìn)行轉(zhuǎn)換。這里，本文利用EGTM的初始值減去其他所有的性能衰減數(shù)據(jù)，相應(yīng)地，失效閾值轉(zhuǎn)換為初始值減去ω。下面，基于轉(zhuǎn)換后的性能衰減數(shù)據(jù)來(lái)驗(yàn)證本文方法的有效性和優(yōu)越性。隨著性能衰減數(shù)據(jù)的積累，采用所提出的參數(shù)估計(jì)方法可以對(duì)模型參數(shù)Θ在每一個(gè)監(jiān)測(cè)點(diǎn)進(jìn)行自適應(yīng)地估計(jì)和更新，結(jié)果如圖2所示。

圖2 模型參數(shù)的自適應(yīng)估計(jì)過(guò)程

圖2的結(jié)果表明，模型參數(shù)可以隨著性能衰減數(shù)據(jù)的積累很快地收斂，并且在獲得一個(gè)新的性能衰減數(shù)據(jù)時(shí)，模型參數(shù)和剩余壽命的PDF都可以進(jìn)行自適應(yīng)地估計(jì)和更新。因此，克服了歷史數(shù)據(jù)和先驗(yàn)信息不足的限制，能夠適用于新研發(fā)航空發(fā)動(dòng)機(jī)的剩余壽命預(yù)測(cè)。

為了進(jìn)一步證明本文方法的優(yōu)越性，選擇以下2種方法進(jìn)行對(duì)比研究：

圖3給出了3種方法在最后6個(gè)監(jiān)測(cè)點(diǎn)預(yù)測(cè)的剩余壽命的PDF曲線，從圖中可以看出，3種方法得到的PDF曲線均包含了實(shí)際的剩余壽命，且隨著循環(huán)次數(shù)的增加逐漸變窄，這意味著預(yù)測(cè)的剩余壽命的不確定性越來(lái)越低。此外，由于本文方法在性能衰減建模中同時(shí)考慮了三源不確定性，并且在每一個(gè)監(jiān)測(cè)點(diǎn)可以對(duì)模型參數(shù)進(jìn)行自適應(yīng)地估計(jì)和在線更新，因此本文方法得到的PDF曲線更加窄而尖銳，說(shuō)明本文方法在剩余壽命預(yù)測(cè)方面具有更高的精度和更好的性能。

圖3 三種方法的比較結(jié)果

為了進(jìn)一步量化比較結(jié)果，本文采用可靠性領(lǐng)域中常用的均方誤差(MSE)指標(biāo)，該指標(biāo)同時(shí)考慮了剩余壽命預(yù)測(cè)的精度和不確定性，在tn時(shí)刻的定義式為

(31)

圖4 所有監(jiān)測(cè)點(diǎn)的MSE

圖5 4個(gè)不同監(jiān)測(cè)點(diǎn)的箱形圖(點(diǎn)劃線表示實(shí)際剩余壽命)

圖4給出了3種方法在所有監(jiān)測(cè)點(diǎn)的MSE比較結(jié)果。由于本文方法和M2方法均不依賴于歷史數(shù)據(jù)和先驗(yàn)信息，且模型的初始參數(shù)是隨機(jī)給定的，因此在監(jiān)測(cè)初期性能衰減數(shù)據(jù)較少時(shí)，剩余壽命預(yù)測(cè)的MSE高于M1，但隨著性能衰減數(shù)據(jù)的積累，本文方法的MSE明顯低于M1和M2。為進(jìn)一步驗(yàn)證，圖5給出了4個(gè)不同監(jiān)測(cè)點(diǎn)的箱形圖。其中，箱子的大小表示預(yù)測(cè)的剩余壽命的不確定性，若箱子越大，則預(yù)測(cè)的不確定性越大，精度越低；點(diǎn)劃線表示當(dāng)前監(jiān)測(cè)點(diǎn)的實(shí)際剩余壽命；箱子中的實(shí)線表示預(yù)測(cè)的剩余壽命的期望，實(shí)線越接近點(diǎn)劃線，說(shuō)明預(yù)測(cè)的精度越高，結(jié)果越準(zhǔn)確。從圖5(a)中可以看出，此時(shí)由于性能衰減數(shù)據(jù)較少，M1優(yōu)于本文方法和M2，但是隨著性能衰減數(shù)據(jù)的積累，本文方法的箱子更小且更加接近實(shí)際的剩余壽命，如圖5(b)～圖5(d)所示。此外，引入相對(duì)誤差(RE)指標(biāo)來(lái)進(jìn)一步量化剩余壽命預(yù)測(cè)的精度，在壽命的20%、45%、70%和95%分位點(diǎn)上分別給出了3種方法RE的比較結(jié)果，如表1所示。

通過(guò)表1的結(jié)果可以看出，本文方法可以有效減少剩余壽命預(yù)測(cè)的RE，進(jìn)而提高了剩余壽命預(yù)測(cè)的精度，尤其在壽命的95%分位點(diǎn)處，剩余壽命預(yù)測(cè)的RE僅為1.68%，明顯小于其他2種方法。以上的實(shí)驗(yàn)結(jié)果表明，本文方法具有更好的預(yù)測(cè)能力，剩余壽命預(yù)測(cè)的精度更高，且不依賴于歷史數(shù)據(jù)和先驗(yàn)信息，更加適用于新研發(fā)的航空發(fā)動(dòng)機(jī)。

表1 不同壽命分位點(diǎn)相對(duì)誤差的比較結(jié)果

5 結(jié) 論

1) 本文在三源不確定性下建立了非線性的性能衰減模型，不僅克服了現(xiàn)有方法的潛在假設(shè)，而且可以較好地建模航空發(fā)動(dòng)機(jī)的性能衰減軌跡，并能夠?qū)⑷床淮_定性的影響納入到預(yù)測(cè)的剩余壽命分布中。

2) 利用Kalman濾波技術(shù)可以對(duì)潛在的性能衰減狀態(tài)進(jìn)行實(shí)時(shí)估計(jì)，并基于首達(dá)時(shí)間的概念，得到了同時(shí)考慮三源不確定性的剩余壽命分布。

3) 在獲得一個(gè)新的性能衰減數(shù)據(jù)時(shí)，基于RTS平滑和ECM算法對(duì)模型參數(shù)進(jìn)行自適應(yīng)地估計(jì)和在線更新，克服了新研發(fā)航空發(fā)動(dòng)機(jī)缺乏歷史數(shù)據(jù)和先驗(yàn)信息的問(wèn)題，降低了參數(shù)估計(jì)的不確定性，顯著提高了剩余壽命預(yù)測(cè)的精度，尤其在壽命的95%分位點(diǎn)處，剩余壽命預(yù)測(cè)的RE僅為1.68%。