王启航，黄金泉，鲁峰

(南京航空航天大学 能源与动力学院，江苏 南京 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种非线性测量噪声统计特性未知的情况下均显著减小了估计的方均根误差，提高了健康参数的估计精度和鲁棒性，实现了复杂噪声背景下航空发动机气路故障诊断。