CN106021826A - 一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法 - Google Patents

Abstract

本发明提出了一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法，该方法着力解决复杂的工况变化掩盖航空发动机真实性能退化趋势致使整机剩余寿命难于预测这个问题，包括如下步骤：首先，识别航空发动机的运行工况；其次，对不同工况下的历史退化数据进行数据标准化；再次，对标准化后的数据进行传感器选择和参数维度约减；接着，使用相似性匹配方法匹配每一台参考发动机和待预测发动机的退化轨迹，获得待预测发动机的剩余寿命估计值和两者之间的相似性距离；最后，通过相似性距离生成加权权重，加权这些寿命估计以获得待预测发动机的剩余寿命。经验证，本发明对不同的发动机测试样本均具有较高的寿命预测精度。

Description

一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法

技术领域

本发明涉及航空发动机剩余寿命预测的技术领域，具体涉及一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法。

背景技术

目前我国航空业正处于蓬勃发展阶段，但依然存在维护和维修成本高等问题。作为飞机核心动力源的航空发动机，是一种安全和可靠性要求极高，设计生产和维修保养花费极大的大型复杂设备。资料显示，自1963年至1975年期间，美国空军的战斗机共发生了飞行事故3824起，其中由于发动机原因导致的事故有1664起，占到总数的43.5％；自1994年至2000年世界航空运输机事故共279，其中因发动机的故障导致飞行事故约20％，因此作为飞机心脏的发动机的维修是航空安全保障的重点。世界范围内航空公司每年大致花费31亿美元用于飞机维修，其中用于航空发动机的维修费用约占总费用的30％，因此需要不断发展与改进现有的航空发动机健康管理技术和维修策略，以减少飞行事故的发生并削减维修花费。随着传感器及监测技术的迅速发展，在航空发动机的使用过程中，积累了大量与寿命信息相关的性能监测数据，采用这类性能退化数据对高可靠和长寿命设备进行寿命预测成为发展趋势。

当前，国内对航空发动机寿命预测的研究还停留在单个零部件或滑油系统上，针对航空发动机整机剩余寿命预测的研究基本没有，并且现有研究都是基于单工况假设下进行的。然而，航空发动机常工作在复杂多变的工况环境下，变化的工况掩盖了系统真实的性能退化趋势，在此情况下现存的技术方法不能进行准确的剩余寿命预计。

发明内容

本发明的目的在于：针对变工况下航空发动机整机剩余寿命难于预测的问题，提出一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法。

本发明采用的技术方案为：一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法，其特征在于：

(1)基于K均值聚类的航空发动机工况识别方法，识别系统在不同时刻所处的工况后将不同工况下的数据标准化到统一的参数空间；

(2)为最大限度的挖掘和借鉴参考发动机退化轨迹的衰退规律，通过时间定位和时间尺度伸缩最优地匹配待预测发动机与参考发动机的性能退化轨迹；

(3)将参考发动机伸缩后轨迹的剩余寿命作为待预测发动机的剩余寿命估计，并使用相似性距离加权这些估计获得待预测发动机的剩余寿命。

本发明与现有技术相比的优点在于：

(1)、本发明通过工况识别方法解决了航空发动机寿命预测领域目前尚未很好解决、而实际中常常存在的变工况问题；

(2)、本发明提出时间定位和时间尺度伸缩来最大化地匹配参考发动机轨迹与待预测发动机轨迹，提高了历史全寿命数据在预测中的参考价值，算法具有较高的预测精度；

(3)、本发明借鉴集成学习中的某些思想，对待预测发动机寿命进行多次预测后加权集成获得待预测的寿命，算法稳健性好，抗干扰能力强。

附图说明

图1为变工况下航空发动机整机剩余寿命预测流程图；

图2为K均值聚类流程图；

图3为工况识别流程图；

图4为变工况工作示意图；

图5为多工况下的退化特性融合示意图，其中，图5(a)为坐标平移示意图，图5(b)为尺度伸缩示意图；

图6为多工况数据标准化流程图；

图7为时间定位示意图；

图8为伸缩匹配原理图，其中，图8(a)为时间尺度伸缩原理图，图8(b)为最优匹配示意图；

图9为涡轮风扇发动机气路部件简图；

图10为工况聚类图；

图11为传感器1，2和3标准化前后的退化数据对比(训练集1#发动机)，其中，图11(a)为标准化前的退化数据，图11(b)为标准化之后的退化数据。

图12为维度约减后的性能退化轨迹(训练集1#发动机)；

图13为测试集1#发动机的3个寿命估计值，其中，图13(a)为使用训练集1#发动机的估计结果，图13(b)为使用训练集2#发动机的估计结果，图13(c)为使用训练集3#发动机的估计结果。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

如图1所示，基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法，具体步骤如下：

1工况识别

复杂的工况变化常常掩盖航空发动机的性能退化，为还原真实的退化趋势，采用先辨识各状态的工况，然后进行数据标准化的策略。本发明采用K均值聚类算法对工况参数进行聚类和分类，从而辨识出系统所处各个状态的运行工况信息。在数据挖掘中，K均值聚类算法应用广泛，是经典聚类算法之一。

设X＝{X₁,X₂,...,X_n}为已知数据集，X中的X₁,X₂,...,X_n是n个数据对象并且每个数据对象都是N维的，即X_i＝(x_i1,x_i2,...,x_iN)。K均值聚类算法就是要找到含有K个聚类中心的集合C＝{C₁,C₂,...,C_K}＝{(c₁₁,c₁₂,...,c_1N),(c₂₁,c₂₂,...,c_2N),...,(c_K1,c_K2,...,c_KN)}使得目标函数：

$J(X,C)=\min \left(\underset{i=1}{\overset{K}{\Σ}}\underset{j=1}{\overset{n_i}{\Σ}}d\right(C_i,X_j\left)\right)---\left(1\right)$

其中，n_i被归为类C_i的数据对象点数，d(C_i,X_j)表示聚类中心与数据对象的欧几里德距离，其定义如下：

$d(C_i,X_j)=\sqrt{{(c_{i1}-x_{j1})}^2+{(c_{i2}-x_{j2})}^2+\mathrm{...}+{(c_{iN}-x_{jN})}^2}---\left(2\right)$

K均值聚类算法的核心思想是把数据集划分成使目标函数达到最小值的K个类。具体步骤如图2所示。

针对训练集数据通过聚类得到各个工况的聚类中心和半径。对于测试数据首先计算实时工况参数与各个聚类中心的距离，距离较近者即为该工况。工况识别的流程如图3所示。

2数据标准化

复杂系统运行在变工况下时，由于测量尺度的不一致性，状态特征监测数据并不能真实反映出系统的退化趋势。如图4所示，假设随着系统的退化，系统的某个特征参数是呈下降趋势的，3条实线(红、黄、蓝)展示了系统在3种运行工况下不同退化的曲线。情况1曲线(一种包含多个工况的运行模式)是上下波动的，完全不能反映出系统随着时间的真实退化趋势；情况2曲线在一定程度上反映了系统的退化趋势，但是是虚假的、不真实的；情况3曲线表象地完全呈现了与真实退化完全相反的趋势。

为实现寿命预测的目的，需要采用数据处理方法，将工况造成的系统监测特征数据的不同尺度特征去掉，使不同工况特征之间具有可比性、呈现一致的退化特征。直观的，就是将图4中所示的3种不同工况下的退化曲线合并(融合映射)到一起，如图5所示。至此，不同运行工况下的特征数据被转化到了同一个尺度之下，具有切实的可比性，并能反映出系统的退化趋势。从图形上看，将三种工况合并的两个核心步骤是：首先，将三条曲线的重心移到一起(与X轴重合)；其次，对坐标平移后的数据进行尺度缩放，归一化到一个可比的特征空间。

多工况数据标准化方法是使用特征层预测方法预测剩余寿命时的预处理环节。使用该方法，不同工况模式下的数据的N个特征将被分别标准化到一个标定的范围，然后形成一个与原始数据特征维数相同的新的N维时间序列。图6展示了多工况数据标准化的流程。

本发明采用正态标准化(Z-Score)方法将多工况特征数据标准化。令在时刻t_i的工况参数为u_i，性能状态参数x_i，合在一起构成系统的状态参数集(t_i,u_i,x_i),i＝1,2,..,n。每个系统状态都将被聚类到P个工况模式O_p之一中，因而形成性能状态样本子集{x}^(p)＝{x_i|u_i∈O_p},i＝1,2,..,n_p,p＝1,2,..,P。每个工况模式的样本均值和样本标准差为：

$\begin{array}{c}{\stackrel{\‾}{x}}^{\left(p\right)}=Mean\left({\left\{x\right\}}^{\left(p\right)}\right)\\ s^{\left(p\right)}=Std\left({\left\{x\right\}}^{\left(p\right)}\right)\end{array}---\left(3\right)$

对属于工况p的n_p个状态参数向量进行标准化：

$y_{ij}=(x_{ij}-{\stackrel{\‾}{x}}_j^{\left(p\right)})/s_j^{\left(p\right)},i=1,2,\mathrm{...},n_p,j=1,2,\mathrm{...},N---\left(4\right)$

式中，n_p为属于工况p的数据状态数，N为性能参数维度。

该标准化方法假定在任一使用循环中，每一个工况发生的概率相等，以确保在同一运行工况下计算得到的样本均值和标准差是给定工况模式下特征均值和标准差的无偏估计。最后，合并不同工况模式下的标准化状态参数得到系统退化N维状态时间序列Y＝{y₁,y₂,...,y_n}。

3状态参数维度约减

高维状态监测参数通常很难处理或处理时计算量较大，实际中通常需要对其进行数据维度的约减。主成分分析(PCA,Principle Component Analysis)是应用最广泛的一种降维方法之一，是一种统计学方法，在信号处理、模式识别、数字图像处理等领域已经得到了广泛的应用。主成分分析方法基本思想是提取出空间原始数据中的主要特征(主元)，减少数据冗余，使得数据在一个低维的特征空间被处理，同时保持原始数据的绝大部分的信息，从而解决数据空间维数过高的瓶颈问题。

设航空发动机退化状态矩阵为X，其列向量X_k＝(x_1k,x_2k,...,x_nk)^T为n维状态参数，发动机的某一性能状态可由x_k描述，x_k的协方差矩阵为：

$R_x=\underset{k=1}{\overset{N}{\Σ}}(X_k-\stackrel{\‾}{x}){(X_k-\stackrel{\‾}{x})}^T---\left(5\right)$

式中，N为退化状态的采样数，为各状态参数的均值：

$\stackrel{\‾}{x}=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}{x}_k---\left(6\right)$

求解R_x的全部特征值λ_i(i＝1,2,...,n)和特征向量v_i，将特征值λ_i按照从大到小的顺序排列：λ₁>λ₂>...>λ_n，则对应的特征向量为v_i(i＝1,2,..,n)。样本x_i投影到特征向量v_i得到该方向对应的主分量为：

$y^{ij}={v_i}^T(x_j-\stackrel{\‾}{x})---\left(7\right)$

所有的特征向量张成一个n维正交空间，x投影到该正交空间得到相应的n维主分量。特征向量所对应的特征值越大，它在重构时的贡献也越大，特征值越小的特征向量在重构时的贡献就越小。设正交空间中前m个主分量为y¹,y²,...,y^m，其累计方差贡献率为：

$h\left(m\right)=\underset{i=1}{\overset{m}{\Σ}}{\mathrm{\λ}}_i/\underset{k=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_k---\left(8\right)$

当前面少数几个主分量的累计方差贡献率足够大，如h(m)>95％，即95％以上的原始数据信息保留在前面几个主分量中，可只取前m(m<n)个主分量来表征原始信息,在保证信息完整的情况下，达到数据降维的目的。

4相似性匹配

(1)时间定位

如图7所示，由于不知道待预测发动机轨迹起点的位置，首先需要进行待预测轨迹的时间定位，找到参考轨迹中与待预测轨迹初始状态参数最接近的时间点。然后将两个轨迹从该时间点截断，取截断点后的两个轨迹进行后续的进一步匹配计算。设L台参考发动机的全寿命循环数分别为{ⁱC},i＝1,2,...,L，待预测发动机的采样点数为n₀，将待预测发动轨迹与各参考发动机进行时间定位得到一个集合{ⁱt₀},i＝1,2,...,L。

(2)时间尺度伸缩匹配

如图8(a)所示，令离散函数x(t),t＝0,1,...,m和y(t),t＝0,1,...,n分别表示参考发动机和待预测发动机的退化轨迹。两个轨迹需要满足以下条件：(A)相同的起点；(B)自变量为正整数。

首先，用一个因子λ去伸展或压缩x(t)的时间尺度，得到一个临时轨迹x'(t)＝x(λt),t＝0,1/λ,...,m/λ；λ∈R⁺。由于临时轨迹的采样点不符合整数采样要求，使用PCHIP插值算法插值序列x'(t)从而获得一个新轨迹x"(t),t＝0,1,...,[m/λ]。

PCHIP算法具体为：令a＝x₀<x₁<...<x_n＝b为区间I＝[a,b]的一个划分，f_k,k＝0,1,..,n为划分点的函数值，令h_k＝x_k+1-x_k,p_k＝f'(x_k)。PCHIP函数I_h(x)满足如下限制：

A.定义域为I＝[a,b]；

B.在每个子区间I_k＝[x_k,x_k+1](k＝0,1,...,n-1)，I_h(x)可以用一个三次样条多项式等式I_k(x)表示；

C.插值条件：I_h(x_k)＝f_k,I_h'(x_k)＝p_k；(k＝0,1,2,...,n)。

在每个子区间[x_k,x_k+1]，PCHIP函数可表示为：

I_k(x)＝f_k+c_k,1(x-x_k)+c_k,2(x-x_k)²+c_k,3(x-x_k)³ (9)

式中：

$\left\{\begin{array}{c}{c}_{k,1}=p_k\\ c_{k,2}=\left(\frac{3}{h_k}\right(f_{k+1}-f_k)-2p_k-p_{k-1})\frac{1}{h_k}\\ c_{k,3}=(p_{k+1}+p_k-2\frac{f_{k+1}-f_k}{h_k})\frac{1}{h_k^2}\\ h_k=x_{k+1}-x_k,p_k=f^{\&prime;}\left(x_k\right)\\ x\&Element;\&lsqb;x_k,x_{k+1}\&rsqb;,k=0,1,...,n-1\end{array}\right.---\left(10\right)$

接着，定义一个相似性距离去量化x"(t)和y(t)之间的相似性：

$D_{\mathrm{\&lambda;}}\left(x\right(t),y(t\left)\right)=\underset{t=0}{\overset{Q}{\&Sigma;}}\left(\frac{{\mathrm{\&alpha;}}^{Q-t}}{\underset{t=0}{\overset{Q}{\&Sigma;}}{\mathrm{\&alpha;}}^{Q-t}}\right|x^{\&prime;\&prime;}\left(t\right)-y\left(t\right){|}^2),Q=min(\&lsqb;m/\mathrm{\&lambda;}\&rsqb;,n)---\left(11\right)$

式中，α是一个0-1之间的平滑参数。在本发明中，α＝0.98。因为近期时间点对未来退化具有更为重要的影响，因而通过给距当期时刻较近的时间点分配更大的权重。[·]表示取整算子。

如图8(b)所示，随着伸缩因子λ的变化，相似性距离随之变化。显然这是一个凸优化问题，也就是说存在一个最优的λ_opt是相似性距离D_λ(x(t),y(t))最小。最优函数值被定义为最优匹配距离D_OM。

$D_{OM}=\underset{\mathrm{\&lambda;}}{\min }{D}_{\mathrm{\&lambda;}}\left(x\right(t),y(t\left)\right)---\left(12\right)$

式中，D_λ(x(t),y(t))为式(11)定义的相似性距离。

5集成加权

L台参考发动机在时间定位点之后的寿命分别记为{^lM},l＝1,2,...,L，其为：

^lM＝^lC-^lt₀,l＝1,2,...,L (13)

式中，^lC为第l个参考发动机的全寿命，^lt₀为待预测发动机与第l个参考发动机的时间定位位置。

通过相似性匹配，得到L个最优伸缩因子和最优匹配距离{^lλ_opt,^lD_OM},l＝1,2,...,L，通过最优伸缩因子计算L个待预测发动机的剩余寿命估计值：

^lR＝^lλ_opt·^lM-n₀,l＝1,2,...,L (14)

式中，n₀为待预测发动机的数据点数，^lλ_opt为最优匹配因子。

至此，得到L个相似性距离和对应的L个RUL估计值，用集合{^lD_OM,^lR},l＝1,2,...,L来表示。为获取待预测发动机的剩余寿命，通过相似性距离生成相似性权重。

$W_l=\frac{S_l}{\underset{l=1}{\overset{L}{\&Sigma;}}S_l},l=1,2,\mathrm{...},L---\left(15\right)$

式中，^lS是第l个参考发动机的相似性的分，其为

加权集成的目的是集成多个寿命估计最终获得待预测发动机的寿命。最简单的集成方法是加权和，待预测发动机RUL的点估计为：

$r=\underset{l=1}{\overset{L}{\&Sigma;}}{R}_l\&CenterDot;W_l---\left(16\right)$

实验验证

本发明使用美国航空航天局NASA提供的涡扇发动机数据集的Dataset 2(包括测试集、训练集)对提出的预测方法进行有效性验证。该数据集由C-MAPSS(Commercial ModularAero-propulsion System Simulation)大型涡轮风扇发动机仿真模型得到。模型通过14个参数的输入来模拟包括风扇、低压压气机、高压压气机、高压涡轮、低压涡轮(Fan,LPC,HPC,HPT,LPT)在内的五大旋转部件失效影响和退化过程，发动机仿真模型主要部件简图如图9所示。

数据由发动机不同运行循环时间点的3个工况参数(飞行高度、飞行马赫数和油门杆角度)21个传感器监测的性能参数组成。传感器监测值是针对研究模块，通过发动机热动力仿真模型获取，且包含噪声。训练数据包括多组独立同一单元的多元时间序列样本，反映了每个单元从初始到失效的全寿命周期的变化过程。每个部件的初始时刻的退化情况是随机且未知的，是否发生退化以及退化程度不同，随着发动机的运行，当性能降级到一个阈值的后整个系统发生失效。用于衡量发动机性能状态的21个监测参数值如表1所示。

表1模型输入参数

序号	描述	符号
			1	燃料流量	Wf(pps)
2	风扇效率参数	fan_eff_mod
			3	风扇流量参数	fan_flow_mod
4	风扇压比参数	fan_PR_mod
			5	低压压气机效率参数	LPC_eff_mod
6	低压压气机流量参数	LPC_flow_mod
			7	低压压气机压比参数	LPC_PR_mod
8	高压压气机效率参数	HPC_eff_mod
			9	高压压气机流量参数	HPC_flow_mod
10	高压压气机压比参数	HPC_PR_mod
			11	高压涡轮效率参数	HPT_eff_mod
12	高压涡轮流量参数	HPT_flow_mod
			13	低压涡轮效率参数	LPT_eff_mod
14	低压涡轮流量参数	LPT_flow_mod

表2监测参数表征量

1)工况识别

使用训练集数据通过K均值聚类得到6个工况的聚类中心和半径分别如表3所示，工况聚类图如图10所示。对于测试数据首先计算实时工况参数与6个聚类中心的距离，距离最近者即判为该工况。

表3聚类中心和半径

2)变工况数据标准化

即使是同一种状态监测参数，不同工况下的数据不具有可比性，使用本发明提出的方法对各工况数据进行标准化处理，去除工况的干扰信息，训练集1#发动机的传感器1，2和3标准化前后的数据如图11所示。

从图11可以看出，传感器2和3的监测参数在标准化后呈现出了一定的趋势性，而传感器1仍然没有明显的退化趋势，说明传感器1完全不能表征任何系统的退化信息。因此，接下来要进行参数的的选择或降维，去除没有趋势性的传感器参数，并对有趋势的参数进行维度的约减。

3)多元参数维度约减

训练集1#发动机状态参数维度约减后的轨迹如图12所示。

4)相似性匹配计算

用训练集1#、2#、3#发动机对测试集1#发动机进行剩余寿命估计的估计结果如图13所示，其中测试机1#发动机的真实剩余寿命为18个循环。

5)预测结果分析

为验证所提出方法更多发动机样本的适用性，使用本发明所提出的方法对测试样本集中随机选择的10台变工况下运行的航空发动机进行算法的进一步验证，预测结果如表4所示。

表4预测结果

序号	真实的剩余寿命	预测的剩余寿命	绝对误差	相对误差(％)
					1	18	15	-3	16.67％
2	79	68	-11	13.92％
					3	79	70	-9	11.39％
4	73	57	-16	21.92％
					5	30	32	2	6.67％
6	67	59	-8	11.94％
					7	99	76	-23	23.23％
8	22	17	-5	22.73％
					9	54	42	-12	22.22％
10	37	31	-6	16.22％
					均值	--	--	-9.1	16.69％

从表4可知，本发明提出的航空发动机整机剩余寿命预测方法的平均预测百分比误差为17％左右，对不同的发动机样本均具有较高的预测精度，提出的预测方法可为航空发动机的整机剩余寿命预测提供一个新的研究思路。

Claims (1)

1.一种基于工况识别和相似性匹配的变工况下航空发动机整机剩余寿命预测方法，其特征在于：

(1)基于K均值聚类的航空发动机工况识别方法，识别系统在不同时刻所处的工况后将不同工况下的数据标准化到统一的参数空间；

(2)为最大限度的挖掘和借鉴参考发动机退化轨迹的衰退规律，通过时间定位和时间尺度伸缩最优地匹配待预测发动机与参考发动机的性能退化轨迹；

(3)将参考发动机伸缩后轨迹的剩余寿命作为待预测发动机的剩余寿命估计，并使用相似性距离加权这些估计获得待预测发动机的剩余寿命。