CN112800616B - 基于比例加速退化建模的设备剩余寿命自适应预测方法 - Google Patents

Abstract

本发明公开了基于比例加速退化建模的设备剩余寿命自适应预测方法，涉及设备剩余寿命预测技术领域。针对加速退化试验条件下单台设备剩余寿命预测问题，首先基于非线性Wiener过程构建比例加速退化模型；其次，在退化模型的基础上建立状态转移方程，并采用卡尔曼滤波KF算法更新设备退化状态；再次，输入设备的性能退化量的观测数据，采用期望最大化‑卡尔曼滤波EM‑KF算法实现对退化模型中未知参数的自适应估计；最后，在退化状态更新和未知参数自适应估计的基础上，基于全概率公式计算出设备剩余寿命的概率密度函数和累积分布函数。利用本发明的方法，实现了单台设备加速退化试验下剩余寿命预测更加准确的效果。

Description

基于比例加速退化建模的设备剩余寿命自适应预测方法

技术领域

本发明涉及设备剩余寿命预测技术领域，特别是涉及基于比例加速退化建模的设备剩余寿命自适应预测方法。

背景技术

出于保证飞行安全和任务完成的需要，机载设备普遍具有高可靠、长寿命等特点，导致传统寿命试验和退化试验难以快速地获取足够的寿命/退化数据来确保剩余寿命预测的准确性和维修决策的科学性。针对传统寿命/退化试验方法存在的不足，加速退化试验逐步兴起，并成为了获取设备退化信息高效且经济的手段。

受环境和自身因素的影响，设备的退化过程呈现出显著的随机性，而加速试验则进一步增强了退化的随机性，因此采用具备时变不确定特征的Wiener过程刻画加速退化过程具备合理性。由于布朗运动在高应力(温度)条件下体现出更强的不确定性，因而现有研究多假设加速应力既影响Wiener过程的漂移系数又影响扩散系数，并已在电缆、LED、加速度计等设备的加速退化建模研究中得到应用和验证。然而，这些研究均认为加速应力仅与Wiener过程漂移系数存在函数关系，而将扩散系数视为不随应力变化的常数，从而忽略了应力越大退化不确定性越强的客观规律，导致此类方法剩余寿命的预测准确性不高。针对传统方法的缺陷，研究人员在加速退化建模过程中将漂移系数与扩散系数均视为应力的函数，并利用加速度计加速退化数据实现了对其剩余寿命的预测，提升了预测的准确性。然而，该方法将漂移系数与扩散系数视作独立变量分别建立加速模型，未能充分考虑漂移系数与扩散系数存在的关联关系，影响了预测精度的提升。为了实现对漂移系数和扩散系数关联特性的准确剖析，研究人员基于加速因子不变原则建立了考虑漂移系数与扩散系数比例关系的退化模型，从而有效降低了建模的不确定性，进一步提升了剩余寿命预测的有效性。然而，该模型无法应用于加速应力场合，且未考虑个体差异和测量误差对退化建模的影响。

在当前上述剩余寿命预测研究中，均需要多个退化试验样本才能确保参数估计与剩余寿命预测的准确性。但在现实环境中，出于节约试验成本的考虑和产品研制进度的制约，在研制阶段参与加速退化试验的机载设备样本数量往往很少，甚至可能只有一个。这种立足于单一设备从初始运行到当前时刻的状态监测数据来进行剩余寿命预测的方法也称为自适应剩余寿命预测方法，目前这方面的研究比较少，也无法应用于加速应力场合。

发明内容

本发明提供的基于比例加速退化建模的设备剩余寿命自适应预测方法，可以解决现有技术中存在的问题。

本发明提供了基于比例加速退化建模的设备剩余寿命自适应预测方法，包括以下步骤：

基于非线性Wiener过程构建比例加速退化模型；

在退化模型的基础上建立状态转移方程，并采用卡尔曼滤波KF算法更新设备退化状态；

输入设备的性能退化量的观测数据，采用期望最大化-卡尔曼滤波EM-KF算法实现对退化模型中未知参数的自适应估计；

在退化状态更新和未知参数自适应估计的基础上，基于全概率公式计算出设备剩余寿命的概率密度函数和累积分布函数。

本发明中的基于比例加速退化建模的设备剩余寿命自适应预测方法，针对加速退化试验条件下单台设备剩余寿命预测问题，首先基于非线性Wiener过程构建比例加速退化模型；其次，在退化模型的基础上建立状态转移方程，并采用卡尔曼滤波KF算法更新设备退化状态；再次，输入设备的性能退化量的观测数据，采用期望最大化-卡尔曼滤波EM-KF算法实现对退化模型中未知参数的自适应估计；最后，在退化状态更新和未知参数自适应估计的基础上，基于全概率公式计算出设备剩余寿命的概率密度函数和累积分布函数。利用本发明的方法，实现了单台设备加速退化试验下剩余寿命预测更加准确的效果。

附图说明

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

图1为实例中行波管加速退化数据；

图2为未知参数自适应估计过程；

图3为本发明方法和对比方法的退化状态更新过程对比；

图4为本发明方法和对比方法的退化状态预测误差对比；

图5为本发明方法和对比方法的剩余寿命预测结果对比；

图6为本发明方法和对比方法的扩散系数更新情况对比。

具体实施方式

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

本发明提供了基于比例加速退化建模的设备剩余寿命自适应预测方法，该方法的主要流程为：首先，基于非线性Wiener过程构建比例加速退化模型；其次，在退化模型的基础上建立状态转移方程，并采用KF(卡尔曼滤波)算法更新设备退化状态；然后，输入设备性能退化量的观测数据，采用EM-KF(期望最大化-卡尔曼滤波)算法实现对退化模型中未知参数的自适应估计；最后，在退化状态更新和未知参数自适应估计的基础上，基于全概率公式计算出设备剩余寿命的概率密度函数和累积分布函数。

下面对本发明方法的每个步骤进行详细说明。

比例加速退化模型

基本线性Wiener退化模型可以表示为：

X(t)＝X(0)+αt+βB(t) (1)

其中，X(0)表示设备在初始时刻的性能退化量，且常令X(0)＝0；α用以描述设备的退化速率，称之为漂移系数，由于设备退化具有差异性，则一般可用正态随机变量来对其进行描述，即β用以刻画设备退化过程的波动大小，称为扩散系数；而B(t)为标准布朗运动，以体现退化过程的时变不确定性，且B(t)～N(0,t)；一般情况下，B(t)与α相互独立。

进一步，考虑到外界环境对设备退化过程的影响，非线性退化模型适用性更强，因此在线性退化模型的基础上，通过对退化过程进行非线性化处理，即可得到更具一般意义的非线性退化模型：

式(2)也被称为时间尺度变换模型。其中，表示时间t的非线性函数，θ表示未知参数。

由于加速应力将同步影响退化模型的漂移系数和扩散系数，且加速应力的大小同退化速率与退化不确定性成正比关系，由此可得：

其中，S₁与S₂分别表示加速应力；是一个仅与应力大小有关的参数；/> 则分别对应不同加速应力条件的漂移系数与扩散系数。

则由式(3)易得：

其中，ρ为比例系数。由于加速应力S₁与S₂具有任意性，则对于式(2)所示退化过程，可知式(4)恒成立。

设备加速模型的一般形式可表示为：

α＝g(S|κ) (5)

其中，κ为未知参数。

则将式(4)与式(5)带入式(2)即可得到比例加速退化模型：

考虑到状态监测过程中由于测量方法、环境影响等对设备退化状态获取产生的不确定性，本发明将测量误差引入比例加速退化模型，可得：

其中，Y(t)表示设备退化量的观测值；v表示测量误差，且满足v～N(0,σ_v ²)；此外，一般认为v独立于B(t)与α。

基于KF的退化状态在线更新

为便于分析，本发明基于Exponential加速模型进行分析，其他加速模型分析过程与该模型相同，在此不再加以赘述。Exponential模型的一般表达式为：

g(S|κ)＝aexp(bS) (8)

其中，κ＝[a,b]；S为电应力。

考虑到不同设备的差异性，本发明令参数a为正态随机变量，且满足则可得：

接下来，本发明基于KF原理对设备的退化状态进行在线更新。假设t_i|S_j为目标设备的第i个监测时刻，且其对应的加速应力为S_j，则Y_i＝Y(t_i|S_j)与X_i＝X(t_i|S_j)分别为对应时刻设备性能退化量观测值和真实值；而Y_1:k＝[Y₁,Y₂,…Y₂]则表示直至t_k时刻已获取的全部退化数据。

联立式(7)与式(8)可得比例加速退化模型的状态转移方程，具体可表示为：

其中，t₀＝0；且易知/>

由于式(10)中存在非线性函数导致传统KF方法无法应用。为此，本发明对其进行线性化处理，可令：

L＝[1,0] (14)

由此可得：

在本发明中，令P_k|k分别表示退化状态的滤波均值与方差，对应的一步预测均值与方差则表示为/>P_k|k-1。而其具体定义式为：

基于上述分析，可给出KF的具体更新过程为：

P_k|k＝P_k|k-1-K_kLP_k|k-1 (21)

其中：

给定退化状态均值与方差的初始值与P_0|0，基于上述KF更新公式，即可实现对设备退化状态的在线更新。

基于EM-KF的参数自适应估计

针对退化模型中的未知参数，本发明在上述基于KF在线更新过程的基础上，利用EM算法对未知参数进行自适应估计。不妨令ψ表示退化模型中的未知参数，则易知

若已知设备性能退化量的观测数据Y_1:k，利用式(15)即可得到未知参数ψ关于设备性能退化状态Z_0:k与观测数据Y_1:k的联合对数似然函数：

其中P(Z_0:k,Y_1:k|ψ)为状态观测数据Y_1:k与设备退化状态的Z_0:k的联合概率密度函数。

基于上述分析可得：

Z_i|Z_i-1,ψ～N(A_iZ_i-1,Q_i) (27)

将式(26)～(29)带入式(25)，并除去常数项可得：

假设第j次迭代后退化模型参数的估计值为基于EM算法，可知第j+1次迭代的计算过程分为E步和M步。

E步：在第j次估计结果ψ^(j)的基础上，对似然函数L(ψ)求解关于隐含状态Z的期望，可得：

M步：对式(31)求极大值。由于式(31)隐含变量较多，无法直接对其进行最大化，为此，本发明采用RTS平滑器对式(31)进行处理。首先，给出RTS后向递归过程：

i＝k,k-1,…,0 (36)

其中，分别表示RTS平滑器的均值、方差以及协方差矩阵；i表示RTS后向平滑的迭代次数；D_i为RTS平滑器的滤波增益。

RTS的基本原理是令的初始值等于KF中系统状态量的最终估计值，进而反向迭代得到参数估计值，由此易知：

利用RTS基本原理可得：

基于上述分析，可将式(31)转为：

其中：

为求E(L(ψ))最大时对应的可令E(L(ψ))关于/>的偏导数等于零，由此可得：

将带入式(41)，并采用MATLAB软件中基于单纯形法的fminsearch函数求解其最大值，即可得到参数估计值b^(j+1),/>

不断对E步和M步进行迭代，直至|ψ^(j+1)-ψ^(j)|小于给定阈值停止迭代，即可实现对退化模型未知参数的自适应估计。

设备剩余寿命自适应预测

设备的寿命T通常被定义为从初始时刻运行直至性能退化量首次超过失效阈值D的时间，其数学表达式为：

T＝inf{t:X(t)≥D|X(0)＜D} (48)

若设备退化过程如式(2)所示，则其寿命的概率密度函数可表示为：

设备的剩余寿命可定义为从当前时刻t_k运行直至性能退化量首次超过失效阈值的时间。参照设备寿命的定义式，可知剩余寿命的定义式为：

L＝inf{l_k:X(t_k+l_k)≥D|X(0)＜D} (50)

若设备的退化过程如式(7)所示，则根据基于时间尺度模型的剩余寿命预测方法，可得其剩余寿命的条件分布函数为：

其中：

φ(l_k)＝Λ(t_k+l_k|v)-Λ(t_k|v) (53)

而φ^-1(·)表示φ(·)的逆函数。

由二维正态分布性质可知，a_k与X_k分别满足如下条件概率分布：

a_k|Y_1:k～N(E(a_k|Y_1:k),D(a_k|Y_1:k)) (55)

则基于全概率公式，即可实现对目标设备剩余寿命的在线预测，对应的剩余寿命概率密度函数与累积分布函数分别为：

实例分析

行波管是机载导航、雷达、电子对抗系统的核心部件，具备高可靠性、高价值、长寿命的特点。本发明基于某型行波管单台加速退化试验数据进行分析(如图1所示)，具体试验条要求为：

1)行波管的性能退化量选用阴极发射电流；

2)加速试验类型为恒定应力加速退化试验，且加速应力为电流密度，该试验中选用加速应力为8A/cm²；

3)当行波管的阴极发射电流下降至初始时刻的10％时可认为该行波管发生失效(对应的真实寿命为7000h)。

由图1可以发现，行波管的退化路径具有明显的非线性特征，且对该行波管性能退化增量进行K-S假设检验后不能拒绝其服从正态分布的假设，因而表明采用非线性Wiener过程对该行波管进行建模分析具有合理性。

实例参数自适应估计

工程经验表明，电子类产品退化过程近似满足幂函数。为此，本发明假设非线性函数设退化模型参数初值为μ_a＝0,/>b＝1,ρ＝1,θ＝0,/>则基于本发明提出的参数自适应估计方法，即可实现对参数的自适应估计，具体估计过程如图2所示。

由图2可知，除其余位置参数均可以较快的收敛至稳定值，表明本发明所提参数自适应估计算法具备较好的收敛性。且该算法运行算总时间约为0.0532s(运行环境：Intel Core I7-9750H处理器，16G内存，Windows7旗舰版操作系统，MATLAB软件)，表明该参数自适应估算方法具备较低的时间复杂度，拥有良好性能。

实例退化状在线更新

为便于分析，记本发明所提剩余寿命自适应预测方法为M0，而将不考虑比例关系的剩余寿命自适应预测方法记为M1。结合上述得到的退化模型参数自适应估计结果，基于KF原理即可对设备的退化状态进行在线更新，具体更新结果如图3所示。

图3中设备退化量的真实值，是将图1中初始时刻退化量设定为0后得到的。由图3可知，M0对应设备退化状态的预测值较M1与真实退化量更为接近，表明考虑漂移系数与扩散系数比例关系的退化模型更能反映设备的真实退化规律，具备更好的模型拟合性。为了更为直观的讨论M0与M1的差异，本发明给出不同方法对应退化状态预测结果的绝对误差，具体如图4所示。

由图4可知，M0的退化状态预测误差要显著小于M1。究其原因，主要是由于M1忽略了漂移系数与扩散系数的比例关系，从而导致该方法对设备退化状态进行估计的不确定性增大，进而产生了较大的误差。因此，有必要在退化建模过程中考虑漂移系数与扩散系数的比例关系。

实例剩余寿命自适应预测

在以上分析的基础上，基于本发明提出的剩余寿命自适应预测方法，即可对设备的剩余寿命进行预测。一般情况下，该型行波管的正常工作应力约为S0＝1A/cm²，其对应的剩余寿命预测曲线如图5所示。

由图5可知，M0与M1对应剩余寿命概率密度函数曲线均可以包含设备的真实剩余寿命，但M0对应曲线较M1明显更为集中，表明在确保剩余寿命准确预测的基础上，M0方法具有更低的预测不确定性，预测精度更高。进一步，本发明令：

β₁ ²＝E(a_k|Y_1:k)×exp(bS₀)×ρ

其中，β₁ ²等价于M0方法在常应力条件下的扩散系数。对应的，不妨设β²表示M1方法在常应力条件下的扩散系数。而β₁ ²与β²的更新过程如图6所示。

由图6可知，在扩散系数的更新全过程，M0对应的扩散系数较M1更小，表明利用比例关系模型可以更为准确的描述退化的时变不确定性，有效降低预测的不确定性，从而显著提升剩余寿预测方法的性能。该结论也进一步验证了以上结论。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (2)

1.基于比例加速退化建模的设备剩余寿命自适应预测方法，其特征在于，包括以下步骤：

基于非线性Wiener过程构建比例加速退化模型；

在退化模型的基础上建立状态转移方程，并采用卡尔曼滤波KF算法更新设备退化状态；

输入设备的性能退化量的观测数据，采用期望最大化-卡尔曼滤波EM-KF算法实现对退化模型中未知参数的自适应估计；

在退化状态更新和未知参数自适应估计的基础上，基于全概率公式计算出设备剩余寿命的概率密度函数和累积分布函数；

其中，基于非线性Wiener过程构建比例加速退化模型包括：

基本线性Wiener退化模型表示为：

X(t)＝X(0)+αt+βB(t) (1)

其中，X(0)表示设备在初始时刻的性能退化量，α用以描述设备的退化速率，称为漂移系数，β用以刻画设备退化过程的波动大小，称为扩散系数，B(t)为标准布朗运动；

通过对退化过程进行非线性化处理，即得到非线性退化模型：

其中，表示时间t的非线性函数，θ表示未知参数；

由于加速应力的大小同退化速率与退化不确定性成正比关系，由此可得：

其中，S₁与S₂分别表示加速应力，是一个仅与应力大小有关的参数，/>分别对应不同加速应力条件的漂移系数与扩散系数；

则由式(3)得：

其中，ρ为比例系数；

设备加速模型的一般形式表示为：

α＝g(S|κ) (5)

其中，κ为未知参数；

则将式(4)与式(5)带入式(2)即得到比例加速退化模型：

将测量误差引入比例加速退化模型，得：

其中，Y(t)表示设备退化量的观测值，v表示测量误差，且满足

其中，在退化模型的基础上建立状态转移方程包括：

基于Exponential加速模型进行分析，Exponential模型的一般表达式为：

g(S|κ)＝a exp(bS) (8)

其中，κ＝[a，b]；S为电应力；

令参数a为正态随机变量，且满足则得：

接下来基于KF原理对设备的退化状态进行在线更新，假设t_i|S_j为目标设备的第i个监测时刻，且其对应的加速应力为S_j，则Y_i＝Y(t_i|S_j)与X_i＝X(t_i|S_j)分别为对应时刻设备性能退化量观测值和真实值；而Y_1：k＝[Y₁，Y₂，…Y₂]则表示直至t_k时刻已获取的全部退化数据；

联立式(7)与式(8)得比例加速退化模型的状态转移方程，具体表示为：

其中，b为Exponential加速模型的形状参数；a_k-1为t_k-1时刻对应的Exponential加速模型的尺度参数；

对非线性函数进行线性化处理，令：

L＝[1，0] (14)

由此得：

令P_k|k分别表示退化状态的滤波均值与方差，对应的一步预测均值与方差则表示为/>P_k|k-1，KF的更新过程为：

P_k|k＝P_k|k-1-K_kLP_k|k-1 (17)

其中：

其中，E(ak-1|Y1:k-1)为a在tk-1时刻状态估计的期望；

给定退化状态均值与方差的初始值与P_0|0，基于公式(16)、(17)、(18)、(19)、(20)，实现对设备退化状态的在线更新；

其中，采用期望最大化-卡尔曼滤波EM-KF算法实现对退化模型中未知参数的自适应估计的方法为：

令ψ表示退化模型中的未知参数，则

输入设备的性能退化量的观测数据Y_1：k，利用式(15)即得到未知参数ψ关于设备性能退化状态Z_0：k与观测数据Y_1：k的联合对数似然函数：

其中P(Z_0：k，Y_1：k|ψ)为状态观测数据Y_1：k与设备退化状态的Z_0：k的联合概率密度函数；

基于上述分析得：

Z_i|Z_i-1，ψ～N(A_iZ_i-1，Q_i) (23)

将式(22)～(25)带入式(21)，并除去常数项得：

假设第j次迭代后退化模型参数的估计值为基于EM算法，对其中的E步和M步进行迭代直至|ψ⁽ⁱ⁺¹⁾-ψ^(j)|小于给定阈值，实现对退化模型未知参数的自适应估计；其中M步对E步得到的隐含状态期望进行最大化时，采用RTS平滑器对隐含状态进行处理，在处理结果上采用卡尔曼滤波KF算法进行反向迭代，得到未知参数的估计值；

其中，基于全概率公式计算出设备剩余寿命的概率密度函数和累积分布函数包括：

剩余寿命的定义式为：

L＝inf{l_k：X(t_k+l_k)≥D|X(0)＜D} (44)

其中，D为失效阈值，若设备的退化过程如式(7)所示，则根据基于时间尺度模型的剩余寿命预测方法，得到剩余寿命的条件分布函数为：

其中：

φ(l_k)＝Λ*(t_k+l_k|v)-Λ(t_k|v) (47)

其中，φ^-1(·)表示φ(·)的逆函数；

由二维正态分布性质，a_k与X_k分别满足如下条件概率分布：

a_k|Y_1k～N(E(a_k|Y_1k)，D(a_k|Y_1k)| (49)

则基于全概率公式，即可实现对目标设备剩余寿命的在线预测，对应的剩余寿命概率密度函数与累积分布函数分别为：

其中，l表示剩余寿命，l_k表示t_k时刻设备对应的剩余寿命。

2.如权利要求1所述的基于比例加速退化建模的设备剩余寿命自适应预测方法，其特征在于，所述EM算法的E步和M步分别为：

E步：在第j次估计结果ψ^(j)的基础上，对似然函数L(ψ)求解关于隐含状态Z的期望，得：

M步：对式(27)求极大值；由于式(27)隐含变量较多，采用RTS平滑器对式(31)进行处理；首先，RTS后向递归过程为：

i＝k，k-1，…，0 (32)

其中，分别表示RTS平滑器的均值、方差以及协方差矩阵；i表示RTS后向平滑的迭代次数；D_i为RTS平滑器的滤波增益；

RTS的基本原理是令的初始值等于KF中系统状态量的最终估计值，进而反向迭代得到参数估计值，由此得：

利用RTS基本原理得：

基于上述分析，将式(27)转为：

其中：

为求E(L(ψ))最大时对应的令E(L(ψ))关于/>的偏导数等于零，由此得：

将带入式(37)，并采用MATLAB软件中基于单纯形法的fminsearch函数求解其最大值，即得到参数估计值b^(j+1)，/>