基于猎人猎物优化算法改进粒子滤波的滚动轴承剩余使用寿命预测技术

张田雨 王庆锋* 舒 悦 肖 旺

(1.北京化工大学机电工程学院,北京 100029;2.北京化工大学高端机械装备健康监控及自愈化北京市重点实验室,北京 100029;3.合肥通用机械研究院有限公司, 合肥 230031;4.国家管网集团联合管道有限责任公司西部分公司,乌鲁木齐 830013)

引 言

滚动轴承广泛应用于旋转类机械设备中,是关键易损件之一。 对智能传感器所采集的滚动轴承实时振动数据的分析,可为预防性维修决策提供依据,延长滚动轴承的使用寿命,大大提高生产效率,因此关于滚动轴承剩余使用寿命(RUL)预测的研究近年来越来越受到人们的重视[1-2]。

目前关于剩余使用寿命预测的方法可归纳为两种[3]:第一种是基于数据驱动的方法,第二种是基于模型驱动的方法。 基于数据驱动的方法是在获取一定数据量的基础上,运用统计分析或机器学习方法对剩余使用寿命进行预测。 该种方法对数据量的大小有要求,所以往往很难符合实际工程化应用需求,存在一定的局限性。 基于模型驱动的方法则是为描述系统的衰退趋势建立数学模型,并根据观测数据对数学模型的参数进行实时更新,进而对剩余使用寿命进行预测。 相较于基于数据驱动的方法,基于模型驱动的方法可以充分利用已知数据包含的故障信息,并得出相对稳定的剩余使用寿命预测结果。

常见的基于模型驱动的方法包括隐马尔可夫模型、卡尔曼滤波、粒子滤波等。 其中粒子滤波预测方法基于贝叶斯滤波算法提出,可改善其他方法在处理长期预测时的不确定性问题,现已广泛应用于剩余使用寿命预测领域。 Li 等[4]提出了基于粒子滤波的预测方法,使用相关矩阵聚类和加权算法计算健康指标,并运用粒子滤波算法预测健康指标的变化趋势,最后得到轴承的剩余使用寿命。 祝志远等[5]提出将粒子滤波算法应用于疲劳裂纹扩展的参数估计和剩余寿命预测上,预测得到了剩余寿命中值和置信区间。 Lei 等[6]提出一种基于模型的滚动轴承剩余使用寿命预测方法,该方法包括健康指标构建和RUL 预测两个模块,在第一个部分构建了一个名为加权最小量化误差的新健康指标,在第二个部分使用粒子滤波预测RUL,减小了预测中的累计误差。 这些采用粒子滤波算法预测剩余寿命的研究可获得更加准确的预测结果,充分体现了粒子滤波算法在解决长期不确定性预测方面的优势,但其本身仍存在多次迭代后粒子权重退化问题。

为解决上述粒子滤波算法存在的缺陷,余臻等[7]提出使用无迹粒子滤波方法预测航天发动机的排气温度,获得了比粒子滤波算法更优良的效果。许雨晨等[8]提出粒子群优化算法改进的粒子滤波预测算法,可以较为准确地预测出滚动轴承的剩余使用寿命。 贺宁等[9]提出使用天牛须搜索算法优化粒子滤波的重采样过程,解决了粒子多样性丧失问题,并将该方法应用于电池的剩余使用寿命预测,获得了接近真实电池寿命的结果,预测中的累计误差明显减小。 徐仁义等[10]提出一种基于均方谐噪比的指标,之后运用改进正则化的粒子滤波算法优化重采样过程,该算法利用基于欧式距离的核函数改进了粒子滤波算法中的重采样过程,在预测滚动轴承的剩余使用寿命中获得了更接近真实剩余使用寿命的结果。

此外,在预测过程中,还会存在特征指标波动性大、无法确定合适的置信区间进而导致预测精度下降的问题,现在通常采用设置95%或99%置信区间的手段,但该方法主观选择的程度较大,无法适用于全部的特征指标,泛化性较差。 为解决上述问题,本文提出一种融合Hodrick-Prescott(HP)趋势滤波-边界线(HPTF-BL)、猎人猎物优化算法(HPO)改进粒子滤波(PF)的滚动轴承剩余使用寿命预测方法。以振动信号为例,首先提取表达设备退化趋势的特征,然后通过HPTF-BL 对特征进行处理,得到特征的上下退化边界线与主要退化趋势,解决了较难选择合适的置信区间且选择主观性大的问题,并减少了预测误差;最后使用HPO-PF 算法对滚动轴承的剩余使用寿命进行预测,解决了粒子滤波算法中粒子权重退化问题,同时提高了预测精度。

1 滚动轴承剩余使用寿命预测模型

1.1 基于小波信息熵的特征提取方法

何正嘉等[11]提出用归一化的特征指标来评判机械设备的健康状态,例如相关系数法、凝聚函数法、小波信息熵法等,这些方法计算出的特征指标均在0 ～1 之间。 由于小波信息熵法使用了熵值这一复杂性衡量指标,该指标可反映出机械设备更多的故障信息,因此对于一台包含滚动轴承的机械设备,适合使用小波信息熵法计算特征指标,来评判其运行状态。

设备运行中产生振动信号,通过小波包分解l次并对每个频带重构得到2l个分解信号xli(k),第i个分解频带信号的能量Eli和相对能量分别为

式中,i=1,2,…,2l;k=1,2,…,n,n∈Z。 相对能量总和

小波信息熵Ent定义为

式中,对数的底取2l,则Ent∈[0,1]。 设备运行状态数增加,不确定性增强,小波信息熵值必然增大,则该时刻旋转机械的运行状态越差;小波信息熵值越小,说明该时刻旋转机械的运行状态越好[12]。

由于滚动轴承的性能退化具有单调不可逆性,可将单调性作为特征指标与滚动轴承性能退化一致性的评判标准[13]。 当特征指标随时间单调递增或者递减时,单调性为1;当特征指标随机波动时,单调性为0。

单调性计算公式如下。

式中,X是时间序列,ε(x)为单位阶跃函数,x≥0时,ε(x)为1;x＜0 时,n为总样本数。

采用在美国辛辛那提大学智能维护系统(IMS)中心实验平台上采集的数据集2(IMS2)的1 号轴承数据[14]作为验证数据,计算当下常用的4 种用于滚动轴承剩余使用寿命预测的特征指标的单调性,4种指标分别为有效值(标签为1)、整流平均值(标签为2)、峭度值(标签为3)和谱距离指标(标签为4),并计算基于小波信息熵(标签为5)的特征指标的单调性,对比结果如图1 所示。

图1 几种特征指标的单调性对比Fig.1 Monotonicity comparison of several indexes

由图1 的对比可知,小波信息熵指标比其他4种特征指标具有更好的单调性,适用于作为后续滚动轴承剩余使用寿命预测的指标。

1.2 基于HPTF-BL 的特征处理方法

1.2.1 HP 滤波

HP 趋势滤波由Hodrick 等[15]提出,是经济学中常用的数据分析方法,可将数据分解为长期趋势项和短期波动项。 由于HP 趋势滤波具有降低数据噪声影响的作用,现常被用于提取各种时间序列的趋势,在数据预测及产品性能退化的可靠性分析等领域应用广泛。

HP 趋势滤波是一种高通滤波器,利用了最小二乘损失函数,采用l2 范数对二次差分矩阵进行计算,可将时间序列xt分解为长期趋势项xk和具有随机波动特性的波动项xc。 其中,趋势项xk被定义为式(5)的最小化函数解。

式中,i为时间序列数据的序号,n为时间序列的样本个数,右侧第一项表示趋势项对原序列的跟踪程度,第二项表示趋势项的光滑程度,λ为平滑参数,用于控制趋势项的平滑程度。 对式(5)求xki序列的一阶偏导,即可获得趋势序列xk。

式中G为系数矩阵,I为单位矩阵。 当λ→0 时,趋势项对序列的跟踪程度达到最大;当λ→+∞时,趋势项序列光滑程度达到最大;当λ=0 时,HP 滤波方法即退化为最小二乘法。 通过上述HP 滤波方法即可将时间序列分解为周期项和波动项。

1.2.2 HPTF-BL 处理方法

本文提出一种针对波动性大的健康指标的处理方法—HP 趋势滤波-边界线(HPTF-BL)方法,可以很好地解决健康指标波动性大、置信区间难以确定的问题,减少最后预测结果的误差。

HPTF-BL 方法是在HP 滤波的基础上对其公式进行形式改写的一种方法,其主要步骤如下:

1)窗口划分,即对获得的已知时间序列xt进行分割,获得每组包含m个数据的多个窗口;

2)对每个窗口内的数据大小进行排序,获得该窗口内的最大值与最小值;

3)用最大值代替窗口内所有的数据,得到上边界xup,用最小值代替窗口内所有数据,得到下边界xlow;

4)用HP 滤波对时间序列xt、步骤3)中得到的上边界xup和下边界xlow进行处理,得到主要趋势项xk,xk的计算公式如式(5)所示。

处理后的健康指标HI如式(7)所示。

1.3 改进粒子滤波算法

1.3.1 粒子滤波算法

粒子滤波是在传统非线性滤波方法如卡尔曼滤波、扩展卡尔曼滤波的基础上发展而来的,在非线性、非高斯系统的模型参数估计中表现出明显的优势,并且已被应用于寿命预测领域[16]。

假设系统在离散时间序列tk=kΔt的状态可以用式(8)的状态传递函数描述。

式中,xk是系统在tk时刻的状态,fk是系统状态传递函数,θk是模型参数向量,ωk是系统噪声。 系统状态值与观测值之间的关系为

式中,zk是系统在tk时刻的观测值,hk是系统的观测函数,vk是观测噪声。 首先根据k-1 时刻的状态对k时刻的先验概率密度函数进行预测

当测得新的观测值后,对状态概率密度函数进行更新,得到k时刻的后验概率密度函数

用重要密度函数的离散采样点和对应权值来描述p(xk/z1:k),则式(11)转化为

式中,δ(·)表示离散采样点,重要性权值wk采用式(13)进行更新。

将式(13)条件概率增加参数项θk

利用新观测值zk实现对系统状态xk和模型参数θk所对应权值的不断更新。

1.3.2 猎人猎物优化算法

猎人猎物优化(hunter-prey optimizer, HPO)算法[17]通过模拟动物猎食过程,在搜索空间中按照一定的规则和策略对种群进行引导和控制,不断更新群体中每个猎人或猎物的位置,并用目标函数评估出新的位置,对一个问题进行寻优。 该算法具有收敛快、寻优能力强的特点。

首先猎人或猎物在搜索范围内按照式(15)随机初始化位置。

式中,xi是猎人或猎物的位置,lb是问题变量的最小值,ub是问题变量的最大值,d是问题变量的维度。式(16)定义了搜索空间的下界和上界。

生成初始种群后,使用目标函数计算每个解的适应度值。 对于猎人的搜索机制,式(17)给出了其数学模型

式中,x(t)是当前猎人位置,x(t+1)是猎人的下一次迭代位置,Ppos是猎物的位置,μ是所有位置的平均值,Z是由式(18)计算出的自适应参数

其中,R1和R3是[0,1]内的随机向量,P是R1＜C的索引值,R2是[0,1]内的随机数,IDX是满足条件(P= =0)的向量R1的索引值,C是探索和开发之间的平衡参数,其值在迭代过程中从1 减小到0.02,具体计算如下。

其中,it是当前迭代次数,Max是最大迭代次数。 计算猎物的位置Ppos,再根据式(20)计算所有位置的平均值μ,然后计算与该平均位置的距离。

距离位置平均值最大的位置被视为猎物位置Ppos。 假设最佳安全位置是最佳全局位置,这将使猎物有更好的生存机会,猎人可能会选择另一个猎物,由式(21)更新猎物位置

式中,x(t)是猎物当前的位置,x(t+1)是猎物的下一次迭代位置,Tpos是全局最优位置,Z是由式(18)计算出的自适应参数,R4是[-1,1]内的随机数;cos 函数及其输入参数允许下一个猎物位置在不同半径和角度的全局最优位置,并提高开发阶段的性能。

为了选择猎人和猎物,结合式(17)和式(21),R5是[0,1]内的随机数,β为调节参数,如果R5值小于β,搜索到的位置被视为猎人,下一个位置将用式(17)更新;如果R5值大于β,搜索到的位置被视为猎物,下一个位置将用式(21)更新。

1.3.3 猎人猎物优化算法改进粒子滤波

利用猎人猎物优化算法改进粒子滤波的重采样过程,克服了粒子经过多次迭代后权重退化的问题。将粒子滤波中的粒子先验状态作为猎人猎物初始种群个体位置,利用迭代寻优过程改善粒子的分布情况,将退化的粒子集优化至高似然值,使大部分粒子都能集中在真实状态附近,解决了常规粒子滤波算法中重采样使用的函数是次优的问题。

HPO-PF 算法的实现步骤如下。

1)设置HPO-PF 算法的参数,包括粒子个数、种群规模、最大迭代次数。

2)初始化粒子集。 从重要性概率密度函数中随机抽取N个粒子。

3)更新粒子位置计算各个粒子的适应度值,将粒子集作为猎人或猎物,更新当前最优位置。 定义适应度函数为

式中,zR当前状态的实际量测值,zP为预测值。

4)种群移动,更新猎人或猎物的位置。 根据式(17)或式(21)计算猎人或猎物当前的位置。

5)判断循环是否停止。 当满足最大迭代次数或达到最优适应度值时,循环结束,否则重复步骤4)。

6)根据式(14)计算各个粒子权重,最后得出估计状态。

1.4 预测模型构建与参数选择

1.4.1 模型构建

本文提出的滚动轴承剩余使用寿命预测方法流程图见图2,方法分为两个阶段:构建退化预测指标以及预测剩余寿命。

图2 剩余使用寿命预测模型Fig.2 Predictive model of RUL

本文使用Paris-Erdogan 裂纹扩展模型[18-19]作为状态空间模型,该状态空间模型能很好地描述轴承的退化过程,并且易于构建,可满足工程上关于预测的可操作性的要求。 其中故障尺寸增长率可表示为

式中,x为故障尺寸;N为材料的疲劳寿命;C0和m为与材料相关的常数;ΔK为应力强度因子,计算式为

式中,β为与材料相关的参数。

将式(24)代入式(23)可得到

假设量测值zk和状态值xk存在线性关系,改写式(25)成如下形式。

式中,A和n为需要根据观测值进行估计的未知参数,vk为观测噪声,ωk为系统噪声,式(26) 与式(27)即为状态空间模型表达式。

预测滚动轴承剩余使用寿命的具体步骤如下。

1)构建退化预测指标

①获取滚动轴承的原始振动信号,计算输入的原始数据的信息熵指标。

②使用1.2 节的HPTF-BL 特征处理方法对计算出的信息熵指标进行处理,得到退化预测指标。

2)预测剩余使用寿命

①使用最小二乘拟合初始化状态空间模型的参数A和n,运用HPO-PF 算法对参数进行更新迭代,轴承失效退化曲线由算法递推预测得到。

②将预测出的失效退化曲线结合预先设定好的失效阈值曲线,得到滚动轴承在该时刻预测的最长寿命、中间寿命以及最短寿命,分别对应预测值上限、参考预测值与预测值下限[20]。

1.4.2 预测结果评估

使用均方根误差(RMSE)以及拟合优度R2作为预测性能的评估指标,来验证预测模型的有效性。均方根误差RMSE 定义如下。

用相关系数(即拟合优度)R2来评估预测结果与滚动轴承的真实退化曲线的拟合程度。R2的值在0 ～1 之间,其值越接近于1,表明预测出的退化曲线的拟合效果越好,预测模型的准确度越高。R2定义如下。

1.4.3 参数选择

在通过HPTF-BL 方法处理特征指标的过程中,需要对λ和m这两个参数进行选择。λ是一个控制时间序列平滑性的指标,在进行特征指标处理时需要获取特征指标的主要退化趋势,根据文献[21],λ选择为16。m表示窗口内数据的多少,至少可选择为3,至多可选择为与已知序列内数据量相同的值,m的值越大,得到的结果越具有代表性,但往往会忽略一些数据信息;m的值越小,得到的结果则更能反映出特征指标的边界信息,具有更好的准确性,本文中m选择为5。

在改进的粒子滤波算法中,本文设置HPO 算法的种群大小为30,最大迭代次数为500。 粒子滤波算法中粒子数的取值对预测结果有较大影响,使用IMS2 的1 号轴承数据[14]作为验证数据,讨论不同粒子数对预测结果的影响。 选用已知数据量为367的数据,使用上述两种评估指标对预测结果进行评估,结果如表1 所示。

表1 不同粒子数对预测结果的影响Table 1 Effect of different particle numbers on the prediction results

通过表1 可知,当粒子数为500 时,粒子滤波可获得最优的预测效果,故粒子数设置为500。

2 实例验证

2.1 实验数据

选用美国辛辛那提大学IMS 中心实验平台上采集的两组数据[14]对本文提出的方法进行验证。两组数据分别如下:数据集1(IMS1)中的3 号轴承运行至失效的数据,该组数据总共包含2 156 个数据文件,3 号轴承在失效实验结束时出现内圈故障;IMS2 中的1 号轴承运行至失效的数据,该组数据共包含984 个数据文件,1 号轴承在失效实验结束时出现外圈故障。 每一个数据文件均对应一个文件序号。

2.2 数据验证及分析

首先,对IMS1 中3 号轴承和IMS2 中1 号轴承x轴方向的原始振动信号进行特征提取,计算出信息熵指标,如图3 所示。

图3 信息熵指标Fig.3 Information entropy index

由文献[22]可知在IMS1 中的数据文件序号为1 608、IMS2 中的数据文件序号为533 时,轴承发生性能退化。 工程应用场景下,一旦滚动轴承发生早期故障,之后得到的数据便可作为剩余使用寿命预测的已知训练数据,故本文分别选取IMS1 在数据文件序号为1 608、IMS2 在数据文件序号为533 之后的轴承数据的信息熵指标,采用HPTF-BL 方法进行处理,处理后的结果如图4 所示。 可以看出一维的健康指标被分为上下边界以及信息熵主要趋势,数据长度与原健康指标的数据长度保持一致,并将原本的健康指标包含其中,消除了原本健康指标存在的局部波动,解决了随机噪声对预测结果的影响,具有单调性好、波动小的特点。 得到的3 条曲线在之后的步骤中同时进行预测,可得到基于信号自身的预测置信区间,增强了在预测趋势步骤中置信区间设置的可解释性,进一步提高了RUL 预测结论的准确性。

图4 处理后的信息熵指标Fig.4 Information entropy index after treatment

当IMS1 的已知数据分别为数据文件序号1 608至数据文件序号1 830、1 930、2 030、2 130 时,使用HPO-PF 算法对已知数据进行训练,迭代预测得到失效退化曲线。 根据文献[11]设定从滚动轴承发生性能退化开始,小波信息熵指标变化值超过0.2时,即为达到运行状态不满意阶段,IMS1 轴承在信息熵指标为0.75 时开始退化,故IMS1 轴承的失效阈值设定为0.95。 IMS1 的3 号轴承在已知数据量不同情况下的预测结果如图5 所示。

图5 IMS1 的3 号轴承预测结果Fig.5 Predicted results for bearing No.3 in data set 1

当IMS2 的已知数据分别为数据文件序号533至数据文件序号750、800、850、900 时,使用HPOPF 算法对已知数据进行训练,迭代预测得到失效退化曲线。 IMS2 轴承在信息熵指标达到0.08 时开始退化,故相应的失效阈值设定为0.3。 IMS2 的1 号轴承在已知数据量不同情况下的预测结果如图6所示。

由图5、6 预测结果可知,随着已知数据文件的增加,HPO-PF 算法预测出的趋势越来越接近真实的信息熵指标,最终该滚动轴承的剩余使用寿命L可用式(30)计算得出。

式中,Np为预测出的轴承失效时刻对应的数据文件序号;NF为开始预测时刻对应的数据文件序号;tper为采样时间,在本文案例中为10 min。 最后可根据式(30)计算出在已知数据量大小不同时轴承的剩余使用寿命,具体结果如表2、3 所示。

由表2 和图5 可知,当IMS1 的已知数据为数据文件序号1 608 ～1 830 时,代入模型训练的数据有向上的趋势,所以预测得出的结果超前真实寿命值800 min 左右;但当已知数据为数据文件序号1 608 ～2 030 时,HPO-PF 算法的预测值上限为2 780 min、参考预测值为2 700 min、预测值下限为2 620 min;当已知数据为数据文件序号1 608 ～2 130时,预测结果与真实剩余寿命仅相差400 min 左右,越来越接近于真实寿命。 所得出的预测值上限、参考预测值和预测值下限可共同为企业预测维修时间点提供建议。

由表3 和图6 可知,当IMS2 的已知数据为数据文件序号533 ～750 和533 ～800 时,代入模型训练的数据略有向上的趋势,但总体上看还是趋于平缓,所以预测得出的结果延迟真实寿命值都在4 000 min 以上,不超过5 000 min;当已知数据为数据文件序号533 ～850 时,HPO-PF 算法的预测结果开始趋近于真实剩余寿命,只相差400 min 左右;当已知数据为数据文件序号533 ～900 时,预测结果与真实剩余寿命仅相差100 min 左右,越来越接近真实寿命。 由以上结果可知,本文预测模型得出的结论可根据已知数据量的变化实时动态改变。

2.3 性能对比

为了进一步验证本文提出的剩余寿命预测方法的优越性,选取5 种方法在IMS2 上进行比较和分析,即使用指数退化模型和Paris-Erdogan 模型作为状态空间函数的粒子滤波算法(PF-E 算法和PF-P算法)、使用指数退化模型和Paris-Erdogan 模型作为状态空间函数的粒子群算法优化粒子滤波算法(PSO-PF-E 算法和PSO-PF-P 算法)以及使用指数退化模型作为状态空间函数的猎人猎物算法优化粒子滤波算法(HPO-PF-E 算法)。 以上5 种方法的基本参数都与本文所提方法保持一致,采用已知数据量为367 的组别进行验证,只使用信息熵指标的主要趋势进行预测,最后通过RMSE 和R2评价预测效果,得到的对比结果如表4 所示。

表4 本文所提方法与其他5 种方法的对比结果Table 4 Comparison between proposed method and five other methods

对比结果表明,使用Paris-Erdogan 模型作为状态空间函数的方法整体上的预测效果都优于使用指数退化模型作为状态空间函数的方法;而在使用Paris-Erdogan 模型的方法中,本文所提方法不论从均方根误差上还是拟合优度(相关系数)上的表现都优于其他方法的预测效果。

3 结论

考虑到实际工厂的剩余使用寿命预测需求,本文提出了一种融合HP 趋势滤波-边界线(HPTFBL)、猎人猎物优化算法改进粒子滤波(HPO-PF)的滚动轴承剩余使用寿命预测方法,该方法从健康指标上下边界刻画、剩余使用寿命预测方法两个方面进行了优化。 实验与对比分析的结果表明,该方法针对不同时刻以前的历史时间序列数据,可预测获得不同时刻下滚动轴承剩余使用寿命的预测值上限、参考预测值和预测值下限,具有预测累计误差小、预测精度高的特点,可为企业工厂的预测性维修决策提供依据。