CN107451392A

CN107451392A - 一种含有多个相关退化过程的剩余寿命预测方法

Info

Publication number: CN107451392A
Application number: CN201710483434.5A
Authority: CN
Inventors: 周东华; 席霄鹏; 陈茂银; 卢晓
Original assignee: Shandong University of Science and Technology
Current assignee: Shandong University of Science and Technology
Priority date: 2017-06-23
Filing date: 2017-06-23
Publication date: 2017-12-08
Anticipated expiration: 2037-06-23
Also published as: CN107451392B

Abstract

本发明公开了一种含有多个相关退化过程的剩余寿命预测方法，属于预测与健康管理中剩余寿命预测领域，包括以下步骤：读入间接的状态监测退化数据，即传感器测量数据，初始化状态空间模型参数；利用序贯卡尔曼滤波方法辨识隐含的退化状态，并在此基础上通过EM算法迭代更新模型的未知参数；基于蒙特卡洛方法对所得退化过程进行外推，在三种不同的系统失效模式下，估计剩余寿命的可靠度函数，并通过数值微分求取对应的概率密度函数；最后，利用均方误差分析模型的拟合效果以及剩余寿命的预测精度。本方法较以往方法具有更强的普适性。

Description

一种含有多个相关退化过程的剩余寿命预测方法

技术领域

本发明属于预测与健康管理中剩余寿命预测领域，具体涉及一种含有多个相关退化过程的剩余寿命预测方法。

背景技术

基于状态监测信息的剩余寿命预测是预测与健康管理领域中的一项重要技术。该技术是将退化过程首次到达预设的失效阈值的时间看作随机过程的首达时问题，其核心思想在于计算首达时间的期望或是概率密度函数，从而为预测维护工作提供合理的先验知识，能够有效地降低维护成本。考虑到预测的不确定性，恰当的概率密度函数更受青睐。

目前，预测与健康管理领域国内外学者提出的剩余寿命预测算法主要包括四大类，分别是物理机理方法，实验方法，数据驱动方法以及复合方法。其中，数据驱动方法只依赖于观测数据，在实际过程中使用较为广泛。对于此类方法，按照信息的获取方式又可分为基于直接监测数据和基于间接监测数据两类。注意到绝大多数工作仅仅针对的是一维退化数据，却不考虑系统的复杂度。也就是说，系统的全部性能指标假定只由单一退化过程进行描述。

但是，实际工业系统通常包含多种不同类型的操作单元，不同单元之间还可能存在一定的相关关系。举例来说，大型高炉由炉喉，炉身，炉腹，炉腰和炉缸构成，其中，炉缸的不断侵蚀会对风口造成很大程度的影响，甚至烧穿。对于这种情况，一个关键的问题在于如何预测含有多个相关退化过程的系统的剩余寿命。

近几年利用Copula函数建模退化过程相关性的方法取得了较为广泛的应用。但是，直接使用Copula函数存在一些明显的局限性：第一，Copula函数只能引入统计意义下的相关系数，并没有封闭的数学形式；第二，利用Copula函数降低维度会丢失大量有用信息，比如协方差矩阵中的交叉项；第三，如何选择合适的Copula函数尚且是一个艰难的问题。

显然，已有的基于Copula函数的退化建模不能较好地描述多个相关的退化过程。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种含有多个相关退化过程的剩余寿命预测方法，设计合理，克服了现有技术的不足，具有良好的效果。

为了实现上述目的，本发明采用如下技术方案：

一种含有多个相关退化过程的剩余寿命预测方法，该方法在处理退化数据时按照以下步骤依次实现：

步骤1：读入M组状态监测退化数据，即传感器测量数据，第j组传感器测量数据记为其中，j＝1,2,…,M；

步骤2：将每组首个数据作为初始测量偏差各个监测时刻下隐含的退化状态记为x_k，针对如下结构初始化模型参数，包括漂移系数λ，非线性系数γ，状态方程的协方差矩阵Q，测量系数矩阵g^(j)以及各个测量的噪声方差

其中，表示第j组传感器在t_k时刻下的测量向量，Φ_k～N(0,Q_k)，τ_k为第k个监测间隔，当Φ_k为对角矩阵时，各个退化过程之间互不相关，模型退化为多个传统的一维模型，而当扩散系数矩阵为一般矩阵时，第i个退化过程x_i(t)不仅仅由对应的第i个扩散项B_i(t)决定，还可能受到B₁(t),…,B_i-1(t),B_i+1(t),…,B_N(t)的影响；

步骤3：基于序贯卡尔曼滤波方法辨识隐含的退化状态x_k|k；

步骤4：利用EM算法迭代更新未知参数λ,γ,Q,的估计；

步骤5：通过蒙特卡洛仿真方法，分别在串联、并联以及加性的失效模式下(记作M₁,M₂,M₃)，计算剩余寿命的可靠度函数R_k(t)，即

其中，L_j表示第j个退化过程x_j(t)的失效阈值，且有j＝1,2,…,N；τ表示蒙特卡洛采样间隔，为第k个监测时刻下退化状态的估计值，Υ_k为第k个监测时刻以前的所有测量数据；

步骤6：利用数值微分求取剩余寿命的概率密度函数，即

步骤7：分别计算三种失效模式下的均方误差MSE，即

其中，为t_k时刻时真实的剩余寿命，根据MSE的值检验模型的预测效果：若MSE随着监测时刻逐渐减小直至趋于零，即说明预测效果良好。按照使用对象的实际需求判断系统的失效模式，一般包括串联失效，并联失效以及加性失效，具体地，若系统各个部件均为关键部件，则为串联失效；若尚有一个部件正常仍能继续工作，则为并联失效；若各个部件对于系统的失效体现为部分失效，且存在累积效应，则为加性失效，最终输出对应的f_k(t)。

本发明所带来的有益技术效果：

本方法能够处理含多个相关隐退化过程的复杂系统，剩余寿命的预测精度大大提高；本方法相对以往剩余寿命预测算法普适性更强，准确度更高；本方法中新的退化模型将单一退化过程升级为多个退化过程，考虑了隐含的相关性，体现在扩散系数矩阵当中，具体为：当扩散系数矩阵为对角矩阵时，各个退化过程之间互不相关，模型退化为多个传统的一维模型；当扩散系数矩阵为一般矩阵时，第i个退化过程x_i(t)不仅仅由对应的第i个扩散项B_i(t)决定，还可能受到B₁(t),…,B_i-1(t),B_i+1(t),…,B_N(t)的影响；对于多个相关的隐退化过程，本方法仍能较为准确的预测系统的剩余寿命。

本发明方法可以应用到以下场合：

1、高速列车，轮船，航天器等的性能退化过程；

2、高炉炼铁等化工生产过程。

附图说明

图1是本发明方法的流程图。

图2是实例一的退化轨迹。

图3是实例一的参数估计结果。

图4是实例一的剩余寿命预测结果。

图5是实例一的均方误差。

图6是实例二的观测数据。

图7是实例二的剩余寿命预测结果。

图8是实例二的均方误差。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

1、一种含有多个相关隐退化过程情况下的剩余寿命预测算法，其流程如图1所示，具体按照以下步骤依次实现：

步骤1：读入M组状态监测退化数据，即传感器测量数据，分别记为y⁽¹⁾(t),y⁽²⁾(t),…,y^(M)(t)，其中，

步骤2：取每组首个数据作为初始测量偏差记各个监测时刻下隐含的退化状态为x_k，针对如下结构初始化模型参数

其中，表示第j组传感器在t_k时刻下的测量向量，Φ_k～N(0,Q_k)，τ_k为第k个监测间隔，

步骤3：基于序贯卡尔曼滤波方法辨识隐含的退化状态x_k|k。首先，针对各组传感器，计算时间更新，即

其中，和分别为第j组传感器在t_k时刻下的退化状态和协方差矩阵的时间更新，λ_k和γ_k分别为t_k时刻下退化状态的漂移系数和非线性系数，Q_k为协方差矩阵的增量。

然后，计算测量更新，即

其中，和分别为第M组传感器在t_k时刻下的退化状态和协方差矩阵测量更新，和P_k+1|k+1分别为t_k时刻下的退化状态和协方差矩阵最终的测量更新，且有

步骤4：利用EM算法迭代更新未知参数λ,γ,Q,的估计。首先是E-步，为保证状态估计信息的可靠性，计算隐退化状态的后向迭代，即

其中，为t_k时刻下后向估计的t_d时刻的退化状态，P_d|k为t_k时刻下后向估计的t_d时刻的协方差矩阵，P_d,d-1|k为t_k时刻下后向估计的t_d时刻和t_d-1时刻的互协方差矩阵，且有

然后，在M-步中，利用数值方法(如MATLAB中的fsolve函数)求解γ_k的估计，即

计算λ_k,Q_k,的估计，即

检验参数估计的收敛性，即

其中，Δ为一个小的正数；

步骤5：初始化蒙特卡洛采样间隔以及离散时间序号l，基于估计的未知参数集通过蒙特卡洛方法进行离散采样，生成大量预测样本即

其中，

分别在串联、并联以及加性的失效模式(记作M₁,M₂,M₃)下，计算剩余寿命的可靠度函数R_k(t)，即

步骤6：当l足够大，基于可靠度数据

利用数值微分求取剩余寿命的概率密度函数，即

步骤7：分别计算三种失效模式下的均方误差MSE，即

其中，为t_k时刻下真实的剩余寿命，根据MSE的值检验模型的预测效果：若MSE随着监测时刻逐渐减小直至趋于零，即说明预测效果良好。按照使用对象的实际需求判断系统的失效模式，一般包括串联失效，并联失效以及加性失效，具体地，若系统各个部件均为关键部件，则为串联失效；若尚有一个部件正常仍能继续工作，则为并联失效；若各个部件对于系统的失效体现为部分失效，且存在累积效应，则为加性失效。并最终输出对应的f_k(t)。

2、实例寿命预测结果

下面结合实例对该发明的具体实施步骤做进一步的具体说明。

下述实例的仿真环境如下：

机型：Intel Core i7-4790(CPU 3.60Ghz，8.00GB RAM)；

操作系统：Windows 10；

软件：Matlab R2016a。

a、实例一：通过一组数值仿真，给出实验步骤的详细说明以及剩余寿命预测结果：

1)考虑二元退化的情况，假定每个退化过程由两个采样速率相同的传感器单独进行监测，即N＝2，M＝2，根据模型(1)生成退化数据，其中，参数设置为τ＝0.1，γ＝0.8，λ＝[0.5；0.6]，Σ＝[0.1,0；0,0.2]，V＝[0.4；0.4]，g⁽¹⁾＝[1.2,0；0,1.1]，g⁽²⁾＝[1,0；0,1]，Ξ⁽¹⁾＝[0.01,0；0,0.04]，Ξ⁽²⁾＝[0.02,0；0,0.03]，两个退化过程的失效阈值分别为L₁＝5.2748和L₂＝6.7175，首达阈值的时间约为第200个监测时刻；

2)取每组首个数据作为初始测量偏差记各个监测时刻下隐含的退化状态为x_k，初始化模型参数

3)利用式(2)-式(8)辨识隐含的退化状态x_k|k，估计的退化轨迹以及真实的退化轨迹如图2所示，可以看出，估计结果能够很好地跟上真实的退化过程，说明了模型的有效性；

4)利用式(9)-式(16)迭代更新未知参数γ_k,λ_k,Q_k,的估计，各个参数的估计结果如图3所示，可以看出，随着观测数据的积累，参数的估计值能够较快地收敛到各自的真值，多源信息对于收敛速度的提升起到了很大的帮助；

5)初始化蒙特卡洛采样间隔以及离散时间序号l，基于参数的估计结果，利用式(17)生成大量预测样本分别在串联，并联以及加性的失效模式下，利用式(18)计算剩余寿命的可靠度函数R_k(t)；

6)利用式(19)求取剩余寿命的概率密度函数f_k(t)，以第170个监测时刻为例，三种失效模式下的剩余寿命预测结果如图4所示，可以看出，真实的剩余寿命较好地落在预测的概率密度函数的下方，此外，加性失效模式下的预测结果相对最为准确，而其它两种模式下的真实值略微偏离预测的概率密度函数的峰值，这是由于加性模式的失效条件相对其它两种模式做出了正确的调整；

7)利用式(20)计算各个失效模式下的MSE，结果如图5所示，可以看出，MSE会随着监测时刻逐渐减小，说明在线预测的概率密度函数较为鲁棒。

b、实例二：以图6为实例，进行实验步骤的详细说明，并给出剩余寿命预测结果：

1)利用C-MAPSS模拟大型商用涡扇发动机的退化，为简便起见，考虑二元退化过程(包含HPC退化和风扇退化)，将传感器3和传感器7分至第一组，同时将传感器4和传感器12分至第二组，以第一个测试单元为例，对应的退化数据如图6所示，数据长度为259，失效阈值分别为L₁＝208.4670和L₂＝-187.7002；

3)利用式(2)-式(8)辨识隐含的退化状态x_k|k；

4)利用式(9)-式(16)迭代更新未知参数γ_k,λ_k,Q_k,的估计；

5)初始化蒙特卡洛采样间隔以及离散时间序号l，基于参数的估计结果，利用式(17)生成大量预测样本以并联失效模式为例，利用式(18)计算剩余寿命的可靠度函数R_k(t)；

6)利用式(19)求取剩余寿命的概率密度函数f_k(t)，在线的剩余寿命预测结果如图7所示，通过和真值进行对比，可以看出，预测结果较为精确；

7)利用式(20)计算各个监测时刻下的MSE，结果如图8所示，可以看出，MSE急剧下降，并最终趋于零，从而验证了所提方法的有效性。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims

1.一种含有多个相关退化过程的剩余寿命预测方法，其特征在于：按照如下步骤依次进行：

步骤1：读入M组状态监测退化数据，即传感器测量数据，第j组传感器测量数据记为其中，j＝1，2，…，M；

<mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>x</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>&lambda;</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <msubsup> <mi>t</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <msub> <mi>&gamma;</mi> <mi>k</mi> </msub> </msubsup> <mo>-</mo> <msubsup> <mi>t</mi> <mi>k</mi> <msub> <mi>&gamma;</mi> <mi>k</mi> </msub> </msubsup> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&Phi;</mi> <mi>k</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>y</mi> <mi>k</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>y</mi> <mn>0</mn> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </msubsup> <mo>+</mo> <msup> <mi>g</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </msup> <msub> <mi>x</mi> <mi>k</mi> </msub> <mo>+</mo> <msubsup> <mi>&psi;</mi> <mi>k</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> </msubsup> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

其中，表示第j组传感器在t_k时刻下的测量向量，Φ_k～N(0,Q_k)，Φ_k为扩散系数矩阵；τ_k为第k个监测间隔，

步骤3：基于序贯卡尔曼滤波方法辨识隐含的退化状态x_k|k；

步骤4：利用EM算法迭代更新未知参数λ,γ,的估计；

步骤5：通过蒙特卡洛仿真方法进行离散采样，分别在串联失效模式M₁、并联失效模式M₂以及加性失效模式M₃下，计算剩余寿命的可靠度函数R_k(t)，即

其中，L_j表示第j个退化过程x_j(t)的失效阈值，且有j＝1,2,…,N；τ表示蒙特卡洛采样间隔；为第k个监测时刻下退化状态的估计值；Υ_k为第k个监测时刻以前的所有测量数据；

步骤6：利用数值微分求取剩余寿命的概率密度函数，即

步骤7：分别计算三种失效模式下的均方误差MSE，即

<mrow> <mi>M</mi> <mi>S</mi> <mi>E</mi> <mo>=</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msubsup> <mo>&Integral;</mo> <mn>0</mn> <mi>&infin;</mi> </msubsup> <msup> <mrow> <mo>(</mo> <mi>t</mi> <mo>-</mo> <msub> <mover> <mi>t</mi> <mo>~</mo> </mover> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <msub> <mi>f</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>t</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> <mo>;</mo> </mrow>

其中，为t_k时刻时真实的剩余寿命，根据MSE的值检验模型的预测效果，并输出对应的f_k(t)。

2.根据权利要求1所述的含有多个相关退化过程的剩余寿命预测方法，其特征在于：在步骤7中，根据MSE的值检验模型的预测效果的具体方法：若MSE随着监测时刻逐渐减小直至趋于零，即说明预测效果良好，按照使用对象的实际需求判断系统的失效模式，具体为若系统各个部件均为关键部件，则为串联失效；若尚有一个部件正常仍能继续工作，则为并联失效；若各个部件对于系统的失效体现为部分失效，且存在累积效应，则为加性失效。