基于改進ARMA模型的火箭發動機穩態工況過程實時故障診斷方法研究

(北京航天動力研究所,北京 100076)

0 引言

作為基于數學模型的故障診斷方法中的一種，時間序列分析方法的優點是具有良好的魯棒性、實用性和實時性。而在時間序列分析方法中，應用于液體火箭發動機故障檢測與診斷中應用的主要是AR模型和ARMA模型。

國內外就時間序列分析方法也做了大量的研究。1990年，文獻[1]以SSME為研究對象，對其運行的全過程都應用了ARMA模型進行故障診斷，并與紅線關機算法進行了比較，體現了其快速性及實用性；1997年，吳建軍等[2]提出了基于ARMA模型計算殘差自相關函數置信區間檢驗的故障檢測策略，并介紹了相關建模方法；2002年，張純良等[3]通過改進ARMA模型，實現了對某空間推進系統的故障特征參數提取；2019年，薛薇等[4]針對重復使用的液體火箭發動機，設計了基于ARMA模型的實時故障仿真系統，對典型故障進行了仿真測試，并通過硬件在回路仿真測試，驗證了算法的適用性。

綜上所述，國內外有關于ARMA的研究有很多，但是絕大部分都是關于ARMA模型的訓練，包括模型識別和參數估計等方面，而關于故障閾值的計算，故障的判斷準則的研究卻幾乎沒有。本文以大推力氫氧補燃循環發動機為研究對象，采用了基于改進ARMA模型進行故障診斷，利用改進的閾值計算方法以及連續準則對它主級工況下故障進行了實時診斷，診斷結果良好。

1 基于ARMA模型的故障診斷

1.1 ARMA模型概述

ARMA(p,q)[5]自回歸滑動平均模型是研究時間序列的一種最為重要的分析方法，該方法認為在一組平穩的時間序列中，某一時刻的觀測值不僅和前p步的觀測值有關，還和前q步的擾動有關，且均為線性關系，它的表達式如式(1)所示：

(1)

其中:φi為自回歸部分系數，θi為滑動平均部分系數；p和q分別為自回歸階數和滑動平均階數；xt為t時刻的觀測值，at為t時刻的偏差。特殊的，當p=0時，模型即為移動平均模型，記為MA(q)，當q=0時，模型即為自回歸模型，記為AR(p)。

ARMA模型所研究的是有序時間序列數據的一個內在聯系，不需要考慮影響數據的各外在因素，通過分析有序數據本身的一個規律，來進行短期的預測。當給定一組系統的時間序列時，通過訓練建立ARMA模型，可以得到短期的預測數據，將預測數據與實際數據相互比較，可以判斷系統是否出現故障。ARMA模型使用有一個重要的前提，即所分析的數據必須是平穩數據。

1.2 ARMA模型的訓練

利用時間序列訓練ARMA的主要步驟如圖1所示，首先是讀入數據，利用傳感器測得所需要的參數數據，然后進行數據的處理，接著進行模型識別和參數估計，最后進行模型的檢驗和預測。

圖1 ARMA模型訓練步驟

首先，判斷原始數據是否平穩，如不平穩，則進行差分處理，使其平穩，這是因為ARMA模型只能訓練平穩的數據；接著進行模型識別，判斷適用于ARMA模型、AR模型還是MA模型；然后通過相關定階及參數估計方法確定模型階數和參數，確定模型具體表達形式；最后進行預測和檢驗，如果不合適，則重新進行模型的訓練過程。

訓練ARMA模型時，最為重要的就是模型定階及參數估計的部分，模型階數和參數的準確性直接影響到檢驗結果的好壞[6]。本文確定模型階數的方法為：首先通過分析計算時間序列的自相關函數ACF和偏自相關函數PACF及它們各自的圖像，利用它們的拖尾性和截尾性實現模型階數的初步確定，得到模型階數的幾種組合情況，接著利用AIC信息準則，選擇AIC值最小的模型作為最優的模型階數。而當模型的階數確定之后，便需要進行參數估計，模型參數估計的方式有很多，最常見的方法有：矩估計法、極大似然估計法、最小二乘法等等，本文利用最小二乘法進行參數估計，最后得到準確的預測模型。

ARMA模型所針對的只是一個監測參數，當需要對多個參數進行判斷時，則需要針對每一個參數單獨進行模型訓練，即有多少個監測參數，就需要訓練多少個ARMA模型。

1.3 基于ARMA模型的發動機穩態工況實時診斷算法設計

當液體火箭發動機工作在主級工況下時，各監測參數都會穩定在一定的動態范圍內，我們利用這些參數來訓練ARMA模型，可以預測參數的一個變化規律。當發動機故障發生，我們用模型預測得到的數據和實際數據就會有很大偏差，但這個殘差到達什么程度則判定為故障的研究卻很稀少，國內外學者們就ARMA模型在故障診斷上的應用有很多的研究，但是大多都是針對如何進行訓練模型而展開，本文基于改進ARMA模型的判斷準則出發，使用了新的閾值計算方法以及連續判定準則，實現了故障的零誤報率及快速報警。

一般而言，基于信號和模型的故障診斷算法，判斷故障是否發生主要依靠于閾值準則，即當預測值和實際值之間的偏差量大于某一閾值時，則判斷故障發生。在本文中，為了避免因為偶然因素導致某一瞬間參數變化過大，從而導致誤判的情況，采用了連續準則，即當預測值和實際值之間的偏差量連續幾次超過某一閾值，則判斷故障發生，反之，則未發生，并通過多次試驗，發現當連續次數設定為二次、三次和四次時，系統均會出現誤報警，而當連續次數設定為五次時，診斷情況最為精確，無誤報警出現，所以設定連續次數為5次。

參數閾值的設定對算法的靈敏度和可靠性的影響很大，如果閾值設定過低，則易出現誤報警，如果閾值范圍太大，又會降低算法的靈敏度。我們的改進閾值確定方法參考自適應相關安全限故障檢測算法閾值的設定準則[7]，在主級工況下，參數的閾值計算公式如下:

Pt=Pave+S*N

(2)

其中:Pt為參數的閾值，Pave為參數的平均值，S為參數的標準偏差，N為帶寬系數。Pave和S由主級工況下各參數正常工作的統計數據求得，而帶寬系數N根據正常參數測量值位于安全帶之外的概率來確定。具體確定閾值的方法如下：

1)首先計算前n個樣本數據的平均值Pave和標準偏差值S；

2)接著先令帶寬系數N值為1，初步確定參數閾值，并記錄前n個樣本數據中參數值超過閾值的參數值P及次數M；

3)利用式(3)計算越界參數的規范化偏差值的平均值N′:

(3)

4)確定閾值公式如式(4)所示：

Pt=Pave+S*(1+N′)

(4)

改進的閾值求解公式計算得到的閾值與傳統的閾值求解算法相比，通過實驗驗證，發現誤報警率有明顯的降低，故障診斷效率有明顯的提高。

2 基于改進ARMA模型判別準則的的發動機典型故障分析

2.1 故障模型和故障數據

圖2 氫氧補燃循環發動機仿真系統圖

本文的故障模型來源于我國新一代大推力氫氧補燃循環液體火箭發動機，它以低溫液氫液氧作為推進劑，采用了單富氫預燃室，燃氣并聯驅動氫氧渦輪泵的補燃循環方案。本文利用了Matlab/Simulink工具構建了它的故障仿真模型，具體結構如圖2所示，它可以仿真發動機各部位多種典型故障，并支持發動機模塊化的建模功能，詳細的設計及仿真分析見文獻[8]。

發動機在主級工況時，各監控參數都穩定在一個合理的動態范圍內，發動機正常工作時，這些波動是可以預測的，變化的這些參數可以通過相應的傳感器測量得到，當發動機有故障發生時，這些參數就會出現較為明顯的變化，傳感器測量的參數也會超出穩態工作范圍。本文主要利用改進ARMA模型來實現發動機主級工況的實時故障診斷研究。

一般而言，給發動機仿真模型添加故障有兩種方式：一種是修改發動機系統仿真模型，直接構建新的故障模型；另一種是可通過添加相應故障故障因子的方法，在發動機模型上直接乘以一個故障因子數，故障因子的設置部位和大小反應了故障發生的部位和故障嚴重程度。考慮到第二種方法的易操作性和可靠性，我們選擇添加故障因子的方法來構建發動機故障模型。本文添加的故障類型為主氧泵汽蝕故障，當發生這類故障時，首先氧泵汽蝕會直接導致氧泵效率的明顯的降低，相應的氧泵的關鍵參數諸如轉速、出口的壓力和流量也會降低；接著根據發動機的工作原理以及氧工質流動情況分析，主氧泵連接了氧預壓泵和預燃室，氧流量互通，壓力相互平衡，從流量關系來看，導致了氧預壓泵流量、推力室氧噴前流量和推力室喉部流量的明顯降低，而從壓力關系來看，主氧泵出口壓力的降低，引起了主氧泵泵前壓力和推力室氧噴前壓力的下降。而其它發動機關鍵參數變化，由于發動機的耦合關系，也會出現相應的變化。

通過上述的故障模式分析，得到了關鍵參數的變化情況，對于我們下一節發動機監測參數的選擇有著重要的參考價值。

2.2 訓練ARMA模型

每一個監測參數都需要訓練對應的ARMA模型，所以監測參數的選擇極為關鍵。根據上一節的故障模式分析結果，并結合了監測參數的可監控性、監測參數對故障的敏感性、監測參數對于監控算法的強魯棒性等方面[9]選擇了氧主泵泵后流量、氧主渦輪泵轉速和氧主泵泵后壓力這三個監測參數，表1列出了所選擇的監測參數以及單位。

表1 用于ARMA模型訓練的氧泵汽蝕監測參數

本文設置仿真總時間為7 s，而故障發生時間則為4 s，設定的采樣步長為0.01 s，并增加了高斯白噪聲。首先通過仿真得到了氧泵汽蝕故障三個監測參數的原始數據圖及零均值之后的圖如圖3所示。

圖3 參數原始數據圖和零均值圖

根據直接觀測，正常工作時候的數據無明顯的周期變化，并且沒有明顯的趨勢性，說明數據具有平穩性，則零均值化之后的數據也是零均值平穩時間序列，滿足了ARMA算法建模的要求，可以用于ARMA算法建模，并進行數據的預測和分析。

接著我們利用這三個監測參數正常工作時的數據來進行ARMA模型的訓練，而用于ARMA模型訓練的數據長度也有一定的要求，訓練長度過長時，會使得計算量增大且易導致模型訓練的困難，而當訓練長度過于短時，則會使得所選數據不能完全表示數據的特征，導致訓練的模型準確度不夠。因此，一般模型數據長度的選擇需要經歷試驗來驗證。這里我們分別使用了正常工作時的前200個、前300個以及前400個數據來進行模型的訓練，發現前200個數據得到的模型差距太大，而前300個和前400個數據得到的訓練模型幾乎一樣，所以在考慮到準確性和為了減少計算量的情況下，選擇利用前300個數據來進行模型訓練。首先進行模型識別步驟，分別計算并畫出前300個數據的自相關函數值和偏自相關函數值，如圖4所示。

圖4 自相關函數圖像和偏自相關函數圖像

由于樣本數據的ACF大多都落在了置信區間內，滿足隨機性要求，并且隨著k增大，函數值逐漸趨近于0，則認為該序列具有平穩性，可以按照ARMA模型進行預測分析，下面分別針對三個監測參數建立各自的ARMA模型。

針對氧主泵泵后流量，分析圖4(a)的圖像：根據偏自相關函數圖的拖尾性，我們可以確定自回歸階數p值為1，而根據自相關函數圖像的拖尾性，我們可以確定滑動平均階數q值為4、5和6。因此我們可以初步得到ARMA(1,4)、ARMA(1,5)和ARMA(1,6)這3個模型。接著利用AIC信息準則，分別計算這3個模型所對應的參數以及AIC值，如表2所示。

根據AIC準則，選擇AIC最小的對應的模型作為最優模型，最終確定預測模型為ARMA(1,4)，因此可以得到氧主泵泵后流量的預測模型如式(5)所示：

xt=-0.77759*xt-1+at+0.35776*at-1+0.206652*at-2+ 0.0569704*at-3+0.057282*at-4

(5)

表2 模型參數的擬合結果

同理，針對氧主渦輪泵轉速，得到預測模型如式(6)所示:

xt=0.855249*xt-1+at+0.324822*at-1-

0.00513586*at-2+0.194542*at-3

(6)

同理，針對氧主泵泵后壓力，得到預測模型如式(7)所示:

xt=-0.780675*xt-1+at+0.201901*at-1+

0.00949483*at-2+ 0.125835*at-3+0.214167*at-4

(7)

2.3 改進ARMA模型判別準則的故障診斷與分析

2.3.1 模型精度分析

三個監測參數都有了預測數學模型，接著根據各自的預測模型，我們可以做出各自的預測曲線，并將它與實際曲線相比較，如圖5所示，左側為預測曲線和實際曲線對比圖形，右側為預測值和實際值得差值圖。

圖5 預測結果圖

根據圖5左側圖形觀察可以知道，預測圖形和實際圖形基本重合，預測效果良好，并且在故障發生后，3個監測參數都有了明顯的下降，這和我們之前預測的故障參數變化趨勢一致，這也證明了搭建的發動機故障仿真模型是正確的；右側圖形可以發現，實際值和預測值得差值相比于實際值較小，證明了預測結果精度較高。而對于時間序列模型的預測效果來說，預測的精度是極為重要的，接下來我們通過相關數值來讓結果更加清晰明了。目前，對預測的精度評定主要是基于誤差理論[10]，即用平均絕對誤差(Mean Absolute Error, MAE)和均方誤差(Mean Squared Error, MSE)來衡量。本文選擇平均絕對誤差MAE作為模型預測精度的判斷標準，通過Matlab編程得到了三個監測參數MAE值如表3所示。

表3 監測參數MAE值

通過表3分析可知，三個監測參數的平均絕對誤差和發生故障后的偏移量相比較小，數值上連5%都不到，所以故障發生之前，監測參數的正常波動并不會影響訓練的ARMA模型對故障發生的判斷，可以用于判斷故障是否發生。

2.3.2 基于改進閾值判別準則的故障診斷

我們利用前300個正常工況下的預測值和實際值的偏差值來進行發動機故障的診斷分析，首先基于改進的參數閾值求解式(2)～(4)，確定了各監測參數預測值和實際值差值的閾值，然后設定連續超過閾值次數為5次時報警，并選取這五次中第一次的時間為發生故障的時間，最后通過仿真分析，得到了氧泵汽蝕故障下各監測參數的故障診斷情況如表4所示，上述過程均利用了Matlab編程[11]實現。

表4 各參數的故障診斷情況

由表4結果分析可知：3個監測參數均未出現提前的誤報警，診斷時間都在幾十毫秒的量級，每一個參數都能迅速地診斷出故障，設定的閾值合理，故障診斷較為迅速。這些均證明了所訓練的ARMA模型以及所使用的改進閾值判別準則和連續準則合理且高效，可以用于發動機故障診斷，這對于以后發動機健康監控系統的設計有著積極的參考意義。

3 結論

本文重點針對新一代大推力氫氧補燃循環液體火箭發動機主級工況進行了實時的診斷算法設計，并通過仿真分析驗證了算法的有效性。首先，通過MATLAB/Simulink工具建立了大推力氫氧補燃循環發動機的故障仿真模型，并通過添加故障因子的辦法得到了對應故障模式的原始故障數據；其次，設計了改進的ARMA模型，并改進了閾值求解算法和故障判別準則；最后，通過仿真分析驗證了所訓練的算法的合理性以及改進閾值判別準則的高效性及適用性，滿足發動機故障診斷系統設計的要求。