基于MUDW和峭度的齿轮故障信号预处理方法

仝 蕊， 康建设， 李宝晨， 张星辉

(1. 陆军工程大学石家庄校区装备指挥与管理系， 河北 石家庄 050003；2. 陆军工程大学科研学术处， 江苏 南京 210007)

作为机械传动系统的核心部件，齿轮的损伤或失效会影响整个系统的稳定运行,而齿轮发生故障时的振动信号呈现非平稳、非线性等特征[1]。为提高故障诊断的准确程度，需要寻求一种更为有效的非线性信号处理方法。

典型的非线性信号处理方法有小波分析、经验模态分解(Empirical Mode Decomposition， EMD)、局域均值分解(Local Mean Decomposition，LMD)和形态学分析(Morphological Analysis，MA)等。小波分析有自身缺陷，在重构采样中会遗漏部分信息[2]。EMD虽具有自适应性，可根据信号特点进行分解，但存在模态混叠现象[3]。LMD与EMD类似，是一种能自适应处理非平稳信号的方法，它在一定程度上改善了EMD易模态混叠的缺陷，但存在端点效应抑制的问题。MA通过结构元素探针在信号中的移动来提取信号的冲击特征，具有抑制脉冲干扰能力较强的优点，但存在盲目选择结构元素和过于依赖相关先验知识的缺陷。GOUTSIAS等[4]针对以上问题提出了形态小波(Morphological Wavelet，MW)分解，将形态滤波方法引入小波多分辨率分解中，在信号处理上具有良好的降噪性能，但信号分解时“隔二抽取”方式会造成逐层信息减半的问题。ZHANG等[5]提出了一种形态非抽样小波(Morphological Un-Decimated Wavelet，MUDW)分解，它省去分解过程下抽样和重构过程上抽样，解决了MW在信号分解时逐层信息减半的问题；但MUDW分解以信号最高分解层的近似信号作为预处理结果，部分故障信息仍会丢失[6]。

基于此，笔者提出一种MUDW分解和峭度融合处理齿轮故障振动信号的方法。首先，应用MUDW对信号进行分解，在此基础上对MUDW初始参数进行分析，并利用网格搜索法对参数组合进行优选，以避免主观经验对信号预处理效果的影响；然后，以峭度作为冲击信号的评估指标来确定包络分析的最优频带，并对各分解层近似信号进行加权融合，以解决传统MUDW分解方法存在的部分信息遗漏问题，更好地提高特征信息比重；最后，利用仿真信号和齿轮故障的实测振动信号对该方法进行验证。

1 基于MUDW的信号预处理

1.1 MUDW分解的基础运算

假设集合Vi为第i层信号空间，Wi为第i层细节空间，T(·)为形态算子，ZHANG等[5]给出了传统MUDW分解方法的一般框架：

(1)

(2)

(3)

MUDW是基于非抽样算法和形态滤波算子构造的，其关键的运算在于形态算子的选择。T(·)通常为腐蚀、膨胀、开运算、闭运算4种形态运算的任意组合形式，常用算子有形态梯度算子、形态差值算子和混合算子[10-11]。因齿轮故障信号中包含着不同尺度的形态特征信息，仅用单尺度结构元素对信号进行MUDW分解不能有效提取多尺度特征信息，需进行多尺度形态差值滤波算子的MUDW变换，同时提取信号中的正负脉冲。所以，MUDW的基础运算可描述为[12]

(4)

f(xi)∘(i+1)g]，

(5)

(6)

式中：f(xi)为原始信号；“∘”和“·”分别表示开运算和闭运算；(i+1)g表示对结构元素g进行i次膨胀操作。

1.2 MUDW参数优选

与线性小波分解类似，MUDW的分解层数N和分解元素长度L也是需要注意的问题，通常是根据待提取的形态特征来选择，但缺乏科学方法指导，可考虑刻画一种指标来衡量能达到信号预处理良好效果的最优参数组合。引用小波能量熵的概念：将小波分析和信息熵结合，反映信号的能量分布信息。将小波能谱[13]表征为E，即信号函数x(t)在尺度为a时的能量值。x(t)经分解后在N层上的小波能量为

E=(E1,E2,…,EN)

，

定义相应的小波能量熵

(7)

定义相对小波能量熵Hw为衡量对特征信息利用程度的指标：

(8)

Hw越小，信号成分越单纯，特征信息越明显，信号处理效果越好[14]。

在此基础上采用网格搜索法[15]对(N，L)参数组合进行优选，具体步骤如下：

1) 初始化网格搜索中的搜索范围和搜索步长。设置N的搜索范围，搜索步长为1；L的最大值不超过信号的脉冲周期，结合参数N可确定L的搜索范围，搜索步长为1；在N和L的坐标系上构造一个二维网格。

2) 计算所有网格点的目标函数值。即相对小波能量熵Hw，根据最小Hw评价当前参数。

3) 将当前参数存放于记忆器中，若满足Hw最小，则搜索结束，将Hw最小的那组(N，L)作为最优解。

2 基于峭度的故障信号融合

2.1 MUDW融合指标构建

通过分析MUDW分解算法可知：分解后的每一层近似信号都不同程度地包含着冲击信息。为了充分利用这些信息，以确保融合过程中能有效提高故障信息含量并减少干扰噪声，需要选择能敏感反映冲击信息的融合指标。

反映冲击程度的常用指标有均方、脉冲、峭度和特征频率幅值。其中，峭度是时域统计指标中无量纲的特征分析值，在机械设备盲信号处理过程中是非高斯性的自然度量指标，它反映了信号概率密度函数分布与高斯分布的偏离程度。峭度对冲击信息比较敏感，能够有效判断故障振动信号的冲击程度，故障越严重(冲击信号越明显)，偏离程度越高，峭度值越大。因此，为了突出信号中的冲击特征，有效筛选出对冲击贡献大的信号频段，笔者采用峭度[16]作为评价指标对各频段进行加权融合。

计算所有频带包络信号的峭度值，最大峭度值对应的频带则被确定为包含冲击信号强度最大的频带。峭度较好地反映了振动信号中冲击能量的大小，是归一化的4阶中心矩，其计算公式为[16-17]

(9)

分解后不同频段信号峭度值反映了本段近似信号的冲击程度，峭度值越大，信号中冲击脉冲信号所占的比重越大[18]，以此衡量含故障信息的贡献度，计算MUDW融合权值。

2.2 MUDW融合权值的计算及信号重构

假设MUDW各分解层近似信号为xi(i=1,2,…,N)，其对应的峭度为Ki，则其融合权值为

ki=Ki∑Ki。

(10)

由Ki的特性分析可知：Ki值越大，其对应的近似信号对特征的贡献度也相对越大，对应的ki也越大。因此，加权融合后重构的信号

(11)

因重构信号包含了各分解层特征信息，因此较融合前信号内包含的特征信息量得到了有效改善，提高了信噪比。

MUDW分解后，利用峭度指标进行加权融合的振动信号预处理流程如图2所示。

3 仿真信号分析

为验证方法的有效性，采用仿真信号模拟齿轮局部异常时的振动信号。设采样频率fs=2 048 Hz，采样时间t=0.5 s，采样长度为1 024的仿真信号

x(t)=x1(t)+x2(t)+n(t)，

(12)

式中：x1(t)模拟故障产生的转频冲击信号，其冲击频率f0=32 Hz；x2(t)模拟正常的啮合振动信号及其他传递到测点的信号，其齿数为20，将正常的啮合振动信号设为sin(2π×640t)，其他部件的振动信号设为0.5(sin(2π×15t)+sin(2π×100t))；n(t)模拟白噪声。

仿真信号x(t)的时域、频域图分别如图3、4所示。由图4可见：信号故障特征频率32 Hz基本淹没在噪声中，故障特征频率和噪声混杂在一起。因此，需利用本文所提出的方法对故障信号进行预处理。

3.1 基于网络搜索法的参数(N，L)优选

首先,确定N的搜索范围。因由MUDW分解得到的信号最小带宽只有大于特征频率的3倍[18]，才能确保其包含的冲击信号足够长，且分解重构后信号与原信号长度保持一致，由此确定其最大分解层数不会过大。由于MUDW的离散性，其分解层数N与信号长度相关，即信号采样长度为2N。分解初始层数由2开始，在分解到最大阶次时伸缩的小波母函数长度不应大于待分析信号的长度，因仿真信号采样长度为1 024，所以N最大值不能超过10，由此设置N的搜索范围为[2，10]。

然后，确定L的搜索范围。每个参数N对应L的一个搜索范围，根据前面分析并参考文献[6]，L搜索范围为[Nmin+1,⎣(fs/f0+N-1)/N」]，其中⎣·」表示向下取整运算，且fs=2 048 Hz，f0=32 Hz，在坐标系上构造(N，L)参数组合的二维网格。

最后，在(N，L)二维网格内，根据式(8)利用最小Hw评价当前参数，将当前最优参数存放于记忆器中。

在(N，L)数值组合的搜索范围中，N可取2～10的任意值，且都能对应得到L的搜索范围，具体变化情况如表1所示。

表1 N和L搜索范围

图5为不同(N，L)参数组合下Hw的数值变化情况。由式(8)计算出每个参数组合下所对应的Hw，N=2，3，…，10时对应的最小Hw分别为0.051 8、0.041 5、0.041 8、0.035 3、0.036 8、0.040 5、0.040 6、0.042 9、0.043 1。计算得到当N=5，L=7时，minHw=0.035 3，说明此时MUDW分解处理效果最好，信号构成成分简单，特征信息成分突出。

3.2 基于MUDW和峭度的振动信号预处理

当N=5，L=7时，设结构元素g0={0，0，0，0，0，0，0}，由式(4)-(6)对信号进行分解。由式(9)计算各层近似信号的峭度分别为：K1=13.301,K2=10.667,K3=11.393,K4=13.805,K5=10.336。根据式(10)计算相应的融合权值分别为：k1=0.223 5，k2=0.179 3，k3=0.191 5，k4=0.232，k5=0.173 7。根据式(11)进行加权融合，可得重构信号

xFinal= 0.223 5x1+0.179 3x2+0.191 5x3+

0.232x4+0.173 7x5。

图6为MUDW融合峭度指标对信号预处理后的频谱图。可以看出：与图4原始信号相比，能够清晰地看到特征频率f0=32 Hz及其2倍频、3倍频，说明噪声及谐波干扰得到了有效抑制。

4 试验验证

为验证本文所提方法的实用性，将该方法应用于齿轮裂纹故障振动信号预处理中。将齿轮齿根裂纹故障加工在输出轴大齿轮齿根上，进行齿轮齿根裂纹故障试验。齿根裂纹采用线切割以角度α进行加工，以深度5 mm进行试验。齿根裂纹加工位置及宽度如图7所示。

在机械振动及故障模拟平台上进行试验，如图8所示，主要设备有二级平行轴变速箱、4 kW三相电磁调速电动机和风冷磁粉制动器。电机输出轴与变速箱输入轴通过联轴器联接，在联轴器的右侧安装了一个转速转矩传感器，轴旋转一转产生60个脉冲。将4个3056B4型压电加速度传感器安装在变速箱上，利用数据采集板及LabVIEW软件将采集到的振动信号存入电脑。

该变速箱齿轮的齿数分别为：低速轴齿轮(故障齿轮)齿数Z1=81，中间轴大齿轮齿数Z2=64，中间轴小齿轮齿数Z3=19，高速轴齿轮齿数Z4=35。输入轴(轴3)转速为800 r/min，传感器位置表示为①、②、③、④，传感器位置及变速箱内部结构如图9所示。试验负载为10 N·m，采样频率为20 kHz，信号采样长度为1 024，计算可得轴3转频f3=800/60=13.33 Hz。设输出轴(轴2)转频为f2，平行轴齿轮传动比i1=Z4/Z2=f2/f3，则计算可得f2=7.29 Hz，一级啮合频率fm1=f1×Z4=f2×Z2=466.6 Hz，轴2转速为f2×60=437 r/min。设故障齿轮轴(轴1)转频为f1，齿轮传动比i2=Z3/Z1=f1/f2，则计算可得f1=1.71 Hz，二级啮合频率fm2=f2×Z3=f1×Z1=138.5 Hz。一般来说，具有局部异常故障的齿轮(如裂纹故障)将以转频为主要频率特征，主要分析其啮合频率及其边频带。齿轮箱含裂齿时，其振动能量将大幅增加，啮合频率及其谐波周围会产生边频带，边频带宽且高。

以传感器③ 6 s采集的振动数据作为原始样本，其所采集的故障振动信号的时域和频域图(功率谱)分别如图10、11所示。可见：故障信号存在明显的调制现象，故障特征频率淹没在噪声干扰中，且信号特征不明显。

采用本文所提方法对齿轮故障信号进行预处理，具体如下：

1) 采用网格搜索法计算分解层数N和分解元素长度L，根据采样长度可得N的搜索范围为[2，10]，N=2，3，…，10时对应的最小Hw分别为0.0573、0.050 8、0.041 3、0.034 8、0.041 7、0.043 1、0.040 1、0.042 1、0.048 2。计算得到当N=5，L=4时，minHw=0.034 8，表明此时MUDW预处理效果最佳，信号构成成分简单，特征信息成分突出。

2)设置结构元素g0={0，0，0，0}，由式(4)-(6)对信号进行分解。由式(9)计算各层信号的峭度分别为：K1=6.027 8，K2=4.226 2，K3=5.131 1，K4=5.036 5，K5=4.028 6。根据式(10)计算相应的融合权值分别为：k1=0.246 5，k2=0.172 8，k3=0.209 9，k4=0.206，k5=0.164 8，根据式(11)进行加权融合，可得重构信号

xFinal= 0.246 5x1+0.172 8x2+0.209 9x3+

0.206x4+0.164 8x5。

基于MUDW和峭度的故障信号预处理效果如图12所示。

由图12可见：通过MUDW分解和采用峭度指标加权融合对故障信号进行预处理后，能较清晰地看到输入轴转频f3=13.33 Hz、故障齿轮的啮合频率fm2=138 Hz及其倍频，以及边频带中含故障特征频率(裂齿转频)的倍频，原振动信号所存在的噪声和谐波干扰现象得到了有效抑制。

由于瞬态冲击能量大，激励起齿轮固有频率，产生了齿轮转频调制现象，出现了以齿轮啮合频率为中心频率、以裂齿所在轴的转频及其高次谐波为调制频率的调制边频带。对比正常齿轮振动信号并对信号进行预处理，效果如图13所示。可以看出：正常齿轮的振动信号经预处理后可以清晰地观察到输入轴转频，由于齿轮未发生故障，因此其他转频及啮合频率并没有明显的冲击信号。对比图12可知：故障齿轮的啮合频率冲击信号较明显(138 Hz)，可以较好地判断出齿轮发生故障的位置。

为了进一步验证加权重构方法的有效性，与传统MUDW信号预处理方法进行对比，效果如图14所示。可见：传统MUDW信号预处理方法受噪声干扰，且部分特征频率不明显或淹没在噪声中，同时遗漏了部分特征信息(如故障轴啮合频率2倍频等)，导致整体含特征信息比重相对较小；而基于MUDW和峭度的信号预处理方法因采用融合处理及参数优选，有效利用了包含在各分解层的特征信息并利用权重突出了冲击信号，比传统MUDW信号预处理方法特征信息比重大且故障特征明显，处理效果更清晰有效。

5 结论

针对齿轮故障振动信号不平稳、非线性强，导致处理效果不理想的问题，笔者在MUDW分解运算的基础框架下，提出了一种基于MUDW和峭度融合处理齿轮故障振动信号的方法，通过对信号各分解层的加权融合，有效提高了特征信息的比重，同时对其中的参数进行了优选，避免了主观经验对信号预处理效果的影响，最后利用仿真信号及齿轮实测信号验证了方法的有效性和实用性。结果表明：该方法对齿轮故障振动信号的预处理效果理想，能有效提取故障信号特征，为后续的故障诊断研究奠定了基础。

参考文献：

[1] 苏勋文,王少萍.齿轮故障演化的非线性动力学机理与实验研究[J].机械传动,2015,39(2):19-24.

[2] SARAVANAN N,RAMACHANDRA K I.Incipient gearbox fault diagnosis using Discrete Wavelet Transform (DWT) for feature extraction and classification using Artificial Neural Network (ANN) [J].Expert systems with applications,2010,37(6):4168-4181.

[3] FREI M G,OSORIO I.Intrinsic time-scale decomposition:time-frequency-energy analysis and real-time filtering of non-stationary signals[C]∥Proceedings of the Royal Society A:Mathematical,Physical & Engineering Sciences.England:The Royal Society,2007,463(2078):321-342.

[4] GOUTSIAS J,HEIJMANS H J A M.Nonlinear multi resolution signal decomposition schemes,part 1:morphological pyramids [J].IEEE transaction on image processing,2000,9(11):1862-1876.

[5] ZHANG J F,SMITH J S,WU Q H.Morphological un-decimated wavelet decomposition for fault location on power transmission lines [J].IEEE transactions on circuits and systems,2006,53(6):1395-1402.

[6] 孙健,李洪儒,王卫国,等.一种基于WMUWD的液压泵振动信号预处理方法[J].振动与冲击,2015,34(21):93-99.

[7] 章立军,阳建宏,徐金梧,等.形态非抽样小波及其在冲击信号特征提取中的应用[J].振动与冲击,2007,26(10):56-59.

[8] COUSTY J,NAJMAN L,DIAS F,et al.Morphological filtering on graphs[J].Computer vision and image understanding,2013,117(4):370-385.

[9] 黄兵峰,沈路,周晓军,等.基于形态非抽样小波分解的滚动轴承故障特征提取[J].农业机械学报,2010,41(2):204-207.

[10] 沈路,周晓军,张杰.形态非抽样小波与灰色关联度的轴承故障诊断[J].电子机械工程,2012,28(5):60-64.

[11] HE W,JIANG Z N,QIN Q.A joint adaptive wavelet filter and morphological signal processing method for weak mechanical impulse extraction[J].Journal of mechanical science and technology,2011,24(8):1709-1716.

[12] 孙健,李洪儒,王卫国,等.基于形态非抽样融合与DCT高阶奇异熵的液压泵退化特征提取[J].振动与冲击,2015,34(22):54-61.

[13] 申弢,黄树红,杨叔子.旋转机械振动信号的信息熵特征[J].机械工程学报,2001,37(6):94-98.

[14] 艾延廷,付琪,田晶,等.基于融合信息熵距的转子裂纹-碰摩耦合故障诊断方法[J].航空动力学报,2013,28(10):2161-2166.

[15] 吕中亮,汤宝平,周忆,等.基于网格搜索法优化最大相关峭度反卷积的滚动轴承早期故障诊断方法[J].振动与冲击,2016,35(15):29-34.

[16] 郭庆丰,王成栋,刘佩森.时域指标和峭度分析法在滚动轴承故障诊断中的应用[J].机械传动,2016,40(11):172-175.

[17] WANG D,TSE P W,TSUI K L.An enhanced kurtogram method for fault diagnosis of rolling element bearings[J].Mechanical systems and signal processing,2013,35:176-199.

[18] ZHANG X H,KANG J S.Rolling element bearings fault diagnosis based on correlated kurtosis kurtogram [J].Journal of vibroengineering,2015,17(6):3023-3034.