自适应谱峭度图及其在滚动轴承故障诊断中的应用

何倩，易彩，吴文逸，林建辉

自适应谱峭度图及其在滚动轴承故障诊断中的应用

何倩，易彩，吴文逸，林建辉

（西南交通大学 牵引动力国家重点实验室，四川 成都 610031）

提出了一种新的自适应谱峭度图方法。首先采用CEEMDAN方法分解信号的傅里叶频谱，然后对得到的低频IMF分量与残差项进行分层累加计算，对每一层的频谱趋势寻找极小值点，使其作为频谱划分边界，将划分后的各频段通过正交滤波器并利用峭度指标对划分后信号中的故障信息进行评估，构造出一种新的自适应谱峭度图。该方法能够克服快速谱峭度图因遵循固定划分规则易使检测出的频带不能包含完整的故障信息或是提取到无效分量的情况，并且分别利用轴承外圈故障的仿真信号和试验信号，验证了该方法的有效性。

故障诊断；自适应谱峭度图；频谱趋势；快速谱峭度图

滚动轴承广泛应用于旋转机械中，在现代机械设备中起着非常重要的作用[1]。在正常情况下，轴承在高速和重载下运行，即使轻微损坏对设备也非常危险，如果不能立即检测到损坏，轴承将不稳定运行，严重时甚至损坏设备。因此，有必要实时监测轴承的运行状态并进行故障诊断。

Dwyer[2]提出了一种非常强大的统计工具——谱峭度（Spectral Kurtosis，SK），可以代替功率谱密度来表示信号的高斯成分及其在频域中的位置。虽然谱峭度的应用领域较广，但其正式的数学定义尚未完善。Antoni等[3]对SK进行了深入而系统的研究，详细阐述了相关理论，并给出数学定义，即能量归一化四阶谱累积量。此外，Antoni[4]证明了上述定义对噪声中非平稳信号的特征检测有很好的效果，提出快速谱峭度图（Fast Kurtogram），并广泛应用于故障诊断领域。

滤波器选择的核心是保证所设计滤波器的边界能够完全覆盖故障共振频带。由于故障信号的谐振频率和带宽可能是任意的，因此能够灵活改变滤波器带宽和中心频率的滤波器是滤波器选择的主要方向。由于快速谱峭度方法中由二叉树或二叉树结构的滤波器组构成的滤波器在每一级的带宽和中心频率是固定的，因此在工程实践中可能无法获得最佳共振带。Barszc等[5]提出了根据固定窄带内信号的平方包络谱（Squared Envelope Spectrum，SES）的峭度选择最佳解调频带，虽然这个过程可修改中心频率，但由于不能调整带宽的缺陷，结果仍然容易出现局部最优解。经验小波变换是从采集到的信号中提取有意义模态的有力方法之一[6]，可将其视为一系列滤波器组，将一个信号分解为若干子信号。尺度空间[7]表示作为一种识别EWT（Empirical Wavelet Transform，经验小波变换）边界的无参数方法，由于其自适应性而被广泛关注，但是相对来说，基于尺度空间的自适应经验小波变换速度慢且不准确，且无效分量较多。

针对上述情况，本文提出基于频谱趋势的自适应划分频带的方法，并将其与峭度指标结合，构造出自适应谱峭度图算法，并通过仿真信号及实测信号证明了该方法的有效性。

1 基本理论

1.1 谱峭度

Dwyer首先将峭度应用于STFT的实部和虚部，并因此引入了频域峭度（Frequency Domain Kurtosis，FDK）的概念[8]。SK最初定义为其频率分量的峭度[9]，并与不同光谱频率的振幅变化进行比较，因此该统计参数表示信号的脉冲性如何随频率变化。相反，Antoni[10]基于Wold-Cramer分解定义了SK，该分解将任何随机非平稳过程描述为因果、线性和时变系统的输出，即：

Antoni已经证明非平稳过程具有非高斯性是其基本特性[10]，然后将SK清楚地表示为CNS（Conditionally Non-Stationary，条件非平稳）过程的能量归一化四阶谱累积量，对应的表达式为：

1.2 快速谱峭度图

Antoni和Randall[2]提出了快速谱峭度图来表示频谱分割平面上信号的峭度值，它可以被认为是中心频率和带宽的不同尺度下峭度值的函数。快速谱峭度的“快”体现在其快速频带分割技术中，该技术使用二元1/3树结构将频带细分为多个合理的宽度。其实现过程如下：

步骤1：在不同层次上分割频带，然后形成一个具有不同带宽和中心频率的窄带二维平面。

步骤2：过滤掉每个窄带信号，并计算过滤信号的峭度值。峭度图是使用这些峭度值填充二维平面生成。

步骤3：峭度值最大的频带位于峭度图中作为最佳共振带，如图1所示，同时可以获得所选共振带的中心频率和带宽。

步骤4：对提取的滤波信号进行包络分析，用于轴承故障诊断。

图1 快速谱峭度图示例结果

2 自适应谱峭度图

2.1 频谱趋势自适应划分频带

许多方法都可以估计给定信号的趋势。传统的方法有线性或多项式回归分析、移动平均、基于趋势估计的傅里叶滤波等，它们的共同局限性是需要提前给出基函数或跨度，且只能在给定情况下有效地估计它们。Wu等[11]提出的基于EMD（Empirical Mode Decomposition，经验模态分解）的趋势是一种数据驱动的、完全自适应的趋势提取方法，首先对信号进行EMD分解，然后将最后几个IMF和残差进行合并得到EMD趋势。EMD可以将复杂信号自适应地分解多个固有模态函数（Intrinsic Mode Function，IMF），但存在严重的模态混叠问题[12]。Torres等[13]提出了一种自适应高斯噪声的集合经验模态分解方法（Complete Ensemble Empirical Mode Decomposition with Adaptive Noise，CEEMDAN），进一步提高了分解信号的准确性和完整性，有效地克服了模态混叠问题。

基于上述分析，本文采用大程度减少模态混叠和不必要迭代筛选且数据分解更完整的自适应噪声的集合经验模态分解（CEEMDAN）来完成频谱趋势的寻找。

基于频谱趋势的自适应频带划分的具体步骤如下：

（1）利用FFT（Fast Fourier Transform，快速傅里叶变换）获得输入信号的频谱；

（2）应用CEEMDAN得到频谱的IMF分量以及残差项；

（3）将最后一项的残差项作为第一项，剩余的IMF按照频率成分从低到高排列，依次从第一项残余项进行累加计算，每累加一个IMF作为一次频谱趋势，即每层频谱趋势；

（4）寻求每次频谱趋势的极小值位置，使其作为边界点，分割输入信号的频谱。

图2显示了原始信号的时域波形和频谱，对频谱进行CEEMDAN分解的结果如图3所示，图4展示了部分频带划分结果，其中蓝色曲线表示残余项与不同个数IMF叠加的结果、红色曲线表示为频谱边界划分的位置。

2.2 自适应谱峭度图

针对快速谱峭度图固定划分所带来的问题，本文提出基于频谱趋势的频带划分方法和峭度结合求解谐振频率和故障频率的方法——自适应谱峭度图，算法流程如图5所示。

自适应谱峭度图的轴承故障诊断的具体步骤如下：

（1）获取轴箱上固定的振动加速度传感器采集得到的振动加速度数据；

（2）计算相对应信号的傅里叶频谱，并对频谱进行自适应噪声的集合经验模态分解，得到各IMF分量及残余项；

（3）将残余项作为第一项，并按照频率成分升序排列IMF分量，依次叠加组合残余项与低频IMF分量，得到每一层相对应的频谱趋势；

图2 振动信号时域波形及频谱

图3 CEEMDAN分解结果

（4）计算每一层频谱趋势的局部极小值作为频谱分割边界，完成频谱划分；

（5）将划分后的各频段通过正交滤波器；

（6）计算每一层各频段对应的峭度值，并将其铺成二维平面，得到自适应谱峭度图；

（7）选取峭度值最大的分量，进行Hilbert解调，进行故障分析，当包络谱中包含有可识别的故障特征频率及其倍频时，则可认为该轴承存在相应的局部故障。

图4 频谱划分结果

3 仿真信号分析

为验证本文所提方法的有效性，构造一组单故障轴承仿真信号[14-15]如式（4）所示。

轴承故障仿真信号的参数如表1所示。信号采样频率为10 kHz，采样时间为1 s。同时将信噪比为-10 dB的高斯白噪声添加在仿真信号中。

仿真信号时域图、频谱图、包络谱如图6所示。从图6（c）可以看出，信号的中心频率位于2500 Hz左右，共振带长大衣位于2000～3000 Hz之间，与设置的仿真信号结果一致。从图6（d）可以看到包含故障噪声的故障特征频率及其倍频，也从侧面印证了通过包络谱可以更为清晰地了解故障情况。

图5 自适应谱峭度图算法流程图

式中：A为第个故障脉冲的振幅；为阻尼比；为真实的故障脉冲周期，＝F/，F为采样频率，为轴承故障特征频率；为随机滑移效应；()为单位阶跃函数；f为共振频率。

表1 轴承故障仿真信号参数

采用基于频谱趋势的自适应谱峭度图法处理该信号，结果如图7所示。从图7（a）可以看到，随着IMF低频分量的增加，每一层极小值点的个数也随之增加，频谱被划分的数量也随之增加。其中峭度值最大的频带位于层数4，中心频率为2422 Hz，带宽为834 Hz，提取此分量的波形如图7（b）所示，对此分量做包络分析，可以看到图7（c）包络谱中包含特征频率，且可以找到2倍、3倍、4倍、5倍等高倍频，故障信息明显，且中心频率f为2422 Hz，与实际中心频率非常接近。考虑到噪声的影响，该误差可以接受。因此，该方法适用于含噪声及周期性冲击信号。

对上述仿真信号进行快速谱峭度图计算，结果如图8所示。从图8（a）可看到峭度值最大的频段位于层数4，中心频率为2343.75 Hz，带宽为312.5 Hz，与实际中心频率相差接近150 Hz。提取此分量，波形如图8（b）所示，对此分量做包络分析，从图8（c）包络谱中可以看到特征频率及2倍频、3倍频，说明快速谱峭度对于轴承故障诊断具有一定的效果。

通过两种方法的对比分析，可以看到虽然两者都能够诊断出故障特征，但从整体上来说，通过使用基于频谱趋势的自适应谱峭度图可以更好地反映出仿真信号信息，即计算所得的中心频率及带宽与实际仿真信号中心频率更为接近，且幅值整体更为突显，也印证了快速谱峭度提因遵循频带固定划分，从而导致得到的最优频带包含不够丰富的故障信息。

4 试验分析

为进一步验证该基于频谱趋势的自适应谱峭度图的有效性，本节采用试验数据对该方法进行验证。

图6 仿真信号时域图、频谱图、包络谱

图7 采用基于频谱趋势的自适应谱峭度图法处理的结果（一）

图8 对仿真信号进行快速谱峭度图计算的结果（二）

试验台运行的轴箱来自高速列车车辆，试验数据来自台架试验轴箱轴承，该轴承为双列圆锥滚子轴承，轴承型号为CEI-2962，其几何尺寸如表2所示，试验台故障类型及位置如图9所示，故障为人为设置，其中外圈损坏的深度为1 mm，宽度分别为1 mm、3 mm、5 mm，采样频率为10 kHz，轴旋转频率f为5.1 Hz，采样周期为10 s。轴承外圈对应的特征频率分别为83.33 Hz。

表2 双列圆锥滚子轴承的基本参数

图10为采集外圈故障信号的时域波形及其傅里叶谱和包络谱。从频谱中可以发现信号有较强的谐波干扰，干扰频率为50 Hz，同时无法从振动信号的包络谱中看到轴承外圈故障特征频率，因此，需对信号进行处理，提取到故障特征信息，从而判断故障类型。

使用基于频谱趋势的自适应谱峭度图处理该信号，结果如图11所示。从图11（a）中可以看到峭度值最大的频带位于层数2，提取此分量的波形图如图11（b）所示，对该频带内的信号进行分析，可以看到图11（c）中包络谱包含特征频率，且可以找到2倍、3倍、4倍等倍频，故障信息明显且丰富，精确地检测到实际的故障特征。

当使用快速谱峭度图处理上述信号时，所得结果如图12所示，图12（a）显示了最大峭度值出现在中心频率625 Hz、带宽416.6667 Hz的频带中，滤波信号包络和包络谱分别如图12（b）和（c）所示，与原始信号的包络谱类似，从该包络谱中无法检测到任何轴承故障信息。因此，快速谱峭度图无法检测到所需的频带，很难确定轴承是否存在外圈故障，效果较差。

图9 试验台及故障类型和位置

图10 外圈故障信号的时域波形及其傅里叶谱和包络谱

5 结论

本文针对快速谱峭度图中采用固定划分容易使检测出的频带不能包含完整的故障信息或是提取到无效分量的问题，提出使用频谱趋势引导故障共振频带探测的方法，构造出基于频谱趋势的自适应谱峭度图。

从频谱趋势理论给出了具体分析过程，可以自适应地得到谱峭度图中每个节点的频带划分，并分别使用本文所提方法与快速谱峭度图方法对轴承故障仿真信号与试验信号进行了验证，表明采用频谱趋势划分频带的方法可以精确地对共振频带进行划分。

图11 使用基于频谱趋势的自适应谱峭度图法处理的结果（二）

图12 对仿真信号进行快速谱峭度图计算的结果（二）

[1]谷玉海，朱腾腾，饶文军，等. 基于EMD二值化图像和CNN的滚动轴承故障诊断[J]. 振动，测试与诊断，2021，41（1）：105-113.

[2]DWYER R F. Fourth-order spectra of Gaussian amplitude-modulated sinusoids[J]. Journal of the Acoustical Society of America，1991，90（2）：918-926.

[3]Antoni J. Randall R B. The spectral kurtosis：Application to the vibratory surveillance and diagnostics of rotating machines[J]. Mechanical Systems & Signal Processing，2006，20（2）：308-331.

[4]Antoni J. Fast computation of the kurtogram for the detection of transient faults[J]. Mechanical Systems & Signal Processing，2007，21（1）：108-124.

[5]T. Barszcz，A. JabŁoński. A novel method for the optimal band selection for vibration signal demodulation and comparison with the Kurtogram[J]. Mechanical Systems and Signal Processing，2011，25（1）：431-451.

[6]J. Gilles. Empirical Wavelet Transform[J]. IEEE Trans. Signal Process，2013，61（16）：3999-4010.

[7]J. Gilles，K. Heal. A parameterless scale-space approach to find meaningful modes in histograms - Application to image and spectrum segmentation[J]. International Journal of Wavelets，Multiresolution and Information Processing，2014，12（6）：1 450 044.

[8]Feng Z P，Liang M，Chu F L. Recent advances in time-frequency analysis methods for machinery fault diagnosis: A review with application examples[J]. Mechanical Systems and Signal Processing，2013，38（1）：165-205.

[9]Dwyer R. Detection of non-Gaussian signals by frequency domain kurtosis estimation[C]. ICASSP'83，IEEE International Conference on Acoustics, Speech, and Signal Processing，IEEE，1983（8）：607-610.

[10]Antoni N J. The spectral kurtosi: A useful tool for characterizing non-stationary signals[J]. Mechanical Systems & Signal Processing，2006，20（2）：282-307.

[11]Wu Z，Huang NE，Long SR，et al. On the trend, detrending, and variability of nonlinear and nonstationary timeseries[J]. Proceedings of the National Academy of Sciences of the United States of America，2007（104）：14889-14894.

[12]孙永鹏，董增寿. 改进EMD在轴承故障诊断中的应用[J]. 太原科技大学学报，2018，39（1）：25-30.

[13]Torres M E，Colominas M A，Schlotthauer G，et al. A complete ensemble empirical mode decomposition with adaptive noise[C]. 2011 IEEE international conference on acoustics, speech and signal processing (ICASSP)，IEEE，2011：4144-4147.

[14]Bard J. A Review of: “Engineering Optimization: Theory and Practice” [J]. Iie Transactions，1997，29（9）：802-803.

[15]McDonald G L，Zhao Q，Zuo M J. Maximum correlated Kurtosis deconvolution and application on gear tooth chip fault detection[J]. Mechanical Systems and Signal Processing，2012（33）：237-255.

Adaptive Spectral Kurtosis Map and Its Application in Fault Diagnosis of Rolling Bearings

HE Qian，YI Cai，WU Wenyi，LIN Jianhui

(State Key Laboratory of Traction Power, Southwest Jiaotong University, Chengdu 610031, China )

A new adaptive spectral kurtosis map method is proposed in this paper. First, the CEEMDAN methods are applied to decompose signal Fourier spectrum. Then the low-frequency IMF component and residuals are obtained and performed layered incremental calculation. The minimum points of frequency spectrum trend of each layer are found and made as spectrum division boundaries. The divided frequency bands pass through orthogonal spectrum filter and the fault information in the signal is evaluated by the kurtosis index to construct a new adaptive spectrum kurtosis map. This method can overcome the problem that the detected frequency band may not contain complete fault information or may extract invalid components by following fixed division rules. And the effectiveness of the method was verified by using the simulation signal and the test signal of bearing outer ring failure respectively.

fault diagnosis；adaptive spectral kurtosis map；spectral trend；fast spectral kurtosis map

TH133.33；TN911.7

A

10.3969/j.issn.1006-0316.2022.11.002

1006-0316 (2022) 11-0009-10

2022-02-25

国家自然科学基金（51905453）；中国博士后面上基金（2019M663899XB）；中央高校基本科研业务费项目（2682020CX50）；牵引动力国家重点实验室自主课题（2020TPL-T14）

何倩（1996－），女，四川南充人，硕士研究生，主要研究方向为机械设备故障诊断与信号分析，E-mail：h243052028@163.com。