基于變分貝葉斯粒子濾波的航空發動機氣路故障診斷方法

王啟航，黃金泉，魯峰

(南京航空航天大學 能源與動力學院，江蘇 南京 210016)

0 引言

航空發動機的性能與可靠性是飛機性能與飛行安全的重要保證。為了降低維護費用，提高飛行安全，必須及時、準確地評估部件的健康狀況，實現發動機故障診斷[1-2]。發動機氣路部件故障占發動機總體故障的90%以上，主要表現在部件的效率、流量等性能參數的變化，其變化會直接導致發動機轉速、溫度、壓力等測量參數的變化。因此發動機氣路故障診斷可以通過傳感器測量參數的變化對健康參數進行估計，分析發動機氣路部件的健康狀況[3]。

目前發動機氣路部件故障診斷方法中基于模型的方法應用最為廣泛，其中卡爾曼濾波(KF)算法簡單，但只適用于線性系統；擴展卡爾曼濾波(EKF)和無跡卡爾曼濾波(UKF)可以解決非線性濾波問題，但三種方法均不適用于非高斯系統，而實際發動機系統的過程噪聲與測量噪聲都比較復雜，不一定滿足高斯分布假設[4]。

粒子濾波(PF)是一種適用于非高斯噪聲的非線性系統狀態估計的濾波方法。粒子濾波通過一組加權粒子的遞推來給出系統狀態后驗概率密度近似形式，理論上對過程噪聲和測量噪聲統計特性無任何限制，是一種適用性極強的貝葉斯濾波方法[5]。粒子濾波精度很大程度上依賴于模型的統計特性，但在實際應用中，測量噪聲統計特性一般是未知且時變的，并表現出非線性的特點。當假設的噪聲特性與實際不符時，PF產生誤差積累，導致狀態估計性能下降，難以滿足測量噪聲特性未知且時變條件下發動機氣路故障診斷對PF濾波精度的要求[6]。

針對以上問題，本文提出基于變分貝葉斯(VB)和核極限學習機(KELM)的邊緣化粒子濾波(VB-KELM-MPF)方法，用學生t分布對測量噪聲建模，引入邊緣化粒子濾波的思想，利用VB對其進行遞推求解得到噪聲參數，同時利用遞推出的測量噪聲和KELM求取測量似然函數用來更新粒子權值，實現非高斯測量噪聲未知參數和發動機健康參數的聯合估計。

1 粒子濾波基本原理

粒子濾波算法由GORDON N J等[5]于1993年提出，其基本思想是從重要性概率密度中抽取隨機樣本，根據測量信息不斷調整權值大小和粒子位置，最后用粒子的加權和代替積分運算，獲得狀態最小方差估計。假設非線性系統的狀態空間模型為：

xk=f(xk-1,uk-1)+ωk-1

(1)

yk=h(xk,uk)+υk

(2)

濾波問題就是求解后驗概率密度函數p(xk|y1:k)。根據最優貝葉斯估計理論有：

狀態預測方程

(3)

狀態更新方程

(4)

(5)

式中δ(·)是Dirac delta函數。

從一個易采樣的重要性密度函數q(xk|xk-1,yk)中抽取一組加權樣本，權值通過重要性采樣法選擇：

(6)

標準粒子濾波算法選擇最易于實現的先驗概率密度作為重要密度函數，即

(7)

則重要性權值化簡為

(8)

將權值歸一化，后驗概率密度可表示為

(9)

可見，當N→∞時，由大數定理即可保證上式逼近真實后驗概率p(xk|y1:k)。

2 非高斯測量噪聲統計特性未知情況下的改進粒子濾波算法

在實際工程應用中，測量數據可能含有非高斯噪聲，甚至出現一定數量的野值，表現出重尾特性。此時，僅用高斯分布來描述噪聲特性的方法將不再適用。因此，測量噪聲采用非高斯學生t分布代替高斯分布，結合VB及KELM實現對噪聲參數和系統狀態的聯合遞推估計。

2.1 非高斯測量噪聲模型建立

對于任一狀態空間模型，假設未知的測量噪聲服從學生t分布：

υk～St(yk|μk,Λk,dfk)=

(10)

式中：μk為均值；尺度矩陣Λk為對角陣Λk=diag(Λk,1,Λk,2,…,Λk,nυ)，dfk為學生t分布的自由度。則測量方程式(2)的統計描述為

yk～p(yk|xk)=St(yk|h(xk),μk,Λk,dfk)

(11)

學生t分布概率密度是一簇關于μk對稱的曲線，dfk的大小決定了曲線的形狀，反映了學生t分布重尾特性的程度。不同自由度下的學生t分布曲線如圖1所示。

圖1 不同自由度下的學生t分布曲線

式(10)不可解析求得，所以給出在機器學習中更為普遍的學生t分布的概率密度形式：

St(yk|h(xk),μk,Λk,uk,dfk)=

(12)

式中uk為引入的服從參數為伽馬分布G(uk|dfk/2,dfk/2)的隱變量。

2.2 噪聲未知參數的變分貝葉斯迭代學習

系統狀態、測量參數和未知噪聲參數的聯合概率密度分布為

p(xk,μk,Λk,uk,dfk|y1:k)∝p(xk|y1:k-1)·

p(yk|xk,μk,Λk,uk)p(μk|Λk)p(Λk)p(uk|dfk)p(dfk)

(13)

考慮到聯合共軛先驗可以表示為高斯分布和伽馬分布的乘積，式(13)右側后驗概率分別表示為：

p(yk|xk,μk,Λk,uk)=N(yk|h(xk)+μk,ukΛk)

(14)

G(Λk,j|ck|k-1,j,dk|k-1,j)}

(15)

p(uk|dfk)=G(uk|dfk/2,dfk/2)

(16)

p(dfk)=G(dfk|ak|k-1,bk|k-1)

(17)

式中：ξk|k-1,j和(βk|k-1Λk,j)-1為高斯分布的均值和方差；βk|k-1為方差系數；a、b、c、d是伽馬分布的參數。

在變分貝葉斯框架下，以上各個后驗分布的變分近似后驗都可以通過計算式(13)的自然對數形式獲得，從而求得學生t分布的未知參數。

首先，由式(14)和式(15)可得均值μk的邊緣VB近似后驗的對數形式：

lnq(μk|Λk)=

Euk,dfk[lnp(yk|xk,μk,Λk,uk)+lnp(μk|Λk)]+

(18)

式中：Ex(·)為概率分布對變量的期望；C為常數。q(μk|Λk)服從N(μk|ξk|k,(βk|kΛk|k)-1)，式 (18)兩邊取指數，得到第一組未知參數的遞推公式：

βk=βk|k-1+E[uk]

(19)

ξk=ξk|k-1+E[uk](βk)-1(yk-h(xk)-ξk|k-1)

(20)

根據貝葉斯原理p(x,y)=p(x|y)p(y)，聯立式(14)、式(15)和式(18)，得到尺度矩陣Λk的VB近似后驗的對數形式：

E[uk]Λk(yk-h(xk)-μk)-

(21)

由共軛先驗的特性可知，q(Λk)服從G(Λk,j|ck|k,j,dk|k,j)，將式(19)和式(20)代入式(21)并兩邊取指數，得第二組未知參數的遞推公式：

(22)

[(yk-h(xk)-ξk|k-1)(yk-h(xk)-ξk|k-1)T]jj

(23)

式中等式右邊最后一項下標表示方陣的主對角元。

類似地，利用VB理論可以推導出自由度dfk和隱變量uk的近似后驗的對數形式。結合式(16)和式(17)，對dfk的VB近似后驗的對數形式兩邊取指數，得第三組未知參數的遞推公式：

(24)

(25)

結合式(14)和式(16)，對uk的VB近似后驗的對數形式兩邊取指數，得第四組未知參數的遞推公式：

(26)

(27)

根據高斯分布和伽馬分布的期望公式可以求出以上各式中的期望值：

(28)

式中Ψ(·)=Γ′(·)/Γ(·)為Digamma函數。

考慮到噪聲參數時變的情況，引入漸消因子ρ∈(0,1)反映噪聲波動，得到各個參數的預測值：

ξk|k-1=ρξk-1,βk|k-1=ρβk-1
ak|k-1=ρak-1,bk|k-1=ρbk-1
ck|k-1,j=ρck-1,j,dk|k-1,j=ρdk-1,j

(29)

2.3 基于VB和KELM的邊緣化粒子濾波算法

本文將統計特性未知的噪聲作為非線性系統的子結構，對其進行邊緣化處理，并利用上節中VB遞推公式近似逼近噪聲參數的后驗分布p(Θk|xk,y1:k)，其中Θk={μk,Λk,dfk}為噪聲參數集。噪聲參數集和系統狀態的聯合后驗概率密度為

p(xk,Θk|y1:k)=p(Θk|xk,y1:k)p(xk|y1:k)

(30)

因為噪聲被邊緣化，p(xk|y1:k)仍通過重要性采樣求取，此時要求每個粒子都要利用上節中VB遞推獲得每個粒子的充分統計量。因此，聯合后驗概率密度可通過如下形式求解：

(31)

式中N為粒子數。粒子的權值為

(32)

VB遞推過程中每一次迭代都要計算N次模型以獲得測量似然分布p(yk|xk)，對于發動機非線性模型而言，該過程將大大增加計算量，導致濾波響應緩慢。針對該問題，本文在VB迭代計算未知噪聲參數時，根據當前時刻粒子狀態和對應的測量值利用KELM建立數據驅動代替發動機模型計算似然分布p(yk|xk)，從而大幅度提高噪聲參數估計速度。本文不再贅述KELM原理，可參考相關文獻[7]。

非高斯噪聲統計特性未知的情況下VB-KELM-MPF算法具體步驟如下。

2)根據式(29)對每個粒子的噪聲參數進行預測。

4)噪聲參數遞推更新。根據式(19)、式(20)、式(22)-式(28)來計算學生t分布的未知參數，每次迭代利用KELM計算p(yk|xk)，直到迭代收斂。

3 基于VB-KELM-MPF的航空發動機氣路故障診斷方法

本文以某型小涵道比渦扇發動機為研究對象驗證VB-KELM-MPF算法的有效性和魯棒性。基于粒子濾波算法的氣路部件故障診斷原理如圖2所示。診斷原理實質就是通過真實發動機的輸出值與模型預測值之間的殘差對部件性能參數變化進行估計，將發動機氣路部件故障診斷轉化為部件性能參數變化的狀態估計問題。

圖2 發動機氣路部件故障診斷原理圖

3.1 航空發動機氣路故障模型

(33)

式中：i=1、2、3、4，分別表示風扇、壓氣機、高低壓渦輪；ηi和Wi分別表示部件效率和流量的實際特性；ηi,d和Wi,d分別表示部件效率和流量的理想特性。SEi,SWi∈(0,1]，取1時表示沒有故障的理想健康狀態。測量參數y=[nL,nH,T22,P22,T3,P3,T43,P43,T5,P5]T，依次為低壓轉子轉速、高壓轉子轉速、風扇出口總溫總壓、壓氣機出口總溫總壓、高壓渦輪出口總溫總壓和低壓渦輪出口總溫總壓。

3.2 仿真與分析

仿真實驗中，設過程噪聲ωx～N(0,Q1)，ωh～N(0,Q2)，其中Q1=0.001 52I2×2，Q2=0.0042I8×8。T5測量噪聲形式：1)0均值高斯分布N(0,0.0032)疊加20%服從均勻分布U(-0.1,0.1)的野值；2)第6s開始噪聲變為均值不為0的高斯分布N(0.01,0.0032)，整個仿真過程疊加20%服從均勻分布U(-0.1,0.1)的野值；3)伽馬分布0.0032×Γ(0.25,0.5)疊加20%服從均勻分布U(-0.1,0.1)的野值。粒子數N=50，仿真步數M=1 000，發動機采樣步長設定為0.02s，在第2s注入故障。PF作為比較算法取R=0.0032I10×10保持不變。

圖3 N(0, 0.0032)+20%U(-0.1, 0.1)噪聲分布下健康參數估計結果

圖4 N(0.03, 0.0032)+20%U(-0.1, 0.1)噪聲分布下健康參數估計結果

圖5 0.0032×Γ(0.25, 0.5)+20%U(-0.1, 0.1)噪聲分布下健康參數估計結果

表1 四種噪聲模式下兩種算法的RMSE

由圖3可以看出，在均值為0、含有野值的測量噪聲情況下，PF的估計結果在野值出現的時刻波動較大；圖4和圖5顯示，由于假設噪聲統計特性與實際測量噪聲不符，導致PF對健康參數的估計出現偏差，甚至發散。由表1可以看出，與已知0均值高斯噪聲分布的診斷結果相比，在非高斯噪聲未知的情況下，PF估計精度均有所下降。VB-KELM-MPF用具有重尾特性的學生t分布建立測量噪聲模型，實時遞推分布中反映均值的參數ξ以及自適應調整自由度df。仿真結果表明：該方法在非0均值噪聲、伽馬噪聲、有野值且噪聲特性未知的情況下均能較準確地估計出健康參數，濾波精度維持在較高的水平，實現了復雜噪聲背景下的發動機故障診斷。

4 結語

本文針對測量噪聲統計特性未知且時變情況下航空發動機氣路故障診斷精度下降甚至發散的問題，提出了基于變分貝葉斯和核極限學習機的粒子濾波算法。該方法對非高斯測量噪聲進行邊緣化處理，將VB原理和KELM方法相結合遞推出噪聲的統計特性，并估計出健康參數。仿真結果表明：相比于標準粒子濾波，本文提出的方法在3種非線性測量噪聲統計特性未知的情況下均顯著減小了估計的方均根誤差，提高了健康參數的估計精度和魯棒性，實現了復雜噪聲背景下航空發動機氣路故障診斷。