基于VMD-SVD能量标准谱-Teager能量算子的轴承故障诊断方法

周 洋，向 阳，黄陈哲，黄进安

（1.武汉理工大学 能源与动力工程学院，武汉 430063；2.武汉理工大学 船舶动力系统运用技术交通行业重点实验室，武汉 430063）

滚动轴承作为减小摩擦和磨损的元件被广泛应用于各种旋转机械设备。当滚动轴承表面出现损伤时，会和配合元件之间发生碰撞，从而激起旋转系统的各阶固有振动，使得产生的振动信号频谱中包含多个共振峰值，使得机械设备产生异常的振动与噪声[1]。然而轻微故障引起的冲击常被强背景噪声淹没，导致难以提取故障特征频率，因此对轴承进行故障诊断，尤其是早期微弱故障诊断一直是学者研究的热点。轴承故障产生的振动信号具有典型的周期性、非平稳特性，包含了丰富的工作状态信息[2]，同时故障信号中常掺杂大量背景噪声，常用的频谱分析方法处理此类信号效果差，难以提取故障特征频率。因此，Huang等[3]提出了经验模态分解（EMD），此方法具有自适应分解特性，处理非线性和非平稳信号具有较高效率。董文智等[4]为解决提取轴承故障特征频率问题，提出了基于集总经验模态分解（EEMD）和奇异值差分谱理论的故障诊断方法。但经验模态分解缺少严格的数学推导，同时存在模态混叠、端点效应和受采样频率影响较大等缺点[5]。基于此，Dragomiretskiy等[6]提出了变分模态分解(VMD)，该方法将信号分解转化为非递归、变分形式的分解模态，具有可靠的理论基础，表现出更好的噪声鲁棒性。Mohanty等[7]利用VMD提取轴承故障特征频率，验证了该方法具有比EMD更好诊断效果。Wan等[8]利用改进的快速谱峭度和VMD相结合的方法，对轴承出现的复合故障进行诊断。刘尚坤等[9]以互信息准则确定VMD分量个数，并与Teager能量算子结合对转子发生局部碰擦和油膜失稳进行时频故障诊断，但在设定互信息的阈值时缺乏理论依据。

针对滚动轴承出现微弱故障时，振动信号表现为调幅调频特性，本文提出了基于VMD-SVD能量标准谱-Teager能量算子联合诊断方法，证明了该方法能有效地提取滚动轴承微弱故障的特征频率。

1 基础理论

1.1 变分模态分解（VMD）

VMD是由Dragomiretskiy等提出信号自适应分解估计方法，其实质是用多个维纳滤波组对信号进行滤波。VMD算法利用了EMD中本征模态函数(IMF)概念，但IMF被重新定义为AM-FM信号

式中：相位φk(t)非递减，φk’(t)≥0；瞬时幅值Ak(t)≥0；并且瞬时幅值Ak(t)和瞬时频率wk(t)=φ’k(t)相对于φk(t)变化缓慢，即在[t-δ,t+δ] (δ=2π∕φk′(t))范围内，uk(t)能视作为幅值为Ak(t)、频率为wk(t)的谐波信号。

VMD是一个变分问题的求解过程，假设每个IMF为具有中心频率的有限带宽，其中心频率和带宽在迭代求解过程中不断更新。变分问题可表示为求K个本征模态函数uk(t)，使得所有uk(t)的估计带宽之和最小，并且各模态函数之和等于分解信号的约束条件。具体分解过程如下：

（1）将各模态函数uk(t)进行希尔伯特(Hilbert)变换，求得解析函数

（2）各模态解析函数分别与预估中心频率e-jωkt相乘，将各模态频谱移至基带

（3）利用H1高斯平滑估计解调信号的带宽，即梯度L2范数平方根，构造受约束变分问题

式中：uk={u1,…,uK}和ωk={ω1,…,ωK}分别为分解得到的各模态分量和中心频率为所有IMF之和；K为迭代次数；δ(t)为冲激函数。

（4）为求解约束变分问题，须引入二次惩罚因子α和Lagrangian乘法算子λ(t（)当信号中混有高斯白噪声时，二次惩罚因子α可确保重构信号的精度；而Lagrangian乘法算子确保约束条件的严格性）将其转变为无约束变分问题，增广Lagrangian表达式为

利用交替方向乘子算法(ADMM)对增广Lagrangian表达式中{uk}、{wk}和λ交替迭代，寻找“鞍点”，获得最优解，即为各模态分量{uk}和中心频率{wk}。

1.2 奇异值能量标准谱

对于一维离散时间序列X={x(1),…,x(N)}，构造m×n阶Hankel矩阵如下

矩 阵 中 ：1＜n＜N，m=N-n+1，延 时 值 为 1，则H∈Rm×m。

存在矩阵U∈Rm×m，V∈Rm×m使得

奇异值标准能量谱定义为[10]

信号中的各成分在奇异值能量标准谱中表现为：有用信号谱线幅值大而陡峭，噪声谱线幅值小而平缓。因此，奇异值能量标准谱线出现的拐点，即为有用信号和噪声的分界点。保留分界点及前k个奇异值，将其后奇异值置为零，得到降噪后矩阵H′

式中：{u1…,uk}为U的前k个列向量；{v1…,vk}为V的前k个列向量。

将矩阵H′中相应位置的元素相加求平均值[11]，即为降噪后信号。

1.3 Teager能量算子

Teager能量算子(Teager Energy Operator，TEO)是Teager在研究非线性语音建模时提出的信号分析算法[12]，记作ψ。设有连续时间信号x(t)，则Teager能量算子定义为

式中：x˙(t)和x¨(t)分别为连续信号的一阶和二阶微分，ψ[x(t)] 能表达振动信号的瞬时能量值。

对于离散信号x(n)，利用差分代替微分形式，则式(10)可以近似表示为

Teager能量算子对离散信号相邻的采样点进行计算时对调制信号具有良好的时频分辨率，能检测出信号中的瞬态冲击成分。

2 算法步骤及流程

针对滚动轴承早期微弱故障信号常受噪声干扰、不易提取特征频率的问题。本文提出了基于VMD-SVD能量标准谱-Teager能量算子的轴承故障诊断方法，其具体步骤和流程如图1所示。

(1)对原始信号采取进行去除趋势项处理，消除基线漂移的干扰；

(2)对预处理后的信号进行VMD分解，设置初始模态数K=2，惩罚因子为2000以及带宽为0。模态数K依次增加1，观察不同模态数K下，各IMF分量的频谱和中心频率，确定过分解时的模态数K，将其减1即为最优模态数K；

(3)计算最优模态数K下各IMF分量的峭度和相关系数，选取包含故障状态信息的敏感IMF分量；

(4)构造敏感IMF分量的Hankel矩阵，并进行SVD分解。利用奇异值能量标准谱，确定拐点，选择有效奇异值个数，对信号进行降噪重构，提高信噪比；

(5)求解重构信号的Teager能量算子并进行FFT变换即得到包络谱；

(6)提取微弱故障特征频率，确定故障原因。

图1 具体步骤和流程

3 实验验证及分析

3.1 实验装置

为了验证本文提出方法的有效性和准确性，将该联合方法应用于实际的滚动轴承故障诊断。本文所使用的数据来源于美国Case Western Reserve大学的轴承中心数据库[13]。该试验平台主要由1台1.5 kW驱动电机、1个扭矩传感器、1台测功仪和电机控制单元组成，试验平台如图2所示。

图2 试验测试平台

驱动端轴承是深沟球轴承，型号为6205-2-RS JEM SKF，具体规格参数如表1所示。

采用电火花技术对电机轴承的内圈、外圈和滚动体进行点蚀来模拟常见故障，采集电机驱动端、风扇端和基座在不同工况下的振动加速度信号。

3.2 滚动轴承内圈单点故障分析

试验中设计4种故障尺寸，分别为0.007、0.014、0.021和0.028 inch(其中1 inch=25.4 mm)；设定4种载荷，分别为0、1、2和3 hp(1 hp=0.75 kW)。选取驱动端振动信号，其中驱动轴转速为1772 r/min，负载为1 hp，故障为微弱尺寸0.007 inch，采样频率fs=12kHz，信号长度为6000个采样点。滚动轴承内圈故障理论公式为

表1 滚动轴承结构参数

式中：fi为内圈故障频率，fn为转频，n为滚动体数量，d为滚动体直径，D为节径，α滚动体接触角。

根据表1中参数和式(12)可求得：fn=29.5 Hz，fi=159.7 Hz

图3(a)、图3(b)和图3(c)分别为内圈单点微弱故障振动加速度信号、频谱和局部频谱图，从图3(a)中可以看出，当轴承出现故障时有明显的冲击成分和强背景噪声干扰；而图3(c)中存在较多峰值点，无法有效确定故障类型。

利用本文提出的方法分析故障信号，首先对原始数据采取去趋势项预处理，消除基线漂移。对预处理后的信号进行VMD分解时应先确定模态数K，通过观察不同模态数K下各IMF的中心频率fc及频谱是否出现重叠来确定最优的模态数K。表2给出了不同模态数K时各IMF的中心频率fc。

从表2可以看出：从K=5开始，求得的IMF之间出现了中心频率相近情况。例如，当K=5时，IMF4和IMF5对应的中心频率仅相差262 Hz，并通过计算这两个IMF分量的频谱，可以判断取K=5时出现了过分解。因此，最优的模态数K=4。对内圈单点故障信号取4个模态分量进行VMD分解，得到4个IMF及其相对应的频谱，如图4(a)和图4(b)所示。

表2 不同K值时各IMF的中心频率

图3 内圈微弱故障振动信号与频谱

图4 各IMF分量的时域波形与频谱

观察图4(a)，可以发现各IMF分量时域波形均存在明显冲击成分，无法直接选取敏感IMF分量；观察图4(b)，可以发现各IMF围绕其中心频率分布，未出现模态重叠现象。

计算各IMF峭度和与预处理后信号之间的相关系数，结果如表3所示。

表3 各IMF峭度与相关系数（K=4）

从表3可以看出IMF3和IMF4的峭度值和相关系数均较大，包含的故障信息也更丰富。因此选取IMF3和IMF4对应位置的元素进行相加，可得到重构信号u(t)，其时域波形和频谱如图5(a)和图5(b)所示。

图5 重构振动加速度信号与频谱

构造信号u(t)的Hankel矩阵并进行SVD分解，求得的奇异值和奇异值能量标准谱分别如图6(a)和图6(b)所示。

从图6(b)中可以看到，振动加速度信号冲击成分主要集中在前面的奇异值能量谱线，当r＞18时，Er＜0.1。为了在保证故障信息不丢失同时最大程度降噪，取前18个奇异值对信号进行降噪重构，得到信号v(t)，图7(a)和图(b)为其时域和频谱。

将图7(a)和图7(b)与图3(a)和图3(b)相比，可以看出降噪后时域波形毛刺明显减少，频谱上的峰值易于区分，便于后续处理。

对降噪信号v(t)进行Teager能量算子解调，图8(a)和图8(b)分别为Teager能量算子时域和频谱。

从图8(b)中可以看出在30 Hz、58 Hz和88 Hz处存在明显峰值，这些数值的最大公约数近似为30，因此可以判断它们分别对应驱动轴的1倍、2倍和3倍转频。由于滚动轴承在运转过程中滚动体和滚道之间会发生相对滑动以及频率分辨率的影响，使得提取的转频和理论转频（29.5 Hz）有一定的差别。同时可以较容易看到160 Hz、318 Hz、478 Hz、638 Hz、798 Hz和956 Hz存在峰值，且这些数值的最大公约数近似为160 Hz（内圈故障特征频率为159.7 Hz），因此可以判断它们分别对应1到6阶内圈故障特征谐波频率。在图8(b)中的218 Hz、260 Hz等处出现峰值，其对应各阶故障特征频率的调制边频带。

图6 奇异值大小与能量标准谱

图7 降噪信号波形与频谱

图8 Teager能量算子波形与包络谱

图9 根据EEMD分解IMF分量的时域波形与包络谱

3.3 与基于EEMD方法对比

采用EEMD算法对与上节相同的故障振动信号进行分解，得到11个IMF分量和1个趋势项。根据EEMD算法的特性[5]，取前3个包含故障信息较多的IMF分量并进行Hilbert变换解调，各IMF分量时域波形和包络谱如图9(a)和图9(b)所示。

从图9(b)能提取1倍转频30 Hz、2倍转频58 Hz和故障特征频率160 Hz，但各IMF分量的包络谱差别不太明显，即验证了EEMD分解存在严重的模态混叠现象，同时在0 Hz附近出现了端点效应，不利于故障特征频率的提取。

根据本文提出的方法，在诊断过程中利用VMD-SVD能量标准谱对预处理后信号进行最大限度的降噪重构，同时采取Teager能量算子增强故障冲击特征，最后对得到的离散时间周期信号进行傅里叶变换，能提取故障特征频率及其谐波频率。与EEMD方法相比，本文提出的方法可以提取故障特征频率及其高阶次谐波频率，利用最大公约数能准确提取故障特征频率，有效地避免了模态混叠和端点效应所导致的不利于提取故障特征频率的问题。

4 结语

(1)将文中所提方法应用于滚动轴承内圈微弱故障信号诊断，根据VMD分解得到各IMF分量中心频率及频谱，能有效确定最优模态数K；利用峭度和相关系数指标能准确选取敏感IMF分量。

(2)利用SVD能量标准谱能够最大限度抑制噪声，对降噪后的信号采取Teager能量算子解调能增

?强故障冲击成分，便于提取转频及高阶故障特征频率。

(3)与EEMD分解方法相比，本文所提方法避免了模态混叠和端点效应所导致的不利于提取故障特征频率的问题。