改进经验小波变换与最小熵解卷积在铁路轴承故障诊断中的应用

乔志城，刘永强，廖英英

(1. 石家庄铁道大学 机械工程学院，石家庄 050043；2. 石家庄铁道大学 省部共建交通工程结构力学行为与系统安全国家重点实验室(筹)，石家庄 050043；3. 石家庄铁道大学 土木工程学院，石家庄 050043)

滚动轴承作为影响旋转机械正常运行的关键部件，对其进行状态监测与故障诊断十分必要。滚动轴承振动信号一般为非线性、非平稳信号，传统的Fourier变换只适用于频率不随时间变化的线性平稳信号，且容易受噪声干扰影响[1-2]。针对滚动轴承故障信号非线性、非平稳特性以及早期故障特征微弱，易受噪声干扰的特点，国内外学者提出经验模态分解(Empirical Mode Decomposition, EMD)的故障特征提取方法。EMD能够将复杂的轴承振动信号自适应分解为若干个本征模态函数分量，实现微弱故障信号与噪声的分离。但EMD本身存在严重的端点效应、模态混叠等缺陷[3-4]。文献[5]提出在信号中加入不同幅值的高斯白噪声的集合经验模态分解(Ensemble Empirical Mode Decomposition, EEMD)方法，利用高斯白噪声的均匀分布统计特性，使不同频率尺度的信号具有连续性，从而有效抑制了模态混叠[6]。文献[7]提出变分模态分解(Varianational Mode Decomposition, VMD)方法，该方法通过迭代方式搜寻变分模型最优解来确定各分量的中心频率与带宽，实现信号频域剖分及各分量的有效分离。相比于EMD，VMD具有完善的理论基础与噪声鲁棒性，能够有效抑制模态混叠[8]。但在使用VMD时需要确定模态分解个数K，K过大或过小都会影响信号分解精度[9-10]。

文献[11]提出了基于小波变换与窄带信号分析理论的经验小波变换(Empirical Wavelet Transform, EWT)。EWT与EMD，VMD等相比[12]，避免了模态混叠与端点效应，且没有迭代过程，不受参数设置的影响，具有计算量小，分解速度快等优点[13]。但EWT在频谱划分过程中，会出现频谱划分过多或者划分不合理等问题。针对这一问题，文献[14]提出根据频谱包络的划分方法，克服了频谱划分集中于高幅值频段的问题，但需要人工对分量个数进行估计，当面对复杂的频谱时仍会因为估计不准并产生分解过度或分解不足的现象；文献[15]提出时频峭度谱的频谱划分方法，但该方法对简单频谱有较好的效果，对复杂频谱的效果有待提升。

本文提出基于互信息的改进EWT频带划分方法，解决传统EWT频带分解过度的问题，根据重构分量的互信息与峭度选择最优分量进行信号重构，但此时的重构信号中包含较多噪声，本文采用最小熵解卷积方法(Minimum Entropy Deconvolution, MED)对重构信号进行降噪处理，然后对降噪信号进行包络分析，能够有效诊断出铁路货车轮对滚动轴承早期微弱故障，具有工程应用价值。

1 理论基础

1.1 经验小波变换

经验小波变换是一种建立在小波理论框架下的自适应小波分析方法。EWT首先对振动信号的频谱进行分割并构造一组具有适应性的小波滤波器组，然后对不同频率成分进行分析，提取具有紧支撑特性的调幅-调频信号[16]。

首先对振动信号进行Fourier变换得到频谱，将频谱定义在[0 π]，然后分割为连续的N段，每个频带可表示为Λn[ωn-1,ωn]。

根据Little wood-Paley和Meyer理论，经验小波的尺度函数φn(ω)和小波函数ψn(ω)在频域的定义为

(1)

(2)

式中：β(x)=x4(35-84x+70x2-20x3)； 0<γ<1；τn=γωn。

经验小波的细节系数和近似系数分别为

(3)

(4)

结合以上公式，得到振动信号的重构公式

(5)

式中： ＾为傅里叶变换； ˇ为傅里叶逆变换； -为复共轭； *为卷积。

EWT频带划分方法-尺度空间法：将频谱划分为N个频带，首先在频谱中确定N+1个边界点，除端点0和π外，还需确定N-1个点。将除0和π外的极大值点按从大到小顺序排列，设找到M个极大值点，此时对应两种情况：①当M≥N时，此时有足够数量的边界点，只保留前N-1个点；②当M≤N时，此时边界点数量不足，保留所有边界点的同时，对N值进行重置。

1.2 最小熵解卷积

最小熵解卷积本质：通过迭代方式，寻找一个最优逆滤波器，使包含故障信息的冲击脉冲信号得到增强，迭代终止条件为峭度值最大[17-18]。假设振动信号出现局部损伤时的振动信号为

y(n)=x(n)*h(n)+e(n)

(6)

式中：y(n)为输出信号；x(n)为输入信号；h(n)为传递函数；e(n)为随机噪声。

根据熵的概念，熵是表征时间序列复杂程度的物理量，当时间序列越复杂，熵值越大；当时间序列规律性越强时熵值越小。因此当包含轴承故障信息的周期冲击脉冲序列增多时，熵值减小[19]。MED的核心就是寻找一个L阶最优滤波器f，使输出信号y(n)最大程度恢复到输入信号x(n)，并使其熵值最小[20]，即

(7)

对两边求导可得

(8)

式中，L为逆滤波器长度。

以解卷积后得到的x(n)的范数作为目标函数去衡量x(n)熵值的大小。

(9)

MED的目的是求目标函数最大时的最优逆滤波器，对式(9)求偏导并让其等于零。

(10)

联立式(8)和式(10)可得

(11)

将式(11)写成矩阵的形式，即

B=A*F

(12)

进行迭代计算得到的逆滤波器的矩阵g为

g=A-1f

(13)

式中：B为输入与输出的互相关矩阵；A为信号y(n)的L×L自相关矩阵。

2 故障诊断方法及步骤

传统EWT通过尺度空间法对频谱进行自适应划分，得到初始分界点。但此时得到的分界点数量较多，将频谱分成过多的频带，对后续分析带来不便。本文提出基于互信息值的初始分界点进行筛选的方法，将分量互信息值大于均值的相邻频带合并为一条频带，分量互信息值小于均值的相邻频带合并为一条频带，以此来减少频带数量。

互信息发展于信息论中熵的概念，定义为两个随机变量间不确定度的差值，能够很好地衡量两随机变量的相关程度，比相关系数更准确[21]。

MI(X,Y)=H(Y)-H(Y|X)

(14)

式中：H(Y)为Y的熵；H(Y|X)为已知X时Y的条件熵。当X与Y相关性越强时，条件熵值H(Y|X)越小，此时互信息值MI(X,Y)越大[22]。图1为轴承故障诊断流程图。改进经验小波变换频谱划分流程：

步骤1对原始振动信号进行FFT(Fast Fourier Transform)，求得振动信号频谱；

步骤2根据尺度空间法确定初始频带分界点，获得初始划分频谱；

步骤3计算根据初始分界点划分频谱得到的各分量互信息值，若该分量互信息大于均值，则其与相邻互信息同大于均值的分量合并，若其相邻分量互信息小于均值，则该分量保持独立；若该分量互信息值小于均值，则将其与相邻互信息同小于均值的分量合并，若其相邻分量互信息大于均值，则该分量保持独立，根据分量合并情况重新确定分界点；

步骤4根据重新确定的分界点，对频谱进行重新划分，从而得到新的分解信号分量；

步骤5最优分量选择——互信息值最大的分量中包含较多与故障信息无关的振荡与噪声成分，因此在选择最优分量时，首先排除互信息最大的分量，当轴承发生故障时峭度值大于3，因此在剩余峭度值大于3的分量中，选择峭度大于均值的分量进行信号重构。

图1 轴承故障诊断流程图Fig.1 Fault diagnosis flow chart

从图1可知，本文提出轴承故障诊断方法的思路，即：首先，通过改进的EWT获得振动信号的分解模态分量，然后根据分量的互信息与峭度优选出最佳分量进行信号重构，为了提高诊断效果，对重构信号进行MED滤波处理，最后对滤波信号进行包络分析，可以明确诊断出轴承是否存在故障。

3 仿真信号分析

为证明本文方法的有效性，由冲击、噪声及谐波三部分构建滚动轴承单点损伤故障模型仿真信号进行验证，如式(15)～式(18)。

冲击信号x1

x1=e-αtAsin(2πf1t)

(15)

噪声信号x2

(16)

谐波信号x3

x3=Csin(2πf2t)

(17)

仿真信号y

y=x1+x2+x3

(18)

式中：A为冲击幅值，取值1.2 m/s2；B为噪声幅值，取值3.0 m/s2；C为谐波幅值，取值1.5 m/s2；α为衰减率，取值800；f1为共振频率，取值3 kHz；f2为谐波频率，取值25 Hz；fs为采样频率，取值10 240 Hz；z为随机数； 采样时长为1 s；fo为故障特征频率，取值64 Hz。图2为仿真信号各组成成分。图3为仿真信号时域图及包络谱。

图2 仿真信号各组成成分Fig.2 Components of simulation signal

图3 仿真信号时域图及包络谱Fig.3 Time-domain diagram and envelope spectrum of simulation signal

从图3(a)能观察到明显的冲击，但图3(b)所示的包络谱，故障特征频率及其倍频成分很不明显，难以确诊轴承故障。利用EWT尺度空间法对频谱进行划分，得到初始频带如图4所示。

图4 初始频带Fig.4 The initial frequency band

从图4可知，由于分界点较多，频谱被划分成53条频带。若采用传统方法，从EWT分解后的53个分量中选择表1中峭度较大的若干分量进行重构。重构信号如图5所示。

图5 重构信号、降噪信号时域图及包络谱Fig.5 Reconstructed signal time domain diagram,denoised signal time domain diagram and envelope spectrum

对图5(a)重构信号进行MED降噪，结果如图5(b)所示，噪声虽然得到削弱，冲击成分变得更加明显，但从图5(c)所示的包络谱中仍然不能确诊轴承故障。因此，根据分量互信息值重新对频带进行划分，如图6所示。

表1 分量峭度值

从图6(a)可知，只有分量1的互信息值大于均值，因此将其余小于均值的频带进行合并，如图6(b)所示。频带数由53减少为2，可知经过频带重划分解决了传统EWT频带划分过多的问题。重构分量见图6(c)，根据表2中的分量峭度与互信息，选择E2分量进行下一步分析。图7为对E2进行MED降噪得到降噪后时域图及包络谱。

图6 分量互信息、重划分频谱及分量时域图Fig.6 Mutual information of component,repartition of frequency spectrum and time domain of reconstructed component

表2 分量互信息与峭度值

从图7(a)可知，分量E2经MED降噪后周期冲击得到明显增强。从图7(b)的包络谱可以明显观察到轴承故障特征频率的倍频f0～7f0，由此可以确定该轴承存在故障。

图7 降噪后E2分量时域图及包络谱Fig.7 Time-domain diagram and envelope spectrum of E2 after denoise

仿真信号分析证明：轴承故障信号经传统EWT分解为较多的分量，根据峭度准则进行信号重构与MED降噪，降噪后的包络谱不能诊断出轴承故障。根据本文方法，在初始边界的基础上，根据互信息值对频带进行重新划分，有效减少频带数量，将谐波信息与故障冲击信息分离，从包络谱中可以诊断出明显的轴承故障，提高了轴承故障诊断的准确率。

4 实验信号

为了证明本文方法在实际应用中的有效性，采用货车轮对轴承实验台进行验证。图8所示为实验台结构图，实验选用CA-YD-188型加速度传感器，采用磁吸座的方式布置在靠近测试轴承的机架上，实验采用197726型双列圆锥滚子轴承，轴承参数如表3所示。

图8 铁路货车轴承轮对实验台Fig.8 Railway freight car wheelset running test bed

表3 双列圆锥型滚动轴承参数

4.1 内圈故障

采用含内圈剥离故障的轮对轴承进行实验，轮对转速为465 r/min，采样频率fs为12.8 kHz，采样时长为1 s。根据理论计算，轴承内圈故障特征频率fi为88.25 Hz，转频fr为7.75 Hz。内圈故障信号时域图及包络谱如图9所示。

从图9(a)能够观察到时域冲击成分，但从图9(b)中仅能观察到内圈故障频率的1倍频，无法判定轴承故障。图10为尺度空间法得到的初始频带。

图9 内圈信号时域图及包络谱Fig.9 Time-domain diagram signal and envelope spectrum of inner race

图10 初始频带Fig.10 The initial frequency band

从图10可知，EWT将信号频谱划分成47个频带，表4为选择的峭度值较大的分量进行信号重构，重构信号如图11所示。

表4 分量峭度值

图11 重构信号、降噪信号时域图及包络谱Fig.11 Reconstructed signal time domain diagram, denoised reconstructed signal time domain diagram and envelope spectrum

对图11(a)重构信号进行MED降噪，得到图11(b)降噪后时域图，从图11(c)中仅能勉强找到内圈故障频率的1倍频f0及转频fr，无法确诊轴承存在内圈故障。图12为根据本文方法对EWT频带进行重新划分。

图12 分量互信息及重划分频谱Fig.12 Band MI and redistributed spectrum

根据图12(a)分量互信息得到图12(b)重划分频谱。频带数量由47减少为8。图13为重划分分量时域图。

图13 分量时域图Fig.13 Time-domain diagrams of components

通过对比表5中各分量的互信息与峭度值，发现分量E3互信息值与峭度值均较大，图14为单独对E3进行MED降噪及包络分析结果。

从图14(a)可知，降噪后故障信号中的噪声得到削减，冲击性增强，但从图14(b)仅能发现内圈故障特征频率的1倍频fi与2倍转频2fr，不能确诊轴承存在内圈故障。

表5 分量互信息与峭度值

图14 分量E3降噪后时域图及包络谱Fig.14 Time domain diagram and envelope spectrum after E3 noise reduction

上述分析证明：互信息值最大的分量与原信号具有最大的相似度，包含较多的谐波与噪声成分，即使峭度值较大也不利于故障的提取，因此在选择分量进行重构时，排除互信息值最大的分量，在剩余分量中选择峭度较大的分量进行重构。

根据最优分量选取原则，排除互信息值最大的E3，结合图15，在剩余峭度值大于3的分量中选择峭度值大于均值的E7，E8进行重构。重构信号及MED降噪结果如图16所示。

图15 分量信号峭度值Fig.15 Kurtosis of component signals

从图16(b)可观察到明显周期冲击，从图16(c)可以观察到内圈故障的多个倍频fi～5fi及以转频fr=7.8 Hz为间隔的边频，据此可判定轴承存在内圈故障。

图16 重构信号、降噪信号时域图及包络谱Fig.16 Reconstructed signal, time-domain diagram and envelope spectrum of denoise signal

内圈故障实验证明：与传统方法对比，本文提出改进频谱重划分方法，能够有效减少频带数量，方便后续处理；该方法能够将包含较多谐波与噪声的互信息值最大分量与其他分量分离，在信号重构时，排除互信息最大的分量，在互信息较小的分量中进行选择，能够优选出包含故障信息的最佳分量。

4.2 外圈故障

采用含外圈剥离故障的货车轮对轴承进行实验验证。采样频率fs为25.6 kHz，转速465 r/min，时长为0.5 s。根据理论计算，轴承外圈故障特征频率为66.75 Hz。图17为外圈信号时域图及包络谱。

图17(a)可观察到故障冲击被淹没在噪声中，图17(b)所示包络谱中找不到外圈故障特征频率及其倍频。证明仅靠时域图和包络谱不能诊断轴承故障。

图17 外圈信号时域图及包络谱Fig.17 Time domain diagram and envelope spectrum of outer race

图18为尺度空间法得到的40个初始频带。在初始分界点基础上，根据互信息对频带进行重新划分，如图19所示。

图18 初始频带Fig.18 The initial frequency band

图19 基于分量互信息重新划分频带Fig.19 Band MI and repartited band

根据图19(a)将互信息值大于或者小于均值的相邻频带进行合并，得到图19(b)重划分频谱，从左至右依次为E1～E7，频带数量由40减少为7。图20为重划分频谱时域图。

图20 分量时域图Fig.20 Time domain diagrams of components

从图20可知，分量E5，E6冲击最明显，依据最优分量选取原则，在表6中排除互信息值最大的E2，根据图21从剩余峭度值大于3的分量中选择峭度大于均值的E5，E6进行重构。图22为重构信号及MED滤波后时域图及包络谱。

表6 分量互信息与峭度值

图21 分量信号峭度值Fig.21 Kurtosis of component signals

图22 重构信号、MED降噪后信号时域图及包络谱Fig.22 Reconstruction signal、time-domain diagram and envelope spectrum after MED denoise

从图22(a)可知，虽然重构信号周期冲击较明显，但仍包含较多噪声。对其进行MED降噪得到图22(b)，降噪后信号冲击性得到增强。从图22(c)可知，外圈故障特征频率的f0～6f0倍频，可以判断轴承存在外圈故障。

为了证明本文方法的优越性，与EEMD方法进行对比。采用EEMD对外圈故障信号进行分解，得到图23的IMF分量时域图(考虑篇幅原因，仅展示前6个IMF分量)。

图23 IMF分量时域图Fig.23 Time domain diagrams of IMFs

根据表7选择峭度最大的IMF2，IMF3与IMF5分量进行重构。图24为重构信号、降噪重构信号时域图及包络谱。

表7 IMF分量峭度值

图24 重构信号、降噪后时域图及包络谱Fig.24 Reconstruction signal、time-domain diagram and envelope spectrum after denoise

对图24(a)重构信号进行MED降噪，图24(b)为降噪后时域图，图24(c)包络谱中不能找到外圈故障特征频率及其倍频，不能确诊轴承故障。

外圈故障实验证明：本文方法能够有效减少EWT分量数量，并能将故障冲击成分与谐波和噪声成分分离，准确诊断出轴承外圈故障。通过与EEMD方法对比，证明本文提出的改进方法具有更高的诊断效率与准确性。

4.3 滚动体故障

采用含滚子剥离故障的货车轮对轴承进行实验，轮对转速为465 r/min，采样频率fs为25.6 kHz，采样时长为0.5 s。根据理论计算，轴承内圈故障特征频率fb为27.08 Hz。滚动体信号时域图及包络谱如图25所示。

图25 滚动体信号时域图及包络谱Fig.25 Time domain diagram and envelope spectrum of ball signal

从图25可知，仅通过频谱分析难以识别出滚动体故障，考虑到篇幅原因，仅对滚动体信号处理结果进行展示。图26为采用本文方法对滚动体故障轴承处理得到的包络谱。

图26 本文方法结果Fig.26 The result of proposed method

从图26可知，明显的滚动体特征频率fb=26.56 Hz及其多个倍频，可以确诊滚动体存在故障。同时证明本文提出方法对滚动体同样适用。

5 结 论

针对传统EWT划分频谱时得到频带过多，为后续最优分量选择带来困难的问题，本文在尺度空间法得到原始分界点的基础上，提出基于互信息的改进频带划分方法，有效减少了频带数量。得到的重构信号中往往包含噪声干扰，使用MED对重构信号进行降噪，提高了信噪比，增强了信号的周期性和冲击性。本文提出方法与EEMD方法对比，取得了更好的诊断效果，有较大的优势，并将其应用于货车轮对轴承的故障诊断。通过改进经验小波变换的频谱划分方法，避免了方法本身的缺点，为EWT的频谱划分问题研究提供了一个新的思路，提高了其可推广应用的范围。