本發(fā)明涉及航空渦扇發(fā)動機(jī)剩余壽命預(yù)測領(lǐng)域，具體是一種多源數(shù)據(jù)融合的航空渦扇發(fā)動機(jī)剩余壽命預(yù)測方法。

背景技術(shù)：

航空發(fā)動機(jī)預(yù)測與健康管理技術(shù)(prognosticsandhealthmanagement，phm)為發(fā)動機(jī)視情維修(condition‐basedmaintenance，cm)提供技術(shù)指導(dǎo)和先驗(yàn)知識，大力推動了航空發(fā)動機(jī)維修技術(shù)逐漸由定時(shí)維修向視情維修技術(shù)的轉(zhuǎn)變，確保飛行安全、降低維修費(fèi)用，因此引起了廣泛的關(guān)注和研究。隨著傳感器技術(shù)和電子設(shè)備的進(jìn)步，發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)能夠被準(zhǔn)確的監(jiān)測和記錄，為phm技術(shù)提供了足夠的數(shù)據(jù)支持。如何根據(jù)已有的發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)，對發(fā)動機(jī)剩余壽命進(jìn)行準(zhǔn)確的預(yù)測，是航空發(fā)動機(jī)phm技術(shù)的研究核心。

基于統(tǒng)計(jì)數(shù)據(jù)驅(qū)動的方法通過挖掘監(jiān)測數(shù)據(jù)包含的潛在信息，以顯示系統(tǒng)內(nèi)部變化以及外部環(huán)境影響所導(dǎo)致系統(tǒng)的性能變化，避免了對復(fù)雜精密度高的航空發(fā)動機(jī)系統(tǒng)機(jī)理進(jìn)行建模，同時(shí)兼顧預(yù)測準(zhǔn)確性，直接實(shí)現(xiàn)渦扇發(fā)動機(jī)的壽命預(yù)測，使得預(yù)測結(jié)果具有模糊性，為今后的研究中，融入專家知識進(jìn)行更精確的預(yù)測提供可行性。統(tǒng)計(jì)數(shù)據(jù)驅(qū)動方法需要對系統(tǒng)退化過程建模，維納過程符合系統(tǒng)非線性非平穩(wěn)退化過程，同時(shí)具有良好的數(shù)學(xué)特性，基于維納過程的退化模型已經(jīng)在航空發(fā)動機(jī)系統(tǒng)退化建模分析和剩余壽命預(yù)測等方面取得了不錯的成效，充分證明了該模型的適用性。

在實(shí)際應(yīng)用中，航空發(fā)動機(jī)的排氣溫度裕度作為一項(xiàng)表征發(fā)動機(jī)換發(fā)依據(jù)和重要的性能指標(biāo)，常被用于建立發(fā)動機(jī)退化模型?；诰€性漂移的維納過程概率分布函數(shù)明確，模型簡單。然而，在現(xiàn)有的航空發(fā)動機(jī)剩余壽命研究中，卻忽略了其他監(jiān)測數(shù)據(jù)中攜帶的部件狀態(tài)信息、發(fā)動機(jī)運(yùn)行環(huán)境和個體之間的差異性以及發(fā)動機(jī)系統(tǒng)運(yùn)行后期退化加速的特點(diǎn)，沒有同時(shí)考慮多源的監(jiān)測數(shù)據(jù)和帶非線性漂移退化過程的結(jié)合。

技術(shù)實(shí)現(xiàn)要素：

針對傳統(tǒng)的航空發(fā)動機(jī)剩余壽命預(yù)測方法的不足，本發(fā)明通過提供一種多源統(tǒng)計(jì)數(shù)據(jù)驅(qū)動的航空渦扇發(fā)動機(jī)剩余壽命預(yù)測方法，以解決傳統(tǒng)統(tǒng)計(jì)數(shù)據(jù)驅(qū)動預(yù)測方法中忽略發(fā)動機(jī)多維監(jiān)測數(shù)據(jù)和設(shè)備運(yùn)行環(huán)境及個體差異的問題。

為實(shí)現(xiàn)上述目的，本發(fā)明的構(gòu)思與技術(shù)解決方案實(shí)現(xiàn)如下：

本發(fā)明的基本構(gòu)思是，充分利用發(fā)動機(jī)傳感器采集的有效監(jiān)測數(shù)據(jù)，在不利于壽命預(yù)測冗余特征的基礎(chǔ)上，進(jìn)行信息融合以提取表征發(fā)動機(jī)運(yùn)行狀態(tài)的健康指標(biāo)和失效閾值，解決傳統(tǒng)預(yù)測方法中監(jiān)測數(shù)據(jù)信息利用不充分的問題。建立基于帶隨機(jī)參數(shù)的非線性維納過程航空發(fā)動機(jī)失效模型，用以表征運(yùn)行環(huán)境的不確定性以及設(shè)備個體間的差異性，模擬發(fā)動機(jī)的非線性退化過程。在此基礎(chǔ)上，對航空發(fā)動機(jī)進(jìn)行實(shí)時(shí)在線壽命預(yù)測。

基于上述基本構(gòu)思，本發(fā)明提供的技術(shù)解決方案為一種多源數(shù)據(jù)融合的航空渦扇發(fā)動機(jī)剩余壽命預(yù)測方法，包括以下步驟：

步驟1：基于多源信息融合的渦扇發(fā)動機(jī)性能退化評估：

步驟1.1：多源監(jiān)測數(shù)據(jù)降維，基于共同主成分分析法提取系統(tǒng)監(jiān)測數(shù)據(jù)的主成分；監(jiān)測數(shù)據(jù)為多維時(shí)間序列，保持時(shí)間維不變，對變量維進(jìn)行共同主成分提?。?/p>

其中u＝[u1,u2,...,ur]^t為變換矩陣，r表示所選取主分量的數(shù)目，ur為第r個主分量對應(yīng)的單位特征向量，xn為第n臺航空發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)，為基于共同主成分分析降維后的第n臺航空發(fā)動機(jī)監(jiān)測數(shù)據(jù)，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)，為第n臺航空發(fā)動機(jī)降維后第mn次采樣的監(jiān)測數(shù)據(jù)，(·)^t表示轉(zhuǎn)置。

步驟1.2：在步驟1.1的基礎(chǔ)上，進(jìn)行發(fā)動機(jī)健康參數(shù)的統(tǒng)計(jì)，研究表明發(fā)動機(jī)的前5％壽命周期被視作沒有發(fā)生退化，故提取訓(xùn)練集中的未退化監(jiān)測數(shù)據(jù)視作發(fā)動機(jī)的健康參數(shù)；

發(fā)動機(jī)總體健康參數(shù)為h＝(h1,h2,...hn,...,hn)^t，hn為第n臺航空發(fā)動機(jī)健康參數(shù)，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)；則航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值m表示如下：

其中m為航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值，sn表示訓(xùn)練集中第n臺航空發(fā)動機(jī)健康參數(shù)的最大采樣次數(shù)，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，∑表示求和；

步驟1.3：在步驟1.2的基礎(chǔ)上，計(jì)算發(fā)動機(jī)運(yùn)行健康指標(biāo)，以評估發(fā)動機(jī)性能退化程度；

每臺發(fā)動機(jī)的健康指標(biāo)計(jì)算公式如下：

其中tj為第n臺發(fā)動機(jī)第j次采樣的時(shí)間，yn(tj)為第j次采樣時(shí)第n臺發(fā)動機(jī)健康指標(biāo)，m為航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)；

步驟1.4：基于發(fā)動機(jī)健康指標(biāo)，確定發(fā)動機(jī)失效閾值；研究表明，航空發(fā)動機(jī)失效過程劃分為四個階段，最后一階段視作發(fā)動機(jī)失效階段；計(jì)算最后一類的聚類中心作為發(fā)動機(jī)失效閾值ω；

步驟2：航空發(fā)動機(jī)退化建模及參數(shù)估計(jì)；

步驟2.1：選取帶非線性漂移的維納過程為表征發(fā)動機(jī)退化過程模型，令初始退化狀態(tài)為零；

y(t)＝y(tǒng)(0)+at^b+σb(t)(7)

其中，y(0)＝0，t為采樣時(shí)間，非線性漂移部分為at^b，擴(kuò)散系數(shù)為σ，b(t)為布朗運(yùn)動；為表征制造偏差和運(yùn)行環(huán)境導(dǎo)致發(fā)動機(jī)個體之間的差異性，a為非線性漂移部分的隨機(jī)變量，滿足其中n(·)表示正態(tài)分布，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差；同時(shí)將b定義為非線性漂移部分的常量表征發(fā)動機(jī)的相同特性；

步驟2.2：基于最大似然法的航空發(fā)動機(jī)退化模型參數(shù)估計(jì)：在建立渦扇發(fā)動機(jī)的退化模型后，需要基于發(fā)動機(jī)健康指標(biāo)對模型進(jìn)行參數(shù)估計(jì)；每臺發(fā)動機(jī)在多個采樣點(diǎn)的健康指標(biāo)服從多維正態(tài)分布，故yn～n(μn,σn)，其中yn為第n臺發(fā)動機(jī)的健康指標(biāo)，μn為第n臺發(fā)動機(jī)的健康指標(biāo)均值向量；σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣；

對于不同發(fā)動機(jī)樣本而言，退化過程獨(dú)立同分布，因此關(guān)于總體發(fā)動機(jī)健康指標(biāo)y的似然函數(shù)的自然對數(shù)形式如公式(12)所示：

其中a為非線性漂移部分的隨機(jī)變量，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，σ為擴(kuò)散系數(shù)，b為非線性漂移部分的常量，y為總體發(fā)動機(jī)健康指標(biāo)，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，|·|表示取模操作，為σn的逆矩陣，yn為第n臺發(fā)動機(jī)的健康指標(biāo)，μn為第n臺發(fā)動機(jī)的健康指標(biāo)均值向量；

通過最大化似然函數(shù)進(jìn)行參數(shù)估計(jì)；

其中a為非線性漂移部分的隨機(jī)變量，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，σ為擴(kuò)散系數(shù)，b為非線性漂移部分的常量，為μa的估計(jì)值，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，tj為第n臺發(fā)動機(jī)第j次采樣的時(shí)間，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，|·|表示取模操作，為σn的逆矩陣，yn為第n臺發(fā)動機(jī)的健康指標(biāo)，y為總體發(fā)動機(jī)健康指標(biāo)；

步驟3：建立航空發(fā)動機(jī)剩余壽命預(yù)測模型；

步驟3.1：基于維納過程的航空發(fā)動機(jī)剩余壽命描述

基于發(fā)動機(jī)性能退化模型具有時(shí)間平移不變性，則退化模型表示為：

其中tc為航空發(fā)動機(jī)當(dāng)前運(yùn)行時(shí)間，li為航空發(fā)動機(jī)剩余壽命，b(·)表示布朗運(yùn)動，y(tc)為航空發(fā)動機(jī)當(dāng)前健康指標(biāo)，y(li+tc)為航空發(fā)動機(jī)失效時(shí)的健康指標(biāo)，a為非線性漂移部分的隨機(jī)變量，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)，y(li)表示發(fā)動機(jī)隨航空發(fā)動機(jī)剩余壽命li退化的維納過程，∫表示積分；可知發(fā)動機(jī)剩余壽命描述仍可表征為帶非線性漂移的維納過程；

步驟3.2：航空發(fā)動機(jī)剩余壽命描述轉(zhuǎn)換為標(biāo)準(zhǔn)維納過程，對于非線性擴(kuò)散過程，若滿足：

其中μ(y,li)表示非線性漂移部分，σ′(li,u)為擴(kuò)散部分，y為健康指標(biāo)，li為航空發(fā)動機(jī)剩余壽命，c1(li)和c2(li)為任意關(guān)于li的函數(shù)，z為維納過程初始值，∫表示積分，表示偏導(dǎo)數(shù)；

步驟4：基于性能衰退的航空發(fā)動機(jī)剩余壽命預(yù)測；

建立航空發(fā)動機(jī)剩余壽命概率密度函數(shù)模型，在步驟1、2、3的基礎(chǔ)上，進(jìn)行實(shí)時(shí)在線壽命預(yù)測；

布朗運(yùn)動首次穿過失效閾值的概率密度函數(shù)表示如下：

其中pb(t)(ω,t)表示布朗運(yùn)動首次穿過失效閾值的概率密度函數(shù)，ω為失效閾值，t為采樣時(shí)間；

故航空發(fā)動機(jī)的剩余壽命概率密度函數(shù)為：

其中表示退化過程首次穿過閾值的概率密度函數(shù)，ω為失效閾值，li為航空發(fā)動機(jī)剩余壽命，tc為航空發(fā)動機(jī)當(dāng)前運(yùn)行時(shí)間，y(tc)為航空發(fā)動機(jī)當(dāng)前健康指標(biāo)，a為非線性漂移部分的隨機(jī)變量，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，f(a)為隨機(jī)變量a的概率密度函數(shù)，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)，exp(·)表示以e為底的指數(shù)函數(shù)；以上結(jié)果滿足航空發(fā)動機(jī)壽命預(yù)測的實(shí)際使用需要。

步驟1.4所述計(jì)算最后一類的聚類中心作為發(fā)動機(jī)失效閾值ω，具體步驟如下：

第一步，對于第n臺發(fā)動機(jī)健康指標(biāo)每次采樣的健康指標(biāo)視作一個樣本；

第二步，隨機(jī)選擇4個初始聚類中心z1(k),z2(k)…z4(k)，k表示迭代次數(shù)，初始值為1，每迭代一次加1；

第三步，逐個將yn(tp),p＝1,2,...,mn按最小距離準(zhǔn)則分配到4個中心為zi(k)的類中；

第三步，分計(jì)算各類的新中心：nj即第j類中的樣本個數(shù)，oj表示屬于第j類的樣本集合；

第四步，如果zj(k+1)≠zj(k)則轉(zhuǎn)到第2步，將模式重新分類，重復(fù)迭代計(jì)算，直到聚類中心收斂，記錄4個聚類中心的最大值ωmax(n)；

第五步，對于每一臺發(fā)動機(jī)，重復(fù)第1、2、3、4步，選取ωmax(n)中的最小值，視作發(fā)動機(jī)失效閾值，即ω＝min(ωmax(n)),n＝1,2,...,100,其中min(·)表示取最小值。

本發(fā)明的有益效果在：

本發(fā)明多源數(shù)據(jù)融合的航空渦扇發(fā)動機(jī)剩余壽命預(yù)測方法，基于多源監(jiān)測數(shù)據(jù)，采用共同主成分分析提取監(jiān)測數(shù)據(jù)的主成分，降低冗余數(shù)據(jù)和噪聲數(shù)據(jù)的影響。多參數(shù)的融合能夠最大限度利用傳感器監(jiān)測的有效信息，從而準(zhǔn)確表征航空發(fā)動機(jī)的退化過程。發(fā)動機(jī)退化后期退化過程加速，為此建立帶非線性漂移的維納過程。建立航空發(fā)動機(jī)剩余壽命預(yù)測模型，實(shí)現(xiàn)對發(fā)動機(jī)的實(shí)時(shí)在線預(yù)測，為視情維修提供先驗(yàn)知識，具備良好的工程應(yīng)用價(jià)值。

附圖說明

圖1為本發(fā)明方法流程圖。

圖2為訓(xùn)練集#1～#10臺航空發(fā)動機(jī)高壓壓氣機(jī)入口總溫度測量數(shù)據(jù)。

圖3為經(jīng)共同主成分分析降維后前六主分量取值。

圖4a為#1號訓(xùn)練發(fā)動機(jī)基于ed的發(fā)動機(jī)健康指標(biāo)曲線。

圖4b為#3號訓(xùn)練發(fā)動機(jī)基于ed的發(fā)動機(jī)健康指標(biāo)曲線。

圖5為航空發(fā)動機(jī)退化過程4個階段。

圖6為#1，#5，#10，#26，#60五臺測試發(fā)動機(jī)剩余壽命概率密度函數(shù)和預(yù)測剩余壽命，以及真實(shí)剩余壽命的對比。

圖7為剩余壽命誤差區(qū)間示意。

圖8為測試機(jī)中100臺發(fā)動機(jī)剩余壽命頻數(shù)分布

圖9為測試機(jī)中100臺發(fā)動機(jī)剩余壽命預(yù)測懲罰函數(shù)得分

具體實(shí)施方式

下面通過一組由nasa于2008年發(fā)布c‐mapss數(shù)據(jù)集來驗(yàn)證本發(fā)明提出的方法。該數(shù)據(jù)集利用渦扇發(fā)動機(jī)仿真模型，在不同操作條件下，模擬產(chǎn)生不同程度的特定失效狀況時(shí)渦扇發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)。數(shù)據(jù)集主要分訓(xùn)練數(shù)據(jù)集、測試數(shù)據(jù)集、剩余壽命數(shù)據(jù)集三部分。其中，訓(xùn)練數(shù)據(jù)集中采集了100臺渦扇發(fā)動機(jī)從健康狀態(tài)運(yùn)行至失效狀態(tài)的全部監(jiān)測數(shù)據(jù)。測試數(shù)據(jù)集包含有100臺發(fā)動機(jī)的監(jiān)測數(shù)據(jù)，均從健康狀態(tài)開始運(yùn)行，在系統(tǒng)失效前的某時(shí)刻停止采集。剩余壽命數(shù)據(jù)集記錄了測試發(fā)動機(jī)的真實(shí)剩余壽命。以#1～#10臺航空發(fā)動機(jī)的高壓壓氣機(jī)入口總溫度測量數(shù)據(jù)為例，訓(xùn)練集監(jiān)測數(shù)據(jù)如圖2所示。

實(shí)施例：

如圖1所示，本發(fā)明基于多源統(tǒng)計(jì)數(shù)據(jù)驅(qū)動的航空渦扇發(fā)動機(jī)剩余壽命預(yù)測方法，包括以下步驟：

1、基于多源信息融合的渦扇發(fā)動機(jī)性能退化評估，采用共同主成分分析和基于歐氏距離(euclideandistance，ed)融合方法，在去除噪聲干擾及不相關(guān)監(jiān)測數(shù)據(jù)的條件下，融合發(fā)動機(jī)有效特征信息，獲取總體發(fā)動機(jī)健康指標(biāo)y。

1.1、多源監(jiān)測數(shù)據(jù)降維，基于共同主成分分析法提取系統(tǒng)監(jiān)測數(shù)據(jù)的主成分。監(jiān)測數(shù)據(jù)為多維時(shí)間序列，保持時(shí)間維不變，對變量維進(jìn)行共同主成分提取。

第n臺航空發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)n＝(1,2,3,..,n)，xni為第n臺航空發(fā)動機(jī)第i次采樣的原始數(shù)據(jù)，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，本實(shí)施例中n＝100，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)。xni＝[xn1,xn2,...,xnj,...xnk]，k表示原始監(jiān)測數(shù)據(jù)的維數(shù)，本實(shí)施例所使用的的數(shù)據(jù)集中原始監(jiān)測數(shù)據(jù)的維數(shù)為21，故k＝21，xnj為第n臺發(fā)動機(jī)第i次采樣數(shù)據(jù)中第j個傳感器測量數(shù)據(jù)。

對于訓(xùn)練數(shù)據(jù)集，基于共同主成分分析法獲取每臺發(fā)動機(jī)的協(xié)方差矩陣，公式如下：

cn＝e[(xn-e(xn))^t(xn-e(xn))](1)

其中cn為第n航空臺發(fā)動機(jī)的協(xié)方差矩陣，xn為第n臺航空發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)，e(·)表示期望，(·)^t表示轉(zhuǎn)置；

對于n臺發(fā)動機(jī)，計(jì)算平均協(xié)方差矩陣如下：

其中為航空發(fā)動機(jī)平均協(xié)方差矩陣，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，cn為第n臺航空發(fā)動機(jī)的協(xié)方差矩陣；λi,i＝1,2,...,k為的第i個特征值，k表示原始監(jiān)測數(shù)據(jù)的維數(shù)，為的單位特征向量矩陣，其中uj為第j個特征值對應(yīng)的單位特征向量；按照特征值從大至小排序，選取r個特征值λ1≥λ2≥...≥λr,r<k記作r個主分量，則u＝[u1,u2,...,ur]^t為變換矩陣，其中r為主分量數(shù)目；原始監(jiān)測數(shù)據(jù)投影到新的特征空間，在獲取更低維度特征變量的同時(shí)，能夠最大限度保留原始監(jiān)測數(shù)據(jù)中包含的有效信息，消除噪聲和不相關(guān)變量的影響。第n臺航空發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)轉(zhuǎn)換公式如下：

其中u＝[u1,u2,...,ur]^t為變換矩陣，r表示所選取主分量的數(shù)目，ur為第r個主分量對應(yīng)的單位特征向量，xn為第n臺航空發(fā)動機(jī)原始監(jiān)測數(shù)據(jù)，xni為第n臺航空發(fā)動機(jī)第i次采樣的原始監(jiān)測數(shù)據(jù)，為基于共同主成分分析降維后的第n臺航空發(fā)動機(jī)監(jiān)測數(shù)據(jù)，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)，為第n臺航空發(fā)動機(jī)降維后第mn次采樣的監(jiān)測數(shù)據(jù)，(·)^t表示轉(zhuǎn)置。

本實(shí)施例中令r＝6，保留前6個主分量，能量所占比例達(dá)到94％以上，各主分量取值如圖3所示，證明所選取的主分量個數(shù)足夠保留原始監(jiān)測數(shù)據(jù)中的信息。對發(fā)動機(jī)訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)中的原始監(jiān)測數(shù)據(jù)均進(jìn)行轉(zhuǎn)換。

1.2、在1.1的基礎(chǔ)上，進(jìn)行發(fā)動機(jī)健康參數(shù)的統(tǒng)計(jì)。研究表明發(fā)動機(jī)的前5％壽命周期可以被視作沒有發(fā)生退化，故提取訓(xùn)練集中的未退化監(jiān)測數(shù)據(jù)視作發(fā)動機(jī)的總體健康參數(shù)h＝(h1,h2,...,hn)^t，其中為第n臺發(fā)動機(jī)的健康參數(shù)，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)，sn,n＝1,2,...,n表示訓(xùn)練集中第n臺航空發(fā)動機(jī)健康參數(shù)的最大采樣次數(shù)，計(jì)算如下：

sn＝[0.05mn](4)

其中sn表示訓(xùn)練集中第n臺航空發(fā)動機(jī)健康參數(shù)的最大采樣次數(shù)，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)，[·]表示向下取整，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)；

1.3、獲取發(fā)動機(jī)的健康數(shù)據(jù)后，首先計(jì)算航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值。航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值m表示如下：

其中m為航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值，sn表示訓(xùn)練集中第n臺航空發(fā)動機(jī)健康參數(shù)的最大采樣次數(shù)，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，∑表示求和；

m的物理含義是指發(fā)動機(jī)在未發(fā)生退化時(shí)的基準(zhǔn)監(jiān)測數(shù)據(jù)。發(fā)動機(jī)的健康指標(biāo)可通過當(dāng)前時(shí)刻監(jiān)測數(shù)據(jù)與發(fā)動機(jī)未退化時(shí)基準(zhǔn)健康數(shù)據(jù)的ed確定。每臺發(fā)動機(jī)的健康指標(biāo)計(jì)算公式如下：

其中tj為第n臺發(fā)動機(jī)第j次采樣的時(shí)間，yn(tj)為第j次采樣時(shí)第n臺發(fā)動機(jī)健康指標(biāo)，m為航空發(fā)動機(jī)未退化監(jiān)測數(shù)據(jù)均值，為第n臺航空發(fā)動機(jī)降維后第j次采樣的監(jiān)測數(shù)據(jù)；故第n臺發(fā)動機(jī)健康指標(biāo)其中mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)。

以#1、#3號訓(xùn)練發(fā)動機(jī)為例，圖4a和圖4b是基于ed的發(fā)動機(jī)健康指標(biāo)曲線，觀察可知發(fā)動機(jī)總體性能呈現(xiàn)下降趨勢，在運(yùn)行初期退化率較小，后期退化速度加快。由此可證明航空發(fā)動機(jī)退化軌跡帶有非線性漂移。

1.4、基于發(fā)動機(jī)健康指標(biāo)，確定發(fā)動機(jī)失效閾值。根據(jù)研究表明，航空發(fā)動機(jī)退化過程如圖5所示，被劃分為四個階段，最后一階段可視作發(fā)動機(jī)失效階段。本研究基于k‐means聚類算法，聚類數(shù)目設(shè)置為4，計(jì)算最后一類的聚類中心作為發(fā)動機(jī)失效閾值ω。具體實(shí)施步驟如下：

第一步，對于第n臺發(fā)動機(jī)健康指標(biāo)每次采樣的健康指標(biāo)視作一個樣本。

第二步，隨機(jī)選擇4個初始聚類中心z1(k),z2(k)…z4(k)，k表示迭代次數(shù)，初始值為1，每迭代一次加1；

第三步，逐個將yn(tp),p＝1,2,...,mn按最小距離準(zhǔn)則分配到4個中心為zi(k)的類中。

第三步，分計(jì)算各類的新中心：nj即第j類中的樣本個數(shù)，oj表示屬于第j類的樣本集合。

第四步，如果zj(k+1)≠zj(k)則轉(zhuǎn)到第2步，將模式重新分類，重復(fù)迭代計(jì)算，直到聚類中心收斂，記錄4個聚類中心的最大值ωmax(n)。

第五步，對于每一臺發(fā)動機(jī)，重復(fù)第1、2、3、4步，選取ωmax(n)中的最小值，視作發(fā)動機(jī)失效閾值，即ω＝min(ωmax(n)),n＝1,2,...,100,其中min(·)表示取最小值。

本實(shí)施例中經(jīng)計(jì)算航空發(fā)動機(jī)的失效閾值ω＝1.27。

2、航空發(fā)動機(jī)退化建模及參數(shù)估計(jì)

2.1、航空發(fā)動機(jī)的退化過程中，初期退化并不顯著，后期由于損耗失效導(dǎo)致退化速度加快，健康指標(biāo)呈現(xiàn)增大趨勢并且整個過程不可逆。因此，選擇帶非線性漂移的維納過程描述渦扇發(fā)動機(jī)的退化過程。對于假設(shè)在采樣時(shí)間t，健康指標(biāo)為y(t)，失效閾值為ω，基于維納過程的航空發(fā)動機(jī)退化過程可以描述為下式：

y(t)＝y(tǒng)(0)+at^b+σb(t)(7)

其中，t為采樣時(shí)間，非線性漂移部分為at^b，擴(kuò)散系數(shù)為σ，b(t)為布朗運(yùn)動；為表征制造偏差和運(yùn)行環(huán)境導(dǎo)致發(fā)動機(jī)個體之間的差異性，a為非線性漂移部分的隨機(jī)變量，滿足其中n(·)表示正態(tài)分布，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差；同時(shí)將b定義為非線性漂移部分的常量表征發(fā)動機(jī)的相同特性；忽略發(fā)動機(jī)原始磨損，對于未運(yùn)行的發(fā)動機(jī)而言，假設(shè)其初始退化狀態(tài)為0，則有y(0)＝0。

發(fā)動機(jī)在某一時(shí)刻健康指標(biāo)等于或者超過失效失效閾值的時(shí)候，系統(tǒng)必須停止運(yùn)行。故航空發(fā)動機(jī)的壽命可以視作維納過程的首次穿越失效閾值的時(shí)間，定義如下：

t＝inf{t:y(t)≥ω|y(0)<ω}(8)

其中：t為發(fā)動機(jī)壽命周期，t為采樣時(shí)間，y(t)為t時(shí)刻健康指標(biāo)，ω為失效閾值，inf{·}表示下確界；

2.2、基于最大似然法的航空發(fā)動機(jī)退化模型參數(shù)估計(jì)。在建立渦扇發(fā)動機(jī)的退化模型后，需要基于發(fā)動機(jī)健康指標(biāo)對模型進(jìn)行參數(shù)估計(jì)。每臺發(fā)動機(jī)在多個采樣點(diǎn)的健康指標(biāo)服從多維正態(tài)分布，故第n臺發(fā)動機(jī)的健康指標(biāo)為：

其中μn為第n臺發(fā)動機(jī)的健康指標(biāo)均值向量，a為非線性漂移部分的隨機(jī)變量，b為非線性漂移部分的常量，tj為在第j次采樣的時(shí)間，e(·)表示期望，(·)^t表示轉(zhuǎn)置；本實(shí)施例中tj＝j(luò)。

其中σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，cov(·)表示協(xié)方差，σ為擴(kuò)散系數(shù)，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，b為非線性漂移部分的常量，tj為第n臺發(fā)動機(jī)第j次采樣的時(shí)間，(·)^t表示轉(zhuǎn)置；

對于不同發(fā)動機(jī)而言，退化過程獨(dú)立同分布，因此退化模型的似然函數(shù)為:

其中l(wèi)(μa,σa,σ,b|y)表示關(guān)于總體發(fā)動機(jī)健康指標(biāo)y的似然函數(shù)，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，σ為擴(kuò)散系數(shù)，b為非線性漂移部分的常量，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)，yn為第n臺發(fā)動機(jī)的健康指標(biāo)，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，μn為第n臺發(fā)動機(jī)的健康指標(biāo)均值向量，|·|表示取模操作，為σn的逆矩陣，exp(·)表示以e為底的指數(shù)函數(shù)，∏表示求積；

其中μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，σ為擴(kuò)散系數(shù)，b為非線性漂移部分的常量，y為總體發(fā)動機(jī)健康指標(biāo)，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，mn為第n臺發(fā)動機(jī)的最大采樣次數(shù)，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，|·|表示取模操作，為σn的逆矩陣，yn為第n臺發(fā)動機(jī)的健康指標(biāo)，μn為第n臺發(fā)動機(jī)的健康指標(biāo)均值向量；

通過最大化似然函數(shù)進(jìn)行參數(shù)估計(jì)。

其中tj為第n臺發(fā)動機(jī)第j次采樣的時(shí)間，表示偏導(dǎo)數(shù)；

其中a為非線性漂移部分的隨機(jī)變量，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，σ為擴(kuò)散系數(shù)，b為非線性漂移部分的常量，為μa的估計(jì)值，n為訓(xùn)練集中發(fā)動機(jī)臺數(shù)，tj為第n臺發(fā)動機(jī)第j次采樣的時(shí)間，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，σn為第n臺發(fā)動機(jī)的健康指標(biāo)協(xié)方差矩陣，|·|表示取模操作，為σn的逆矩陣，yn為第n臺發(fā)動機(jī)的健康指標(biāo)，y為總體發(fā)動機(jī)健康指標(biāo)；

基于單純形法尋找式(15)的最小值，繼而獲得μa估計(jì)值。

本實(shí)施例最終參數(shù)估計(jì)結(jié)果如表1所示：

表1航空發(fā)動機(jī)退化模型參數(shù)

3、建立航空發(fā)動機(jī)剩余壽命預(yù)測模型：

3.1、基于維納過程的航空發(fā)動機(jī)剩余壽命描述?；诎l(fā)動機(jī)性能退化模型具有時(shí)間平移不變性，則退化模型可以被表示為：

其中tc為航空發(fā)動機(jī)當(dāng)前運(yùn)行時(shí)間，li為航空發(fā)動機(jī)剩余壽命，b(·)表示布朗運(yùn)動，y(tc)為航空發(fā)動機(jī)當(dāng)前健康指標(biāo)，y(li+tc)為航空發(fā)動機(jī)失效時(shí)的健康指標(biāo)，a為非線性漂移部分的隨機(jī)變量，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)，y(li)表示發(fā)動機(jī)隨航空發(fā)動機(jī)剩余壽命li退化的維納過程，∫表示積分；可知發(fā)動機(jī)剩余壽命描述仍可表征為帶非線性漂移的維納過程；此時(shí)維納過程首次穿越ω-y(tc)的時(shí)間即為發(fā)動機(jī)的剩余壽命。

3.2、根據(jù)3.1可知航空發(fā)動機(jī)的剩余壽命模型可用帶非線性漂移的維納過程進(jìn)行描述。為估計(jì)發(fā)動機(jī)的剩余壽命，首先應(yīng)將帶非線性漂移的維納過程轉(zhuǎn)換為布朗運(yùn)動。航空發(fā)動機(jī)剩余壽命描述轉(zhuǎn)換為標(biāo)準(zhǔn)維納過程，對于非線性擴(kuò)散過程，若滿足：

其中μ(y,li)表示非線性漂移部分，σ′(li,u)為擴(kuò)散部分，y為健康指標(biāo)，li為航空發(fā)動機(jī)剩余壽命，c1和c2為任意關(guān)于時(shí)間的函數(shù)，z為初始退化量，∫表示積分，表示偏導(dǎo)數(shù)；

則可經(jīng)公式(19)將非線性擴(kuò)散過程轉(zhuǎn)換布朗運(yùn)動。

其中y為航空發(fā)動機(jī)健康指標(biāo)，為轉(zhuǎn)換后的航空發(fā)動機(jī)健康指標(biāo)，li為航空發(fā)動機(jī)剩余壽命，為轉(zhuǎn)換后的剩余壽命，ki，(i＝1,2,3)為任意系數(shù)且k1>0，t0為初始時(shí)間，ti,(i＝0,1,2)為初始時(shí)間且滿足ti≥0，z為初始退化量，c1和c2為任意關(guān)于時(shí)間的函數(shù)；

對于發(fā)動機(jī)剩余壽命模型而言，公式(18)中的擴(kuò)散部分σ(y,li)＝σ²，漂移部分μ(y,li)＝ab(li+tc)^b-1。當(dāng)航空發(fā)動機(jī)剩余壽命li＝0時(shí)，y(li)＝0，據(jù)此可知初始時(shí)間t0＝0，維納過程初始值z＝0。因此，可得下式：

其中a為非線性漂移部分的隨機(jī)變量，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)，tc為航空發(fā)動機(jī)當(dāng)前運(yùn)行時(shí)間，li為航空發(fā)動機(jī)剩余壽命；令k1＝1,ki＝0(i＝2,3)，tj＝0(j＝1,2)，則轉(zhuǎn)換過程表示如下：

其中y為航空發(fā)動機(jī)健康指標(biāo)，為轉(zhuǎn)換后的航空發(fā)動機(jī)健康指標(biāo)，li為航空發(fā)動機(jī)剩余壽命，為轉(zhuǎn)換后的剩余壽命，a為非線性漂移部分的隨機(jī)變量，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)；

轉(zhuǎn)換后的失效閾值：其中為ω發(fā)動機(jī)失效閾值；

4、基于性能衰退的航空發(fā)動機(jī)剩余壽命預(yù)測：

建立航空發(fā)動機(jī)剩余壽命概率密度函數(shù)模型，在步驟1、2、3的基礎(chǔ)上，進(jìn)行實(shí)時(shí)在線壽命預(yù)測。

布朗運(yùn)動首次穿過失效閾值的概率密度函數(shù)表示如下：

pb(t)(ω,t)表示布朗運(yùn)動首次穿過閾值的概率密度函數(shù)，ω為失效閾值，t為采樣時(shí)間；

因此航空發(fā)動機(jī)剩余壽命概率密度函數(shù)為：

其中：px(li)(ω,li)表示退化過程首次穿過失效閾值的概率密度函數(shù)，li為航空發(fā)動機(jī)剩余壽命，為轉(zhuǎn)換后的剩余壽命，ω為發(fā)動機(jī)失效閾值，為轉(zhuǎn)換后的失效閾值，tc為航空發(fā)動機(jī)當(dāng)前運(yùn)行時(shí)間，y(tc)為航空發(fā)動機(jī)當(dāng)前健康指標(biāo)，a為非線性漂移部分的隨機(jī)變量，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)，表示偏導(dǎo)數(shù)；以上結(jié)果滿足航空發(fā)動機(jī)壽命預(yù)測的實(shí)際使用需要。

由于非線性漂移部分的隨機(jī)變量則a的概率密度函數(shù)為：

其中f(a)表示非線性漂移部分的隨機(jī)變量a的概率密度函數(shù)，μa和分別為隨機(jī)變量a的均值和方差；

根據(jù)全概率公式，剩余壽命概率密度函數(shù)模型可以被修正為：

其中ω為失效閾值，li為航空發(fā)動機(jī)剩余壽命，tc為航空發(fā)動機(jī)當(dāng)前運(yùn)行時(shí)間，y(tc)為航空發(fā)動機(jī)當(dāng)前健康指標(biāo)，a為非線性漂移部分的隨機(jī)變量，μa為隨機(jī)變量a的均值，σa為隨機(jī)變量a的標(biāo)準(zhǔn)差，f(a)為隨機(jī)變量a的概率密度函數(shù)，b為非線性漂移部分的常量，σ為擴(kuò)散系數(shù)，exp(·)表示以e為底的指數(shù)函數(shù)；

發(fā)動機(jī)預(yù)測剩余壽命計(jì)算如下：

其中，li為航空發(fā)動機(jī)剩余壽命，表示退化過程首次穿過閾值的概率密度函數(shù)；

圖6按照已運(yùn)行周期大小隨機(jī)選擇#1，#5，#10，#26，#60五臺測試發(fā)動機(jī)，繪制其剩余壽命概率密度函數(shù)并且計(jì)算出預(yù)測剩余壽命，以及真實(shí)剩余壽命的對比。

發(fā)動機(jī)剩余壽命預(yù)測準(zhǔn)確度的評價(jià)指標(biāo)包括及時(shí)預(yù)測數(shù)目、確定系數(shù)r2以及懲罰得分s。剩余壽命誤差區(qū)間如圖7所示，預(yù)測誤差其中rul為發(fā)動機(jī)真實(shí)剩余壽命，為發(fā)動機(jī)預(yù)測剩余壽命；確定系數(shù)的計(jì)算方法如下式：

其中：ruli為第i臺發(fā)動機(jī)真實(shí)剩余壽命，為第i臺發(fā)動機(jī)預(yù)測剩余壽命，為全部發(fā)動機(jī)真實(shí)剩余壽命的平均值；

懲罰得分的計(jì)算方法如下式：

其中e為預(yù)測誤差，為預(yù)測誤差為負(fù)值的發(fā)動機(jī)數(shù)目，為預(yù)測誤差為正值的發(fā)動機(jī)數(shù)目；

確定系數(shù)r越接近1表明預(yù)測模型性能越好，懲罰函數(shù)值越低說明預(yù)測模型的性能越好。

本實(shí)施例測試機(jī)中100臺發(fā)動機(jī)剩余壽命頻數(shù)分布及懲罰函數(shù)得分如圖8和圖9所示，表2為剩余壽命預(yù)測結(jié)果。

表2預(yù)測性能指標(biāo)結(jié)果

通過表2中的結(jié)果可知，本發(fā)明提出的方法可以提高航空發(fā)動機(jī)剩余壽命預(yù)測的準(zhǔn)確度，考慮不相關(guān)監(jiān)測數(shù)據(jù)的影響提高模型預(yù)測性能。

本發(fā)明與現(xiàn)有方法的最大區(qū)別之處在于對多維監(jiān)測數(shù)據(jù)進(jìn)行共同主成分提取，消除對于壽命預(yù)測冗余的特征；基于多維監(jiān)測數(shù)據(jù)進(jìn)行融合取代單變量預(yù)測；同時(shí)，引入帶非線性的維納過程描述發(fā)動機(jī)退化率非線性的特點(diǎn)，得到航空發(fā)動機(jī)剩余壽命預(yù)測的概率分布，進(jìn)行發(fā)動機(jī)剩余壽命點(diǎn)估計(jì)。實(shí)驗(yàn)結(jié)果驗(yàn)證了本方法的實(shí)用性。