基于自適應(yīng)擬合建模的航空發(fā)動(dòng)機(jī)健康參數(shù)估計(jì)*

顧嘉輝， 魯 峰， 黃金泉， 強(qiáng)子健

(南京航空航天大學(xué) 能源與動(dòng)力學(xué)院 江蘇省航空動(dòng)力系統(tǒng)重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210016)

0 引 言

航空發(fā)動(dòng)機(jī)長期在高溫、高壓和高轉(zhuǎn)速的惡劣環(huán)境中工作，會(huì)發(fā)生緩慢的性能退化現(xiàn)象，影響飛行安全[1，2]。通過收集商用航空發(fā)動(dòng)機(jī)的機(jī)載傳感器數(shù)據(jù)可以對(duì)發(fā)動(dòng)機(jī)的狀態(tài)進(jìn)行監(jiān)視與故障診斷[3～5]，分析性能退化的程度幫助制定維修策略，從而提高飛行安全和降低航班延誤率。

航空發(fā)動(dòng)機(jī)的氣路健康狀態(tài)常用各個(gè)部件的流量和效率的相對(duì)變化量(即健康參數(shù))來表示，然而由于其無法直接測量，只能根據(jù)機(jī)載傳感器的測量值，如兩個(gè)轉(zhuǎn)子轉(zhuǎn)速和各個(gè)截面的壓力和溫度通過算法進(jìn)行估計(jì)。常用的方法有各種濾波算法[5～8]。

在狀態(tài)估計(jì)方面，卡爾曼濾波算法研究得較多[6～10]，然而發(fā)動(dòng)機(jī)部件級(jí)模型計(jì)算量較大，使得非線性濾波算法目前難以滿足機(jī)載實(shí)時(shí)性要求，同時(shí)發(fā)動(dòng)機(jī)線性狀態(tài)變量模型(state variable model,SVM)存在建模誤差，使得線性卡爾曼濾波算法對(duì)健康參數(shù)的估計(jì)精度較低。

常用的發(fā)動(dòng)機(jī)狀態(tài)變量模型建模方法有：擾動(dòng)法[11]、擬合法[12]、改進(jìn)擬合法[13]等。但系數(shù)矩陣僅通過對(duì)發(fā)動(dòng)機(jī)模型在某些特定的性能退化下擾動(dòng)或擬合得到，對(duì)于這些特定情況以外的性能退化情況，由于發(fā)動(dòng)機(jī)非線性特性強(qiáng)，線性模型精度較低，因此健康參數(shù)的估計(jì)誤差會(huì)隨著飛行循環(huán)數(shù)的增加，逐漸變大。本文提出自適應(yīng)擬合法來建立含健康參數(shù)的狀態(tài)變量模型。該方法利用上一個(gè)采樣周期內(nèi)健康參數(shù)的估計(jì)值實(shí)時(shí)重新計(jì)算健康參數(shù)相關(guān)系數(shù)矩陣，以提高線性模型在各種退化程度下的精度，從而使得加載了卡爾曼濾波器的狀態(tài)變量模型能夠在全壽命期限內(nèi)較準(zhǔn)確地估計(jì)健康參數(shù)。

1 基于自適應(yīng)擬合法建模的航空發(fā)動(dòng)機(jī)狀態(tài)變量模型

1.1 航空發(fā)動(dòng)機(jī)非線性模型

本文使用的航空發(fā)動(dòng)機(jī)非線性模型為某型大涵道比渦扇發(fā)動(dòng)機(jī)部件級(jí)模型[14]，共有5個(gè)旋轉(zhuǎn)部件：風(fēng)扇(fan)、低壓壓氣機(jī)(low-pressure compressor,LPC)、高壓壓氣機(jī)(high-pressure compressor,HPC)、高壓渦輪(high-pressure turbine,HPT)和低壓渦輪(low-pressure turbine,LPT)。

假設(shè)該發(fā)動(dòng)機(jī)部件級(jí)模型可以表示為

(1)

式中 狀態(tài)量x=[N1N2]T為風(fēng)扇轉(zhuǎn)速與高壓轉(zhuǎn)速,輸入量u=[WfVSV]T為燃燒室燃油流量與高壓壓氣機(jī)導(dǎo)葉角度,飛行條件v=[HMa]T為飛行高度和飛行馬赫數(shù),輸出量y=[N1N2T25P25T3P3T495P13]T為兩轉(zhuǎn)子轉(zhuǎn)速、低壓壓氣機(jī)出口總溫和總壓、高壓壓氣機(jī)出口總溫和總壓、排氣溫度以及外涵道出口總壓。

1.2 航空發(fā)動(dòng)機(jī)狀態(tài)變量模型

在發(fā)動(dòng)機(jī)某穩(wěn)態(tài)點(diǎn)(x0,u0,v0)附近對(duì)式(1)進(jìn)行泰勒級(jí)數(shù)展開，忽略二階及以上分量，就可以得到發(fā)動(dòng)機(jī)在該穩(wěn)態(tài)點(diǎn)附近的小偏差線性狀態(tài)變量模型

(2)

式中Δx為發(fā)動(dòng)機(jī)狀態(tài)量的增量，Δu為輸入量的增量，Δy為發(fā)動(dòng)機(jī)輸出量的增量，A,B,C,D為相應(yīng)維數(shù)的系數(shù)矩陣。考慮航空發(fā)動(dòng)機(jī)在服役一段時(shí)間后，其氣路部件性能會(huì)不可避免地發(fā)生退化，從而導(dǎo)致發(fā)動(dòng)機(jī)工作狀態(tài)偏離其額定狀態(tài)。因此在狀態(tài)變量模型中增加了能夠表示發(fā)動(dòng)機(jī)部件性能退化程度的健康參數(shù)，描述真實(shí)發(fā)動(dòng)機(jī)動(dòng)態(tài)特性的狀態(tài)變量模型應(yīng)為

(3)

參照文獻(xiàn)[15]給出的某型航空發(fā)動(dòng)機(jī)的每個(gè)健康參數(shù)在不同飛行循環(huán)數(shù)下的大小可知，在飛行循環(huán)數(shù)6 000后，低壓渦輪流量和效率的變化依然非常小，絕對(duì)值都在1 %左右。同時(shí)，考慮到卡爾曼濾波算法中狀態(tài)量的可觀性要求，式(3)中，Δh為不包含低壓渦輪流量和效率的8個(gè)健康參數(shù)組成的列向量

Δh=[SW1SE1SW2SE2SW3SE3SW4SE4]T

(4)

式中SW,SE分別為對(duì)應(yīng)流量和效率健康參數(shù)、下標(biāo)1～4分別對(duì)應(yīng)風(fēng)扇、低壓壓氣機(jī)、高壓壓氣機(jī)和高壓渦輪。式(3)中Δu和Δh對(duì)于狀態(tài)變量模型的作用相同，在建模時(shí)都可以作為控制量來考慮，因此式(3)可以改寫成

(5)

1.3 自適應(yīng)擬合建模法

自適應(yīng)擬合法首先參考文獻(xiàn)[5]發(fā)動(dòng)機(jī)參數(shù)性能退化的范圍對(duì)每個(gè)旋轉(zhuǎn)部件健康參數(shù)的可能變化范圍進(jìn)行劃分，如表1。

表1 各部件健康參數(shù)劃分范圍 %

可知，在進(jìn)行健康參數(shù)估計(jì)前：不考慮低壓渦輪流量和效率健康參數(shù)，則可能出現(xiàn)的氣路健康參數(shù)的組合有4×4×5×4=320種。將所有可能的健康參數(shù)組合記為

Δh(i)=[SW1(i)SE1(i)SW2(i)SE2(i)SW3(i)SE3(i)SW4(i)SE4(i)],i=1,2,…,320

(6)

將320組健康參數(shù)組合注入部件級(jí)模型中得到由于性能退化引起的穩(wěn)態(tài)狀態(tài)偏差和輸出偏差，分別記為

(7)

(8)

在開始健康參數(shù)估計(jì)前，將輸出偏差，即8×320數(shù)組和所有系數(shù)矩陣都以單浮點(diǎn)精度(MATLAB中為Single)形式保存在mat.文件下。

機(jī)載進(jìn)行健康參數(shù)估計(jì)時(shí)，根據(jù)前一時(shí)刻健康參數(shù)的估計(jì)值從mat.文件中選擇最接近的、線性無關(guān)的8組健康參數(shù)組合Δh(si),i=1,2,…，8(行向量)以及對(duì)應(yīng)的Δy(si),Δx(si),i=1,2,…，8(列向量)，利用于更新L,M矩陣：對(duì)i=1且不考慮控制量，當(dāng)系統(tǒng)穩(wěn)定后有

(9)

類似的方程組共有8組，對(duì)應(yīng)選中的8組健康參數(shù)組合。將8組方程組進(jìn)行整理，對(duì)式(9)的第一式有

(10)

式中 ΔH1(si)為2×16維矩陣，L(j,:)為L矩陣第j行所有元素組成的行向量，0為1×8維零矩陣。該式中未知量僅為L(1,:),L(2,:)，即L矩陣可解。對(duì)式(9)的第二式有

(11)

式中 ΔH2(si)為6×48維的矩陣，M(j,:)為M矩陣的第j行所有元素組成的行向量，0為1×8維的零矩陣，()3:8為該列向量的第3到第8行。同理，M矩陣可解。將該時(shí)刻計(jì)算得到的L,M矩陣代替式(8)中原來的L,M矩陣，再進(jìn)行健康參數(shù)估計(jì)，然后下一時(shí)刻重復(fù)該步驟。

關(guān)于線性無關(guān)樣本的選擇方法，按各部件健康參數(shù)組合的歐氏距離進(jìn)行排序，以風(fēng)扇部件為例說明如下：假設(shè)前一時(shí)刻風(fēng)扇流量和效率健康參數(shù)的估計(jì)值為：S1=-3.5 %,S1=-3.4 %，對(duì)應(yīng)圖1中菱形所示，圖中另外4個(gè)離散圖形對(duì)應(yīng)表1中風(fēng)扇健康參數(shù)的4組樣本。

圖1 風(fēng)扇健康參數(shù)樣本分布

計(jì)算4組樣本與當(dāng)前估計(jì)值對(duì)應(yīng)的風(fēng)扇健康參數(shù)的歐氏距離的平方r1(i)(下標(biāo)1～4分別對(duì)應(yīng)風(fēng)扇、低壓壓氣機(jī)、高壓壓氣機(jī)和高壓渦輪，下同)

r1(i)=(S1-SW1(i))2+(S1-SE1(i))2

(12)

將r1(i)(i=1,2,3,4)中最小距離平方對(duì)應(yīng)的樣本序號(hào)(即表1中的序號(hào))記為q1。類似得到q2～q4，則8組線性無關(guān)樣本的各部件健康參數(shù)對(duì)應(yīng)的序號(hào)可以選取為

(13)

相比改進(jìn)擬合法，自適應(yīng)擬合法不再更新式(8)中的A,B,C,D,Kst矩陣，僅由前一時(shí)刻健康參數(shù)的估計(jì)值從mat.文件中選取8組健康參數(shù)組合，按照式(10)、式(11)簡單計(jì)算L,M矩陣，以此提高線性狀態(tài)變量模型在各種性能退化情況下的模型精度，使得結(jié)合了線性卡爾曼濾波器的該線性模型能夠在全壽命期限內(nèi)準(zhǔn)確估計(jì)健康參數(shù)。

2 設(shè)計(jì)點(diǎn)數(shù)字仿真驗(yàn)證

2.1 性能緩慢退化的仿真驗(yàn)證

本文在設(shè)計(jì)飛行點(diǎn)(H=10 668 m,Ma=0.785)從慢車以上(N1=70 %)到最大狀態(tài)(N1=110 %)采用自適應(yīng)參數(shù)擬合法建立13組離散的系數(shù)矩陣Aa,Ba,Ca,Da,Kst與320種輸出偏差Δy(i)組成的高維數(shù)組，對(duì)于該飛行點(diǎn)的其他狀態(tài)以風(fēng)扇換算轉(zhuǎn)速為索引，采用線性插值的方法獲得系數(shù)矩陣與高維數(shù)組。

數(shù)字仿真時(shí)分別采用改進(jìn)擬合法與自適應(yīng)擬合法建模，使用線性卡爾曼濾波器對(duì)健康參數(shù)進(jìn)行估計(jì)，以驗(yàn)證自適應(yīng)擬合法建模的有效性和精度。對(duì)于改進(jìn)擬合法建模時(shí)各健康參數(shù)的階躍量統(tǒng)一選取絕對(duì)值大小為5 %。數(shù)字仿真時(shí)間為60 s，采樣周期為25 ms(即2 400個(gè)采樣周期)，性能退化從第二個(gè)采樣周期開始線性緩慢變化，到第60 s時(shí)退化到對(duì)應(yīng)6 000飛行循環(huán)數(shù)[5]。此外，雖然狀態(tài)變量模型中沒有包含低壓渦輪的流量和效率健康參數(shù)，但數(shù)字仿真中依然參考文獻(xiàn)[5]參數(shù)退化表中飛行循環(huán)數(shù)6 000參數(shù)退化情況，加入了低壓渦輪的性能退化。同時(shí)假設(shè)輸出傳感器含高斯白噪聲，滿足vk～N(0,Rk)，其中Rk=0.0022I8×8，以及加入PID控制器對(duì)風(fēng)扇物理轉(zhuǎn)速閉環(huán)控制[17]。本文的所有數(shù)字仿真都是在個(gè)人筆記本電腦Windows 10系統(tǒng)的32位MATLAB R2012a下完成，其配置為：i5—6300HQ @ 2.30 GHz，4 GB RAM。

圖2分別為兩種建模方法在設(shè)計(jì)飛行點(diǎn)的巡航狀態(tài) 下模擬性能緩慢退化的健康參數(shù)估計(jì)結(jié)果。

圖2 2種建模方法下健康參數(shù)的估計(jì)結(jié)果

由圖可知，基于改進(jìn)擬合法建模的LKF隨著飛行循環(huán)數(shù)的增加，健康參數(shù)的估計(jì)誤差逐漸加大，主要原因是改進(jìn)擬合法只能在建模時(shí)選取的健康參數(shù)階躍量附近的小范圍內(nèi)模型誤差較小，其他退化情況下模型精度都難以保證。相比之下，自適應(yīng)擬合法建立的線性模型由于實(shí)時(shí)更新L,M矩陣，模型精度始終保持較高，因此對(duì)健康參數(shù)的估計(jì)誤差明顯小于前者。

本文健康參數(shù)的個(gè)數(shù)為8個(gè)，故改進(jìn)擬合法建模樣本的個(gè)數(shù)為8，而自適應(yīng)擬合法由于樣本個(gè)數(shù)的增加，模型精度得到提高，但增加了存儲(chǔ)空間、降低了機(jī)載實(shí)時(shí)性。圖3為自適應(yīng)擬合建模法不同樣本數(shù)下，對(duì)文獻(xiàn)[5]退化情況表的健康參數(shù)估計(jì)精度、存儲(chǔ)空間和計(jì)算時(shí)間的關(guān)系。圖中健康參數(shù)的估計(jì)精度用均方根誤差總和sum(RMSE)表示

(14)

圖3 樣本個(gè)數(shù)對(duì)估計(jì)精度、存儲(chǔ)空間和計(jì)算時(shí)間的影響

圖3中存儲(chǔ)空間為該設(shè)計(jì)點(diǎn)下從慢車到最大的離散13組系數(shù)矩陣以及對(duì)應(yīng)的健康參數(shù)組合引起的輸出偏差數(shù)組的總大小。此外，參考文獻(xiàn)[18]，該線性模型還加入BP神經(jīng)網(wǎng)絡(luò)(用于補(bǔ)償由于發(fā)動(dòng)機(jī)制造、安裝等造成的個(gè)體性能差異以及通過相似換算擴(kuò)展到其他飛行點(diǎn)后帶來的穩(wěn)態(tài)偏差)，故計(jì)算時(shí)間為個(gè)人筆記本上包含神經(jīng)網(wǎng)絡(luò)的線性模型的總運(yùn)行時(shí)間。本文折衷考慮，故選擇自適應(yīng)擬合法建模樣本的個(gè)數(shù)為4×4×5×4=320。

2.2 含氣路突變故障的仿真驗(yàn)證

考慮到發(fā)動(dòng)機(jī)個(gè)體性能退化的差異以及可能出現(xiàn)氣路突變故障，參考文獻(xiàn)[5]參數(shù)退化表模擬另外2種含氣路突變故障的性能退化情況，如表2，表中括號(hào)內(nèi)為在性能退化的基礎(chǔ)上在3 000飛行循環(huán)時(shí)注入的氣路突變故障的大小，故障參考文獻(xiàn)[6]，選取Level 1故障，以檢驗(yàn)突變故障的診斷能力。

表3 2種發(fā)動(dòng)機(jī)健康參數(shù)退化情況

自適應(yīng)擬合建模法對(duì)2種退化情況的估計(jì)結(jié)果見圖4(a1)和圖4(b1)。針對(duì)氣路突變故障的診斷，采用一定時(shí)間內(nèi)健康參數(shù)的變化量mk作為檢測依據(jù)

(15)

設(shè)定閾值，當(dāng)|mk|≥0.003時(shí)，說明發(fā)生氣路突變故障，如圖4(a2)和圖4(b2)，此時(shí)需要結(jié)合健康參數(shù)的估計(jì)結(jié)果進(jìn)行故障診斷。對(duì)于本文，一個(gè)飛行架次內(nèi)只采樣一組數(shù)據(jù)，而真實(shí)飛行時(shí)，閾值的設(shè)置需綜合考慮采樣周期、傳感器特性和卡爾曼增益矩陣大小。

圖4 自適應(yīng)擬合法建模的LKF估計(jì)結(jié)果

3 飛行區(qū)塊數(shù)字仿真及結(jié)果分析

由于真實(shí)的航空發(fā)動(dòng)機(jī)在每個(gè)飛行循環(huán)內(nèi)的巡航狀態(tài)時(shí)，其飛行高度H、飛行馬赫數(shù)Ma和風(fēng)扇轉(zhuǎn)速N1都不一定是設(shè)計(jì)值。因此對(duì)于設(shè)計(jì)飛行點(diǎn)附近的其他飛行狀態(tài)點(diǎn)先采用相似換算、再按照風(fēng)扇換算轉(zhuǎn)速插值系數(shù)矩陣和輸出偏差數(shù)組、然后用自適應(yīng)擬合法重新計(jì)算L,M矩陣的方法來估計(jì)氣路健康參數(shù)。為了確定相似原理的適用范圍，參考文獻(xiàn)[19]，計(jì)算全飛行包線內(nèi)各個(gè)飛行點(diǎn)距離設(shè)計(jì)飛行點(diǎn)的廣義距離

(16)

式中T2ds與P2ds為設(shè)計(jì)飛行點(diǎn)的風(fēng)扇進(jìn)口總壓和總溫。設(shè)定閾值d≤0.07，得到相似原理的適用范圍大致由以下4個(gè)飛行點(diǎn)：1)H=9 300 m,Ma=0.60;2)H=11.1 km,Ma=0.90;3)H=9.9 km,Ma=0.60;4)H=11.7 km,Ma=0.90組成的平行四邊形飛行區(qū)塊。

選擇該飛行區(qū)塊的4個(gè)頂點(diǎn)，每個(gè)飛行點(diǎn)下選取13個(gè)不同的風(fēng)扇換算轉(zhuǎn)速，每個(gè)飛行狀態(tài)都參考文獻(xiàn)[5]發(fā)動(dòng)機(jī)參數(shù)退化表的6 000飛行循環(huán)數(shù)模擬性能緩慢退化，記錄每個(gè)飛行狀態(tài)點(diǎn)的8個(gè)健康參數(shù)均方根誤差總和(sum(RMSE))，得到圖5。

圖5 2種建模方法下4個(gè)飛行點(diǎn)健康參數(shù)估計(jì)誤差總和

可知，在該飛行區(qū)塊4個(gè)頂點(diǎn)的各個(gè)風(fēng)扇換算轉(zhuǎn)速下，采用自適應(yīng)擬合法建模的線性卡爾曼濾波(linear Kalman filtering,LKF)相比改進(jìn)擬合法建模的LKF對(duì)健康參數(shù)的估計(jì)誤差明顯減小，尤其對(duì)于N1=95 %以上的狀態(tài)點(diǎn)，均方根誤差總和減小40 %以上。

4 結(jié)束語

1)相比改進(jìn)擬合的建模方法，基于此模型的卡爾曼濾波器對(duì)健康參數(shù)的估計(jì)誤差大大減小；

2)僅增加2個(gè)矩陣簡單計(jì)算，故實(shí)時(shí)性較高，數(shù)字仿真時(shí)間60 s下計(jì)算時(shí)間小于20 s；

3)通過采用相似換算，整個(gè)飛行區(qū)塊內(nèi)從慢車到最大的所有矩陣和數(shù)組的存儲(chǔ)空間小于200 kB，考慮到該發(fā)動(dòng)機(jī)為商用，其巡航包線較小，總存儲(chǔ)空間不大于1 MB。