基于谱峭度和Teager能量算子的轴承故障特征增强

刘建春，胡俊锋，熊国良，唐晓红，张 龙

(1.江西现代职业技术学院 信息工程分院，南昌 330095；2.华东交通大学 机电与车辆工程学院，南昌 330013；3.南昌铁路局 科学技术研究所，南昌 330002)

滚动轴承是机械设备中应用最广泛、亦是最易出现故障的零部件，其运行状态直接影响到整个机械设备的性能[1]。滚动轴承故障会引起设备的异常振动与噪声，严重干扰设备的正常运行，轻则造成巨大的经济损失，重则威胁人员安全。因此，开展滚动轴承的状态检测与故障诊断研究显得十分必要，受到国内外学者的广泛关注。当滚动轴承出现表面损伤类故障时，其损伤部位将会和与之接触的轴承元件发生碰撞，产生周期性冲击成分。然而在现实工程中，受工作环境、振动传输路径和多振动源相互耦合等因素影响，最终采样得到的数据往往包含众多噪声干扰，表征着故障信息的周期性冲击成分往往被噪声掩盖，这无疑加大了滚动轴承故障诊断的难度。寻求行之有效的故障特征提取方法、及时发现滚动轴承故障一直是故障诊断的热点问题[2–3]。

共振解调法是目前应用最广泛的故障诊断方法之一，其关键在于寻求一个合适的带通滤波器对信号进行噪声消除，然而带通滤波器的参数（即中心频率和带宽）选择往往取决于研究人员的经验知识和海量数据分析，在获得满意结果之前需要进行多次调试，费时费力。谱峭度能够通过计算信号每根谱线峭度值的方式发现发生冲击的频段，以此寻找最优中心频率和带宽，构造最优带通滤波器对信号进行处理，实现信号中的噪声消除[4]。谱峭度拥有十分优异的消噪性能，在故障诊断领域的应用极为广泛。丁康等结合Morlet小波的滤波特性，用于改善谱峭度滤波器的参数选择问题，提高了故障诊断效果[5]。程军圣等针对齿轮故障信号的非平稳调制特性，提出基于局部均值分解时频分析的谱峭度分析方法，成功提取出齿轮故障振动信号特征[6]。苏文胜等以EMD作为预处理器突出信号中的高频共振分量，而后对滤波信号谱峭度处理进一步减少噪声干扰，用于解决早期故障信号中的故障冲击特征难以提取的问题[7]。然而谱峭度在对背景噪声较大的低信噪比信号进行处理时，由于表征着故障信息的循环冲击成分被噪声掩盖严重，处理后的信号噪声虽然得到了极大的消除，但是带内噪声残留仍然较为严重，极大地影响了后期诊断的准确性。以上研究均未对谱峭度的带内噪声问题进行深入探讨。针对该问题，本文提出利用Teager能量算子[8]能够有效追踪系统能量的特性，在信号能量的角度进一步消除谱峭度处理后的带内噪声，使故障特征信息得到进一步增强。

1 方法介绍

如图1所示为本文提出的SK-Teager能量算子故障诊断模型，其诊断步骤如下：首先对振动信号进行采集得到原始振动信号，由于受到噪声干扰较大，故障冲击成分被噪声淹没，因此需要对该信号进行处理以消除噪声干扰，突出故障脉冲信息；谱峭度能够根据信号特征处理得到谱峭度图，从谱峭度图中可得到最优中心频率及带宽，以此构成最优滤波器对信号进行滤波处理得到SK滤波信号，初步消除信号中的噪声干扰；而后通过Teager能量算子跟踪该滤波信号的瞬时总能量E，从系统能量的角度进一步消除SK处理过程中残余的带内噪声；最后进行包络谱分析提取出包络谱中的故障特征频率，通过对故障特征频率理论值的对比分析获得最终诊断结果。

1.1 谱峭度

谱峭度最早由 Dwyer[9]提出，而后 Antoni[10]在此基础上基于FIR带通滤波器提出了快速谱峭度算法，谱峭度的定义如下

其中Sn（f）为信号的n阶谱距，其数值等于E〈|X(f，t)|n〉，E〈·〉为均值计算，|·| 为取模计算，X(f，t)为信号x(t)处于频率f的复包络。

快速谱峭度是一种基于1/3二叉树滤波器结构的滤波方法，算法的详细解释如下：

图1 SK-Teager能量算子故障诊断模型

(1)构造一截止频率fc=1/8+ε(ε≥0)的标准低通FIR滤波器h(n)（此处采样频率为归一化频率，即fs=1）。同时，基于h(n)分别构造准解析低通滤波器h0(n)（分析频带[0，1/4]）和高通滤波器h1(n)（分析频带[1/4，1/2]）

(2)如图2所示，分别以低通滤波器h0(n)和高通滤波器h1(n)对滤波结果进行2倍降采样，依次迭代，便可得到如图3所示滤波器树及其滤波结果。其中为第k层第i个滤波器的处理结果。

(3)根据式(1)对以上每一个滤波结果进行峭度计算

(4)将所有的谱峭度汇总，便可得到谱峭度图。

图2 低通和高通滤波器分解图

图3 二叉树滤波器组及其滤波结果

为提高分析精度，谱峭度对频带进行进一步划分，即将3×2k个滤波器加入k+1和k+2级之间，类似于上述方法，以构建三个准解析带通滤波器gj(n（)频带分别是：[0,1/6]、[1/6,1/3]、[1/3,1/2]）的形式对每个进行分解，最终获得三个子代信号

1.2 Teager能量算子

原始振动信号x(t)经过谱峭度处理后获得滤波信号c(t)，则该信号的能量算子Ψ定义为

假设质量为m的质量块和刚度为k的弹簧组成一个线性无阻尼振动系统，则该振动系统的运动方程

式中c(t)为质量块位移量；A为振动幅值；ω为固有（圆）频率；θ为初始相位角。

振动系统处于任意时刻的瞬时总能量E的取值总等于质量块的动能及弹簧的势能之和，即

根据上式可得出结论：振动系统的瞬时总能量E与振动幅值A及固有（圆）频率ω的平方成正比。结合公式(4)和公式(5)

对比式(6)与式(7)，发现信号的能量算子Ψ和系统瞬时总能量E之间仅相差一个m/2，为常数。因此得出结论：Teager能量算子Ψ能够有效跟踪振动系统总能量E。

公式(4)为连续时间信号的Teager能量算子定义，对于离散时间信号x(n)而言，其能量算子Ψ将变换为

由上式可知：若要得到任意时刻下离散时间信号c(n)的能量算子Ψ，只需知晓其中三个采样点便可。因此Teager能量算子对于信号中的瞬态成分十分敏感，拥有十分优异的时间分辨率。

2 仿真数据分析

为验证所提方法在滚动轴承故障诊断中的有效性，根据滚动轴承真实运动状态构造存在外圈故障的仿真信号，采样频率12 kHz，外圈故障特征频率为80 Hz。其时域波形如图5(a)所示，由于本文针对的是强噪声状态下的低信噪比信号，因此加入高斯白噪声能量较大，原始信号中的循环冲击成分被噪声掩盖较为严重。图4是对该信号谱峭度分析后获得的谱峭度图，从中可以看出其最大峭度值处于第3层（颜色最深），获得最优滤波器参数：中心频率1 875 Hz，带宽750 Hz。

图4 仿真信号谱峭度图

根据该参数构造带通滤波器，并将原始信号图5(a)通过滤波器处理，得到SK滤波信号如图5(b)，经过谱峭度处理后，原始信号中的噪声分量得到了极大程度的消除，信噪比得到了一定的提升，但是仍然存在带内噪声，且对于循环冲击分量的干扰效果较为严重，突出效果并不明显。利用Teager能量算子对图5(b)的瞬时值及其微分非线性组合进行计算，估计整个信号源产生动态信号时所需的总能量，获得图5(c)所示结果。显而易见，图5(c)在图5(b)的基础上进一步消除了噪声干扰，且被抑制到较小程度，从信号能量的角度二次增强了循环冲击分量的凸显效果，故障特征提取效果极为优异。其包络谱中能够轻易提取出故障特征频率80 Hz及其倍频成分，地毯噪声极小，可判别该信号存在严重外圈故障。

图5 仿真外圈故障诊断

3 实验数据分析

3.1 外圈故障

为了进一步验证所提方法在实际故障状态下的轴承信号处理效果，搭建如图6所示试验平台，对型号为N205EM的圆柱滚子轴承以线切割的加工方式模拟轴承的内圈故障、外圈故障和滚动体故障，故障直径分别为0.05 mm、0.17 mm和0.45 mm。考虑到传感器分布造成的采样差异，在测试轴承座的3点钟、12点方向和底座上均安置有加速度传感器。本节采用底座加速度传感器（测量距离最远）测得的直径为0.05 mm(故障最小)的外圈故障轴承数据进行分析。采样过程中主轴转速1 188 r/min，采样频率12 kHz。理论外圈故障频率为108.82 Hz。

图6 旋转机械振动分析及故障诊断试验平台

由图8(a)可知，由于故障尺寸较小，且测点距离故障点位置较远，因此原始振动信号噪声干扰严重，几乎完全掩盖了故障冲击成分，从时域波形图中无法直观提取等间隔故障冲击分量。对该信号谱峭度分析，得到谱峭度图如图7所示。

图7 实验室外圈故障信号谱峭度图

可知谱峭度最大值处于第2层第3个频带（中心频率3 750 Hz，带宽1 500 Hz）。图8(b)为由最优参数构造的带通滤波器对原始信号处理后的时域波形图，从中可以看到谱峭度能够将循环冲击分量和噪声分量有效分离，降低信号噪声的同时，也致使时域波形图中出现了周期性冲击成分，但是仍然存在如1 200和4 200数据点附近的冲击特征突出不明显，难以确定冲击位置等现象。图8(c)相较图8(b)而言则在信号能量的角度二次消除了噪声分量，且将图8(b)中被掩盖的不确定冲击分量进一步凸显，这有利于后期的故障判别。在图8(d)包络谱中出现了明显的等间隔谱峰，且其值和理论故障特征频率及其倍频成分一一对应，表明本文方法已成功提取隐藏于强噪声干扰下的微弱故障循环冲击成分。

3.2 内圈故障

图8 实验室外圈故障诊断

实验数据来源于美国Case Western Reserve University轴承数据中心[11]，试验台的整体结构如图9所示，包含电机、编码器、功率计及内部控制等四个主要部分。测试轴承有两种，分别安装于电机驱动端和风机端。于电机驱动端、风机端和底座上安装了加速度传感器以记录轴承的运行状态。由于考虑强噪声问题，选取传递路径较远的风机端加速度传感器采样获得的驱动端内圈故障轴承数据进行分析，已知轴承型号为SKF6205-2RS，故障尺寸0.356 mm，采样频率12 kHz，主轴转速1 797 r/min，内圈故障特征频率162.19 Hz。

图9 滚动轴承故障模拟实验台

由于选取的故障尺寸较小，且传感器测量位置与故障轴承的传递距离较远。因此图11(a)中的时域波形中夹杂着较大的背景噪声干扰，异常冲击分量难以直观提取出来。根据图10谱峭度图寻觅出的最优频带参数构造出的带通滤波器对原始信号处理后，相比原始信号而言。

图11(b)中的噪声得到了一定程度的压制，隐藏于原始信号的大部分冲击点位置已经可以基本辨别，但仍存在部分幅值能量较小的冲击成分被残存的带内噪声压制现象。运用Teager能量算子追踪图11(b)的系统总能量后，其循环冲击特征成分在能量的角度进一步被加强，噪声消除效果极为明显。其包络谱中的基频162.1 Hz与故障特征频率理论值接近，且倍频成分及转频调制十分明显，包络谱中呈现出典型的内圈故障特征。因此可判定该轴承存在内圈故障，与实际情况相符合。可见本文方法在早期微弱故障特征提取中的有效性。

图10 实验室内圈故障信号谱峭度图

图11 实验室内圈故障诊断

3.3 混合故障

以图6所示试验平台为载体，通过在N205EM圆柱滚子轴承的内圈和外圈分别加工出直径为0.05 mm的故障以模拟滚动轴承的混合故障，采样频率12 kHz，主轴转速1 175 r/min（即转频19.58 Hz），计算出轴承外圈故障特征频率约104.09 Hz，内圈故障特征频率约152 Hz。

通过安装于测试轴承座12点钟方向的加速度传感器采样得到如图13(a)所示数据，由于加工的故障尺寸较小，原始信号中的噪声成分能量较大，大部分故障脉冲信息被掩盖。通过谱峭度对图13(a)处理后，从图12谱峭度图中可得到最优带通滤波器参数：中心频率3 500 Hz，带宽1 000 Hz。

图12 实验室混合故障信号谱峭度图

原始信号经最优带通滤波器处理后获得如图13(b)结果，从中可发现原始信号中的噪声被大幅消除，信噪比得到了较大程度的提升。图13(c)则以系统能量的角度在图13(b)的基础上进一步突出了故障特征信息。图13(d)包络谱中的峰值频率及意义如表1所示，包络谱中103.6 Hz、208.6 Hz、311.2 Hz分别对应外圈故障特征频率的1～3倍频；151.6 Hz、304.7 Hz、456.2 Hz对应内圈故障特征频率的1～3倍频；且内圈故障频率的周边存在间隔为转频19.58 Hz的边频成分。以此可判定该数据存在内圈和外圈混合故障。最终结论与实际情况相符，可见本文方法在处理混合故障时亦是完全适用的。

图13 实验室混合故障诊断

表1 图13(d)各峰值频率成分数值及意义

4 结 语

(1)谱峭度能够以最大峭度值为原则，根据信号特征选取最优参数，构成最优带通滤波器对信号进行有效处理。解决了共振解调法带通滤波器参数难以选取的问题。

(2)谱峭度对数据进行处理时，会出现带内噪声残余问题，严重影响后续诊断的准确性。在此利用teager能量算子能够有效跟踪信号能量的特性，通过信号能量的角度进一步消除带内噪声，增强故障冲击特征，有利于后期的轴承故障识别。

(3)经过实验室数据分析表明：本文方法不仅适用于滚动轴承单一故障的诊断，对于轴承的混合故障也同样适用。

[1]王建东，马增强，李延忠，等.一种自适应调整滤波参数的轴承故障诊断方法[J].噪声与振动控制，2017，37(1)：137-141.

[2]夏均忠，苏涛，马宗坡，等.基于EMD的滚动轴承故障特征提取方法[J].噪声与振动控制，2013，33(2)：123-127.

[3]张文义，于德介，陈向民.基于信号共振稀疏分解与能量算子解调的轴承故障诊断方法[J].中国电机工程学报，2013(20)：111-118.

[4]马朝永，张学飞，胥永刚.形态分量分析和谱峭度在轴承故障诊断中的应用[J].噪声与振动控制，2016，36(4)：169-173.

[5]丁康，黄志东，林慧斌.一种谱峭度和Morlet小波的滚动轴承微弱故障诊断方法[J].振动工程学报，2014，27(1)：128-135.

[6]程军圣，杨怡，杨宇.基于LMD的谱峭度方法在齿轮故障诊断中的应用[J].振动与冲击，2012，31(18)：20-23.

[7]苏文胜，王奉涛，张志新，等.EMD降噪和谱峭度法在滚动轴承早期故障诊断中的应用[J].振动与冲击，2010，29(3)：18-21.

[8]王冰，李洪儒，许葆华.基于LMD和平滑Teager能量算子解调的电机滚动轴承故障特征提取[J].机械传动，2012，36(09)：89-92.

[9]DWYER R F.Detection of non-gaussian signals by frequency domain kurtosis estimation[J].Acoustic，SpeechandSignalProcessing.Boston:IEEE International Conference on ICASSP,1983,83:607-610.

[10]ANTONI J.Fast computation of the kurtogram for the detection of transient faults[J].Mechanical Systems and Signal Processing,2007,21(1):108-124.

[11]Case Western Reserve University.Bearing data center EB/OL.2014-1-6.http://csegroups.case.edu/bearing data center/ pages/ welcome- case- western- reserveuniversity-bearing-data-center-website.