CN105653851A - 基于分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法 - Google Patents

Abstract

本发明公开了一种基于分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法，其特征包括如下步骤：1)采集滚动轴承运行状态振动信号，从中提取特征数值；2)设置报警阈值、分阶段点阈值、停车阈值分别作为预测初始点、第二阶段和第三阶段分界点、预测截止点。3)设置粒子滤波算法初始值；4)建立不同劣化阶段滚动轴承初始模型；5)设置模型初始参数数值；6)应用最小二乘法更新模型参数；7)应用上述不同劣化阶段物理模型及粒子滤波算法进行滚动轴承剩余寿命预测；8)滚动轴承剩余寿命预测输出。该方法具有实时性、有效性、准确率较高等优点，适合于多种设备的滚动轴承剩余寿命预测。

Description

基于分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法

技术领域

本发明属于设备故障诊断领域，涉及针对常用滚动轴承剩余寿命预测模型建立的方式方法。

背景技术

滚动轴承是旋转机械中应用广泛的零部件之一，其寿命的长短影响设备的正常运转。准确预测滚动轴承的剩余寿命，可以合理安排滚动轴承的性能检查及维护更换工作，避免停车过早影响生产的进行或停车过晚造成设备损坏，影响生产安全。目前如何有效进行滚动轴承剩余寿命预测是目前研究的一大难点，众多工作人员及科研人员都在这一研究上进行了很多探索。将物理方法和数据驱动方法有效结合，按照裂纹劣化不同阶段分别进行滚动轴承剩余寿命预测是一大研究亮点。

目前滚动轴承故障导致轴承寿命缩短的原因有很多，主要原因如下：滚动轴承长期受交变载荷作用引起疲劳剥落；滚动轴承零部件之间因为相对滑动儿引起的磨损；在外力和环境温度作用下发生塑性变形；由于环境介质引起的滚动轴承零件金属表面发生化学反应造成的腐蚀；滚动轴承由于过载或缺陷导致断裂等突发性事故及保持架损坏等多种故障。

滚动轴承剩余寿命预测方法目前有很多种，主要包括：基于Paris疲劳寿命预测模型的轴承寿命预测方法，应用概率理论进行轴承的寿命预测，基于人工智能技术的轴承疲劳寿命的研究及基于状态监测系统的疲劳寿命预测研究等。Paris理论公式符合第二个区域的劣化规律。在第三区域，裂纹劣化速度明显加快，时间很短，但是对于低周期的疲劳部件，时间不算很短，且第三部分的准确预测也能提高滚动轴承剩余寿命预测的准确性，便于精准计算滚动轴承更换时间。为提高剩余寿命预测结果的准确性，将滚动轴承运行状态的差异分为第二阶段及第三阶段两种情况，则不同的运行阶段对应不同的寿命预测模型。

近年来，随着科技的发展进步，新兴的人工智能诊断技术及状态监测系统正在得到逐步发展和应用，此种方法适用于解决物理规律复杂、不确定性影响因素较多的滚动轴承剩余寿命预测情况。物理模型主要从故障机理出发，模型劣化规律与实际工况相似。将二者相结合可以很大程度上提高滚动轴承剩余寿命预测结果的准确性。本文依据上述内容提出了将人工智能技术及物理模型相结合的分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法。

发明内容

本发明的目的是考虑滚动轴承不同劣化阶段的特点，提出不同劣化阶段的预测模型，应用其进行滚动轴承剩余寿命预测，很大程度上提高了滚动轴承剩余寿命预测结果的准确性。该方法结合现代智能算法，为滚动轴承寿命预测提供了一套智能的、实时的、准确的滚动轴承剩余寿命预测方法。

本发明公开了一种基于分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法，具体步骤为：

步骤(1)采集滚动轴承运行状态振动信号，从中提取特征数值：应用信号采集系统采集滚动轴承运行状态信息，将其输入特征数值公式，求取特征数值作为滚动轴承剩余寿命预测的输入值。

步骤(2)设置报警阈值、分阶段点阈值、停车阈值分别作为预测初始点、第二阶段和第三阶段分界点、预测截止点：将运行状态振动区域按照标准分为A、B、C、D4个区域，设置报警阈值为B区域上限的1.25倍，设置分阶段点阈值为C区域上限数值，设置停车阈值为C区域上限的1.25倍。

步骤(3)设置粒子滤波算法初始值；

1)取预测时刻点前100～200个数据点进行预测；

2)特征数值初始值为采集到的振动波形数据计算出的第一个特征数值；

3)按照公式z_t+1＝z_t+randn获取服从高斯分布的单个时间点的N个粒子，randn为服从高斯分布的0～1之内的随机数；z_t即为t时刻的粒子,z_t+1即为t+1时刻粒子。

4)设置粒子的初始权重均为

5)设置过程噪声w(t)及观测噪声v(t)的数值分别为0～1之间的数值。

步骤(4)建立不同劣化阶段滚动轴承初始模型。

本次预测针对滚动轴承中晚期运行状态进行预测。

通过对Paris理论公式及Foreman理论公式进行改进，获取物理模型和状态空间模型相结合的分阶段滚动轴承剩余寿命预测模型，预测模型如下：

第二阶段：

状态方程： $\frac{V(t+1)}{S-S_{or}}=k_1\×q^n\×{\sqrt{\mathrm{\π}}}^n\×{(V_{dq}-V_{nr})}^n+w\left(t\right)$ (10)

观测方程： $y(t+1)=\frac{V(t+1)}{S-S_{or}}\×(S-S_{or})+v\left(t\right)$ (11)

第三阶段：

状态方程： $\frac{V(t+1)}{S-S_{or}}=(k_1\×q^n\×{\sqrt{\mathrm{\π}}}^n\×{\left(V_{dq}-V_{nr}\right)}^n+w(t\left)\right)\×\frac{x\left(t\right)}{V_{ed}-x\left(t\right)}$ (12)

观测方程： $y(t+1)=\frac{V(t+1)}{S-S_{or}}\×(S-S_{or})+v\left(t\right)$ (13)

其中，V_dq表示当前振动数值，本次预测振动数值采用包络谱峰值，但实际预测振动数值可采用多种振动特征数值，例如振动峰值、有效值等；V_nr表示正常运转振动数值；S_or指代开始劣化时应力循环次数；S指代应力循环次数；V(t+1)为轴承t+1时刻的特征数值大小；w(t)在t时刻状态转移的扰动，为轴承运转过程中产生的过程噪声向量数值；v(t)为t时刻轴承运转过程中产生的观测噪向量值；y(t+1)即为最终状态预测数值；k₁是第一材料系数，n是第二材料系数，二者依据材料特性和实验环境发生变化，是变量；q指代载荷；V_ed为振动停车阈值，V_max为目前采集到的振动数值的最大值，即为x(t)；其余参数同第二阶段参数。

方程式中k₁、n的确定方法为：

1)取预测时刻点前100～200个特征数值；

2)将上述特征数值带入滚动轴承剩余寿命预测模型；

3)应用最小二乘法拟合特征数值及应力循环次数之间的非线性关系，即在满足最小均方误差的情况下，求取出上述状态方程中的k₁、n；此处申明，k₁、n是随数据的变化实时更新的。

步骤(5)设置模型初始参数数值；

设置预测步长为10，设置动态状态空间模型材料常数k₁、n初始值分别为2.0；1×10^-3，正确设置模型参数数值可以缩短运算时间，增加预测结果的准确性。

步骤(6)当特征数值大于等于报警阈值，将预测时刻点前100～200个特征数值带入滚动轴承寿命预测算法；应用最小二乘法更新初始模型参数，从而使模型符合滚动轴承实际运行状态。当特征数值大于等于分阶段点阈值，更新第三阶段模型参数，否则，更新第二阶段模型参数。

步骤(7)应用上述分阶段物理模型及粒子滤波算法进行滚动轴承剩余寿命预测。

预测过程中，当特征数值大于等于报警阈值小于分阶段点阈值时，应用参数更新后的第二阶段模型及粒子滤波算法进行滚动轴承运行状态评估；当特征数值大于等于分阶段点阈值小于停车阈值时，应用参数更新后的第三阶段模型及粒子滤波算法进行滚动轴承运行状态评估。

预测时每隔100分钟应用最小二乘法更新模型参数，修正模型轨迹，从而使模型符合滚动轴承实际运行状态。

步骤(8)当特征数值小于停车阈值时，预测循环进行，否则，预测循环停止，输出滚动轴承运行状态特征数值，并依据下述公式求取滚动轴承剩余寿命：

$L=(D-d)\×\frac{f_n}{f_s}\×t_{\min }$ (14)

L指代剩余寿命

D指代预测的特征数值点数与输入的特征数值的点数的总和

d指代数据提供的输入的特征数值的点数

f_n指代采样点数

f_s指代采样频率

t_min指代输入的两个特征数值之间的时间间隔；

本发明可以进行滚动轴承的实时预测，使滚动轴承在工作状态下即可进行工作，适用性强，可适用于多种工作场合；

本发明公开了一种基于分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法，其中状态方程的数据是可以实时更新替换的，状态方程中的未知参数应用最小二乘法原理求取。

附图说明

图1：特征值数据变化趋势图

图2:滚动轴承寿命预测结果图

图3:滚动轴承寿命预测流程

具体实施方式

下面将结合附图对本发明的具体的寿命预测流程做进一步说明。

如图3所示，本发明的具体流程如下所示：

1、采集滚动轴承运行状态振动信号，从中提取特征数值：应用信号采集系统采集滚动轴承运行状态信息，将其输入特征数值公式，求取特征数值作为滚动轴承剩余寿命预测的输入值。

2、ISO2372及ISO7919标准可知，振动标准一共分为A、B、C、D4个区域。设置报警阈值、分阶段点阈值、停车阈值分别作为预测初始点、第二阶段和第三阶段分界点、预测截止点：将运行状态振动区域按照标准分为A、B、C、D4个区域，设置报警阈值为B区域上限的1.25倍，设置分阶段点阈值为C区域上限数值，设置停车阈值为C区域上限的1.25倍。

3、设置粒子滤波算法初始值；

1)取预测时刻点前100～200个数据点进行预测；

2)特征数值初始值为采集到的振动波形数据计算出的第一个特征数值；

3)按照公式z_t+1＝z_t+randn获取服从高斯分布的单个时间点的N个粒子，randn为服从高斯分布的0～1之内的随机数；z_t即为t时刻的粒子,z_t+1即为t+1时刻粒子。

(4)设置粒子的初始权重均为

(5)设置过程噪声w(t)及观测噪声v(t)的数值分别为0～1之间的数值。

4、建立不同劣化阶段滚动轴承初始模型。

本次预测针对滚动轴承中晚期运行状态进行预测。

通过对Paris理论公式及Foreman理论公式进行改进，获取物理模型和状态空间模型相结合的分阶段滚动轴承剩余寿命预测模型，预测模型如下：

第二阶段：

状态方程： $\frac{V(t+1)}{S-S_{or}}=k_1\×q^n\×{\sqrt{\mathrm{\π}}}^n\×{(V_{dq}-V_{nr})}^n+w\left(t\right)$ (15)

观测方程： $y(t+1)=\frac{V(t+1)}{S-S_{or}}\×(S-S_{or})+v\left(t\right)$ (16)

第三阶段：

状态方程： $\frac{V(t+1)}{S-S_{or}}=(k_1\×q^n\×{\sqrt{\mathrm{\π}}}^n\×{\left(V_{dq}-V_{nr}\right)}^n+w(t\left)\right)\×\frac{x\left(t\right)}{V_{ed}-x\left(t\right)}$ (17)

观测方程： $y(t+1)=\frac{V(t+1)}{S-S_{or}}\×(S-S_{or})+v\left(t\right)$ (18)

其中，V_dq表示当前振动数值，本次预测振动数值采用包络谱峰值，但实际预测振动数值可采用多种振动特征数值，例如振动峰值、有效值等；V_nr表示正常运转振动数值；S_or指代开始劣化时应力循环次数；S指代应力循环次数；V(t+1)为轴承t+1时刻的特征数值大小；w(t)在t时刻状态转移的扰动，为轴承运转过程中产生的过程噪声向量数值；v(t)为t时刻轴承运转过程中产生的观测噪向量值；y(t+1)即为最终状态预测数值；k₁是第一材料系数，n是第二材料系数，二者依据材料特性和实验环境发生变化，是变量；_q指代载荷；V_ed为振动停车阈值，V_max为目前采集到的振动数值的最大值，即为x(t)；其余参数同第二阶段参数。

方程式中k₁、n的确定方法为：

1)取预测时刻点前100～200个特征数值；

2)将上述特征数值带入滚动轴承剩余寿命预测模型；

3)应用最小二乘法拟合特征数值及应力循环次数之间的非线性关系，即在满足最小均方误差的情况下，求取出上述状态方程中的k₁、n；此处申明，k₁、n是随数据的变化实时更新的。

5、设置模型初始参数数值；

设置预测步长为10，设置动态状态空间模型材料常数k₁、n初始值分别为2.0；1×10^-3，正确设置模型参数数值可以缩短运算时间，增加预测结果的准确性。

6、当特征数值大于等于报警阈值，将预测时刻点前100～200个特征数值带入滚动轴承寿命预测算法；应用最小二乘法更新初始模型参数，从而使模型符合滚动轴承实际运行状态。当特征数值大于等于分阶段点阈值，更新第三阶段模型参数，否则，更新第二阶段模型参数。

7、应用上述分阶段物理模型及粒子滤波算法进行滚动轴承剩余寿命预测。

预测过程中，当特征数值大于等于报警阈值小于分阶段点阈值时，应用参数更新后的第二阶段模型及粒子滤波算法进行滚动轴承运行状态评估；当特征数值大于等于分阶段点阈值小于停车阈值时，应用参数更新后的第三阶段模型及粒子滤波算法进行滚动轴承运行状态评估。

预测时每隔100分钟应用最小二乘法更新模型参数，修正模型轨迹，从而使模型符合滚动轴承实际运行状态。

8、当特征数值小于停车阈值时，预测循环进行，否则，预测循环停止，输出滚动轴承运行状态特征数值，并依据下述公式求取滚动轴承剩余寿命：

$L=(D-d)\×\frac{f_n}{f_s}\×t_{min}$ (19)

L指代剩余寿命

D指代预测的特征数值点数与输入的特征数值的点数的总和

d指代数据提供的输入的特征数值的点数

f_n指代采样点数

f_s指代采样频率

t_min指代输入的两个特征数值之间的时间间隔；

这个循环过程可以实时在线进行，即在轴承的工作过程中可以同步完成轴承的寿命预测。

Claims (1)

1.基于分阶段物理模型及粒子滤波的滚动轴承剩余寿命预测方法具体分为以下几个步骤：

(1)采集滚动轴承运行状态振动信号，从中提取特征数值：应用信号采集系统采集滚动轴承运行状态信息，将其输入特征数值公式，求取特征数值作为滚动轴承剩余寿命预测的输入值；

(2)设置报警阈值、分阶段点阈值、停车阈值分别作为预测初始点、第二阶段和第三阶段分界点、预测截止点：将运行状态振动区域按照标准分为A、B、C、D4个区域，设置报警阈值为B区域上限的1.25倍，设置分阶段点阈值为C区域上限数值，设置停车阈值为C区域上限的1.25倍；

(3)设置粒子滤波算法初始值；

1)取预测时刻点前100～200个数据点进行预测；

2)特征数值初始值为采集到的振动波形数据计算出的第一个特征数值；

3)按照公式z_t+1＝z_t+randn获取服从高斯分布的单个时间点的N个粒子，randn为服从高斯分布的0～1之内的随机数；z_t即为t时刻的粒子,z_t+1即为t+1时刻粒子；

4)设置粒子的初始权重均为

5)设置过程噪声w(t)及观测噪声v(t)的数值分别为0-1之间的数值；

(4)建立不同劣化阶段滚动轴承初始模型；

本次预测针对滚动轴承中晚期运行状态进行预测；

通过对Paris理论公式及Foreman理论公式进行改进，获取物理模型和状态空间模型相结合的分阶段滚动轴承剩余寿命预测模型，预测模型如下：

第二阶段：

状态方程： $\frac{V(t+1)}{S-S_{or}}=k_1\×q^n\×{\sqrt{\mathrm{\π}}}^n\×{(V_{dq}-V_{nr})}^n+w\left(t\right)$ (1)

观测方程： $y(t+1)=\frac{V(t+1)}{S-S_{or}}\×(S-S_{or})+v\left(t\right)$ (2)

第三阶段：

状态方程： $\frac{V(t+1)}{S-S_{or}}=(k_1\×q^n\×{\sqrt{\mathrm{\π}}}^n\×{\left(V_{dq}-V_{nr}\right)}^n+w(t\left)\right)\×\frac{x\left(t\right)}{V_{ed}-x\left(t\right)}---\left(3\right)$

观测方程： $y(t+1)=\frac{V(t+1)}{S-S_{or}}\×(S-S_{or})+v\left(t\right)$ (4)

其中，V_dq表示当前振动数值；V_nr表示正常运转振动数值；S_or指代开始劣化时应力循环次数；S指代应力循环次数；V(t+1)为轴承t+1时刻的特征数值大小；w(t)在t时刻状态转移的扰动，为轴承运转过程中产生的过程噪声向量数值；v(t)为t时刻轴承运转过程中产生的观测噪向量值；y(t+1)即为最终状态预测数值；k₁是第一材料系数，n是第二材料系数，二者依据材料特性和实验环境发生变化，是变量；q指代载荷；V_ed为振动停车阈值，V_max为目前采集到的振动数值的最大值，即为x(t)；其余参数同第二阶段参数；

(5)设置模型初始参数数值；

设置预测步长为10，设置动态状态空间模型材料常数k₁、n初始值分别为2.0；1×10^-3；

(6)当特征数值大于等于报警阈值，将预测时刻点前100～200个特征数值代入滚动轴承寿命预测算法；应用最小二乘法更新模型参数当特征数值大于等于分阶段点阈值，更新第三阶段模型参数，否则，更新第二阶段模型参数；

(7)应用上述分阶段物理模型及粒子滤波算法进行滚动轴承剩余寿命预测；

预测过程中，当特征数值大于等于报警阈值小于分阶段点阈值时，应用参数更新后的第二阶段模型及粒子滤波算法进行滚动轴承运行状态评估；当特征数值大于等于分阶段点阈值小于停车阈值时，应用参数更新后的第三阶段模型及粒子滤波算法进行滚动轴承运行状态评估；

预测时每隔100分钟应用最小二乘法更新模型参数，修正模型轨迹，从而使模型符合滚动轴承实际运行状态；

(8)当特征数值小于停车阈值时，预测循环进行，否则，预测循环停止，输出滚动轴承运行状态特征数值，并依据下述公式求取滚动轴承剩余寿命：

$L=(D-d)\×\frac{f_n}{f_s}\×t_{min}$ (5)

L指代剩余寿命

D指代预测的特征数值点数与输入的特征数值的点数的总和

d指代数据提供的输入的特征数值的点数

f_n指代采样点数

f_s指代采样频率

t_min指代输入的两个特征数值之间的时间间隔。