基于谐波匹配补偿和无键相阶次跟踪的轴承故障诊断

武杰， 卢振连，， 马洪儒， 朱艳芳， 吴耀春， 薛晓峰， 姜阔胜

（1. 安阳工学院 机械工程学院，河南 安阳 455000；
2. 安徽理工大学 机械工程学院，安徽 淮南 232000；
3. 运城学院 机电工程系，山西 运城 044000）

轴承作为煤矿机械设备关键运动部件，直接影响设备的健康状态[1-2]。因此，对煤矿机械设备轴承进行故障诊断具有重要意义。传统的轴承故障诊断方法主要分为经典的时域统计分析方法和以傅里叶理论为核心的全局域变换方法。由于煤矿机械设备运行过程中振动大、冲击强，轴承振动信号呈现出瞬态非平稳的特性[3-4]，传统的轴承故障诊断方法难以识别故障特征。

传统阶次跟踪方法是通过有键相硬件设备获取非平稳信号的瞬时频率，再计算出瞬时相位[5-7]，进一步通过重采样获取角度域的循环平稳信号，从而实现非平稳信号到循环平稳信号的转变。但煤矿机械设备工作环境恶劣，不利于有键相硬件设备的安装，导致瞬时频率获取困难，阶次跟踪难以进行，严重影响了对轴承的故障诊断。无键相阶次跟踪方法可通过振动信号直接获得瞬时频率[8-11]，克服了对有键相硬件设备的依赖，能够很好地实现对轴承故障特征的提取。但该方法运用的前提是通过振动信号获取准确的瞬时频率，进而计算出瞬时相位。然而，煤矿机械设备振动剧烈，转速呈现出强烈的非线性特性，传统的基于短时傅里叶变换（Short Time Fourier Transform，STFT）的时频分析方法在估计该工况下的瞬时频率时鲁棒性较差。为解决该问题，文献[12]提出了一种基于Chirplet 变换的瞬时频率估计方法，该方法利用Chirplet 变换来逼近信号，从而实现对瞬时频率的细化，在一定程度上提高了瞬时频率的估计精度，但由于振动信号的非线性调制，在速度波动剧烈的情况下，时频分辨率较低，导致计算的瞬时相位精度不高[13-14]。文献[15]针对剧烈的非线性调制，提出了瞬时故障特征频率的概念，从轴承包络信号的时频谱中提取出瞬时频率，但在转速急剧变化的情况下，包络信号的转频由于能量较弱，无法准确估计。

本文提出了一种基于谐波匹配补偿（Harmonic Matching Compensation，HMC）的时频分析方法，可实现对轴承剧烈非平稳运行工况下瞬时频率的准确估计；
进一步采用无键相阶次跟踪方法，实现对轴承故障的准确识别。

基于HMC 和无键相阶次跟踪的轴承故障诊断流程如图1 所示。首先，采用基于HMC 的时频分析方法对振动信号进行处理，准确估计出瞬时频率。其次，由于在瞬态非平稳工况下，信号中混叠的各阶谐波分量随转速一起波动，通过Vold-Kalman 滤波方法提取谐波分量信号。再次，计算出瞬时相位，进而实现信号的角度域重采样。最后，对重采样的信号进行快速傅里叶变换，从而进行阶次谱分析，识别出轴承故障特征。

图1 基于HMC 和无键相阶次跟踪的轴承故障诊断流程Fig. 1 Bearing fault diagnosis process based on harmonic matching compensation and keyless phase order tracking

1.1 基于HMC 的瞬时频率估计

在瞬态非平稳工况下，对于一个有限的能量信号x(t)，其HMC 被定义为

式中：t，f分别为时间和频率；
S(t)为瞬时频率的解调函数；
hT(t)为 窗函数，其中心为τ、宽度为T，为了便于区分信号中的谐波分量，选择具有低混淆特性的汉宁窗。

对式（1）进一步变形，可得

在HMC 计算过程中，对于有限能量信号x(t)，先通过STFT 对瞬时频率进行粗略估计，之后采用泰勒公式将粗略估计的瞬时频率f(t)展开：

式中：n为阶数；
Rn(t)为泰勒余项。

在式（3）中，余项的大小取决于余项中t的阶次。当t处于高阶次时，余项值不断减小，同时无限逼近瞬时频率f(t)。因此，可以用关于t的负指数幂的形式来近似余项，同时通过调节参数c实现对瞬时频率的匹配：

式中：f0为幂指数函数的常数项；
d为指数因子。

由式（4）可知，粗略估计出的瞬时频率曲线为指数幂形式。因此，瞬时频率的解调函数S(t)可以被定义为

则式（2）中解调操作因子 e xp(−j2πS(t))可表示为

为了方便计算，d可设定为0.5。对于和c，可通过将瞬时频率曲线上任意2 个时频点 (t1，f1)和(t2，f2)代入式（5）进行求解。因此，可在每个汉宁窗中，分别获得相应的解调函数S(t)和解调操作因子exp(−j2πS(t))。

将 解 调 函 数S(t)和 解 调 操 作 因 子exp(−j2πS(t))分别转换成向量形式和。因此，经过HMC 变换解调后的信号为

式中 real[·]为实数运算操作。

1.2 基于Vold-Kalman 滤波的谐波分量提取

Vold-Kalman 滤波器是一种可有效跟踪信号波动变化的时变滤波器，能跟随信号特性自适应分离和提取谐波[16-17]。相对于传统的带通滤波器，Vold-Kalman 滤波器的中心带宽可以跟随转速的波动自适应调整，这种性质使其很适用于对煤矿机械设备轴承信号的跟踪滤波。在提取第k次谐波分量的同时，为有效抑制其他相邻谐波分量，Vold-Kalman 滤波器的瞬时带宽Ck需满足如下条件：

式中Ok为第k次谐波和其相邻谐波之间的距离。

1.3 瞬时相位计算

当第k次谐波分量e(n)的时域信号从原始振动信号中被提取出来后，采用Hilbert 变换计算第k次谐波的瞬时相位：

式中：
unwrap[·]为 相位解卷积操作，可将e(n)的瞬时相 位 从 区 间 [−π，π]映 射 到 相 位 累 积 空 间；
(n)为e(n)经Hilbert 变换后的结果。

理论上，第k次谐波信号的瞬时相位 φk(n)与原始信号瞬时相位 φ(n)之间严格遵循线性比例关系。但信号在实际的传递过程中，受振动传递路径的影响，瞬时相位会发生偏移。因此，式（9）可修正为

式中：a为 φk(n) 和 φ(n)之 间的线性比例系数；

φs(n)为由振动传递路径引入的相位偏移。

由于振动传递路径引入的相位偏移较大，需要对相位进行修正。修正后的实际瞬时相位为

1.4 信号角度域重采样

在角度域重采样过程中，为了避免产生频率混淆， ∆θ需要满足如下约束：

式中fs为信号采样频率。

通过时间域和角度域的映射关系，以 ∆θ为采样间隔，可实现原始时间域信号在角度域的等角度重采样，从而获得角度域重采样信号。

根据煤矿机械设备轴承的实际瞬态非平稳运行工况，设置轴承振动仿真信号[18-19]：

假设轴承外圈固定在轴承座中，内圈随转轴一起旋转，且在轴承外圈上存在1 处局部剥落缺陷。因此，轴承的故障脉冲信号s(t)可设置为[19]

假设轴承故障信号中包含2 个关于转频的谐波分量，其幅值分别为B1=0.005g（g为重力加速度），B2=0.01g， 初始相位分别为 φ1=π/6， φ2=−π/3。在测试过程中轴承经历了剧烈的转速上升和转速下降过程，因此瞬时频率[18]为

将轴承故障特征频率对外圈转频归一化处理[20]，设轴承的外圈故障特征频率为转轴转频的4.35 倍，则故障脉冲的平均角度间隔为360°/4.35≈82.76°。设置信号的采样频率为20 kHz，采样持续时间为2 s，噪声干扰为−3 dB 的白噪声。

轴承瞬时转速模拟曲线如图2 所示，可模拟煤矿机械设备作业过程中转速剧烈波动工况。

图2 轴承瞬时转速模拟曲线Fig. 2 Simulation curve of bearing instantaneous speed

仿真信号的时域波形如图3 所示。由于转速剧烈波动，各谐波分量及噪声干扰混淆，难以发现由轴承缺陷所激发的故障冲击特征。此外，受轴承转速剧烈波动影响，传统的时域和频域统计指标分析不再适用。

图3 仿真信号时域波形Fig. 3 Time domain waveform of simulated signal

采用本文方法对仿真信号进行处理。首先，通过快速谱峭度分析法[21]获得仿真信号的谱峭度，如图4 所示，进而识别出轴承的最佳共振频带。然后，根据最佳共振频带设计中心频率为1 875 Hz、带宽为416.67 Hz 的带通滤波器，提取出仿真信号的共振解调序列，如图5（a）所示。最后，采用Hilbert 变换得到仿真信号的包络，如图5（b）所示。

图4 仿真信号谱峭度Fig. 4 Spectral kurtosis of simulated signal

图5 仿真信号的共振解调序列和包络Fig. 5 Resonance demodulation sequence and envelope of simulated signal

为验证基于HMC 的时频分析方法在转速剧烈波动工况下对瞬时频率估计的有效性，分别采用HMC 方法和STFT 方法对瞬时频率进行估计，得到的时频谱如图6 所示。可看出采用HMC 方法得到的时频谱中时频聚集性显著且分辨率高，而采用STFT 方法得到的时频谱在剧烈转速波动和噪声双重干扰下较模糊，且时频聚集性和分辨率较差。

图6（a）中有2 个能量突出的谐波分量，对谐波分量Ⅰ进行瞬时频率估计，结果如图7 所示。可看出估计的瞬时频率曲线与实际的瞬时频率曲线（通过式（16）计算所得）高度重合，瞬时频率估计值与实际值之间的最大相对误差不超过1%。

图6 不同方法下仿真信号时频谱Fig. 6 Time-frequency spectrum of simulated signal under different methods

图7 仿真信号瞬时频率估计结果Fig. 7 Instantaneous frequency estimation result of simulated signal

以HMC 方法精确估计出的瞬时频率为中心频率，利用Vold-Kalman 滤波器对谐波分量Ⅱ进行自适应提取，结果如图8 所示。

图8 Vold-Kalman 滤波器提取的谐波分量Fig. 8 Harmonic components extracted by Vold-Kalman filter

从图8 可看，出提取出的谐波分量波形平滑，幅值为0.01g，与设置的仿真信号幅值B2=0.01g相同。根据式（8），取Vold-Kalman 滤波器的瞬时带宽为 0.5阶次。

依据式（9）和式（11），估计出瞬时相位，并与实际瞬时相位对比，结果如图9 所示。可看出估计的瞬时相位曲线与实际瞬时相位曲线高度重合，瞬时相位估计值与实际值的最大相对误差不超过2%。

图9 仿真信号瞬时相位估计结果Fig. 9 Instantaneous phase estimation result of simulated signal

根据估计出的瞬时相位，通过式（12）获得时间域和角度域的映射关系，进而完成角度域的重采样。再将重采样的角度域信号进行快速傅里叶变换，得到包络阶次谱，如图10 所示。可看出本文方法和有键相阶次跟踪方法得到的包络阶次谱高度吻合，通过本文方法可清晰找到轴承外圈故障特征阶次（4.35 阶次）及相关倍频特征阶次，证明轴承外圈存在故障。

图10 仿真信号包络阶次谱Fig. 10 Envelope order spectrum of simulated signal

为验证本文方法对煤矿机械设备轴承故障诊断的有效性，搭建轴承模拟试验系统，如图11 所示。

图11 轴承模拟试验系统Fig. 11 Bearing simulation test system

将轴承外圈通过过盈配合固定于轴承座内，轴承内圈与轴一起旋转。电动机与伺服驱动系统提供动力输入，驱动轴承转动。采用3 个灵敏度为50 mV/g的单轴振动加速度传感器采集轴承振动信号。转速计为光电脉冲编码器，用来采集轴承转频信息。信号采集系统主要设备为Coco−80 数字采集仪。为了模拟煤矿机械设备运行时轴承运行工况，设置轴承经历瞬间启动、载荷及转速迅速上升和瞬间停止3 个过程。设置采样频率为25.6 kHz。

轴承几何参数见表1。试验设置轴承外圈存在人为设置的磨损故障缺陷，如图12 所示。

表1 轴承几何参数Table 1 Bearing geometry parameters

图12 轴承外圈故障Fig. 12 Outer ring fault of bearing

通过归一化操作[20]可得轴承外圈故障特征阶次：

式中：z为滚动体个数；
fr为外圈转频；
r为节圆半径；
D为滚动体直径；

ϕ为接触角。

将表1 数据代入式（17），计算出轴承外圈故障特征阶次为3.072 7。

采集的原始振动信号、转速和截取的试验信号如图13 所示。由于转速的剧烈瞬态波动，从图13（a）和图13（c）中无法识别出任何关于轴承故障的冲击特征。

图13 试验采集信号Fig. 13 Test acquisition signal

通过快速谱峭度分析法获得截取的试验信号的谱峭度，如图14 所示，进而可识别出轴承的最佳共振频带。根据获得的最佳共振频带设计中心频率为1 718.8 Hz、带宽为312.5 Hz 的带通滤波器，进而提取出截取的试验信号的共振解调序列，并采用Hilbert 变换得到截取的试验信号的包络，如图15所示。

图14 试验信号谱峭度Fig. 14 Spectral kurtosis of test signal

图15 试验信号的共振解调序列和包络Fig. 15 Resonance demodulation sequence and envelope of test signal

分别采用HMC 方法和STFT 方法处理试验信号，得到时频谱，如图16 所示。可看出采用STFT 方法处理后的时频谱中识别不出任何频率波动的趋势，而由于HMC 方法在瞬态条件下的优异性能，可清楚识别出瞬时频率的变化趋势。

图16（a）中有2 个谐波能量较为显著。为了确定这2 个谐波分量的位置，在0.058 38 s 处对其进行检测，其中心频率分别为 975.3 Hz 和953.6 Hz，表明2 个谐波分量之间的频率间隔为21.7 Hz。因此，确定2 个谐波分量分别为试验信号的第45（975.3 Hz/21.7 Hz≈45）个 和 第44（953.6 Hz/21.7 Hz≈44）个 分量。本文通过第45 个谐波分量估计试验信号的瞬时频率，结果如图17 所示。可看出瞬时频率估计曲线与实际的瞬时频率曲线高度重合，瞬时频率估计值与实际值之间的最大相对误差不超过1%。

图16 不同方法下试验信号时频谱Fig. 16 Time-frequency spectrum of test signal under different methods

图17 试验信号瞬时频率估计结果Fig. 17 Instantaneous frequency estimation result of test signal

采用HMC 方法准确估计出瞬时频率后，通过Vold-Kalman 滤波器对第45 个谐波分量进行自适应提取。根据式（9）和式（10），估计出瞬时相位，并与实际瞬时相位对比，结果如图18 所示。可看出估计的瞬时相位曲线与实际的瞬时相位曲线高度重合，瞬时相位估计值与实际值之间的最大相对误差不超过2%。

图18 试验信号瞬时相位估计结果Fig. 18 Instantaneous phase estimation result of test signal

根据估计出的瞬时相位，通过式（12）获得时间域和角度域的映射关系，进而完成角度域重采样。对重采样的包络信号进行快速傅里叶变换，获得包络阶次谱，如图19 所示。可看出本文方法和有键相阶次跟踪方法得到的包络阶次谱高度吻合，通过无键相阶次跟踪方法可清晰找到轴承外圈故障特征阶次（3.072 7 阶次）及相关倍频特征阶次，由此可明确推断出轴承外圈存在故障。

图19 试验信号包络阶次谱Fig. 19 Envelope order spectrum of test signal

为进一步验证本文方法的有效性，对轴承进行内圈故障加工，如图20 所示。应用本文方法得到轴承内圈故障振动信号的包络阶次谱，如图21 所示。可看出本文方法和有键相阶次跟踪方法所获得的故障特征阶次相同，但本文方法获得的故障特征阶次幅值明显高于有键相阶次跟踪方法获得的故障特征阶次幅值，更有利于识别故障特征。

图20 轴承内圈故障Fig. 20 Inner ring fault of bearing

图21 轴承内圈故障振动信号的包络阶次谱Fig. 21 Envelope order spectrum of vibration signal of bearing inner ring fault

1） 基于HMC 的时频分析方法不依赖传统的局部平稳假设要求，可自适应地匹配调整瞬时频率信号，并对瞬时频率信号进行自适应补偿，进而降低了信号的非线性度，避免了传统基于STFT 的时频分析方法在瞬态非平稳运行工况下时频谱模糊的问题，从而精确估计出瞬时频率。

2） 将基于HMC 的时频分析方法与无键相阶次 跟踪相结合，应用于轴承故障诊断，能够准确估计瞬时频率，获得较为精确的瞬时相位，通过包络阶次谱可准确且明显地表征出轴承故障特征阶次，从而有效识别故障特征。