CN110851980B - 一种设备剩余寿命预测方法及系统 - Google Patents

Abstract

本发明公开一种设备剩余寿命预测方法及系统，该方法包括：建立基于非线性扩散过程的设备退化数学模型；获取加速应力下设备退化参数的估计值，为第一数据；根据第一数据计算正常工作应力下的退化参数值，得到第二数据；通过拟合优度检验确定第二数据的分布类型；根据第二数据的分布类型，得到第二数据的后验分布函数和期望值；根据后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；根据全概率公式修正第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；根据第二剩余寿命概率密度函数预测设备剩余寿命。本发明能够提高设备剩余寿命的预测精度。

Description

一种设备剩余寿命预测方法及系统

技术领域

本发明涉及可靠性工程技术领域，特别是涉及一种设备剩余寿命预测方法及系统。

背景技术

近年来，在航空航天和军工等领域出现了越来越多长寿命、高可靠性的设备，对这些设备进行准确定寿是保持设备性能和做好预防性维修工作重要前提。对于这类高可靠性设备来说，即便是在加速寿命试验中也难以发生失效，传统的基于失效数据的方法无法有效预测此类设备的寿命。某些产品的部分性能指标会随试验时间发生退化，对退化过程中的性能退化量进行测量并统计分析，无需等待产品失效便可预测产品的寿命信息，并且当提高产品的所处环境的某种应力水平时，会加速设备的退化失效过程。因此，利用超过正常应力水平来加速设备性能退化，采集退化数据并对正常应力水平下工作的设备进行寿命预测的加速退化试验方法被广泛应用于设备的寿命预测和可靠性评估中。

加速退化试验通过抽样的方法对设备整体的寿命信息进行估计，无法准确描述个体的寿命状态。对于设备个体而言，在测试或贮存状态中积累的性能退化数据，往往因为数据量较少而无法对个体寿命信息进行准确预测。针对这类问题，基于Bayesian理论的剩余寿命预测方法得到学者们的深入研究。这类方法通常以批量设备在退化试验中的数据作为先验信息，以个体设备工作或定期测试中的数据作为样本信息，利用Bayesian方法得到个体设备剩余寿命预测模型中参数的后验信息，对其剩余寿命进行准确、实时的预测。目前，缺少将加速退化数据作为先验信息进行，利用状态监测数据进行统计推断的相关研究。

发明内容

本发明的目的是提供一种设备剩余寿命预测方法及系统，能够提高对设备剩余寿命预测的精度。

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

一种设备剩余寿命预测方法，包括：

建立基于非线性扩散过程的设备退化数学模型；

获取加速应力下设备退化参数的估计值，为第一数据；

根据所述第一数据计算正常工作应力下的退化参数值，得到第二数据；

通过拟合优度检验确定所述第二数据的分布类型；

根据所述第二数据的分布类型，得到第二数据的后验分布函数和期望值；

根据所述后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；

根据全概率公式修正所述第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；

根据所述第二剩余寿命概率密度函数预测设备剩余寿命。

可选的，所述设备退化数学模型为：

其中，a表征不同设备间的差异性，b表征同类设备的共性特征，σ_B为扩散系数，ω为设备的失效阈值。

可选的，所述获取加速应力下设备退化参数的估计值，包括：

根据维纳过程建立设备性能退化增量与时间增量的极大似然函数

其中，X_ijk为第j个样品在第k个加速应力下的第i次测量值，t_ijk为第j个样品在第k个加速应力下的第i次的测量时刻，其中i＝1,2,……,n₁，j＝1,2,……,n₂，k＝1,2,……,n₃，ΔX_ijk＝X_ijk-X_(i-1)jk为性能退化量的增量，Δt_ijk＝t_ijk-t_(i-1)jk为时间增量，参数a_jk和

均表示第j个样品在第k个加速应力下的参数值；

将b作为已知量求得参数a_jk和

关于b的解析表达式；

将所述解析表达式代入所述极大似然函数，得到估计值

将所述估计值

代入所述a_jk和

关于b的解析表达式，得到加速应力下设备退化参数的估计值

和

可选的，所述根据所述第一数据计算正常工作应力下的退化参数值，包括：

根据Nelson假设，确定加速因子；

获取阿伦尼斯模型的线性表达式；

根据所述加速应力下设备退化参数的估计值、加速因子和线性表达式得到正常工作应力下的退化参数值。

可选的，所述第一剩余寿命概率密度函数为：

其中，R(t_k)为t_k时刻的可靠度函数，

可选的，所述第二剩余寿命概率密度函数为：

一种设备剩余寿命预测系统，包括：

模型建立模块，用于建立基于非线性扩散过程的设备退化数学模型；

参数获取模块，用于获取加速应力下设备退化参数的估计值，为第一数据；

计算模块，用于根据所述第一数据计算正常工作应力下的退化参数值，得到第二数据；

拟合优度检验模块，用于通过拟合优度检验确定所述第二数据的分布类型；

后验分布确定模块，用于根据所述第二数据的分布类型，得到第二数据的后验分布函数和期望值；

第一剩余寿命概率密度函数确定模块，用于根据所述后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；

第二剩余寿命概率密度函数确定模块，用于根据全概率公式修正所述第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；

设备剩余寿命预测模块，用于根据所述第二剩余寿命概率密度函数预测设备剩余寿命。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明通过建立基于非线性扩散过程的设备退化过程的模型，并建立模型参数与加速应力之间的关系，利用加速退化数据确定模型参数的先验分布类型和相关超参数，利用状态监测数据对模型参数的后验分布进行更新，并采用基于吉布斯采样的马尔科夫链蒙特卡洛方法进行数值求解，最后，基于首达时间的概念，通过马尔科夫链蒙特卡洛方法获得考虑模型参数随机性的剩余寿命分布，能够提高设备剩余寿命的预测精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明设备剩余寿命预测方法流程图；

图2为本发明设备剩余寿命预测系统模块图；

图3为本发明实施例各加速应力下电连接器接触电阻退化曲线图；

图4为本发明实施例工作应力下电连接器接触电阻退化曲线图；

图5为本发明实施例电连接器剩余寿命概率密度函数示意图；

图6为本发明实施例电连接器剩余寿命的均方误差示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种设备剩余寿命预测方法及系统，能够提高对设备剩余寿命预测的精度。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明设备剩余寿命预测方法流程图，如图1所示，一种设备剩余寿命预测方法，包括：

步骤101：建立基于非线性扩散过程的设备退化数学模型；

步骤102：获取加速应力下设备退化参数的估计值，为第一数据；

步骤103：根据所述第一数据计算正常工作应力下的退化参数值，得到第二数据；

步骤104：通过拟合优度检验确定所述第二数据的分布类型；

步骤105：根据所述第二数据的分布类型，得到第二数据的后验分布函数和期望值；

步骤106：根据所述后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；

步骤107：根据全概率公式修正所述第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；

步骤108：根据所述第二剩余寿命概率密度函数预测设备剩余寿命。

其中，步骤101具体包括：

设X(t)表示样品设备在t时刻的性能退化量，则基于扩散过程的退化过程{X(t),t≥0}可以表示为：

其中，X(0)表示初始时刻的性能退化量，通常X(0)＝0，μ(t；θ)为漂移系数，是时间t的非线性函数，用以表示模型的非线性特征，其中参数向量θ＝(a,b)为未知参数。这里a表征不同设备间的差异性，b表征同类设备的共性特征，σ_B称为扩散系数，B(t)为标准Brownian运动，且B(t):N(0,t)。

对于不同函数形式的μ(t；θ)可以描述不同形式的非线性随机退化过程。显然，若μ(t；θ)＝μ，则式(1)则转化为线性随机退化模型，即Wiener过程。因此，线性随机退化模型是本文所讨论的非线性随机退化模型的特例。不失一般性，本文主要考虑X(0)＝0的情况。

对于(1)式给出的随机退化过程，在首达时间意义下，设备的寿命T可定义为：

T＝inf{t:X(t)≥ω|X(0)＜ω} (2)

其中，寿命T为随机变量，ω为设备的失效阈值，通常由工业标准和专家经验确定。

对于式(1)出的随机退化过程，当不考虑参数a和σ_B的随机性时，设备寿命的概率密度函数可以近似为：

其中，

其中，S(t)是一个表达式，无具体含义。

幂函数和指数函数这两种典型形式在描述金属疲劳和机械设备的磨损退化数据方法应用较为广泛，通常考虑μ(t:θ)为幂函数形式和指数函数形式。本发明以μ(t；θ)＝abτ^b-1为例进行说明，当b＝1时，式(1)描述的是Wiener过程，当b≠1时，式(1)描述的是应用更广泛的扩散过程。将μ(t；θ)＝abτ^b-1代入式(3)，可化简为：

其中，a表征不同设备间的差异性，b表征同类设备的共性特征，σ_B为扩散系数，ω为设备的失效阈值，通常由工业标准和专家经验确定。

步骤102具体包括：

假设对样品设备进行恒定应力加速退化试验，S₀为工作应力水平，S_k为第k个加速应力水平，X_ijk为第j个样品设备在第k个加速应力下的第i次测量值，t_ijk为对应的测量时刻，其中i＝1,2,……,n₁；j＝1,2,……,n₂；k＝1,2,……,n₃。ΔX_ijk＝X_ijk-X_(i-1)jk为性能退化量的增量，Δt_ijk＝t_ijk-t_(i-1)jk为时间增量。根据Wiener过程的特性，ΔX_ijk服从如下的正态分布：ΔX_ijk：

即ΔX_ijk：

可根据每个样品设备n₁次测量数据(ΔX_ijk,Δt_ijk)建立如下形式的极大似然函数：

其中，参数a_jk和

均表示第j个样品在第k个加速应力下的参数值。

然后采用两步法对参数a_jk,

b进行估计：首先将参数b视为已知量，分别求得不同应力下不同样品设备的参数a_jk和

关于b的解析表达式。然后将所有的解析式代入极大似然函数，并利用加速退化数据可以得到估计值

将估计值

代入a_jk和

关于b的解析表达式，即可得到各应力下各样品设备的参数估计值

和

例如：以第k个加速应力下，第j个样品设备为例，共有n₁个加速退化数据。对式(5)两端同时取对数，可得：

分别求取式(6)于a_jk和

的一阶偏导数，可得：

令式(7)为0，可解得：

令式(8)为0，并将式(9)代入，可得：

共可得到n₂×n₃个a_jk和

关于b的解析表达式，将其全部代入式(11)，利用加速退化数据可以得到估计值

将

分别代入式(9)、(10)，可得到各应力下每个样品设备基于扩散过程的漂移参数和扩散参数平方的估计值

即第一数据。

步骤103具体包括：

根据Nelson假设，给出了加速因子的一种定义：假设F_p(t_p)和F_q(t_q)分别表示产品在应力S_p和S_q下的累积失效概率，若存在F_p(t_p)＝F_q(t_q)，则可以将应力S_p相对于S_q的加速因子定义为A_pq＝t_q/t_p。

样品设备在各加速应力下的失效机理不变是保证加速因子A_pq与可靠度(时间)无关的充要条件，即：

F_p(t_p)＝F_q(A_pqt_p)(12)

对式(12)两端同时对时间t求导，可得：

将式(4)代入式(13)，化简后可得：

由于b表征同类设备的共性特征，故在不改变设备失效机理的情况下，可认为b_p＝b_q，记为b。则式(14)可化简为：

为求得加速因子A分别与a和

的关系，需保证A_pq与时间t无关，则可得到

即

常用的加速模型有阿伦尼斯模型、逆幂律模型和艾林模型。在以温度应力作为影响设备性能退化的主要应力时，通常采用阿伦尼斯模型来描述性能退化参数与温度应力的关系，其一般表达式为：

η＝Mexp[E_a/kT](17)

其中，η为设备在温度加速应力水平作用下的性能退化参数；T为温度加速应力水平；m是与失效模式、加速试验类型以及其他因素相关的常数；E_a为激活能，其具体大小与发生失效模式的材料有关；k为Boltzmann常数，k＝8.6171×10^-5eV/K。

将加速模型式(17)进行线性化处理可表示为：

其中，m＝lnM；h＝-E_a/k；

为加速应力函数。m和h均为待估计常数。

将a和

分别代入式(18)可得：

可根据不同加速应力下的估计值

并利用最小二乘法得到m_a、h_a、

和

的估计值，其中，m_a、h_a、

和

分别为m和h关于参数a和

的估计值。

由加速模型式(18)可以得到A_pq的另一种表达式为：

将h_a和

分别代入式(21)，可进一步可计算出各加速应力S_k相对于工作应力S₀的加速因子A_k0，即可将S_k下的参数值

折算到正常工作应力S₀下，即

即为第二数据；其中，正常工作应力是设备在正常工作状态下所受的环境应力。比如，电连接器要求正常工作时，环境温度保持20℃，则20℃即为该电连接器的正常工作应力。

步骤104具体包括：

为简化后续表达式，定义D₀＝1/σ²B₀。由于漂移参数a₀和扩散参数D₀不服从共轭先验分布，那么首先要确定参数a₀和D₀的分布类型，进而才能对超参数进行估计和Bayesian更新。在得到折算后工作应力下a₀和D₀的折算值后，即可使用拟合优度检验的方法确定a₀和D₀的最优拟合分布类型。

由于Anderson-Darling拟合优度检验方法具有良好的统计特性，故选用此方法来确定参数值的分布类型。Anderson-Darling统计量可描述数据服从特定分布类型的程度，数据与分布拟合的越好，Anderson-Darling统计量的AD值越小。其计算公式为：

其中，F_n(x)为经验分布函数，F(x)为累积分布函数，n为数据个数。

为了便于超参数的先验估计和后续计算，令

设a_i:f_a(θ_a)，D_i:f_D(θ_D)，π(a,D)为a和D的联合概率密度函数，其中θ_a和θ_D分别为参数a和D的超参数向量。假定f_a(θ_a)与f_D(θ_D)不相关，则有π(a,D)＝f_a(θ_a)·f_D(θ_D)。此时，对参数分布的确定就是依据加速退化数据确定f_a(θ_a)和f_D(θ_D)。

其中，参数a和D是为了说明二者之间的关系，后面出现的a₀和D₀表示的是经过折算后，在工作应力下的参数值。

假设通过AD拟合优度检验得到

D₀～Ga(α,β)，其中，超参数μ_a和

分别为正态分布的均值和方差，超参数α和β分别为伽玛分布的形状参数和尺度参数。

步骤105具体包括：

根据第二数据的各自分布的概率密度函数建立极大似然方程对超参数值进行估计。

因为参数a和D具体的分布类型需要根据估计出的数据才能确定，这里为了便于方法的叙述，先假设

D₀～Ga(α,β)

，即a₀服从正态分布，D₀服从Gamma分布。

已知

D₀～Ga(α,β)，则可得到a₀和D₀的联合先验密度函数π(a₀,D₀)为：

利用Bayesian公式，可推导出参数a₀和D₀的联合后验密度函数π(a*,D*|X)为：

其中，π(a₀,D₀)为a₀和D₀的联合先验密度函数，a*和D*分别为a₀和D₀的后验参数，L(X|a₀,D₀)为极大似然函数。此时极大似然函数L(X|a₀,D₀)和联合先验密度函数π(a₀,D₀)分别为：

其中，X＝[X(t₁),X(t₂),……,X(t_m)]为m个现场退化数据。

由于直接求取参数的联合后验密度函数π(a*,D*|X)较为困难，难以得到解析解，所以考虑采用马尔科夫链蒙特卡洛(Markov Chain Monte Carlo,MCMC)方法进行求解。采用MCMC方法求取参数的后验分布可以选择多种抽样方法，这里给出Gibbs抽样方法的一般步骤。

1)分别给出参数a₀和D₀的初始值

和

2)假设第k次迭代开始时参数分别为

和

则按照下面的子步骤得到第k次迭代后的参数值

和

(i)由式(28)，即满条件分布

抽取

(ii)由式(29)，即满条件分布

抽取

在实际问题中，可以采用上述抽样方法得到Monte Carlo样本，进而求得参数的后验分布，也可以借助于WinBUGS软件、Matlab MCMC工具箱等方法进行计算。

已知利用工作应力下的退化数据进行Bayesian更新后的参数a*和D*的联合后验概率密度函数π(a*,D*|X)，则参数a*和D*后验分布的边缘密度函数为：

进一步可得后验参数的期望值为：

步骤106具体包括：

将式(30)、(31)得到的结果

和

代入式(4)，即可得到产品寿命的概率密度函数

鉴于此，根据剩余寿命L_k的定义式L_k＝{l_k:T-t_k|T＞t_k}，任意时刻t_k的第一剩余寿命概率密度函数可以表示为：

其中，R(t_k)为t_k时刻的可靠度函数，

注意到，此时由式(34)得到的产品剩余寿命的概率密度函数，仅仅利用了后验参数的期望值，并未考虑到参数的随机性。

步骤107具体包括：

根据式(4)和全概率公式可以得到在考虑参数a₀和D₀随机性的情况下，产品寿命的概率密度函数为：

由于式(35)难以获得解析表达式，并且其结果依赖于参数的分布形式，不具有普适性。本文利用MCMC方法，充分考虑参数a₀和D₀的随机性，以得到随机参数作用下的产品寿命概率密度函数。

利用参数a₀和D₀的联合后验概率密度函数π(a*,D*|X)，以此作为平稳分布建立Markov链来得到π(a*,D*|X)的样本(a₀,D₀)⁽ⁱ⁾，基于这些样本，即可进行统计推断。若得到样本(a₀,D₀)⁽¹⁾,……,(a₀,D₀)⁽ⁿ⁾，易知其相互独立，则根据大数定理可以得到式(35)的估计值：

故利用式(36)可以得到考虑参数a₀和D₀随机作用下的产品寿命概率密度函数，将其代入式(34)，即可得到产品的第二剩余寿命概率密度函数

根据第二剩余寿命概率密度函数能够预测设备的剩余寿命。

实施例：

以某型军用电连接器为例，验证本发明所提方法的有效性，验证对象包括本发明所提出的基于扩散过程的融合加速退化数据和实际退化数据的设备剩余寿命预测方法(记为M₁)、基于Wiener过程的融合加速退化数据和实际退化数据的剩余寿命预测方法(记为M₂)和基于对数变换数据lnX(t)的剩余寿命预测方法(记为M₃)。

样本的正常工作应力S₀＝20℃，随机抽取24个样品，分别在三种加速应力S₁＝60℃，S₂＝80℃，S₃＝100℃下进行加速退化试验，各应力下分配8个样品。假设在应力S₁下以时间间隔48小时进行30次测量；在应力S₂下以时间间隔36小时进行25次测量；在应力S₃下以时间间隔24小时进行20次测量。根据该型电连接器的国家标准，样品接触阻值的失效阈值ω＝5mΩ，各应力下样品的退化数据如图3所示。

首先利用Anderson-Darling统计量以置信度90％对ΔX_ijk进行假设检验，证明样品的退化过程服从Wiener过程。将加速退化数据(ΔX_ijk,Δt_ijk)代入式，并利用两步法对参数

b进行估计，得到

以及不同应力下每个样品的

的估计值，如表1所示。

表1 不同加速应力下各样品

估计值

利用最小二乘法和

的估计值，对式(19)、(20)中的未知参数h进行估计，得到

由此可得到各加速应力折算到工作应力下的加速因子值，如表2所示。

表2 加速因子值

利用表2中的A_i0的值，可将各加速应力下的参数估计值折算到工作应力下，结果如表3所示。

表3

折算到工作应力下的参数值

利用Anderson-Darling统计量分别对a₀和D₀可能的分布类型进行最优拟合优度检验。选择常见的指数分布、极值分布、正态分布、对数正态分布、伽玛分布和威布尔分布作为可能的分布类型，具体结果如表4所示。

表4 a₀和D₀在各分布类型下的AD值

Anderson-Darling统计量的AD值越小，说明待拟合数据与该分布类型拟合得越好，由表4可知，a₀最优拟合于正态分布，D₀最优拟合于伽玛分布。依据式和式，可以得到超参数μ_a、

α和β的极大似然估计值。由此可确定a₀和D₀的先验分布为a₀～N(3.272×10^-6,0.405×10^-6)，D₀～Ga(6,3.0279×10⁴)。

对某一正常工作应力使用中的同型号电连接器每隔1500小时进行一次接触电阻值测量，共获得10个退化数据，如图4所示。由图4可知，在进行第8次测量后，退化轨迹首次达到失效阈值。因时，可认为该样品的实际寿命约为L＝1.26×10⁴小时。

在获取到第i(1≤i≤8)次测量值后，对产品进行剩余寿命预测。三种方法的预测结果如表5所示。

表5 三种方法进行剩余寿命预测的结果

从表5中可以看出，在剩余寿命预测的初始阶段，三种方法均存在一定程度的偏差。利用不断得到的工作应力下的退化数据对漂移系数和扩散系数进行更新，M₁方法得到的剩余寿命预测值越来越接近剩余寿命的真实值，其相对误差逐渐减小，在第6次参数更新后，稳定在2.5％左右。M₂方法不考虑退化数据的非线性，以Wiener过程为模型进行剩余寿命寿命预测，随着数据的非线性特征的凸显，预测值的相对误差明显增大。M₃方法利用对数变换将非线性的加速退化数据进行了变换，再利用线性的Wiener过程进行剩余寿命预测。由于M₃方法将数据的非线性考虑在内，其预测值的相对误差较M₂方法有所改善，但相比于M₁方法则有较大差距。该结果说明，线性退化模型无法应用于非线性退化产品的寿命预测中，同时对数变换对该加速退化数据线性化的能力在这个实例中的作用是有限的，也验证了并不是所有的数据仅通过数据变换就能转化为线性数据。

为了直观的比较三种方法的区别，分别绘制出三种方法在各测量时刻的剩余寿命概率密度函数、平均剩余寿命和实际剩余寿命，如图5所示。

在图5中，实际剩余寿命用红色圆点标注，而剩余寿命预测的平均值用蓝色圆点标注。由图5可知，M₁方法得到的剩余寿命预测的概率密度函数和均值都要明显优于M₂和M₃方法的结果。此外，M₃方法的预测结果要优于M₂方法，但仍与实际结果相差较大。相比之下，M₁方法预测的剩余寿命的均值在后期与实际剩余寿命几乎重叠。验证结果表明，在退化数据具有非线性特征时，考虑非线性随机模型对数据进行建模并实现剩余寿命预测是非常必要的。通过式可以确定三种方法在各测量点处的均方误差，如图6所示。

从图6可以看出，利用M₁方法得到的剩余寿命的均方误差小于M₂和M₃方法，进一步验证了本发明提出的融合加速退化数据和实测退化数据的剩余寿命预测方法的有效性和优越性。

本发明还公开了一种设备剩余寿命预测系统，如图2所示，该系统包括：

模型建立模块201，用于建立基于非线性扩散过程的设备退化数学模型；

参数获取模块202，用于获取加速应力下设备退化参数的估计值，为第一数据；

计算模块203，用于根据所述第一数据计算正常工作应力下的退化参数值，得到第二数据；

拟合优度检验模块204，用于通过拟合优度检验确定所述第二数据的分布类型；

后验分布确定模块205，用于根据所述第二数据的分布类型，得到第二数据的后验分布函数和期望值；

第一剩余寿命概率密度函数确定模块206，用于根据所述后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；

第二剩余寿命概率密度函数确定模块207，用于根据全概率公式修正所述第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；

设备剩余寿命预测模块208，用于根据所述第二剩余寿命概率密度函数预测设备剩余寿命。

本发明还公开了如下技术效果：

本发明通过融合加速退化数据和状态监测数据，解决了高可靠性设备的剩余使用寿命的估计问题。首先，提出了一种基于非线性扩散过程的模型来表征设备的退化过程。然后，通过建立模型参数和加速应力之间的关系，利用加速退化数据确定模型参数的先验分布类型和相关超参数。为了实现加速退化数据和状态监测数据(正常工作应力下的退化数据)的融合，在新的状态监测数据可用时，采用贝叶斯方法对模型参数的后验分布进行更新，并采用基于吉布斯采样的马尔科夫链蒙特卡洛方法进行数值求解。最后，基于首达时间的概念，通过马尔科夫链蒙特卡洛方法获得考虑模型参数随机性的近似剩余寿命分布。本发明具有更高的剩余寿命估计精度，具有一定的工程实用价值。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (4)

1.一种设备剩余寿命预测方法，其特征在于，包括：

建立基于非线性扩散过程的设备退化数学模型；

获取加速应力下设备退化参数的估计值，为第一数据；

包括：

根据维纳过程建立设备性能退化增量与时间增量的极大似然函数

其中，X_ijk为第j个样品在第k个加速应力下的第i次测量值，t_ijk为第j个样品在第k个加速应力下的第i次的测量时刻，i＝1,2,……,n₁，j＝1,2,……,n₂，k＝1,2,……,n₃，ΔX_ijk＝X_ijk-X_(i-1)jk为性能退化量的增量，Δt_ijk＝t_ijk-t_(i-1)jk为时间增量；参数a_jk和

均表示第j个样品在第k个加速应力下的参数值；

将b作为已知量求得参数a_jk和

关于b的解析表达式；

将所述解析表达式代入所述极大似然函数，得到估计值

将所述估计值

代入所述a_jk和

关于b的解析表达式，得到加速应力下设备退化参数的估计值

和

根据所述第一数据计算正常工作应力下的退化参数值，得到第二数据；

通过拟合优度检验确定所述第二数据的分布类型；

根据所述第二数据的分布类型，得到第二数据的后验分布函数和期望值；

根据所述后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；

所述第一剩余寿命概率密度函数为：

其中，R(t_k)为t_k时刻的可靠度函数，

根据全概率公式修正所述第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；

所述第二剩余寿命概率密度函数为：

根据所述第二剩余寿命概率密度函数预测设备剩余寿命。

2.根据权利要求1所述的设备剩余寿命预测方法，其特征在于，所述设备退化数学模型为：

其中，a表征不同设备间的差异性，b表征同类设备的共性特征，σ_B为扩散系数，ω为设备的失效阈值。

3.根据权利要求1所述的设备剩余寿命预测方法，其特征在于，所述根据所述第一数据计算正常工作应力下的退化参数值，包括：

根据Nelson假设，确定加速因子；

获取阿伦尼斯模型的线性表达式；

根据所述加速应力下设备退化参数的估计值、加速因子和线性表达式得到正常工作应力下的退化参数值。

4.一种设备剩余寿命预测系统，其特征在于，包括：

模型建立模块，用于建立基于非线性扩散过程的设备退化数学模型；

参数获取模块，用于获取加速应力下设备退化参数的估计值，为第一数据；

包括：

根据维纳过程建立设备性能退化增量与时间增量的极大似然函数

其中，X_ijk为第j个样品在第k个加速应力下的第i次测量值，t_ijk为第j个样品在第k个加速应力下的第i次的测量时刻，i＝1,2,……,n₁，j＝1,2,……,n₂，k＝1,2,……,n₃，ΔX_ijk＝X_ijk-X_(i-1)jk为性能退化量的增量，Δt_ijk＝t_ijk-t_(i-1)jk为时间增量；参数a_jk和

均表示第j个样品在第k个加速应力下的参数值；

将b作为已知量求得参数a_jk和

关于b的解析表达式；

将所述解析表达式代入所述极大似然函数，得到估计值

将所述估计值

代入所述a_jk和

关于b的解析表达式，得到加速应力下设备退化参数的估计值

和

计算模块，用于根据所述第一数据计算正常工作应力下的退化参数值，得到第二数据；

拟合优度检验模块，用于通过拟合优度检验确定所述第二数据的分布类型；

后验分布确定模块，用于根据所述第二数据的分布类型，得到第二数据的后验分布函数和期望值；

第一剩余寿命概率密度函数确定模块，用于根据所述后验分布函数、期望值和设备退化数学模型得到第一剩余寿命概率密度函数；

所述第一剩余寿命概率密度函数为：

其中，R(t_k)为t_k时刻的可靠度函数，

第二剩余寿命概率密度函数确定模块，用于根据全概率公式修正所述第一剩余寿命概率密度函数，得到第二剩余寿命概率密度函数；

所述第二剩余寿命概率密度函数为：

设备剩余寿命预测模块，用于根据所述第二剩余寿命概率密度函数预测设备剩余寿命。

Images