基于Vold-Kalman广义解调的变转速轴承和齿轮复合故障诊断

2019-04-03 01:17赵德尊李建勇程卫东温伟刚

振动与冲击 2019年6期

赵德尊, 李建勇, 程卫东, 温伟刚

(1.清华大学机械工程系，北京 100084； 2.北京交通大学机械与电子控制工程学院,北京 100044; 3.载运工具先进制造与测控技术教育部重点实验室,北京 100044)

轴承和齿轮是机械结构中用于支撑和传递力矩的重要部件[1]。由于复杂的工作环境和运行状态，轴承和齿轮也是最易受损的零部件。实际工况中，轴承和齿轮往往同时发生故障。受转速波动的影响，多个故障特征相互影响、彼此干扰的同时其自身重复频率也表现出时变非平稳性，给故障诊断带来困难。因此，探究变转速条件下轴承和齿轮复合故障特征的提取具有重要意义。

对于故障轴承，随着故障点与其表面间的接触，测量的振动信号中将产生一系列短时、快速衰减的冲击脉冲。冲击脉冲的重复频率即故障特征频率是确定轴承故障点位置的重要指标，另外，由轴承故障引起的高频共振将调制故障特征频率[2]。对于故障齿轮，有文献指出，其振动信号中除了含有具有幅值优势的啮合频率外还包含由齿轮故障激起的共振带，该共振带同样会调制齿轮转频[3]。转频可作为有效指标诊断齿轮故障[4]。据此，计算最优的滤波参数，区分提取轴承和齿轮分别激起的共振带，是其复合故障特征提取的关键步骤，也引起了诸多学者的关注。

基于谱峭度对振动信号中冲击脉冲的敏感性，Antoni[5]提出了快速 Kurtogram 算法，该算法可以有效确定最优带通滤波参数，通过带通滤波提取由机械部件局部故障引起的高频共振。为进一步提高快速Kurtogram算法抗干扰能力，Wang等[6]提出了基于二进制小波包分解的快速 Kurtogram增强算法，并成功的应用于滚动轴承故障诊断。Combet等[7]提出了利用谱峭度值确定最优滤波参数的齿轮早期故障特征提取方法。基于谱峭度的最优滤波参数计算方法得到诸多学者认可和广泛应用的同时，其对于复合共振带提取的无能无力也引起了学者的重视。Wang等[8]改进了谱峭度滤波算法，即通过设定啮合指标确定只含有齿轮共振带的啮合共振谱峭度图，继而与包含轴承和齿轮共振带信息的Kurtogram图对比分析完成轴承和齿轮复合共振带的区分。然而，该方法只适用于恒转速条件，对于变转速复合故障特征的提取无能为力。

变转速模式普遍存在于旋转机械的实际运行中[9]。变工况下，表征轴承故障的故障特征频率和表征齿轮故障的转频将表现出时变非平稳的特点。现有的以FFT(Fast Fourier Transform)为核心的故障诊断方法将因频谱模糊现象而失去作用。针对传统的计算阶比分析方法在计算效率和精度方面的不足[10]，探索无需角域重采样的复合故障特征提取方法对提高故障诊断技术的有效性和实用性具有重要意义。

Vold-Kalman阶比跟踪是Vold等[11]基于Kalman滤波器提出的可以用于估计振动信号中某一谐波成分的算法。与传统的带通滤波或经验模式分解算法相比，该算法可在时域中直接提取某一特定的近似、交叉的谐波分量，避免了由时域至频域变换带来的相位偏差和能量泄露[11-12]，因此受到了诸多学者关注。其基本原理从第一代发展到第二代，具备了成熟的理论体系，并广泛的应用于商业软件中。然而现有的Vold-Kalman滤波算法均是与阶比分析相结合处理变转速旋转机械振动信号[13-14]。广义解调算法是Olhede等[15]提出的一种可以将时频分布是倾斜、非线性的非平稳信号转换成时频分布是线性的且平行于时间轴的平稳信号的分析方法。基于广义解调算法对于调幅-调频信号的处理优势，很多学者将其应用于旋转机械设备的故障诊断[16]。

综上所述，针对传统的以谱峭度为核心的带通滤波和计算阶比分析算法的不足，结合Vold-Kalman滤波和广义解调变换的优势，提出了滚动轴承和齿轮时变非平稳复合故障特征提取方法。该方法舍弃了传统的高频共振带提取和角域重采样的故障诊断策略，直接利用Vold-Kalman滤波提取时变非平稳的故障特征成分，其次通过广义解调实现提取的时变非平稳故障特征成分的平稳化重置，最后基于频谱分析完成重置成分的定量表达。仿真和实验分析结果证实了本文算法提取、重置和定量表达故障特征成分的有效性；同时避免了传统带通滤波算法带来的相位偏差和能量泄露，频谱中故障特征成分具有更明显的幅值优势。

1 理论基础

1.1 Vold-Kalman滤波

旋转机械中产生的调制信号可表示为幅值和载波信号的乘积

(1)

式中：Ak为第k个谐波成分的幅值；φk(t)为载波成分，其表达式如式(2)所示

(2)

由于机械系统的固有惯量，调制信号中的幅值包络有别于载波信号以较低的频率缓慢变化，因此可以用低阶多项式表示。对于离散信号其状态方程可表示为

(3)

设定s=2，式(3)可改写为

Ak(n-1)-2Ak(n)+Ak(n+1)=εk(n)

(4)

将式(4)展开为矩阵

(5)

因此，不同阶次的多项式具有相同的矩阵形式

MA=ε

(6)

式中：M为N×N矩阵。

实测信号包含各谐波成分和噪声或误差成分，其观测方程可表示为

(7)

式中：y(n)为实测振动信号；δ(n)为噪声或误差项。

为提取特定的谐波分量xk(t)，系统的观测方程可表示为

(8)

其矩阵形式为

Y-CA=δ

(9)

根据系统状态方程式(6)和观测式(9)以及预先设定的目标分量的载波矩阵C，应用最小二乘法，使得非齐次项ε(n)和噪声或误差δ(n)的平方和最小[17]，即

(10)

式中：C*为C的共轭转置；A*为A的共轭转置；r为加权因子，具体取值可参考文献[18]。

进一步求得目标分量的幅值包络矩阵A

A=(r2MTM+E)-1C*Y

(11)

通过载波矩阵C和幅值包络A即可重构目标分量

X=AC

(12)

1.2 广义解调变换

广义解调变换的本质为广义傅里叶变换。对于任意的解析信号x(t)，其广义傅里叶变换的定义为

(13)

式中：s0(t)为随时间t变化的实值函数，即相位函数；广义傅里叶变换实际上是对x(t)e-2πjs0(t)做标准的傅里叶变换。同样对XG(f)进行逆傅里叶变换可得

(14)

根据式(14)，如果XG(f)=δ(f-f0)，x(t)可改写为

x(t)=e-2πj[f0t+s0(t)]

(15)

通过式(15)可知，解析信号中某一特定时变非平稳成分的相位函数s0(t)只要满足式(16)，其即可被转换成线性的、平行时间轴的平稳成分。

(16)

式中：f(t)为特定时变非平稳成分的频率方程；f0为频率方程对应的频率点；s0(t)为频率方程对应的相位函数；ds0(t)/dt为对s0(t)求导。

基于上述广义解调变换的理论介绍，给出对原始信号x(t)中感兴趣的某一时变非平稳成分的解调过程：①对原始信号x(t)进行希尔伯特变换，获得解析信号y(t)=x(t)+jH[x(t)]，其中H[x(t)]是x(t)的希尔伯特变换结果；②根据感兴趣的时变非平稳成分的频率方程和式(16)计算其对应的相位函数s0(t)，基于相位函数得到变换结果d(t)=y(t)e-2πjs0(t)；③对变换结果d(t)进行希尔伯特变换得到新的解析信号z(t)=d(t)+jH[d(t)]，即原始信号x(t)的广义解调信号。

2 本文算法

故障特征频率和转频分别是判定轴承和齿轮健康状况的关键性指标。本文利用Vold-Kalman滤波算法和广义解调变换对时变非平稳的故障特征成分进行提取和平稳化重置，进而通过频率谱对重置的故障特征成分进行定量表达，最终基于理论频率点和频谱中的峰值判断轴承和齿轮的健康状况。计算故障特征成分方程及其对应的相位函数和频率点是所提算法的关键。本文利用转频方程、机械系统中各转频间的比例参数以及目标轴承故障特征系数对上述参量进行计算。本文方法的流程图如图1所示。具体的步骤总结如下：

图1 算法流程图Fig.1 Flowchart of proposed diagnosis strategy

步骤1同步测取振动信号和转速脉冲信号，并根据转速信号计算转频方程。

步骤2基于转频方程、机械结构传动比、目标轴承故障特征系数，计算轴承和齿轮故障特征频率方程，并根据式(16)计算频率方程对应的相位函数和频率点。

步骤3通过故障特征频率方程计算该故障特征对应的载波矩阵C，并根据载波矩阵分别对原始信号的包络进行Vold-Kalman滤波，即利用式(11)求取载波矩阵C对应的幅值包络矩阵A，通过式(12)实现时变非平稳故障特征成分的重构。

步骤4基于各个滤波分量对应的相位函数，分别对滤波分量进行广义解调，得到解调信号。

步骤5对各个解调成分进行频谱分析。通过频谱峰值和理论频率点的对比，初步确定故障点位置。

步骤6基于初步诊断结果，提取其二倍和三倍谐波成分。通过频率谱的定量表达进一步验证步骤5诊断结果的正确性。

3 仿真验证

为验证本文算法的有效性，构造了变转速条件下滚动轴承和齿轮复合故障振动信号。

变转速故障轴承信号的仿真模型为

(17)

式中：Ai,b为第i个故障冲击幅值；βb为结构衰减系数；ωb为由轴承引起的共振带中心频率；τj为由滚动体滑移带来的故障冲击间隔误差，其取值一般为0.01～0.02；μ(t)为单位阶跃函数；ti,b为第i个冲击发生时刻，可由式(18)计算

(18)

式中：fcc为轴承故障特征系数；fb(t)为轴承转频。

变转速下故障齿轮仿真模型如下

(19)

式中：第一行为齿轮啮合振动；第二行为由齿轮故障引起的高频共振；n为齿轮啮合频率序号；本文取值1，2和3；An,g为第n个啮合频率幅值；z为齿轮齿数；Fg(t)为齿轮转频的积分函数；Am,g,βg,ωg和tm,g的含义等同于式(17)中的Ai,b,βb,ωb和ti,b。其中tm,g的计算公式为

(20)

式中：fg(t)为齿轮转频。

因此，变转速条件下轴承和齿轮复合故障振动信号的仿真模型可表示为

x(t)=xb(t)+xg(t)+n(t)

(21)

式中：n(t)为高斯白噪声。仿真模型的参数见表1。

表1 变转速下轴承和齿轮复合故障仿真模型参数

根据仿真模型得到的信号的时域表达，如图2所示。

图2 仿真信号Fig.2 Simulated signal

根据式(16)计算仿真模型中预设的转频方程fb(t)对应的相位函数sb(t)=1.875t2和频率点fpb=39 Hz。基于转频方程fb(t)及其对应的相位函数sp(t)和频率点fpb、传动比和轴承故障特征系数分别计算轴承和齿轮故障特征频率方程及其对应的相位函数和频率点，详见表2。此外，为说明本文算法不会因不存在的特征成分，在分析结果中出现相应峰值以造成误判，随机选取了故障特征系数fcc1=2.8作为对比分析。

对原始信号进行希尔伯特变换得到包络信号。首先，分别利用表2中的故障特征频率方程对包络信号进行Vold-Kalman滤波，得到四个滤波信号。其次，根据表中的相位函数分别对相应的滤波信号进行广义解调，得到四个解调信号。再次，对解调信号分别做FFT得到频率谱依次如图3中四个小图所示。最后，通过频率谱中峰值与表2中频率点的对比确定故障点位置。图3中只有两张图中存在明显峰，即只有分别利用表2中第一行和第三行参数获得的频率谱中存在明显峰值，且两图中的峰值144.3和25.98分别与表2中频率点fp0和fpd1近似或相同。因此可以初步判定该仿真模型中轴承存在与故障特征系数3.7关的故障，以及主动轮存在故障。而分别利用表2中第二行和第四行参数获得的频率谱中在相应的频率点fpc和fpd2处并没有出现峰。说明上述信号中不存在与故障特征系数fcc1=2.8有关的故障特征，而被动轮也不存在故障。

表2 仿真复合故障特征频率方程及其对应的相位函数和频率点

图3 Vold-Kalman广义解调成分频率谱Fig.3 Spectra of Vold-Kalman generalized demodulation components

为进一步验证上述诊断结果，分别对图3(a)和图3(c)中峰的二倍及三倍频进行提取、重置以及定量表达。其中Vold-Kalman滤波用到的频率方程分别为2f0(t)，3f0(t)，2fd1(t)和3fd1(t)。广义解调用到的相位函数分别为2s0(t)，3s0(t)，3sd1(t)和3sd1(t)。结合初步验证结果得到完整频率谱分别如图4和图5所示。图4中第二、第三个峰值分别为第一个峰值的2倍、3倍。同样，图5中的峰值也存在上述比例关系。因此得出结论：仿真信号中轴承存在与故障特征系数3.7相关的故障以及主动齿轮存在故障。诊断结果与仿真模型相一致。

4 实验验证

本节在机械故障仿真试验台上测取了减速条件下轴承和齿轮复合故障特征振动信号，进一步验证所提算法的有效性。试验台如图6所示，其中目标轴承存在内圈故障，其故障为人为破坏形成的裂纹故障，故障程度为轻度(宽：0.2 mm 深：0.4 mm)。齿轮箱中被动轮存在断齿故障，即其中一个齿被人为破坏，出现断裂。试验台机械结构参数见表3。本试验中采样频率设置为48 000 Hz，截取2 s时长的信号(见图7)进行分析。上述信号对应的转频曲线如图8所示。该转频曲线由转速计测取的轴的转速脉冲信号计算。对转频曲线进行拟合得到拟合方程：fa(t)=-0.042t3+0.029t2-2.767t+30.86。基于式(16)，对上述方程进行改写，计算其对应的相位函数sa(t)=-0.011t4+0.01t3-1.38t2和频率点fpa=30.86 Hz。

图4 轴承频谱Fig.4 Bearing spectrum

图5 齿轮频谱Fig.5 Gear spectrum

根据转频地fa(t)及其对应的相位函数sa(t)和频率点fpa，以及表3中的参数计算所有故障特征成分方程及其相位函数和频率点如表4所示。

图6 机械故障仿真试验台Fig.6 Machinery fault simulator

参数值轴承类型ER-10k滚珠直径d/mm7.94节圆直径D/mm33.5滚珠数量nb8接触角α0fcco3.052fcci4.95fccb1.99主动轮齿数zg18带传动比ib2.56齿轮传动比ig1.5

首先根据该工况下外圈故障时对应的故障特征频率方程和相位函数对原始信号的包络进行处理。即，基于fo(t)对原始信号的包络进行Vold-Kalman滤波获取表征轴承外圈故障的时变非平稳成分，基于相位函数so(t)对滤波成分进行广义解调，最后利用FFT计算解调信号的频率谱。同样分别利用表4中内圈、滚珠、主动轮和被动轮故障时对应参数获得相应的频率谱。所有频谱中，只有利用轴承内圈和主动轮故障时对应的参数获得的频谱中存在明显峰，且峰值153.1和8.057分别与频率点fpi和fpdrive近似，如图9所示。因此可以初步判定该轴承内圈存在故障，齿轮箱中被动轮存在故障。

图7 实测振动信号Fig.7 Measured signal

图8 实测振动信号对应的转频曲线Fig.8 Rotational frequency trend

分别利用方程2fi(t)，3fi(t)，2fdriven(t)和相位函数2si(t)，3si(t)，2sdriven(t)对原始信号的包络进行Vold-Kalman广义解调。整合上一步骤中获得的频率谱，得到轴承和被动轮完整频率谱分别如图10和11所示。图中峰值近似等于频率点fpi，2fpi，3fpi，fpdriven和2fpdriven，进一步验证了诊断结果的正确性。

为验证本文算法在避免能量泄露和消除无关项干扰方面的优越性，利用带通滤波与广义解调相结合的算法对本节所测量的原始振动信号进行分析。处理过程如下：首先基于原始信号的时频表达确定带通滤波参数提取轴承故障激起的高频共振带；其次利用相位函数Si(t)对原始信号的包络进行广义解调；最后通过FFT获得频率谱如图12所示。本文方法获得相对应的频谱图如图9(a)所示。对比分析图9(a)和图12可知本文算法获得频谱图更加简洁、无其他干扰成分。另外，本文算法获得的频谱峰值(7.6×106)相比图12中的(1.4×105)更有优势。