匹配压缩脊线提取在齿轮箱故障诊断中的应用*

2019-06-25王箫剑涂晓彤李鸿光李富才包文杰

振动、测试与诊断 2019年3期

王箫剑，涂晓彤，李鸿光，李富才，包文杰

(上海交通大学机械系统与振动国家重点实验室上海，200240)

引言

行星齿轮箱是由太阳轮、行星轮、大齿圈以及行星架等元件组成的机构，具有降低转速比、放大电机扭转力等功能。因其传递平稳、承载力大、体积小巧等优点，广泛应用于船舶、航天、汽车等领域中。然而此类行业对行星齿轮箱的安全平稳要求极高，所以对行星齿轮箱的工作情况进行实时监测是十分必要的。目前，越来越多的行星齿轮箱在变转速的情况下工作，比如汽车的升降速、飞机的起落过程等等。传统的频域信号处理方法仅能用于平稳信号的分析，所以一种既能定量、又能直观地描述信号故障特征的非平稳信号处理方法是极其必要的。时频分析对于提取非平稳信号的时变特性很有优势。通常使用一些传统的方法，例如连续小波变换或者短时傅里叶变换，在时频平面上呈现非平稳信号特征[1-2]。但是，由于不确定性定理，它们并不能很好地保证时频平面的分辨率。为了处理这一问题，匹配解调技术[3-5]和匹配同步压缩技术[6-7]应运而生，然而，每一时刻瞬时频率的预提取仍然是目前面临的问题。对于这些瞬时频率(时频图谱脊线)提取问题，国内外研究很多[8-10]，但在强噪声、多分量情况下，因噪声、无关分量对目标分量的干扰，对于弱分量提取较为困难。

笔者提出了匹配压缩脊线提取(demodulated synchrosqueezing ridge extraction，简称DSRE)，该方法适用于处理强噪声、多分量信号，并可应用于变转速旋转机械故障诊断与在线监测。DSRE可分为两个步骤：a.根据时频图进行脊线提取；b.根据所提脊线，进行匹配解调变换和二阶同步压缩得到新的时频图和新的脊线。为了处理多分量信号，提出了时变滤波方法，用于分量剔除与提取，得到每一阶分量的瞬时频率。将该方法应用于行星齿轮箱的故障诊断中，能够直观地提取基频及故障频率的特征，便于行星齿轮箱实时监测与故障诊断。

1 DSRE基本原理

1.1 脊线提取技术与包络线滤波技术

对于非平稳信号，由于每一时刻的瞬时频率不同，需要对其进行时频变换，例如小波变换、短时傅里叶变换等。

对于任意一个信号x(t)=A(t)ei2π(f0t+φ(t))，f0为其载波频率，φ(t) 为其瞬时相位角，对应φ′(t) 为其瞬时频率。其短时傅里叶变换可以表示为

(1)

通常选用标准差为σ的高斯窗函数

(2)

由于短时傅里叶变换的时域窗长度保持恒定，因而无法有效处理时变信号。

在时频图上每一时刻具有最大能量的点的序列，被定义为脊线。脊线的纵坐标，对应着信号的瞬时频率。

对于脊线提取的方式有很多种[8-10]，笔者用于脊线提取的过程[8]如下。

1) 对信号经过时频变换得到的图谱，计算每一时间点幅值旳极值Qm(t)及其所对应的瞬时频率νm(t)

(3)

其中：TFD(t,ω)为t时间点ω频率处对应的时频变换幅值；Np(t)对应t时刻的最大极值点个数；m为t时刻对应的第m个最大极值点。

2) 定义路径优化函数

{mc(τ1),mc(τ2),…,mc(τN)}=

{νm1(τ1),…,νmN(τN)]

(4)

F由当前时刻的时频谱幅值Qm(t)(当前时刻的能量)，及当前时刻的频率与上一时刻的频率之差νm(tn)-νk(tn-1) (时频图谱脊线的连续性)共同决定

F[Qm(tn),νm(tn),νk(tn-1)]=

logQm(tn)+w1(νm(tn)-νk(tn-1))

(m∈{1,2,…,NP(tn)}；k∈{1,2,…，NP(tn-1)})

(5)

其中：w为衡量脊线连续性的因子。

w1(Δξ)=-σ2fs|Δξ|

(6)

3) 对上述路径优化函数进行求解，得到的设计变量即为每一时刻点对应的瞬时频率。

此路径优化函数的优点是既考虑了脊线当前点的能量，又考虑了脊线的连续性。在处理含噪声的信号时，从脊线连续性角度，可以有效地消除噪点对所提信号信息的干扰。

对于多分量信号，可经过包络线滤波技术去除所提脊线分量信号。在时频图上，找到对应所提脊线的上下包络线，将上下包络线中间部分幅值均重置为零。

1.2 旋转算子滤波-匹配解调技术

多项式匹配解调算法的思路是在处理信号每个时间点邻域的时候，消除信号脊线在时间点邻域的高次项，从而在每一个时间点的邻域，频率都不随时间改变，即在所有高斯窗里，信号频率均时不变。

将上文所提取的脊线，作为瞬时频率参考，据此提出了旋转算子dR(t)和平移算子dS(t,u)

其中：φ′(t)为信号的瞬时频率，由脊线经多项式拟合得到；u为t时刻邻域对应时间点。

旋转算子的作用是去除整个时域内瞬时频率脊线的1阶及其以上导数，而平移算子的作用是将每个时间点邻域赋予恒等同于该时间点的瞬时频率，如图1所示。

图1 匹配解调过程旋转变换及平移变换Fig.1 Rotating step and shifting step in the moduling process

经过旋转平移变换之后，信号Sd(t,u)为

SdR(t)=x(t)dR(t)=A(t)ei2πf0t

(9)

Sd(t,u)=sdR(t)dS(t)=A(t)ei2π[f0t+φ′(u)t]

(10)

此信号经过短时傅里叶变换之后，即得到较为精确的时频图。

匹配解调方法的表达式为

(11)

对于多分量信号，可通过旋转算子滤波得到目标分量。该分量脊线提取旋转算子后，对信号旋转变换，信号将近似于时不变信号，此时滤波，即可在时域信号上去除其他分量和噪声，如图2所示。

图2 匹配解调-旋转算子滤波Fig.2 The filtering step in the moduling process

对滤波之后的信号进行平移变换和短时傅里叶变换，即可得到所需信号的时频图。

1.3 2阶匹配同步压缩技术

匹配同步压缩技术就是把信号的各时频点在信号的瞬时频率处进行重排，从而使能量汇集在信号的瞬时频率附近。本研究匹配同步压缩的目的是使匹配解调所得的时频变换图谱上面的能量更集中，从而有利于接下来脊线的提取。

时频信号的压缩方法如下。

以一个信号的短时傅里叶变换为例

(12)

首先，将此短时傅里叶变换经过从u-t到v的换元可得

(13)

对于一个恒定幅值的信号，将信号瞬时相位泰勒级数展开，并保留二次项系数

x(v+t)=Ae(a+b(t+v)+0.5c(t+v)2)i

(14)

其中：a，b，c分别为泰勒展开常数项、一次项、二次项系数。

再假设信号频率1阶导数远大于2阶导数的情况下，将式(14)经过对时间的求导，可得

R(∂tSTFT(t,w))=

(15)

傅里叶变换的导数可以表示为

R(∂tSTFT(t,w))=i(b+ct)STFT(t,ω)+

icTSTFT(t,ω)

对比上面两式发现

(16)

然后，将短时傅里叶图谱对于式(16)解得的瞬时频率再赋值，将每一时间点全部能量聚集到瞬时频率周围，可以得到压缩之后的图谱

(17)

其中：Ssd(t,f)为1.2节变换得到的时频分布；g*(0)为信号能量加权因子；γ为能量阈值。

1.4 多分量信号提取故障分量方法实施

对于多分量信号，首先，提取能量最大信号分量的脊线，时频变换得到所需信息后，去除该分量，并提取下一阶信号分量；其次，依次提取余下信号分量。具体步骤如下：

1) 进行短时傅里叶变换，用脊线提取法估计出能量最大分量的脊线，同时采用包络法估计脊线的上下边界；

2) 用多项式拟合所得脊线，得到多项式各参数，然后采用匹配解调技术，经过旋转算子变换后滤波，并进行平移变换得到新的时频图；

3) 根据得到的瞬时频率进行2阶同步压缩，得到该分量的时频图，并进行第2次脊线提取提出修正之后的脊线及其上下边界；

4) 在时频图上根据脊线及其上下边界去除能量最大的分量，重复步骤1～3，直到提取出故障信号。

对每一阶信号分量的提取，该法的计算代价为O(3NfN)+O(M2N)，其中：N为信号时域采样点；Nf为频域点数目；M为1.1节中Np(t)最大值。

2 仿真分析

2.1 仿真信号构造

为了检验该信号处理方法，采用含有高斯白噪声的非平稳仿真信号，仿真信号f(t)如下

(18)

该仿真信号的采样频率为1 024 Hz，由3种瞬时频率随时间改变的信号分量叠加而成，同时加入高斯白噪声η(t)，该白噪声平均值为0，信噪比为-1 dB。时域波形如图3(a)所示，频谱图如图3(b)所示。

图3 仿真信号特征Fig.3 Characteristics of the simulation signal

2.2 能量最大信号分量提取

对信号进行时频变换并提取脊线，所得结果如图4(a)所示。图4(a)为采用短时傅里叶变换得到的时频图，红线代表采用文献[8]提取的脊线。如图所示，所提脊线虽然连续，并能反应信号特征，但由于噪声的存在，脊线存在一定波动，和真实情况相比具有误差。

对信号进行匹配解调旋转变换，得到时频图见图4(b)。如图所示，信号经过匹配解调的旋转变换后，其第1阶分量的瞬时频率基本保持不变。

对所得信号进行傅里叶变换及带通滤波，滤去噪声信号和其他阶分量信号，如图5所示。

为了验证所提脊线的准确性，比较本研究方法所提脊线的精确程度与文献[8]中所述方法的精确程度，如图6所示。其中：点划线为本研究方法所提脊线；虚线为文献[8]所提脊线；实线为真实脊线。由图可以看出，本研究算法所提脊线与真实脊线基本重合，而文献[8]中所提脊线在真实脊线附近存在波动。此波动误差会随着所提分量的阶数的增高而逐渐积累，从而影响弱能量分量提取的精度。

图4 旋转变换前后时频图谱Fig.4 Time-frequency transform before and after the rotating process

图5 旋转变换信号滤波前后傅里叶变换图谱Fig.5 The Fourier transform of the rotated signal before and after the filtering process

图6 第1阶分量本研究所提脊线与真实脊线对比Fig.6 Comparation between the extracted 1 st ridge curves and true 1 st ridge curve

笔者采用范数来定量表现所提脊线精确程度

(19)

其中：EIF为估计瞬时频率；TIF为真实瞬时频率。

经计算可得本研究算法所得误差为0.001 0，而文献[8]所述算法误差为0.003 0，故本算法所提脊线精度比文献[8]脊线准确。

2.3 能量次大信号分量提取

在提取出第1阶信号瞬时频率之后，在时频图上将第1阶信号频率成分滤去，并进行第2次脊线提取，如图7所示。

图7 去除1阶信号分量后时频变换图谱及脊线预提取Fig.7 The time-frequency transform and previous ridge extraction after removing the first signal

对于该多分量信号，图8为第2阶分量的脊线对比图。其中：点划线为本研究所提脊线；虚线为采用文献[8]并结合本研究的脊线流程所提脊线；实线为真实脊线。由图可以看出，本研究所提脊线拟合程度更高。

对于该脊线，本算法范数为0.019 1，文献[8]算法范数为0.048 4，故本算法精确程度更高。

图 8 第2阶信号分量所提脊线与真实脊线对比Fig.8 Comparation between the extracted 2st ridge curves and true 2st ridge curve

3 行星齿轮箱故障诊断

本试验的振动数据采集自行星齿轮箱故障模拟试验台。如图9所示，试验台包括驱动电机、交流电机、行星齿轮箱、固定轴齿轮箱和制动器。采用转速传感器和加速度传感器分别采集电机轴转速和行星齿轮箱振动信号，电机轴转速可由驱动电机自由调节。本试验中电机为变转速，包括一段升速(1 800～2 700 r/min)和一段减速过程(2 700～1 800 r/min)。信号采样频率为12 800 Hz，采集时长为16 s。

图9 行星齿轮试验台Fig.9 The test rig of the planetary gearbox

表1为行星齿轮箱的齿轮参数，表2为行星齿轮箱内各零件的特征阶次。

表1 行星齿轮箱齿轮参数

表2 行星齿轮箱故障阶次

Tab.2 Fault orders of the planetary gearbox

10.1672.500.83

对振动信号原始数据时频变换，信号的时域波形如图10(a)所示，时频变换图谱如图10(b)所示。由于外界噪声和信号多重分量的叠加，难以从时域信号或时频变换图像中发现故障特征。信号能量最集中的分量是齿轮加速度信号的三倍频分量，即使采用传统方法提取能量最集中的信号分量，再进行阶次分析的方式也难以实现。

图10 行星齿轮箱振动信号特征Fig.10 Characteristics of the acceleration signal of vibration of a planetary gearbox

首先，采用本方法提取基频。在这个信号中，基频并不是能量最大的分量，能量最大的分量是三倍频，所以先提取信号的三倍频，如图11(a)所示。提取三倍频并滤去三倍频信号分量之后，即可提取信号一倍频特征频率，如图11(b)所示。

以此类推，将倍频成分全部去除后，继续进行脊线提取，这样即可提取出故障频率，如图12所示。实线为采用文献[8]方法并结合本研究提脊线流程所提脊线；虚线为本研究提取脊线。由图可知，本研究所提脊线方法可行，且较为准确。

将所提故障分量信号与基频信号分量作商，即可得到故障信号分量，如图13(a)所示。由图可知，故障频率为2.5阶，为太阳轮故障。将所得故障信号分量与基频信号分量于一张时频图上呈现，亦可直观地表现故障信号的变化过程及其强弱，如图13(b)所示，其中：上方线为故障分量；下方线为基频分量。若时间足够长，可通过时频图的明暗观察故障的变化情况。