一種基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法
【專利摘要】本發(fā)明基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法包括步驟：多源監(jiān)測參數(shù)去噪處理與特征提取；對多源監(jiān)測時間序列進(jìn)行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例；對多源參數(shù)進(jìn)行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進(jìn)行樣本訓(xùn)練，獲得融合模型在多階段內(nèi)的參數(shù)；結(jié)合訓(xùn)練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI；利用Kalman濾波算法，對發(fā)動機(jī)從性能完好到性能失效全過程進(jìn)行最佳擬合，并最小化預(yù)測模型的誤差；結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI；利用Kalman濾波算法，對預(yù)測模型時變參數(shù)進(jìn)行實時估計；確定預(yù)測模型，引入時間機(jī)制，實時估計發(fā)動機(jī)的失效時間。
【專利說明】
一種基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法

【技術(shù)領(lǐng)域】
:
[0001]本發(fā)明涉及發(fā)動機(jī)剩余壽命預(yù)測方法，特別涉及一種基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法。

【背景技術(shù)】
:
[0002]作為飛機(jī)的核心部件，航空發(fā)動機(jī)的結(jié)構(gòu)極其復(fù)雜，因此，單一監(jiān)測參數(shù)并不能準(zhǔn)確表征其性能狀況。為綜合利用觀測信息以準(zhǔn)確描述發(fā)動機(jī)的實時健康狀況，信息融合技術(shù)在發(fā)動機(jī)的健康管理中得到了廣泛的應(yīng)用，例如基于信息融合的故障診斷、性能評估以及性能趨勢預(yù)測等。其中，衰退趨勢的預(yù)測對發(fā)動機(jī)的預(yù)知維修決策尤為重要，是目前計劃維修轉(zhuǎn)向預(yù)知維修過程中的研究熱點。
[0003]信息融合技術(shù)是用數(shù)學(xué)方法和技術(shù)工具綜合不同源信息，從而得到高品質(zhì)的有用信息。因此，存在各種不同種類、不同等級的融合，如數(shù)據(jù)融合、圖像融合、特征融合、決策融合、傳感器融合、分類器融合等。隨著計算機(jī)技術(shù)和網(wǎng)絡(luò)通信技術(shù)的飛速發(fā)展，信息融合的應(yīng)用領(lǐng)域不斷擴(kuò)展，從開始誕生的軍事領(lǐng)域，逐漸向其它領(lǐng)域滲透，如:智能機(jī)器人與智能車輛領(lǐng)域、醫(yī)學(xué)圖象處理與診斷、氣象預(yù)報、地球科學(xué)、農(nóng)業(yè)應(yīng)用領(lǐng)域、現(xiàn)代制造領(lǐng)域和經(jīng)濟(jì)商業(yè)領(lǐng)域等。在發(fā)動機(jī)的健康管理中，信息融合技術(shù)顯得尤為重要，其功能表現(xiàn)在:充分利用信息的冗余性與互補(bǔ)性，提高時間或空間的分辨率，增加目標(biāo)特征矢量的維數(shù)，降低信息的不確定性，改善信息的置信度，從而降低推理的模糊程度，提高系統(tǒng)決策能力。
[0004]在電子設(shè)備、電氣設(shè)備、控制設(shè)備、航空發(fā)動機(jī)等復(fù)雜設(shè)備健康管理中，以健康狀態(tài)預(yù)測為核心的健康管理技術(shù)是復(fù)雜裝備健康管理的核心技術(shù)。對民航發(fā)動機(jī)的健康狀態(tài)預(yù)測主要是對發(fā)動機(jī)性能下降階段的健康狀態(tài)預(yù)測。對民航發(fā)動機(jī)衰退趨勢的準(zhǔn)確預(yù)測，可以制定更加有效的維修計劃，能夠在保證發(fā)動機(jī)可用性和可靠性的前提下有效降低運行成本，實現(xiàn)安全與經(jīng)濟(jì)的雙贏。從方法體系的角度可以將現(xiàn)有的民航發(fā)動機(jī)健康狀態(tài)變化趨勢預(yù)測方法劃分為模型驅(qū)動和數(shù)據(jù)驅(qū)動兩種方法。前者往往需要基于發(fā)動機(jī)工作原理、材料屬性、工況等因素建立解析模型，其建模代價一般比較高。后者則無需先驗假設(shè)，可以直接基于發(fā)動機(jī)的健康狀態(tài)數(shù)據(jù)進(jìn)行建模，其建模代價相對較小。民航發(fā)動機(jī)狀態(tài)監(jiān)視技術(shù)的發(fā)展為數(shù)據(jù)驅(qū)動建模提供了豐富的數(shù)據(jù)，使得數(shù)據(jù)驅(qū)動建模預(yù)測方法具有良好的應(yīng)用前景。
[0005]常見的基于融合進(jìn)行發(fā)動機(jī)狀態(tài)預(yù)測的過程大多沒有考慮發(fā)動機(jī)各個參數(shù)與發(fā)動機(jī)整體性能衰退的突變性。因此，已有研究并未對發(fā)動機(jī)的衰退過程做分階段處理，從而在融合中不能準(zhǔn)確反映發(fā)動機(jī)的性能衰退非平穩(wěn)性，預(yù)測模型的樣本訓(xùn)練也會因為單階段處理弓丨起參數(shù)估計不準(zhǔn)確的問題。


【發(fā)明內(nèi)容】

:
[0006]為解決常用發(fā)動機(jī)剩余壽命預(yù)測方法所不能解決的問題，本發(fā)明對現(xiàn)有發(fā)動機(jī)健康狀態(tài)預(yù)測技術(shù)進(jìn)行改進(jìn)，實現(xiàn)多源數(shù)據(jù)的分階段融合，基于多階段信息融合結(jié)果建立預(yù)測模型并利用多階段訓(xùn)練樣本估計模型的參數(shù)，引入實時監(jiān)測數(shù)據(jù)，采用智能算法實現(xiàn)多階段模型參數(shù)的更新與預(yù)測，最終預(yù)測發(fā)動機(jī)的性能衰退趨勢，得到剩余使用壽命的準(zhǔn)確估計。
[0007]本發(fā)明采用如下技術(shù)方案:一種基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其包括如下步驟:
[0008](a)多源監(jiān)測參數(shù)去噪處理與特征提??；
[0009](b)在步驟(a)的基礎(chǔ)上，對多源監(jiān)測時間序列進(jìn)行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例；
[0010](C)在步驟(b)的基礎(chǔ)上，對多源參數(shù)進(jìn)行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進(jìn)行樣本訓(xùn)練，獲得融合模型在多階段內(nèi)的參數(shù)；
[0011](d)在步驟(C)的基礎(chǔ)上，結(jié)合訓(xùn)練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI ；
[0012](e)在步驟(d)的基礎(chǔ)上，利用Kalman濾波算法，對發(fā)動機(jī)從性能完好到性能失效全過程進(jìn)行最佳擬合，以確定健康指標(biāo)趨勢預(yù)測模型結(jié)構(gòu)，并最小化預(yù)測模型的誤差；
[0013](f)在步驟(C)的基礎(chǔ)上，結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI，此時的HI是預(yù)測對象；
[0014](g)在步驟(f)的基礎(chǔ)上，利用Kalman濾波算法，對預(yù)測模型的時變參數(shù)進(jìn)行實時估計；
[0015](h)在步驟(e)和步驟(g)的基礎(chǔ)上，確定預(yù)測模型，引入時間機(jī)制，從而實時估計發(fā)動機(jī)的失效時間。
[0016]本發(fā)明具有如下有益效果:本發(fā)明所采用的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法是一種數(shù)據(jù)驅(qū)動方法，從發(fā)動機(jī)的多源監(jiān)測參數(shù)中獲得退化軌跡，而不需要對發(fā)動機(jī)性能衰退趨勢的形態(tài)進(jìn)行假設(shè)。多階段的融合過程不僅能夠?qū)崿F(xiàn)數(shù)據(jù)的融合而且能夠?qū)崿F(xiàn)特征的融合，從而切實表征發(fā)動機(jī)性能退化的非平穩(wěn)性或突變性。因此，本發(fā)明能夠利用較為全面的監(jiān)測信息，基于發(fā)動機(jī)的多階段性能退化特征，借助人工智能的算法獲得滿足要求的預(yù)測結(jié)果，具備實際應(yīng)用價值。

【專利附圖】

【附圖說明】
:
[0017]圖1為本發(fā)明的結(jié)構(gòu)框架示意圖。
[0018]圖2為本發(fā)明研究實體在全壽命周期內(nèi)的低壓壓縮機(jī)出口總溫(T24)監(jiān)測時間序列。
[0019]圖3為發(fā)動機(jī)多階段融合步驟示意圖。
[0020]圖4為發(fā)動機(jī)全壽命周期內(nèi)的健康指標(biāo)HI的退化趨勢。
[0021]圖5為發(fā)動機(jī)壽命預(yù)測步驟示意圖。
[0022]圖6為單階段融合與多階段融合的結(jié)果比較示意圖。
[0023]圖7為訓(xùn)練集下融合得到的HI時間序列。
[0024]圖8為訓(xùn)練集下的HI預(yù)測過程。

【具體實施方式】
:
[0025]請參照圖1所示，本發(fā)明基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法主要包含信息融合與衰退趨勢預(yù)測兩大模塊。信息融合模塊包括:數(shù)據(jù)處理，用于去除工頻噪聲、背景噪聲以及隨機(jī)脈沖干擾等因素對純凈信號的影響；平穩(wěn)性分析，用于分析監(jiān)測時間序列的平穩(wěn)性，并準(zhǔn)確定位引起性能衰退非平穩(wěn)的突變點；階段劃分與特征提取，用于對非平穩(wěn)的監(jiān)測時間序列進(jìn)行階段劃分，實現(xiàn)階段內(nèi)的平穩(wěn)性，從而達(dá)到對時間序列的階段線性化處理；多階段分析，用于多源數(shù)據(jù)間的關(guān)聯(lián)分析，建立多源的、具有階段平穩(wěn)性的監(jiān)測數(shù)據(jù)與發(fā)動機(jī)隱含的健康狀況間的對應(yīng)關(guān)系。衰退趨勢預(yù)測包含:模型建立，在多階段信息融合的基礎(chǔ)上，采用智能算法基于歷史性能數(shù)據(jù)構(gòu)建性能衰退趨勢預(yù)測模型；多階段參數(shù)估計，采用分階段樣本訓(xùn)練的方法估計預(yù)測模型的參數(shù)；參數(shù)更新與預(yù)測，結(jié)合實時監(jiān)測數(shù)據(jù)，采用智能算法對模型參數(shù)進(jìn)行更新與預(yù)測，最后計算發(fā)動機(jī)剩余使用壽命。
[0026]請參照圖1至圖5所示，本發(fā)明基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法的具體步驟如下:
[0027](a)多源監(jiān)測參數(shù)去噪處理與特征提取。
[0028](b)在步驟(a)的基礎(chǔ)上，對多源監(jiān)測時間序列進(jìn)行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例。
[0029](C)在步驟(b)的基礎(chǔ)上，對多源參數(shù)進(jìn)行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進(jìn)行樣本訓(xùn)練，獲得融合模型在多階段內(nèi)的參數(shù)。
[0030](d)在步驟(C)的基礎(chǔ)上，結(jié)合訓(xùn)練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI。
[0031](e)在步驟(d)的基礎(chǔ)上，利用Kalman濾波算法，對發(fā)動機(jī)從性能完好到性能失效全過程進(jìn)行最佳擬合，以確定健康指標(biāo)趨勢預(yù)測模型結(jié)構(gòu)，并最小化預(yù)測模型的誤差。
[0032](f)在步驟(C)的基礎(chǔ)上，結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI，此時的HI是預(yù)測對象。
[0033](g)在步驟(f)的基礎(chǔ)上，利用Kalman濾波算法，對預(yù)測模型的時變參數(shù)進(jìn)行實時估計。
[0034](h)在步驟(e)和步驟(g)的基礎(chǔ)上，確定預(yù)測模型，引入時間機(jī)制，從而實時估計發(fā)動機(jī)的失效時間。
[0035]步驟(a)中，進(jìn)行多源監(jiān)測參數(shù)的去噪處理。
[0036]由于多源監(jiān)測數(shù)據(jù)具有非平穩(wěn)非線性的特征，采用經(jīng)驗?zāi)B(tài)分解方法(EmpiricalMode Decomposit1n, EMD)對多源監(jiān)測數(shù)據(jù)進(jìn)行去噪處理，EDM將原始監(jiān)測數(shù)據(jù)分解為從高頻到低頻的多個本征模函數(shù)(Intrinsic Mode Funct1n, IMF)與殘余項之和。其中，IMF還包括信號起主導(dǎo)作用的模態(tài)與噪聲起主導(dǎo)作用的模態(tài)，只需從MF中剔除噪聲起主導(dǎo)作用的模態(tài)，將信號起主導(dǎo)作用的模態(tài)與殘余項進(jìn)行疊加，即可進(jìn)行部分重構(gòu)得到降噪之后的監(jiān)測時間序列。EMD分解的公式如下:

η
[0037]z(0 = X/.V/f；.+Γη(I)

/-1
[0038]其中，ζ (t)為監(jiān)測時間序列，頂F1, IMF2, MFn分別代表從高到低不同頻率段的信號成分，殘余項^代表信號的平均趨勢。
[0039]分解過程中，IMF必須滿足以下條件:
[0040](I)在整個時間范圍內(nèi)，極值點與過零點的數(shù)目相等或最多相差一個；
[0041](2)所有極大值形成的上包絡(luò)線與所有極小值形成的下包絡(luò)線的均值始終為O。EMD分解步驟:
[0042](I)找出z(t)的所有局部極大值，對該極值進(jìn)行三次樣條插值，得到由所有極大值點構(gòu)成的上包絡(luò)線，記為a(t)；
[0043](2)找出z(t)的所有局部極小值，對該極值進(jìn)行三次樣條插值，得到由所有極小值點構(gòu)成的下包絡(luò)線，記為b(t)；
[0044](3)計算上下包絡(luò)線的均值記為c (t)，c{t) = a(t)+2b{t)；
[0045](4)計算d(t) = z(t)_c(t),當(dāng)d(t)滿足上述IMF的兩條性質(zhì)時,d(t)即第一個MF分量；否則，以d(t)作為輸入量，重復(fù)步驟(I)?(3)，直到第一個MF產(chǎn)生，記第一個IMF分量為IMF1⑴；
[0046](5)計算rjt) = zahMFjOdfna)作為輸入量，重復(fù)步驟⑴?(4)，得到第二個IMF分量IMF2 (t);計算1*2(1:) = z (t)-1MF2(t),重復(fù)上述步驟,直到余項rn(t)為一個單調(diào)序列或小于給定閾值，EMD分解即結(jié)束。
[0047]分解結(jié)束之后，需根據(jù)最小連續(xù)均方誤差(consecutive mean squareerror, CMSE)的準(zhǔn)則區(qū)分MF各個分量中哪些部分對信號起主導(dǎo)作用，哪些部分對噪聲其主導(dǎo)作用。CMSE計算公式如下:
[0048]CMSE Qk,zk+l)?)]2 =去￡[.Λ.)]2⑵
/V tN
[0049]其中，知為讓項相加重構(gòu)出的去噪監(jiān)測時間序列，包含k_l個IMF分量和一個余項。當(dāng)CMSE的值最小時，說明在k處MF分量出現(xiàn)轉(zhuǎn)折，而MF首次轉(zhuǎn)折的地方即為噪聲起主導(dǎo)作用模態(tài)與信號起主導(dǎo)作用模態(tài)的分界，去噪之后的監(jiān)測時間序列由以下公式進(jìn)行部分重構(gòu):頂Fk+1 (t) +IMFk+2 (t) +...+IMFn^1 (t) +rn(t)。
[0050]步驟(b)中，對多源監(jiān)測時間序列進(jìn)行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例。
[0051](I)采用以下啟發(fā)式算法對時間序列的平穩(wěn)性進(jìn)行分析并劃分時間窗:
一I
,_[(^V1-1)X-1)XS2Ofl1「I , 1-
0ΓΛ = - X--1--
r ?"N, +Ν,-2N, iY9,，、
[0052]L12」L.? 2」.()'W)
[0053]其中，Np N2分別表示i點左右兩部分的時間點個數(shù)W1Q)、U2 (i)分別表示i點左右兩部分的均值A(chǔ)ahS2Q)分別表示i點左右兩部分的標(biāo)準(zhǔn)偏差；SD(i)為合并偏差；τα)為檢驗統(tǒng)計值，τα)越大，表明該點左右兩部份統(tǒng)計特性的差異越大。該算法通過計算τα)將非平穩(wěn)時間序列每個時間點的左右兩部分進(jìn)行比較。τα)值最大時對應(yīng)的i點即為該非平穩(wěn)時間序列的分割點，可將整個時間窗一分為二。當(dāng)?shù)谝粋€分割點獲得之后，可用同樣的方法對劃分好的兩個時間窗分別進(jìn)行分割即可獲得四個時間窗，以此類推。時間窗的數(shù)量應(yīng)當(dāng)依據(jù)實際應(yīng)用情況和子窗內(nèi)時間序列的平穩(wěn)性而定。
[0054](2)劃分時間窗的點可作為時間序列的分割點，但要判斷分割點是否為性能突變點，還需要計算慢變量:
c Δ5 Sj-Si
[0055]~= T-14)

At tj - ?-
[0056]其中，Sv為慢變量，Δ t為時間間隔，j > i, i, j = I, 2,..., η, η為時間點的個數(shù)，
ΔS為監(jiān)測值在相應(yīng)時間間隔內(nèi)的變化量。理論上，時間間隔越小越能反映詳細(xì)的慢變量變化過程，但是由于噪聲影響，該過程的波動幅度會隨著時間間隔的縮小而變大?？紤]到最終要獲取的是慢變量的相對變化(從一個穩(wěn)態(tài)到另一個穩(wěn)態(tài)的轉(zhuǎn)變)，可適當(dāng)調(diào)整時間間隔。通過時間窗口 At的滑動可求得整個時間窗內(nèi)的慢變量序列。慢變量的突變即體現(xiàn)參數(shù)退化的突變點。在圖4給出的時間序列中，分割點有3個，但通過計算慢變量可確定其第二個分割點為該時間序列的實際突變點。
[0057](3)各個監(jiān)測參數(shù)時間序列的突變點確定后，即可計算監(jiān)測參數(shù)平均突變點與平均變化比例:

I 1
[0058]ρα,,=?ΥΑ,(5)

1 i=l

1、二
[0059]m λ?=-Z^m in(6)

1 i=l
[0060]其中，Pita為第n個平均突變點；mAn為第η個平均突變點處的平均變化比例；Pin為第i個監(jiān)測參數(shù)的第η個突變點，min為第i個監(jiān)測參數(shù)從退化初到第η個突變點的變化量在整個退化過程中的變化量的比例，i = 1,2,...,1,1為監(jiān)測參數(shù)的個數(shù)，η為突變點的個數(shù)。min的公式如下所示:
I 5^丨 η 12j jI Λ_
minX SiJ)-- X SiJ)
「λλ/^?)?=1^ t=n-2 / ) Z=I^ t=N-4,-X
_1]5奸2/ 5⑴
Σ Md Σ U
? 二Ir^n-2 I ? 二It 二 Ν-4
[0062]其中，f是第i個監(jiān)測參數(shù)在退化初期的5個循環(huán)時的監(jiān)測值的平均值，

-5 ?=1
I u.?'2I Λ
7 Σ Sij是第i個監(jiān)測參數(shù)在第η個突變點附近5個循環(huán)時的監(jiān)測值的平均值，-JSi t
J t=n-2D i=W_4
是第i個監(jiān)測參數(shù)在性能退化末期的5個循環(huán)時的監(jiān)測值的平均值，N是第i個監(jiān)測參數(shù)的監(jiān)測值個數(shù)。
[0063]多組監(jiān)測時間序列的突變點不嚴(yán)格在一個點上，但是由于監(jiān)測參數(shù)退化突變的相似性，例如，前半段退化緩慢后半段退化迅速，各個監(jiān)測參數(shù)的突變點會集中在某個點附近，因此將各個監(jiān)測時間序列突變點的平均值作為融合時間序列的突變點。這樣才能確保多源監(jiān)測時間序列的突變性與融合時間序列(HI時間序列)的突變性相一致。例如(基于圖4的數(shù)據(jù))，計算得到突變點在第217個循環(huán)時附近，退化量在整個退化過程中的比例是
0.507，說明在突變點附近的HI應(yīng)為0.493。
[0064]步驟(c)中，對多源參數(shù)進(jìn)行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進(jìn)行樣本訓(xùn)練，獲得融合模型在多階段內(nèi)的參數(shù)。
[0065](I)基于步驟(b)獲得的突變點，建立分階段融合模型:
vl — + XA^
[0066],-v2 = fl20+X4C8；
yn = aM +
[0067]其中，yl表示第I個突變點前的HI ;y2表示第I個突變點到第2個突變點間的HI ；yn表示第η-1個突變點到第η個突變點間的HI ；Χ表示多源監(jiān)測參數(shù)U1, x2，...，X1)，I表示性能參數(shù)的個數(shù)；a1(l, 4表示第I個突變點前的模型系數(shù)A1 = (an, a12,...，an) ；a20, Αζ表示第I個突變點到第2個突變點間的模型系數(shù)，A2 = (a21, a22, a21) ；an0, <表示第n_l個突變點到第n個突變點間的模型系數(shù)，An = (anl, an2, anl)。
[0068](2)融合模型的系數(shù)需基于訓(xùn)練樣本集和樣本數(shù)據(jù)通過多元線性回歸的方法獲得。假設(shè)前5個循環(huán)時對應(yīng)的HI為I ;突變點附近的5個循環(huán)時對應(yīng)的HI為kAn，kAn =1-Hiita ;后5個循環(huán)時對應(yīng)的HI為O。訓(xùn)練樣本集可表示為:
Xt kl

kAl
[0069]M=1:(9)
_XF k0 _
[0070]樣本集M中，XT、XY1...Xyn^P Xf分別表示HI為l、kA1...kAn和O時對應(yīng)的性能參數(shù)監(jiān)測值樣本。以樣本訓(xùn)練集M進(jìn)行監(jiān)測樣本訓(xùn)練求解系數(shù)an(l，4。例如，計算得到融合模型的系數(shù)為:
[0071][a10, aj = [335.848，—0.117，—0.003，—0.013，0.784，0.128，-0.148，0.003]
[0072][a20, a2i] = [0, -0.035,0.001, -0.020,38.560,0.037, -0.339,0.001]
[0073]步驟(d)中，將多階段的融合模型系數(shù)an(l，尤代入回歸模型中，同時引入訓(xùn)練集中發(fā)動機(jī)從性能完好到失效的全壽命周期多源監(jiān)測數(shù)據(jù)，可得到全壽命周期健康指標(biāo)HI時間序列，如圖4所示。總之，多源監(jiān)測參數(shù)通過去噪處理產(chǎn)生的多組信號是通過多元回歸的方法進(jìn)行融合得到一組新的退化時間序列，即HI時間序列。
[0074]步驟(e)中，利用Kalman濾波算法，對發(fā)動機(jī)從性能完好到性能失效全過程進(jìn)行最佳擬合，以確定健康指標(biāo)趨勢預(yù)測模型結(jié)構(gòu)，并最小化預(yù)測模型的誤差。
[0075](I)通過狀態(tài)空間方法建立健康指標(biāo)預(yù)測模型:
[0076]xt = F.xt_1+wt (10)
[0077]yt=H +mt(11)
[0078]式中:xt為t時刻對應(yīng)的狀態(tài)向量Xi =[^%?；\...為"];4，0?,3,...為"是預(yù)



1......0"
測模型的時變參數(shù)；兔為HI預(yù)測值；F為狀態(tài)轉(zhuǎn)移矩陣F= i J 0 ； ; H為觀測矩陣



O*...*.1-1

LJ/2+1
H = ^t°,tl,t2,t3,t4,,.,,tn^ ； Wt ?N(0, Qt), mt ?N(0, Rt), E[wtXmt = 0] ；wt 為狀態(tài)向量噪聲；mt
為觀測過程噪聲?；谝陨厦枋觯A(yù)測模型可表示為:
[0079]% 二 a0t+心…+ " +St(12)
[0080]其中，δ t表示預(yù)測誤差，m是預(yù)測模型的階數(shù)，決定了預(yù)測模型的結(jié)構(gòu)。δ t與m的值可根據(jù)最佳擬合(預(yù)測值與實際值的方差最小)確定，最佳擬合由Kalman濾波算法實現(xiàn)。另外，考慮到預(yù)測模型的時變參數(shù)《Κ?,...X隨時間變化的特性，采用具有實時遞歸運算能力的Kalman濾波算法對時變參數(shù)進(jìn)行估計。
[0081](2)利用以下Kalman濾波算法公式實現(xiàn)時變參數(shù)估計與最佳擬合:
? 之卜I =F名—1卜I,..
[0082]TCJ3)
I^l i=F-Pt-1\t-1-F +Qt
? ^t\t = %—ι + KtSt
[0083]\C14)
\p,r Pt^l-KtHp^l
[0084]其中，χφ-ι為根據(jù)t_l時刻及該時刻之前的監(jiān)測信息估計t時刻的后驗狀態(tài)，￡[-T/--Vi] = 0，Xt為t時刻實際的后驗狀態(tài);ptlH為狀態(tài)估計誤差的協(xié)方差，
，Xt為t時刻的實際狀態(tài)；Kt表示卡爾曼增益，Kt = PtIt^1HVst, St =
HptlHH1+!^ vt為預(yù)測的偏差值,4 = |八-火|，yt為實際值1力預(yù)測值。
[0085]其中時變參數(shù)估計與檢測時間序列最佳擬合的算法流程:
[0086]1.時變參數(shù)集狀態(tài)集為 X (t) = {x1； X2, , xn};
[0087]I〗.對于t = 1:n，進(jìn)行以下操作:
[0088]( i )建立狀態(tài)方程xt = F.xt_!+wt,描述t-Ι時刻的狀態(tài)與t時刻狀態(tài)的關(guān)系；建立觀測方程A = + ，描述t時刻的狀態(tài)與t時刻測量值的關(guān)系；其中
= [“? ?'u;, β卜…，& ]，



I …...0
丑= [iV,iV/,...,i"]，wt ?N(0，Qt)，mt ?N(0，Rt)，E[wtXmt = O],及F=: Q...；;
-0......1-H+1
[0089]( ii )預(yù)測階段:根據(jù)t_l時刻的狀態(tài)與t_l時刻及t_l時刻前的觀測信息預(yù)測t時刻的狀態(tài):=F'^-l|r-l ;根據(jù)t_l時刻的狀態(tài)協(xié)方差與t_l時刻及t-ι時刻前的觀測信息預(yù)測t時刻的狀態(tài)協(xié)方差;Ρ?ΙΗ = F.Ρ?_ιΜ.FT+Qt,其中Ε[_χ,~毛卜1 ] =。，Ρ,\,-ι=cov(χ,—毛卜!);
[0090](iii)更新階段:結(jié)合卡爾曼增益更新狀態(tài)與協(xié)方差的值，作為下一步預(yù)測的基礎(chǔ)值；更新狀態(tài)更新協(xié)方差=PtIt = PtIt-1-KtHptlH= PtIt^1HVStjSt =
Hpt|t-1HT+Rt ;
[0091]II1.時變參數(shù)的預(yù)測:將上一步驟中所得到的I作為根據(jù)t時刻及以前的觀測信息與t-ι時刻的狀態(tài)信息的時變參數(shù)狀態(tài)預(yù)測值；引入時間因子t即可計算得到HI的預(yù)測值；
[0092]IV.最佳擬合的檢驗:計算不同預(yù)測模型階數(shù)下的HI預(yù)測誤差序列^ = |只-丸|，
將最小預(yù)測誤差對應(yīng)的階數(shù)作為最終預(yù)測模型的階數(shù)，即確定A+4中的m
值，其中St取誤差均值。
[0093]需要注意的是，融合多源監(jiān)測參數(shù)獲得HI時間序列的方法是多元回歸，而Kalman濾波算法一方面用于訓(xùn)練集下的HI最佳擬合與模型結(jié)構(gòu)(δ t與m的值)的確定，另一方面用于測試集下的HI預(yù)測模型時變參數(shù)的預(yù)測。
[0094]例如，最佳擬合的情況下，預(yù)測模型的結(jié)構(gòu)為丸=0?+4 + ?>2+^，其中δ =±0.138。此預(yù)測模型則可以用于測試集中的HI預(yù)測。
[0095]步驟(f)中，將測試集中的實時多源監(jiān)測數(shù)據(jù)代入回歸模型中，回歸模型的系數(shù)為步驟(C)中計算得到的an(l，4"通合多源監(jiān)測參數(shù)，得到測試集下的健康指標(biāo)HI時間序列，此時的HI是需要預(yù)測的對象。
[0096]該步驟的實質(zhì)是將步驟(d)的訓(xùn)練數(shù)據(jù)替換為測試數(shù)據(jù)，引入多階段融合模型中，即可得到當(dāng)前發(fā)動機(jī)的實時HI時間序列。例如，計算得到訓(xùn)練集下的HI時間序列如圖7所示。
[0097]步驟(g)中，利用Kalman濾波算法(詳見步驟(e)中⑵中的Kalman算法流程)，對預(yù)測模型的時變參數(shù)進(jìn)行估計。將步驟(f)中得到的HI時間序列引入HI預(yù)測模型V, = Oi0 + a)t + aft2 + δ，δ = ±0.138。采用第步驟(e)中的Kalman算法流程,從而實現(xiàn)時變參數(shù)心aj, W的實時預(yù)測。
[0098]步驟(h)中，在確定的HI預(yù)測模型中引入時間機(jī)制，從而實時估計發(fā)動機(jī)的失效時間。結(jié)合步驟(e)中確定的模型結(jié)構(gòu)與步驟(g)中估計的預(yù)測模型時變參數(shù)，代入時間，預(yù)測HI的變化趨勢，最終可計算發(fā)動機(jī)剩余壽命。預(yù)測的HI等于O時對應(yīng)的時間即為發(fā)動機(jī)性能失效的時間。預(yù)測結(jié)果如圖8所示。該發(fā)動機(jī)的實際失效時間是98個循環(huán)時之后，預(yù)測的剩余使用壽命是90，因此具有較好的預(yù)測精度。
[0099]當(dāng)利用Kalman濾波算法對預(yù)測模型的時變參數(shù)進(jìn)行預(yù)測時，基于以下對比可說明基于多階段信息融合的預(yù)測方法的優(yōu)越性:
[0100](I)在融合多源監(jiān)測參數(shù)計算HI時，如果不采用分階段的融合方法，則會出現(xiàn)如圖6虛線所示的情況。由于監(jiān)測參數(shù)不分階段處理，用于融合各個監(jiān)測參數(shù)的回歸模型的回歸系數(shù)在整個監(jiān)測時間窗內(nèi)是不變的，這將導(dǎo)致融合得到的HI時間序列不能準(zhǔn)確反映發(fā)動機(jī)的實際退化特性。即使Kalman濾波方法適用于非平穩(wěn)數(shù)據(jù)的預(yù)測，并考慮的模型參數(shù)的時變特性，由于輸入的HI不準(zhǔn)確，得到的剩余壽命預(yù)測結(jié)果也將是不準(zhǔn)確的。
[0101](2)當(dāng)融合多源監(jiān)測參數(shù)計算HI時，采用分階段的融合方法，則能達(dá)到圖6點劃線的預(yù)測效果。由于考慮的退化過程的突變特性，采用分階段的方法對多源監(jiān)測參數(shù)在不同時間窗分別進(jìn)行融合，此時的回歸模型的系數(shù)在整個檢測時間窗是分階段的，因此計算得到的HI可以較為準(zhǔn)確地反映發(fā)動機(jī)的實際退化過程。此時，由于Kalman濾波算法的輸入量HI較為準(zhǔn)確，結(jié)合Kalman濾波自身的優(yōu)點(非平穩(wěn)數(shù)據(jù)的預(yù)測適用性以及對時變參數(shù)的預(yù)測更新能力)，預(yù)測得到的剩余使用壽命自然更接近真實值。
[0102]以上所述僅是本發(fā)明的優(yōu)選實施方式，應(yīng)當(dāng)指出，對于本【技術(shù)領(lǐng)域】的普通技術(shù)人員來說，在不脫離本發(fā)明原理的前提下還可以作出若干改進(jìn)，這些改進(jìn)也應(yīng)視為本發(fā)明的保護(hù)范圍。
【權(quán)利要求】
1.一種基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:包括如下步驟 (a)多源監(jiān)測參數(shù)去噪處理與特征提?。? (b)在步驟(a)的基礎(chǔ)上，對多源監(jiān)測時間序列進(jìn)行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例； (C)在步驟(b)的基礎(chǔ)上，對多源參數(shù)進(jìn)行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進(jìn)行樣本訓(xùn)練，獲得融合模型在多階段內(nèi)的參數(shù)； (d)在步驟(C)的基礎(chǔ)上，結(jié)合訓(xùn)練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI ； (e)在步驟(d)的基礎(chǔ)上，利用Kalman濾波算法，對發(fā)動機(jī)從性能完好到性能失效全過程進(jìn)行最佳擬合，以確定健康指標(biāo)趨勢預(yù)測模型結(jié)構(gòu)，并最小化預(yù)測模型的誤差； (f)在步驟(C)的基礎(chǔ)上，結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標(biāo)HI，此時的HI是預(yù)測對象； (g)在步驟(f)的基礎(chǔ)上，利用Kalman濾波算法，對預(yù)測模型的時變參數(shù)進(jìn)行實時估計； (h)在步驟(e)和步驟(g)的基礎(chǔ)上，確定預(yù)測模型，引入時間機(jī)制，從而實時估計發(fā)動機(jī)的失效時間。
2.如權(quán)利要求1所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:步驟(a)中，多源監(jiān)測參數(shù)去噪處理與特征提取，采用經(jīng)驗?zāi)B(tài)分解方法EMD對多源監(jiān)測數(shù)據(jù)進(jìn)行去噪處理，EDM將原始監(jiān)測數(shù)據(jù)分解為從高頻到低頻的多個本征模函數(shù)MF與殘差分量之和，頂F包括信號起主導(dǎo)作用的模態(tài)與噪聲起主導(dǎo)作用的模態(tài)，從MF中剔除噪聲起主導(dǎo)作用的模態(tài)，將信號起主導(dǎo)作用的模態(tài)與殘余項進(jìn)行疊加，即可進(jìn)行部分重構(gòu)得到降噪之后的監(jiān)測時間序列，EMD分解的公式如下:
其中，Z (t)為監(jiān)測時間序列，IMF1, IMF2,IMFn*別代表從高到低不同頻率段的信號成分，殘余項rn代表信號的平均趨勢。
3.如權(quán)利要求2所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:步驟(b)中包括如下: (1)采用以下啟發(fā)式算法對時間序列的平穩(wěn)性進(jìn)行分析并劃分時間窗:
其中,NpN2分別表示i點左右兩部分的時間點個數(shù)^1Qhu2Q)分別表示i點左右兩部分的均值；Si(i)、S2⑴分別表示i點左右兩部分的標(biāo)準(zhǔn)偏差；SD(i)為合并偏差；T(i)為檢驗統(tǒng)計值； (2)將劃分時間窗的點作為時間序列的分割點，則需判斷分割點是否為性能突變點，計算慢變量:
其中，Sv為慢變量，Δ t為時間間隔，j > i, i, j = I, 2,..., η, η為時間點的個數(shù)，Δ S為監(jiān)測值在相應(yīng)時間間隔內(nèi)的變化量，通過時間窗口 At的滑動求得整個時間窗內(nèi)的慢變量序列，慢變量的突變即體現(xiàn)參數(shù)退化的突變點； (3)各個監(jiān)測參數(shù)時間序列的突變點確定后，即計算監(jiān)測參數(shù)平均突變點與平均變化比例:
其中，Pto為第η個平均突變點:1?為第η個平均突變點處的平均變化比例；Pin為第i個個監(jiān)測參數(shù)的第η個突變點，min為第i個監(jiān)測參數(shù)從退化初到第η個突變點的變化量在整個退化過程中的變化量的比例，i = 1,2,...,1,1為監(jiān)測參數(shù)的個數(shù)，η為突變點的個數(shù)。
4.如權(quán)利要求3所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:步驟(C)中包括如下: (1)基于步驟(b)中獲得的突變點，建立分階段融合模型:
其中，yl表示第I個突變點前的HI ；y2表示第I個突變點到第2個突變點間的HI ；yn表示第η-1個突變點到第η個突變點間的HI ；Χ表示多源監(jiān)測參數(shù)Oq，x2,..., X1), I表示性能參數(shù)的個數(shù)；a1(l, 4表示第I個突變點前的模型系數(shù)A1 = (an, a12,...，an) ；a20,為1'表示第I個突變點到第2個突變點間的模型系數(shù)，A2 = (a21, a22, a21) ；an0,尤表示第n_l個突變點到第η個突變點間的模型系數(shù)，An = (anl, an2, anl)； (2)模型系數(shù)通過訓(xùn)練樣本集和樣本數(shù)據(jù)通過多元線性回歸的方法獲得，假設(shè)前5個循環(huán)時對應(yīng)的HI為I ;突變點附近的5個循環(huán)時對應(yīng)的HI為kAn，= 1-1nto ;后5個循環(huán)時對應(yīng)的HI為0，訓(xùn)練樣本集可表示為:
樣本集M中，XT、XY1...Xyn和Xf分別表示HI為l、kA1...kAn和O時對應(yīng)的性能參數(shù)監(jiān)測值樣本，以樣本訓(xùn)練集M進(jìn)行監(jiān)測樣本訓(xùn)練求解系數(shù)an(l，ξ。
5.如權(quán)利要求4所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:步驟(d)中，將多階段的融合模型系數(shù)an(l，4代入回歸模型中，同時引入訓(xùn)練集中發(fā)動機(jī)從性能完好到失效的全壽命周期多源監(jiān)測數(shù)據(jù)，得到全壽命周期健康指標(biāo)HI時間序列。
6.如權(quán)利要求5所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:步驟(e)中包括如下: (1)通過狀態(tài)空間方法建立健康指標(biāo)預(yù)測模型:
Xt = F * xt-1+wtV/ = I l.x{ + m , 式中:xt為t時刻對應(yīng)的狀態(tài)向量Λ.,是預(yù)測模型的時變參數(shù)；SHI預(yù)測值；F為狀態(tài)轉(zhuǎn)移矩陣F= | J ° ； ; H為觀測矩陣H =


-0......1-?+I[t°, t1，t2，t3, t4，…，tn] ；wt ?N (0, Qt)，mt ?N (0, Rt)，E [wt X mt = 0] ；wt 為狀態(tài)向量噪聲；mt為觀測過程噪聲，基于以上描述，預(yù)測模型表示為:yt = a^t +a)t.., + a'tnf + St 其中，\表示預(yù)測誤差，m是預(yù)測模型的階數(shù)，決定了預(yù)測模型的結(jié)構(gòu)，31與111的值可根據(jù)最佳擬合即預(yù)測值與實際值的方差最小確定，最佳擬合由Kalman濾波算法實現(xiàn)； (2)利用以下Kalman濾波算法公式實現(xiàn)時變參數(shù)估計與最佳擬合:\ρ\ ! =ρ.ρ,-ι\,-ι.ρτ +Qt? ^\,=^,\,-χ+Κ>δ, 其中，11 I為根據(jù)t-1時刻及該時刻之前的監(jiān)測信息估計t時刻的后驗狀態(tài)，4-1] = 0，Xt為t時刻實際的后驗狀態(tài)；PtlH為狀態(tài)估計誤差的協(xié)方差，^|/-1=cov(^-Vi)，xt為t時刻的實際狀態(tài)；Kt表示卡爾曼增益，Kt = PtIt^1HVSt, St =HptIwH^Rt, vt為預(yù)測的偏差值，1^=Iit-Ah Yt為實際值，兔為預(yù)測值。
7.如權(quán)利要求6所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:所述時變參數(shù)估計與時間序列最佳擬合檢測的算法流程包括: I.時變參數(shù)集狀態(tài)集為X (t) = {χ1； X2,, xj ； II.對于t = 1:n，進(jìn)行以下操作: (i )建立狀態(tài)方程xt = F.XH+wt，描述t-1時刻的狀態(tài)與t時刻狀態(tài)的關(guān)系；建立觀測方程免=〃4+叫，描述t時刻的狀態(tài)與t時刻測量值的關(guān)系；其中Xi =[0--為3,..,,<]， H = [t°，t1，t2，t3，t4，...，tn], wt ?N(0，Qt)，mt ?N(0，Rt)，E[wtXmt = 0]，及:1 0:F=；:0 ； 0I L......J?+l (? )預(yù)測階段:根據(jù)t-1時刻的狀態(tài)與t-1時刻及t-1時刻前的觀測信息預(yù)測t時刻的狀態(tài):為M =F^t-1\t-1;根據(jù)t-1時刻的狀態(tài)協(xié)方差與t-Ι時刻及t-Ι時刻前的觀測信息預(yù)測 t 時刻的狀態(tài)協(xié)方差；pt|H = F.FT+Qt，其中五[A- ,J = Oj ^iM=Cov(X1-Xihl); (iii)更新階段:結(jié)合卡爾曼增益更新狀態(tài)與協(xié)方差的值，作為下一步預(yù)測的基礎(chǔ)值；更新狀態(tài):~=x+-1+ji:A;更新協(xié)方差-.Ptlt = Ptit-1-KtHptI^1 ;其中，Kt = Pt,^1HVSt, St =Hptlt-PRt ； II1.時變參數(shù)的預(yù)測:將上一步驟中所得到的％作為根據(jù)t時刻及以前的觀測信息與t-ι時刻的狀態(tài)信息的時變參數(shù)狀態(tài)預(yù)測值；引入時間因子t即可計算得到HI的預(yù)測值； IV.最佳擬合的檢驗:計算不同預(yù)測模型階數(shù)下的HI預(yù)測誤差序列K= !J't-Λ?，將最小預(yù)測誤差對應(yīng)的階數(shù)作為最終預(yù)測模型的階數(shù)，即確定λ = <+在中的m值，其中St取誤差均值。
8.如權(quán)利要求7所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:所述步驟(f)中，將測試集中的實時多源監(jiān)測數(shù)據(jù)代入回歸模型中，回歸模型的系數(shù)為步驟(c)中計算得到的4『，融合多源監(jiān)測參數(shù)，得到測試集下的健康指標(biāo)HI時間序列，此時的HI是需要預(yù)測的對象。
9.如權(quán)利要求8所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:所述步驟(g)中，利用Kalman濾波算法對預(yù)測模型的時變參數(shù)進(jìn)行估計，將步驟(f)中得到的HI時間序列引入HI預(yù)測模型，采用步驟(e)中的Kalman算法流程實現(xiàn)時變參數(shù)的實時預(yù)測。
10.如權(quán)利要求9所述的基于多階段信息融合的航空發(fā)動機(jī)剩余壽命預(yù)測方法，其特征在于:所述步驟(h)中結(jié)合步驟(e)中確定的模型結(jié)構(gòu)與步驟(g)中的預(yù)測模型時變參數(shù)，代入時間，預(yù)測HI的變化趨勢，最終計算出發(fā)動機(jī)剩余壽命。
【文檔編號】G06F19/00GK104166787SQ201410341589
【公開日】2014年11月26日 申請日期:2014年7月17日 優(yōu)先權(quán)日:2014年7月17日
【發(fā)明者】劉君強(qiáng), 張馬蘭, 左洪福, 謝吉偉 申請人:南京航空航天大學(xué)