CN104166787A - 一种基于多阶段信息融合的航空发动机剩余寿命预测方法 - Google Patents

Abstract

本发明基于多阶段信息融合的航空发动机剩余寿命预测方法包括步骤：多源监测参数去噪处理与特征提取；对多源监测时间序列进行平稳性分析，计算各个参数监测时间序列突变点，计算突变点处参数退化比例；对多源参数进行多阶段划分，建立回归融合模型，利用历史监测数据进行样本训练，获得融合模型在多阶段内的参数；结合训练集中的监测数据，融合多源监测参数，得到健康指标HI；利用Kalman滤波算法，对发动机从性能完好到性能失效全过程进行最佳拟合，并最小化预测模型的误差；结合测试集中的实时监测数据，融合多源监测参数，得到健康指标HI；利用Kalman滤波算法，对预测模型时变参数进行实时估计；确定预测模型，引入时间机制，实时估计发动机的失效时间。

Description

一种基于多阶段信息融合的航空发动机剩余寿命预测方法

技术领域：

本发明涉及发动机剩余寿命预测方法，特别涉及一种基于多阶段信息融合的航空发动机剩余寿命预测方法。

背景技术：

作为飞机的核心部件，航空发动机的结构极其复杂，因此，单一监测参数并不能准确表征其性能状况。为综合利用观测信息以准确描述发动机的实时健康状况，信息融合技术在发动机的健康管理中得到了广泛的应用，例如基于信息融合的故障诊断、性能评估以及性能趋势预测等。其中，衰退趋势的预测对发动机的预知维修决策尤为重要，是目前计划维修转向预知维修过程中的研究热点。

信息融合技术是用数学方法和技术工具综合不同源信息，从而得到高品质的有用信息。因此，存在各种不同种类、不同等级的融合，如数据融合、图像融合、特征融合、决策融合、传感器融合、分类器融合等。随着计算机技术和网络通信技术的飞速发展，信息融合的应用领域不断扩展,从开始诞生的军事领域，逐渐向其它领域渗透，如：智能机器人与智能车辆领域、医学图象处理与诊断、气象预报、地球科学、农业应用领域、现代制造领域和经济商业领域等。在发动机的健康管理中，信息融合技术显得尤为重要，其功能表现在：充分利用信息的冗余性与互补性，提高时间或空间的分辨率，增加目标特征矢量的维数，降低信息的不确定性，改善信息的置信度，从而降低推理的模糊程度，提高系统决策能力。

在电子设备、电气设备、控制设备、航空发动机等复杂设备健康管理中，以健康状态预测为核心的健康管理技术是复杂装备健康管理的核心技术。对民航发动机的健康状态预测主要是对发动机性能下降阶段的健康状态预测。对民航发动机衰退趋势的准确预测，可以制定更加有效的维修计划，能够在保证发动机可用性和可靠性的前提下有效降低运行成本，实现安全与经济的双赢。从方法体系的角度可以将现有的民航发动机健康状态变化趋势预测方法划分为模型驱动和数据驱动两种方法。前者往往需要基于发动机工作原理、材料属性、工况等因素建立解析模型，其建模代价一般比较高。后者则无需先验假设，可以直接基于发动机的健康状态数据进行建模，其建模代价相对较小。民航发动机状态监视技术的发展为数据驱动建模提供了丰富的数据，使得数据驱动建模预测方法具有良好的应用前景。

常见的基于融合进行发动机状态预测的过程大多没有考虑发动机各个参数与发动机整体性能衰退的突变性。因此，已有研究并未对发动机的衰退过程做分阶段处理，从而在融合中不能准确反映发动机的性能衰退非平稳性，预测模型的样本训练也会因为单阶段处理引起参数估计不准确的问题。

发明内容：

为解决常用发动机剩余寿命预测方法所不能解决的问题，本发明对现有发动机健康状态预测技术进行改进，实现多源数据的分阶段融合，基于多阶段信息融合结果建立预测模型并利用多阶段训练样本估计模型的参数，引入实时监测数据，采用智能算法实现多阶段模型参数的更新与预测，最终预测发动机的性能衰退趋势，得到剩余使用寿命的准确估计。

本发明采用如下技术方案：一种基于多阶段信息融合的航空发动机剩余寿命预测方法，其包括如下步骤：

(a)多源监测参数去噪处理与特征提取；

(b)在步骤(a)的基础上，对多源监测时间序列进行平稳性分析，计算各个参数监测时间序列突变点，计算突变点处参数退化比例；

(c)在步骤(b)的基础上，对多源参数进行多阶段划分，建立回归融合模型，利用历史监测数据进行样本训练，获得融合模型在多阶段内的参数；

(d)在步骤(c)的基础上，结合训练集中的监测数据，融合多源监测参数，得到健康指标HI；

(e)在步骤(d)的基础上，利用Kalman滤波算法，对发动机从性能完好到性能失效全过程进行最佳拟合，以确定健康指标趋势预测模型结构，并最小化预测模型的误差；

(f)在步骤(c)的基础上，结合测试集中的实时监测数据，融合多源监测参数，得到健康指标HI，此时的HI是预测对象；

(g)在步骤(f)的基础上，利用Kalman滤波算法，对预测模型的时变参数进行实时估计；

(h)在步骤(e)和步骤(g)的基础上，确定预测模型，引入时间机制，从而实时估计发动机的失效时间。

本发明具有如下有益效果：本发明所采用的基于多阶段信息融合的航空发动机剩余寿命预测方法是一种数据驱动方法，从发动机的多源监测参数中获得退化轨迹，而不需要对发动机性能衰退趋势的形态进行假设。多阶段的融合过程不仅能够实现数据的融合而且能够实现特征的融合，从而切实表征发动机性能退化的非平稳性或突变性。因此，本发明能够利用较为全面的监测信息，基于发动机的多阶段性能退化特征，借助人工智能的算法获得满足要求的预测结果，具备实际应用价值。

附图说明：

图1为本发明的结构框架示意图。

图2为本发明研究实体在全寿命周期内的低压压缩机出口总温(T24)监测时间序列。

图3为发动机多阶段融合步骤示意图。

图4为发动机全寿命周期内的健康指标HI的退化趋势。

图5为发动机寿命预测步骤示意图。

图6为单阶段融合与多阶段融合的结果比较示意图。

图7为训练集下融合得到的HI时间序列。

图8为训练集下的HI预测过程。

具体实施方式：

请参照图1所示，本发明基于多阶段信息融合的航空发动机剩余寿命预测方法主要包含信息融合与衰退趋势预测两大模块。信息融合模块包括：数据处理，用于去除工频噪声、背景噪声以及随机脉冲干扰等因素对纯净信号的影响；平稳性分析，用于分析监测时间序列的平稳性，并准确定位引起性能衰退非平稳的突变点；阶段划分与特征提取，用于对非平稳的监测时间序列进行阶段划分，实现阶段内的平稳性，从而达到对时间序列的阶段线性化处理；多阶段分析，用于多源数据间的关联分析，建立多源的、具有阶段平稳性的监测数据与发动机隐含的健康状况间的对应关系。衰退趋势预测包含：模型建立，在多阶段信息融合的基础上，采用智能算法基于历史性能数据构建性能衰退趋势预测模型；多阶段参数估计，采用分阶段样本训练的方法估计预测模型的参数；参数更新与预测，结合实时监测数据，采用智能算法对模型参数进行更新与预测，最后计算发动机剩余使用寿命。

请参照图1至图5所示，本发明基于多阶段信息融合的航空发动机剩余寿命预测方法的具体步骤如下：

(a)多源监测参数去噪处理与特征提取。

(b)在步骤(a)的基础上，对多源监测时间序列进行平稳性分析，计算各个参数监测时间序列突变点，计算突变点处参数退化比例。

(c)在步骤(b)的基础上，对多源参数进行多阶段划分，建立回归融合模型，利用历史监测数据进行样本训练，获得融合模型在多阶段内的参数。

(d)在步骤(c)的基础上，结合训练集中的监测数据，融合多源监测参数，得到健康指标HI。

(e)在步骤(d)的基础上，利用Kalman滤波算法，对发动机从性能完好到性能失效全过程进行最佳拟合，以确定健康指标趋势预测模型结构，并最小化预测模型的误差。

(f)在步骤(c)的基础上，结合测试集中的实时监测数据，融合多源监测参数，得到健康指标HI，此时的HI是预测对象。

(g)在步骤(f)的基础上，利用Kalman滤波算法，对预测模型的时变参数进行实时估计。

(h)在步骤(e)和步骤(g)的基础上，确定预测模型，引入时间机制，从而实时估计发动机的失效时间。

步骤(a)中，进行多源监测参数的去噪处理。

由于多源监测数据具有非平稳非线性的特征，采用经验模态分解方法(EmpiricalMode Decomposition,EMD)对多源监测数据进行去噪处理，EDM将原始监测数据分解为从高频到低频的多个本征模函数(Intrinsic Mode Function,IMF)与残余项之和。其中，IMF还包括信号起主导作用的模态与噪声起主导作用的模态，只需从IMF中剔除噪声起主导作用的模态，将信号起主导作用的模态与残余项进行叠加，即可进行部分重构得到降噪之后的监测时间序列。EMD分解的公式如下：

$z\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{IMF}}_i+r_n---\left(1\right)$

其中，z(t)为监测时间序列，IMF₁,IMF₂,...,IMF_n分别代表从高到低不同频率段的信号成分，残余项r_n代表信号的平均趋势。

分解过程中，IMF必须满足以下条件:

(1)在整个时间范围内，极值点与过零点的数目相等或最多相差一个；

(2)所有极大值形成的上包络线与所有极小值形成的下包络线的均值始终为0。EMD分解步骤：

(1)找出z(t)的所有局部极大值，对该极值进行三次样条插值，得到由所有极大值点构成的上包络线，记为a(t)；

(2)找出z(t)的所有局部极小值，对该极值进行三次样条插值，得到由所有极小值点构成的下包络线，记为b(t)；

(3)计算上下包络线的均值记为c(t)，

(4)计算d(t)＝z(t)-c(t)，当d(t)满足上述IMF的两条性质时，d(t)即第一个IMF分量；否则，以d(t)作为输入量，重复步骤(1)～(3)，直到第一个IMF产生，记第一个IMF分量为IMF₁(t)；

(5)计算r₁(t)＝z(t)-IMF₁(t)，将r₁(t)作为输入量，重复步骤(1)～(4)，得到第二个IMF分量IMF₂(t)；计算r₂(t)＝z(t)-IMF₂(t)，重复上述步骤，直到余项r_n(t)为一个单调序列或小于给定阈值，EMD分解即结束。

分解结束之后，需根据最小连续均方误差(consecutive mean square error,CMSE)的准则区分IMF各个分量中哪些部分对信号起主导作用，哪些部分对噪声其主导作用。CMSE计算公式如下：

$\mathrm{CMSE}({\hat{z}}_k,{\hat{z}}_{k+1})=\frac{1}{N}\underset{i}{\overset{N}{\mathrm{\Σ}}}{[{\hat{z}}_k\left(t_i\right)-{\hat{z}}_{k+1}\left(t_i\right)]}^2=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left[{\mathrm{IMF}}_k\left(t_i\right)\right]}^2---\left(2\right)$

其中，为k项相加重构出的去噪监测时间序列，包含k-1个IMF分量和一个余项。当CMSE的值最小时，说明在k处IMF分量出现转折，而IMF首次转折的地方即为噪声起主导作用模态与信号起主导作用模态的分界，去噪之后的监测时间序列由以下公式进行部分重构：IMF_k+1(t)+IMF_k+2(t)+...+IMF_n-1(t)+r_n(t)。

步骤(b)中，对多源监测时间序列进行平稳性分析，计算各个参数监测时间序列突变点，计算突变点处参数退化比例。

(1)采用以下启发式算法对时间序列的平稳性进行分析并划分时间窗：

$\left\{\begin{array}{c}{S}_D={\left[\frac{(N_1-1)\×S_1{\left(i\right)}^2+(N_2-1)\×S_2{\left(i\right)}^2}{N_1+N_2-2}\right]}^{\frac{1}{2}}\×[\frac{1}{N_1}+\frac{1}{N_2}]\\ T\left(i\right)=\left|\frac{u_1\left(i\right)-u_2\left(i\right)}{S_D\left(i\right)}\right|\end{array}\right.---\left(3\right)$

其中，N₁、N₂分别表示i点左右两部分的时间点个数；u₁(i)、u₂(i)分别表示i点左右两部分的均值；S₁(i)、S₂(i)分别表示i点左右两部分的标准偏差；S_D(i)为合并偏差；T(i)为检验统计值，T(i)越大，表明该点左右两部份统计特性的差异越大。该算法通过计算T(i)将非平稳时间序列每个时间点的左右两部分进行比较。T(i)值最大时对应的i点即为该非平稳时间序列的分割点，可将整个时间窗一分为二。当第一个分割点获得之后，可用同样的方法对划分好的两个时间窗分别进行分割即可获得四个时间窗，以此类推。时间窗的数量应当依据实际应用情况和子窗内时间序列的平稳性而定。

(2)划分时间窗的点可作为时间序列的分割点，但要判断分割点是否为性能突变点，还需要计算慢变量：

$S_V=\frac{\mathrm{\ΔS}}{\mathrm{\Δt}}=\frac{S_j-S_i}{t_j-t_i}---\left(4\right)$

其中，S_V为慢变量，Δt为时间间隔，j＞i，i,j＝1,2,...,n，n为时间点的个数，ΔS为监测值在相应时间间隔内的变化量。理论上，时间间隔越小越能反映详细的慢变量变化过程，但是由于噪声影响，该过程的波动幅度会随着时间间隔的缩小而变大。考虑到最终要获取的是慢变量的相对变化(从一个稳态到另一个稳态的转变)，可适当调整时间间隔。通过时间窗口Δt的滑动可求得整个时间窗内的慢变量序列。慢变量的突变即体现参数退化的突变点。在图4给出的时间序列中，分割点有3个，但通过计算慢变量可确定其第二个分割点为该时间序列的实际突变点。

(3)各个监测参数时间序列的突变点确定后，即可计算监测参数平均突变点与平均变化比例：

$P_{\mathrm{An}}=\frac{1}{l}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{P}_{\mathrm{in}}---\left(5\right)$

$m_{\mathrm{An}}=\frac{1}{l}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{m}_{\mathrm{in}}---\left(6\right)$

其中，P_An为第n个平均突变点；m_An为第n个平均突变点处的平均变化比例；P_in为第i个监测参数的第n个突变点，m_in为第i个监测参数从退化初到第n个突变点的变化量在整个退化过程中的变化量的比例，i＝1,2,...,l，l为监测参数的个数，n为突变点的个数。m_in的公式如下所示：

$\begin{array}{c}{m}_{\mathrm{in}}=(\frac{1}{5}\underset{t=1}{\overset{5}{\mathrm{\Σ}}}{S}_{i\_t}-\frac{1}{5}\underset{t=n-2}{\overset{n+2}{\mathrm{\Σ}}}{S}_{i\_t})/(\frac{1}{5}\underset{t=1}{\overset{5}{\mathrm{\Σ}}}{S}_{i\_t}-\frac{1}{5}\underset{t=N-4}{\overset{N}{\mathrm{\Σ}}}{S}_{i\_t})\\ =(\underset{t=1}{\overset{5}{\mathrm{\Σ}}}{S}_{i\_t}-\underset{r=n-2}{\overset{n+2}{\mathrm{\Σ}}}{S}_{i\_t})/(\underset{t=1}{\overset{5}{\mathrm{\Σ}}}{S}_{i\_t}-\underset{t=N-4}{\overset{N}{\mathrm{\Σ}}}{S}_{i\_t})\end{array}---\left(7\right)$

其中，是第i个监测参数在退化初期的5个循环时的监测值的平均值，是第i个监测参数在第n个突变点附近5个循环时的监测值的平均值，是第i个监测参数在性能退化末期的5个循环时的监测值的平均值，N是第i个监测参数的监测值个数。

多组监测时间序列的突变点不严格在一个点上，但是由于监测参数退化突变的相似性，例如，前半段退化缓慢后半段退化迅速，各个监测参数的突变点会集中在某个点附近，因此将各个监测时间序列突变点的平均值作为融合时间序列的突变点。这样才能确保多源监测时间序列的突变性与融合时间序列(HI时间序列)的突变性相一致。例如(基于图4的数据)，计算得到突变点在第217个循环时附近，退化量在整个退化过程中的比例是0.507，说明在突变点附近的HI应为0.493。

步骤(c)中，对多源参数进行多阶段划分，建立回归融合模型，利用历史监测数据进行样本训练，获得融合模型在多阶段内的参数。

(1)基于步骤(b)获得的突变点，建立分阶段融合模型：

$\left\{\begin{array}{c}y1=a_{10}+{\mathrm{XA}}_1^T\\ y2=a_{20}+{\mathrm{XA}}_2^T\\ .\\ .\\ .\\ \mathrm{yn}=a_{n0}+{\mathrm{XA}}_n^T\end{array}\right.---\left(8\right)$

其中，y1表示第1个突变点前的HI；y2表示第1个突变点到第2个突变点间的HI；yn表示第n-1个突变点到第n个突变点间的HI；X表示多源监测参数(x₁,x₂,...,x_l)，l表示性能参数的个数；a₁₀,表示第1个突变点前的模型系数A₁＝(a₁₁,a₁₂,...,a_1l)；a₂₀,表示第1个突变点到第2个突变点间的模型系数，A₂＝(a₂₁,a₂₂,...,a_2l)；a_n0,表示第n-1个突变点到第n个突变点间的模型系数，A_n＝(a_n1,a_n2,...,a_nl)。

(2)融合模型的系数需基于训练样本集和样本数据通过多元线性回归的方法获得。假设前5个循环时对应的HI为1；突变点附近的5个循环时对应的HI为k_An，k_An＝1-m_An；后5个循环时对应的HI为0。训练样本集可表示为：

$M=\left[\begin{array}{cc}{X}_T& k^1\\ .& .\\ .& .\\ .& .\\ X_{Y1}& k_{A1}\\ .& .\\ .& .\\ .& .\\ X_{\mathrm{Yn}}& k_{\mathrm{An}}\\ .& .\\ .& .\\ .& .\\ X_F& k^0\end{array}\right]---\left(9\right)$

样本集M中，X_T、X_Y1...X_Yn和X_F分别表示HI为1、k_A1...k_An和0时对应的性能参数监测值样本。以样本训练集M进行监测样本训练求解系数a_n0,例如，计算得到融合模型的系数为：

[a₁₀,a_1i]＝[335.848，-0.117，-0.003，-0.013，0.784，0.128，-0.148，0.003]

[a₂₀,a_2i]＝[0，-0.035，0.001，-0.020，38.560，0.037，-0.339，0.001]

步骤(d)中，将多阶段的融合模型系数a_n0,代入回归模型中，同时引入训练集中发动机从性能完好到失效的全寿命周期多源监测数据，可得到全寿命周期健康指标HI时间序列，如图4所示。总之，多源监测参数通过去噪处理产生的多组信号是通过多元回归的方法进行融合得到一组新的退化时间序列，即HI时间序列。

步骤(e)中，利用Kalman滤波算法，对发动机从性能完好到性能失效全过程进行最佳拟合，以确定健康指标趋势预测模型结构，并最小化预测模型的误差。

(1)通过状态空间方法建立健康指标预测模型：

x_t＝F·x_t-1+w_t   (10)

${\hat{y}}_t=H\·x_t+m_t---\left(11\right)$

式中：x_t为t时刻对应的状态向量是预测模型的时变参数；为HI预测值；F为状态转移矩阵H为观测矩阵w_t～N(0,Q_t)，m_t～N(0,R_t)，E[w_t×m_t＝0]；w_t为状态向量噪声；m_t为观测过程噪声。基于以上描述，预测模型可表示为：

${\hat{y}}_t=a_t^0+a_t^{1}t...+a_t^m{t}^m+{\mathrm{\δ}}_t---\left(12\right)$

其中，δ_t表示预测误差，m是预测模型的阶数，决定了预测模型的结构。δ_t与m的值可根据最佳拟合(预测值与实际值的方差最小)确定，最佳拟合由Kalman滤波算法实现。另外，考虑到预测模型的时变参数随时间变化的特性，采用具有实时递归运算能力的Kalman滤波算法对时变参数进行估计。

(2)利用以下Kalman滤波算法公式实现时变参数估计与最佳拟合：

$\left\{\begin{array}{c}{\hat{x}}_{t|t-1}=F\·{\hat{x}}_{t-1|t-1}\\ p_{t|t-1}=F\·p_{t-1|t-1}\·F^T+Q_t\end{array}\right.---\left(13\right)$

$\left\{\begin{array}{c}{\hat{x}}_{t|t}={\hat{x}}_{t|t-1}+K_t{\mathrm{\δ}}_t\\ p_{t|t}=p_{t|t-1}-K_t{\mathrm{Hp}}_{t|t-1}\end{array}\right.---\left(14\right)$

其中，为根据t-1时刻及该时刻之前的监测信息估计t时刻的后验状态，x_t为t时刻实际的后验状态；p_t|t-1为状态估计误差的协方差，x_t为t时刻的实际状态；K_t表示卡尔曼增益，K_t＝p_t|t-1H^T/S_t，S_t＝Hp_t|t-1H^T+R_t。v_t为预测的偏差值，y_t为实际值，为预测值。

其中时变参数估计与检测时间序列最佳拟合的算法流程：

Ⅰ.时变参数集状态集为X(t)＝{x₁,x₂,...,x_n}；

Ⅱ.对于t＝1：n，进行以下操作：

(ⅰ)建立状态方程x_t＝F·x_t-1+w_t，描述t-1时刻的状态与t时刻状态的关系；建立观测方程描述t时刻的状态与t时刻测量值的关系；其中 $x_t={[a_t^0{,a}_t^1,a_t^2,a_t^3,...,a_t^n]}^T,$ w_t～N(0,Q_t)，m_t～N(0,R_t)，E[w_t×m_t＝0],及

(ⅱ)预测阶段：根据t-1时刻的状态与t-1时刻及t-1时刻前的观测信息预测t时刻的状态：根据t-1时刻的状态协方差与t-1时刻及t-1时刻前的观测信息预测t时刻的状态协方差；p_t|t-1＝F·p_t-1|t-1·F^T+Q_t，其中

(ⅲ)更新阶段：结合卡尔曼增益更新状态与协方差的值，作为下一步预测的基础值；更新状态：更新协方差：p_t|t＝p_t|t-1-K_tHp_t|t-1；其中，K_t＝p_t|t-1H^T/S_t，S_t＝Hp_t|t-1H^T+R_t；

Ⅲ.时变参数的预测：将上一步骤中所得到的作为根据t时刻及以前的观测信息与t-1时刻的状态信息的时变参数状态预测值；引入时间因子t即可计算得到HI的预测值；

IV.最佳拟合的检验：计算不同预测模型阶数下的HI预测误差序列将最小预测误差对应的阶数作为最终预测模型的阶数，即确定中的m值，其中δ_t取误差均值。

需要注意的是，融合多源监测参数获得HI时间序列的方法是多元回归，而Kalman滤波算法一方面用于训练集下的HI最佳拟合与模型结构(δ_t与m的值)的确定，另一方面用于测试集下的HI预测模型时变参数的预测。

例如，最佳拟合的情况下，预测模型的结构为其中δ＝±0.138。此预测模型则可以用于测试集中的HI预测。

步骤(f)中，将测试集中的实时多源监测数据代入回归模型中，回归模型的系数为步骤(c)中计算得到的a_n0,融合多源监测参数，得到测试集下的健康指标HI时间序列，此时的HI是需要预测的对象。

该步骤的实质是将步骤(d)的训练数据替换为测试数据，引入多阶段融合模型中，即可得到当前发动机的实时HI时间序列。例如，计算得到训练集下的HI时间序列如图7所示。

步骤(g)中，利用Kalman滤波算法(详见步骤(e)中(2)中的Kalman算法流程)，对预测模型的时变参数进行估计。将步骤(f)中得到的HI时间序列引入HI预测模型δ＝±0.138。采用第步骤(e)中的Kalman算法流程，从而实现时变参数的实时预测。

步骤(h)中，在确定的HI预测模型中引入时间机制，从而实时估计发动机的失效时间。结合步骤(e)中确定的模型结构与步骤(g)中估计的预测模型时变参数，代入时间，预测HI的变化趋势，最终可计算发动机剩余寿命。预测的HI等于0时对应的时间即为发动机性能失效的时间。预测结果如图8所示。该发动机的实际失效时间是98个循环时之后，预测的剩余使用寿命是90，因此具有较好的预测精度。

当利用Kalman滤波算法对预测模型的时变参数进行预测时，基于以下对比可说明基于多阶段信息融合的预测方法的优越性：

(1)在融合多源监测参数计算HI时，如果不采用分阶段的融合方法，则会出现如图6虚线所示的情况。由于监测参数不分阶段处理，用于融合各个监测参数的回归模型的回归系数在整个监测时间窗内是不变的，这将导致融合得到的HI时间序列不能准确反映发动机的实际退化特性。即使Kalman滤波方法适用于非平稳数据的预测，并考虑的模型参数的时变特性，由于输入的HI不准确，得到的剩余寿命预测结果也将是不准确的。

(2)当融合多源监测参数计算HI时，采用分阶段的融合方法，则能达到图6点划线的预测效果。由于考虑的退化过程的突变特性，采用分阶段的方法对多源监测参数在不同时间窗分别进行融合，此时的回归模型的系数在整个检测时间窗是分阶段的，因此计算得到的HI可以较为准确地反映发动机的实际退化过程。此时，由于Kalman滤波算法的输入量HI较为准确，结合Kalman滤波自身的优点(非平稳数据的预测适用性以及对时变参数的预测更新能力)，预测得到的剩余使用寿命自然更接近真实值。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (10)

1.一种基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：包括如下步骤 

(a)多源监测参数去噪处理与特征提取； 

(b)在步骤(a)的基础上，对多源监测时间序列进行平稳性分析，计算各个参数监测时间序列突变点，计算突变点处参数退化比例； 

(c)在步骤(b)的基础上，对多源参数进行多阶段划分，建立回归融合模型，利用历史监测数据进行样本训练，获得融合模型在多阶段内的参数； 

(d)在步骤(c)的基础上，结合训练集中的监测数据，融合多源监测参数，得到健康指标HI； 

(e)在步骤(d)的基础上，利用Kalman滤波算法，对发动机从性能完好到性能失效全过程进行最佳拟合，以确定健康指标趋势预测模型结构，并最小化预测模型的误差； 

(f)在步骤(c)的基础上，结合测试集中的实时监测数据，融合多源监测参数，得到健康指标HI，此时的HI是预测对象； 

(g)在步骤(f)的基础上，利用Kalman滤波算法，对预测模型的时变参数进行实时估计； 

(h)在步骤(e)和步骤(g)的基础上，确定预测模型，引入时间机制，从而实时估计发动机的失效时间。 

2.如权利要求1所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：步骤(a)中，多源监测参数去噪处理与特征提取，采用经验模态分解方法EMD对多源监测数据进行去噪处理，EDM将原始监测数据分解为从高频到低频的多个本征模函数IMF与残差分量之和，IMF包括信号起主导作用的模态与噪声起主导作用的模态，从IMF中剔除噪声起主导作用的模态，将信号起主导作用的模态与残余项进行叠加，即可进行部分重构得到降噪之后的监测时间序列，EMD分解的公式如下： 

其中，z(t)为监测时间序列，IMF₁,IMF₂,...,IMF_n分别代表从高到低不同频率段的信号 成分，残余项r_n代表信号的平均趋势。 

3.如权利要求2所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：步骤(b)中包括如下： 

(1)采用以下启发式算法对时间序列的平稳性进行分析并划分时间窗： 

其中，N₁、N₂分别表示i点左右两部分的时间点个数；u₁(i)、u₂(i)分别表示i点左右两部分的均值；S₁(i)、S₂(i)分别表示i点左右两部分的标准偏差；S_D(i)为合并偏差；T(i)为检验统计值； 

(2)将划分时间窗的点作为时间序列的分割点，则需判断分割点是否为性能突变点，计算慢变量： 

其中，S_V为慢变量，Δt为时间间隔，j＞i，i,j＝1,2,...,n，n为时间点的个数，ΔS为监测值在相应时间间隔内的变化量，通过时间窗口Δt的滑动求得整个时间窗内的慢变量序列，慢变量的突变即体现参数退化的突变点； 

(3)各个监测参数时间序列的突变点确定后，即计算监测参数平均突变点与平均变化比例： 

其中，P_An为第n个平均突变点；m_An为第n个平均突变点处的平均变化比例；P_in为第i个个监测参数的第n个突变点，m_in为第i个监测参数从退化初到第n个突变点的变化 量在整个退化过程中的变化量的比例，i＝1,2,...,l，l为监测参数的个数，n为突变点的个数。 

4.如权利要求3所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：步骤(c)中包括如下： 

(1)基于步骤(b)中获得的突变点，建立分阶段融合模型： 

其中，y1表示第1个突变点前的HI；y2表示第1个突变点到第2个突变点间的HI；yn表示第n-1个突变点到第n个突变点间的HI；X表示多源监测参数(x₁,x₂,...,x_l)，l表示性能参数的个数；a₁₀,表示第1个突变点前的模型系数A₁＝(a₁₁,a₁₂,...,a_1l)；a₂₀,表示第1个突变点到第2个突变点间的模型系数，A₂＝(a₂₁,a₂₂,...,a_2l)；a_n0,表示第n-1个突变点到第n个突变点间的模型系数，A_n＝(a_n1,a_n2,...,a_nl)； 

(2)模型系数通过训练样本集和样本数据通过多元线性回归的方法获得，假设前5个循环时对应的HI为1；突变点附近的5个循环时对应的HI为k_An，k_An＝1-m_An；后5个循环时对应的HI为0，训练样本集可表示为： 

样本集M中，X_T、X_Y1...X_Yn和X_F分别表示HI为1、k_A1...k_An和0时对应的性能参数监测值样本，以样本训练集M进行监测样本训练求解系数a_n0,。

5.如权利要求4所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：步骤(d)中，将多阶段的融合模型系数a_n0,代入回归模型中，同时引入训练 集中发动机从性能完好到失效的全寿命周期多源监测数据，得到全寿命周期健康指标HI时间序列。 

6.如权利要求5所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：步骤(e)中包括如下： 

(1)通过状态空间方法建立健康指标预测模型： 

x_t＝F·x_t-1+w_t

式中：x_t为t时刻对应的状态向量是预测模型的时变参数；为HI预测值；F为状态转移矩阵H为观测矩阵H＝[t⁰,t¹,t²,t³,t⁴,...,tⁿ]；w_t～N(0,Q_t)，m_t～N(0,R_t)，E[w_t×m_t＝0]；w_t为状态向量噪声；m_t为观测过程噪声，基于以上描述，预测模型表示为： 

其中，δ_t表示预测误差，m是预测模型的阶数，决定了预测模型的结构，δ_t与m的值可根据最佳拟合即预测值与实际值的方差最小确定，最佳拟合由Kalman滤波算法实现； 

(2)利用以下Kalman滤波算法公式实现时变参数估计与最佳拟合： 

其中，为根据t-1时刻及该时刻之前的监测信息估计t时刻的后验状态， x_t为t时刻实际的后验状态；p_t|t-1为状态估计误差的协方差， x_t为t时刻的实际状态；K_t表示卡尔曼增益，K_t＝p_t|t-1H^T/S_t， S_t＝Hp_t|t-1H^T+R_t，v_t为预测的偏差值，y_t为实际值，为预测值。 

7.如权利要求6所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：所述时变参数估计与时间序列最佳拟合检测的算法流程包括： 

Ⅰ.时变参数集状态集为X(t)＝{x₁,x₂,...,x_n}； 

Ⅱ.对于t＝1：n，进行以下操作： 

(ⅰ)建立状态方程x_t＝F·x_t-1+w_t，描述t-1时刻的状态与t时刻状态的关系；建立观测方程描述t时刻的状态与t时刻测量值的关系；其中 H＝[t⁰,t¹,t²,t³,t⁴,...,tⁿ],w_t～N(0,Q_t)，m_t～N(0,R_t)，E[w_t×m_t＝0],及

(ⅱ)预测阶段：根据t-1时刻的状态与t-1时刻及t-1时刻前的观测信息预测t时刻的状态：根据t-1时刻的状态协方差与t-1时刻及t-1时刻前的观测信息预测t时刻的状态协方差；p_t|t-1＝F·p_t-1|t-1·F^T+Q_t，其中

(ⅲ)更新阶段：结合卡尔曼增益更新状态与协方差的值，作为下一步预测的基础值；更新状态：更新协方差：p_t|t＝p_t|t-1-K_tHp_t|t-1；其中，K_t＝p_t|t-1H^T/S_t，S_t＝Hp_t|t-1H^T+R_t； 

Ⅲ.时变参数的预测：将上一步骤中所得到的作为根据t时刻及以前的观测信息与t-1时刻的状态信息的时变参数状态预测值；引入时间因子t即可计算得到HI的预测值； 

IV.最佳拟合的检验：计算不同预测模型阶数下的HI预测误差序列将最小预测误差对应的阶数作为最终预测模型的阶数，即确定中的m值，其中δ_t取误差均值。 

8.如权利要求7所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特 征在于：所述步骤(f)中，将测试集中的实时多源监测数据代入回归模型中，回归模型的系数为步骤(c)中计算得到的融合多源监测参数，得到测试集下的健康指标HI时间序列，此时的HI是需要预测的对象。 

9.如权利要求8所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：所述步骤(g)中,利用Kalman滤波算法对预测模型的时变参数进行估计,将步骤(f)中得到的HI时间序列引入HI预测模型，采用步骤(e)中的Kalman算法流程实现时变参数的实时预测。 

10.如权利要求9所述的基于多阶段信息融合的航空发动机剩余寿命预测方法，其特征在于：所述步骤(h)中结合步骤(e)中确定的模型结构与步骤(g)中的预测模型时变参数，代入时间，预测HI的变化趋势，最终计算出发动机剩余寿命。