基于EMD多尺度威布尔分布与HMM的轴承性能退化评估方法

2022-02-22周雁翔

振动与冲击 2022年3期

王冉, 周雁翔, 胡雄, 陈进

(1.上海海事大学物流工程学院，上海 201306；2.上海交通大学机械系统与振动国家重点实验室，上海 200240)

滚动轴承是旋转机械中一种不可或缺的重要部件，对滚动轴承采用基于状态的预测性维护策略是提高旋转设备安全性与可靠性的重要手段[1-2]。滚动轴承从正常运行到最终失效往往要经历一系列性能退化的过程。与故障诊断技术不同，性能退化评估更侧重于动态地揭示设备在运行过程中健康状态退化的规律，因此更有利于开展预测性维护，避免机械设备严重故障的发生。

滚动轴承性能退化评估的主要步骤包括性能退化特征提取以及评估模型的建立，其核心目标在于构建一种能够有效反映轴承性能退化趋势的性能指标，且该指标需具有良好的单调性、敏感性与鲁棒性的性能指标[3]。常用的性能退化特征包括时域统计特征、频域特征、时频域特征等[4]。在评估模型方面，常用的评估模型包括支持向量机、模糊C均值[5-6]、支持向量数据描述[7]、隐马尔可夫模型(hidden Markov model，HMM)[8-10]等。其中，HMM的双随机链结构能够较为准确地描述轴承在退化过程中隐含的衰退状态与观测变量之间的关系，因此被广泛应用于轴承的性能退化评估中。例如：Jiang等[11]将冗余属性投影与HMM相结合建立了性能退化评估模型。Yu[12]建立了一种基于自适应HMM的在线学习框架用于量化评估轴承的性能退化状态。这些方法验证了HMM在轴承性能退化评估上的有效性与优越性。然而，上述方法在鲁棒性以及对早期故障的敏感性方面仍然存在着一些不足。在实际工程应用中，轴承振动难免受到环境噪声、退化过程的随机性等因素的干扰，使得时、频域退化特征出现随机波动，从而影响轴承性能退化评估结果的稳定性和准确性。尤其在轴承的早期退化阶段，轴承故障信息往往淹没在复杂的噪声中[13]，导致难以从中提取出有效的退化特征。

退化特征提取是轴承性能退化评估的关键。与常用的时、频域等特征提取方法不同，从数据的统计分布出发来提取设备的退化特征更适用于复杂的工程数据，具有更强的鲁棒性。在常见的统计分布模型中，威布尔分布是一种适用于机械零部件磨损累计失效的概率分布模型，更接近实际工程信号的长尾分布特点，在可靠性工程及寿命预测领域中被广泛应用[14]。Ben Ali等[15]利用威布尔分布对均方根、峭度等轴承时域特征进行拟合，以抑制原始时域特征中的随机波动和异常干扰，然后把威布尔分布处理后的特征输入神经网络实现轴承的剩余寿命预测。陈昌等[16]将威布尔分布与最小二乘支持向量机结合对滚动轴承的退化趋势进行了预测。侯美慧等[17]将威布尔分布用于岸桥铰点的退化特征提取中，用实际工况下的岸桥全寿命铰点振动数据验证了威布尔分布在复杂工程应用中的有效性和可行性。然而，上述方法均采用原始信号的威布尔分布参数作为衡量设备性能状态变化的指标，反映的故障特征信息不够充分，且原始数据包含噪声干扰，有必要对其进行预处理以增强故障特征信息。

滚动轴承的振动信号是一种典型的非平稳信号，对其进行多尺度分解有助于从不同尺度上挖掘轴承隐藏的性能退化特征。经验模态分解(empirical mode decomposition，EMD)是一种自适应的多尺度信号分解方法，能够把复杂信号分解为一系列本征模态分量(intrinsic mode function，IMF)分量之和，非常适用于处理非平稳信号，在轴承及齿轮箱的故障诊断中已得到成功应用[18-19]。

综上，为了提高轴承性能退化结果的准确性与稳定性，针对轴承振动信号非平稳性的特点，并且鉴于威布尔分布参数对性能退化的敏感性以及HMM在性能退化评估中的有效性，本文提出一种基于EMD多尺度数据威布尔分布与HMM的轴承性能退化评估方法。该方法通过EMD对轴承振动信号进行多尺度分解，并筛选出包含明显故障信息的IMF分量，然后从各IMF分量中提取威布尔形状参数作为退化特征，最后结合HMM模型达到准确评估轴承性能退化状态的目的。

1 理论基础

1.1 经验模态分解

EMD是一种适用于非平稳、非线性信号的多尺度信号分解方法。它可以将信号自适应地分解为若干个本征模态函数之和。分解得到的IMF分量须满足以下两个条件：

(1) 该函数局部极值点以及过零点的数量最多只能相差一个。

(2) 任意时刻，该函数的上包络线与下包络线关于x轴局部对称。

对任一信号x(t)使用EMD方法进行分解主要包括以下步骤。

步骤1确定该信号的局部极大值以及极小值点，并分别拟合形成上包络线及下包络线。

步骤2计算上下包络线的平均值m1，将均值拟合，得到均值包络线。

具体过程为：首先，计算中间信号

h1=x(t)-m1

(1)

然后，判断h1是否满足上述IMF分量的两个条件，若满足，则h1为原始信号x(t)的第一个IMF分量；若不满足则以该中间信号作为初始信号，然后重复步骤1～步骤2得到新的IMF分量。

步骤3当函数hn满足单调序列或是常数时，循环结束，记r=hn，r称为残余分量。最终，原始信号可表示为多个IMF分量及残余分量的线性叠加。

1.2 威布尔分布

常用的三参数威布尔分布的概率密度函数f(x)及累计密度函数F(x)可表示为

(2)

(3)

式中:β为形状参数；η为尺度参数；γ为位置参数。在威布尔分布的3个参数中，形状参数决定了分布曲线的基本形状，尺度参数可以将曲线进行放大和缩小，而位置参数决定了函数坐标系中的起始位置。由此可见，位置参数与轴承性能退化趋势无关，因此主要分析数据威布尔分布的形状参数与尺度参数。

本文采用极大似然估计方法对威布尔分布的主要参数进行估计。给定观测序列X=x1,…,xN，则在形状参数β与尺度参数η下观测序列的极大似然函数为

(4)

式中，N为数据长度。根据极大似然原理及牛顿迭代法可求解得到形状参数β及尺度参数η。

1.3 隐马尔可夫模型

隐马尔可夫模型是一种双重随机过程，不仅状态转移过程是随机的，且每一状态的观测值也是随机的。其中，描述状态转移的马尔科夫过程不能被直接观测，但可利用观测过程进行估计。

一个典型的离散HMM可用以下参数描述

(1)N：模型中隐含状态的数目。假定状态集合为S={S1,…,SN}，记t时刻模型所处状态为qt，则qt∈{S1,…,SN}。

(2)M：各个状态对应的观测值数目。记观测值集合为V={v1,…,vM}，若ot为t时刻的观测值，则ot∈V。

(3) π：初始状态的概率向量π={π1,…,πN}，其中πi=P(q1=Si)，1≤i≤N。

(4)A：状态转移概率矩阵A=(aij)N×N，其中aij=P(qt+1=Sj|qt=Si),1≤i,j≤N。

(5)B：观测值概率矩阵B=(bjk)N×N，其中bjk=P(ot=vk|qt=Sj),1≤j≤N,1≤k≤M。

典型的HMM可由以上5个参数表示，常被简记为λ=(π,A,B)。在HMM的训练过程中，常采用期望最大法对模型参数进行估计。然后，在模型参数已知的条件下，采用前向-后向算法对观测序列O={o1,…,oT}的概率P(O|λ)进行计算。在前向算法中，定义前向变量为1～t时刻观测序列与t时刻模型处于状态Si的联合概率，即

αt(i)=P(o1,…,ot,qt=Si|λ)

(5)

记前向变量的初始值α1(i)=πibi(o1)，则前向变量的迭代公式为

(6)

迭代结束后，观测序列的概率可表示为

(7)

为了防止数据下溢，通常对式(7)的概率取对数，采用对数似然概率log(P(O|λ))作为给定模型参数情况下观测值出现的概率，简写为LLP。

需要说明的是：上述HMM的观测值序列及相应的观测概率矩阵B是离散的，然而轴承的振动信号是连续的，需要用连续的概率密度函数来描述观测值的概率分布。由于混合高斯模型可以无限逼近任意分布，因此常采用混合高斯模型来拟合任意状态Sj下观测值ot的概率密度函数，即

(8)

式中:Mj为状态Sj的高斯分量数目;wjm为第m个高斯模型的权重;μjm与Vjm分别为第m个高斯分布的均值向量与协方差矩阵。

2 基于EMD多尺度威布尔分布与HMM的轴承性能退化评估方法

本文提出了基于EMD多尺度威布尔分布与HMM的轴承性能退化评估方法。该方法首先采用经验模态分解将轴承振动数据分解到不同尺度的IMF分量中，并以峭度为标准选取IMF分量；然后提取每个IMF分量的威布尔分布形状参数构建退化特征集，将其输入HMM模型中进行训练，最终得到性能退化评估曲线。该方法分为离线训练与在线评估两个阶段，具体流程如图1所示，主要步骤包括以下。

图1 轴承性能退化评估流程图

步骤1利用EMD对采集的轴承振动数据进行多尺度分解。

步骤2分别计算各个IMF分量的峭度值，筛选得到故障信息明显的多个IMF分量。

步骤3对各个IMF分量分别对一个滑动时间窗口内的数据进行威布尔分布拟合，提取威布尔分布的形状参数构建多维性能退化特征集。

步骤4将轴承正常状态下的数据通过步骤1～步骤3得到的退化特征集作为训练数据输入HMM模型，得到性能退化评估模型λ。

步骤5在线评估阶段，将待测数据的退化特征输入评估模型λ，利用前向算法中的式(7)计算得到待测数据满足正常模型λ的对数似然概率LLP。若LLP的值较大，说明待测数据处于正常状态的概率较高；反之，若LLP的值变小，则说明轴承当前待测状态开始偏离正常状态。因此，可根据LLP的变化趋势及时发现轴承的早期故障，评估其退化状态。

3 试验验证

本文使用的滚动轴承全寿命周期振动信号来自杭州轴承试验中心的ABLT-1A轴承寿命强化试验台。试验装置如图2(a)所示，主要包括台架、加速度传感器、测试轴承、采集系统等，可同时进行4个轴承的加速疲劳试验。4个测试轴承及3个加速度传感器(通道1～3)的安装位置如图2(b)所示。

(a) 试验测试装置

试验轴承采用型号为6307的单列深沟球轴承采样频率为25.6 kHz，每一分钟采集一次，采样时间为0.8 s。本文采用最先失效的B12轴承对所提出的方法进行验证。试验从开始采集直至轴承完全失效共经历2 469 min，得到2 469组振动数据。

轴承全寿命信号的时域波形及有效值曲线如图3所示。从图3可以区分出轴承的正常阶段及失效阶段，但对轴承早期故障的识别较为困难。同时，如图3(b)所示，有效值指标存在波动和异常值，不宜直接作为轴承的性能指标。

(a) 轴承全寿命数据

3.1 EMD分解与威布尔分布校验

首先，对轴承振动信号进行EMD分解，然后根据峭度值的大小筛选出主要的IMF分量，得到高于各分量平均峭度值的6个IMF分量，即IMF2、IMF3、IMF10、IMF11、IMF13、IMF17。对各IMF分量进行威布尔分布校验，检验结果表明各分量数据均近似满足威布尔分布。其中，IMF3分量整体数据的威布尔分布校验结果如图4所示。

图4 IMF3威布尔分布校验

3.2 特征提取

以IMF3为例，取窗口大小为10 min，步长为1 min，每分钟含20 480个数据点，对该分量分别提取威布尔尺度参数、形状参数等特征指标，结果如图5所示。可见，尺度参数不随故障发展呈单调变化趋势，且无法有效反映轴承早期故障，不宜作为评估指标。相比之下，形状参数与轴承退化趋势较为一致，单调性较好，且在轴承的退化初期有所变化，因此本文选取形状参数作为轴承的退化特征。

(a) IMF3尺度参数

对6个IMF分量分别提取威布尔形状参数，组成多维退化特征集W=[W1,W2,W3,W4,W5,W6]。

3.3 HMM模型训练与性能退化评估

在HMM训练阶段，选择轴承正常状态下的前300组数据作为训练集，将其退化特征集作为HMM的输入，建立HMM模型进行训练。试验中，混合高斯模型数M取2，隐含状态数N取5。

在评估阶段，将轴承全寿命数据作为测试集进行性能退化评估，得到的性能退化曲线如图6所示。可见，第1 291 min前HMM输出概率值较高且较为平稳，表明轴承运行状况良好；在1 291 min后，性能指标开始下降，表明轴承早期故障的出现；在2 330 min后，曲线出现剧烈下降，表明轴承性能大幅退化，进入故障失效阶段。

图6 基于EMD多尺度数据威布尔分布与HMM的性能退化评估结果

为了验证所提出方法的准确性，对两个变化点即第1 291 min和2 330 min的数据进行包络谱分析。本试验中，轴承转频fr=50 Hz，内圈故障特征频率fi=246 Hz。1 291 min与2 330 min的时域波形及包络谱如图7所示。如图7(a)所示，1 291 min的时域波形中存在微弱的冲击成分。同时，如图7(c)所示，从其包络谱中可以观测到转频fr及内圈故障频率fi成分，说明此时轴承出现早期故障，发生轻微点蚀。在第2 330 min，图7(b)的时域波形振幅明显增大，冲击信号增多。同时，如图7(d)所示，包络谱中故障特征频率分量fi较为明显，且受转频调制影响，在fi附近出现明显的边频带，表明轴承出现内圈严重故障。上述分析结果验证了评估结果的正确性与有效性。

(a) 1 291 min时域波形

为了进一步验证方法的优越性，将所提方法与以下6种方法进行对比分析。采用常用的11个时域特征指标、时域特征与主成分分析(principal components analysis，PCA)降维、时域特征与局部保持投影(locality preserving projections，LPP)、EMD后各分量的能量指标、原始单一尺度信号的威布尔形状参数、EMD多尺度数据的威布尔尺度参数等退化特征提取方法分别与HMM相结合，得到的性能退化曲线如图8所示。可见，图8(a)～图8(d)的性能指标曲线在轴承正常阶段均存在异常波动，容易导致误判，且无法识别出轴承的早期故障。图8(e)的性能退化曲线在正常阶段虽然比较平稳，但是仍然难以识别轴承的早期故障。图8(f)所采用的方法与本文最为相近，从图8可以看出该方法虽然对轴承的早期故障有一定反映，但是整体趋势不明显，波动较大，说明本文采用的威布尔分布的形状参数比图8(f)采用的尺度参数更适合作为性能退化特征。根据图8划分的具体退化阶段评估结果如表1所示。可见，仅有本方法能及时识别出轴承的早期故障，并能准确区分出轴承的早期故障与严重失效阶段。