基于线调频小波路径追踪和逐步解调滤波的滚动轴承故障诊断

2019-06-21 07:46:42刘东东程卫东温伟刚

振动与冲击 2019年11期

刘东东，程卫东，温伟刚

(北京交通大学机械与电子控制工程学院，北京 100044)

滚动轴承是旋转机械最关键的部件之一，由于其常处于变载荷﹑变转速的工况中，工作环境恶劣，极易发生故障[1-2]。阶比跟踪技术是最常用的处理时变工况下振动信号的方法[3]。阶比分析中鉴相时标的确定需要转速信息，转速可以通过转速计测取，也可以在信号中提取[4]。由于测取转速受到安装成本以及安装空间的限制，如何提取转速引起学者的研究。文献[5]提出了基于多尺度线调频基的稀疏信号分解方法，自适应地选取适当尺度对振动信号投影分解，提取齿轮包络信号的瞬时频率。但是处理滚动轴承振动信号时，高幅值的故障特征频率(Fault Characteristic Frequency, FCF)会干扰转频的提取。文献[6]使用FCF代替转频对信号重采样，根据故障特征阶比模板判断轴承故障。文献[7]根据估计转频的范围构造一组带通滤波器对特征频率进行滤波，根据依靠等效转频采样得到的故障相角域判断轴承状态。虽然根据等效转频进行阶比分析可以恢复信号的周期性，但是等角度重采样过程复杂，计算效率较低[8]。文献[9]也证明阶比分析存在包络畸变现象。

广义解调时频分析算法[10]可以将时频分布为倾斜、非线性的频率分量转换成平行于时间轴的频率。由于其非常适用于处理非稳态的调幅-调频信号，已经开始应用于振动信号的处理。文献[11]使用广义解调时频分析，将齿轮振动信号分解为若干个单分量信号，再对各个单分量信号进行包络分析，得到包络阶次谱。文献[12]通过对齿轮振动信号进行基于多尺度线调频基的稀疏信号分解，估计啮合频率，再对信号进行广义解调和FFT变换，得到齿轮的FCF。然而，在轴承的振动信号中，由于调制转频分量太小，因此很难直接提取，用于估计相位函数。文献[13]将广义解调算法应用到滚动轴承的故障诊断中，但是相位函数的估计以实测转速为基础，在实际工况中，会受到安装成本以及安装位置的制约，此算法的应用也会受到限制。

从以上分析可知，以阶比分析为基础的轴承故障诊断法会有效率以及误差等问题，而广义解调算法的使用会依赖转速计等辅助设备。当工况变化时，轴承故障引起的高频共振不仅会调制FCF，而且FCF也会调制转频。因此，理论上在包络时频谱(Time-Frequency Representation, TFR)中会有转频分量，但是，转频的幅值一般较小，很难直接提取[14]。瞬时故障特征频率(Instantaneous Fault Characteristic Frequency, IFCF)具有幅值优势[15]，一般容易提取。文献[16]提出了线调频小波路径追踪(Chirplet Path Pursuit, CPP)算法，该算法与常用提取频率的峰值搜索算法相比，抗噪声能力强，拟合精度较高。文献[17]使用CPP算法估计齿轮啮合倍频和轴承IFCF，判断含有齿轮振动下的轴承运行状态。文献[18]根据CPP算法估计的IFCF和转频对应关系，判断轴承故障位置，但是该算法以实测转速为参考。

本文利用CPP算法的优良特性，使用该算法在降采样的包络信号中估计IFCF趋势线，据此计算IFCF的相位函数。然而，仅仅依靠IFCF在解调频谱中的峰值，无法判断故障位置，甚至出现误判。为此，本文利用估计的IFCF，根据轴承各部分故障特征因子，采用重复估计转频的方式，计算潜在的转频。根据IFCF和转频的相位函数，构造逐步解调滤波算法，对信号各频率进行逐步解调和滤波。通过对仿真信号和试验信号的处理，证明了该算法的有效性。

1 算法部分

为了在不使用转速测量设备的情况下，对轴承故障位置进行辨别，应对轴承信号的特点进行分析，以提取必要瞬时频率，估计相位函数。若滚动轴承出现局部故障，故障点与配合面会形成一系列冲击，该冲击会激起系统的高频共振。转速不变时，冲击间隔相等，间隔的倒数就是故障特征频率[19]。理论上FCF只与转速、故障位置和轴承参数相关，外圈、内圈和滚珠存在故障时的故障特征频率fo,fi和fb的公式

(1)

(2)

(3)

式中:fr为转速，N为滚动体个数，D为节圆直径，d为滚动体直径，φ为接触角。

1.1 线调频小波路径追踪

CPP算法是近年提出的一种检测瞬时频率的算法, 它的原理是将信号划分为一系列的动态时间支撑区，从建立的Chirplet原子库中，按照相关性最大的原理，选择各支撑区间的原子，根据最佳路径连接原则对选定原子进行连接，依据各个支撑区间的原子的频率估计信号的瞬时频率。

建立Chirplet原子库

(4)

式中:I为动态时间支撑区间，I=[kN2-j,(k+1)N2-j]×Δt；k为区间的序号，k=0,1,…,2j-1；N为采样长度；j为区间尺度系数，j=0,1,…,log2N-5；1I(t)为矩形窗函数，即t∈I时，1I(t)=1；t∉I时，1I(t)=0；aμ为线调频小波的调频率系数；bμ为频率的偏置系数；按照采样定理，aμt+bμ应小于fs/2.

定义信号的chirplet变换为信号与原子库中所有原子的内积

(5)

式中：t=n×Δt，n为t时刻点数，Δt为采样时间间隔。

在原子库中挑选一组原子进行连接，使分析信号y(t)在时间支撑区I上的具有最大投影系数

(6)

设时间支撑区I内的最大投影系数kI表示的信号分量为dI(t)

dI(t)=|kI|exp[-i(at2/2+bt)]1I(t)

(7)

连接投影系数kI表示的信号分量dI(t)，使得信号d具有最大的总能量

(8)

式中：Π包含信号的整个时间长度，而且不重叠。

该算法的具体连接方法如下：

(1) 初始化，设置e(i)=0，pred(i)=0，其中，e(i)为选择的i-1个动态时间支撑区的信号的总能量，pred(i)为第i个时间支撑区的前置支撑区的序号。

(2) 对动态时间支撑区集合{Ii|i∈z}，寻找相邻的下一个动态支撑区的集合{Ik|k∈z}，如果e(i)+d(i)>e(k)，则有

(9)

该连接算法能保证Chirplet原子的组合与原信号的相似度最高。因此，CPP算法是通过连接各个支撑区的瞬时频率为线性的原子逼近信号的瞬时频率。

1.2 广义解调算法

由于转速的变化使轴承故障引起的冲击失去了周期性，信号的频谱出现频率模糊现象，因此，有必要将时变频率转化为周期频率。广义解调算法是一种新的处理时变信号的方法，可以根据相位函数将非周期的瞬时频率转换为周期频率。该算法的本质为广义傅里叶变换，对于信号x(t)，其广义傅里叶变换的定义为

(10)

式中：s0(t)为与时间t相关的实值函数，这实际是对x(t)e-2πjs0(t)做标准的傅里叶变换。

2 轴承的故障诊断方法

对于某一特定时变频率，广义解调算法可以根据该频率的相位函数将其转变为周期频率。因此，在频谱中该频率由模糊转化为突出的峰值。利用CPP算法可以在包络信号中估计IFCF，因此，根据估计的IFCF可以实现对其解调。但是仅仅依靠该峰值，无法判断轴承故障及其位置。这是因为转子失衡或者安装原因也会产生随转速变化的频率，如果利用CPP算法估计的瞬时频率为该频率，对其解调也会出现峰值。因此，如果判断该频率是轴承故障引起的，还要辨别故障位置，必须需要借助参考频率。

在变转速的工况下，高频共振不仅会调制FCF，而且FCF还会调制转频。由于转频的幅值较低，一般很难在信号中直接估计。假设估计的IFCF是外圈故障引起的，根据FCF与转频的比例fr=fcf/Co(Co为外圈故障特征系数)，可以计算潜在转频。如果在解调滤波后的转频位置不出现峰值，那么该故障就不是外圈故障，然后再假设内圈故障，对转频进行估计。为了减少FCF以及幅值较小的转频被噪声污染，可以对解调后的信号进行带通滤波，滤除背景噪声和解调造成的频率干扰。这里将根据CPP算法估计的轴承IFCF，依据特征系数计算其他频率，根据得到的频率对信号进行逐步解调，恢复各频率的周期性，然后通过带通滤波滤除噪声，这种处理方式称为逐步解调滤波算法。该算法主要包括：① 使用CPP算法估计IFCF；② 利用估计的IFCF计算谐频及转频；③ 根据各频率估计相位函数；④ 根据相位函数对信号进行解调；⑤ 对解调信号逐步带通滤波。

为了直观描述算法原理，图1给出了一个仿真信号的处理实例。假设f(t)为IFCF，fr(t)为转频，通过提取的IFCF趋势线可以估计IFCF的相位函数，假设外圈故障，便可以通过式(1)计算转频，进而估计出转频的相位函数。使用带通滤波器对解调后信号进行滤波，然后对滤波信号求和，便可以求得IFCF和转频解调滤波信号。如果假设正确，在频谱的fo和fr的位置会出现峰值。

(a) 原始信号时频谱

(b) 解调信号的时频谱

(d) 解调滤波信号的时频谱

(e) 解调滤波信号f′(t)的时频谱

(f) 逐步解调滤波信号时频谱

(g) 原始信号的频谱

(h) 逐步解调滤波信号的频谱

图1 逐步解调滤波算法流程图

Fig.1 Flowchart of stepwise demodulation filter algorithm

基于CPP和逐步解调滤波的滚动轴承故障诊断方法的具体步骤如下：

步骤1根据轴承的几何参数，结合式(1)～式(3)，计算轴承各个部分的故障特征系数Co=fo/fr，Ci=fi/fr，Cb=fb/fr。

其中，Co,Ci和Cb分别表示外圈、内圈和滚动体的故障特征系数。

步骤2对原始信号进行Hilbert变换得到包络信号，使用CPP算法在降采样包络信号中估计IFCF，并使用多项式对其进行拟合f1(t)，计算IFCF的倍频

(11)

其中，f1(t)为IFCF的拟合方程，an,…,a1,a0均为常数。

步骤3假设外圈、内圈和滚动体出现故障，并估计对应转频fro1(t)=f1(t)/Co，fri1(t)=f1(t)/Ci和frb1(t)=f1(t)/Cb以及倍频。

其中，fro1(t),fri1(t)和frb1(t)分别表示外圈、内圈和滚动体出现故障时的转频。

步骤4构造解调函滤波的相位函数

(12)

并求得解析信号y(t)=x(t)+jH[x(t)]，其中H[x(t)]是x(t)的Hilbert变换，对y(t)进行广义解调后得到解析滤波分量dk(t)=yk(t)e-2πjsk(t)。其中：fk(t)为f1(t)及其倍频；fro1(t),fri1(t)和frb1(t)及其倍频；fk(0)为对应频率的初始值。

步骤5以拟合函数fk(t)的初始值fk(0)为中心频率，fwidth为带宽，构造窄带带通滤波器。使用该滤波器对解调后的信号滤波，得到解调滤波信号mfk(t)。对滤波后的解调信号mfk(t)求和便可以得到逐步解调滤波的信号

mf(t)=∑mfk(t)

(13)

步骤6对逐步解调滤波信号mf(t)进行FFT变换得到解调滤波信号的频谱，根据频谱完成故障诊断。

3 仿真分析

为了验证算法的有效性，构造了变转速轴承仿真模型x(t)并对其进行分析。仿真模型x(t)由轴承故障引起冲击x1(t)和高斯噪声n(t)构成

x(t)=x1(t)+n(t)

(14)

故障轴承的冲击成分x1(t)的表达式

(15)

式中：N为冲击数量，Am为第m个冲击的幅值，由于轴承转速越大，冲击幅值越大，设置Am随着转速线性增加，β为结构衰减系数，u(t)为单位阶跃函数，wr为轴承故障激发的系统的共振频率，tm表示第m个冲击出现的时间，tm由递推公式(16)确定

(16)

式中：f(t)为滚动轴承转频；τ为滚动体的滑移引起的冲击间隔误差，一般为0.01～0.02；n代表轴承转动一周产生的冲击数。

设置轴承转速变化规律为f(t)=10t+5(r/s)，设置外圈和内圈故障特征因子为Co=3.5和Ci=5，其他参数如表1所示。

表1 仿真信号参数表

处理轴承外圈发生故障的仿真数据，截取3～6 s的数据进行分析。根据设置转速f(t)=10t+5(r/s)，该时间段内转速线性增加，变换范围为35～65 Hz。

图2为仿真信号的时域波形图，观察可知，随着转速升高，波形幅值变大。图3为仿真信号频谱。使用CPP算法在降采样的包络信号中估计，如图4所示的IFCF趋势线。按照逐步解调滤波算法的步骤，对信号进行处理，得到如图5所示的频谱。根据相位函数的构造方式，如果是外圈故障，解调频谱理论出现峰值的位置应该为转频的初始值35 Hz，外圈FCF的初始值35×Co=122.5 Hz及其倍数处。如果是内圈故障，则在转频，内圈FCF的初始值35×5=175 Hz及其倍数

图2 仿真信号的时域波形图

图3 仿真信号的频谱

出现峰值。观察解调频谱，出现峰值位置符合理论计算的外圈出现故障的位置，由此可判断轴承外圈出现故障。在如图6所示的包络频谱中无突出的峰值，证明了该算法的可以提高频谱质量。

图4 提取的IFCF趋势线

图5 逐步解调滤波信号的频谱

图6 原始信号的包络频谱

4 试验验证

利用振动试验台测取变转速工况下的振动信号，通过处理试验信号，对算法的有效性进一步验证。其中，为了模拟轴承故障，使用电火花对轴承的外圈和内圈加工裂纹，试验使用轴承如图8所示。试验台的结构如图7所示。轴承的参数见表2。

图7 试验台的结构

根据试验滚动轴承的具体参数得到试验轴承的外圈、内圈以及滚动体的故障特征系数分别为Co=2.55，Ci=4.45和Cb=1.7。

图8 试验滚动轴承

轴承型号滚珠数n滚珠直径d/mm节圆直径D/mm接触角α600074.817.650

对升速状态下外圈故障滚动轴承信号进行处理，图9为振动信号的时域波形图。图10为根据转速计测得脉冲拟合的转频。使用CPP算法在如图11所示的降采样的包络信号中，估计IFCF趋势线如图12所示。利用多项式拟合IFCF趋势线，并按照逐步解调滤波算法的步骤处理包络信号，得到逐步解调滤波信号。对逐步解调滤波信号进行FFT变换，得到了如图13所示的频谱图。观察可知，峰值位置符合试验轴承的外圈故障出现故障时的计算值，因此外圈存在故障。原始包络信号的频谱如图14所示。图14无突出峰值，通过对比证明该算法的有效性。