CN108563806A - 基于相似性的发动机气路参数长期预测方法及系统 - Google Patents

Abstract

本发明涉及一种基于相似性的发动机气路参数长期预测方法及系统，其中方法包括：通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离利用得到的统计距离和历史轨迹样本，针对每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合；通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。本发明相对于自回归滑动平均、反向传播神经网络和传统的基于相似性的预测方法相比具有更高的预测精度。

Description

基于相似性的发动机气路参数长期预测方法及系统

技术领域

本发明涉及航空发动机性能预测技术领域，尤其涉及一种基于相似性的发动机气路参数长期预测方法及系统。

背景技术

作为一类高价值、高可靠性复杂装备，航空发动机的工作性能随着新的技术的应用而得到不断提升，相应的，其结构复杂度、采购价格和维护价格也水涨船高。为了适应新时代的航空业发展需求，航空发动机又面临着来自经济性和安全性方面的压力。与此同时，当今的能效管理与环境保护标准也对发动机的性能水平和可靠性水平也提出了更高的要求。为了应对维护难度、维修成本、运营性能等因素给发动机维修方面所带来的多重压力，学者们提出了视情维修(Condition Based Maintenance,CBM)和预测与健康管理(Prognostic and Health Management,PHM)的理念。这些理念旨在利用机内测试传感器所提供的数据特征来进行异常探测、故障诊断和性能衰退趋势预测，以此实现预防性维修，降低预计外故障的发生概率。顾名思义，对发动机衰退特征趋势的预测是PHM中的一项重要内容，这项技术对预防性维修的实施具有决定性意义，是PHM领域的研究热点之一。

一般情况下，发动机的性能状态特征具有多元时间序列的形式。现今，能够适用于多元参数时间序列预测的理论框架有很多，例如滑动自组织平均模型(Auto-RegressiveMoving Average,ARMA)、灰色模型(Grey Model,GM)、蒙特卡洛类方法、模糊逻辑、人工神经网络(Artificial Neural Network,ANN)、卡尔曼滤波(Kalman Filter,KF)类方法。这些方法有的具有确定形式的输出，有的具有概率形式的输出，在实际的性能状态短期预测实验当中都体现出了较高的精度和鲁棒性。然而，不同于短期预测领域的繁荣，针对发动机乃至于广义的复杂设备的性能状态的长期预测的研究和理论却极为有限。虽然具有不同的侧重点，和短期预测一样，长期预测对于设备的PHM同样具有重要意义：一方面，因为几乎所有的维护维修计划都需要在得知设备尽可能远的预期性能的前提下进行制定，长期预测给出的设备性能长期变化趋势对维修范围制定环节提供了不可或缺的信息支撑；另一方面，因为发动机的长期的性能衰退受到更多的随机事件的影响(例如小范围的维修调整，以及微小结构损伤)，算法输出的预测结果应该具有概率形式，以实现预测结果的置信程度的准确表达。然而，上述所提到的方法的长期预测性能远远不能满足实用要求，而且大多数的方法只能够给出确定性的预测结果，难以支撑后续的管理决策环节。即使部分算法，如粒子滤波和卡尔曼滤波器能够给出性能特征的长期分布，但是这些方法仍然不能提供性能特征在未来各个时间点的概率密度函数的准确估计。因此，作为一种非线性、动态的、高度随机的过程，性能参数的长期衰退趋势预测一直被各种预测类研究当作一件在讨论范围之外的议题而被避而不谈。

着眼于这个问题，有研究学者提出了基于相似性的预测(Similarity BasedPrediction,SBP)理论，进而发展出了一系列的适用于长期性能预测的方法。SBP提出了一种假设，这种假设规定若待预测衰退样本的历史特征轨迹与某些历史样本的对应时间段的轨迹相似，那么其以后的发展趋势一定与这部分历史样本的后续发展轨迹相似。基于这个假设，SBP首先评估待预测轨迹与历史数据库中各样本轨迹的相似性，而后利用得到的相似程度和历史轨迹的后续发展轨迹实现待预测轨迹在未来各个时间点的概率密度函数的重构。这个技术框架使SBP能够有效地利用历史样本提供的同构信息进行目标样本未来发展趋势的刻画，而且，只要目标样本的待预测时间点能够被历史样本的长度所覆盖，SBP便能够对其进行预测。

观察SBP的技术框架可以看出，SBP的预测性能取决于其在相似性评估步骤所采用的度量方法和概率密度函数重构步骤所采用的重构方法。

SBP首次被提出时，见文献1“J.Liu,D.Djurdjanovic,J.Ni,N.Casoetto,J.Lee,Similarity based method for manufacturing process performance prediction anddiagnosis,Computers in industry,58(2007)558-566”，衰退轨迹之间的相似程度是用维度的匹配矩阵来度量的，匹配矩阵所包含的元素是两个衰退序列之间的逐点交叉计算得到的Mahalanobis距离。其中T_objective表示目标轨迹的长度，表示第i个历史轨迹的长度。匹配矩阵法有效地解决了计算两个不同长度的时间序列之间的相似程度的度量问题，但是如果衰退轨迹的长度达到数千或者更长，那么其计算负担将会达到不可接受的程度。文献2“A.Bleakie,D.Djurdjanovic,Analytical approach to similarity-based prediction of manufacturing system performance,Computers in Industry,64(2013)625-633”选择对等长的多元时间序列进行相似性评估，并对逐点计算得到的Mahalanobis距离进行平均，凭此提高了计算速度。文献3“T.Wang,Trajectory similaritybased prediction for remaining useful life estimation,University ofCincinnati,2010”进一步改进了相似性评估的方法，其加入了时间迟滞和时间放缩因子的考量，比较适用于具有不同衰退速率和衰退起始点的轨迹之间的相似性的对比，但是这种方法的迟滞因子和放缩因子的搜索速度较慢，不适用于具有较多预测参数和预测样本点的预测过程。上述文献2和文献3同时还加入了时间久远程度的考量，在距离加权平均的过程中降低了距当前时刻较远的逐点Mahalanobis距离的权重，以加强距离较近的数据对模型的影响。所有传统的SBP方法对应于不同的状态参数(即归一化之前的特征元素)都应用了统一的协方差矩阵，对所获得的逐点距离都进行了加权平均，且都选择了通过s_oj＝exp(-βd_oj)形式的指数映射实现距离到相似度的转化。在实际应用场合中，这些运算过程将带来有不可忽视的隐患：统一的协方差矩阵对于样本轨迹噪声大小的差异性的区分能力不足，对逐点距离的加权平均运算消除了两条轨迹之间的不相似程度的统计性信息，而相似度转换过程引入了主观因素，容易造成对原始信息的扭曲。因此，上述所提到的相似程度度量方法极易造成对轨迹相似程度的误判。

对于概率密度函数估计环节，文献1和文献2采用了GMM、加权似然估计方法(Weighted Likelihood Estimation,WLE)、期望最大算法(Expectation-Maximization,EM)算法，而文献3采用了核密度估计(Kernel Density Estimation,KDE)方法。WLE能够使目标轨迹的预期发展趋势向具有更高相似度的历史轨迹的方向偏离，从而有效地提高了预测结果的可靠性。但是值得注意的是，WLE和KDE并没有考虑到离群样本的预测情况，而当随着性能状态特征维度的升高，离群样本在整个样本群体中所占的比例会显著提升。对于离群样本而言，其相对于历史样本群体所表现出的普遍不相似性，使得概率密度估计过程难以找到可信的参考样本，最终导致预测结果失去意义。与此同时，WLE和KDE给预测方法的泛化能力带来了较大的损害，而泛化能力方面的讨论被之前的SBP方法研究所普遍忽视，如图1所示。如图1所示，在不失一般性的前提下，若将历史衰退轨迹简化到极端情形，可以认为所有的历史衰退轨迹都被约束在典型衰退模式所驱动的轨迹束的内，如衰退模式A和衰退模式B。显然，当目标轨迹处于轨迹束A和轨迹束B内部时，利用WLE和KDE方法能够有效排除不相似历史轨迹的干扰，从而自然而然地获得准确的预测结果。然而，对于处于边缘地带的目标轨迹，例如图1中的红色轨迹，根据该条轨迹与轨迹束A和轨迹束B所包含的样本的相似性信息，WLE和KDE的内在的加权平均机制会不变地认为该目标轨迹处于轨迹束A和轨迹束B之间(如图中的红色虚线所示)，而接下来，其未来的发展趋势也会相应地在轨迹束A和轨迹束B之间得到重构。换句话来说，由加权平均方法的固有特性所决定，WLE和KDE所生成的轨迹一定会不可避免地处于被作为基使用的典型历史衰退轨迹在高维特征空间所形成的包络面之内，这种特性势必会造成预测方法具有难以消除的系统误差，从而给模型的泛化能力带来较大的损害。

发明内容

本发明要解决的技术问题在于，针对现有长期的发动机预测方法的上述缺陷，提出了一种基于相似性的发动机气路参数长期预测方法及系统，利用GMM模型对任意概率密度函数的无限逼近性能来实现概率密度函数的重构。

为了解决上述技术问题，本发明第一方面，提供了一种基于相似性的发动机气路参数长期预测方法，包括以下步骤：

统计距离计算步骤、通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离；

高斯元生成步骤、利用得到的统计距离和历史轨迹样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合；

降序聚合步骤、通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。

本发明第二方面，提供了一种基于相似性的发动机气路参数长期预测系统，包括：

统计距离计算模块，用于通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离；

高斯元生成模块，用于利用得到的统计距离和历史轨迹样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合；

降序聚合模块，用于通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。

实施本发明的基于相似性的发动机气路参数长期预测方法及系统，具有以下有益效果：经对机队的长期衰退趋势的预测实验验证，本发明相对于自回归滑动平均、反向传播神经网络和传统的基于相似性的预测方法相比具有更高的预测精度；通过在离群样本上的预测实验，也验证了本发明所提出的DOA方法相对于核密度估计等加权聚合方法能够有效地增强算法的泛化能力。

附图说明

图1为加权似然估计方法和核密度估计方法的泛化能力缺陷示意图；

图2为本发明的DBSA-GMM方法结构的概念图；

图3为根据本发明优选实施例的基于相似性的发动机气路参数长期预测方法的流程图；

图4为根据本发明优选实施例的基于相似性的发动机气路参数长期预测方法中降序聚合步骤的DOA算法流程图；

图5为根据本发明优选实施例的基于相似性的发动机气路参数长期预测系统的模块框图；

图6为衰退轨迹样例图；

图7为六种预测方法对DEGT进行预测的平均相对误差结果图；

图8为样本群体之间的统计距离示意图；

图9为样本4及10条相似样本的示意图；

图10为所有SBP方法对样本4的DEGT的预测结果图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

为了弥补SBP方法的不足，本发明提出了一种基于相似性的发动机气路参数长期预测方法，该方法基于距离降序聚合的高斯混合模型(Distance Based SequentialAggregation with Gaussian Mixture Model,DBSA-GMM)实现，简称为DBSA-GMM。与传统的SBP方法所不同的是，DBSA-GMM以另一种方式来利用GMM模型对任意概率密度函数的无限逼近性能来实现概率密度函数的重构：首先，DBSA-GMM对各个历史轨迹和目标轨迹的噪声大小分别进行评估，通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征，并借助得到的逐点地距离特征序列评估目标轨迹与各历史轨迹之间的统计距离；而后，与WLE和KDE方法所不同的是，DBSA-GMM跳过了距离到相似度的转化步骤，直接利用得到的统计距离和历史样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合；最后，通过本发明中新提出的降序聚合(Descending Order Aggregation,DOA)方法将所有的高斯元聚合为一个GMM模型，从而采用一种投票的方式，让所有历史样本以相等的权重参与到目标轨迹的概率密度函数的重构过程当中。

在介绍本发明的DBSA-GMM的具体方法之前，首先明确如下的术语：

循环——发动机从起飞—爬升—巡航—降落的一次飞行任务的整个工作过程；

特征——经过标准化处理的表征发动机性能状态的多元性能状态参数向量；

特征元素——特征的某个单一维度，表示经过标准化的单一性能状态参数；

轨迹——在发动机衰退过程中，特征在多元参数空间中的演化路径。

与所有的SBP方法一样，DBSA-GMM方法的适用条件如下：

1、样本库中含有足够多的样本，并且衰退轨迹应该具有足够多的时间点的特征采样数据；

2、每一条衰退轨迹的噪声随发动机性能衰退过程不应该有太大的变化，特征元素的噪声可以认为是高斯噪声且之间相互独立；

3、由于较高纬度的向量之间的相似性难以度量，性能状态特征的维度不能过高(一般不应超过10维)。

图2给出了本发明的DBSA-GMM方法结构的概念图。其中，纵轴为性能特征衰退轨迹在单维参数坐标轴上的投影。如图所示，该方法将在目标轨迹涵盖了足够多的观测数据时被激活。图中S1为统计距离计算步骤，计算目标轨迹与历史轨迹之间的统计距离；图中S2为高斯元生成步骤，S3为降序聚合步骤，共同实现基于统计距离的高斯元聚合。DBSA-GMM首先计算目标轨迹与历史轨迹之间的统计距离d_i(X_o,X_j)，而后采用DOA方法，将每个历史轨迹所提供的假想高斯元整合起来形成最终的概率密度函数。当目标轨迹的新的观测值出现后，DBSA-GMM会对其预测进行刷新预测结果，以不断地提高预测精度。当目标轨迹的走到终点时，其会被归并入历史数据库当中，作为历史样本支持新的目标轨迹的预测，从而形成一种样本循环使用方式，以增量化的形式不断提升DBSA-GMM的预测性能。

请参阅图3，为根据本发明优选实施例的基于相似性的发动机气路参数长期预测方法的流程图。如图3所示，该实施例提供的基于相似性的发动机气路参数长期预测方法包括以下步骤：

首先，在步骤S1中，执行统计距离计算步骤，通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离。

随后，在步骤S2中，执行高斯元生成步骤，利用得到的统计距离和历史轨迹样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合。

最后，在步骤S3中，执行降序聚合步骤，通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。

优选地，统计距离计算步骤S1包括：

(1)通过以下公式计算目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列：

其中，为距离特征量，Φ(X_o,X_i)为距离特征的时间序列；其中M维的目标轨迹和第i个历史轨迹样本分别为 和分别为和的第p个特征元素，和分别为高斯分布的标准差；

(2)将Φ(X_o,X_i)作为服从M自由度的非中心卡方分布，通过求解以下非中心卡方分布的概率密度函数中作为统计距离：

式中，z为时间序列的长度，v为卡方分布的自由度，λ为表示所有高斯元平方和的卡方分布非中心参数,I_v/2-1(·)表示修正的第一类贝塞尔函数。更优选地，该统计距离计算步骤中通过以下公式求解统计距离

式中，L(Φ|M,λ)′为经样本权重重新分配后的对数似然函数，为距离特征，w_t为的权重，M为自由度。优选地，高斯元生成步骤中通过以下函数生成目标轨迹的假想的高斯函数形式的概率密度估计：

其中，为样本i生成的假想估计；和分别为成对的两个高斯元；和分别表示两个高斯元和的权重；N(·)表示高斯函数，为样本i的第p个元素， 为的标准差，为根据X_i所获取的的标准差。

下面对本发明方法各个步骤的具体实现过程进行详细介绍。

1、统计距离计算步骤

1.1距离度量定义

式中，N(·)表示高斯分布。由此可以推得：

众所周知，欧氏距离是一种被广泛应用于表达高维空间中点之间的距离的度量。参照欧氏距离的数学形式

将式(2)的两边同时平方，并沿维度方向对获得的各个元素进行求和运算，可以得到一个具有单一维度的距离特征量。当沿时间轴方向对逐点的数据进行同样的运算，便可以得到距离特征的时间序列，这个距离特征序列传递了两个多元轨迹的逐点距离信息。

式中为距离特征量，Φ(X_o,X_i)为距离特征的时间序列。

结合之前提出的噪声独立假设，若X_o≡X_i，则Φ(X_o,X_i)服从M自由度的中心卡方分布，否则，Φ(X_o,X_i)将服从M自由度的非中心卡方分布。

非中心卡方分布的概率密度函数如下：

式中，v为卡方分布的自由度，λ为表示所有高斯元平方和的卡方分布非中心参数,I_v/2-1(·)表示修正的第一类贝塞尔函数，具有如下形式

其中，Γ(·)表示欧拉伽马方程，其表达式为：

观察式(3)至式(5)可以看出，非中心参数λ实际上在统计意义上传递了两个具有高斯噪声的多元时间序列之间的距离特征序列的平方的期望信息。因此，本发明定义为两个多元时间序列之间的统计距离。

1.2统计距离与其他距离度量的比较

作为度量多元时间序列相似程度的关键要素，距离度量的选择乃至构造方法一直是该领域的研究热点。除了被最广泛采用的欧氏距离和马氏距离之外，皮尔逊(Pearson)相关系数以及相关距离，均方根距离，短期序列距离等度量方式也在各个场合得到了应用并表现出了不同的适应范围。和这些距离度量方式相比，本发明中提出的统计距离对样本的信噪比之间的差异性具有更高的区分度，同时也满足标准距离定义所包含的非负性、对称性、自反性的特性。本发明所提出的统计距离不满足三角不等式，然而，根据现有文献所证明的，在对多元时间序列相似程度的排序与欧氏距离具有等价意义，而作为的统计形式，也具有类似的数理特性。从数学角度上来说，统计距离实际上是一种伪距离，其能够根据样本噪声的大小动态判断样本所包含的信息的可参考程度，从而弱化含有较大噪声的样本的影响。

1.3统计距离求解

为了根据获取的特征距离序列Φ求取非中心参数λ，本发明对多种方法进行了参考，事实证明，当λ和ν皆未知时，λ的求解需要经历一个非常耗时的逼近-寻优的迭代过程。在本发明中，因为已知ν＝M，提出了一种λ的简化求解方法，该方法为无导数优化过程，极大地降低了求解运算量：

采用极大似然估计的方式对λ进行求解，λ的对数似然函数为

已知ν＝M，则ln(L(Φ|v,λ))相对于λ的导数为

通过对式(6)至式(7)的观察，可以发现，f_X(·)具有球对称的特征，即是说，f_X(·)的形状对于各个高斯元的非中心程度的贡献不具备区分能力，而仅仅与它们的贡献的和相关。因此可以将卡方分布的元素重新分解为M个正态分布，并暂且认为：

式中第一项为平方后符合具有一个自由度的非中心卡方分布的正态分布，该正态分布的均值为第2～M项为标准正太分布。因为卡方分布的和仍然为卡方分布，所以f_X(·)可以看作一个由第一项平方生成的自由度为1的非中心卡方分布与第2～M项平方和生成的自由度为M-1的中心卡方分布的和，其中，根据概率密度函数的数学意义，前者的表达式可以写为

式中φ(·)表示标准正态分布。而后，将两个卡方分布的概率密度函数做卷积，f_X(z；M,λ)即可重构为：

式中，当M已知时，C为常数。与之对应地，式(9)相对于λ的导数可以写作

式中，

这里值得一提的是，在时间序列分析研究当中，通常认为时间距离较近的特征值具有更强的代表性，因此，随着时间久远程度的提升，往往会赋予观测值相对较小的权重。为了照顾到这一特性，本发明中对的权重加以调整，令其权重为：

式中，ρ是用以调整权重衰减速率的控制参数。相应的，式(13)被改写为：

通过观察式(15)，可以发现求和运算当中的第二项相对于λ来说是单调递减函数(证明过程较为简单，此处不予给出)。因此，本发明采用斐波那契方法，便能够完成λ的快速求解。

2、高斯元生成步骤

为了表述方便，用Θ＝{X_i},i∈{1,2,…,N}表示整个历史样本群体，表示包含目标轨迹与所有历史轨迹之间的统计距离的距离矢量。其中N为历史样本的个数，X_o表示目标样本并且T表示X_o的长度且所有的统计距离都是根据t∈[1,T]时刻观察到的数据得到的。而后，为了根据第i个历史轨迹估算目标轨迹的第p个特征元素在t^*＞T时刻的均值和方差(假设第p个特征元素在t^*＞T时刻符合高斯分布)，可以根据前文所做出的噪声稳定假设，列出下式：

式中，和分别为根据t∈[1,T]观测值所得到的X_i和X_o第p个特征元素的标准差的无偏估计，为的标准差。和都可以通过标准的噪声幅值估计手段获得。是根据X_i所获取的的标准差。根据每个X_i的噪声变化的不同，其所对应的也会发生不同的变化。

因此，在已知，并且设定与λ_T(X_i,X_o)成正比例的前提下，如果假定与整个样本群体所组成的轨迹束的分散度成比例变化，便可以推得：

式中，分母M表示算法默认各个特征元素对统计距离的贡献度相等，S_p(·)的表达式为：

用于表达整体样本库所组成的轨迹束的分散度。式中，z表示时间序列的长度。

因此，以为例，从第i个历史样本的观察角度来说，通过λ_T(X_i,X_o)和可以获得一个包含两个等权重的高斯元的的假想概率密度分布函数：

式中，和分别表示两个高斯元的权重，而在没有其他外界条件信息的情况下，可以认为同样的，根据其他的历史轨迹样本和与之所对应的统计距离，并分散各个高斯元的权重，可以获取一个包含2N个等权重高斯元的集合观察式(19)可以得知，即使省略进一步的优化和验证等因素的考虑，的概率密度函数也仍然可以被粗略地估计为：

而只要Θ包含足够多的样本，就会具有足够的可靠性。

3、高斯元的降序聚合步骤

观察式(19)可以得知，每一个X_i都会产生两个关于其自身对称的成对高斯元分支。然而，从直觉上来说，在实际场景中X_o的位置一定是唯一的，因而不可能存在其处于另一个高斯元两侧的情形。因此可以推断，每个成对高斯元其中的一个一定是算法的在非中心参数开方过程所生成的虚构分支。为此本发明提出了DOA方法用以剔除式(20)中所包含的虚假高斯元，使得到的概率密度函数更为紧凑，并提高预测精度。DOA通过计算各个高斯元与排除其自身所在的历史样本的剩余样本群体生成的GMM分布的相容性来实现成对高斯元之间的可置信度的对比。剩余样本群体所生成的GMM分布为：

单独高斯元与剩余样本群体生成的GMM的相容性计算方法根据McLachlan所提出的概率密度函数在积分空间中的相关系数计算公式推导得来：

DOA的流程图在图4中给出。如图4所示，优选地，本发明中降序聚合步骤S3包括：

首先，在步骤S301中，流程开始；

随后，在步骤S302中，输入数据，包括：高斯元生成步骤中所获得的序列、各个所对应的与目标轨迹的统计距离

随后，在步骤S303中，将高斯元生成步骤中所获得的序列按照它们与目标轨迹的统计距离进行降序排列；

随后，在步骤S304中，初始化，令i＝1；

随后，在步骤S305中，分别计算各个和相对应的排除自身的混合高斯估计

随后，在步骤S306中，分别计算相容系数：

随后，在步骤S307中，根据相容系数更新成对高斯元的权重：

随后，在步骤S308中，检测是否i＝N，是则转步骤S310，否则转步骤S309；

在步骤S309中，令i＝i+1，并转步骤S305，继续执行直到所有被更新；

在步骤S310中，根据以下公式输出的最终概率密度函数：

最后，在步骤S311中，流程结束。

DOA算法的本质是：对应于更大统计距离的高斯元对的相对距离更大，因而具有更高的分离度而变得容易取舍。对其进行优先区分，能够及早地排除具有更大偏离程度的高斯元，并为后续的对应于较低统计距离的高斯元对的相容性评估过程提供更为可靠的参考值。

通过DOA算法，本发明提供的DBSA-GMM方法能够给出概率形式的预测结果，而为了能够支持需要确定性预测结果的PHM业务，也为了方便与其他方法进行性能对比，DBSA-GMM方法的确定性格式的预测结果可以以所有时刻点的概率密度函数的顶点所连成的轨迹为准，与此同时，采用中位数或者是期望值也是可以可行的。

请参阅图5，为根据本发明优选实施例的基于相似性的发动机气路参数长期预测系统的模块框图。如图5所示，该实施例提供的基于相似性的发动机气路参数长期预测系统100包括：统计距离计算模块110、高斯元生成模块120和降序聚合模块130。

其中，统计距离计算模块110用于通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离。该统计距离计算模块110与前述方法中统计距离计算步骤S1的实现过程相同。

高斯元生成模块120用于利用得到的统计距离和历史轨迹样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合。该高斯元生成模块120与前述方法中高斯元生成步骤S2的实现过程相同。

降序聚合模块130用于通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。该降序聚合模块130与前述方法中降序聚合步骤S3的实现过程相同。

优选地，统计距离计算模块110包括：第一处理单元和第二处理单元。其中：

第一处理单元用于通过以下公式计算目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列：

其中，为距离特征量，Φ(X_o,X_i)为距离特征的时间序列；其中M维的目标轨迹和第i个历史轨迹样本分别为 和分别为和的第p个特征元素，和分别为高斯分布的标准差；

第二处理单元用于将Φ(X_o,X_i)作为服从M自由度的非中心卡方分布，通过求解以下非中心卡方分布的概率密度函数中作为统计距离：

式中，z为时间序列的长度，v为卡方分布的自由度，,λ为表示所有高斯元平方和的卡方分布非中心参数,I_v/2-1(·)表示修正的第一类贝塞尔函数。更优选地，该第二处理单元通过以下公式求解统计距离

式中，L(Φ|M,λ)′为经样本权重重新分配后的对数似然函数，为距离特征，w_t为的权重，M为自由度。

优选地，高斯元生成模块120通过以下函数生成目标轨迹的假想的高斯函数形式的概率密度估计：

其中，为样本i生成的假想估计；和分别为成对的两个高斯元；和分别表示两个高斯元和的权重；N(·)表示高斯函数，为样本i的第p个元素， 为的标准差，为根据X_i所获取的的标准差。

优选地，降序聚合模块130包括：降序排列单元和迭代输出单元；其中，降序排列单元用于将高斯元生成步骤中所获得的序列按照与目标轨迹的统计距离进行降序排列。迭代输出单元用于从i＝1至T，依次计算所有高斯元的权重和并根据以下公式输出的最终概率密度函数：该迭代输出单元的具体迭代流程如图4中步骤S304～S310。

下面通过试验来比较本发明方法与其他几种方法在航空发动机性能预测方面的应用。

1、试验对象

为了验证本发明DBSA-GMM方法的优越性，本发明采用来自于一个由同型号发动机组成的机队的数据来作为验证实验的样本。该型号发动机为民航大涵道比双涵道发动机。在飞行过程中，航空发动机通过进气道吸入空气，经低压压气机(Low PressureCompressor,LPC)、高压压气机(High Pressure Compressor,HPC)对空气进行增压后，进入燃烧室(Combustion Chamber,CC)。高压空气在燃烧室中与雾化油液进行混合并燃烧，产生的高压高温燃气依次经过高压涡轮(High Pressure Turbine,HPT)和低压涡轮(LowPressure Turbine,LPT)。高压涡轮从燃气中抽取功率以驱动高压压气机，以形成核心机的稳定工作循环。低压涡轮使用抽取出的功率驱动同样安装在低速轴上的风扇(Thrust Fan，FAN)和低压压气机，完成外涵道气体的压缩和加速，而从低压涡轮排出的气体具有较高速度，同样也能提供少部分推力。从内涵道和外涵道喷出的气体经过压缩后会具有较高的喷射速度，使得发动机能够借助气体的反冲力来获得推力。

发动机的部件通常工作在极为严酷的条件下。举例来说：发动机的高压压气机的叶片通常为细长的悬臂梁结构，在其对空气进行压缩的过程中需要承受巨大的气动弯矩和强烈的气流颤振，这些因素极易造成叶片的疲劳损伤；高压涡轮的叶片需要在超过10000rpm和超过1200℃的条件下工作，在这样的条件下，即使不考虑气动弯矩的影响，巨大的离心力和金属的热疲劳也极易导致叶片的蠕变、烧蚀、掉块和变形。逐渐累积的部件损伤和变形会使发动机整机的性能逐渐退化，使得发动机整体推力下降、耗油率升高。

为了实现性能状态的实时监控，现今发动机都装有机内测试系统。机内测试系统能够调用传感器测试发动机的关键的气路参数，将其整理并通过无线传输手段试试返回到地面管理中心。所监控的气路参数一般为温度(T)，压力(总压P和静压P_s)，转速(N)等。

2、实验数据及实验安排

因为航空发动机的自身的复杂性，其气路参数的原始值通常具有很高的非线性特征并且与其实际的运行工况强烈地耦合在一起。相比较直接使用原始值的做法而言，当今的航空发动机性能判定方法通常都会采用GPA方法对原始数据进行预处理，并采用预处理后提取的特征来作为性能判定依据。GPA方法的工作过程如下：首先，GPA方法会以健康状态的发动机在整个飞行包线内的气路参数测量值为基准，绘制出表现发动机的各个气路参数相对于所有的运行工况和控制指令的性能“基线”。在进行发动机性能评估时，GPA将待评测发动机的实测气路参数与性能基线相对比，并将得到的偏差值作为最终评定发动机健康程度的特征参数。借鉴全球三大航空发动机生产商通用电气公司(General ElectricCompany,GE)、罗尔斯.罗伊斯公司(Rolls Royce,RR)、普拉特惠特尼公司(Pratt&Whitney,PW)所共同采用的评价指标选取方案，本发明采用发动机排气温度偏差值(Delta ExhaustGas Temperature,DEGT)、燃油流量偏差值(Delta Fuel Flow,DFF)、核心机转速偏差值(Delta Core Speed,DN2)来组成航空发动机的健康状态的特征向量。DEGT、DFF、DN2的单位分别为℃、百分比和百分比。

最初选取的样本总共包含120条完整的发动机气路参数偏差值衰退轨迹。为了保证所有的衰退过程都处于同一起点，本发明首先对样本进行了筛选，择取了自新服役状态或大修后发动机开始衰退的历史样本，最终获得了77条衰退轨迹。这些轨迹平均长度约为4000个飞行循环，照此估算，整个历史样本库总共包含77条轨迹和300000条以上的飞行循环数据，对于维度为3的特征空间内的轨迹相似性比对属充足比较充足的研究案例。图6给出了其中的两条轨迹的气路偏差值数据作为说明范例。颜色条所标注的是数据采集所位于的时间点相对于其自身所在的整个衰退轨迹的时间长度的比值。如图6所示，航空发动机的性能特征值的演变轨迹并非是平滑的曲线，而往往呈现出有向点云的形式。与此同时，不同的衰退轨迹的曲折程度、噪声大小也具有较大的差异性。这种类型的长期预测任务对几乎所有的传统预测方法都是几乎不可能完成的，对SBP类方法也具有相当的难度，而DBSA-GMM方法却比较擅长处理此类问题。在开始预测之前，文中首先采用高斯核平滑方法对各条轨迹加以平滑，提取出的平滑轨迹将用于预测精度的计算，提取出的数据噪声将用来支持DBSA-GMM的进一步计算。与此同时，为了验证本发明DBSA-GMM所做出的高斯噪声假设是否得到满足，本发明还采用了Kolmogorov-Smirnov方法来验证数据噪声是否服从高斯分布。

作为性能对比，自回归滑动平均(Auto-Regressive Moving Average,ARMA)，反向传播神经网络(Back Propagation Artificial Neural Network,BP-ANN)同样被用来预测航空发动机的健康状态参数；与此同时，参照文献，将统计距离和KDE方法相结合，本发明构造了SD-KDE方法并投入到发动机性能状态预测实验，专以用来验证DOA方法能够加强预测算法泛化能力的设想，SD-KDE方法给出的特征元素的在未来各时间点的概率密度函数为

为了验证文中所提出的统计距离相对于传统SBP方法所采用的距离度量的优越性，也同时为了验证DBSA-GMM方法相对于传统SBP方法的整体优越性，文献1和文献2当中的方法也同样被应用于预测实验当中用作算法性能对比。

文中采用交叉实验的方法，即77个样本轮流作为目标样本，被剩余的76个作为历史样本的衰退轨迹所训练的模型预测的方式进行试验。预测算法共执行77次。鉴于最短的衰退轨迹仅包含1916个数据点，文中采用前1200个样本点来训练各个预测模型，以后续的716个样本点作为验证样本。

3、试验结果

本发明只给出了DEGT的预测效果图，预测结果如图7所示。而对DEGT、DN2和DFF的预测精度以数字的形式在表1中给出。图7和表1当中的数字和曲线都是77组预测实验的均值。表1中的粗斜体数字表示所有预测方法针对单一参数进行预测的最优结果。同时，为了对比方法的泛化能力，本发明对所有的历史衰退轨迹两两之间的统计距离进行了统计，统计结果如图8所示。根据图中每一列的数值的平均值，选取了具有最大的特异性的10个样本作为离群样本，并对各类方法在离群样本群体上的预测结果加以统计以评估方法的泛化能力。对于离群样本的预测结果也在图7和表1中给出。

表1所有预测方法的平均相对误差

正如之前的SBP预测方法研究案例所体现的一样，图7中，ARMA和BP-ANN方法在预测初始的30个循环内都取得了较高的预测精度，然而，当滚动预测继续执行时，预测误差开始急剧上升，迅速地超过了SBP类预测方法。这是因为ARMA和BP-ANN将趋势分析焦点过分地集中于时间序列的短期动态特性而只能进行短期预测所导致的。由于ARMA和BP-ANN需要借助于滚动预测过程来预测数据的长期变化趋势，致使误差的持续累积，最终突破了算法的收敛范围，使算法最终预测失败。反观之，得益于相似的基本算法结构，无论是传统的SBP方法还是SD-KDE和DBSA-GMM，都给出了贴合度较高的预测结果。不仅如此，文中所提出的统计距离也确实体现出了比传统的相似性指标更好的轨迹相似行判断效果：根据图7可以得知，传统的SBP方法确实能够给出较为平稳的预测趋势，但是当预测的时间范围较短时，传统的SBP预测方法并没有取得比预测时间范围较长时更好的预测结果，而当对离群样本进行预测时，传统SBP方法给出的预测误差随预测时间范围的验证呈现出下降趋势。根据预测时间范围越长，随机因素的影响越大的假设，预测误差应该随预测时间范围的变长而上升，而传统SBP方法得出的结果却与之相反。这说明传统SBP方法对衰退轨迹的相似程度产生了一定的误判，从而导致预测结果出现了较大的偏差。另一个值得注意的现象是，文献2考虑到了时间久远程度对逐点Mahalanobis权重的影响，因而取得了好于文献1的结果。相较而言，SD-KDE和DBSA-GMM都给出了更为精确的预测结果，而且给出的随预测范围缓慢升高的误差也更符合之前做出的理论假设。两种方法在整个样本群体上的预测精度极为相近，而本发明的DBSA-GMM方法却在离群样本的预测实验中取得了更为精确的预测结果，这验证了DOA对预测算法的泛化能力提升产生的实际作用。

为了进一步验证DOA算法的优势，本发明根据图8选取了特异性最强的离群样本(样本4)进行了单样本预测，此次试验仍然采用初始的1200循环的数据进行模型训练，并将预测范围延伸到了预测起始点之后的2800循环。图9给出了与目标轨迹最为接近的10个历史样本。从图中可以看出，除了被安装在同一架飞机上的另一台发动机产生的衰退轨迹之外，其他9条衰退轨迹与目标轨迹的相似程度都很低。由于相似样本的不足，传统SBP方法难以发挥作用，而统计距离和DOA驱动的DBSA-GMM则能够较好地应对此类问题。

图10给出了所有SBP方法对样本4的DEGT的预测结果。图中包含了DEGT的原始值、平滑值、所有SBP方法的预测值，以及DBSA-GMM方法给出的95％置信区间。由图可知，文献1和文献2中给出的方法在较短的预测范围内表现出了较大的预测误差，而这个现象也与图7所绘制出的平均相对误差曲线相吻合。这是因为对于传统的SBP方法而言，最相似的轨迹将被分配为近乎为1的权重，而其他轨迹的权重近乎为0，这使得概率密度函数的重构过分倚重于单一历史样本，从而造成了预测结果的偏离。与此同时，由于数据噪声的差异性，传统的相似性度量方法比较容易造成具有不同噪声大小的轨迹之间相似程度的误判，同样也拉低了算法的预测精度。对比传统SBP方法，SD-KDE和DBSA-GMM的预测精度更高。在较短的预测范围内，SD-KDE与DBSA-GMM的预测精度相当，但是值得注意的是由于KDE过程采用了加权平均计算过程，所以输出的预测结果具有明显的偏差，而当预测范围逐渐升高时，SD-KDE的预测偏差变的更为明显，导致了误差的迅速升高。与SD-KDE不同的是，DBSA-GMM的预测结果与实际的DEGT发展趋势一直具有较高的贴合度，并且其给出的95％置信区间也对实际曲线具有较高的覆盖率。这是因为DBSA-GMM所采用的DOA能够以等权重的方式综合各个历史样本提供的信息以完成概率密度函数的重构，从而最大化地利用了历史样本所能提供的信息。从图7的SD-KDE和DBSA-GMM的性能对比曲线可以看出，二者在整个样本群体的预测精度相近，而对于离群样本DBSA-GMM具有更高的预测精度，这说明对于非离群样本而言，SD-KDE方法的预测精度要高于DBSA-GMM，而这样的现象也符合第一节中所提及的泛化能力讨论部分所提出的SBP方法的预测精度假设。因此，在执行预测任务时，可以先计算目标样本与各个历史样本的统计距离，而后根据其特异性的强弱酌情选择SD-KDE或DBSA-GMM进行预测。

与所有的SBP方法一样，SD-KDE和DBSA-GMM不依赖算法的外延过程来实现滚动预测，因此能够支持并行计算、避免误差积累、节省计算时间，使得SD-KDE和DBSA-GMM能够应用于更广泛的工业场景，例如对化学反应过程和对多种具有相关性的产品的市场预期的预测。

综上所述，为了预测航空发动机的气路参数偏差值的长期衰退趋势，本发明提出一种针对含有较大噪声的衰退轨迹预测的DBSA-GMM方法。通过卡方分析，DBSA-GMM首先估计衰退轨迹与历史衰退轨迹之间的统计距离，而后采用新提出的DOA方法对历史数据加以聚合，最终得到目标衰退轨迹的未来较长时间范围内的概率密度函数。在航空发动机的性能趋势预测实验当中，DBSA-GMM相对于ARMA，ANN和原始SBP方法给出了更为精确的预测结果，而其泛化能力也在传统SBP方法的基础上得到了强化。

为了追求更高的计算效率，DBSA-GMM在重构目标轨迹的概率密度函数时采用了将统计距离平均分配在各个特征维度上的分解策略，省略了特征之间联合概率密度分布的估算。然而，当预测任务具有更高的精度需求，或在计算资源充足时，这部分损失的精度仍然能够通过期望最大算法(Expectation-Maximization)加以恢复。由于航空发动机的健康状态的衰退实际上是局部的结构损伤所导致的，因此可以预计，结合结构寿命分析和DBSA-GMM方法的剩余寿命预测应该能够进一步提高算法的预测精度。

Claims (10)

1.一种基于相似性的发动机气路参数长期预测方法，其特征在于，包括以下步骤：

统计距离计算步骤、通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离；

高斯元生成步骤、利用得到的统计距离和历史轨迹样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合；

降序聚合步骤、通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。

2.根据权利要求1所述的基于相似性的发动机气路参数长期预测方法，其特征在于，所述统计距离计算步骤包括：

(1)通过以下公式计算目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列：

其中，为距离特征量，Φ(X_o,X_i)为距离特征的时间序列；其中M维的目标轨迹和第i个历史轨迹样本分别为 和分别为和的第p个特征元素，和分别为高斯分布的标准差；

(2)将Φ(X_o,X_i)作为服从M自由度的非中心卡方分布，通过求解以下非中心卡方分布的概率密度函数中作为统计距离：

式中，z为时间序列的长度，v为卡方分布的自由度，,λ为表示所有高斯元平方和的卡方分布非中心参数,I_v/2-1(·)表示修正的第一类贝塞尔函数。

3.根据权利要求2所述的基于相似性的发动机气路参数长期预测方法，其特征在于，所述统计距离计算步骤中通过以下公式求解统计距离

式中，L(Φ|M,λ)′为经样本权重重新分配后的对数似然函数，为距离特征，w_t为的权重，M为自由度。

4.根据权利要求1所述的基于相似性的发动机气路参数长期预测方法，其特征在于，所述高斯元生成步骤中通过以下函数生成目标轨迹的假想的高斯函数形式的概率密度估计：

其中，为样本i生成的假想估计；和分别为成对的两个高斯元；和分别表示两个高斯元和的权重；N(·)表示高斯函数，为样本i的第p个元素，

为的标准差，为根据X_i所获取的的标准差。

5.根据权利要求4所述的基于相似性的发动机气路参数长期预测方法，其特征在于，所述降序聚合步骤包括：

(1)将高斯元生成步骤中所获得的序列按照与目标轨迹的统计距离进行降序排列；

(2)初始化，令i＝1；

(3)分别计算各个和相对应的排除自身的混合高斯估计

(4)分别计算相容系数：

(5)根据相容系数更新成对高斯元的权重：

(6)重复步骤(3)至(5)直到所有被更新；

(7)根据以下公式输出的最终概率密度函数：

6.一种基于相似性的发动机气路参数长期预测系统，其特征在于，包括：

统计距离计算模块，用于通过计算得到目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列，并利用得到的逐点距离特征的时间序列评估目标轨迹与各个历史轨迹之间的统计距离；

高斯元生成模块，用于利用得到的统计距离和历史轨迹样本，针对于每个预测时间点上的单个特征元素，使每个历史轨迹样本都生成一个目标轨迹的假想的高斯函数形式的概率密度估计，成为一个假想高斯元集合；

降序聚合模块，用于通过降序聚合方法对获得的假想高斯单元集合进行聚合，得到目标特征的高斯混合模型。

7.根据权利要求6所述的基于相似性的发动机气路参数长期预测系统，其特征在于，所述统计距离计算模块包括：

第一处理单元，用于通过以下公式计算目标轨迹和各个历史轨迹之间的逐点距离特征的时间序列：

其中，为距离特征量，Φ(X_o,X_i)为距离特征的时间序列；其中M维的目标轨迹和第i个历史轨迹样本分别为 和分别为和的第p个特征元素，和分别为高斯分布的标准差；

第二处理单元，用于将Φ(X_o,X_i)作为服从M自由度的非中心卡方分布，通过求解以下非中心卡方分布的概率密度函数中作为统计距离：

式中，z为时间序列的长度，v为卡方分布的自由度，,λ为表示所有高斯元平方和的卡方分布非中心参数,I_v/2-1(·)表示修正的第一类贝塞尔函数。

8.根据权利要求7所述的基于相似性的发动机气路参数长期预测系统，其特征在于，所述第二处理单元通过以下公式求解统计距离

式中，L(Φ|M,λ)′为经样本权重重新分配后的对数似然函数，为距离特征，w_t为的权重，M为自由度。

9.根据权利要求6所述的基于相似性的发动机气路参数长期预测系统，其特征在于，所述高斯元生成模块通过以下函数生成目标轨迹的假想的高斯函数形式的概率密度估计：

其中，为样本i生成的假想估计；和分别为成对的两个高斯元；和分别表示两个高斯元和的权重；N(·)表示高斯函数，为样本i的第p个元素，

为的标准差，为根据X_i所获取的的标准差。

10.根据权利要求9所述的基于相似性的发动机气路参数长期预测系统，其特征在于，所述降序聚合模块包括：

降序排列单元，用于将高斯元生成步骤中所获得的序列按照与目标轨迹的统计距离进行降序排列；

迭代输出单元，用于从i＝1至T，依次计算所有高斯元的权重和并根据以下公式输出的最终概率密度函数：