基于多尺度线调频基稀疏信号分解的齿轮故障诊断方法

摘要

本发明公开了一种基于多尺度线调频基稀疏信号分解的齿轮故障诊断方法。本发明方法将齿轮振动信号的时间跨度在各尺度下等分为动态时间支撑区，然后在动态时间支撑区上建立线性调频基函数，形成多尺度线性调频基函数库，计算每个动态时间支撑区上投影系数最大的线性调频基函数，通过对最大系数线性调频基函数的连接，自适应的形成与齿轮振动信号具有最大相关系数且瞬时频率具有物理意义的分解信号，进而可以获得非平稳转速下故障齿轮的幅值调制频率和啮合频率随时间的变化情况，判断故障位置，进行故障诊断。

说明书

技术领域

本发明涉及一种齿轮故障诊断方法，特别涉及一种基于多尺度线调频基稀疏信号分解的齿轮故障诊断方法。

背景技术

当齿轮箱发生故障且转速恒定时，通常会形成以齿轮啮合频率及其高次谐波为载波频率，以齿轮所在轴转频及其高次谐波为调制频率的调制现象。但当齿轮箱转速波动时，由于转频的波动会同时导致调制频率和载波频率的波动，因此其振动信号为非平稳信号，对齿轮箱振动信号进行FFT变换，将难以对调制边频带进行识别，进而难以依据调制边频带诊断齿轮故障。常用的齿轮振动信号时频分析方法有小波变换和EMD方法。小波变换由于采用可变的时频窗函数，素有“数学显微镜”之称，但是由于时频不确定原理，无法同时在时域与频域取得较高分辨率。此外，小波变换的分解尺度只与信号的采样率有关，而与信号本身无关，所以本质上小波变换不是一种自适应的信号分解方法。EMD方法将多分量信号自适应的分解为若干个瞬时频率具有物理意义的IMF(Intrinsic Mode Function)分量之和，进一步采用Hilbert变换求出每个IMF分量的瞬时频率和瞬时幅值，从而实现对复杂信号的解调。但是实际上每个IMF分量并非单分量信号，仍为多分量信号，这导致其瞬时频率出现无法解释的不规则性，所以EMD方法不适合于窄带多分量信号的分解。同时，EMD方法在理论上还存在诸于过包络、欠包络、模态混淆和端点效应等问题，需要进一步研究解决。

发明内容

为了解决现有齿轮故障诊断方法存在的上述技术问题，本发明提供一种基于多尺度线调频基稀疏信号分解的齿轮故障诊断方法。

本发明解决上述技术问题的技术方案包括以下步骤：

1)利用加速度振动探头对齿轮箱进行测量，获得加速度振动信号，信号长度为2的整数次方；

2)将振动加速度信号长度以N/2^j为长度进行等分，形成动态时间支撑区I＝[kN2^-j～(k+1)N2^-j]，I为动态动态时间支撑区，j为分析尺度系数，j＝0，1，...，logN-1，N为振动加速度信号长度，k＝0，1，...，2^j-1；

3)在动态时间支撑区上定义多尺度线性调频基函数库

4)将加速度振动信号对多尺度线性调频基函数库中基函数进行投影，计算每个动态时间支撑区上的最大投影系数及其对应的基函数：

5)根据最大投影系数及其对应的线性调频基函数得到该动态时间支撑区内的分解信号c_I(t)，

$c_I\left(t\right)=\mathrm{abs}\left(2{\beta }_I\right)e^{-i(a_{\beta }t+b_{\beta }{t}^2-\mathrm{angle}\left(2{\beta }_I\right))}{1}_I\left(t\right)$

β_I为最大投影系数，1I(t)为矩形窗函数，当t∈I时为1，当$t\notin I$时为0，I为最大投影系数所对应线调频基函数的动态时间支撑区，a_β为对应线调频基函数的频率偏置系数，b_β为对应线调频基函数的频率斜率；

6)连接动态时间支撑区下的分解信号c_I(t)，形成覆盖加速度振动信号长度的信号分量，保留能量最大的信号分量作为本次分解的信号分量；

7)从加速度信号中减去分解信号分量，形成残余信号；

8)将残余信号能量与加速度振动信号能量之比与终止阈值比较，若大于终止阈值，则将残余信号作为新的分解信号重复4～8步；如果小于终止阈值则停止分解，转到步骤9；

9)根据分解信号分量判断齿轮故障及其位置。

本发明的技术效果在于：本发明方法将齿轮振动信号的时间跨度在各尺度下等分为动态时间支撑区，然后在动态时间支撑区上建立线性调频基函数，形成多尺度线性调频基函数库，计算每个动态时间支撑区上投影系数最大的线性调频基函数，通过对最大系数线性调频基函数的连接，自适应的形成与齿轮振动信号具有最大相关系数且瞬时频率具有物理意义的分解信号，获得非平稳转速下故障齿轮的幅值调制频率和啮合频率随时间的变化情况，可对非平稳转速下齿轮故障进行较准确的诊断。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的基于多尺度线调频基的稀疏信号分解流程图。

图2为本发明中连接动态时间支撑区下的分解信号流程图。

图3为本发明中断齿齿轮振动信号时域波形图。

图4为本发明中断齿齿轮振动信号频谱图。

图5为本发明中断齿齿轮分解信号时频图。

图6为本发明中2倍转频曲线与分解所得调制频率时频图。

图7为本发明中正常齿轮分解信号时频图。

图8为本发明中1、2倍转频曲线与分解所得调制频率时频图。

具体实施方式

图1为本发明的基于多尺度线调频基的稀疏信号分解流程图。下面结合流程图对基于多尺度线调频基的稀疏信号分解方法原理进行详细说明。

1)利用加速度振动探头对齿轮箱进行测量，获得振动加速度信号，采样长度定为2的整数次方，根据轴承转速和齿轮齿数设定采样频率；

2)划分动态时间支撑区；将振动加速度信号长度N以2的j次方进行划分，形成动态时间支撑区，I＝[kN2^-j～(k+1)N2^-j]，j为分析尺度系数，j＝0，1，...，logN-1，N为采样长度，k＝0，1，...，2^j-1；

3)在动态时间支撑区上定义多尺度线性调频基函数库

$D\left(h_{a_{\mu },b_{\mu },I}\right)=\{h_{a_{\mu },b_{\mu },I}\left(t\right)=K_{a_{\mu }.b_{\mu },I}{e}^{-i(a_{\mu }t+b_{\mu }{t}^2)}{1}_I\left(t\right)\}$

要求$\left|\right|h_{a_{\mu },b_{\mu },I}\left|\right|=1$则

$\left|\right|h_{a_{\mu },b_{\mu },I}\left(t\right)\left|\right|=K_{a_{\mu },b_{\mu },I}\times N$

$K_{a_{\mu },b_{\mu },I}=1/N$

上式中：D为基函数库；为多尺度线性调频基函数；I为动态时间支撑区；为归一化系数，使得$\left|\right|h{a}_{\mu },b_{\mu },I\left|\right|=1;$a_μ为频率偏置系数，b_μ为频率斜率。根据采样定理a_μ+2b_μ应该小于f_s/2，f_s为采样率；1_I(t)为矩形窗函数，当t∈I时为1，当$t\notin I$时为0。

在尺度系数j下，当振动加速度信号长度为N时，每个动态时间支撑区包括N/2^j个采样点，采样长度被分割为2^j个动态时间支撑区。尺度系数j取0～logN-1，因此本发明方法要求采样长度为2的整数次方。定义的多尺度线性调频基函数在动态分析时间段内的瞬时频率为a_μ+2b_μ。

4)将加速度振动信号对多尺度线性调频基函数库中基函数进行投影，求出每个动态时间支撑区上的最大投影系数及其对应的基函数：

根据信号分析理论，任意信号f(t)可以展开为一组基函数的线性组合，即

$f\left(t\right)=\underset{n\in Z}{\Sigma }{a}_n{h}_n$

如果该组基函数为正交基，则可用内积计算它们的展开系数，即：

a_μ＝<f(t)，h_n>/‖h_n‖

a_n的大小反映了f(t)与基函数的相似程度。

支撑区I内的最大投影系数(即展开系数)β_I计算公式为：

${\beta }_I=\underset{I}{\max }<f\left(t\right),h_{a_{\mu },b_{\mu },I}\left(t\right)>$

最大投影系数β_I中包含了分解信号的幅值和初始相位信息。推导过程如下：

若分析信号为：

θ(t)为分析信号相位函数，为初始相位，则有：

$<f\left(t\right),h_{a_{\mu },b_{\mu },I}\left(t\right)>$

周期函数在一个周期内的积分为0，即：

$\underset{T_0}{\int }{e}^{-i2\pi /T_0\times t}\mathrm{dt}=0$

T₀为周期信号周期，所以：

当θ(t)-a_μt-b_μt²＝0时不等式成立，有：

所以最大投影系数β_I中包含了分解信号的幅值(r)和初始相位信息，对应的基函数中则包含了分解信号的频偏信息(a_μ)和频率斜率信息(b_μ)；

5)根据最大投影系数和对应的线性调频基函数定义该动态时间支撑区下的分解信号c_I(t)：

$c_I\left(t\right)=\mathrm{abs}\left(2{\beta }_I\right)e^{i(a_{\beta }t+b_{\beta }{t}^2-\mathrm{angle}\left({2\beta }_I\right))}{1}_I\left(t\right)$

c_I(t)为动态时间支撑区I内最大投影系数对应的分解信号；在动态分析时间段I内

f_I(t)＝c_I(t)+r_I(t)

r_I(t)为分解的残余信号，因为r_I(t)与c_I(t)对应的基函数正交，所以‖f_I(t)‖²＝‖c_I(t)‖²+‖r_I(t)‖²

投影系数的最大化保证了残余信号能量的最小。通过公式可以计算出每个动态分析时间段内的最大投影系数和相应的基函数。

6)连接动态时间支撑区下的分解信号c_I(t)，形成覆盖加速度振动信号长度的信号分量，不同的连接方法形成不同的信号分量，保留能量最大的信号分量为本次分解的分解信号分量；

为了使整个分析时间段内的残余信号能量最小，需要采用合适的动态时间支撑区连接方法，在该连接方法下满足在整个分析时间段内使分解信号的总能量最大，即：

$\mathrm{Max}\left(\underset{I{\in \Pi }^n}{\Sigma }\right|\left|c_I\left(t\right)\right|{|}^2)$${\Pi }^n=\{I_1^n,I_2^n,...\}\in \left\{I\right\}$

式中n代表第n次分解，且∏ⁿ覆盖整个振动信号时间长度，不重叠，其对应的最大投影系数和基函数分别为

${\beta }^n=\{{\beta }_{I_1}^n,{\beta }_{I_2}^n,...\}$

$H^n=\{h_{a_{{\mu }_1},b_{{\mu }_1},I_1}^n,h_{a_{{\mu }_2},b_{{\mu }_2},I_2}^n,...\}$

∏ⁿ的连接方法应保证在本次分解中分解信号总能量最大，连接算法如下

1、初始化。对步骤2划分的动态时间支撑区进行编号，形成动态时间支撑区集合{I_m，m∈Z}。d(m)为连接到第m个动态时间支撑区之前分解信号的总能量，pre(m)为与第m个动态时间支撑区连接的前一个动态时间支撑区序号，e(m)为第m个动态时间支撑区内与最大投影系数对应的分解信号的能量，初始化时，置d(m)＝0，pre(m)＝0；

2、对于动态时间支撑区集合{I_m，m∈Z}中的每一个元素I_m，查找出与I_m相邻的所有下一个动态时间支撑区，形成集合{I_n}，即{I_n}中所有元素的起始时间与I_m相邻。稀疏信号分解方法是根据投影系数大小进行分解，投影系数大的先分解，小的后分解，但当多个分量具有相同的幅值时，其具有相同的投影系数，出现交叉分解的现象，为解决等振幅分解的问题，本发明在支撑区连接算法中引入了保留系数δ，即保留投影系数相差不大的基函数参与连接，如果

d(m)+e(m)＞d(n)×δ

有

d(n)＝d(m)+e(m)

pre(n)＝m

通常δ取1，但如果分量信号频率出现交叉解调的现象，则将δ的值逐次减小后重新分解，直到分量信号既光滑连续又能将等分量信号分解出来。

为了保证分解信号分量频率变化的连续性和光滑性，相邻动态时间支撑区在连接点的频偏相差在一定阈值之内，即$a_{{\mu }_m}-a_{{\mu }_n}<ϵ,$ε为偏差阈值，且I_m，I_n相邻，这样可以保证在分析信号有多个幅值相同的分量时单次分解只会选择一个分量成分进行连接。ε的取值影响分解信号分量频率曲线的光滑性，取值越大频率曲线的光滑性越差，反之则越光滑，因为本发明方法采用线调频基作为基函数，其频率曲线为倾斜的直线，所以相互连接的两个动态时间支撑区的频率成分不连续，会出现一定的频率跃变的情况，ε过小会导致基函数之间无法连接，因此ε应取能光滑分解的最小值。通常初始时，ε取10Hz进行分解。如果能连续分解，则逐次减小ε的取值，直到能光滑分解出频率连续的分量信号，优先选择相邻动态分析时间段频率变化较平坦的频率成分。连接动态时间支撑区下的分解信号的流程图如图2所示。

7)从加速度信号中减去分解信号分量，形成残余信号；原始信号与初次分解信号的关系为：

$f\left(t\right)=\underset{I_i\in \Pi }{\Sigma }[c_{I_i}\left(t\right)+r_{I_i}\left(t\right)]$

$r^l=\underset{I_i\in {\Pi }^l}{\Sigma }{r}_{I_i}\left(t\right)=f\left(t\right)-\underset{I_i\in {\Pi }^l}{\Sigma }{c}_{I_i}\left(t\right)$

$c^l=\underset{I_i\in {\Pi }^l}{\Sigma }{c}_{I_i}\left(t\right)$

至此便完成了信号的第一次分解，r^l为第一次分解的残余信号，c^l为分解信号。基于多尺度线调频基的稀疏信号分解方法是一个逐次分解的过程，本次分解的残余信号可以作为下轮分解的分析信号进行进一步的分解，直到残余信号能量小于一定的阈值便停止分解。分解到第n次，有

$r^n=\underset{I_i\in {\Pi }^n}{\Sigma }{r}_{I_i}\left(t\right)=r^{n-1}-\underset{I_i\in {\Pi }^n}{\Sigma }{c}_{I_i}\left(t\right)$

$c^n=\underset{I_i\in {\Pi }^n}{\Sigma }{c}_{l_i}\left(t\right)$

可以证明随着n的增大，残余信号rⁿ的能量迅速衰减到0。分解方法可以逐次分解出包含在信号中的信号分量，剔除了FFT变换中谐波成分的影响，也不会如二次时频表示一样出现交叉干扰项，能很好的分解出多分量非平稳信号。

8)判断残余信号能量与分析信号能量之比是否小于终止分解阈值，若大于，则将残余信号作为新的分解信号重复4～8步，如果小于终止阈值则停止分解，转到步骤9；

9)根据分解信号分量判断齿轮故障及其位置。实际齿轮箱系统中，当齿轮出现故障时，调幅调频效应通常是同时存在的，会出现啮合频率被转频或其高次倍频调制的调制现象，考虑多倍频调制的情况为

式中，A为信号幅值；B为调幅的调制指数；k为倍频系数，k＝1、2、3...；Z为齿轮齿数；f_n为非平稳转速函数；为频率调制函数；

根据积化和差公式，f(t)可以分解为三个信号分量

采用本发明方法，可以分解出f(t)的上述三个信号分量，从而得到信号的三个频率成份为

从上述三式可得到随时间变化的调幅调制频率kf_n(t)为：

kf_n(t)＝f₂(t)-f₁(t)

因此本发明方法可以获得非平稳转速下故障齿轮的幅值调制频率随时间的变化情况，而当齿轮发生故障时，其振动信号幅值调制频率函数通常为故障齿轮所在轴的转频或其高次倍频，所以通过对比分解得到的幅值调制频率函数kf_n(t)与变速箱各轴转频及其高次倍频，可以判断出故障齿轮所在轴及其调制倍频系数k，通常啮合频率(Zf_n)远远大于频率调制函数的频率所以将分解得到的第一个频率成份f₁(t)除以f₂(t)将近似得到齿轮的齿数Z：

Z≈f₁(t)k/[f₂(t)-f₁(t)]

因此本发明方法不仅可以判断故障齿轮所在轴，还可以获得故障齿轮的近似齿数，非常适合于非平稳转速下的齿轮箱故障诊断。

附图3为断齿齿轮振动信号时域波形图。将齿轮箱故障试验台上的主动齿轮人为切割一个齿，模拟齿轮断齿故障，输入轴齿轮齿数55，输出轴齿轮齿数75，因此啮合频率为转频的55倍，当齿轮箱发生断齿故障时会形成以齿轮啮合频率及其高次谐波为载波频率，以齿轮所在轴转频及其高次谐波为调制频率的齿轮啮合频率调制现象，从而产生调制边频带。采集齿轮箱振动加速度信号，采样频率为4096Hz，采样时长为1.9998秒，在非恒定转速下采集一组断齿振动信号和一组正常齿轮振动信号，正常齿轮和断齿齿轮参数相同。

从时域波形图上看出，在断齿的地方存在冲击现象，但冲击的时间间隔并不均匀，说明了在转速波动下齿轮振动冲击信号的非平稳性。

图4为经FFT得到的信号频谱。从图4无法直接识别调制边频带，从而无法判断齿轮的故障类型和位置。

图5、图6为对断齿振动信号进行基于多尺度线调频基的稀疏信号分解结果。取尺度系数0～4，则最小分析点数为8192/2⁴，即512点。因为信号的最大频率成分不超过采样率的一半，所以分析频偏范围取0～2048Hz，搜寻频偏分辨率为1Hz，搜寻调频率范围取-500～500Hz，搜寻分辨率为1Hz。经第一次分解，得到信号的啮合频率分量，如图5曲线1所示，曲线2为通过转速度计测量后乘以齿数获得的啮合频率估计，图中曲线1和曲线2基本重合，验证了基于多尺度线调频基的稀疏信号分解方法对非平稳信号分解的有效性。对残余信号进行第二次分解，分解信号的频率曲线如图5曲线3所示，曲线3在啮合频率曲线1的上方，且其几乎平行于啮合频率，图5中曲线3与曲线1频率之差与转轴的2倍频重合，如图6所示，图6中曲线1为转速的两倍频曲线，曲线2为图5中曲线3与曲线1频率之差，即调制频率曲线，图6中调制频率与2倍转频重合，进一步根据公式，得到故障齿轮的齿数为55.0349，与故障齿轮齿数55近似相等，所以通过本发明方法可以获得故障齿轮的转频和故障齿轮的齿数，非常适合于变速齿轮箱故障诊断。

图7、图8为对正常齿轮进行基于多尺度线调频基的稀疏信号分解得到的结果图，分解参数与断齿齿轮分解参数相同。经过第一次分解将得到信号的啮合频率分量，如图7曲线1所示，曲线2为通过转速度计测量后乘以齿数获得的啮合频率估计，图中曲线1和曲线2基本重合，验证了基于多尺度线调频基的稀疏信号分解方法对非平稳信号分解的有效性，同时也验证了正常齿轮箱振动信号的主要频率成份为齿轮的啮合频率分量。对残余信号进行第二次分解，分解信号的频率曲线如图7曲线3所示，图中曲线1与曲线3频率之差即图8中曲线1，图8曲线2为转速的两倍频曲线，曲线3为转速的一倍频曲线，图8中曲线1前端与转频的两倍频重合，这是由于信号中噪声的影响，而信号中端与末端与转频的一倍频和两倍频都不重合，所以可以判断齿轮的啮合频率没有出现转频及其高次频的调制现象，齿轮为正常齿轮。

基于以上实验分析，证明了本发明方法非常适合于齿轮箱故障诊断，可以获得故障轴的转频及故障齿轮的齿数，也进一步验证了当齿轮出现故障时会出现啮合频率被转频及其高次倍频调制的调制现象，而正常齿轮信号的主要成份为啮合频率，无频率调制现象，可以据此判断齿轮箱是否正常。