火箭發(fā)動機渦輪泵故障檢測及預(yù)測算法研究

王姝淇，劉梓琰，王 冠，陳海寶

(1.上海交通大學(xué)，上海 200240；2.北京宇航系統(tǒng)工程研究所，北京 100076)

0 引言

在液體火箭發(fā)動機中，渦輪泵是用于輸送推進劑的核心部件，被稱為液體火箭發(fā)動機的心臟，火箭發(fā)動機的多數(shù)故障都與它有關(guān)，振動則是渦輪泵故障的重要起因之一。渦輪泵工作轉(zhuǎn)速高、工作環(huán)境復(fù)雜惡劣，易發(fā)生組件斷裂、燒蝕等問題，造成嚴重后果。因此有必要研究高效的火箭發(fā)動機故障檢測、故障識別方法，保障火箭的安全性與可靠性。

由于火箭發(fā)動機系統(tǒng)的復(fù)雜性，構(gòu)建基于物理模型的方法十分困難，而近年來，隨著深度學(xué)習(xí)及計算能力的發(fā)展，數(shù)據(jù)驅(qū)動的方法取得了廣泛的應(yīng)用。Zhang等[1]設(shè)計了一種基于端到端的深度學(xué)習(xí)模型，利用振動信號實現(xiàn)軸承故障診斷分析；馬俊添等[2]通過深度殘差網(wǎng)絡(luò)結(jié)合時空卷積網(wǎng)絡(luò)構(gòu)建健康因子完成軸承剩余壽命預(yù)測；Hong等[3]使用小波包分解提取衰減過程，并使用經(jīng)驗?zāi)B(tài)分析算法提取波形趨勢，共同處理振動信號的時頻域特征，最后用自組織神經(jīng)網(wǎng)絡(luò)進行健康預(yù)測；Guo等[4]利用卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Network，CNN)構(gòu)建了健康因子；Ren 等[5]結(jié)合了時域特征和頻域特征，并提出頻譜分割求和方法提取了退化趨勢。

為了進一步揭示機械的退化過程，通?；诒O(jiān)測信號來構(gòu)建健康指標(Health Index，HI)，用以表征機械的健康狀態(tài)。HI的構(gòu)建有不同的方式，如Huang等[6]、Hu等[7]通過RMS進行了滾動元件軸承的剩余使用壽命(Remaining Useful Life，RUL)預(yù)測；Li等[8]通過計算信號的數(shù)學(xué)形態(tài)模式譜構(gòu)建HI；Hu等[9]通過計算旋轉(zhuǎn)諧波振幅的平均對數(shù)比來監(jiān)測齒輪的損傷過程。另外，通過融合多個物理特征構(gòu)建HI模型的方法取得了廣泛的應(yīng)用，如Guo等[10]基于循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Net-work，RNN)構(gòu)建了用于軸承的RNN-HI，并進行了后續(xù)的RUL預(yù)測；Malhotra等[11]通過一個包含編碼器-解碼器的長段時間記憶網(wǎng)絡(luò)(Long Short Term Memory，LSTM)構(gòu)建HI。

基于上述綜述方法，本文首先構(gòu)建時頻混合特征集及精簡特征集，提出了基于特征集的故障預(yù)測方法和HI模型，以及基于時頻譜圖異常譜線的故障檢測方法。

1 數(shù)據(jù)描述

本文使用的數(shù)據(jù)為某型液體火箭發(fā)動機試車數(shù)據(jù)。全部試車數(shù)據(jù)覆蓋共計9臺發(fā)動機29次試車試驗，每次試車提供一個包含大量測量參數(shù)的數(shù)據(jù)集，參數(shù)可分類為緩變數(shù)據(jù)(如氧渦輪出口壓力、燃氣發(fā)生器室壓等)、速變數(shù)據(jù)(如氧渦輪軸向振動、氧泵軸向振動等)。各次試車的基本狀態(tài)統(tǒng)計如表1所示，其中試驗編號方式為發(fā)動機編號-該發(fā)動機試車次數(shù)，如09-1代表9號發(fā)動機第一次試車(后文對發(fā)動機統(tǒng)一表述為xx號發(fā)動機，對試車試驗統(tǒng)一表述為xx-yy次試車)，時長代表該次試車試驗點火時間長度，健康狀態(tài)與故障模式代表該次試車后發(fā)動機的狀態(tài)，由試車試驗后拆機檢查得到，健康狀態(tài)為無法判斷的則代表未拆機檢查，可能正常也可能故障。

由于振動信號具有響應(yīng)快、敏感、特征豐富的特點，本文主要對其中的振動數(shù)據(jù)進行分析。數(shù)據(jù)中振動參數(shù)共有9個：氧泵殼體軸向振動、氧泵殼體徑向振動、氧泵殼體切向振動、氧渦輪軸向振動、氧渦輪徑向振動、氧渦輪切向振動、泵后支座軸向振動、泵后支座徑向振動、泵后支座切向振動。

實際數(shù)據(jù)記錄中并非所有參數(shù)均有記錄，考慮到希望覆蓋更多次試車，本文中的特征集構(gòu)建采用各次試車的共有振動參數(shù)：氧泵殼體軸向振動、氧泵殼體徑向振動、氧泵殼體切向振動，而故障預(yù)測方法和HI模型方法均基于特征集完成；而時頻譜圖的異常檢測采用具有異常譜線特點的氧渦輪軸向振動、氧渦輪徑向振動、氧渦輪切向振動數(shù)據(jù)。

2 特征集構(gòu)建及精簡

航天發(fā)動機是一個高度復(fù)雜的精密系統(tǒng)，而實際測量的振動數(shù)據(jù)包含大量無用信息，混雜著固有振動頻率、零部件的振動、故障缺陷引起的振動、環(huán)境噪聲，因此需提取混合信號中的有效特征。下面首先構(gòu)建時域頻域統(tǒng)計特征集[12]，并通過特征矢量降維對時頻域特征集進行精簡。

2.1 時頻混合特征集

從時域、頻域兩方面分別進行統(tǒng)計特征集的計算，混合后形成完整的時頻混合特征集。

2.1.1 時域特征集

時域分析根據(jù)傳感器采集到的振動信號在時間域上的統(tǒng)計特征作為輸入特征。通常通過對時域信號切片，計算切片窗口內(nèi)的統(tǒng)計特征來代表該切片的特征，以此進行數(shù)據(jù)降維。目前，常用的時域特征及相應(yīng)計算公式如表2所示。

表2 時域特征集構(gòu)成Tab.2 Time-domain feature set

另外，不同監(jiān)測數(shù)據(jù)表示不同的物理特性，具有不同的量綱和數(shù)量級，為消除數(shù)據(jù)不規(guī)范對效果的影響，進一步對特征集進行歸一化處理。

2.1.2 頻域特征集

振動信號包含各類要素信息，不同要素通常對應(yīng)于不同頻段的信號，可通過分析頻域特征抽取所需的故障要素。信號經(jīng)過傅里葉變換后，得到的各頻率的幅值大小，對該頻域信息做統(tǒng)計特征的分析作為頻域特征提取，常用的頻域特征如表3所示。與時域特征集構(gòu)建類似，完成頻域特征集構(gòu)建后，同樣須要對其進行歸一化處理。

表3 頻域特征集構(gòu)成Tab.3 Frequency-domain feature set

分別按照表2與表3所示的特征，構(gòu)建氧渦輪泵軸向、徑向與切向振動3個參數(shù)的特征集，每個參數(shù)11個時域特征，12個頻域特征共計23個特征。將上述3個參數(shù)特征混合在一起，形成69維時頻特征集，用于生成第3章算法所需的精簡特征集以及第5章HI模型構(gòu)建算法的輸入。

2.2 精簡特征集

上述69維時頻特征集雖然可以從不同的方面反映發(fā)動機的狀態(tài)，但特征總維度大，包含冗余信息。為了降低數(shù)據(jù)維度，采用主成分分析(Principal Component Analysis，PCA)方法來進行特征降維，從原始數(shù)據(jù)集中精簡得到顯著特征集。

PCA通過正交變換將可能存在相關(guān)性的變量重新組合成一組新的線性不相關(guān)的綜合變量，以實現(xiàn)采用較少的綜合變量來盡可能多地反映原來變量大部分信息的目的。2.1節(jié)中構(gòu)建的69維特征集，經(jīng)PCA計算后成為69維的正交特征，按每個特征的方差大小排序，可得到表4。其中單個成分方差貢獻率為某主成分的方差與全部主成分方差和的比值，累積貢獻率為前k個主成分方差貢獻率的累加。PCA降維的有效性常用降維后主成分的累積貢獻率來表示，如果前k個主成分的方差累計貢獻率超過了85%，那么說明用前k個主成分去代替原來的n個變量后，不能解釋的方差不足15%。

表4 69維特征集PCA主成分方差貢獻率Tab.4 Variance contribution of principal component for 69 features

在氧泵軸徑切的69維時頻混合特征集中，前10個成分的方差貢獻率及累積貢獻率如表4所示。前4個主成分的累積貢獻率達到了95%，而隨著主成分的增加，累積貢獻率上升逐漸緩慢。因此僅保留前4個主成分即可保留足夠的信息量，將PCA后的前4個主成分作為精簡特征集，用于故障預(yù)測方法。

3 基于精簡特征集的故障預(yù)測方法

考慮火箭射前發(fā)動機短程試車的場景需求，因此本章開展故障預(yù)測方法研究，實現(xiàn)通過上一次試車數(shù)據(jù)預(yù)測下一次試車是否會出現(xiàn)異常的功能。該問題屬于時間序列分類(Time Series Classification，TSC)問題。

根據(jù)Fawaz等[13]對傳感器讀數(shù)、光譜圖、模擬數(shù)據(jù)等不同領(lǐng)域的數(shù)據(jù)集的實驗結(jié)果，ResNet在時間序列分類上有出色的表現(xiàn)。ResNet通過調(diào)整模型結(jié)構(gòu)，引入快捷連接的方式構(gòu)建輸出為F(x)+x的殘差塊，由該殘差塊疊加的網(wǎng)絡(luò)更容易優(yōu)化，能使網(wǎng)絡(luò)層數(shù)大大加深。該結(jié)構(gòu)如圖1所示，上文中F(x)對應(yīng)圖中殘差塊輸出中的W2·ReLU(W1·x)，其中W1和W2分別為兩網(wǎng)絡(luò)層的權(quán)重矩陣。

圖1 ResNet網(wǎng)絡(luò)殘差塊結(jié)構(gòu)Fig.1 Structure of residual block

鑒于ResNet在時間序列分類問題的良好效果，這里也采用ResNet網(wǎng)絡(luò)來完成，多維時間序列采用精簡特征集，通過訓(xùn)練網(wǎng)絡(luò)完成該特征序列到下一次試車預(yù)測結(jié)果的映射。網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示。

圖2 殘差結(jié)構(gòu)的時間序列分類網(wǎng)絡(luò)Fig.2 Classification network with residual structure

本章使用的數(shù)據(jù)為氧泵殼體軸向、氧泵殼體徑向、氧泵殼體切向振動數(shù)據(jù)的精簡特征集，將下一次試車的“正?！被颉爱惓！鼻闆r作為當前試車數(shù)據(jù)的訓(xùn)練標簽。在全部29次試車的數(shù)據(jù)中，以異常試車結(jié)束的發(fā)動機由于沒有維修，下一次試車必然也是異常的，因此標簽為“異?！?；而以正常試車結(jié)束的發(fā)動機，包括09-1、14-6、28A-6、33-2次試車，其下一試車狀態(tài)未知，從而標簽“未知”。將上述4個標簽“未知”的試車去除后，剩余可用的數(shù)據(jù)共計25次試車。將25個可用試車隨機劃分為19個訓(xùn)練集試車和6個測試集試車。

為了增加樣本量、并挖掘局部的特征，提升預(yù)測結(jié)果的穩(wěn)定性，采用了時間切片處理和投票機制。時間切片和標簽設(shè)置如圖3所示，切片時長設(shè)置為100 s，切片的標簽與當前試車的預(yù)測標簽一致，即每個時間切片樣本都對應(yīng)著下一次試車的分解檢查情況(正常：0；異常：1)，進而作為網(wǎng)絡(luò)監(jiān)督訓(xùn)練的輸入信息。經(jīng)過時間切片后的數(shù)據(jù)樣本量分別為：訓(xùn)練集77個樣本，測試集30個樣本。

圖3 時間切片和標簽設(shè)置示意圖Fig.3 Diagram of time slice and label setting

對于一個包含多個時間切片的試車，最終預(yù)測結(jié)果采用投票機制：當預(yù)測異常的切片數(shù)目多于正常時，認為綜合預(yù)測結(jié)果為異常；當預(yù)測正常的切片數(shù)目多于異常時，認為綜合預(yù)測結(jié)果為正常。最終預(yù)測結(jié)果以各個切片的綜合結(jié)果來決定，以此提高預(yù)測的準確率和可靠性。

測試集預(yù)測效果如表5所示，對于測試集的6次試車預(yù)測均準確，最終效果準確率100%。值得指出，這里由于隨機劃分的因素，測試集均是500 s的試車長度，因此都是5個時間切片，實際上對于任意長度的試車可同理預(yù)測。

表5 測試集預(yù)測結(jié)果Tab.5 Prediction results on test dataset

4 基于時頻譜圖異常譜線檢測的故障檢測方法

前期研究表明氧渦輪泵試驗中32號發(fā)動機的異常試車在氧渦輪軸向、氧渦輪徑向、氧渦輪切向振動參數(shù)的時頻譜圖5 kHz附近呈現(xiàn)傾斜異常譜線的特點，而正常試車在該頻段并沒有出現(xiàn)傾斜譜線，且該情況在多臺設(shè)備出現(xiàn)，因此推測時頻譜圖異常譜線，是判斷氧渦輪泵異常與否的重要判據(jù)。

圖5 氧渦輪徑向振動時頻譜圖Fig.5 Time-spectrum diagrams of oxygen turbo pump radial vibration

本章針對時頻譜圖的異常譜線檢測展開研究，采用的振動參數(shù)是氧渦輪軸向、氧渦輪徑向、氧渦輪切向，采用32-1、32-2、32-3、32-4、32-5、32-6次試車。通過時頻域分析繪制時頻譜圖，采用神經(jīng)網(wǎng)絡(luò)方法對時頻譜圖中的異常譜線進行檢測，由此確定故障檢測的目標。

4.1 時頻譜圖繪制

上述時域特征、頻域特征是單獨在時域或頻域的分析，但考慮到信號的頻率特性可能隨著時間變化，為了捕獲這一時變特性，時頻域分析提供一種在時域和頻域均有分辨率的分析方法。短時傅里葉變換是時頻域分析的主要方法，其定義如下

(1)

式中，x(m)為輸入信號，ω(m)為窗函數(shù)。通過短時傅里葉變換獲得時間和頻率的二維函數(shù)X(n，ω)，利用二維圖像表示，即得到信號時頻譜圖。

前期研究表明，氧渦輪泵試驗異常試車在時頻譜圖5 kHz附近呈現(xiàn)傾斜異常譜線的特點，因而重點關(guān)注5 kHz附近的頻譜，同時放大頻段4.5～5.5 kHz 繪制時頻譜圖。圖4、圖5、圖6分別是32號發(fā)動機6次試車的氧渦輪軸向、氧渦輪徑向、氧渦輪切向振動參數(shù)在5 kHz附近的時頻譜圖，由左至右，由上至下分別為32-1至32-6次試車。

圖6 氧渦輪切向振動時頻譜圖Fig.6 Time-spectrum diagrams of oxygen turbo pump tangential vibration

圖4～圖6中，32-1至32-4均表現(xiàn)出高亮豎紋特征，而32-5、32-6則表現(xiàn)出高亮斜紋特征(紅圈標記處)，可見32號發(fā)動機的氧渦輪軸向、徑向、切向振動數(shù)據(jù)在狀態(tài)異常時均表現(xiàn)出相近的異常譜線特征，這為后續(xù)異常譜線識別分析檢測奠定基礎(chǔ)。

4.2 基于ResNet的故障檢測方法

ResNet網(wǎng)絡(luò)在圖像特征提取方面具有較好的優(yōu)勢，本文將采用ResNet對時頻譜圖檢測。為了挖掘小時間范圍的異常譜線，將頻段為4.5～5.5 kHz的時頻譜圖做每1 s的切片，如圖7所示。

圖7 局部時頻譜圖獲取方法Fig.7 Method of obtaining local time-spectrum diagram

根據(jù)局部時頻譜圖中有無斜向譜線作為第一個訓(xùn)練標簽，考慮到不同工況下的不同情況，使用工況作為第二個訓(xùn)練標簽。類別標簽，即有無異常，用1表示異常，用0表示正常；工況標簽分別用0，1，2表示低、中、高工況。利用32-5和32-6次試車的氧渦輪軸向振動參數(shù)的時頻譜圖切片作為數(shù)據(jù)集，隨機劃分為訓(xùn)練集和測試集。將數(shù)據(jù)的局部時頻譜圖按秒切片，將每個切片圖像轉(zhuǎn)換成 (3，224，224)維度的數(shù)組，作為ResNet網(wǎng)絡(luò)輸入，以對應(yīng)的類別標簽與工況標簽作為網(wǎng)絡(luò)輸出的期望。圖8、圖9分別展示了ResNet網(wǎng)絡(luò)訓(xùn)練和測試的結(jié)果中，由每個局部時頻譜圖的1 s 時間切片樣本經(jīng)ResNet網(wǎng)絡(luò)預(yù)測得到的類別與工況與其真實狀態(tài)的對比情況，其中橫坐標的值代表樣本的序號，每個樣本代表1 s時間切片的時頻譜圖。

圖8 訓(xùn)練效果Fig.8 Training effect

圖9 測試效果Fig.9 Testing effect

圖9表明了算法的測試結(jié)果，在每個橫坐標對應(yīng)的樣本上，若橙色點與藍色點重合則代表準確，反之則不準確，經(jīng)過統(tǒng)計，在200余個測試樣本中，異常譜線檢測的準確率約77%，工況識別的準確率約為99%，證明了利用小樣本數(shù)據(jù)進行故障診斷的可行性。

5 基于時頻混合特征集的HI模型構(gòu)建

為了為氧轉(zhuǎn)子運行提供一個能夠及時發(fā)現(xiàn)隱患的評價工具，應(yīng)用上述構(gòu)建的69維時頻混合特征集來構(gòu)建HI，該特征集的特征量較多，從而可以對HI有更好的擬合。本節(jié)使用包含19、32、21、28A、33、14、09號發(fā)動機共26次試車的數(shù)據(jù)，將試車數(shù)據(jù)按2.1節(jié)中方法構(gòu)建的氧泵殼體軸向、氧泵殼體徑向、氧泵殼體切向振動數(shù)據(jù)的69維完整時頻混合特征集。

將同一個發(fā)動機不同次試車在去除不穩(wěn)定的啟動段后連接起來，模擬HI問題的連續(xù)運行到失效的數(shù)據(jù)形式，構(gòu)建線性回歸的HI模型。在目標HI退化曲線中最常用的、最成功的形式是指數(shù)衰減[14]，由于指數(shù)衰減模式的前期平緩、后期加速的模式更符合實際情況，如下式所示

ht=

(2)

式中，β代表起始為1的部分占總體的比例(也是末尾為0的部分占總體的比例)，L為數(shù)據(jù)長度。

在對HI進行擬合時，最常用的是線性回歸方法[14-15]，邏輯回歸會因扭曲系統(tǒng)原有的退化模式(如指數(shù)模式)而不適用[11]。為了保留信號/特征中的原始模式，采用線性回歸模型對指數(shù)衰減進行擬合，線性回歸模型為

(3)

式中，f(x)代表HI的值，x=(x1，x2，…，xN)是N維特征矢量，(α，β)=(α，β1，β2，…，βN)是N+1個模型參數(shù)，該模型完成從N維矢量到HI值的映射。根據(jù)表1，僅有19號發(fā)動機的6次連續(xù)試車是確定滿足條件的從正常到失效的試車(19-1正常，19-2正常，19-3正常，19-4正常，19-5無法判斷，19-6異常)，因此基于19號發(fā)動機建立損傷模型，擬合效果如圖10所示。對存在異常的試車驗證，其數(shù)據(jù)如表6所示。

圖10 19號發(fā)動機擬合效果Fig.10 Fitting effect of No.19 engine

表6 異常試車數(shù)據(jù)詳情Tab.6 Data information of abnormal operation

在跨發(fā)動機驗證時，模型輸出結(jié)果體現(xiàn)出擾動較大、毛刺較多的情況，為了更好地判斷曲線趨勢，還對模型輸出結(jié)果曲線進行平滑降噪處理，通過S-G濾波降低噪聲干擾，放大趨勢性的變化，降噪后體現(xiàn)為結(jié)果圖中的橙色曲線。

對于異常的32和21號發(fā)動機驗證如圖11所示。由圖11可以看出，32號發(fā)動機前面健康值一直較高，但在末尾體現(xiàn)出明顯的衰退趨勢，表明32號發(fā)動機在末尾有明顯退化趨勢，對應(yīng)著32-6次試車可能存在異常，這與32-6次試車的實際標簽相符合；21號發(fā)動機在末尾也存在明顯的健康衰退趨勢，因此認為其存在異常，這與21-2次試車的實際標簽符合。另外，參考表1，24與34號發(fā)動機，各自僅有一次時長為100 s的試車，且事后檢查存在故障。對于這類發(fā)動機該方法并不適用，這是由于經(jīng)過數(shù)據(jù)預(yù)處理和特征提取后的數(shù)據(jù)長度非常小，構(gòu)建的HI曲線可能完整地處于異常狀態(tài)，難以觀察到退化趨勢，而各類試車故障模式，對應(yīng)故障狀態(tài)的HI不同，沒有趨勢僅觀測數(shù)值難以進行判斷。

(a)32號發(fā)動機

對于正常發(fā)動機的HI趨勢曲線如圖12所示，分別是28A、33、14、09號發(fā)動機，可見均未表現(xiàn)出健康值的下降趨勢，沒有顯著的退化特征，其中33號發(fā)動機構(gòu)建的HI曲線時間長度短，但計算得到的HI很高，因此也可以表征該發(fā)動機狀態(tài)正常。

綜上所述，在跨發(fā)動機驗證中，異常發(fā)動機能夠在HI曲線中表現(xiàn)出明顯的下降趨勢，表征了健康狀態(tài)的衰退；而正常發(fā)動機未表現(xiàn)出健康值的下降趨勢，沒有顯著的退化特征。證明該方法可以對渦輪泵健康狀態(tài)進行有效評估。

(a)28A號發(fā)動機

6 總結(jié)

本文提出了針對某型液體火箭發(fā)動機渦輪的故障預(yù)測、故障檢測與健康狀態(tài)評估方法，并在真實火箭發(fā)動機數(shù)據(jù)集上進行了測試驗證。在故障預(yù)測方面，基于精簡特征集和ResNet實現(xiàn)了對下一次試車正常或異常的預(yù)測，引入時間切片的方法，最終達到100%的預(yù)測準確率。在故障檢測方面，通過時頻域分析對于時頻譜圖中的重點關(guān)注頻段體現(xiàn)出的異常譜線與異常之間存在的關(guān)聯(lián)問題，設(shè)計ResNet網(wǎng)絡(luò)來實現(xiàn)對于異常譜線的識別，達到了較好的識別準確率。在健康狀態(tài)評估方面，通過設(shè)立的經(jīng)驗退化趨勢函數(shù)和時頻混合特征集擬合試車數(shù)據(jù)實現(xiàn)健康評估，并在發(fā)動機渦輪泵數(shù)據(jù)進行測試，驗證該方法可以量化表征渦輪泵健康狀態(tài)。