基于优化VMD与改进阈值降噪的滚动轴承早期故障特征提取

2021-07-14赵小强

振动与冲击 2021年13期

陈鹏, 赵小强,2,3

(1. 兰州理工大学电气工程与信息工程学院，兰州 730050; 2. 甘肃省工业过程先进控制重点实验室, 兰州 730050) 3. 兰州理工大学国家级电气与控制工程实验教学中心, 兰州 730050)

滚动轴承作为旋转类机械设备的关键零部件，其运行状态直接关系到整个设备的安全状况[1]。在其旋转过程中，由于与其它元件接触摩擦引起损伤，当反复经过损伤点时，采集的振动信号会产生周期性脉冲[2]，该脉冲是判断滚动轴承发生故障的主要依据。但是轴承的运行环境复杂，常受到传输路径和强噪声的干扰，导致早期故障的微弱信号淹没在噪声中而难以提取[3]。因此，如何有效地分析滚动轴承振动信号，提取出故障引起的冲击特征进行故障辨识是整个故障诊断的关键。

对于滚动轴承的振动信号分析而言，应用信号处理方法对振动信号进行降噪的意义是显而易见的。近年来，在小波降噪的基础上，学者们提出了许多自适应信号处理降噪方法，并应用于滚动轴承的故障特征提取[4]。如经验模态分解(empirical mode decomposition，EMD)[5]、集总经验模态分解(ensemble empirical mode decomposition， EEMD)[6]、互补集合经验模态分解(complete ensemble empirical mode decomposition,CEEMD)[7]、局部均值分解(local mean decomposition,LMD)[8]、局部特征尺度分解(local characteristic-scale decomposition， LCD)[9]和固有时间尺度分解(intrinsic time-scale decomposition，ITD)[10]等。这些方法都在不同程度上实现了滚动轴承振动信号的降噪，但是这些方法都是基于递归分解原理，没有从本质上解决其模态混叠和端点效应的缺点。而Dragomiretskiy等[11]提出了一种新的信号自适应分解方法——变分模态分解(variational modal decomposition, VMD)。由于其通过引入变分约束来实现信号的分解，克服了上述方法存在的模态混叠和端点效应问题，被广泛应用于旋转机械的振动信号降噪[12]。然而VMD存在模态分解个数k和二次罚因子α需要进行预设置且不同参数值对于分解效果影响大的问题。对于此问题，如通过手动调节较难获得最优参数。因此，如何自适应地实现VMD最优参数选择，对于信号的分解效果和后续特征提取至关重要。对此，王奉涛等[13]提出了基于能量的模态个数k选择；Lian等[14]提出了预设置k的范围，在设置的范围内分解信号，计算不同分解模态个数下的能量损失来获得最佳k。两位学者采取预设一定的范围值，对k进行了优化，自适应性不够。唐贵基等[15]利用粒子群优化算法优化组合参数，建立了基于最小包络熵的适应度函数；刘建昌等[16]详细分析了组合参数对VMD分解效果的影响，并通过遗传变异粒子群算法对其进行优化，适应度函数采用了唐贵基等人提出的最小包络熵；Yan等[17]提出了基于最小包络熵的CSA优化VMD分解方法。学者们通过智能优化算法实现VMD组合参数的自适应选择，建立了能够表征信号稀疏性的包络熵适应度函数，实现了信号的分解降噪。可见，通过选择全局和局部搜索能力较强的优化算法，建立合理的适应度函数，可实现VMD参数的自适应选择分解降噪而获得较好的效果。

Mirjalili等[18]在2016年提出了WOA算法，证明具有较强的局部和全局搜索能力。峭度由于可以检测脉冲的瞬变性而应用于选择模态分量，但是对异常噪声敏感、稳定性不高。L-峭度[19]是一种描述信号分布特性的指标，是不同阶次统计的线性组合，相比峭度只是样本的四阶统计量而言，稳定性更高，能够消除异常值和噪声信号的影响，可以更好的检测周期性脉冲[20]。相关系数可以表达模态分量和原始信号的相关性，但对噪声敏感。

基于以上分析，本文构建L-峭度与相关系数的自适应度函数，利用WOA优化VMD参数实现振动信号的自适分解降噪，通过构建的准则选取最优模态分量，再通过改进阈值降噪技术对选取的最优模态分量进行进一步降噪，实现振动信号故障特征的有效提取和故障辨识。

1 VMD算法

VMD的本质是通过迭代搜寻变分模型，把信号x(t)分解为具有一定稀疏性的k个模态分量，每个模态分量以中心频率ωk进行分解，最终得到约束变分问题为

(1)

式中:{uk}={u1,…,uk}为分解后的k个模态分量；{ωk}={ω1,…,ωk}为分解后每个模态分量的频率中心。为了求解式(1)，利用二次罚函数项和Lagrange乘子将式(1)转化为如下形式的无约束问题

L({uk},{ωk},λ)=

(2)

式中:α为惩罚因子；λ(t)为Lagrange乘子。

(3)

VMD的具体实现过程如下：

步骤2重复n←n+1；

(4)

更新ωk

(5)

对于所有ω，双上升步长

(6)

直至满足迭代终止条件

(7)

式中,ε为判别精度，ε>0。

2 优化VMD参数与改进阈值降噪

2.1 构建L-峭度和相关系数准则

传统的均值、偏差和峭度等统计量被用于选择有效模态，这些统计量识别异常的方法只是剔除标准差和均值之间偏离甚远的数据，这对于小样本数据而言是不准确的，而恰恰旋转机械设备的故障数据都属于小样本数据，因此常用的以上统计量影响样本数据统计的正确性。而L-峭度是一种描述信号分布特性的指标，是不同阶次统计的线性组合，相比峭度只是样本的四阶统计量而言，稳定性更高，能够消除异常值和噪声信号的影响，可以用来很好地确定最佳模态分量。L-峭度的表达式为式(8)

(8)

如果信号的分布形状扁平，则λ2的值将趋于很大的值，这意味着分布发散，其与标准偏差表征一致，而L-峭度将成为一个很小的值，反之亦然。此外，L-峭度值也由λ4确定，它可以衡量两个极值(即EX4:4-EX1:4)和两个中心值(即EX3：4-EX2：4)之间的距离。如果分布是平坦的，则信号分布通常是均匀的，并且λ4接近于零。另一方面，如果分布具有尖峰态，则两个极值之间的距离通常较大，而两个中心值之间的距离较小，且λ4较大。因此，L-峭度，也可以代表分布特征。

同时，相关系数的表达式为(9)，其可以表征分解后模态分量与原始信号的相关程度，但对噪声敏感。

(9)

因此，考虑两者的优点构建L-峭度和相关系数的适应度函数如式(10)，该适应度函数具备筛选出与原始信号相关性最大的分量和能够稳定检测信号中周期性故障脉冲并有效消除噪声干扰的能力。

fitness=LKCI=(LKI·CI)max

(10)

2.2 WOA优化VMD

设初始最优位置为当前猎物的位置，其它猎物在最初以该位置目标进行移动，并不断利用式(11)和(12)更新位置

D=|C·X*(t)-X(t)|

(11)

X(t+1)=X*(t)-A·D

(12)

A=2a·rand-a

(13)

C=2·rand

(14)

式中：rand表示产生的[0,1]之间随机数，随着循环迭代的增加，a作为收缩因子将从2线性减小到0，其表达式为

a=(2-2t/Tmax)

(15)

式中，Tmax表示设置的最大运行迭代次数。鲸鱼通过式(13)和(14)来实现猎物搜索过程中包围区域的不断缩小。同时，每个鲸鱼又通过螺旋式来更新其与目标位置的距离，其模拟螺旋式更新的数学模型可以表示为如下

X(t+1)=D′·ebl·cos(2πl)+X*(t)

(16)

式中：D′=|X*(t)-X(t)|表示第i只鲸鱼和当前最优位置之间的距离；b为常量系数，通常为b=1;l为[1,-1]之间的随机数。

为了收缩包围和螺旋更新同步进行，通过概率实现对两种方式的选择更新，其数学表达式如下

(17)

当|A|≥1时，将随机选择鲸鱼远离当前最优鲸鱼来寻找最优的鲸鱼，可以提高算法的全局搜索能力，该数学表达式如下

D=|C·Xrand-X(t)|

(18)

X(t+1)=Xrand-A·D

(19)

式中：Xrand表示随机选取鲸鱼的位置向量。

通过WOA对VMD的参数k和α进行组合优化，如图1所示。其步骤如下

图1 WOA优化VMD流程

(1) WOA初始化参数k和α的取值范围为k=[2,15]，α=[200,8 000]以避免不完整的信号分解。

(2) 对振动信号进行VMD分解。

(3) 计算每个组合参数对应位置的适应度值。更新适应度值大于当前值时的最佳适应值。

(4) 确定是否终止迭代。如果t小于Tmax，让t=t+ 1，并且鲸鱼的位置会更新。否则，迭代终止，并保存最佳参数。

2.3 改进阈值算法

硬阈值和软阈值是常用的降噪技术，在各种应用中得到了广泛应用[21]，其数学表达式分别为式(20)和(21)

(20)

(21)

(22)

式中：ki是通过试验设定;σi是第i分量的噪声标准偏差；N是信号的长度。噪声标准偏差通过式(23)获得

(23)

硬阈值函数在阈值处是不连续的，会产生严重的波动。在软阈值函数中，不连续问题得到了解决，但该函数在阈值位置是不可微的，这将导致对于随机噪声信号不能准确地估计。在本文中，引用改进的阈值降噪[22]公式为

(24)

其中u=α(|ci(t)|-λi),0<αandμ≤1。

在式(24)中，当α接近0时，阈值式(24)中的函数具有硬阈值函数功能，当α接近1时，则该阈值函数的功能接近于软阈值功能。可见，该函数具有一定的灵活性，参数α一般设置为0.05。

3 设计的故障特征提取方法步骤

滚动轴承运行工况复杂、环境噪声污染大导致早期故障信号中的故障特征难以提取。本文构建基于L-峭度和相关系数准则，利用WOA优化VMD参数实现振动信号的自适应分解降噪，再选取模态分量中噪声信息少且与原始信号最大相关的信号进行改进阈值降噪，突出故障特征频率。其故障特征提取步骤如图2所示。

图2 故障提取步骤

(1) 对采集得到的振动信号进行WOA优化VMD自适应分解；

(2) 计算分解之后各分量的LKCI指标；

(3) 选择最大LKCI指标分量并对最大指标分量进行改进阈值降噪；

(4) 对降噪后的信号进行希尔伯特分析，获得故障特征频率。

4 试验验证

4.1 试验1

为了验证所提故障特征提取方法的准确性和有效性，构建一组周期性脉冲序列信号并添加高斯白噪声，其仿真信号模型如下

y(t)=s(t)+n(t)

(25)

式中：s(t)为模拟故障冲击的周期性冲击分量;位移常数A0=0.6;阻尼系数g=0.04;重复周期T=0.01 s;固有频率fn=4 000 Hz;采样点数N=4 096;采样频率为20 kHZ;n(t)为高斯白噪声。产生的周期性冲击信号如图3所示，对该冲击性信号进行加噪获得强噪声干扰的混合仿真信号如图4所示。

图3 周期性冲击信号

图4 仿真信号波形

对仿真信号进行直接包络如图5所示，可见在频率500～1 000 Hz，故障特征频率出现了模糊，其余特征频率被噪声淹没而没有得到突出。再对仿真信号按照本文建立的适应度函数，通过WOA优化VMD进行分解，然后进行Teager能量算子解调得到结果如图6所示。其中Teager能量算子解调相比直接对信号进行包络的故障特征更加明显，但在频率800～1 000 Hz之间出现了故障特征淹没。继续通过建立最小包络熵适应度函数，采用WOA优化VMD分解得到包络熵最小的最佳模态分量，然后进行改进阈值降噪，得到包络分析结果如图7所示，显而易见其故障特征频率和其它噪声频率相混合而不易观察，诊断结果较差。最后，应用本文所提出的方法对仿真故障振动信号进行分析如图8所示，从中可以明显的看见所有故障特征频率都得以突出，取得了较好的效果。由上述分析可见，本文所提出方法的可行性和优越性。

图5 仿真信号包络

图6 Teager能量算子

图7 包络熵优化

图8 本文方法故障特征提取

4.2 试验2

工程试验数据选择来自美国西储大学轴承试验中心在6205-2RSJEMSKF型深沟球滚动轴承上采集的振动加速度信号，该轴承上的损伤故障采用电火花加工技术设置。本文以故障直径为0.018 mm、负载为0、转速为1 730 r/min来模拟滚动轴承的早期微弱故障。在采样频率为12 kHZ下进行轴承的加速度信号采集。由滚动轴承的故障特征计算理论，可以计算得到内圈故障特征频率fi=156.1 Hz、外圈故障特征频率fo=103.4 Hz、转频为Fr=29.95 Hz。首先，对内圈进行故障特征提取，其内圈的加速度振动信号时域图如图9所示。

图9 内圈振动信号时域波形

对内圈振动信号应用本文建立的适应度函数通过WOA优化VMD进行分解，其WOA优化VMD的迭代曲线如图10所示。

图10 WOA优化VMD

根据最优组合参数进行VMD分解，得到分解后的IMF分量如图11所示。

(a) IMF1

其中两个IMF分量的LKCI值如图12所示。

图12 IMF分量的LKCI值

选择最优的IMF2分量进行Teager能量算子解调如图13 所示，可见转频、故障特征频率的1倍频和2倍频较为突出，3倍频已较弱。为了进一步对比分析本文方法的有效性，以最小包络熵为自适应度函数，建立WOA优化VMD选取最佳模态分量，通过改进阈值降噪技术进行降噪得到最后的结果如图14所示，从图14中可以看到只有转频特征较突出，对于故障特征频率和2倍频实现了较微弱的识别，其余倍频都没有得到有效的识别。最后，通过本文提出的方法进行分析，其中设定改进阈值降噪中的ki=0.3，得到结果如图15所示。由图15可见，内圈故障特征频率以及2倍、3倍和4倍频得到了显著的突出，清楚的判断出属于内圈故障。综上所分析，本文提出的方法得到了较优的特征提取效果。

图13 Teager能量算子

图14 包络熵优化

图15 本文方法故障特征提取

对外圈进行故障特征提取分析，外圈故障振动信号的时域波形如图16所示。同理，应用本文建立的适应度函数，采用WOA优化VMD进行分解，然后进行Teager能量算子解调如图17所示，可以观察到1倍频到4倍频得到突出，但在高频部分的故障特诊频率被淹没而难以辨识。其次，建立最小包络熵适应度函数，采用WOA优化VMD后分解后，选取最佳模态分量进行改进阈值降噪后的结果如图18所示，可以观察到较微弱的1倍频和2倍频故障特征频率，特征提取效果较差。