CN112949204A - 一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法 - Google Patents

Abstract

一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法，先建立建立与滚动轴承累计旋转圈数有关的状态空间模型，再确定三个阈值，并基于阈值确定起始退化时间和起始预测时间；然后基于起始退化时间之后的观测值序列，由极大似然法更新模型参数，利用贝叶斯信息准则自动选择最优状态空间模型；最后根据选择的模型及其更新后的参数，计算失效时间的概率密度函数，得到剩余寿命预测结果；本发明能够由滚动轴承自身状态监测数据驱动而无需参考任何其他滚动轴承监测数据，并将监测数据实时匹配最优状态空间模型，有效表征工业实际中滚动轴承运行过程中的退化情况，提高了滚动轴承的剩余寿命预测精度。

Description

一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法

技术领域

本发明属于滚动轴承剩余寿命预测技术领域，具体涉及一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法。

背景技术

近年来，随着军事技术的迅速发展和现代工业水平的整体提升，对机械设备安全运行的要求也越来越高。设备中的核心零部件极易发生不同程度的故障，导致设备失效，造成经济损失和人员伤亡。因此，保证机械设备关键零部件安全服役，对保障国防安全和促进现代工业生产意义重大。滚动轴承在绝大多数滚动轴承系统中扮演着重要角色，被称作“工业的关节”，其健康状态的好坏直接影响到机械设备的运转效率。实现对滚动轴承剩余寿命的准确预测，提前更换剩余寿命不足的轴承，从而保证机械设备的安全具有重要意义。

现如今大部分的数据驱动剩余寿命预测方法依赖退化过程的经验知识和历史失效数据来训练预测模型，这些方法有两个前提假设：1、能够提前获取大量失效数据；2、不同滚动轴承的退化行为是相似的。然而，对于实际工程设备，往往只有少量失效数据或没有失效数据；即使积累了足够的失效数据，由于制造差异、环境、工况等因素，不同滚动轴承之间的退化过程也往往有所差异，很难保证训练好的模型在不同设备上的适应性。因此，上述两个假设在在一些实际情况中是不合理的，导致基于此假设的数据驱动剩余寿命预测方法效果不好。

发明内容

为了克服现有技术的以上缺点，本发明的目的在于提供了一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法，滚动轴承剩余寿命根据其自身的状态监测数据而被预测，而没有来自其他任何滚动轴承的监测数据，减少了模型对训练数据相似性和数量的要求，并能够使监测数据实时匹配最优状态空间模型，提高了滚动轴承剩余寿命的预测精度。

为了达到上述目的，本发明采取的技术方案为：

一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法，包括以下步骤：

第1步，建立与滚动轴承累计旋转圈数有关的状态空间模型：

y(t)＝x(t，0)+ωx(t，θ) (2)

其中，公式(1)为状态方程，公式(2)为观测方程，x(t，θ)为t时刻的健康状态，x₀为滚动轴承初始健康状态，

是t时刻的累计退化量，退化函数

根据实际情况设为线性函数、幂函数、指数函数或对数函数，θ是退化函数

的参数向量，r(t)代表t时刻滚动轴承的累计旋转圈数，其计算公式为

其中s(τ)代表轴承转速；

y(t)是t时刻的观测值，能够反应状态退化的特征指标；ωx(t，θ)代表观测噪声，ω服从正态分布N(0，σ²)，因此，y(t)服从正态分布N(x(t，θ)，σ²(x(t，θ))²)；

第2步，确定三个阈值：预警阈值D₁，报警阈值D₂和失效阈值D₃，这三个阈值将观测值划分为4个不同的级别；

2.1)预警阈值D₁是级别I和II的分界线，分别代表健康阶段和初始退化阶段；在滚动轴承早期运行阶段，将记录的观测值作为健康阶段下的随机样本，计算这些观测值的3σ区间，并选择其上边界或下边界作为预警阈值D₁，选择3σ区间的上边界用于观测值随着退化过程增加的情况，而3σ区间的下边界用于观测值随着退化过程减少的情况；

2.2)失效阈值D₃是等级III和IV之间的边界，即严重退化阶段和失效阶段，需要人为根据实际情况进行预先指定；

2.3)报警阈值D₂是级别II和III之间的边界，定义为预警阈值D₁和失效阈值D₃的平均值；

第3步，实时监测并采集滚动轴承的振动信号，计算观测值序列，确定首次退化时间t_FDT和首次预测时间t_FPT；当观测值达到预警阈值D₁，将这个时间点定义为滚动轴承的起始退化时间FDT；当观测值达到报警阈值D₂，将这个时间点定义为起始预测时间FPT；

第4步，状态空间模型参数更新：

状态空间模型的未知参数由起始退化时间FDT之后的一系列观测值Y＝(y_FDT，y_FDT+1，…，y_k)′确定，随着观测值序列不断增加，状态空间模型参数不断更新；

由极大似然估计法估计模型参数θ和σ，未知参数θ和σ²的估计值是使下面对数似然函数最大化的值：

其中，θ和σ²为未知参数，Y为观测值序列，k代表Y的观测值数量，

是t_i时刻轴承的累计退化量，y_i代表t_i时刻的观测值，初始健康状态x₀设置为预警阈值D₁，r_i是滚动轴承从起始退化时间t_FDT到时刻t_i的累计转动量，计算公式如下：

其中，FDT代表起始退化时间，s_j代表t_j时刻轴承转速；

计算公式(3)关于σ²的偏导数，并且令

则σ²的估计值由θ表示，公式如下：

将公式(5)带入公式(3)，将对数似然函数简化为仅含未知参数θ的形式：

基于公式(6)，使用下山单纯形法求解θ的估计结果，然后，通过将θ的估计输入到公式(5)中来计算σ²的估计结果；

第5步，最优状态空间模型自动选择：

伴随着观测值序列的不断增加，基于贝叶斯信息准则(BIC)评估建立的状态空间模型的拟合性能，自动选择当前最优状态空间模型，BIC计算公式如下：

其中，

代表似然函数

的最大值，N_P＝N_θ+1代表未知参数的个数，N_θ代表θ中参数的个数；选择获得最小BIC分数的模型作为当前最优状态空间模型；

第6步，计算失效时间的概率密度函数：

失效时间T看作是一个时间变量，定义为当观测值y(t)超过失效阈值D₃的时刻，其累计分布函数表示为：

F(t)＝Pr(T≤t)＝Pr(y(t)≥D₃) (8)

其中，Pr()表示某事件发生的概率；

由式(2)知观测值y(t)服从正态分布N(x(t，θ)，σ²(x(t，θ))²)，令

则z服从标准正态分布，其概率密度函数和累计分布函数分别表示为ψ(z)和Ψ(z)，失效时间的累计分布函数表示如下：

其中

x(t，θ)为t时刻的健康状态，Ψ(g(t))表示t时刻的累计分布函数；

失效时间T不能为负，T＜0的概率需要从失效时间的累计分布函数中截去，截去后的累计分布函数表示如下：

其中，Ψ(g(0))表示t＝0时刻的累计分布函数；

失效时间T的概率密度函数表示如下：

其中：

是t时刻轴承的累计退化量；

第7步，计算剩余寿命：

将第5步中确定的最优状态空间模型的表达式及其更新后的参数带入公式(11)得到失效时间T的概率密度函数，剩余寿命预测结果L计算如下：

L＝T-t_k+t_FDT (14)

其中，T表示失效时间，t_k代表当前时刻，t_FDT代表首次退化时间。

本发明的有益效果为：

本发明提出了一种数据与模型自适应匹配的方法来预测滚动轴承的剩余寿命，由滚动轴承自身状态监测数据驱动而无需参考任何其他滚动轴承监测数据，考虑了不同滚动轴承退化过程之间的差异性，降低了数据驱动方法对训练数据量的依赖性。本发明方法能有效表征工业实际中滚动轴承运行过程中的退化情况，使监测数据实时匹配最优状态空间模型，提高了滚动轴承的剩余寿命预测精度。

附图说明

图1为本发明方法流程图。

图2为三个阈值和起始退化时间、起始预测时间示意图。

图3为四组模拟退化轨迹示意图。

图4为本发明方法对四种模拟退化轨迹在160天时的剩余寿命预测结果。

图5为本发明方法对四种模拟轨迹在全预测周期的预测结果。

具体实施方式

下面结合附图和实施例对本发明做进一步阐述。

参照图1，一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法，包括以下步骤：

第1步，建立与滚动轴承累计旋转圈数有关的状态空间模型：

y(t)＝x(t，θ)+ωx(t，θ) (2)

其中，公式(1)为状态方程，公式(2)为观测方程，x(t，θ)为t时刻的健康状态，x₀为滚动轴承初始健康状态，

是t时刻的累计退化量，退化函数

根据实际情况设为线性函数、幂函数、指数函数或对数函数，θ是退化函数

的参数向量，r(t)代表t时刻滚动轴承的累计旋转圈数，其计算公式为

其中s(τ)代表轴承转速；

y(t)是t时刻的观测值，能够反应状态退化的特征指标；ωx(t，θ)代表观测噪声，ω服从正态分布N(0，σ²)，因此，y(t)服从正态分布N(x(t，θ)，σ²(x(t，θ))²)；

第2步，参照图2，确定三个阈值：预警阈值D₁，报警阈值D₂和失效阈值D₃，这三个阈值将观测值划分为4个不同的级别；

2.1)预警阈值D₁是级别I和II的分界线，分别代表健康阶段和初始退化阶段；在滚动轴承早期运行阶段，将记录的观测值作为健康阶段下的随机样本，计算这些观测值的3σ区间，并选择其上边界或下边界作为预警阈值D₁，选择3σ区间的上边界用于观测值随着退化过程增加的情况，而3σ区间的下边界用于观测值随着退化过程减少的情况；

2.2)失效阈值D₃是等级III和IV之间的边界，即严重退化阶段和失效阶段，需要人为根据实际情况进行预先指定；

2.3)报警阈值D₂是级别II和III之间的边界，定义为预警阈值D₁和失效阈值D₃的平均值；

第3步，实时监测并采集滚动轴承的振动信号，计算观测值序列，确定首次退化时间t_FDT和首次预测时间t_FPT；当观测值达到预警阈值D₁，将这个时间点定义为滚动轴承的起始退化时间FDT；当观测值达到报警阈值D₂，将这个时间点定义为起始预测时间FPT；

第4步，状态空间模型参数更新：

状态空间模型的未知参数由起始退化时间FDT之后的一系列观测值Y＝(y_FDT，y_FDT+1，…，y_k)′确定，随着观测值序列不断增加，状态空间模型参数不断更新；

由极大似然估计法估计模型参数θ和σ，未知参数θ和σ²的估计值是使下面对数似然函数最大化的值：

其中，θ和σ²为未知参数，Y为观测值序列，k代表Y的观测值数量，

是t_i时刻轴承的累计退化量，y_i代表t_i时刻的观测值，初始健康状态x₀设置为预警阈值D₁，r_i是滚动轴承从起始退化时间t_FDT到时刻t_i的累计转动量，计算公式如下：

其中，FDT代表起始退化时间，s_j代表t_j时刻轴承转速；

计算公式(3)关于σ²的偏导数，并且令

则σ²的估计值可以由θ表示，公式如下：

将公式(5)带入公式(3)，将对数似然函数简化为仅含未知参数θ的形式：

基于公式(6)，使用下山单纯形法求解θ的估计结果，然后，通过将θ的估计输入到公式(5)中来计算σ²的估计结果；

第5步，最优状态空间模型自动选择：

伴随着观测值序列的不断增加，基于贝叶斯信息准则(BIC)评估建立的状态空间模型的拟合性能，自动选择当前最优状态空间模型，BIC计算公式如下：

其中，

代表似然函数

的最大值，N_P＝N_θ+1代表未知参数的个数，N_θ代表θ中参数的个数；选择获得最小BIC分数的模型作为当前最优状态空间模型；

第6步，计算失效时间的概率密度函数：

失效时间T看作是一个时间变量，定义为当观测值y(t)超过失效阈值D₃的时刻，其累计分布函数表示为：

F(t)＝Pr(T≤t)＝Pr(y(t)≥D₃) (8)

其中，Pr()表示某事件发生的概率；

由式(2)知观测值y(t)服从正态分布N(x(t，θ)，σ²(x(t，θ))²)，令

则z服从标准正态分布，其概率密度函数和累计分布函数分别表示为ψ(z)和Ψ(z)，失效时间的累计分布函数表示如下：

其中

x(t，θ)为t时刻的健康状态，Ψ(g(t))表示t时刻的累计分布函数；

失效时间T不能为负，T＜0的概率需要从失效时间的累计分布函数中截去，截去后的累计分布函数表示如下：

其中，Ψ(g(0))表示t＝0时刻的累计分布函数；

失效时间T的概率密度函数表示如下：

其中：

是t时刻轴承的累计退化量；

第7步，计算剩余寿命：

将第5步中确定的最优状态空间模型的表达式及其更新后的参数带入公式(11)得到失效时间T的概率密度函数，剩余寿命预测结果L计算如下：

L＝T-t_k+t_FDT (14)

其中，T表示失效时间，t_k代表当前时刻，t_FDT代表首次退化时间。

为了证明本发明方法的有效性，模拟了4种不同退化轨迹的振动信号有效值(RMS)，结合本发明方法开展剩余寿命预测。

为了产生两阶段的滚动轴承退化过程，将退化过程假设为分段函数，如公式(15)所示：

其中x₀为5m/s²，转速r(t)恒定为1200rpm，起始退化时间t_FDT为100天，

假设为四种常用的函数，分别是线性函数、幂函数、指数函数和对数函数，图3显示了四组不同的模拟退化轨迹。

四组模拟退化轨迹观测值连同工作时间一起输入到本发明方法的状态空间模型中，每种信号的剩余寿命概率密度函数由其自身的观测值进行实时预测，为了使预测过程可视化，四组模拟退化轨迹在第160天的剩余寿命预测结果如图4所示。

将每组模拟退化轨迹的观测值输入到本发明方法的状态空间模型中，将会自动选择最优状态空间模型，并根据每个轨迹的退化趋势自适应地估计模型参数。可以看出，状态空间模型对预测轨迹的期望与160天后观测值的趋势一致，这表明已获得正确的状态空间模型和准确的参数；预测轨迹的95％置信区间也与160天后观测值重叠。由于观测存在随机噪声，观测值会在一个时间间隔内围绕故障阈值波动，将观察值第一次超过阈值与观察值低于阈值的最后时间的时间区间定义为失效时段(FTI)。从图4中可以看出，预测故障时间的概率密度函数覆盖了失效时段FTI，这表明了预测结果的准确性。

为了评估本发明方法在整个预测过程中的性能，计算了预测的剩余寿命概率密度函数的期望值和95％置信区间，将其与有效预测区间进行比较，并提供了四个模拟退化轨迹的状态空间模型选择结果，如图5所示，有效预测区间EPB＝FTI-FDT；可以看出，预测的剩余寿命的概率密度函数与起始预测时间FPT的有效预测区间EPB高度吻合，并且在退化过程中每个时刻的最优状态空间模型也已正确选择，这表明本发明方法在剩余寿命预测中的有效性。

应当指出，在不脱离本发明构想的前提下，对本发明方法所做的调整和变形，也应视为本发明的保护范围。

Claims (1)

1.一种数据与模型自适应匹配的滚动轴承剩余寿命预测方法，其特征在于，包括以下步骤：

第1步，建立与滚动轴承累计旋转圈数有关的状态空间模型：

其中，公式(1)为状态方程，公式(2)为观测方程，x(t，θ)为t时刻的健康状态，x₀为滚动轴承初始健康状态，

是t时刻的累计退化量，退化函数

根据实际情况设为线性函数、幂函数、指数函数或对数函数，θ是退化函数

的参数向量，r(t)代表t时刻滚动轴承的累计旋转圈数，其计算公式为

其中s(τ)代表轴承转速；

y(t)是t时刻的观测值，能够反应状态退化的特征指标；ωx(t，θ)代表观测噪声，ω服从正态分布N(0，σ²)，因此，y(t)服从正态分布N(x(t，θ)，σ²(x(t，θ))²)；

第2步，确定三个阈值：预警阈值D₁，报警阈值D₂和失效阈值D₃，这三个阈值将观测值划分为4个不同的级别；

2.1)预警阈值D₁是级别I和II的分界线，分别代表健康阶段和初始退化阶段；在滚动轴承早期运行阶段，将记录的观测值作为健康阶段下的随机样本，计算这些观测值的3σ区间，并选择其上边界或下边界作为预警阈值D₁，选择3σ区间的上边界用于观测值随着退化过程增加的情况，而3σ区间的下边界用于观测值随着退化过程减少的情况；

2.2)失效阈值D₃是等级III和IV之间的边界，即严重退化阶段和失效阶段，需要人为根据实际情况进行预先指定；

2.3)报警阈值D₂是级别II和III之间的边界，定义为预警阈值D₁和失效阈值D₃的平均值；

第3步，实时监测并采集滚动轴承的振动信号，计算观测值序列，确定首次退化时间t_FDT和首次预测时间t_FPT；当观测值达到预警阈值D₁，将这个时间点定义为滚动轴承的起始退化时间FDT；当观测值达到报警阈值D₂，将这个时间点定义为起始预测时间FPT；

第4步，状态空间模型参数更新：

状态空间模型的未知参数由起始退化时间FDT之后的一系列观测值Y＝(y_FDT，Y_FDT+1，…，y_k)′确定，随着观测值序列不断增加，状态空间模型参数不断更新；

由极大似然估计法估计模型参数θ和σ，未知参数θ和σ²的估计值是使下面对数似然函数最大化的值：

其中，θ和σ²为未知参数，Y为观测值序列，k代表Y的观测值数量，

是t_i时刻轴承的累计退化量，y_i代表t_i时刻的观测值，初始健康状态x₀设置为预警阈值D₁，r_i是滚动轴承从起始退化时间t_FDT到时刻t_i的累计转动量，计算公式如下：

其中，FDT代表起始退化时间，s_j代表t_j时刻轴承转速；

计算公式(3)关于σ²的偏导数，并且令

则σ²的估计值可以由θ表示，公式如下：

将公式(5)带入公式(3)，将对数似然函数简化为仅含未知参数θ的形式：

基于公式(6)，使用下山单纯形法求解θ的估计结果，然后，通过将θ的估计输入到公式(5)中来计算σ²的估计结果；

第5步，最优状态空间模型自动选择：

伴随着观测值序列的不断增加，基于贝叶斯信息准则(BIC)评估建立的状态空间模型的拟合性能，自动选择当前最优状态空间模型，BIC计算公式如下：

其中，

代表似然函数

的最大值，N_P＝N_θ+1代表未知参数的个数，N_θ代表θ中参数的个数；选择获得最小BIC分数的模型作为当前最优状态空间模型；

第6步，计算失效时间的概率密度函数：

失效时间T看作是一个时间变量，定义为当观测值y(t)超过失效阈值D₃的时刻，其累计分布函数表示为：

F(t)＝Pr(T≤t)＝Pr(y(t)≥D₃) (8)

其中，Pr()表示某事件发生的概率；

由式(2)知观测值y(t)服从正态分布N(x(t，θ)，σ²(x(t，θ))²)，令

则z服从标准正态分布，其概率密度函数和累计分布函数分别表示为ψ(z)和Ψ(z)，失效时间的累计分布函数表示如下：

其中

x(t，θ)为t时刻的健康状态，Ψ(g(t))表示t时刻的累计分布函数；

失效时间T不能为负，T＜0的概率需要从失效时间的累计分布函数中截去，截去后的累计分布函数表示如下：

其中，Ψ(g(0))表示t＝0时刻的累计分布函数；

失效时间T的概率密度函数表示如下：

其中：

是t时刻轴承的累计退化量；

第7步，计算剩余寿命：

将第5步中确定的最优状态空间模型的表达式及其更新后的参数带入公式(11)得到失效时间T的概率密度函数，剩余寿命预测结果L计算如下：

L＝T-t_k+t_FDT (14)

其中，T表示失效时间，t_k代表当前时刻，t_FDT代表首次退化时间。

Images