基于最优窗函数Gabor变换的变工况行星齿轮箱故障诊断*

胡瑞杰 庞学博 佘彩青 刘艾强 胡少梁 李宏坤

（1.大连理工大学机械工程学院；2.中车大连机车研究所有限公司；3.西安航天动力研究所）

0 引言

风能是一种清洁无公害的可再生能源，蕴含巨大能量，目前日益受到世界各国重视，风电已经成我国第三大主力电源[1]。但由于风力发电机服役环境恶劣，长期处于变转速、变载荷及强瞬时冲击的工况下，其行星齿轮箱极易发生故障[2]。当行星齿轮箱发生故障时，其振动信号表现出明显的非平稳、非线性特点，且受故障、转速和传递路径等多种因素影响表现出调频、调幅和调相等特点[3]，给特征提取造成困扰。传统的时域与频域分析方法只适合提取恒定转速下平稳信号的特征，无法提取变转速下非平稳信号特征。时频分析是提取非平稳信号特征的一个重要方法，通过联合信号的时间与频率，对其故障特征进行表示[4]。近些年也有学者从自适应分解方法的角度出发研究变转速机械设备故障诊断问题[5]。

传统的时频分析方法有短时傅里叶变换(Shorttime Fourier Transform，STFT)、小波变换(Wavelet transform，WT)、Gabor 变换和魏格纳变换(Wigner-Ville Distribution，WVD)等等[6]，广泛应用于机械设备故障诊断、状态监测以及寿命预测等领域[7]。以STFT、WT、Gabor 变换为代表的线性时频分布通过内积来处理信号，从而确定各频率成分随时间变化的趋势。但由于其窗函数固定，缺乏自适应性。根据Heisenberg不确定性原理，这类方法的时间分辨率与频率分辨率不能同时达到最优[8]。以WVD 为代表的双线性时频分布根据能量守恒原理[9]，是一种利用时域与频域的概率密度函数来表示信号的能量二次型时频分布。这类分布改善了时频分辨率，但存在交叉干扰项影响特征提取[10]。关于交叉干扰项的抑制是目前时频分布的热点研究方向[10-11]。

在齿轮与轴承故障分析过程中，需要特别注意升降速阶段的信号。因为升降速过程比平稳过程包含了更多故障特征信息，分析该阶段信号更容易提取到难以发现的故障特征频率[12-13]。行星齿轮箱升、降速阶段输入轴转速信号可近似为线性上升或下降的直线，因此其振动信号可近似为多分量的线性调频信号。分数阶傅里叶变换(Fractional Fourier Transform，FrFT)作为傅里叶变换的广义形式，具有独特的时频旋转特性，可将时频平面旋转至信号自项具有最小宽度的方向，从而大幅提高时频分辨率，该方法十分适合处理线性调频信号[14]。

本文将FrFT 的优点引入到传统时频变换方法中，在分数域优化Gabor变换的时频窗函数，得到最小时间带宽积对应的窗函数，从而得到高质量的时频变换图像。利用该方法对行星齿轮箱升速阶段振动包络信号进行分析，若齿轮发生故障时，则在其时频图像上可提取明显的调制频率成分。仿真分析与实际信号分析都表明该方法具有良好的可行性与噪声鲁棒性，适合处理升降速阶段低信噪比振动信号。

1 FrFT变换及滤波原理

1.1 FrFT原理

一般的，信号x(t)的p阶FrFT 可以表示为Xp(u)或者FPx(t)，FPx(t)可以视为算子FP作用于信号x(t)，其结果在u域上的表示。FrFT的定义如下：

式中，Kp(t,u)为FrFT变换核。

式中，α=pπ/2为FrFT的旋转角度；n为整数。对应阶数的逆变换为：

FrFT可以理解为对任意信号在时频域内绕原点以α角度进行旋转构成的分数域上的表示，α可以为任意角度，坐标轴如图1(a)所示。由式(2)可知α只出现在三角函数位置，所以参数p∈(-2,2]。当p=0时，FrFT的结果为原函数，p=1时是普通的傅里叶变换。

齿轮箱加速或减速的振动信号非常接近多分量线性调频信号（Linear Frequency Modulation，LFM）[14]，FrFT非常适合对LFM信号进行自适应滤波。设某含有两个LFM分量的信号时频图像如图1(b)所示，其中某一分量的时频图像与时间轴的夹角为β。

图1 FrFT旋转特性Fig.1 FrFT rotation characteristics

由图1(b)可知，当α与β正交时，该LFM信号在u轴上的投影有最好的聚集性，而其他分量不聚集，此时设计窄带滤波可实现LFM分量与其他分量的分离。由于噪声分量的能量在时频域内均匀分布，在任何变换阶次FrFT 的结果中都不聚集。因此经过FrFT 滤波后信号具有很好的噪声鲁棒性。为了得到最佳的变换结果，首要任务是找到FrFT的最佳变换阶次p0。

1.2 寻找最优变换阶次p0

根据图1（b）可知，当已知LFM 信号的线性调频率k0时，最佳变换角α的计算公式为[15]：

实际采集信号往往无法得知线性调频率等参数。工程中常采用峰值搜索思想确定FrFT 的最佳阶次，即按照一定步长Δp对信号进行连续的FrFT，然后在分数阶域u与变换阶次p构成的(p,u)平面上寻找峰值。由于LFM 信号在分数域良好的聚集性，在(p,u)平面上会出现明显的峰值，如图2所示。此时按照阈值在平面上进行二维搜索，就可以确定LFM 信号对应的最佳变换阶次p0与聚集中心u0。

在实测信号处理过程中，变工况下采集的行星齿轮箱振动信号通常是多分量调频信号。若齿轮存在分布式或局部式故障，在啮合过程中振动信号会有明显的幅值和频率调制现象并且伴有冲击成分。这种多分量信号的FrFT 幅值谱中峰值众多，一个峰值附近往往伴随着其他突出的峰值，而且弱分量信号容易被强分量信号淹没难以提取特征频率，所以根据峰值搜索法寻找最佳变换阶次容易出现较大误差。

文献[16]提出采用转速信号来确认最佳变换阶次的方法，根据输入轴转速信号计算出齿轮各级齿轮啮合频率及其转频，对转频及各级齿轮啮合频率分量进行最小二乘法拟合，计算出各分量的调频率，此时可采用式(4)计算目标分量的最佳旋转角。由于输入转速不受振动和噪声的干扰，且齿轮箱内各级齿轮传动比恒定，所以根据该方法计算得到的调频率及变换阶次具有较高精度，并且具有较高鲁棒性。但是在实际信号采集过程中由于设备尺寸及空间限制，通常不容易获得输入轴转速信号。本文引入无量纲参数峭度，其表达式为：

式中：μ为信号x的均值；σ为信号x的标准差。

峭度是归一化4阶中心矩，它反映了信号的分布特性[17]。峭度对信号中非正态分布成分很敏感，可以用来检测冲击成分。当齿轮正常运行时，振动信号的幅值分布近似服从正态分布，峭度值约为3。当齿轮出现局部式或分布式故障时，振动信号幅值会偏离正态分布，峭度值也随之增大，且峭度值越大说明信号幅值越偏离正态分布，冲击成分所占的比重越大，而故障信息往往隐藏在这些冲击成分较多的调制信号中。因此选择峭度K作为选择最佳变换阶次的指标，计算FrFT 所有变换结果的峭度并利用峰值搜索则可以确定p0，如图3所示。某齿轮实测信号中冲击成分的峭度突出，在FrFT 分量峭度谱中搜索峰值可得最佳变换阶次p0为1.03。

图3 某齿轮实测信号FrFT分量峭度谱Fig.3 FrFT component kurtosis spectrum of a gear measured signal

1.3 FrFT与Gabor变换的联系

Gabor变换基本定义如下：

式中：f(τ)为原始信号；h(t)为Gabor变换窗函数。

本文中的定义为[18]:

在上式定义中，Gabor变换的幅值具有FrFT的旋转特性，但相位却不具有该性质：

对Gabor 变换定义进行修改使其相位也满足FrFT旋转特性：

文献[18]已经证明，当信号f(t) 的FrFT 结果为FP f(t)，f(t) 的Gabor 变换结果为的Gabor变换结果是时满足：

式中：Rα为时间与频率轴逆时针旋转因子。

由上式可以证明：信号f(t)的Gabor 变换结果等价于将FP f(t)进行Gabor 变换后再逆时针旋转α(α=pπ/2)角度。

2 分数域最优窗函数确定

2.1 广义时间带宽积

时间带宽积（Time-bandwidth product，TBP）是时间与带宽的乘积，是评估时频聚集性的重要指标，它能表示时频域内有效信号的支撑区域。时频图像的聚集性随着TBP的值下降而提高。Gabor变换的时间带宽积是指其加权窗函数的时间带宽积[19]，定义如下：

式中，h表示Gabor变换窗函数，Tx·h,Bx·h分别表示加权窗函数的时宽与带宽，Tx,Bx分别表示原始信号的时宽与带宽。Th,Bh分别表示窗函数信号的时宽与带宽。ηt,ηf分别表示原始信号时域与频域的均值。

为了获得更高质量的时频图像，Durak等提出了广义时间带宽积（Generalized Time-bandwidth product，GTBP）的概念，其定义式为[20]：

对于线性调频信号来说，TBP是其时频域内时间与带宽乘积的区域，如图4(a)所示。GTBP 的原理是将FrFT 引入到TBP 的计算过程中，利用FrFT 的旋转特性在其分数域时频轴上计算TBP，如图4(b)所示。这样改进后可以显著降低TBP，获得更时频聚集性更好的时频图像，为后续故障特征提取做好准备。

图4 GTBP与TBP对比Fig.4 The comparison of GTBP and TBP

2.2 确定最优窗函数

文献[20]已经证明，最小时间带宽积对应的Gabor变换窗函数表达式为：

结合本文1.3 节提出的Gabor 变换与FrFT 的关系，对Gabor变换窗函数进行优化，求解出在最小广义带宽积下的最优窗函数。首先对式(13)进行扩展：

因此，可推导出GTBP 原则下Gabor 变换的最优窗函数为：

式中,hminTBP(xp),p(t) 为hminTBP(xp)(t) 的p阶FrFT 变换结果;hGTBP(t)为GTBP原则下Gabor变换的最优窗函数。

3 算法流程

上述章节详细叙述了在分数域优化Gabor 变换窗函数的方法，进而获得能量聚集性与时频分辨率较好的时频分布。利用该时频图像对采集的行星齿轮箱升速过程中的振动信号进行故障特征频率提取，可以判断行星齿轮箱的健康状态。并且该方法具有良好的噪声鲁棒性，可以从严重的背景噪声中提取出故障特征频率，时频图像受噪声干扰程度很小。具体算法流程（图5）如下：

1）对采集的振动信号进行去均值、包络解调等预处理，得到包络信号。

2）对包络信号采用两级步长法进行p∈(-2,2]阶FrFT变换，并计算FrFT结果的峭度值，利用峰值搜索法寻找最大峭度值，对应的旋转因子即为p0。

3）计算加权窗函数时间带宽积TBP与广义时间带宽积GTBP。

4）利用式(24)计算出GTBP 原则下Gabor 变换的最优窗函数hGTBP(t)。

对包络信号进行分数域最优窗函数Gabor变换，得到能量聚集性较好的时频图像，利用该图像进行行星齿轮箱健康状态判断。

图5 本算法流程图Fig.5 Algorithm flow chart of this paper

有两点需要特别注意，第一是基于峰值搜索时信号FrFT峭度谱是否具有唯一峰值。与实信号的FFT结果一样，实信号的FrFT 变换结果也存在对称项。对称项会导致信号FrFT 峭度谱中的峰值不唯一，从而无法搜索到正确的p0。因为包络信号具有保留正频率部分、消除负频率部分的性质，因此，提出基于包络信号消除对称项的方法。即先由实信号得到包络信号，再对包络信号进行p∈[0,2]阶FrFT变换，即可消除对称项。第二是峭度搜索法确定p0时的步长问题。由于LFM信号只聚集在(p,u)平面有限范围内，考虑到运行时间与计算精度，可以采用两级步长的搜索方法确定p0。即先在p∈[0,2]范围内进行大步长Δp1=0.01 的连续FrFT，并对计算结果的峭度谱进行峰值搜索，初步确定(p01,u01)的位置，然后在区间p∈[p01-nΔp1,p01+nΔp1]内再进行Δp2=0.001 的连续FrFT，并对结果峭度谱进行峰值搜索，从而确定(p02,u02)。采用两级步长峰值搜索法不仅能保证精度，又能大幅提高效率。

4 仿真分析

为了验证所提算法在变转速工况下行星齿轮箱齿轮故障诊断中的有效性，现构建一太阳轮故障仿真信号。当行星齿轮箱中太阳轮发生故障时，会形成被太阳轮所在轴转频以及故障特征频率同时调制的信号。分析下式所示的太阳轮故障仿真信号[21]:

式中，fr=15+1.5t模拟行星齿轮箱升速过程中太阳轮转频；fs表示太阳轮断齿故障特征频率；fm表示啮合频率；调制系数A=0.9，B=1。初相位φ=θ=0。为了模拟实际运行情况下采集信号背景噪声的影响，在仿真信号中加入白噪声δ(t)，信噪比为-1dB。采样频率为1024Hz，采样时间为4s。

根据卷积理论，原始信号包含9 个时变成分[22-23]：fm、fm±fs、fm±fr、fm±fs±fr。其中fs=(10/3)fr，fm=(50/3)fr。主要调制成分为fs，fr及它们的组合fr±fs。采用本文提出算法处理仿真信号，结果如图6所示。

由图6(a)可知信号受噪声干扰严重，时域特征湮没在背景噪声里。对包络信号进行时频分析得到图6(d)、(e)、(f)。图6(d)是对包络信号采用传统Gabor 变换的结果，可以看出时频分布结果分辨率较差、时频聚集性较差，无法识别调制频率。图6(e)是对包络信号采用WVD变换的结果，时频聚集性较Gabor 变换结果有所提升，可以识别到调制频率随时间变化的趋势，但频率分辨率较差，时频模糊严重。且存在无法消除的干扰项，对特征提取造成影响。图6(f)为对包络信号采用本文算法得到的时频变换结果，可以看到该算法降噪效果显著，仿真信号中主要的调制频率fr，fs，fs+r，fs-r都可以提取出来，这些频率分量的最佳变换阶次p0依次为：1.6257，1.8748，1.9033，1.8232。图像的时频聚集性与时频分辨率比较传统的时频分布方法有较大的改善。因为图6(e)中的主要调制频率成分都与太阳轮故障特征频率fs有关，所以可以判断该行星齿轮箱中太阳轮发生了故障，此结果与所模拟的故障相一致，从而验证了本文算法在变转速行星齿轮箱故障诊断过程中的可行性和优越性。

图6 仿真信号分析结果Fig.6 Simulation signal analysis results

5 试验与工程信号分析

5.1 试验设置

为了验证本文算法在实际变转速工况下行星齿轮箱故障诊断中的可行性，以太阳轮断齿和行星轮缺齿为例，在行星齿轮箱故障诊断试验台上进行分析。该试验台如图7所示，主要包括驱动电机、变频调速器、一级行星齿轮箱、磁粉制动器、加速度传感器(型号为美国 DYTRAN 公司生产的3035B 型传感器，灵敏度为100mV/g)和安装NI9234采集卡的工控机组成。行星齿轮箱结构参数如表1 所示。通过线切割的方式分别加工太阳轮断齿故障和某一行星轮缺齿故障，如图8 所示。采样测点选择在行星齿轮箱外壳的轴向与径向的位置，采样频率设置为12800Hz，采样时间为10s，输入轴转频fr变化范围为15～30Hz线性上升。通过计算可得行星齿轮箱各齿轮故障特征频率[24]，如表2所示。考虑到运算速度与占用内存，对原始信号进行降采样处理，处理后信号采样频率为1024Hz。

图7 行星齿轮箱试验台Fig.7 Planetary gearbox test bench

图8 行星齿轮箱故障件Fig.8 Faulty parts of planetary gearbox

表1 行星齿轮箱结构参数Tab.1 Structural parameters of planetary gearbox

表2 行星齿轮箱故障特征频率Tab.2 Characteristic frequency of planetary gearbox fault

5.2 太阳轮断齿故障分析

对太阳轮断齿故障实测信号进行分析，结果如图9所示。由图9(a)～(d)可知原始信号受噪声干扰严重，无法从时域信号中提取特征，且变工况下故障特征频率随时间变化，无法从频谱图中提取出目标分量。采用时频分析方法对实测包络信号进行分析，如图9(f)～(h)所示。图9(f)为采用传统Gabor 变换结果，可以看出时频分辨率较差、时频聚集性较差，无法识别调制频率。图9(g)为采用WVD变换结果，时频聚集性较Gabor变换结果有所提升，可以识别到调制频率随时间变化的趋势，但频率分辨率较差，受噪声干扰严重。且存在无法消除的干扰项，对特征提取造成影响。图9(h)为采用本文方法处理包络信号得到的结果，可以看到该方法抑制噪声效果显著，时频分辨率与时频聚集性较传统时频分析方法有较大改善，提取到的主要调制频率成分为fr，2fr，3fr，(7/3)fs，2fs，这些频率成分的最佳变换阶次p0依次为：1.6257，1.7952，1.8608，1.9161 和1.9280。频率成分随时间的变化趋势和太阳轮转频相似，且出现了比较明显的太阳轮故障特征频率fs的倍频。这种现象说明行星齿轮箱太阳轮出现了故障，符合试验设置。

图9 太阳轮断齿信号分析结果Fig.9 The analysis results of the broken tooth signal of the sun gear

5.3 行星轮缺齿故障分析

对行星轮缺齿故障实测信号进行分析，结果如图10 所示。从原始信号与包络信号的时域、频域图中无法提取出明显的故障特征。图10(f)为对包络信号采用Gabor变换结果，时频分辨率较差、时频聚集性较差，无法识别调制频率。图10(g)为WVD变换结果，可以识别出部分调制频率成分随着时间变换，时频聚集性提高，但受噪声干扰严重，时频分辨率差，且存在无法去除的交叉干扰项，很难从时频图像中提取出有效成分。图10(h)为采用本文算法分析结果，可以看到时频聚集性与时频分辨率有了较大的提升，该算法抑制了大部分背景噪声，提取到的主要调制频率成分为fp与fr的组合：fr,6fp，9fp，5fr-(2/3)fp，6fr-(2/3)fp，16fp，7fr+fp，8fr+fp，9fr+fp，这些频率成分的最佳变换阶次p0依次为：1.6257，1.8350，1.8886，1.9107，1.9261，1.9369，1.9429，1.9497 和1.9550。频率成分随时间的变化趋势和太阳轮转频相似，除转频外所有频率成分都与行星轮故障特征频率fp相关，且出现了5fr-(2/3)fp，6fr-(2/3)fp等fp与fr的倍频组合的成分，振动信号受转频与故障特征频率同时调制的特点明显，可以判定是行星齿轮箱中某行星轮发生故障，符合实验设置。在图中也观察到行星齿轮箱固有频率fn，当各个频率成分增大至共振带时，振动能量明显增大。

图10 行星轮缺齿信号分析结果Fig.10 Planetary gear missing tooth signal analysis results

5.4 正常齿轮分析

为了验证所提算法的有效性，采集同试验台行星齿轮箱正常太阳轮变转速工况下振动信号进行对比实验，采样频率为12800Hz，采样时间为5s，输入轴转频fr变化范围是10～25Hz 线性上升。对原始信号进行降采样处理，处理后信号采样频率为1024Hz，分析结果如图11 所示。由图可知信号噪声严重，无法从时域图中提取出有效成分，且因为变转速的原因频谱图中频率模糊严重，无法提取到转频fr与啮合频率fm。对包络信号进行时频分析得到图11(f)～(h)。Gabor 变换与WVD的噪声鲁棒性差，包络时频图像中无法观察到任何调制频率分量，而本文算法分析结果中可以在低频段提取到一阶转频fr与三阶转频3fr，最佳变换阶次p0依次为1.7956和1.93，无明显的故障特征频率。考虑到行星齿轮箱的制造与安装误差，出现转频及其倍频对振动信号的调制属于正常现象。

图11 正常齿轮信号分析结果Fig.11 Healthy gear signal analysis results

6 总结

本文在充分研究FrFT 的基础上，将其处理线性调频信号的优势引入到时频分析方法中，在分数域优化传统Gabor变换的时频窗函数，创新性的利用峭度峰值寻找最优变换阶次p0，获得最小广义时间带宽积原则下最优窗函数，达到提高时频图像能量聚集性与时频分辨率的目的。最后，采用该算法对太阳轮故障仿真信号以及实测行星齿轮箱中太阳轮断齿及行星轮缺齿振动信号进行分析，结果表明本文方法可以在强背景噪声下提取出非平稳信号的时变特征，获得准确的诊断结果，验证了所提出方法的可行性与优越性。