曾　力,　龍　偉,　李炎炎

(1. 四川大學空天科學與工程學院, 四川 成都 610065; 2. 四川大學制造科學與工程學院, 四川 成都 610065)

燃氣輪機氣路主要由進氣道、導流葉片、壓氣機、燃燒室、渦輪等組成,這些部件在工作期間,受到高溫、高壓、高應力、吸入物、粉塵等的影響,性能不斷退化甚至出現部件的突發(fā)性失效,導致災難性后果[1-3].因此,如何對氣路部件的工作狀態(tài)進行監(jiān)控,評估燃氣輪機的健康狀況,為發(fā)動機維修提供合理的依據是當前研究的一個重點.通常健康狀況可用健康參數表示,例如風扇、高低壓壓氣機的效率系數和流量系數等,這些參數不可或難以直接測量.由于健康狀況的變化可以通過一些可測量參數(高低壓轉子轉速,壓氣機進氣口壓力、溫度,燃燒室進氣壓力、溫度等)的變化得以體現.因此,可以通過傳感器獲取這些可測量參數,建立這些參數到系統(tǒng)健康參數的映射模型,然后利用估計算法推導相應的健康參數.

無跡卡爾曼濾波器(unscented Kalman filter, UKF)是以卡爾曼濾波器作為框架,通過無跡變換(unscented transformation, UT)來逼近非線性函數的樣本和均值,沒有忽略高階項.UKF是對分布的概率密度函數的近似,避免求解Jacobian矩陣,其精度要高于擴展卡爾曼濾波器(extended Kalman filter, EKF).目前,UKF已應用于燃氣輪機氣路部件故障問題的診斷.Sebastien等[4]提出一種基于卡爾曼濾波與主成分分析的混合算法,通過把測量數據分解到主元空間和殘差空間,在兩個空間中分別構造兩個統(tǒng)計量,結合卡爾曼濾波實現對燃氣輪機性能的監(jiān)測與故障診斷.張鵬[5]針對計算誤差與噪聲引起的協(xié)方差矩陣負定的問題提出根據殘差的變化來改進濾波收斂速度的辦法.李冬等[6]提出通過主元分析對測量參數進行分解,提取反映性能衰減的成分,據此建立卡爾曼濾波模型.胡宇等[7]針對發(fā)動機氣路故障診斷模型建模不準確的情況,提出采用高斯回歸過程(Gauss process regression, GPR)對訓練數據進行建模,利用超球體單行采樣和平方根濾波來提高濾波的穩(wěn)定性.

燃氣輪機在工作中,隨著使用周期的增加導致的氣路部件老化或是突發(fā)性狀況,可能出現狀態(tài)突變的情況[8-9].這種情況下,普通的UKF無法根據測量值的變化準確估計出健康參數.針對這個問題,本文提出基于殘差相似性的漸消無跡卡爾曼濾波法(fading unscented Kalman filter based on residual similarity, FUKF-RS),通過強制殘差序列保持正交來確定漸消因子.針對遺忘因子ρ通過經驗值獲取而導致影響估值精度的問題,利用殘差陣的相似度代替遺忘因子計算誤差陣.另外,在構造濾波器時采用最小偏度采樣.此時采樣點個數為n+2(n為狀態(tài)空間維數),小于對稱采樣的樣本點數2n+1,縮短了運算時間.實驗證明此方案能夠在確保精度的情況下,增加系統(tǒng)的魯棒性,加強對突變狀態(tài)的跟蹤能力,其計算復雜度較低.

1　發(fā)動機健康參數估計算法設計

1.1　對突變狀態(tài)的估計

UKF算法將卡爾曼濾波的使用條件從線性場合推廣到非線性場合,并具有較高的估計精度.但是對于一些強非線性場合,例如被估計對象的狀態(tài)發(fā)生大范圍突變的時候,UKF的估值精度會發(fā)生大幅度下降.圖1描述某型燃氣輪機的高壓渦輪在不同時刻發(fā)生跳變時,利用UKF對渦輪流量系數進行估計得到的曲線.

圖1　基于UKF的流量系數突變時的估計值Fig.1　Estimation of the discharge coefficientduring state mutation based on UKF

圖1中,在第1 000 s和1 450 s時,高壓渦輪的狀態(tài)發(fā)生了突變,導致流量系數發(fā)生了2次跳變,但UKF的濾波曲線較為平滑,與真實值相比,無法準確跟蹤其值的變化.造成這種情況的主要原因是狀態(tài)的突變導致量測值殘差突然變化,而由于濾波器中當前數據所占比重較小,使得殘差變化值引起的濾波增益不夠大,無法通過估計值來體現出狀態(tài)的突變.利用漸消卡爾曼濾波(fading unscented Kalman filter, FUKF)算法可以提高濾波器對參數突變的跟蹤能力,其思想是給預測協(xié)方差乘以一個加權因子,通過改變加權因子的大小來調節(jié)歷史數據和當前數據所占的比重來達到跟蹤突變狀態(tài)的作用.將其應用在燃氣輪機氣路部件故障診斷分析中,利用濾波器由系統(tǒng)的測量參數估計出健康參數.

1.2　估值算法的設計

本文對健康參數的估計采用FUKF,結合燃氣輪機的結構點,算法如下[10]:

(1) 采集樣本點并確定其相應的權值

根據前面所述,按照最小偏度采樣的辦法采集樣本點,并計算其相應的權值,具體計算公式參見文獻[11].時刻k健康參數為χi,k,其中,i=1,2,…,10,代表第i個樣本.χi,k由高壓渦輪(h)、低壓渦輪(l)、風扇(f)、高壓壓氣機(p)的效率系數(下標e)和流量系數(下標d)組成,即χi,k=(he,hd,le,ld,fe,fd,pe,pd)T.

(2) 健康參數一步預測

基于燃氣輪機模型實現健康參數的一步預測,第i列健康參數在時刻k+1的預測值為

(1)

wk為燃氣輪機的系統(tǒng)噪聲.

(3) 健康參數均值和方差的計算

(2)

(3)

(4) 測量參數一步預測

基于燃氣輪機模型實現測量參數的一步預測,第i列測量參數在時刻k+1的預測值為

(4)

式中:f(·)為發(fā)動機的非線性測量模型;

vk為測量系統(tǒng)的噪聲.

(5) 測量參數均值和方差的計算

(5)

(6)

式中:Gk+1=diag (r1,(k+1),r2,(k+1),…,r8,(k+1))為漸消陣;

rg,(k+1)為第g個漸消因子,其值用于調節(jié)第g個測量參數的方差(g=1,2,…,8).

(6) 健康參數與測量參數間協(xié)方差的計算

時刻k+1健康參數與測量參數之間的協(xié)方差為

(7)

(7) 健康參數修正

健康參數修正主要包括時刻k+1卡爾曼增益Mk+1計算(如式(8))、健康參數估計值χk+1的修正以及協(xié)方差Pk+1的更新.

(8)

修正后的健康參數最優(yōu)估計值為

(9)

式中:yi,k+1為時刻k+1的傳感器實時測量值;

χk+1為健康參數一步預測值.

時刻k+1健康參數協(xié)方差陣修正值為

(10)

1.3　漸消因子求解

使FUKF穩(wěn)定工作的條件是在確定卡爾曼增益Mk+1時,滿足時刻k及時刻j測量參數殘差之間的協(xié)方差為0[9]，即

(11)

定義時刻k、j測量參數殘差的協(xié)方差為

(Pχy,k-Mk+1Ck),

(12)

式中：ψ(·)為健康參數最優(yōu)估計值隨采樣時刻的變化函數;

Pχy,k為時刻k的健康參數與測量參數之間的協(xié)方差；

Ck為時刻k測量參數殘差的方差.

通過式(11)強制每一步殘差序列保持正交可以確保序列中的信息全部被提取出來用于當前狀態(tài)估計,同時濾波器具有自適應跟蹤狀態(tài)變化的能力.當濾波器保持穩(wěn)定輸出時,此時增益最優(yōu),Ckj=0;如果系統(tǒng)模型誤差較大時,Ckj≠0,增益值并非最優(yōu),此時可通過調節(jié)漸消因子Gk+1使等式Pχy,k-Mk+1Ck=0成立,那么可認為增益最優(yōu),此時

(13)

式中:I為對角元素為1,其余元素為0的單位矩陣.

利用先驗知識可以確定漸消因子之間的比值為

r1,(k+1)∶r2,(k+1)∶…∶r8,(k+1)=a1∶a2∶…∶a8,

(14)

式中：ag為漸消因子的比值.

令ag=vgck+1，其中ck+1為比例公共因子.

根據前面描述的漸消無跡卡爾曼濾波器成立的條件及漸消陣Gk+1可得

Ck+1=ck+1diag (a1,a2,…,a8)×

(15)

式中:Rk+1為測量噪聲方差.

對式(15)兩邊求跡,整理可得

(16)

式中：tr(·)為矩陣求跡運算符.

從式(17)可得單個漸消因子的求解公式.

(17)

殘差陣由前后時刻殘差值按不同比例組成,如式(18)所示.

(18)

可以通過先驗知識確定出a1,a2,a3,…,a8之間的比值，計算出因子ck+1便可以得到漸消因子的值.ρ值越大,越能突出當前殘差的影響,ρ值越小,越突出歷史信息的影響.ρ的取值通常是通過人為經驗確定的,并沒有精確的量化取值辦法,過大或過小都會對狀態(tài)估計結果造成影響.另外,根據濾波器為殘差約束充分條件知道,基于殘差約束的卡爾曼濾波器的漸消因子必須保證前后殘差正交.如果系統(tǒng)狀態(tài)發(fā)生突變,則殘差陣將發(fā)生波動,無法確保正交.因此，本文提出利用前后估值時刻的殘差陣的對角元素來構造向量，通過前后時刻向量的相似度來表征殘差陣之間的相似性，并用該相似度值來替代ρ.

2　基于殘差相似性的誤差陣求解

λt=dk,vk+1/(|dk|·|vk+1|).

(19)

根據突出當前信息,兼顧歷史數據的原則,Ck+1改為

(20)

式(20)確保了當前信息前的系數始終大于歷史數據的系數.

3　實驗仿真

針對某型航空燃氣輪機氣路部件性能漸變性退化和突變性退化,通過改變燃氣輪機健康參數的值來模擬實際工作中出現的故障.通過漸變和突變的方式,改變健康參數,從而達到模擬部件性能退化的目的,分別利用FUKF和FUKF-RS對健康參數進行估計,通過變化量的形式對健康參數的退化值進行表示.另外,為了與圖1對比,分別利用FUKF和FUKF-RS對圖1條件下的狀態(tài)突變進行估計.

仿真環(huán)境為高度H=0,馬赫數為0,分別施加系統(tǒng)噪聲Q和測量噪聲R,分別為0. 0022In×n和0. 0022Im×m.在初始時刻，樣本點取值均為1，即he=hd=Le=Ld=fe=fd=Pe=Pd=1.設定整體仿真計算時間為3 600 s,每隔60 s采集一次傳感器的測量值,同時估計此時對應的健康參數值.仿真中根據故障原理設定風扇、壓氣機的效率系數和流量系數均下降,高、低壓渦輪的效率系數下降,流量系數升高.其中,fe、fd分別下降3%和3.5%,pe、pd分別下降2.5%和3%,he、he分別下降3%和4%,hd、ld均上升2%.本文只列出高壓壓氣機的效率系數pe與流量系數pd的估計值，利用FUKF和FUKF-RS對pe與pd估值分別如圖2、3所示.

圖3　基于FUKF-RS估計高壓壓氣機的健康參數Fig.3　Estimation of the health parameter of high pressure compressor based on FUKF-RS

由圖2可以看出，FUKF能估計出高壓壓氣機健康參數的變化趨勢,但與理論值相比,存在較大的誤差,這是由于遺忘因子的取值不佳造成的.

相比圖2、3的估值更接近于理論值.這是由于FUKF-RS使當前采樣時刻測量參數殘差的信息在進行測量值協(xié)方差更新過程中所占比例較大,更為接近實際測量工況.

針對圖1中高壓渦輪流量系數突變時估值不準的問題，分別利用FUKF和FUKF-RS對該工況下的流量系數估值，如圖4、5所示.

圖4　基于FUKF-RS的高壓渦輪流量系數突變時的估計值Fig.4　Estimation of the discharge coefficient of high pressure turbine during state mutation based on FUKF-RS

由圖4可以看出,FUKF-RS具有較強的狀態(tài)突變跟蹤能力,相比圖1的結果有了較大的提高.

圖5　基于FUKF的高壓渦輪流量系數突變時的估計值Fig.5　Estimation of the discharge coefficient of high pressure turbine during state mutation based on FUKF

由圖5可以看出,FUKF同樣具有一定的突變跟蹤能力,但相比FUKF-RS,其跟蹤能力較低.

穩(wěn)態(tài)時，FUKF-RS的狀態(tài)估計精度較高，狀態(tài)突變時，有較強的狀態(tài)跟蹤能力.主要是因為通過強制每一步估計的殘差正交,確保殘差序列中的信息被完全提取出來，這樣使得濾波器具有對狀態(tài)突變的跟蹤能力.同時，利用表征殘差相似性的參數來替代漸消算法中的ρ，避免了因為ρ取值不當造成的估值誤差.

4　結　論

針對燃氣輪機在長期工作過程中由于不斷退化導致部件性能突變,引起健康參數在極短時間內出現突變,導致難以準確估計健康參數值的問題.

(1) 本文采用FUKF對健康參數估計,針對FUKF中求解漸消因子時存在的遺忘因子取值不當可能影響健康參數的估值精度的問題,提出FUKF-RS算法.

(2) 利用采樣前后時刻測量值協(xié)方差矩陣的對角元素所組成向量的余弦值來表征向量之間的相似度,利用該相似度替代遺忘因子,實現了求解過程中的量化取值.

(3) 針對發(fā)動機部件出現故障,健康參數出現較大范圍波動時,FUKF-RS相比FUKF,其參數跟蹤精度提高了3%左右.而UKF則無法對健康參數進行準確估計.

(4) 發(fā)動機在正常工作,健康參數緩慢退化的情況下,FUKF-RS的估值精度相比FUKF提高了2%左右,且估值曲線更平滑.

致謝:四川大學實驗技術立項項目資助(20170135).