基于谱峭度的滚动轴承故障包络阶比跟踪分析

刘亭伟，郭 瑜，李 斌，高 艳

（昆明理工大学 机电工程学院，昆明 650093）

旋转机械设备的变速运行工况是一种普遍现象，机械变速运行过程中的振动信号往往包含有更为丰富的振动特征信息［1］，但该工况下的振动信号是非平稳的，一般具有幅值、频率时变性，其振动特征通常也是时变的，因此建立在稳速运行工况下的传统分析方法一般不易在变速运行过程中实现对故障的有效诊断，这已成为目前旋转机械故障诊断的瓶颈［2］。同时轴承故障初始阶段引起的振动通常为微弱振动信号，机械系统中大量干扰源的存在常使故障信号被湮没在强干扰等背景噪声中。如何在变速运行工况下对干扰进行有效抑制和在噪声环境中对微弱故障信号进行提取是目前机械故障诊断学科中的重要研究内容之一。

包络分析（又称为共振解调分析）目前已被广泛应用于滚动轴承的故障诊断，它可以有效地提取出滚动轴承故障初期的微弱冲击信号。包络分析中的关键在于滤波频率和带宽（共振区）的确定。传统包络分析方法一般是根据经验选择，降低了包络分析的灵活性和准确性。近年来，基于谱峭度的包络分析方法［3-5］被提出，它可以根据谱峭度值自适应地确定共振中心频率和滤波频带，进而通过共振解调算法把微弱的冲击振动从低频干扰和强背景噪声中提取出来。在近来的研究中，该方法通常是对旋转机械稳速运行工况的振动信号，该方法是否适用于变速运行工况下的振动信号分析是一值得进一步探讨的问题。在研究中发现，虽然旋转机械轴承故障引起的微弱冲击序列在变速运行过程中丧失了其在平稳运行下的周期性，但该冲击序列的冲击性本质并无改变，即持续续时间短，具有很宽的频带，其带宽涵盖了轴承及相关结构的共振频带并可引起相关的固有振动，因此包络提取应同样适合变速运行工况。但应注意的是直接对变速工况下的振动信号进行频谱分析会造成频率模糊现象［6］，以至于得不到故障特征，解决该问题的途径是使用针对旋转机械变速运行工况振动分析的阶比跟踪分析（Order Tracking Analysis）法［7］，其利用旋转机械振动信号与转速的阶比关系，通过等角度重采样方法将变速工况下的频变信号转化为阶比恒定的角度域信号，对等角度信号进行频谱分析，即可得到清晰的阶比谱，进而实现对故障特征的准确提取。

本文将谱峭度分析方法和阶比分析方法相结合，提出了基于谱峭度的包络阶比分析方法，用该方法成功实现了对旋转机械变速运行工况下的滚动轴承故障特征提取。

1 基于谱峭度的包络提取

峭度（Kurtosis）是描述波形尖峰度的一个量纲，对冲击信号敏感，可用于评价振动信号中冲击成分的强弱，但易受噪声干扰而无法表征冲击成分。近年来提出了谱峭度［4］（Spectral Kurtosis，SK）用其来反映原始信号在某个频率成分上峭度值大小。其计算公式表示为［3-5］：

谱峭度快速计算方法采用1/3-二叉树结构构建一系列带通滤波器组实现谱峭度的快速计算［5］。其主要计算步骤为：

（1）获取各级复包络。采用与小波包分解算法相似的处理方法。首先构建两个具有线性相位的准解析低通滤波器h0（n）和高通滤波器h1（n）：

式中：h（n）为截止频率fc=1/8+ε，ε≥0的标准低通FIR滤波器（此处为归一化的频率，即1代表采样频率fs，一般fc取 0.3［5］）；h0（n）为由h（n）频移 1/8 后的准解析低通滤波器，其带宽为［0，1/4］；h1（n）为由h（n）频移3/8后的准解析高通滤波器，其带宽为［1/4，1/2］。然后利用以上两个滤波器h0（n）和h1（n）对原信号进行滤波，2倍降采样，以此方式迭代，将原信号进行K级分解，把原信号分解成不同的子带信号（n），其中系数k=0，…，K-1为分解级数，系数i=0，…，2k-1为子代信号位置索引。分解算法可用以下递归形式实现：

式中：当k=0时，c0（n）≡x（n）。对于子代信号（n），可以认为它是原始信号x（t）经中心频率fi=（i+2-1）·2-k-1和带宽Δfk=2-k-1的窄带滤波器滤波后降采样得到的信号［5］。由式（2）不难看出滤波器的实部相当于一个零相移滤波器，虚部相当于一个90°相移滤波器，因此分解后得到的子代信号（n）的实部和虚部相位差为90°，根据包络检波原理，子代信号的模即为信号的包络（n）称为在中心频率fi=（i+2-1）2-k-1和带宽Δfk=2-k-1处的复包络信号。

为了提高分析的精度，该方法对频带做进一步细分，在k+1和k+2级之间插入3（2k个滤波器，和以上方法类似，通过构建三个准解析的带通滤波器gj（n），j=0，1，2，它们的频带分别为［0，1/6］、［1/6，1/3］和［1/3，1/2］，对每个复包络信号（n）进行分解，得到三个子代信号

（2）计算各复包络信号（n）的谱峭度。根据式（1）计算复包络信号（n）在中心频率fi和带宽Δfk处的谱峭度值。因此式（1）可改写为：

（3）获得最佳包络参数和最佳复包络信号。当式（4）谱峭度达到最大时，其对应的中心频率和带宽即为优化包络中心频率fo和带宽Δfo，其对应的复包络即为优化复包络信号co（n）。即：

式中；argmax表示取最大值的参数，即获得谱峭度值为最大时对应的优化中心频率fo和带宽Δfo，以及对应的优化复包络信号co（n）。

2 阶比分析简介

如上所述旋转机械变速运行工况下的振动信号为频率时变信号，直接对其进行频谱分析会产生所谓的频率模糊现象［6］，为解决该问题产生了阶比分析方法。阶比定义为每转的循环次数，即振动信号频率与对应转轴转频的倍数［8］，可用以下公式表示：

式中：f为振动信号的频率；n为对应转轴的转速，单位为r/min，则n/60为转频；l为阶次。

滚动轴承由内圈、外圈、滚动体和保持架四部分组成。滚动轴承故障特征频率可由文献［9］给出，将滚动轴承故障频率代入公式（6）可得滚动轴承的故障特征阶比，分别表示如下，

内圈故障特征阶比：

外圈故障特征阶比：

滚动体故障特征阶比：

保持架故障特征阶比：

式中：Z为滚动体的个数；d为滚动体直径；D为轴承节径；α为接触角；nb为轴承所在轴的转速；n为所选取的参考轴的转速。

阶比分析根据升降速阶段转速与故障频率的阶比对应关系，把时域等时间间隔（非平稳）采样信号通过等角度重采样变为角度域准平稳信号。阶比分析的关键在于如何实现等角度采样，即阶比跟踪采样。

研究中采用了基于转速计脉冲的计算阶比跟踪［10-11］技术实现等角度采样，其原理可简述为：① 计算阶比跟踪首先对原始振动信号和键相信号进行常规等时间间隔同步采样；② 利用键相脉冲信号进行转速估计，并作为振动相位的参考基准，并由转速估计得到等角度采样的键相时标；③ 由得到的等角度采样时标对原振动信号进行插值重采样，得到振动信号的等角度采样信号；④ 最后以等角度采样信号为数据通过FFT实现阶比谱分析，见文献［10-11］。

3 基于谱峭度的滚动轴承包络阶比分析

基于谱峭度的滚动轴承包络阶比分析原理框图如图1所示，其基本分析步骤为：

（1）基于谱峭度的包络提取：首先按式（3）对原始变速振动信号x（n）进行分解，得到不同中心频率和带宽下的复包络信号（n），然后根据式（4）计算各复包络信号的谱峭度值，最后根据式（5）得到优化的包络中心频率fo和带宽Δfo，及优化的复包络信号co（n），实现对信号包络的自适应提取，把湮没于强背景干扰中的微弱轴承故障振动信号提取出来。

图1 基于谱峭度的滚动轴承包络阶比分析流程图Fig.1 The envelope order analysis flowsheet of rolling element bearing based on SK

（2）等角度采样时标的获取：为了实现等角度采样必须获得等角度采样时标，计算阶比跟踪技术采用以下方式实现。一般情况下，在一段很小的时间间隔内可以认为转轴角加速度恒定［11］，可用一元二次方程来描述转角与时间的关系，即：

式中：a0、a1、a2为待定系数。把同步采集的键相脉冲信号s（t）连续的三个脉冲发生时刻t1、t2、t3和对应的转角0、ΔΦ、2ΔΦ（每圈一个脉冲时取 ΔΦ=2π）代入式（11），求解a0、a1、a2，用矩阵形式可表示为：

将求得的系数a0、a1、a2代入式（13）得到任意角度所对应的时刻t为：

对转角按等角度采样间隔离散化，上式可写为：

式中：Tn为第n点等角度采样时标；Δθ为等角度采样间隔。

（3）包络信号等角度重采样：等角度重采样是根据已知的离散时域振动信号，本研究中为第1步得到的优化包络信号co（n），在离散的等角度采样时标Tn上插值，得到各等角度时刻的振动幅值大小，最终得到等角度信号co（Tn）。插值重采样方法有多种，在高采样率条件下，最简单的线性插值即可取得较好的效果，在仿真和试验中运用拉格朗日线性插值［12］，其公式为：

式中：ni和ni+1为离Tn最近的两个相邻的等时间间隔采样时刻，co（n）为优化的包络信号。co（Tn）即为所求的等角度采样包络信号。

（4）故障特征分析：对得到的等角度信号co（Tn）进行FFT得到阶比谱，将其与理论计算获得的滚动轴承故障特征阶比进行对比，实现对滚动轴承故障的辩识。

4 仿真试验分析

为了验证本方法的有效性，进行了仿真试验，仿真过程如下：

仿真一个固有频率（系统固有频率）fr=7 000 Hz幅值为0.25的单一微弱冲击衰减振动用于仿真轴承故障的每转内的一次冲击：

仿真一阶数分别为1×、2×和3×的降速强低频干扰：

式中：f（τ）为降速过程的瞬时频率，其起始频率f0=25 Hz（对应起始转速为1 500 r/min），终止频率f1=5 Hz（对应终止转速为300 r/min），持续时间T为4 s，并有：

最后仿真一幅值为0.1的高斯白噪声干扰n（t）。则仿真的轴承振动信号为：

式中：R（t）为每转重复5次r（t）的变频（按照瞬时频率曲线变化）仿真故障振动信号，即仿真一滚动轴承冲击故障特征阶比为5×及其谐波的冲击信号。仿真中取采样频率fs=20 kHz，仿真的滚动轴承振动信号如图2所示（图中显示时间从3.1～3.3 s）不难发现冲击振动成分在整个信号中并不明显。

图2 滚动轴承振动仿真信号Fig.2 Simulated vibration signal of rolling elements bearing

在没经过基于谱峭度包络分析，直接对原始仿真信号进行阶比跟踪得到的阶比功率谱如图3（横坐标为阶次l）所示，在图3中只见到了仿真的强低频干扰信号y（t）的阶比，而没有见到仿真的故障振动信号R（t）的阶比。这说明故障振动信号非常微弱，被强的低频干扰信号所淹没。

图3 没经过包络分析处理的仿真信号阶比谱Fig.3 Order spectrum of simulated signal without envelope analysis

按论文介绍的方法对振动仿真信号进行分解（仿真中最大分解级数K=6），并计算在各中心频率fi和带宽Δfk处复包络信号cj k（n）的谱峭度值（参见第3节步骤1），得到的谱峭度图如图4所示。图4中横坐标为频率，纵坐标为尺度k，尺度不同带通滤波器的中心频率和带宽不同，尺度越大，带宽越窄，图中方块颜色深浅代表谱峭度值大小，灰度值越大谱峭度值越大。从图中可看出在中心频率fo为7 083 Hz，带宽Δfo为833 Hz时谱峭度值最大，此参数即为优化的共振解调参数，这与仿真的共振频率7 000 Hz相吻合。其对应的包络信号即为优化包络信号co（n），但由于其为变频冲击信号，故对co（n）进行计算阶比跟踪等角度重采样（参见第3节步骤2和3）得到等角度包络信号co（Tn），对co（Tn）进行FFT，得到的阶比功率谱如图5所示，对比图3我们可以发现在经过谱峭度的包络分析后，可清晰地得到仿真的5×故障特征阶比及其谐波，这说明强的低频干扰被抑制甚至被去除，微弱的故障信号被分离出来。仿真结果表明本方法的可行性。

5 测试试验分析

研究中以ZJS50综合设计型机械设计试验台为测试对象，如图6所示。

在中间齿轮箱输入轴位置处安装有一故障滚动轴承，进行滚动轴承故障试验，以外圈故障为例，试验参数如下：故障轴承型号为N1007；安装轴转速从1 200 r/min（对应转频为20 Hz）开始降速，以降速工况的振动信号为分析对象，按式（7）计算得滚动轴承外圈故障的特征阶比约为6.4×；数据采集设备为四通道NIUSB9215A采集卡；采样频率为20 kHz；加速度传感器型号为4530，灵敏度为69.7 pC/g，安装位置为轴承座；转速脉冲测量传感器为电涡流位移传感器，型号为4842，灵敏度为-3.645 V/mm。图7为采集的原始振动信号时间波形和转速脉冲信号（注：图中转速脉冲信号幅值单位：V，由于要同时显示两个通道的信号，图中只标出了振动信号的幅值单位）。

对测试振动信号不进行包络分析，只进行阶比跟踪得到的阶比功率谱如图8所示，图8中阶比功率谱错乱无章，找不到有用信息。

图8 未经过包络分析处理的测试数据阶比功率谱Fig.8 Order pectrum of the test data without envelope analysis

对测试振动信号进行基于谱峭度的包络提取（见第3节步骤1），得到的谱峭度图如图9所示，获得的优化共振解调频率为5 833 Hz，优化带宽为1 666 Hz。

图9 测试数据的谱峭度图Fig.9 Fast kurtogram of the test data

按照第3节步骤2和3对提取的优化包络时域信号进行计算阶比等角度重采样（取阶比跟踪触发转速为900 r/min，每转重采样点数为56，即最大分析阶比为28×，数据长度为2 048），对应的阶比功率谱如图10所示。

图10 经包络处理的测试数据阶比功率谱Fig.10 Order spectrum of the test data with envelope analysis

对比图8，从图10中可以清晰地看到6.272×的滚动轴承外圈故障特征阶比及其谐波，与滚动轴承外圈故障的理论特征阶比6.4×及其谐波基本吻合（一般滚动体和内外圈之间存在1～2%转频的滑动［13］，以及阶比谱分辨率等导致了实测故障频率和理论故障频率存在一定误差）。分析结果表明本方法对实际测试有效。

6 结论

在变速工况下，由于滚动轴承故障对应的冲击信号仍然可以激起轴承及其周围结构的固有振动，因此共振解调方法在此同样适用。运用基于谱峭度的包络提取方法可以根据谱峭度值自适应的确定共振解调的优化滤波相关参数，进而把微弱的滚动轴承故障冲击信号提取出来。阶比分析可通过等角度采样实现对旋转机械变速运行工况的有效分析。本文将以上两种方法相结合，实现了对变速运行工况下滚动轴承故障特征的提取。仿真和测试试验验证了本方法的有效性，对发展旋转机械变速运行工况下的滚动轴承故障特征提取和诊断方法有较好的研究意义。

［1］褚福磊，彭志科，冯志鹏，等.机械故障诊断中的现代信号处理方法［M］.北京：科学出版社，2009.

［2］李志农，丁启全，吴昭同，等.旋转机械升降速过程的双谱-FHMM识别方法［J］.振动工程学报，2003，16（2）：171-174.

［3］Sawalhi N，Randall R B，Endo H.The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis［J］.Mechanical Systems and Signal Processing，2007，21（6）：2616-2633.

［4］ Antoni J，Randall R B.The spectral kurtosis：application to the vibratory surveillance and diagnostics of rotating machines［J］.Mechanical Systems and Signal Processing，2006，20（2）：308-331.

［5］Antoni J.Fast computation of the kurtogram for the detection of transient faults［J］.Mechanical Systems and Signal Processing，2006，21（1）：108-124.

［6］郭 瑜，秦树人.基于瞬时频率估计及时频滤波的阶比分量提取［J］.中国机械工程，2003，14（17）：1506-1509.

［7］郭 瑜，秦树人.旋转机械非稳定信号的伪转速跟踪阶比分析［J］.振动与冲击，2004，23（1）：61-69.

［8］郭 瑜，秦树人，梁玉前.时频分析阶比跟踪技术［J］.重庆大学学报，2002，25（5）：17-21.

［9］陈 进.机械设备振动监测与故障诊断［M］.上海：上海交通大学出版社，1997.

［10］杨炯明，秦树人，季 忠.旋转机械阶比分析技术中阶比采样实现方式的研究［J］.中国机械工程，2005，16（3）：249-253.

［11］ Fyfe K R，Munck E D S.Analysis of computed order tracking［J］.Mechanical Systems and Signal Processing，1997，13（4）：667-641.

［12］郭 瑜，秦树人，汤宝平，等.基于瞬时频率估计的旋转机械阶比跟踪［J］.机械工程学报，2003，39（3）：32-36.

［13］ Sawalhi N，Randall R B.Simulating gear and bearing interactions in the presence of faults partⅠ.the combined gear bearing dynamic model and the simulation of localised bearing faults［J］. Mechanical Systems and Signal Processing，2008，22（8）：1924-1951.