一種基于工況識別和相似性匹配的變工況下航空發動機整機剩余壽命預測方法
【專利摘要】本發明提出了一種基于工況識別和相似性匹配的變工況下航空發動機整機剩余壽命預測方法，該方法著力解決復雜的工況變化掩蓋航空發動機真實性能退化趨勢致使整機剩余壽命難于預測這個問題，包括如下步驟：首先，識別航空發動機的運行工況；其次，對不同工況下的歷史退化數據進行數據標準化；再次，對標準化后的數據進行傳感器選擇和參數維度約減；接著，使用相似性匹配方法匹配每一臺參考發動機和待預測發動機的退化軌跡，獲得待預測發動機的剩余壽命估計值和兩者之間的相似性距離；最后，通過相似性距離生成加權權重，加權這些壽命估計以獲得待預測發動機的剩余壽命。經驗證，本發明對不同的發動機測試樣本均具有較高的壽命預測精度。
【專利說明】
-種基于工況識別和相似性匹配的變工況下航空發動機整機 剩余壽命預測方法
技術領域
[0001] 本發明設及航空發動機剩余壽命預測的技術領域，具體設及一種基于工況識別和 相似性匹配的變工況下航空發動機整機剩余壽命預測方法。
【背景技術】
[0002] 目前我國航空業正處于蓬勃發展階段，但依然存在維護和維修成本高等問題。作 為飛機核屯、動力源的航空發動機，是一種安全和可靠性要求極高，設計生產和維修保養花 費極大的大型復雜設備。資料顯示，自1963年至1975年期間，美國空軍的戰斗機共發生了飛 行事故3824起，其中由于發動機原因導致的事故有1664起，占到總數的43.5%;自1994年至 2000年世界航空運輸機事故共279,其中因發動機的故障導致飛行事故約20%，因此作為飛 機屯、臟的發動機的維修是航空安全保障的重點。世界范圍內航空公司每年大致花費31億美 元用于飛機維修，其中用于航空發動機的維修費用約占總費用的30%，因此需要不斷發展 與改進現有的航空發動機健康管理技術和維修策略，W減少飛行事故的發生并削減維修花 費。隨著傳感器及監測技術的迅速發展，在航空發動機的使用過程中，積累了大量與壽命信 息相關的性能監測數據，采用運類性能退化數據對高可靠和長壽命設備進行壽命預測成為 發展趨勢。
[0003] 當前，國內對航空發動機壽命預測的研究還停留在單個零部件或滑油系統上，針 對航空發動機整機剩余壽命預測的研究基本沒有，并且現有研究都是基于單工況假設下進 行的。然而，航空發動機常工作在復雜多變的工況環境下，變化的工況掩蓋了系統真實的性 能退化趨勢，在此情況下現存的技術方法不能進行準確的剩余壽命預計。
【發明內容】
[0004] 本發明的目的在于:針對變工況下航空發動機整機剩余壽命難于預測的問題，提 出一種基于工況識別和相似性匹配的變工況下航空發動機整機剩余壽命預測方法。
[0005] 本發明采用的技術方案為:一種基于工況識別和相似性匹配的變工況下航空發動 機整機剩余壽命預測方法，其特征在于：
[0006] (1)基于K均值聚類的航空發動機工況識別方法，識別系統在不同時刻所處的工況 后將不同工況下的數據標準化到統一的參數空間；
[0007] (2)為最大限度的挖掘和借鑒參考發動機退化軌跡的衰退規律，通過時間定位和 時間尺度伸縮最優地匹配待預測發動機與參考發動機的性能退化軌跡；
[000引（3)將參考發動機伸縮后軌跡的剩余壽命作為待預測發動機的剩余壽命估計，并 使用相似性距離加權運些估計獲得待預測發動機的剩余壽命。
[0009] 本發明與現有技術相比的優點在于：
[0010] (1)、本發明通過工況識別方法解決了航空發動機壽命預測領域目前尚未很好解 決、而實際中常常存在的變工況問題；
[0011] (2)、本發明提出時間定位和時間尺度伸縮來最大化地匹配參考發動機軌跡與待 預測發動機軌跡，提高了歷史全壽命數據在預測中的參考價值，算法具有較高的預測精度；
[0012] (3)、本發明借鑒集成學習中的某些思想，對待預測發動機壽命進行多次預測后加 權集成獲得待預測的壽命，算法穩健性好，抗干擾能力強。
【附圖說明】
[0013] 圖1為變工況下航空發動機整機剩余壽命預測流程圖；
[0014] 圖2為K均值聚類流程圖；
[0015] 圖3為工況識別流程圖；
[0016] 圖4為變工況工作示意圖；
[0017] 圖5為多工況下的退化特性融合示意圖，其中，圖5(a)為坐標平移示意圖，圖5(b) 為尺度伸縮示意圖；
[0018] 圖6為多工況數據標準化流程圖；
[0019] 圖7為時間定位示意圖；
[0020] 圖8為伸縮匹配原理圖，其中，圖8(a)為時間尺度伸縮原理圖，圖8(b)為最優匹配 示意圖；
[0021 ]圖9為滿輪風扇發動機氣路部件簡圖；
[0022] 圖10為工況聚類圖；
[0023] 圖11為傳感器1，2和3標準化前后的退化數據對比（訓練集1#發動機），其中，圖11 (a)為標準化前的退化數據，圖11(b)為標準化之后的退化數據。
[0024] 圖12為維度約減后的性能退化軌跡(訓練集1#發動機）；
[0025] 圖13為測試集1#發動機的3個壽命估計值，其中，圖13(a)為使用訓練集1#發動機 的估計結果，圖13(b)為使用訓練集2#發動機的估計結果，圖13(c)為使用訓練集3#發動機 的估計結果。
【具體實施方式】
[0026] 下面結合附圖W及【具體實施方式】進一步說明本發明。
[0027] 如圖1所示，基于工況識別和相似性匹配的變工況下航空發動機整機剩余壽命預 測方法，具體步驟如下：
[002引 1工況識別
[0029] 復雜的工況變化常常掩蓋航空發動機的性能退化，為還原真實的退化趨勢，采用 先辨識各狀態的工況，然后進行數據標準化的策略。本發明采用K均值聚類算法對工況參數 進行聚類和分類，從而辨識出系統所處各個狀態的運行工況信息。在數據挖掘中，K均值聚 類算法應用廣泛，是經典聚類算法之一。
[0030] 設X=陸，X2, . . .，Xn}為已知數據集，X中的Xl，X2, . . .，Xn是n個數據對象并且每個 數據對象都是N維的，即Xi= (Xii，Xi2, ...，XinKK均值聚類算法就是要找到含有K個聚類中 屯、的集合C= (Cl ,〇2 , . . . ,Ck} = { (C11 , C12 , . . . , C1N) , (C21 , C22 , . . . , C2N) , . . . , (CK1 , CK2 ,..., CKN)}使得目標函數：
[0031]
(I)
[0032] 其中，m被歸為類Cl的數據對象點數，d(Ci，、)表示聚類中屯、與數據對象的歐幾里 德距離，巧巧女化下：
[0033]
(2)
[0034] K均值聚類算法的核屯、思想是把數據集劃分成使目標函數達到最小值的K個類。具 體步驟如圖2所示。
[0035] 針對訓練集數據通過聚類得到各個工況的聚類中屯、和半徑。對于測試數據首先計 算實時工況參數與各個聚類中屯、的距離，距離較近者即為該工況。工況識別的流程如圖3所 /J、- O
[0036] 2數據標準化
[0037] 復雜系統運行在變工況下時，由于測量尺度的不一致性，狀態特征監測數據并不 能真實反映出系統的退化趨勢。如圖4所示，假設隨著系統的退化，系統的某個特征參數是 呈下降趨勢的，3條實線(紅、黃、藍)展示了系統在巧中運行工況下不同退化的曲線。情況1曲 線(一種包含多個工況的運行模式)是上下波動的，完全不能反映出系統隨著時間的真實退 化趨勢;情況2曲線在一定程度上反映了系統的退化趨勢，但是是虛假的、不真實的；情況3 曲線表象地完全呈現了與真實退化完全相反的趨勢。
[0038] 為實現壽命預測的目的，需要采用數據處理方法，將工況造成的系統監測特征數 據的不同尺度特征去掉，使不同工況特征之間具有可比性、呈現一致的退化特征。直觀的， 就是將圖4中所示的3種不同工況下的退化曲線合并（融合映射巧Ij 一起，如圖5所示。至此， 不同運行工況下的特征數據被轉化到了同一個尺度之下，具有切實的可比性，并能反映出 系統的退化趨勢。從圖形上看，將=種工況合并的兩個核屯、步驟是：首先，將=條曲線的重 屯、移到一起(與X軸重合）；其次，對坐標平移后的數據進行尺度縮放，歸一化到一個可比的 特征空間。
[0039] 多工況數據標準化方法是使用特征層預測方法預測剩余壽命時的預處理環節。使 用該方法，不同工況模式下的數據的N個特征將被分別標準化到一個標定的范圍，然后形成 一個與原始數據特征維數相同的新的腺隹時間序列。圖6展示了多工況數據標準化的流程。
[0040] 本發明采用正態標準化(Z-Score)方法將多工況特征數據標準化。令在時刻ti的 工況參數為111，性能狀態參數化，合在一起構成系統的狀態參數集(1：1，111，義〇4 = 1，2，..，]1。 每個系統狀態都將被聚類到P個工況模式Op之一中，因而形成性能狀態樣本子集lx} W = Ixi UiG〇p}，i = l，2，..，化，p = l，2，. .，P。每個工況模式的樣本均值和樣本標準差為：
[0041]
(3)
[0042] 對屬于工況P的化個狀態參數向量進行標準化：
[0043]
(4)
[0044] 式中，np為屬于工況P的數據狀態數，N為性能參數維度。
[0045] 該標準化方法假定在任一使用循環中，每一個工況發生的概率相等，W確保在同 一運行工況下計算得到的樣本均值和標準差是給定工況模式下特征均值和標準差的無偏 估計。最后，合并不同工況模式下的標準化狀態參數得到系統退化N維狀態時間序列Y = {yi,Y2,. ..,yn}〇
[0046] 3狀態參數維度約減
[0047] 高維狀態監測參數通常很難處理或處理時計算量較大，實際中通常需要對其進行 數據維度的約減。主成分分析(PCAiPrinciple Component Analysis)是應用最廣泛的一種 降維方法之一，是一種統計學方法，在信號處理、模式識別、數字圖像處理等領域已經得到 了廣泛的應用。主成分分析方法基本思想是提取出空間原始數據中的主要特征（主元），減 少數據冗余，使得數據在一個低維的特征空間被處理，同時保持原始數據的絕大部分的信 息，從而解決數據空間維數過高的瓶頸問題。
[004引設航空發動機退化狀態矩陣為X，其列向量Xk=(Xlk，X2k，. . .，xnk)%n維狀態參數， 發動機的某一性能狀態可由Xk描述，Xk的協方差矩陣為：
[0049]
(5)
[0化0]式中，N為退化狀態的采樣數，玄為各狀態參數的均值：
[0化1]
(6)
[0052]求解Rx的全部特征值、Q = I，2, ...，n)和特征向量Vi,將特征值、按照從大到小 的順序排列:MH2〉. . .Hn,則對應的特征向量為ViQ = I，2, ..，n)。樣本Xi投影到特征向量 Vi得到該方向對應的主分量為：
[005；3]
(7)
[0054]所有的特征向量張成一個n維正交空間，X投影到該正交空間得到相應的n維主分 量。特征向量所對應的特征值越大，它在重構時的貢獻也越大，特征值越小的特征向量在重 構時的貢獻就越小。設正交空間中前m個主分量為yi，y2,...，ym，其累計方差貢獻率為：
[0 化 5]
(8)
[0056]當前面少數幾個主分量的累計方差貢獻率足夠大，如h(m)〉95%，即95% W上的原 始數據信息保留在前面幾個主分量中，可只取前m(m<n)個主分量來表征原始信息，在保證 信息完整的情況下，達到數據降維的目的。
[0化7] 4相似性匹配
[005引 （1)時間定位
[0059] 如圖7所示，由于不知道待預測發動機軌跡起點的位置，首先需要進行待預測軌跡 的時間定位，找到參考軌跡中與待預測軌跡初始狀態參數最接近的時間點。然后將兩個軌 跡從該時間點截斷，取截斷點后的兩個軌跡進行后續的進一步匹配計算。設L臺參考發動機 的全壽命循環數分別為{1C}，i = l，2, ...，L，待預測發動機的采樣點數為no,將待預測發動 軌跡與各參考發動機進行時間定位得到一個集合{Ho}，i = l，2,...，L。
[0060] (2)時間尺度伸縮匹配
[0061] 如圖8(a)所示，令離散函數x(t)，t = 0，l, . . .，m和y(t)，t = 0，l, . . .，n分別表示參 考發動機和待預測發動機的退化軌跡。兩個軌跡需要滿足W下條件：（A)相同的起點；（B)自 變量為正整數。
[0062] 首先，用一個因子A去伸展或壓縮x(t)的時間尺度，得到一個臨時軌跡x'（t)=x(A t)，t = 0，lA，. . .，mA;AGr。由于臨時軌跡的采樣點不符合整數采樣要求，使用PCHIP插 值算法插值序列x'（t)從而獲得一個新軌跡x"(t)，t = 0，l,. . .，[mA]。
[0063] PCHIP算法具體為:令a = xo<:?a<. . .<xn = b為區間I = [a,b]的一個劃分，fk，k = 0, 1，. .，n為劃分點的函數值，令hk=xk+i-xk，pk = f'（xk)JCHIP函數Ih(X)滿足如下限制：
[0064] A.定義域為I = [a，b];
[0065] B.在每個子區間Ik=[祉，祉+1]化=0，1，...，11-1)，14^)可^用一個；次樣條多項 式等式Ik(X)表示；
[0066] C.插值條件：Ih(xk)=fk，Ih'（xk)=pk;(k = 0,l,2,...,n)。
[0067] 在每個子區間[Xk，祉+1 ]，PCHIP函數可表示為：
[006引 Ik(X) =fk+Ck,l(X-Xk)+Ck,2(X-Xk)2+Ck,3(X-Xk)3 (9)
[0069] 式中；
[0070]
(10)
[0071] 接著，定義一個相似性距離去量化x"(t)和y(t)之間的相似性：
[0072]
'IP
[0073] 式中，a是一個0-1之間的平滑參數。在本發明中，a = 〇.98。因為近期時間點對未來 退化具有更為重要的影響，因而通過
^距當期時刻較近的時間點分配更大的 權重。[?]表示取整算子。
[0074] 如圖8(b)所示，隨著伸縮因子A的變化，相似性距離隨之變化。顯然運是一個凸優 化問題，也就是說存在一個最優的^DPt是相似性距離化(x(t)，y(t))最小。最優函數值被定 義為最優匹配距離Dom。
[00 巧]
U2)
[0076] 式中，DA(x(t)，y(t))為式（11)定義的相似性距離。
[0077] 5集成加權
[007引L臺參考發動機在時間定位點之后的壽命分別記為{iM}，1 = 1，2, . . .，L，其為：
[0079] Im=Ic-Ho,1 = 1,2, . . . ,L (13)
[0080] 式中，Ic為第1個參考發動機的全壽命，Ho為待預測發動機與第1個參考發動機的 時間定位位置。
[0081] 通過相似性匹配，得到L個最優伸縮因子和最優匹配距離{iADpt，iD〇M}，l = l， 2,...，L，通過最優伸縮因子計算L個待預測發動機的剩余壽命估計值：
[0082] Ir=I入。Pt ? iM-n〇，l = l，2, . . .，L (14)
[008引式中，no為待預測發動機的數據點數，iAopt為最優匹配因子。
[0084] 至此，得至化個相似性距離和對應的L個抓L估計值，用集合{1d〇m，1R}，1 = 1，2,...， L來表示。為獲取待預測發動機的剩余壽命，通過相似性距離生成相似性權重。
r U5)
[0085]
[0086] 式中，Is是第1個參考發動機的相似性的分，其夫
[0087] 加權集成的目的是集成多個壽命估計最終獲得待預測發動機的壽命。最簡單的集 成方法是加權和，待預測發動機RUL的點估計為：
[008引
（16)
[0089] 實驗驗證
[0090] 本發明使用美國航空航天局NASA提供的滿扇發動機數據集的化taset 2(包括測 試集、訓練集）對提出的預測方法進行有效性驗證。該數據集由C - M A P S S (C O m m e r C i a 1 Modular Aero-propulsion System Simulation)大型滿輪風扇發動機仿真模型得到。模型 通過14個參數的輸入來模擬包括風扇、低壓壓氣機、高壓壓氣機、高壓滿輪、低壓滿輪(Fan, LPC，HPC，HPT，LPT)在內的五大旋轉部件失效影響和退化過程，發動機仿真模型主要部件簡 圖如圖9所示。
[0091] 數據由發動機不同運行循環時間點的3個工況參數(飛行高度、飛行馬赫數和油口 桿角度)21個傳感器監測的性能參數組成。傳感器監測值是針對研究模塊，通過發動機熱動 力仿真模型獲取，且包含噪聲。訓練數據包括多組獨立同一單元的多元時間序列樣本，反映 了每個單元從初始到失效的全壽命周期的變化過程。每個部件的初始時刻的退化情況是隨 機且未知的，是否發生退化W及退化程度不同，隨著發動機的運行，當性能降級到一個闊值 的后整個系統發生失效。用于衡量發動機性能狀態的21個監測參數值如表1所示。
[0092] 表1模型輸入參數 「nnno1
[0097] I)工況識別
[0098] 使用訓練集數據通過K均值聚類得到6個工況的聚類中屯、和半徑分別如表3所示， 工況聚類圖如圖10所示。對于測試數據首先計算實時工況參數與6個聚類中屯、的距離，距離 最近者即判為該工況。
[0099] 表3聚類中屯、和半徑
[0100]
[0101] 2)變工況數據標準化
[0102] 即使是同一種狀態監測參數，不同工況下的數據不具有可比性，使用本發明提出 的方法對各工況數據進行標準化處理，去除工況的干擾信息，訓練集1#發動機的傳感器1，2 和3標準化前后的數據如圖11所示。
[0103] 從圖11可W看出，傳感器2和3的監測參數在標準化后呈現出了一定的趨勢性，而 傳感器1仍然沒有明顯的退化趨勢，說明傳感器1完全不能表征任何系統的退化信息。因此， 接下來要進行參數的的選擇或降維，去除沒有趨勢性的傳感器參數，并對有趨勢的參數進 行維度的約減。
[0104] 3)多元參數維度約減
[0105] 訓練集1#發動機狀態參數維度約減后的軌跡如圖12所示。
[0106] 4)相似性匹配計算
[0107] 用訓練集1#、2#、3#發動機對測試集1#發動機進行剩余壽命估計的估計結果如圖 13所示，其中測試機1#發動機的真實剩余壽命為18個循環。
[0108] 5)預測結果分析
[0109] 為驗證所提出方法更多發動機樣本的適用性，使用本發明所提出的方法對測試樣 本集中隨機選擇的10臺變工況下運行的航空發動機進行算法的進一步驗證，預測結果如表 4所示。
[0110] 表4預測結果 「01111
[0112]從表4可知，本發明提出的航空發動機整機剩余壽命預測方法的平均預測百分比 誤差為17%左右，對不同的發動機樣本均具有較高的預測精度，提出的預測方法可為航空 發動機的整機剩余壽命預測提供一個新的研究思路。
【主權項】
1. 一種基于工況識別和相似性匹配的變工況下航空發動機整機剩余壽命預測方法，其 特征在于： (1) 基于K均值聚類的航空發動機工況識別方法，識別系統在不同時刻所處的工況后將 不同工況下的數據標準化到統一的參數空間； (2) 為最大限度的挖掘和借鑒參考發動機退化軌跡的衰退規律，通過時間定位和時間 尺度伸縮最優地匹配待預測發動機與參考發動機的性能退化軌跡； (3) 將參考發動機伸縮后軌跡的剩余壽命作為待預測發動機的剩余壽命估計，并使用 相似性距離加權這些估計獲得待預測發動機的剩余壽命。
【文檔編號】G06F17/50GK106021826SQ201610543483
【公開日】2016年10月12日
【申請日】2016年7月11日
【發明人】劉紅梅, 李連峰, 呂琛, 馬劍
【申請人】北京航空航天大學