一種基于多階段信息融合的航空發(fā)動機剩余壽命預測方法
【專利摘要】本發(fā)明基于多階段信息融合的航空發(fā)動機剩余壽命預測方法包括步驟：多源監(jiān)測參數(shù)去噪處理與特征提??；對多源監(jiān)測時間序列進行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例；對多源參數(shù)進行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進行樣本訓練，獲得融合模型在多階段內(nèi)的參數(shù)；結(jié)合訓練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI；利用Kalman濾波算法，對發(fā)動機從性能完好到性能失效全過程進行最佳擬合，并最小化預測模型的誤差；結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI；利用Kalman濾波算法，對預測模型時變參數(shù)進行實時估計；確定預測模型，引入時間機制，實時估計發(fā)動機的失效時間。
【專利說明】
一種基于多階段信息融合的航空發(fā)動機剩余壽命預測方法
【技術(shù)領(lǐng)域】
:
[0001]本發(fā)明涉及發(fā)動機剩余壽命預測方法，特別涉及一種基于多階段信息融合的航空發(fā)動機剩余壽命預測方法。
【背景技術(shù)】
:
[0002]作為飛機的核心部件，航空發(fā)動機的結(jié)構(gòu)極其復雜，因此，單一監(jiān)測參數(shù)并不能準確表征其性能狀況。為綜合利用觀測信息以準確描述發(fā)動機的實時健康狀況，信息融合技術(shù)在發(fā)動機的健康管理中得到了廣泛的應用，例如基于信息融合的故障診斷、性能評估以及性能趨勢預測等。其中，衰退趨勢的預測對發(fā)動機的預知維修決策尤為重要，是目前計劃維修轉(zhuǎn)向預知維修過程中的研究熱點。
[0003]信息融合技術(shù)是用數(shù)學方法和技術(shù)工具綜合不同源信息，從而得到高品質(zhì)的有用信息。因此，存在各種不同種類、不同等級的融合，如數(shù)據(jù)融合、圖像融合、特征融合、決策融合、傳感器融合、分類器融合等。隨著計算機技術(shù)和網(wǎng)絡通信技術(shù)的飛速發(fā)展，信息融合的應用領(lǐng)域不斷擴展，從開始誕生的軍事領(lǐng)域，逐漸向其它領(lǐng)域滲透，如:智能機器人與智能車輛領(lǐng)域、醫(yī)學圖象處理與診斷、氣象預報、地球科學、農(nóng)業(yè)應用領(lǐng)域、現(xiàn)代制造領(lǐng)域和經(jīng)濟商業(yè)領(lǐng)域等。在發(fā)動機的健康管理中，信息融合技術(shù)顯得尤為重要，其功能表現(xiàn)在:充分利用信息的冗余性與互補性，提高時間或空間的分辨率，增加目標特征矢量的維數(shù)，降低信息的不確定性，改善信息的置信度，從而降低推理的模糊程度，提高系統(tǒng)決策能力。
[0004]在電子設備、電氣設備、控制設備、航空發(fā)動機等復雜設備健康管理中，以健康狀態(tài)預測為核心的健康管理技術(shù)是復雜裝備健康管理的核心技術(shù)。對民航發(fā)動機的健康狀態(tài)預測主要是對發(fā)動機性能下降階段的健康狀態(tài)預測。對民航發(fā)動機衰退趨勢的準確預測，可以制定更加有效的維修計劃，能夠在保證發(fā)動機可用性和可靠性的前提下有效降低運行成本，實現(xiàn)安全與經(jīng)濟的雙贏。從方法體系的角度可以將現(xiàn)有的民航發(fā)動機健康狀態(tài)變化趨勢預測方法劃分為模型驅(qū)動和數(shù)據(jù)驅(qū)動兩種方法。前者往往需要基于發(fā)動機工作原理、材料屬性、工況等因素建立解析模型，其建模代價一般比較高。后者則無需先驗假設，可以直接基于發(fā)動機的健康狀態(tài)數(shù)據(jù)進行建模，其建模代價相對較小。民航發(fā)動機狀態(tài)監(jiān)視技術(shù)的發(fā)展為數(shù)據(jù)驅(qū)動建模提供了豐富的數(shù)據(jù)，使得數(shù)據(jù)驅(qū)動建模預測方法具有良好的應用前景。
[0005]常見的基于融合進行發(fā)動機狀態(tài)預測的過程大多沒有考慮發(fā)動機各個參數(shù)與發(fā)動機整體性能衰退的突變性。因此，已有研究并未對發(fā)動機的衰退過程做分階段處理，從而在融合中不能準確反映發(fā)動機的性能衰退非平穩(wěn)性，預測模型的樣本訓練也會因為單階段處理弓丨起參數(shù)估計不準確的問題。
【發(fā)明內(nèi)容】
:
[0006]為解決常用發(fā)動機剩余壽命預測方法所不能解決的問題，本發(fā)明對現(xiàn)有發(fā)動機健康狀態(tài)預測技術(shù)進行改進，實現(xiàn)多源數(shù)據(jù)的分階段融合，基于多階段信息融合結(jié)果建立預測模型并利用多階段訓練樣本估計模型的參數(shù)，引入實時監(jiān)測數(shù)據(jù)，采用智能算法實現(xiàn)多階段模型參數(shù)的更新與預測，最終預測發(fā)動機的性能衰退趨勢，得到剩余使用壽命的準確估計。
[0007]本發(fā)明采用如下技術(shù)方案:一種基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其包括如下步驟:
[0008](a)多源監(jiān)測參數(shù)去噪處理與特征提取；
[0009](b)在步驟(a)的基礎(chǔ)上，對多源監(jiān)測時間序列進行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例；
[0010](C)在步驟(b)的基礎(chǔ)上，對多源參數(shù)進行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進行樣本訓練，獲得融合模型在多階段內(nèi)的參數(shù)；
[0011](d)在步驟(C)的基礎(chǔ)上，結(jié)合訓練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI ；
[0012](e)在步驟(d)的基礎(chǔ)上，利用Kalman濾波算法，對發(fā)動機從性能完好到性能失效全過程進行最佳擬合，以確定健康指標趨勢預測模型結(jié)構(gòu)，并最小化預測模型的誤差；
[0013](f)在步驟(C)的基礎(chǔ)上，結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI，此時的HI是預測對象；
[0014](g)在步驟(f)的基礎(chǔ)上，利用Kalman濾波算法，對預測模型的時變參數(shù)進行實時估計；
[0015](h)在步驟(e)和步驟(g)的基礎(chǔ)上，確定預測模型，引入時間機制，從而實時估計發(fā)動機的失效時間。
[0016]本發(fā)明具有如下有益效果:本發(fā)明所采用的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法是一種數(shù)據(jù)驅(qū)動方法，從發(fā)動機的多源監(jiān)測參數(shù)中獲得退化軌跡，而不需要對發(fā)動機性能衰退趨勢的形態(tài)進行假設。多階段的融合過程不僅能夠?qū)崿F(xiàn)數(shù)據(jù)的融合而且能夠?qū)崿F(xiàn)特征的融合，從而切實表征發(fā)動機性能退化的非平穩(wěn)性或突變性。因此，本發(fā)明能夠利用較為全面的監(jiān)測信息，基于發(fā)動機的多階段性能退化特征，借助人工智能的算法獲得滿足要求的預測結(jié)果，具備實際應用價值。
【專利附圖】
【附圖說明】
:
[0017]圖1為本發(fā)明的結(jié)構(gòu)框架示意圖。
[0018]圖2為本發(fā)明研究實體在全壽命周期內(nèi)的低壓壓縮機出口總溫(T24)監(jiān)測時間序列。
[0019]圖3為發(fā)動機多階段融合步驟示意圖。
[0020]圖4為發(fā)動機全壽命周期內(nèi)的健康指標HI的退化趨勢。
[0021]圖5為發(fā)動機壽命預測步驟示意圖。
[0022]圖6為單階段融合與多階段融合的結(jié)果比較示意圖。
[0023]圖7為訓練集下融合得到的HI時間序列。
[0024]圖8為訓練集下的HI預測過程。
【具體實施方式】
:
[0025]請參照圖1所示，本發(fā)明基于多階段信息融合的航空發(fā)動機剩余壽命預測方法主要包含信息融合與衰退趨勢預測兩大模塊。信息融合模塊包括:數(shù)據(jù)處理，用于去除工頻噪聲、背景噪聲以及隨機脈沖干擾等因素對純凈信號的影響；平穩(wěn)性分析，用于分析監(jiān)測時間序列的平穩(wěn)性，并準確定位引起性能衰退非平穩(wěn)的突變點；階段劃分與特征提取，用于對非平穩(wěn)的監(jiān)測時間序列進行階段劃分，實現(xiàn)階段內(nèi)的平穩(wěn)性，從而達到對時間序列的階段線性化處理；多階段分析，用于多源數(shù)據(jù)間的關(guān)聯(lián)分析，建立多源的、具有階段平穩(wěn)性的監(jiān)測數(shù)據(jù)與發(fā)動機隱含的健康狀況間的對應關(guān)系。衰退趨勢預測包含:模型建立，在多階段信息融合的基礎(chǔ)上，采用智能算法基于歷史性能數(shù)據(jù)構(gòu)建性能衰退趨勢預測模型；多階段參數(shù)估計，采用分階段樣本訓練的方法估計預測模型的參數(shù)；參數(shù)更新與預測，結(jié)合實時監(jiān)測數(shù)據(jù)，采用智能算法對模型參數(shù)進行更新與預測，最后計算發(fā)動機剩余使用壽命。
[0026]請參照圖1至圖5所示，本發(fā)明基于多階段信息融合的航空發(fā)動機剩余壽命預測方法的具體步驟如下:
[0027](a)多源監(jiān)測參數(shù)去噪處理與特征提取。
[0028](b)在步驟(a)的基礎(chǔ)上，對多源監(jiān)測時間序列進行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例。
[0029](C)在步驟(b)的基礎(chǔ)上，對多源參數(shù)進行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進行樣本訓練，獲得融合模型在多階段內(nèi)的參數(shù)。
[0030](d)在步驟(C)的基礎(chǔ)上，結(jié)合訓練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI。
[0031](e)在步驟(d)的基礎(chǔ)上，利用Kalman濾波算法，對發(fā)動機從性能完好到性能失效全過程進行最佳擬合，以確定健康指標趨勢預測模型結(jié)構(gòu)，并最小化預測模型的誤差。
[0032](f)在步驟(C)的基礎(chǔ)上，結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI，此時的HI是預測對象。
[0033](g)在步驟(f)的基礎(chǔ)上，利用Kalman濾波算法，對預測模型的時變參數(shù)進行實時估計。
[0034](h)在步驟(e)和步驟(g)的基礎(chǔ)上，確定預測模型，引入時間機制，從而實時估計發(fā)動機的失效時間。
[0035]步驟(a)中，進行多源監(jiān)測參數(shù)的去噪處理。
[0036]由于多源監(jiān)測數(shù)據(jù)具有非平穩(wěn)非線性的特征，采用經(jīng)驗模態(tài)分解方法(EmpiricalMode Decomposit1n, EMD)對多源監(jiān)測數(shù)據(jù)進行去噪處理，EDM將原始監(jiān)測數(shù)據(jù)分解為從高頻到低頻的多個本征模函數(shù)(Intrinsic Mode Funct1n, IMF)與殘余項之和。其中，IMF還包括信號起主導作用的模態(tài)與噪聲起主導作用的模態(tài)，只需從MF中剔除噪聲起主導作用的模態(tài)，將信號起主導作用的模態(tài)與殘余項進行疊加，即可進行部分重構(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)所有極大值形成的上包絡線與所有極小值形成的下包絡線的均值始終為O。EMD分解步驟:
[0042](I)找出z(t)的所有局部極大值，對該極值進行三次樣條插值，得到由所有極大值點構(gòu)成的上包絡線，記為a(t)；
[0043](2)找出z(t)的所有局部極小值，對該極值進行三次樣條插值，得到由所有極小值點構(gòu)成的下包絡線，記為b(t)；
[0044](3)計算上下包絡線的均值記為c (t)，c{t) = a(t)+2b{t)；
[0045](4)計算d(t) = z(t)_c(t),當d(t)滿足上述IMF的兩條性質(zhì)時,d(t)即第一個MF分量；否則，以d(t)作為輸入量，重復步驟(I)?(3)，直到第一個MF產(chǎn)生，記第一個IMF分量為IMF1⑴；
[0046](5)計算rjt) = zahMFjOdfna)作為輸入量，重復步驟⑴?(4)，得到第二個IMF分量IMF2 (t);計算1*2(1:) = z (t)-1MF2(t),重復上述步驟,直到余項rn(t)為一個單調(diào)序列或小于給定閾值，EMD分解即結(jié)束。
[0047]分解結(jié)束之后，需根據(jù)最小連續(xù)均方誤差(consecutive mean squareerror, CMSE)的準則區(qū)分MF各個分量中哪些部分對信號起主導作用，哪些部分對噪聲其主導作用。CMSE計算公式如下:
[0048]CMSE Qk,zk+l)?)]2 =去￡[.Λ.)]2⑵
/V tN
[0049]其中，知為讓項相加重構(gòu)出的去噪監(jiān)測時間序列，包含k_l個IMF分量和一個余項。當CMSE的值最小時，說明在k處MF分量出現(xiàn)轉(zhuǎn)折，而MF首次轉(zhuǎn)折的地方即為噪聲起主導作用模態(tài)與信號起主導作用模態(tài)的分界，去噪之后的監(jiān)測時間序列由以下公式進行部分重構(gòu):頂Fk+1 (t) +IMFk+2 (t) +...+IMFn^1 (t) +rn(t)。
[0050]步驟(b)中，對多源監(jiān)測時間序列進行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例。
[0051](I)采用以下啟發(fā)式算法對時間序列的平穩(wě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點左右兩部分的標準偏差；SD(i)為合并偏差；τα)為檢驗統(tǒng)計值，τα)越大，表明該點左右兩部份統(tǒng)計特性的差異越大。該算法通過計算τα)將非平穩(wěn)時間序列每個時間點的左右兩部分進行比較。τα)值最大時對應的i點即為該非平穩(wěn)時間序列的分割點，可將整個時間窗一分為二。當?shù)谝粋€分割點獲得之后，可用同樣的方法對劃分好的兩個時間窗分別進行分割即可獲得四個時間窗，以此類推。時間窗的數(shù)量應當依據(jù)實際應用情況和子窗內(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)測值在相應時間間隔內(nèi)的變化量。理論上，時間間隔越小越能反映詳細的慢變量變化過程，但是由于噪聲影響，該過程的波動幅度會隨著時間間隔的縮小而變大?？紤]到最終要獲取的是慢變量的相對變化(從一個穩(wěn)態(tài)到另一個穩(wěn)態(tài)的轉(zhuǎn)變)，可適當調(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)測時間序列的突變點不嚴格在一個點上，但是由于監(jiān)測參數(shù)退化突變的相似性，例如，前半段退化緩慢后半段退化迅速，各個監(jiān)測參數(shù)的突變點會集中在某個點附近，因此將各個監(jiān)測時間序列突變點的平均值作為融合時間序列的突變點。這樣才能確保多源監(jiān)測時間序列的突變性與融合時間序列(HI時間序列)的突變性相一致。例如(基于圖4的數(shù)據(jù))，計算得到突變點在第217個循環(huán)時附近，退化量在整個退化過程中的比例是
0.507，說明在突變點附近的HI應為0.493。
[0064]步驟(c)中，對多源參數(shù)進行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進行樣本訓練，獲得融合模型在多階段內(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ù)需基于訓練樣本集和樣本數(shù)據(jù)通過多元線性回歸的方法獲得。假設前5個循環(huán)時對應的HI為I ;突變點附近的5個循環(huán)時對應的HI為kAn，kAn =1-Hiita ;后5個循環(huán)時對應的HI為O。訓練樣本集可表示為:
Xt kl
kAl
[0069]M=1:(9)
_XF k0 _
[0070]樣本集M中，XT、XY1...Xyn^P Xf分別表示HI為l、kA1...kAn和O時對應的性能參數(shù)監(jiān)測值樣本。以樣本訓練集M進行監(jiā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，尤代入回歸模型中，同時引入訓練集中發(fā)動機從性能完好到失效的全壽命周期多源監(jiān)測數(shù)據(jù)，可得到全壽命周期健康指標HI時間序列，如圖4所示?？傊?，多源監(jiān)測參數(shù)通過去噪處理產(chǎn)生的多組信號是通過多元回歸的方法進行融合得到一組新的退化時間序列，即HI時間序列。
[0074]步驟(e)中，利用Kalman濾波算法，對發(fā)動機從性能完好到性能失效全過程進行最佳擬合，以確定健康指標趨勢預測模型結(jié)構(gòu)，并最小化預測模型的誤差。
[0075](I)通過狀態(tài)空間方法建立健康指標預測模型:
[0076]xt = F.xt_1+wt (10)
[0077]yt=H +mt(11)
[0078]式中:xt為t時刻對應的狀態(tài)向量Xi =[^%?；\...為"];4，0?,3,...為"是預
1......0"
測模型的時變參數(shù)；兔為HI預測值；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
為觀測過程噪聲。基于以上描述，預測模型可表示為:
[0079]% 二 a0t+心…+ " +St(12)
[0080]其中，δ t表示預測誤差，m是預測模型的階數(shù)，決定了預測模型的結(jié)構(gòu)。δ t與m的值可根據(jù)最佳擬合(預測值與實際值的方差最小)確定，最佳擬合由Kalman濾波算法實現(xiàn)。另外，考慮到預測模型的時變參數(shù)《Κ?,...X隨時間變化的特性，采用具有實時遞歸運算能力的Kalman濾波算法對時變參數(shù)進行估計。
[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為預測的偏差值,4 = |八-火|，yt為實際值1力預測值。
[0085]其中時變參數(shù)估計與檢測時間序列最佳擬合的算法流程:
[0086]1.時變參數(shù)集狀態(tài)集為 X (t) = {x1； X2, , xn};
[0087]I〗.對于t = 1: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 )預測階段:根據(jù)t_l時刻的狀態(tài)與t_l時刻及t_l時刻前的觀測信息預測t時刻的狀態(tài):=F'^-l|r-l ;根據(jù)t_l時刻的狀態(tài)協(xié)方差與t_l時刻及t-ι時刻前的觀測信息預測t時刻的狀態(tài)協(xié)方差;Ρ?ΙΗ = F.Ρ?_ιΜ.FT+Qt,其中Ε[_χ,~毛卜1 ] =。，Ρ,\,-ι=cov(χ,—毛卜!);
[0090](iii)更新階段:結(jié)合卡爾曼增益更新狀態(tài)與協(xié)方差的值，作為下一步預測的基礎(chǔ)值；更新狀態(tài)更新協(xié)方差=PtIt = PtIt-1-KtHptlH= PtIt^1HVStjSt =
Hpt|t-1HT+Rt ;
[0091]II1.時變參數(shù)的預測:將上一步驟中所得到的I作為根據(jù)t時刻及以前的觀測信息與t-ι時刻的狀態(tài)信息的時變參數(shù)狀態(tài)預測值；引入時間因子t即可計算得到HI的預測值；
[0092]IV.最佳擬合的檢驗:計算不同預測模型階數(shù)下的HI預測誤差序列^ = |只-丸|，
將最小預測誤差對應的階數(shù)作為最終預測模型的階數(shù)，即確定A+4中的m
值，其中St取誤差均值。
[0093]需要注意的是，融合多源監(jiān)測參數(shù)獲得HI時間序列的方法是多元回歸，而Kalman濾波算法一方面用于訓練集下的HI最佳擬合與模型結(jié)構(gòu)(δ t與m的值)的確定，另一方面用于測試集下的HI預測模型時變參數(shù)的預測。
[0094]例如，最佳擬合的情況下，預測模型的結(jié)構(gòu)為丸=0?+4 + ?>2+^，其中δ =±0.138。此預測模型則可以用于測試集中的HI預測。
[0095]步驟(f)中，將測試集中的實時多源監(jiān)測數(shù)據(jù)代入回歸模型中，回歸模型的系數(shù)為步驟(C)中計算得到的an(l，4"通合多源監(jiān)測參數(shù)，得到測試集下的健康指標HI時間序列，此時的HI是需要預測的對象。
[0096]該步驟的實質(zhì)是將步驟(d)的訓練數(shù)據(jù)替換為測試數(shù)據(jù)，引入多階段融合模型中，即可得到當前發(fā)動機的實時HI時間序列。例如，計算得到訓練集下的HI時間序列如圖7所示。
[0097]步驟(g)中，利用Kalman濾波算法(詳見步驟(e)中⑵中的Kalman算法流程)，對預測模型的時變參數(shù)進行估計。將步驟(f)中得到的HI時間序列引入HI預測模型V, = Oi0 + a)t + aft2 + δ，δ = ±0.138。采用第步驟(e)中的Kalman算法流程,從而實現(xiàn)時變參數(shù)心aj, W的實時預測。
[0098]步驟(h)中，在確定的HI預測模型中引入時間機制，從而實時估計發(fā)動機的失效時間。結(jié)合步驟(e)中確定的模型結(jié)構(gòu)與步驟(g)中估計的預測模型時變參數(shù)，代入時間，預測HI的變化趨勢，最終可計算發(fā)動機剩余壽命。預測的HI等于O時對應的時間即為發(fā)動機性能失效的時間。預測結(jié)果如圖8所示。該發(fā)動機的實際失效時間是98個循環(huán)時之后，預測的剩余使用壽命是90，因此具有較好的預測精度。
[0099]當利用Kalman濾波算法對預測模型的時變參數(shù)進行預測時，基于以下對比可說明基于多階段信息融合的預測方法的優(yōu)越性:
[0100](I)在融合多源監(jiān)測參數(shù)計算HI時，如果不采用分階段的融合方法，則會出現(xiàn)如圖6虛線所示的情況。由于監(jiān)測參數(shù)不分階段處理，用于融合各個監(jiān)測參數(shù)的回歸模型的回歸系數(shù)在整個監(jiān)測時間窗內(nèi)是不變的，這將導致融合得到的HI時間序列不能準確反映發(fā)動機的實際退化特性。即使Kalman濾波方法適用于非平穩(wěn)數(shù)據(jù)的預測，并考慮的模型參數(shù)的時變特性，由于輸入的HI不準確，得到的剩余壽命預測結(jié)果也將是不準確的。
[0101](2)當融合多源監(jiān)測參數(shù)計算HI時，采用分階段的融合方法，則能達到圖6點劃線的預測效果。由于考慮的退化過程的突變特性，采用分階段的方法對多源監(jiān)測參數(shù)在不同時間窗分別進行融合，此時的回歸模型的系數(shù)在整個檢測時間窗是分階段的，因此計算得到的HI可以較為準確地反映發(fā)動機的實際退化過程。此時，由于Kalman濾波算法的輸入量HI較為準確，結(jié)合Kalman濾波自身的優(yōu)點(非平穩(wěn)數(shù)據(jù)的預測適用性以及對時變參數(shù)的預測更新能力)，預測得到的剩余使用壽命自然更接近真實值。
[0102]以上所述僅是本發(fā)明的優(yōu)選實施方式，應當指出，對于本【技術(shù)領(lǐng)域】的普通技術(shù)人員來說，在不脫離本發(fā)明原理的前提下還可以作出若干改進，這些改進也應視為本發(fā)明的保護范圍。
【權(quán)利要求】
1.一種基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:包括如下步驟 (a)多源監(jiān)測參數(shù)去噪處理與特征提?。? (b)在步驟(a)的基礎(chǔ)上，對多源監(jiān)測時間序列進行平穩(wěn)性分析，計算各個參數(shù)監(jiān)測時間序列突變點，計算突變點處參數(shù)退化比例； (C)在步驟(b)的基礎(chǔ)上，對多源參數(shù)進行多階段劃分，建立回歸融合模型，利用歷史監(jiān)測數(shù)據(jù)進行樣本訓練，獲得融合模型在多階段內(nèi)的參數(shù)； (d)在步驟(C)的基礎(chǔ)上，結(jié)合訓練集中的監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI ； (e)在步驟(d)的基礎(chǔ)上，利用Kalman濾波算法，對發(fā)動機從性能完好到性能失效全過程進行最佳擬合，以確定健康指標趨勢預測模型結(jié)構(gòu)，并最小化預測模型的誤差； (f)在步驟(C)的基礎(chǔ)上，結(jié)合測試集中的實時監(jiān)測數(shù)據(jù)，融合多源監(jiān)測參數(shù)，得到健康指標HI，此時的HI是預測對象； (g)在步驟(f)的基礎(chǔ)上，利用Kalman濾波算法，對預測模型的時變參數(shù)進行實時估計； (h)在步驟(e)和步驟(g)的基礎(chǔ)上，確定預測模型，引入時間機制，從而實時估計發(fā)動機的失效時間。
2.如權(quán)利要求1所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:步驟(a)中，多源監(jiān)測參數(shù)去噪處理與特征提取，采用經(jīng)驗模態(tài)分解方法EMD對多源監(jiān)測數(shù)據(jù)進行去噪處理，EDM將原始監(jiān)測數(shù)據(jù)分解為從高頻到低頻的多個本征模函數(shù)MF與殘差分量之和，頂F包括信號起主導作用的模態(tài)與噪聲起主導作用的模態(tài)，從MF中剔除噪聲起主導作用的模態(tài)，將信號起主導作用的模態(tài)與殘余項進行疊加，即可進行部分重構(gòu)得到降噪之后的監(jiān)測時間序列，EMD分解的公式如下:
其中，Z (t)為監(jiān)測時間序列，IMF1, IMF2,IMFn*別代表從高到低不同頻率段的信號成分，殘余項rn代表信號的平均趨勢。
3.如權(quán)利要求2所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:步驟(b)中包括如下: (1)采用以下啟發(fā)式算法對時間序列的平穩(wěn)性進行分析并劃分時間窗:
其中,NpN2分別表示i點左右兩部分的時間點個數(shù)^1Qhu2Q)分別表示i點左右兩部分的均值；Si(i)、S2⑴分別表示i點左右兩部分的標準偏差；SD(i)為合并偏差；T(i)為檢驗統(tǒng)計值； (2)將劃分時間窗的點作為時間序列的分割點，則需判斷分割點是否為性能突變點，計算慢變量:
其中，Sv為慢變量，Δ t為時間間隔，j > i, i, j = I, 2,..., η, η為時間點的個數(shù)，Δ S為監(jiān)測值在相應時間間隔內(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ā)動機剩余壽命預測方法，其特征在于:步驟(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ù)通過訓練樣本集和樣本數(shù)據(jù)通過多元線性回歸的方法獲得，假設前5個循環(huán)時對應的HI為I ;突變點附近的5個循環(huán)時對應的HI為kAn，= 1-1nto ;后5個循環(huán)時對應的HI為0，訓練樣本集可表示為:
樣本集M中，XT、XY1...Xyn和Xf分別表示HI為l、kA1...kAn和O時對應的性能參數(shù)監(jiān)測值樣本，以樣本訓練集M進行監(jiān)測樣本訓練求解系數(shù)an(l，ξ。
5.如權(quán)利要求4所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:步驟(d)中，將多階段的融合模型系數(shù)an(l，4代入回歸模型中，同時引入訓練集中發(fā)動機從性能完好到失效的全壽命周期多源監(jiān)測數(shù)據(jù)，得到全壽命周期健康指標HI時間序列。
6.如權(quán)利要求5所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:步驟(e)中包括如下: (1)通過狀態(tài)空間方法建立健康指標預測模型:
Xt = F * xt-1+wtV/ = I l.x{ + m , 式中:xt為t時刻對應的狀態(tài)向量Λ.,是預測模型的時變參數(shù)；SHI預測值；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為觀測過程噪聲，基于以上描述，預測模型表示為:yt = a^t +a)t.., + a'tnf + St 其中，\表示預測誤差，m是預測模型的階數(shù)，決定了預測模型的結(jié)構(gòu)，31與111的值可根據(jù)最佳擬合即預測值與實際值的方差最小確定，最佳擬合由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為預測的偏差值，1^=Iit-Ah Yt為實際值，兔為預測值。
7.如權(quán)利要求6所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:所述時變參數(shù)估計與時間序列最佳擬合檢測的算法流程包括: I.時變參數(shù)集狀態(tài)集為X (t) = {χ1； X2,, xj ； II.對于t = 1: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 (? )預測階段:根據(jù)t-1時刻的狀態(tài)與t-1時刻及t-1時刻前的觀測信息預測t時刻的狀態(tài):為M =F^t-1\t-1;根據(jù)t-1時刻的狀態(tài)協(xié)方差與t-Ι時刻及t-Ι時刻前的觀測信息預測 t 時刻的狀態(tài)協(xié)方差；pt|H = F.FT+Qt，其中五[A- ,J = Oj ^iM=Cov(X1-Xihl); (iii)更新階段:結(jié)合卡爾曼增益更新狀態(tài)與協(xié)方差的值，作為下一步預測的基礎(chǔ)值；更新狀態(tài):~=x+-1+ji:A;更新協(xié)方差-.Ptlt = Ptit-1-KtHptI^1 ;其中，Kt = Pt,^1HVSt, St =Hptlt-PRt ； II1.時變參數(shù)的預測:將上一步驟中所得到的％作為根據(jù)t時刻及以前的觀測信息與t-ι時刻的狀態(tài)信息的時變參數(shù)狀態(tài)預測值；引入時間因子t即可計算得到HI的預測值； IV.最佳擬合的檢驗:計算不同預測模型階數(shù)下的HI預測誤差序列K= !J't-Λ?，將最小預測誤差對應的階數(shù)作為最終預測模型的階數(shù)，即確定λ = <+在中的m值，其中St取誤差均值。
8.如權(quán)利要求7所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:所述步驟(f)中，將測試集中的實時多源監(jiān)測數(shù)據(jù)代入回歸模型中，回歸模型的系數(shù)為步驟(c)中計算得到的4『，融合多源監(jiān)測參數(shù)，得到測試集下的健康指標HI時間序列，此時的HI是需要預測的對象。
9.如權(quán)利要求8所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:所述步驟(g)中，利用Kalman濾波算法對預測模型的時變參數(shù)進行估計，將步驟(f)中得到的HI時間序列引入HI預測模型，采用步驟(e)中的Kalman算法流程實現(xiàn)時變參數(shù)的實時預測。
10.如權(quán)利要求9所述的基于多階段信息融合的航空發(fā)動機剩余壽命預測方法，其特征在于:所述步驟(h)中結(jié)合步驟(e)中確定的模型結(jié)構(gòu)與步驟(g)中的預測模型時變參數(shù)，代入時間，預測HI的變化趨勢，最終計算出發(fā)動機剩余壽命。
【文檔編號】G06F19/00GK104166787SQ201410341589
【公開日】2014年11月26日 申請日期:2014年7月17日 優(yōu)先權(quán)日:2014年7月17日
【發(fā)明者】劉君強, 張馬蘭, 左洪福, 謝吉偉 申請人:南京航空航天大學