应用加权马氏距离与退化模型预测丝杠副寿命①

文娟

(浙江工业大学机械工程学院 杭州 310023)

(恒丰泰精密机械股份有限公司 温州 325000)

0 引言

作为数控机床的关键单元,滚珠丝杠副的性能关系到整个数控机床系统的可靠性与安全性。一旦滚珠丝杠副出现异常,则会降低整个数控机床的性能,影响加工质量,更甚者可能会造成机床停机,带来灾难。因此,适时监测滚珠丝杠副的性能状态,并根据状态监测结果指导维修计划对于提高整个制造系统可靠性与降低生产成本具有重要意义。状态维修主要包括故障诊断和故障预警2 个部分[1-4]。其中,故障预警的主要任务为判断设备的当前状态,并根据当前状态对其剩余服役寿命进行预测。准确的故障预警能够在确保安全的情况下节约维修费用,因而受到了越来越多研究者的关注。

寿命预测方法一般分为2 类:数据驱动方法和基于模型的方法。其中,数据驱动方法利用人工神经网络和支持向量机等机器学习方法建立设备寿命与状态监测数据间的关系。针对丝杠副的寿命预测,高宏力等人[5]提出了一种基于多模型融合技术和B 样条模糊神经网络的丝杠寿命预测技术。赵敏等人[6]利用多变量灰色模型建立丝杠寿命预测模型,实现丝杠的寿命预测。Li 等人[7]采用高斯过程回归预测滚珠丝杠副的剩余寿命,并研究了传感器富余与传感器不足时的不同寿命预测策略。此类方法需要大量历史失效数据来构建准确的寿命预测模型。然而在实际应用中,由于丝杠副具有较高的可靠性,因此很难获得足够的滚珠丝杠副性能退化数据。另一方面,基于模型的寿命预测方法试图建立能够反映设备退化过程的数学模型[8-9],并结合设备的实时状态监测信息更新模型参数,从而达到寿命预测的目的。与数据驱动方法不同,基于模型的方法需要的状态监测数据较少,且能充分利用设备的先验知识,因而更适用于滚珠丝杠副的寿命预测,但此类方法中存在着一些问题。

基于模型的寿命预测技术中的一个重要问题就是如何构建合适的健康指标,一个合适的健康指标能够简化退化建模并提高寿命预测的准确度。利用信号处理技术提取信号特征作为健康指标是一种非常广泛的健康指标构建方法。Yan 等人[10]利用时域特征反映轴承的健康状态。Kim 等人[11]提取信号的频域特征来识别轴承退化状态。Singleton 等人[12]应用Choi-Williams 分布对信号进行处理,提取时频域特征作为轴承的健康指标。以上特征均能从不同的角度对信号进行表述,但一个特定的特征一般只对一个特定退化阶段的特定损伤类型敏感。由于滚珠丝杠副存在多种损伤,一个单一的特征很难表征其整体退化状态。因此,应当利用多个特征的优势,有机结合多信息域特征,得到一个能够反映其退化趋势的健康指标。Nie 等人[13]提出一种利用马氏距离融合多个特征构建轴承健康指标的方法,但融合过程中,由于有些原始特征可能与退化过程相关性小,会对融合结果造成负作用。为了解决这一问题,本文提出了一种基于加权马氏距离(weighted Mahalanobis distance,WTMD)的特征构建方法。在特征融合前,对原始特征进行选择,并赋予其不同的权值,从而得到一种质量更高的健康指标——加权马氏距离,用于滚珠丝杠副的剩余寿命预测。

基于模型的寿命预测技术中的另一大问题就是如何构建合适的模型描述设备的退化过程。Peng和Tseng[14]提出了一种通用的线性模型,用于预测产品的平均失效时间。Gebraeel 等人[8]利用指数模型描述轴承的退化过程,并预测轴承的剩余寿命。司小胜[9]提出了一种基于非线性模型的寿命预测方法,并将其应用于电池的寿命预测中。这些研究中,指数模型是一类比较常用的模型,特别适用于指数型增长的退化模式,因而在寿命预测领域得到了较为广泛的应用[15-16]。基于以上研究,本文试图将指数模型应用到滚珠丝杠副剩余寿命预测中,利用指数模型描述丝杠副的退化过程;通过粒子滤波(particle filter,PF)算法结合健康指标与退化模型,估计其健康状态并更新模型参数,完成剩余寿命预测。

1 PF 算法

PF 是一种利用一群具有特定权值的粒子近似状态估计值的方法。与传统的卡尔曼滤波(Kalman filtering,KF)方法不同,PF 方法不受限于高斯假设,特别适用于非线性非高斯应用场景,已经在寿命预测领域得到了许多应用[17-19]。

1.1 贝叶斯估计

大部分的动态系统都可以由2 个方程描述,即状态转移方程和观测方程,具体如下:

其中,ft是上一时刻到当前时刻的状态转移函数,ht是当前状态xt的观测函数,nt是独立同分布过程噪声,vt是独立同分布测量噪声。

分别使用x0:t={x0,…,xt} 和y1:t={y1,…,yt} 表示当前所有状态和所有观测值,并假设系统的状态服从一阶马尔科夫过程,即p(xt|x0:t-1)=p(xt|xt-1)。利用贝叶斯理论,可以利用式(3)与(4)递归地计算时刻t的后验概率密度函数p(xt|y1:t)。

其中,p(xt|y1:t-1) 为时刻t的先验概率密度函数,p(xt|xt-1) 为由状态转移函数定义的转移密度,p(yt|xt) 为观测函数定义的似然,p(yt|y1:t-1) 为常数项。

1.2 PF 算法

实际应用中很难求得式(3)与(4)的解析解,而PF 算法就是一种近似求解后验概率密度函数的方法。利用PF 算法,系统状态的后验概率密度函数可以表示为

(1) 初始化:令t=0,从初始分布p(x0) 中采样得到N个粒子

(3) 权值更新:根据式(7)利用当前时刻t的测量值更新各粒子所对应的权值。

并根据式(8)对权值进行归一化处理。

(4) 重采样:为了增加有效粒子的数量,提高计算效率,需要移除权值极小的粒子,并复制权值较大的粒子。重采样后得到新的粒子集及其对应权值

(5) 状态估计:利用重采样得到的粒子及其权值估计当前时刻t的状态。

(6) 然后返回至步骤(2),在下一个测量时间点重复步骤(2)～(5)。

2 滚珠丝杠副剩余寿命预测方法

本文提出的基于模型的滚珠丝杠副剩余寿命预测方法如图1 所示。该方法包括2 部分,即健康指标构建与剩余寿命预测。

图1 提出方法流程图

为了验证本文提出方法,必须开展滚珠丝杠副性能退化实验,获取足够的状态监测信息。为了在较短周期内获得滚珠丝杠副的全寿命周期数据,在滚珠丝杠副加速性能退化实验台上开展滚珠丝杠副性能退化实验,并采集其退化过程中的状态监测信息,本文所采集的数据为振动信号。

在健康指标构建部分,首先采用多信息域信号分析方法分别对滚珠丝杠副的退化数据进行处理,提取多个原始特征。然后,评估每个原始特征对于滚珠丝杠副性能退化的敏感程度,从中选取较为敏感的特征。为了去除所选特征中的冗余信息,节约计算资源,采用相关聚类[21]对其进行聚类,选择每个类中的典型特征作为加权马氏距离方法的输入,得到健康指标WTMD。在寿命预测部分,采用指数模型描述滚珠丝杠副的退化过程,通过PF 算法结合实时监测信息与模型,利用实时WTMD 值更新滚珠丝杠副退化状态与模型参数。最后,利用当前时刻的状态估计值与参数估计值预测滚珠丝杠副之后的退化趋势,并预测其剩余寿命。第3 节与第4 节将分别对这2 部分进行详细叙述。

3 健康指标构建

3.1 特征选择

采用不同的信号分析方法提取滚珠丝杠副振动信号的原始特征,通常对设备退化过程敏感的原始特征应具有单调退化趋势,即与时间有较好的相关性[21]。为了衡量每个特征的单调性,计算特征与时间的斯皮尔曼系数作为该特征的趋势性指数,反映该特征的单调性[22]:

其中,{Tk}k=1:K为滚珠丝杠副全寿命周期中的信号采集时刻点,为从振动信号中提取的原始特征,K为退化过程中信号采集时刻点的总数,与分别为{Tk}k=1:K与的平均值。如果一个特征随时间单调递增或递减,则其趋势性指数值ρi为1;反之,如果一个特征保持不变或者随机变化,则ρi为0。因此,较高的趋势性指数值意味着较好的单调性。

由式(10)得出的特征趋势性指数值位于区间[0,1]中,为了选择单调性较好的特征,选择趋势性指数值大于0.5 的原始特征进行后续处理。然而这些原始特征中可能存在冗余信息,因为有些原始特征所代表的退化趋势十分相似。为了识别并移除这些冗余信息,采用相关性聚类算法来对这些特征进行分类[21],将相似特征归为一类。为了确定特征的类别数,使用一个衡量聚类算法有效性的指数——PBM 指数[23]评估不同类别数时聚类算法的有效性,选择PBM 指数最高的类别数作为特征的类别数进行分类。PBM 指数的定义如下:

确定特征的类别数为H*后,将所选择的特征分成H*类,同一类特征则认为其蕴含相似的退化信息,从每一类中选取趋势性指数最高的特征作为典型特征。

3.2 特征融合

为了充分利用每个特征的优势,结合典型特征中所包含的有用信息,采用加权马氏距离方法对得到的典型特征进行融合,计算滚珠丝杠副退化过程中振动信号特征与原始健康状态信号特征的WTMD,将其作为滚珠丝杠副的健康指标。

马氏距离(Mahalanobis distance,MD)计算是一种简单有效的特征融合方法[13],该方法计算给定特征向量与设备正常状态下特征向量集的MD 作为设备的健康指标。传统MD 计算方法对于每个特征的处理方式一样,从而每个特征对于MD 计算结果的影响相同。而在滚珠丝杠副的退化过程中,每个特征对于设备状态变化的敏感程度不同。为了更好地表征滚珠丝杠副的健康状态,敏感程度较高的特征应该在特征融合中发挥较大的作用。因此,本文提出一种基于加权马氏距离的特征融合方法,如下:

其中,zih=[zi1,zi2,…,ziH*] 为H*个典型特征;W=diag(w1,w2,…,wH*);wh是第h个典型特征的归一化权值,wh=,反映了该特征对于滚珠丝杠副退化的敏感程度;C为滚珠丝杠副正常状态下典型特征向量集的协方差矩阵;μ为其均值。式(12)中,敏感程度较高的特征具有较大的权值,对计算结果WTMD 的影响较大,因此,通过该方法得到的健康指标对滚珠丝杠副的性能退化更为敏感。

4 剩余寿命预测

4.1 退化模型

数控机床的工况复杂多变,在整个制造和操作过程中都存在较大的随机性。因此,认为滚珠丝杠副的退化过程是一个随机过程。考虑随机误差的影响,滚珠丝杠副的退化过程可以采用如下指数形式表述:

式中,ϕ表示滚珠丝杠副的正常状态;β是随机变量;ε(t)=σW(t) 是一个中心布朗运动,其均值为0、方差为σ2t,代表随机过程中的随机误差。通过两边取对数,该模型可以变换为线性形式:

其中,β服从正态分布N(μ1,)。

建立退化过程的模型后,则需要结合状态监测信息对滚珠丝杠副的状态与模型参数进行估计。假设时刻tk已有其健康指标的对数观测值L1:k={L(t1),L(t1),…,L(tk)},并已知参数β,则L1:k的联合分布可以表示为

式中,Δt=tj-tj-1是一个恒定的时间区间。当给定β的先验分布时,则可以由贝叶斯理论获得tk时刻β的后验分布p(β|L1:k)。

4.2 剩余寿命预测

得到tk时刻的状态和参数估计值后,接着需要确定剩余寿命概率密度函数,即确定健康指标达到预定阈值γ所需时间的概率分布。根据式(14),可以将模型变换为

令uk=t-tk,则上式变为

式中,V(uk)=B(uk+tk)-B(tk) 是一个标准的布朗运动。利用健康状态估计与参数更新结果,可以通过式(18)得到tk时刻的剩余寿命概率密度函数。

5 案例分析

为了阐述提出方法在实际中的应用并证明其有效性,本研究开展了滚珠丝杠副加速性能退化实验,采用实验中采集的退化数据对提出方法进行验证。

5.1 实验数据

正常条件下,滚珠丝杠副退化需要很长的时间,因此,很难获取其全寿命周期数据。为了解决这一问题,在如图2 所示滚珠丝杠副加速性能退化实验台上开展加速寿命实验,通过增加负载,加速退化速度,在较短时间内使滚珠丝杠副失效,并采集其退化过程中的状态监测信息,用于剩余寿命预测研究。

图2 滚珠丝杠副加速性能退化实验台

滚珠丝杠副的工作转速一般低于3000 r/min,即50 Hz,结合实际中关于机床振动的经验,实验中,设定滚珠丝杠副振动信号的采样频率为5 kHz。在开展全寿命周期实验前,对滚珠丝杠副进行了加速磨损实验,设置磨损过程中滚珠丝杠副的轴向负载为1563 N,转速为300 r/min,并持续采集其振动数据。分析振动信号得出其磨损速度较慢,短时间内信号的变化不大。通过分析信号的变化情况,并结合关于滚珠丝杠副退化的实际经验,确定每隔30 min 采集一次信号较为合适。因此,在全寿命周期实验中,每隔30 min 采集一次振动信号。

由于靠近电机端轴承上安装的加速度传感器所采集到的振动信号能够较好地反映滚珠丝杠副的状态变化,因此,本文选用该传感器采集到的振动信号进行分析,构建健康指标,进行寿命预测。当电机端轴承振动信号超过10 g 时,出现了明显的异常噪声。为了避免滚珠丝杠副损伤蔓延至整个测试平台,引发安全问题,因此在振动信号超过10 g 时停止实验,认为滚珠丝杠副失效。磨损退化过程中,共采集了145 次信号后滚珠丝杠副失效。由于经济与时间成本的限制,仅采集了一组滚珠丝杠副退化过程的监测数据,用于验证文中提出方法。

图3 为滚珠丝杠副整个退化过程中靠近电机端轴承上安装的加速度传感器所采集到的振动信号,其横坐标为振动信号的每一个采样点。由图可知,最初,滚珠丝杠副工作较为稳定,振动信号的变化较小。随着时间的增加,其磨损增加,振动信号出现小幅度的变化。进入末期后,滚珠丝杠副退化加快,振动信号迅速增大,直至失效。

图3 滚珠丝杠副振动信号

5.2 滚珠丝杠副健康指标构建

本文选用10 个时域特征,28 个时频域特征以及2 个三角函数特征[24]作为原始特征,如表1 所示。时频域特征的提取过程中,采用小波包变换对信号进行3 层分解,得到8 个小波包系数,将各小波包系数的能量及其与总能量的比值作为特征。此外,采用经验模态分解算法对信号进行分解,将前6个本征模态函数(intrinsic mode function,IMF)的瞬时能量及其与总瞬时能量的比值作为特征。三角函数特征的提取过程中,先采用三角函数对原始信号进行变换,然后提取变换后信号的特征,将其作为原始信号的特征[24]。本研究分别采用反双曲线正弦函数与反正切函数对原始信号进行变换,再分别提取其标准差(standard deviation,SD)作为特征,具体计算公式如下:

表1 滚珠丝杠副振动信号原始特征

式中,X=为采集到的数据序列,xn为每一个采样点的信号幅值,N为数据序列中采样点的总数。

采用式(10)计算40 个原始特征的趋势性指数值,计算结果如图4 所示。选择趋势性指数大于0.5的特征,共选出22 个特征。将所选择的特征输入相关聚类算法中,计算类别不同时的PBM 指数,如图5所示。由图可知,当H=3 时,PBM 指数值最大,因此将这22 个特征分为3 类。从每个类中选取趋势性指数值最大的特征作为这一类的典型特征。为了便于观察分析每个典型特征的特点,对选取出来的3 个典型特征分别进行归一化处理,得到结果如图6 所示,由图可知,这3 个特征展现了整个寿命周期中3种不同的退化趋势,各自蕴含了不同的退化信息。Y5 在滚珠丝杠副工作初期变化较大,但关于中后期的退化情况信息量较少;Y12 描述了寿命中期的波动情况;Y32 直到寿命末期才表现出一定的退化趋势。

图5 不同类别个数的PBM 指数

图6 滚珠丝杠副的典型特征

为了充分利用每个典型特征的优势,采用加权马氏距离方法对3 个典型特征进行融合,计算退化过程中滚珠丝杠副振动信号典型特征与原始健康状态下典型特征集的WTMD,作为其健康指标。为了说明加权马氏距离方法的有效性,对于40 个原始特征,采用马氏距离方法计算其健康指标MD。对2个指标分别进行归一化,如图7 所示。由图可知,WTMD 的退化趋势更明显,蕴含更加丰富的性能退化信息,更适用于退化建模。

图7 滚珠丝杠副退化过程中的WTMD 值与MD 值

为了量化地比较2 种健康指标的优劣,采用式(10)分别计算WTMD 与MD 的趋势性指数,结果如表2 所示。由表可知,WTMD 的趋势性指数高于MD 的趋势性指数。因此,WTMD 对滚珠丝杠副的损伤发展更为敏感,因而更加适用于其剩余寿命预测。

表2 两种不同健康指标的趋势性指数值

5.3 滚珠丝杠副寿命预测

将加权马氏距离方法得到的健康指标WTMD输入状态空间模型中,利用PF 算法估计滚珠丝杠副的健康状态与更新模型参数,并预测其剩余寿命。利用该方法得到的滚珠丝杠副退化状态估计结果如图8 所示。由图可知,提出方法能够比较准确地估计WTMD 值,追踪滚珠丝杠副的退化路径。

图8 滚珠丝杠副退化状态估计结果

为了说明指数模型在滚珠丝杠副剩余寿命预测中的有效性,采用线性维纳过程模型[14]、非线性维纳过程模型[9]与其对比,图9 为3 种模型的滚珠丝杠副剩余寿命预测结果。由图可知,对于指数模型而言,在状态更新刚开始的时候,由于观测到的WTMD 值较少,其剩余寿命预测结果存在较大误差,准确度低。但随着观测数据的增加,其预测误差减小,在利用了83 组数据进行健康状态估计与参数更新后,指数模型的结果收敛于真实值,但其他2 种模型的预测结果依然与真实值有一定的误差。对比3 种模型的预测结果,线性模型的预测误差最大,这是因为滚珠丝杠副的损伤并不以线性形式发展。非线性模型的预测准确率高于线性模型,但该模型的收敛速度比指数模型慢,直到退化过程的末期才能比较准确地预测滚珠丝杠副的剩余寿命。图10 为不同模型在时间t=62.5 h 时预测的失效时间概率密度函数曲线,由图可知,指数模型的分布中心最接近真实失效时间,但其概率密度分布宽度比线型模型宽,而非线性模型的分布最宽,即误差范围最大。

图9 滚珠丝杠副剩余寿命预测结果

图10 不同模型预测的失效时间概率密度函数

为了更加直观地比较以上3 种模型的寿命预测效果,采用以下公式计算时刻t的预测误差:

利用式(20)分别计算3 种模型在各个信号采集时刻点的预测误差。为了便于比较,将每种模型的预测误差值分成3 个区间,并计算处于每个区间的信号采集时刻点占总时刻点数的百分比,如表3所示。由表可知,指数模型预测结果中,52.08%的时刻点的预测误差不超过30%,而非线性模型预测结果中,预测误差不超过50%的时刻点仅有30.56%,线性模型预测结果中,仅有4.86%的时刻点的预测误差不超过50%。因此,指数模型的预测误差最小。这是因为滚珠丝杠副的损伤整体上呈指数形式增长,所以在这3 种模型中,指数模型最适合于描述其退化过程,并具有最优的寿命预测性能。

表3 三种模型的预测误差比较

为了说明粒子滤波的有效性,使用卡尔曼滤波算法进行对比。图11为采用指数模型描述退化过程时,粒子滤波与卡尔曼滤波的寿命预测结果对比。由图可知,在预测开始阶段,卡尔曼滤波的效果要好于粒子滤波,但随着退化的加深,在滚珠丝杠副寿命后期,粒子滤波的预测结果较卡尔曼滤波更接近真实值。对于设备寿命预测而言,为避免其失效带来的安全隐患与经济损失,其寿命后期的预测准确率更为重要。因此,粒子滤波的预测性能更佳,更适用于滚珠丝杠副的寿命预测。

图11 不同滤波算法的寿命预测结果

本文所建立的模型能够用于实际安装于机床的滚珠丝杠副的寿命预测。应用过程中,同样采用加速度传感器采集滚珠丝杠副的状态监测信息,然后利用采集到的信号更新模型参数,得到实际工作状态下的模型参数估计结果,进而预测其寿命。本文提出的滚珠丝杠副寿命预测方法目前已经应用到某企业的数控机床状态监测系统中,用于估计实际工作中的滚珠丝杠副的剩余寿命。

6 结论

本文以滚珠丝杠副为对象,针对其工作过程中的性能退化问题,基于实验中采集的全寿命周期数据,利用加权马氏距离方法与指数模型估计其健康状态并预测其剩余寿命。

为了说明提出方法的优势,将所提出方法的实验结果与对照方法的实验结果进行对比。比较马氏距离方法和加权马氏距离方法得到的健康指标,结果表明WTMD 的趋势性指数值高于MD,对滚珠丝杠副的损伤发展更为敏感,更适用于其剩余寿命预测。同时,指数模型能够很好地描述滚珠丝杠副的退化过程,利用该模型能够有效地估计设备健康状态并预测其剩余寿命,性能优于线性维纳过程模型与非线性维纳过程模型。此外,由于提出方法对于状态监测数据数量的要求低于数据驱动方法,因而更适用于实际应用。

由于本文采集的数据量有限,对于提出方法的验证存在一定的局限性,后续研究中将采集更多的样本以进一步验证提出方法的泛化性。