基于线调频小波路径追踪阶比循环平稳解调的齿轮故障诊断

1、基于线调频小波路径追踪阶比循环平稳解调的齿轮故障诊断    基于线调频小波路径追踪阶比循环平稳解调的齿轮故障诊断#陈向民，于德介，罗洁思*基金项目：国家自然科学基金(50875078)、国家高技术研究发展计划(863 计划): (2009AA04Z414)和高等学校博士学科点专项科研基金（20090161110006）作者简介：陈向民，（1984-），男，博士研究生，主要研究方向：信号处理与机械故障诊断通信联系人：于德介，（1957-），男，教授，博士研究生导师，主要研究方向：机械系统动力学、设备状态监测与故障诊断、动态信号分析与处理、振动与噪声控制.

2、E-mail:  调，既可以解幅值调制信号，也可以解频率调制信号，近年被引入到机械故障诊断领域，可以有效提取淹没在噪声中的周期性故障特征2-7。但在变转速情况下，如在风力发电齿轮箱中8，齿轮箱中转轴的转频随着风力强弱的改变而变化，信号中转轴的转动周期为一时变函数，因而其自相关函数也为时变的非周期性函数，不满足二阶平稳性条件（信号的相关函数45 随时间周期平稳变化），因此，需预先对信号进行周期平稳化处理。工程实际中常采用阶比跟踪方法来将非周期平稳信号平稳化，其基本思想是通过等角度采样技术将时域非周期平稳信号转换为角域周期平稳信号9-11。常用的阶比跟踪方法有硬件阶比跟踪法、计算阶比跟踪

3、法、基于瞬时频率估计的阶比跟踪法。相对于硬件阶比跟踪法和计算阶比跟踪法，基于瞬时频率估计的阶比跟踪法不需要使用专门的角度编码盘、跟踪滤波50 器和转速计等硬件以及上述硬件设备的安装条件，因而得到了广泛的应用。峰值跟踪法是目前常用的基于瞬时频率估计的阶比跟踪法12-13，其基本思想是先通过时频分析（如Gabor变换、小波变换等）将分析信号表示为时间和频率的联合函数，进而保留每个时刻联合函数最大值所对应的频率值以获取瞬时频率估计，最终依据该瞬时频率曲线对待分析信号进行插值重采样以实现等角度采样。但由于峰值跟踪法在信噪比较低的情况下估计效果不太理想，55 同时还存在时频分析法的优化选择、非阶比信号和

4、高次阶比谐波的影响等问题需要解决,其应用范围有限14。Emmanue 等近年提出了线调频小波路径追踪算法15，该算法通过对线调频小波图中的线调频小波原子进行连接，自适应的获得频率呈曲线变化的信号分量，己在地震引力波的分析中得到应用。最近，线调频小波路径追踪算法已被引入到机械故障诊断领域16-17，该方法60 能够精确地估计齿轮的啮合频率，且具有很高的抗噪能力。相对于峰值跟踪算法，线调频小波路径追踪算法表现出更高的频率拟合精度和更高的抗噪能力。针对变转速齿轮箱的实际工况，本文将线调频小波路径追踪算法与循环平稳解调算法相结合，提出基于线调频小波路径追踪的阶比循环平稳解调方法，并将其应用于变转速工况

5、下的齿轮故障诊断中。该方法首先通过线调频小波路径追踪算法求取信号中齿轮的啮合频率，65 用齿轮的啮合频率除以齿数便得到转轴的转速；然后，依据该转速对原始信号进行等角度重采样，得到角域周期平稳信号；再求取角域重采样信号的循环自相关函数，并在特征循环阶比处对循环自相关函数进行切片解调，便得到阶比解调谱；最后，根据阶比解调谱进行故障诊断。算法仿真和应用实例表明，本文方法能有效地提取齿轮故障特征。1 循环平稳解调70 如果非平稳随机信号x(t)的统计特性（相关函数、功率谱等）随时间周期或多周期（各个周期不能通约）平稳变化，则称x(t)为循环平稳信号。考虑信号x(t) 非对称形式的自相关函数1：( ,

6、) ( ) *( ) x R t =E x t x t+ (1)假定信号满足时间遍历性，对信号用采样周期0 T 进行采样，则( , ) x R t 可用时间平均表75 示为：*0 0( , ) ( ) ( )lim 1 ( ) ( )2 1xNN n NR t E x t x tx t nT x t nTN = += + + + (2)由式(2)可知，相关函数( , ) x R t 是关于时间t （ 固定）周期为0 T 的函数，因此，可将 ( , ) x R t 用傅里叶级数展开：( , ) ( ) j2 tx x R t R e = (3)80 相应的傅里叶系数为：/2 * 2/2

7、( ) lim 1 ( ) ( ) T j tx T TR xt x t e dtT = + (4)式中T 为整个采样区间，即0 T=(2N+1)T， 为循环频率（使( ) 0 x R 的 值），( ) x R 称为循环自相关函数，是关于循环频率 和时延 的函数。循环频率 包括零循环频率和非零循环频率，零循环频率刻画信号的平稳部分，非零循环频率对应信号的循环平稳85 性。对于调幅信号：( ) 1 cos(2 )cos(2 ) n z xt =A +B ft ft (5)式中n f 为调制频率， z f 为载波频率。将式(5)代入式(4)，积分运算时利用式(6)：1 2/2/2 1 2lim 1

8、 0 ( ) T j t j tT Te e dtT = (6)90 整理得到的积分结果为：2 2222 222 2222 2 ( )2 22 ( )cos(2 )1 cos(2 ) 02 2cos(2 )(1 )2cos(2 ) 2( ) 81 cos(2 ) 24 2 ( 2 )8(216nnzz n zn zz nj fz nj fz nxj fn zj f j f fn zj f fnA f B fA B f e fA B f e fRAe B f fA Be e f fA Be f ±±±± ± ±± ±

9、+ =+ =±= ±=+ =±+ =± ±= ± ±2 ) z f(7)从式(7)中可以看出，循环自相关函数的非零值只存在于循环频率等于调制频率及其倍频、二倍频载波频率、调制频率和载波频率的倍数和差等地方，其它频率值处均为零。对式(7)中的循环自相关函数，取定一个 值，即可得到一个自变量为时延 的切片，对该切片进95 行解调分析即可得到信号的故障特征。2 线调频小波路径追踪算法线调频小波路径追踪算法采用的多尺度线调频基元函数库如下：( ) ( ) 1 ( ) ( ), , , , , ,2 D h h t K e t Ii

10、a t b ta b I a b I a b I + = = (8)式中：100 D 为基元函数库；( ) , , h t a b I为多尺度线调频基元函数;I 为动态分析时间段,I =kN2 j (k +1)N2 j ，其中N 为分析信号的采样长度，j 为分析尺度系数， 0，1，.，log（2 1） j = N ， k = 0,1,L,2 j 1； Ka , b , I 为归一化系数，使得ha ,b ,I =1 ； a 为频率偏置系数， b 为调频率;根据采样定理a b t + 2 应该小于2 s f ；1 (t) I 105 为矩形窗函数，当t I 时为1，当t I 时为0。式(

11、8)定义的多尺度线性调频基函数在动态分析时间段内的瞬时频率为a b t + 2 。通过多尺度线性调频基函数对信号进行逐段投影分析，计算获得每个时间分析段I 内的最大投影系数和对应的线调频基元函数，该基元函数即为在时间分析段I 中与分析信号最为相似的频110 率成分。线调频基函数的多尺度特性使得它具有动态匹配分析信号的特性，而基函数中包含的调频率信息则使得其适合分析频率呈曲线变化的非平稳信号。当信号与多尺度线性调频基函数越相似时，其投影系数也越大，基元函数的能量也越大，因此要求找到一种动态分析时间段连接方法，在该连接方法下满足在整个分析时间内使连接的所有基元函数信号的总能量最大，即：1 ( )

12、)( )(2 2, , +Ia b I I ti a t b tMax K e , ,. 1 2 115 I I I (9) 覆盖整个分析时间段，不重叠，其对应的最大投影系数和基元函数分别为： , , 1 2I I L = , , 1 1 1 2 , 2 2 a ,b ,I a b , I L H h h = (10) 的连接方法应保证在投影中使连接的基函数在整个分析时间段内的总能量最大，线调频小波路径追踪算法提出的连接算法如下：120 (1) 初始化。以i 为时间支持区序号， d(i) 为第i 个时间支持区之前分解信号的总能量，pre(i) 为连接到第i 个时间支持区的前置时间支持区序号，e

13、(i) 为第i 个时间支持区最大投影系数对应的分解信号的能量，初始化时，置d(i) = 0 , pre(i) = 0 ；(2) 对于动态分析时间段集合I ,i Z i 中的每一个元素i I ，查找出与其相邻的所有下一个动态分析时间段集合 j I ,即 j I 中所有元素的起始时间与i I 相邻。如果：125 d(i) + e(i) > d( j) (11)有：d( j) = d( j) + e(i) pre( j) = i (12) 的连接方法可以保证在整个分析时间段内基元函数组合形成的信号与分析信号最为相似，而基元函数在动态分析时间支持区 , ,. 1 2 I I 内的瞬时频率为i a

14、 b t i i + 2 ， i i t I ，130 对应的线性直线连接形成的频率曲线则是对分解所得信号的瞬时频率估计。3 基于线调频小波路径追踪的阶比循环平稳解调算法阶比跟踪算法需要利用转速信号对振动信号进行等角度重采样，得到角域周期平稳信号。基于线调频小波路径追踪的阶比循环解调算法的具体步骤如下：(1) 利用线调频小波路径追踪算法估计信号中齿轮的啮合频率，用该啮合频率除以齿数135 得到转速，用三阶多项式对转速进行拟合，则得到转速表达式：f (t)=at3+bt2+ct+d (13)(2) 确定阶次跟踪的最大分析阶次Dmax ；(3) 计算角度重采样的角度间隔 ，根据采样定理，采样率至少

15、应为最大分析阶次的两倍，所以： max max2 2D D140 = = (14)(4) 计算重采样后数据的长度N ：02 ( ) T N f t dt= (15)式中T 为时域采样的总时间；(5) 从下式计算等角度重采样的键相时标n T ：4 3 2 4 3 20 0 0 0max( 1,2, )4 3 2 2 4 3 2 n n n naT bT cT dT n aT bT cT dTn ND+ + + = + + + + = L145 (16)式中0 T 为时域采样开始时间；(6) 对时域等间隔采样信号进行等角度重采样根据所求出的键相时标n T ，利用Langrange 线性插

16、值公式对振动信号进行插值，求出振动信号在角域里对应于键相时标n T 的幅值。对于给定的插值节点n T 及对应的幅值( ) n x T ，Langrange 线性插值公式为：111( ) ( ) (i ) (i)( )n i n i i n ii ixT xt x t x t T t t T tt t+= + 150 (17)( ) n x T 即为角域重采样信号；(7) 由式(4)求取( ) n x T 的循环自相关函数( ) x R ；(8) 根据式(7)，在特征循环阶比 处对循环自相关函数( ) x R 进行切片,再对切片进行解调分析，得到切片解调谱，根据切片解调谱进行故障诊断。155 4

17、 算法仿真变转速齿轮箱振动信号为多分量非平稳信号，当齿轮箱出现故障时，会形成被故障齿轮所在轴转频或其高次倍频调制的调制信号，取仿真信号如下：( ) (1 cos(2 (20 1.429 cos(1.4 )cos(2 10 (20 1.429 cos(1.4 )x t t tt t = + × × × × × × ×× × × × × × ×(18)信号转频为：160 fc= 20+2× ×sin(1.4× ×t) (

18、19)模拟齿轮齿数为10，调制频率为1 倍转频。对仿真信号进行采样，采样频率为1024Hz，采样点数为1024，信号时长为1 秒，加入信噪比为0db 的高斯白噪声，信号波形如图1 所示：0 0.2 0.4 0.6 0.8-4-2024时间 t/s幅值 A165 图1 信号波形对图1 所示信号运用线调频小波路径追踪算法估计啮合频率，然后由啮合频率除以模拟齿数得到信号转速，求出的转速曲线如图2 中的曲线2 所示，从图2 可看出，线调频小波路径追踪算法能很好地逐段拟合频率变化曲线。 0 0.2 0.4 0.6 0.8152025时间 t/s频率 f/Hz211701精确转频曲线2基于线调频

19、小波路径追踪估计的转频曲线图2 仿真信号瞬时转频曲线175 由于信号的模拟齿轮齿数为10，对应角域阶比为10，根据角域采样定理，最大分析阶比max D 应大于信号的最高阶比，因而本文选取信号的最大分析阶比max D 为30 阶，利用已求取的瞬时转频对原始信号进行等角度重采样，得到的角域周期平稳信号如图3 所示：0 20 40 60 80 100 120-202弧度 rad幅值 A图3 角域周期平稳信号180求取角域周期平稳信号的循环自相关函数，因为转频对应的阶比为1，所以在循环阶比 =1处对循环自相关函数进行切片，得到的切片解调谱如图4 所示，可看出，在接近阶比1 处有明显的峰值，说明信号被1

20、 倍转频调制，与实际相符，从而验证了本文方法的有效性。0 5 10 15 20 2500.0050.010.0150.02阶比 n幅值 A0.9961185 图4 循环阶比 =1时的切片解调谱5 应用实例为进一步验证本文方法的有效性，分别对断齿、裂纹和正常齿轮振动信号进行分析。变转速齿轮箱试验台简图如图5 所示。现将主动齿轮切割一个齿，并在无负载条件下模拟齿轮190 断齿故障，输入轴齿数为55，输出轴齿数为75。实验中利用加速度计拾取信号，采样频率为4096Hz，采样点数为8192，采样时长为2 秒。在变转速下采集到的断齿齿轮振动信号时域波形如图6 所示： 图5 模拟实验台简图0 0

21、.5 1 1.5-0.500.5时间 t/s幅值 a/(m.s-2)195图6 断齿齿轮原始振动信号基于线调频小波路径追踪算法估计得到的断齿齿轮主动轴转频曲线如图7 所示：0 0.5 1 1.521.52222.52323.5时间 t/s频率 f/Hz200 图7 基于线调频小波路径追踪的断齿齿轮转频曲线利用图7 所示转速对信号等角度重采样，得到的角域周期平稳信号如图8 所示。求取该角域周期平稳信号的循环自相关函数，在循环阶比 =1处对循环自相关函数进行切片解调，得到的切片解调谱如图9 所示，从图可看出，在阶比1、阶比2 处有明显的峰值，说明该主205 动轴齿轮出现了故障。0 20 40 60

22、 80 100 120-0.4-0.200.20.4弧度 rad幅值 a/(m.s-2)图8 断齿齿轮角域周期平稳信号0 5 10 15 2000.51x 10-4阶比 n幅值 a/(m.s-2)12210 图9 断齿齿轮循环阶比 =1时的切片解调谱 同样，用一对模数相等，齿数均为37 的齿轮进行实验，并将主动齿轮齿根处用激光切割深度约为1.5mm 的裂纹以模拟齿轮裂纹故障。实验采样频率为5120Hz，采样点数为8192，在变转速下采集到的裂纹齿轮振动加速度信号如图10 所示。图11 为基于线调频小波路径追215 踪算法估计得到的主动轴转频曲线，利用该转频曲线对原始振动信号进行等角度重采样，得到的角域周期平稳信号如图1