本发明涉及系统可靠性领域，更具体的涉及复杂设备的剩余寿命预测方法。
背景技术：
由于设备不同部件的多个退化过程包含了相同的退化信息，这多个退化过程往往具有一定的依赖关系。因此对于具有多个部件的大型复杂系统，如何同时对其多个部件的退化情况进行监测和状态估计，并且合理的描述多个退化过程之间的依赖关系，从而实现基于多退化指标的设备剩余寿命(remainingusefullife,rul)预测具有重要的工程意义。
通过对目前的文献检索发现，现有技术大多仅对单一的退化指标进行分析，仅有少部分文献同时分析多个退化指标。大部分基于多退化指标的预测方法将多个退化指标之间的依赖关系建立在多个退化过程的内部随机性当中，如xi等人在《remainingusefullifepredictionformulticomponentsystemswithhiddendependencies》中为每个退化指标建立维纳过程的退化模型，并且使用一个扩散系数矩阵(diffusioncoefficientmatrix)将多个退化指标之间的依赖关系建立在维纳过程的随机项中；peng等人所著的《jointonlinerulpredictionformultivariatedeterioratingsystems》一文，以及fang等人所著的《copula-basedreliabilityanalysisofdegradingsystemswithdependentfailures》一文中，均为每个退化指标建立随机过程的退化模型，并且使用copula函数将多个指标之间的依赖关系建立在随机过程模型的随机项中。
由于不同设备的多个退化过程之间可能具有不同的依赖结构，例如，某些设备的不同部件的退化速度之间可能存在一定相关性，这种相关性可能是导致不同退化指标依赖关系的原因。当为一个退化指标建立对应的随机过程模型时，模型中的退化轨迹参数在一定程度上反映了该指标的退化速率，因此针对这种依赖结构，将依赖关系体现在随机过程模型的参数相关性上更为合理。然而当前的基于多指标的预测方法均将多指标之间的依赖关系建立在随机过程模型的随机项中，对于这种依赖关系体现为不同指标退化速率相关性的多部件复杂设备，尚未出现一种能够合理解释依赖关系的多退化指标设备剩余寿命预测方法，使得对于这类设备的健康状态估计与剩余寿命预测不能满足工程应用的要求。
技术实现要素：
为了克服现有技术的不足，本发明提供一种基于多退化指标的复杂设备剩余寿命预测方法。对多个退化指标建立多元退化模型，将该多元退化模型中的模型参数根据其物理意义分划分外部参数和退化轨迹参数，并且在该多元退化模型中使用协方差矩阵建立多个退化指标的退化轨迹参数之间的相关关系。使用最大似然估计对外部参数进行离线估计，且使用基于核平滑的粒子滤波(kernelsmoothing–particlefilter,ks-pf)方法对模型中的退化轨迹参数与多个指标的退化状态进行在线估计，从而实现基于多退化指标的设备剩余寿命预测。此外，本发明还采用最优调整(optimaltuning,ot)的方法对ks-pf中的核参数进行在线的自适应调整。
本发明解决其技术问题所采用的技术方案的详细步骤如下：
步骤1：边缘退化模型的建立；
假设设备有k个退化指标，为每个退化指标建立维纳过程的边缘退化模型：
xk(t)＝μk(θk,t) σkb(t)(1)
公式(1)为状态转移方程，公式(2)为量测方程，其中，k＝1,2…k代表第k个退化指标，xk(t)是第k个退化指标在t时刻的退化状态，μk(θk,t)是一个单调可微函数，被称为退化轨迹，为线性或非线性函数；θk＝[ak,bk,…]t为退化轨迹参数，决定指标的退化速率，个体不确定性指的是同一批次的不同设备之间的退化多样性，通过对退化指标的状态值进行回归得到，b(t)是一个标准布朗运动过程，σk是波动系数；σkb(t)具有正态增量，即有
令xk,j＝xk(tj)，yk,j＝yk(tj)，δt＝tj-tj-1，其中0≤j≤m，m为待预测设备中退化指标量测点的个数，仅考虑第k个退化指标，则边缘退化模型的公式(1)和公式(2)被转化为粒子滤波方法需要的状态转移方程和量测方程，其中，状态转移方程为f(xk,j|xk,j-1)：
量测方程h(yk,j|xk,j)为：
当公式(3)和(4)被用于多个退化指标时，令xj＝[x1,j,…,xk,j]t表示所有退化指标在tj时刻的状态值，令yj＝[y1,j,…,yk,j]t代表所有退化指标在tj时刻的量测值，则状态转移方程和量测方程转化为f(xj|xj-1)和h(yj|xj)；此时的多元退化模型f(xj|xj-1)和h(yj|xj)未建立多个退化指标之间的依赖关系；
在多个退化指标同时退化的情况中，假设设备的失效时刻是先到达对应失效阈值的某退化指标到达该退化指标的失效阈值的时刻，因此对一个具有k个退化指标的复杂设备，每个退化指标的失效阈值定义为hk,k＝1,2,…,k，该设备在tj时刻开始预测时的剩余使用寿命为：
步骤2：外部参数离线估计；
外部参数
其中，q＝[min{ti,tj}]1≤i,j≤m，im是一个m阶单位矩阵；
外部参数的似然函数如下：
为了估计外部参数，首先将退化状态
步骤3：测试集设备退化轨迹参数与退化指标状态的同时在线估计；
步骤3.1：基于参数相关性的多元退化模型建立；
假设模型(1)中的退化轨迹参数θ1,θ2,…,θk之间线性相关，则不同退化指标x1(t),x2(t),…,xk(t)之间的复杂依赖关系由线性相关的退化轨迹参数以及θk和xk(t)之间的线性或非线性的映射共同构建和描述；首先估计θ1,θ2,…,θk之间的皮尔逊相关系数，假设退化轨迹被设置为μk(θk,t)＝akexp(bk·t)，θk＝[ak,bk]t，则
使用基于核平滑的粒子滤波方法对测试集设备的多个退化指标及其对应的退化轨迹参数进行在线估计，首先建立基于参数间相关性的状态转移方程，通过在状态转移方程中构建一个协方差矩阵来实现建立参数之间相关性的目的，从一个多元高斯分布中进行采样，得到多个退化指标的退化轨迹参数的粒子群，而该多元高斯分布的协方差矩阵体现了多个退化指标参数之间的相关关系；
步骤3.2：ks-pf实现退化轨迹参数与多退化指标的同时在线估计
使用粒子滤波中的贯序重要性重采样(sequentialimportanceresampling,sir)框架，按照状态转移方程(10)和量测方程(4)，对测试集设备的多指标退化状态x＝[x1(t),…,xk(t)]t和退化轨迹参数θk,j＝[ak,j,bk,j,…]t,k＝1,…,k进行在线估计；
步骤3.3：最优调整ot方法实现核参数s的自适应在线调整；
采用最优调整(optimizedtuning,ot)方法，在每个状态估计时刻tj，均选择一个最优的核参数sj；
在tj时刻，设备多个退化指标的多元退化状态的先验分布qj＝p(xj|y1:j-1)表示为公式(11)，后验分布pj＝p(xj|y1:j)表示为公式(15)，先验分布与后验分布的kl散度计算为：
贝叶斯滤波中的更新公式为：
将公式(18)和公式(11)带入公式(17)，并且使用狄拉克函数δ(·)的挑选性质，得到：
将公式(13)带入公式(14)，得到：
最后，将公式(20)带入公式(19)，得到：
因此，在每个tj时刻，在sir框架中用到的状态转移方程(10)里的核参数s均使用最优的核参数：
用状态转移方程(10)和量测方程(4)，按照粒子滤波的sir框架，即在每个tj时刻均使用公式(13)、(14)得到更新后的粒子集
步骤4：多退化指标的复杂设备剩余寿命预测；
在每个tj时刻均得到粒子集
在tj时刻开始预测时，将估计得到的每个粒子代表的参数和系统状态均代入模型(3)中，则通过设备失效定义得到基于每个粒子预测的剩余使用寿命；基于第i个粒子预测得到的设备剩余使用寿命为：
在tj时刻开始预测时，预测得到的设备的剩余使用寿命的概率密度函数(probabilitydensityfunction,pdf)为：
即求得剩余使用寿命的概率密度函数，即式(25)中的p(lj|y1:j)。
所述步骤3.1的基于参数相关性的多元退化模型建立的具体步骤如下：
令δuj＝[δμ1(tj),…,δμk(tj)]t，ωj＝[ω1,j,…,ωk,j]t，并且令
其中，
当粒子
所述步骤3.2中ks-pf实现退化轨迹参数与多退化指标的同时在线估计的具体步骤如下：
1)初始化：
令j＝0，设置t0时刻的设备的多指标退化状态和退化轨迹参数的初始分布p(x0)和p(θ0)；
2)状态估计：
令j加1，若存在量测向量yj，顺序执行以下步骤，否则终止该程序；
a)预测：从建议分布p(xj,θj|xj-1,θj-1)中按照状态转移方程(10)采样得到粒子集
b)更新：使用tj时刻新到来的退化指标量测向量yj按照量测方程(4)对粒子权值进行更新：
并且该权值被归一化为：
c)重采样：有效粒子数被计算为
d)迭代执行步骤a)-步骤c)，即得到状态估计。
所述步骤2中，退化状态
所述
本发明的有益效果在于通过为多个退化指标建立随机过程的多元退化模型，并将多个退化指标之间的依赖关系建立在模型退化轨迹参数的相关性上，实现了对多退化指标的合理建模；使用最大似然方法实现了模型中外部参数的离线估计；构建了一种最优调整核参数的基于核平滑的粒子滤波(otks-pf)方法，从有噪声的量测数据中对多退化指标的真实退化状态以及退化轨迹参数进行在线估计；最后基于最后一个时刻的多退化指标状态估计的概率分布，预测得到设备的剩余使用寿命。
附图说明
图1是本发明的实现框架图。
图2是第一组实验中训练集第1号发动机的两个退化指标的退化状态及退化状态量测值。其中图2(a)为训练集第1号发动机的量测值及退化指标1，图2(b)为训练集第1号发动机的量测值及退化指标2。
图3是第一组实验中所有训练集发动机的退化轨迹参数箱线图。
图4是第一组实验中训练集第1号和8号发动机的退化状态及模型拟合曲线。
图5是第一组实验中所有训练集发动机的退化轨迹参数散点图。
图6是第一组实验中测试集第5号发动机的退化轨迹参数的在线估计结果。其中图6(a)为采用退化指标1的估计结果，图6(b)为采用退化指标2的估计结果。
图7是第一组实验中测试集第19号发动机分别使用单指标和多指标联合预测得到的rul预测结果。其中图7(a)显示了rul点预测结果及概率密度函数预测结果；图7(b)显示了rul点预测结果及置信区间预测结果。
图8是第一组实验中测试集第68号发动机分别使用单指标和多指标联合预测得到的rul预测结果。其中图8(a)显示了rul点预测结果及概率密度函数预测结果；图8(b)显示了rul点预测结果及置信区间预测结果。
图9是第一组实验所有20个测试集发动机的rul平均预测结果。其中图9(a)显示了预测的平均误差结果，图9(b)显示了预测的平均95％置信区间宽度结果。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
本发明的实现框架如图1所示，由离线模型训练和模型的在线估计两部分构成。在离线模型训练部分，使用一组训练集设备对退化模型进行参数估计，以得到能够反映这组设备的共同特征的模型参数(模型的外部参数)。而在线估计部分，则针对特定的测试集设备，对其进行个性化。具体的，利用该设备的量测信息去在线估计其特有的退化信息(退化轨迹参数和退化指标的状态值)，从而实现对该测试集设备的rul预测。
本发明的具体步骤如下：
步骤1：边缘退化模型的建立
假设设备有k个退化指标，为每个退化指标建立维纳过程的边缘退化模型：
xk(t)＝μk(θk,t) σkb(t)(1)
公式(1)为状态转移方程，公式(2)为量测方程，其中，k＝1,2…k代表第k个退化指标，xk(t)是第k个退化指标在t时刻的退化状态，μk(θk,t)是一个单调可微函数，被称为退化轨迹，为线性或非线性函数；θk＝[ak,bk,…]t为退化轨迹参数，决定指标的退化速率，μk(θk,t)体现了退化过程的个体不确定性以及非线性不确定性；个体不确定性指的是同一批次的不同设备之间的退化多样性，通过对退化指标的状态值进行回归得到，b(t)是一个标准布朗运动过程，σk是波动系数；σkb(t)具有正态增量，即有
假设
令xk,j＝xk(tj)，yk,j＝yk(tj)，δt＝tj-tj-1，其中0≤j≤m，m为待预测设备中退化指标量测点的个数，仅考虑第k个退化指标，则边缘退化模型的公式(1)和公式(2)被转化为粒子滤波方法需要的状态转移方程和量测方程，其中，状态转移方程为f(xk,j|xk,j-1)：
量测方程h(yk,j|xk,j)为：
当公式(3)和(4)被用于多个退化指标时，令xj＝[x1,j,…,xk,j]t表示所有退化指标在tj时刻的状态值，令yj＝[y1,j,…,yk,j]t代表所有退化指标在tj时刻的量测值，则状态转移方程和量测方程转化为f(xj|xj-1)和h(yj|xj)；但此时的多元退化模型f(xj|xj-1)和h(yj|xj)未建立多个退化指标之间的依赖关系；
在多个退化指标同时退化的情况中，假设设备的失效时刻是：先到达对应失效阈值的某退化指标到达该退化指标的失效阈值的时刻，因此对一个具有k个退化指标的复杂设备，每个退化指标的失效阈值定义为hk,k＝1,2,…,k，该设备在tj时刻开始预测时的剩余使用寿命为：
步骤2：外部参数离线估计；
外部参数
其中，q＝[min{ti,tj}]1≤i,j≤m，im是一个m阶单位矩阵；
外部参数的似然函数如下：
训练集设备中
步骤3：测试集设备退化轨迹参数与退化指标状态的同时在线估计；
步骤3.1：基于参数相关性的多元退化模型建立；
对于不同退化指标之间的依赖关系是由其退化速率之间的相关性导致的复杂设备，假设模型(1)中的退化轨迹参数θ1,θ2,…,θk之间线性相关，则不同退化指标x1(t),x2(t),…,xk(t)之间的复杂依赖关系由线性相关的退化轨迹参数以及θk和xk(t)之间的线性或非线性的映射(公式(1)中的线性或非线性函数μk(θk,t))共同构建和描述。因此，首先需要估计θ1,θ2,…,θk之间的皮尔逊相关系数，假设退化轨迹被设置为μk(θk,t)＝akexp(bk·t)，θk＝[ak,bk]t，则
使用基于核平滑的粒子滤波方法对测试集设备的多个退化指标及其对应的退化轨迹参数进行在线估计，首先需要建立合理的基于参数间相关性的状态转移方程，通过在状态转移方程中构建一个协方差矩阵来实现建立参数之间相关性的目的。直接从一个多元高斯分布中进行采样，得到多个退化指标的退化轨迹参数的粒子群，而该多元高斯分布的协方差矩阵体现了多个退化指标参数之间的相关关系。
令δuj＝[δμ1(tj),…,δμk(tj)]t，ωj＝[ω1,j,…,ωk,j]t，并且令
其中：
从公式(10)中看到，当粒子
步骤3.2：ks-pf实现退化轨迹参数与多退化指标的同时在线估计
使用粒子滤波中的贯序重要性重采样(sequentialimportanceresampling,sir)框架，按照状态转移方程(10)和量测方程(4)，对测试集设备的多指标退化状态x＝[x1(t),…,xk(t)]t和退化轨迹参数θk,j＝[ak,j,bk,j,…]t,k＝1,…,k进行在线估计，具体步骤如下：
1)初始化：
令j＝0，设置t0时刻的设备的多指标退化状态和退化轨迹参数的初始分布p(x0)和p(θ0)；
2)状态估计：
令j加1，若存在量测向量yj，顺序执行以下步骤，否则终止该程序；
2.1)预测：从建议分布p(xj,θj|xj-1,θj-1)中按照状态转移方程(10)采样得到粒子集
2.2)更新：使用tj时刻新到来的退化指标量测向量yj按照量测方程(4)对粒子权值进行更新：
并且该权值被归一化为：
2.3)重采样：有效粒子数被计算为
2.4)迭代执行步骤2.1)，即得到状态估计；
步骤3.3：最优调整ot方法实现核参数s的自适应在线调整；
公式(10)中的参数
粒子滤波方法的在线估计准确度与建议分布的选择有关，一个最优的建议分布应该使得更新后的粒子权值方差最小，这也可以理解为：一个好的建议分布应该使得先验分布与后验分布尽可能的相似。kl散度(kullback-leiblerdivergence)是一种度量两个概率分布相似程度的指标，因此ot方法在每个tj时刻，均选择使得先验分布与后验分布的kl散度最小的核参数sj。
在tj时刻，设备多个退化指标的多元退化状态的先验分布qj＝p(xj|y1:j-1)表示为公式(11)，后验分布pj＝p(xj|y1:j)表示为公式(15)，先验分布与后验分布的kl散度计算为：
贝叶斯滤波中的更新公式为：
将公式(18)和公式(11)带入公式(17)，并且使用狄拉克函数δ(·)的挑选性质，得到：
将公式(13)带入公式(14)，得到：
最后，将公式(20)带入公式(19)，得到：
因此，在每个tj时刻，在sir框架中用到的状态转移方程(10)里的核参数s均使用最优的核参数：
这样，使用该部分介绍的otks-pf方法，用状态转移方程(10)和量测方程(4)，按照粒子滤波的sir框架，即可在每个tj时刻均使用公式(13)、(14)得到更新后的粒子集
步骤4：多退化指标的复杂设备剩余寿命预测
根据步骤3.2和步骤3.3描述的“最优调整核参数的核平滑(otks)的基于核平滑的粒子滤波(otks-pf)实现多个退化指标的状态在线估计方法，在每个tj时刻均可得到粒子集
在tj时刻开始预测时，将估计得到的每个粒子代表的参数和系统状态均代入模型(3)中，则可以通过设备失效定义得到基于每个粒子预测的剩余使用寿命；基于第i个粒子预测得到的设备剩余使用寿命为：
在tj时刻开始预测时，预测得到的设备的剩余使用寿命的概率密度函数(probabilitydensityfunction,pdf)为：
即求得剩余使用寿命的概率密度函数，即上式中的p(lj|y1:j)。
1.发动机退化数据介绍与数据预处理
使用nasa发布的c-mapss软件仿真的飞机涡轮发动机的退化数据进行多退化指标方法的验证和分析。该发动机退化数据集包含100个同系列单工况单故障模式的发动机的退化数据，每个发动机均是从某个特定状态一直运行到系统故障。每个发动机的退化状态均通过21个传感器得到的时间序列量测数据来反映。
(1)数据预处理及多退化指标的构造
首先，由于发动机原始的测量信号量纲不一样，因此对原始测量数据进行z-score标准化。之后需要构造发动机的退化指标，原始数据中有21个传感器的量测数据，对原始数据使用主成分分析(principalcomponentsanalysis,pca)进行降维处理以去除噪声数据保留有用信息。
经过pca处理及分析，发现仅需二个维度即可保留原始数据87.6％的信息量，因此将原始数据降为二维，并且将降维得到的二维数据作为退化指标1和退化指标2来共同反映该发动机的退化情况。
对每个退化指标建立如公式(1)和(2)的边缘退化模型，由于在该发动机数据中使用两个退化指标进行设备的rul预测，因此公式(1)和(2)中的k＝1,2代表第k个退化指标。
这样，经过pca降维后得到数据即作为该发动机的退化指标1和退化指标2的量测值yk(t),k＝1,2，再将降维后得到的退化指标1和退化指标2的退化状态量测值均进行简单滑动平均滤波，得到相对更平滑的曲线，滑动窗口大小设置为25；将滤波后的曲线作为发动机的两个退化指标的退化状态值xk(t),k＝1,2。图2展示了1号发动机通过pca降维后得到的两个退化指标的量测值yk(t)及经过滑动平均后得到的退化状态值xk(t)。
(2)失效阈值及设备失效的定义
定义hk,k＝1,2分别为两个退化指标的失效阈值。这里，对退化指标1设定失效阈值h1＝1.22，对退化指标2设定失效阈值h2＝0.85。定义一个发动机的真实失效时刻teol为：发动机发生退化时，先退化至其对应的失效阈值的指标的退化状态值xk(t)到达其阈值的时刻。
在实验验证部分，取一部分发动机作为训练集，其他发动机作为测试集，训练集发动机的全寿命的退化数据均视为已知，而测试集发动机则视为仅已知部分退化数据。对于所有100个发动机，采用k-fold来验证我们的方法，k值取5，即将100个发动机随机分为5组，每组20个发动机。将实验重复五次，称为5组实验，每次取其中一组的20个发动机数据做测试集，其余80个发动机为训练集。在每组实验中，均使用训练集数据进行参数估计，并对测试集发动机在不同时刻进行多次的剩余寿命预测。由于篇幅原因，以下分析均以第一组实验为例，我们也给出了所有5组实验的预测结果，并进行了结果对比分析。
2.边缘退化模型的建立及外部参数估计；
对于建立的公式(1)和(2)的边缘退化模型，首先要确定退化轨迹μk(θk,t)的形式。通过对两个退化指标状态值的退化轨迹进行观察，看到其呈现指数形式的退化，因此本文为两个退化指标建立如下的指数模型：
μk(θk,t)＝ak bk·exp(ck·t)(47)
对训练集发动机已知的两个退化指标的量测值yk(t),k＝1,2进行滑动窗口为25的滑动平均处理，得到其退化状态值xk(t),k＝1,2。图3为使用公式(26)的模型对所有训练集发动机的两个退化指标的状态值进行最小二乘拟合后，对拟合的退化轨迹参数进行统计得到的箱线图。因此对参数b1,b2取其各自的中位数作为固定值：b1＝0.085,b2＝0.056。这样，模型(26)变为：
图4展示了第一组实验中的2个训练集发动机(第1号和第8号发动机)使用公式(27)进行拟合时的模型拟合效果。从图4中可以看到，使用公式(27)的退化轨迹模型对两个退化指标的退化状态进行拟合，得到的拟合结果较好。之后，使用训练集数据以及公式(8)和(9)，对退化模型的外部参数
表1第一组实验的外部参数估计结果
3.测试集发动机退化轨迹参数在线估计及剩余寿命预测
按照公式(10)的形式，建立两个退化指标的状态转移方程f(xj,cj|xj-1,cj-1)如下：
其中，k＝1,2代表第k个退化指标。
b1＝0.085，b2＝0.056(对第一组实验)，
从公式(28)可以看到在状态转移的过程中仅包含了状态的增量，因此虽然此时θk＝[ak,ck]t，但在otks-pf的状态参数同估计时，仅需估计参数ck。参数c1和c2之间的皮尔逊相关系数由训练集计算为
使用otks-pf方法对多个退化指标的状态值x1(t),x2(t)及退化轨迹参数c1,c2进行在线估计，状态和参数的初始分布被设置为高斯分布，即：x0～n(μx,σx)，c0～n(μc,σc)，其中
图6展示了在第一组实验中测试集第5号发动机使用该方法进行双指标联合预测及单指标预测的参数估计效果。可以看到，该基于参数相关性的双指标联合预测方法与单指标的预测方法效果相近；在otks-pf的状态参数同估计中，随着已知观测点的增多，估计得到的参数均能逐渐收敛到参数的真实值。
图7和图8分别展示了第一组实验中测试集第19号发动机和测试集第68号发动机的寿命预测结果。其中，位于“rul”轴与“开始预测时刻”轴组成的平面内的黄线代表设备的真实rul值，(a)图中的蓝线，黑线，红线分别代表使用退化指标1，退化指标2以及使用基于参数的联合预测方法预测得到的rul的概率密度曲线，而蓝点，黑点和红点则是对应的预测得到的rul期望值。
图7和图8中的(a)图展示了分别使用退化指标1、退化指标2以及使用基于参数的联合预测方法预测得到的rul的概率密度曲线以及预测值；(b)图展示了基于参数相关性的多指标联合估计的rul预测结果及95％置信区间。观察图7和图8可以得到以下分析及结论：
(1)基于参数相关性的粒子滤波多退化指标设备rul预测方法可以同时估计两个退化指标的状态值，并且跟踪到先到达其失效阈值的退化指标，从而预测得到准确的设备rul的后验概率。从图7和图8的(a)图中可以看到，对于19号发动机，联合估计方法得到的预测pdf随着已知量测点的增加，不断地跟踪到退化指标2预测得到的rul的pdf；而对于第68号发动机，联合估计的预测pdf不断的跟踪到指标1预测得到的rul的pdf。
这是因为对于第19号发动机的两个退化指标，退化指标2首先退化至其对应的失效阈值，该时刻是第19号发动机的真实失效时刻。而对于第第68号发动机，退化指标1首先退化至其对应的失效阈值，该时刻是第68号发动机的真实失效时刻。因此，对于同一批设备，无论其哪个退化指标先退化至其对应的失效阈值，本发明提出的联合预测的方法均能跟踪到该退化指标的预测结果，这证明了本发明提出的基于参数相关性的多指标联合预测方法的有效性。
(2)基于参数相关性的粒子滤波多退化指标设备rul预测方法随着开始预测时刻的增长，其rul预测的精确度不断提升。从图7和图8的(b)图可以看到，随着已知观测点的增多，基于参数相关性的联合估计方法预测得到的rul的95％置信区间的宽度不断减小，预测的精度不断提高，并且预测的准确度也较高。这符合使用粒子滤波实现rul预测的逻辑，因为已知的量测值越多，则得到的先验知识越多，因此预测也会更准确，更精确。这说明了本发明所提出的方法具有实际的应用价值。
为了更直观的对比使用单个退化指标实现rul预测与双指标联合估计实现rul预测的预测效果，在每组实验中，对所有的测试集发动机从其寿命终止时刻teol的三分之二时刻开始预测，之后每隔1个时刻点预测一次。这样，每个测试集发动机均能得到若干个寿命预测的结果，再对该结果取平均。图9展示了第一组实验的预测误差与95％置信区间的平均值。由图9可以看出，在第一组实验的20个测试发动机中，基于参数相关性的联合估计的预测误差基本低于退化指标1和退化指标2的单个估计的预测误差，这符合分析的结果。并且联合估计的95％置信区间宽度也均低于单指标估计得到的95％置信区间宽度。
表2五组实验的单指标预测与多指标联合预测结果
对于每一组实验，对所有20个发动机的平均预测误差和平均95％置信区间宽度算平均值，表2展示了所有5组实验的预测结果对比情况。从表2中可以看到，当考虑同一批次具有个体差异性的不同测试集设备时，本专利提出的基于参数相关性的多退化指标设备rul预测方法在预测准确度(由预测误差体现)和预测精度(由95％置信区间体现)上都比单指标的预测方法有较大的优势，这说明本发明提出的方法能够对具有多个退化指标的大型复杂设备进行有效且准确的rul预测。
1.一种基于多退化指标的复杂设备剩余寿命预测方法，其特征在于包括下述步骤：
步骤1：边缘退化模型的建立；
假设设备有k个退化指标，为每个退化指标建立维纳过程的边缘退化模型：
xk(t)＝μk(θk,t) σkb(t)(1)
公式(1)为状态转移方程，公式(2)为量测方程，其中，k＝1,2…k代表第k个退化指标，xk(t)是第k个退化指标在t时刻的退化状态，μk(θk,t)是一个单调可微函数，被称为退化轨迹，为线性或非线性函数；θk＝[ak,bk,…]t为退化轨迹参数，决定指标的退化速率，个体不确定性指的是同一批次的不同设备之间的退化多样性，通过对退化指标的状态值进行回归得到，b(t)是一个标准布朗运动过程，σk是波动系数；σkb(t)具有正态增量，即有
令xk,j＝xk(tj)，yk,j＝yk(tj)，δt＝tj-tj-1，其中0≤j≤m，m为待预测设备中退化指标量测点的个数，仅考虑第k个退化指标，则边缘退化模型的公式(1)和公式(2)被转化为粒子滤波方法需要的状态转移方程和量测方程，其中，状态转移方程为f(xk,j|xk,j-1)：
量测方程h(yk,j|xk,j)为：
当公式(3)和(4)被用于多个退化指标时，令xj＝[x1,j,…,xk,j]t表示所有退化指标在tj时刻的状态值，令yj＝[y1,j,…,yk,j]t代表所有退化指标在tj时刻的量测值，则状态转移方程和量测方程转化为f(xj|xj-1)和h(yj|xj)；此时的多元退化模型f(xj|xj-1)和h(yj|xj)未建立多个退化指标之间的依赖关系；
在多个退化指标同时退化的情况中，假设设备的失效时刻是先到达对应失效阈值的某退化指标到达该退化指标的失效阈值的时刻，因此对一个具有k个退化指标的复杂设备，每个退化指标的失效阈值定义为hk,k＝1,2,…,k，该设备在tj时刻开始预测时的剩余使用寿命为：
步骤2：外部参数离线估计；
外部参数
其中，q＝[min{ti,tj}]1≤i,j≤m，im是一个m阶单位矩阵；
外部参数的似然函数如下：
为了估计外部参数，首先将退化状态
步骤3：测试集设备退化轨迹参数与退化指标状态的同时在线估计；
步骤3.1：基于参数相关性的多元退化模型建立；
假设模型(1)中的退化轨迹参数θ1,θ2,…,θk之间线性相关，则不同退化指标x1(t),x2(t),…,xk(t)之间的复杂依赖关系由线性相关的退化轨迹参数以及θk和xk(t)之间的线性或非线性的映射共同构建和描述；首先估计θ1,θ2,…,θk之间的皮尔逊相关系数，假设退化轨迹被设置为μk(θk,t)＝akexp(bk·t)，θk＝[ak,bk]t，则
使用基于核平滑的粒子滤波方法对测试集设备的多个退化指标及其对应的退化轨迹参数进行在线估计，首先建立基于参数间相关性的状态转移方程，通过在状态转移方程中构建一个协方差矩阵来实现建立参数之间相关性的目的，从一个多元高斯分布中进行采样，得到多个退化指标的退化轨迹参数的粒子群，而该多元高斯分布的协方差矩阵体现了多个退化指标参数之间的相关关系；
步骤3.2：ks-pf实现退化轨迹参数与多退化指标的同时在线估计
使用粒子滤波中的贯序重要性重采样框架，按照状态转移方程(10)和量测方程(4)，对测试集设备的多指标退化状态x＝[x1(t),…,xk(t)]t和退化轨迹参数θk,j＝[ak,j,bk,j,…]t,k＝1,…,k进行在线估计；
步骤3.3：最优调整ot方法实现核参数s的自适应在线调整；
采用最优调整(optimizedtuning,ot)方法，在每个状态估计时刻tj，均选择一个最优的核参数sj；
在tj时刻，设备多个退化指标的多元退化状态的先验分布qj＝p(xj|y1:j-1)表示为公式(11)，后验分布pj＝p(xj|y1:j)表示为公式(15)，先验分布与后验分布的kl散度计算为：
贝叶斯滤波中的更新公式为：
将公式(18)和公式(11)带入公式(17)，并且使用狄拉克函数δ(·)的挑选性质，得到：
将公式(13)带入公式(14)，得到：
最后，将公式(20)带入公式(19)，得到：
因此，在每个tj时刻，在sir框架中用到的状态转移方程(10)里的核参数s均使用最优的核参数：
用状态转移方程(10)和量测方程(4)，按照粒子滤波的sir框架，即在每个tj时刻均使用公式(13)、(14)得到更新后的粒子集
步骤4：多退化指标的复杂设备剩余寿命预测；
在每个tj时刻均得到粒子集
在tj时刻开始预测时，将估计得到的每个粒子代表的参数和系统状态均代入模型(3)中，则通过设备失效定义得到基于每个粒子预测的剩余使用寿命；基于第i个粒子预测得到的设备剩余使用寿命为：
在tj时刻开始预测时，预测得到的设备的剩余使用寿命的概率密度函数为：
即求得剩余使用寿命的概率密度函数，即式(25)中的p(lj|y1:j)。
2.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述步骤3.1的基于参数相关性的多元退化模型建立的具体步骤如下：
令δuj＝[δμ1(tj),…,δμk(tj)]t，ωj＝[ω1,j,…,ωk,j]t，并且令
其中，
当粒子
3.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述步骤3.2中ks-pf实现退化轨迹参数与多退化指标的同时在线估计的具体步骤如下：
1)初始化：
令j＝0，设置t0时刻的设备的多指标退化状态和退化轨迹参数的初始分布p(x0)和p(θ0)；
2)状态估计：
令j加1，若存在量测向量yj，顺序执行以下步骤，否则终止该程序；
a)预测：从建议分布p(xj,θj|xj-1,θj-1)中按照状态转移方程(10)采样得到粒子集
b)更新：使用tj时刻新到来的退化指标量测向量yj按照量测方程(4)对粒子权值进行更新：
并且该权值被归一化为：
c)重采样：有效粒子数被计算为
d)迭代执行步骤a)-步骤c)，即得到状态估计。
4.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述步骤2中，退化状态
5.根据权利要求1所述的基于多退化指标的复杂设备剩余寿命预测方法，其特征在于：
所述
