最大相关峭度解卷积结合稀疏编码收缩的齿轮微弱故障特征提取

唐贵基， 王晓龙

(华北电力大学能源动力与机械工程学院， 河北 保定 071003)

最大相关峭度解卷积结合稀疏编码收缩的齿轮微弱故障特征提取

唐贵基， 王晓龙

(华北电力大学能源动力与机械工程学院， 河北 保定 071003)

针对强背景噪声环境下齿轮早期故障诊断问题，提出了最大相关峭度解卷积(Maximum Correlated Kurtosis Deconvolution, MCKD)结合稀疏编码收缩的微弱故障特征提取方法。由于最大相关峭度解卷积算法的处理结果同时受滤波器长度参数及解卷积周期参数的影响，为自适应地实现最佳的解卷积效果，利用粒子群算法优良的寻优特性，对最大相关峭度解卷积算法的最佳影响参数组合进行搜索。原故障信号经影响参数优化的最大相关峭度解卷积算法处理后，冲击特征会明显增强，为剔除剩余噪声，对所获解卷积信号做进一步稀疏编码收缩降噪处理，并通过分析降噪信号的包络谱来识别故障特征频率成分。实例分析结果验证了该方法的有效性和可靠性。

故障诊断； 齿轮； 微弱特征； 粒子群优化； 最大相关峭度解卷积

引 言

齿轮是机械设备中广泛应用的传动部件,但故障发生概率高、容易损坏，准确有效地提取齿轮早期故障微弱特征对于避免发生严重故障、 减少经济损失意义重大。在实际工程应用中，由于振动传输路径复杂多变、工作环境噪声干扰严重、多振动源的激励和响应互相耦合等诸多因素的影响，使得齿轮早期微弱故障特征提取相对困难, 这也一直是故障诊断领域的一个研究热点[1-2]。

针对该问题，不少学者进行了深入研究，文献[3]利用经验模态分解方法对精轧机齿轮故障信号进行自适应分解，并对得到的固有模态分量做Hilbert包络解调来获取每个分量的幅频信息，最终实现强噪声背景下齿轮早期裂纹故障的准确诊断；文献[4]利用小波包的高分辨率分解及重构能力，将齿轮微弱故障信号分解到不同频段，然后选择有效频段进行信号重构，顺利提取出微弱故障特征信息；文献[5]将共振稀疏分解与包络谱相结合，采用共振稀疏分解将微弱冲击从齿轮振动信号中分离出来，然后对冲击序列做包络分析，可对齿轮早期故障做出准确判断；文献[6]提出离散小波变换与神经网络相结合的锥齿轮早期故障智能识别方法，首先通过离散小波变换提取故障特征，继而利用决策树来筛选故障特征，最后通过神经网络实现故障分类；针对齿轮微弱故障特征提取困难的问题,文献[7]提出基于稀疏信号分解的多阶分数阶傅里叶变换自适应滤波方法，该方法能很好地从频率呈曲线变化的多分量信号中剥离出啮合频率包络调制信号，对滤波后的信号做进一步解调分析，可顺利提取出微弱故障特征，上述各方法在齿轮早期故障诊断应用中均取得了一定的效果。

最大相关峭度解卷积[8]可通过循环迭代实现解卷积运算，突出信号中包含的连续脉冲，但该算法的处理结果受滤波器长度及解卷积周期两个参数的影响，如果任何一个参数设置不合适，都难以达到理想的分析效果。为了实现参数的自适应筛选，本文提出利用粒子群算法对MCKD算法的最佳影响参数组合进行搜寻，并将影响参数优化的MCKD算法与稀疏编码收缩算法相结合，用于处理齿轮早期故障信号，提出了最大相关峭度解卷积结合稀疏编码收缩的微弱故障特征提取方法，实例分析结果表明该方法可实现齿轮早期故障的精确诊断。

1 基于粒子群优化的MCKD算法

1.1 MCKD算法

零均值信号yn(n=1,2,…,N)的相关峭度表达式为

(1)

相关峭度是在峭度基础上提出的概念，当解卷积周期参数T=0时，相关峭度即退化为峭度，然而峭度对少数异常突出脉冲过于敏感，相关峭度则充分考虑了冲击成分的连续性，与峭度相比能更准确地衡量信号中周期性脉冲序列所占的比重。

假设y为一个连续冲击信号，x为冲击信号y经周围环境及路径传输后实际采集到的响应信号，上述过程可通过如下公式表示为

(2)

为便于分析，在此先不考虑噪声e的影响，则MCKD算法本质是寻找一个FIR滤波器，通过输出信号x恢复输入信号y，即

(3)

为突出解卷积运算结果中的连续尖脉冲，该算法以信号的相关峭度为评定标准，并将相关峭度最大化作为目标函数，即

(4)

上述优化求解问题等价于求解方程

(5)

由公式(3)～(5)可归纳得到如下表达式

(6)

以矩阵的形式将上式重新表述为

(7)

整理式(7)可得

(8)

由于已知

(9)

则滤波器系数可通过下式得到

(10)

将滤波器系数带入式(3)，即可得到实际采集信号x的解卷积信号y(原冲击信号)。

1.2 最佳影响参数组合的自适应选取

MCKD算法可有效增强信号中隐含的冲击脉冲，适合处理信噪比低并具有周期性冲击特性的信号，然而滤波器长度参数L和解卷积周期参数T的选择对解卷积结果有着至关重要的影响。如果保持一个影响参数不变，仅以另一个影响参数为优化对象，讨论该参数对运算结果的影响，这种局部寻优方式由于忽略了两个参数间的交互作用，因此得到的也只是相对最优影响参数。粒子群算法[9]作为一种群体智能优化算法，具有良好的全局寻优能力，本文采用粒子群算法对MCKD算法的两个影响参数进行并行搜索优化，可实现滤波器长度参数L和解卷积周期参数T的自适应筛选。

假设在一个D维空间中，由M个粒子组成种群X=(X1,X2,…,XM)，其中第i个粒子表示一个D维向量Xi=(xi1,xi2,…,xiD)，代表第i个粒子在D维搜索空间中的位置(即优化问题的潜在解)。第i个粒子的速度为Vi=(vi1,vi2,…,viD)，其个体局部极值为Pi=(pi1,pi2,…,piD)，种群全局极值为G=(g1,g2,…,gD)，每个粒子通过个体局部极值和种群全局极值迭代更新自身的速度和位置，公式表述如下[10]

(11)

式中ω为惯性权重；d=1,2,…,D；i=1,2,…,M；k为当前迭代次数；c1和c2为加速度因子；η为介于[0,1]间的随机数。

利用粒子群算法搜寻影响参数时，需确定一个适应度函数，粒子每次更新位置时计算一次该位置对应的适应度值，通过适应度值对比进行更新。由于包络谱熵能有效衡量信号中的周期性冲击成分[11]，本文便将解卷积信号的包络谱熵作为寻优过程中的适应度函数，用于评定MCKD算法解卷积运算结果的优劣。零均值信号y(j)(j=1,2,…,N)的包络谱熵Ep可表示成如下形式

(12)

式中pj为a(j)的归一化形式，a(j)为信号y(j)包络谱的离散化形式。

以任意粒子Xi的位置为MCKD算法的滤波器长度参数和解卷积周期参数，设定好这两个影响参数后，利用MCKD算法对原信号进行处理并计算所获解卷积信号的包络谱熵，熵值Ep(Xi)即为该粒子位置的适应度值。当解卷积结果中周期性脉冲成分突出时，包络谱熵值较小，而解卷积效果欠佳、规律性冲击成分不明显时，包络谱熵值则相对较大，因此包络谱熵值最小化即为最终搜索目标。利用粒子群算法搜寻优化MCKD算法最佳影响参数组合的具体步骤如下：

(1)初始化粒子群算法的各项参数。

(2)初始化粒子种群，以影响参数组合[L,T]作为粒子的位置，每个位置对应一对影响参数，随机产生一定数量的参数组合作为粒子的初始位置，随机初始化每个粒子的移动速度。

(3)以粒子的位置为MCKD算法的两个影响参数，并对原信号做解卷积运算，计算每个粒子Xi的位置对应的适应度值Ep(Xi)。

(4)对比适应度值大小并更新个体局部极值和种群全局极值。

(5)利用公式(11)更新粒子的速度和位置。

(6)循环迭代, 转至步骤(3)，直至迭代次数达到最大设定值，输出最佳适应度值及粒子的位置。

2 稀疏编码收缩

稀疏编码收缩算法可从高斯信号中推测得到隐含的非高斯成分，Hyvarinen提出如下模型来表示非高斯信号的概率密度函数[12]

(13)

式中d为信号的标准差，α为控制概率密度函数稀疏性的参数，对于具有冲击特性的非高斯信号来说，参数α可设定为0.1[13]。

假设x为原始信号，其统计特征呈现出非高斯性质，v是均值为零的高斯白噪声信号，则观测信号可表示为

y=x+v

(14)

(15)

3 诊断流程

分析实际信号时，选择合适的影响参数,有效突出信号中隐含的规律性故障冲击，是利用MCKD算法处理齿轮早期故障信号的关键所在。本文利用粒子群算法自适应地筛选最佳的影响参数组合，可避免参数设定时人为主观盲目选择带来的弊端。齿轮故障信号经影响参数优化的MCKD算法处理后，信噪比会显著提高，原信号中被噪声所淹没的周期性脉冲能够有效突显出来，但仍然会存在较多噪声干扰。由于齿轮故障冲击成分通常具有非高斯性质，而噪声成分则呈现出高斯分布特性，稀疏编码收缩算法可实现高斯信号和非高斯信号的有效分离，因此利用该算法对所获的解卷积信号做进一步处理，保留冲击成分的同时可有效剔除剩余噪声。综合上述分析，本文提出最大相关峭度解卷积结合稀疏编码收缩的齿轮微弱故障特征提取方法，有望实现齿轮早期故障的准确识别，诊断流程如图1所示，具体实现过程如下：

图1 诊断流程图Fig.1 Diagnosis flow chart

(1)设定粒子群算法各项参数并搜寻MCKD算法的两个影响参数。在粒子群参数设置过程中，如果最大进化代数及种群规模这两个参数值设置的太大，则计算负担过重，如果设置的过小，则无法保证搜寻到最佳的影响参数组合，参考文献[10]中的参数设定值以及多次试验对比分析，本文设定的粒子群算法各项参数值如表1所示。

表1 粒子群算法各项参数

(2)粒子群寻优结束后会得到一组最佳影响参数组合[Lo,To]，设定MCKD算法的滤波器长度参数为Lo，解卷积周期参数为To，并利用影响参数优化的MCKD算法对故障信号做预处理。

(3)利用稀疏编码收缩算法对所得的解卷积信号做进一步降噪处理。

(4)对降噪信号做Hilbert包络解调运算，并计算包络谱。

(5)将齿轮故障特征频率的理论计算值与包络谱中的谱线进行对比，实现故障最终判别。

4 应用实例

4.1 点蚀故障

在齿轮箱故障试验台从动齿轮的齿面上加工出一个微小凹痕来模拟齿轮轻微点蚀故障，图2为试验平台的结构简图，在齿轮箱端盖靠近轴承处布置4只加速度传感器，设置的采样频率为5 120 Hz，采样点数为8 192点。其中，电机转速为850 r/min,输入轴主动轮齿数为55，输出轴从动轮齿数为75，通过计算可知输入、输出轴转频分别为fi=14.17 Hz，fo=10.39 Hz。

图2 试验平台Fig.2 Experimental platform

实测信号波形如图3所示，时域波形中隐约可见少量冲击,但毫无规律可寻。对实测信号直接做包络谱分析，结果如图4所示。没有发现幅值突出的频率成分，由此可见传统包络谱分析手段对于该故障信号来说显得无能为力。

图3 实测信号波形Fig.3 Waveform of measured signal

图4 实测信号包络谱Fig.4 Envelope spectrum of measured signal

图5 本文提出方法的实测信号分析结果Fig.5 Analysis results of measured signal by proposed method

下面利用本文提出方法对实测信号进行分析，结果如图5所示。图5(a)为包络谱熵随进化代数变化的关系曲线，包络谱熵最小值7.149 5出现在第9代进化种群中，粒子群算法寻优所得最佳影响参数组合为[56,486]，由此设定MCKD算法的滤波器长度参数为56，解卷积周期参数为486，利用影响参数优化的MCKD算法对实测信号做预处理，结果如图5(b)所示。通过分析解卷积信号的波形可发现，原本被强烈噪声所掩盖的冲击成分被有效突显出来。图5(c)是采用稀疏编码收缩算法对解卷积信号做进一步处理所获的降噪信号，解卷积信号经过降噪后，大量冗余干扰噪声被剔除，而与齿轮故障相关的冲击成分则被保留下来，图5(d)是降噪信号的包络谱分析结果，图中主要成分为输出轴转频及其倍频fo～9fo，由此判断输出轴上的从动齿轮存在局部损伤，分析结果与实际情况相符，实例分析表明本文提出方法可顺利识别出齿轮早期故障的微弱特征频率信息。

下面分别单独应用参数优化MCKD算法和稀疏编码收缩算法对实测信号进行处理来验证二者结合的必要性。图6为原信号经MCKD算法处理所得解卷积信号的包络谱(MCKD算法的滤波器长度及解卷积周期参数与上述分析相同，分别为56，486)，虽然谱图中存在多个与故障特征相关的频率成分，但由于存在较多干扰谱线，特征频率基频及5倍频成分处谱峰并不突出。图7为原信号经稀疏编码收缩算法处理所得降噪信号的包络谱，图中仅特征频率基频处谱线幅值偏高，但峰值仍不够明显。通过对比可发现，将MCKD算法和稀疏编码收缩算法结合起来分析故障信号，可有效剔除原信号中的强烈噪声干扰成分，最大限度突出微弱故障信息。

图6 解卷积信号包络谱Fig.6 Envelope spectrum of deconvolution signal

图7 降噪信号包络谱Fig.7 Envelope spectrum of denoising signal

为验证所述方法优势，利用文献[3]所用的经验模态分解(Empirical Mode Decomposition, EMD)方法和文献[4]所用的小波包分解方法分别对实测信号进行处理，并做包络解调分析，各取其最佳包络谱(所有分量均做包络谱，取效果最佳的分量)与本文方法进行对比。图8为实测信号EMD处理结果，自适应分解后共获得7个IMF分量，其中IMF1分量的包络解调效果最佳，图9为该分量的时域波形及包络谱。对实测信号做2层小波包分解，小波函数为文献[4]采用的db5小波，结果如图10所示，d20～d23是小波包分解所获4个节点系数的重构信号，其中d21节点重构信号的包络解调效果最佳，重构信号波形及包络谱如图11所示。对比分析后发现，图9(b)、图11(b)所示包络谱中，齿轮故障相关频率成分的幅值水平以及倍频成分的连续性与图5(d)相比，均存在较大差距，由此表明本文方法在齿轮微弱故障特征提取上具有一定优势。

图8 EMD处理的结果Fig.8 Processing results of EMD

图9 IMF1分量的波形及包络谱Fig.9 Waveform and envelope spectrum of IMF1

图10 小波包处理的结果Fig.10 Processing results of wavelet packet

图11 重构信号的波形及包络谱Fig.11 Waveform and envelope spectrum of reconstructed signal

4.2 磨损故障

试验数据来自美国学者J Rafiee的个人网站[14]，试验装置包括齿轮箱、电机、加载机构、脉动分析仪、三轴加速计、转速计及减震器。电机实测转速为1 397 r/min，输入轴主动轮齿数为29，输出轴从动轮齿数为24，输入、输出轴转频分别为fi=23.29 Hz，fo=28.14 Hz。加速度计位于齿轮箱外部靠近输入轴位置处，采样频率为16 384 Hz。试验过程采集了主动齿轮发生轻微磨损、中度磨损和断齿3种不同程度故障的振动信号，本文选择对轻微磨损阶段的实测信号进行分析，分析点数为16 384点。

图12所示实测信号波形中，基本没有出现明显冲击脉冲,对该信号做进一步包络谱分析，结果如图13所示。图中存在多个谱峰，除了在输入轴二倍转频2fi处找到一根幅值偏高谱线外，再没有发现其它故障相关频率成分，仅通过这一有限信息很难准确判断齿轮箱是否出现故障。

图12 实测信号波形Fig.12 Waveform of measured signal

图13 实测信号包络谱Fig.13 Envelope spectrum of measured signal

图14是利用本文方法对实测信号进行分析所得结果。图14(a)中，粒子种群仅通过3代进化就搜寻到了滤波器长度参数和解卷积周期参数的最佳组合[212,91]。通过观察图14(b)解卷积信号波形可发现，经影响参数优化的MCKD算法处理后，实测信号中隐含的规律性冲击被顺利剥离出来，对解卷积信号做进一步稀疏编码收缩降噪处理后，噪声成分被剔除殆尽，脉冲序列的周期特征变得更加明显。图14(d)所示降噪信号包络谱中，输入轴转频及其谐波fi～6fi处谱线峰值十分突出，基本不存在其他干扰成分，谱图相对干净。通过综合分析降噪信号的波形及包络谱可知，冲击脉冲出现的频率与输入轴转频相同，由此断定输入轴上的主动齿轮发生了故障。

图14 本文提出方法的实测信号分析结果Fig.14 Analysis results of measured signal by proposed method

同样，分别利用经验模态分解和db5小波包分解对实测信号进行处理(篇幅所限，所获各分量时域波形省略)，并取最佳分量的包络谱与本文方法进行对比。经验模态分解处理后共获得6个IMF分量，其中IMF3分量的包络谱效果最佳，图15是该分量的波形及包络谱。2层小波包分解重构后，发现d23节点重构信号的包络分析结果最佳，图16是重构信号的波形及包络谱。通过对比可知，图15(b)和图16(b)均没有图14(d)效果明显，齿轮轻微磨损故障信号对比分析结果再次验证了本文方法的优势。

图15 IMF3分量的波形及包络谱Fig.15 Waveform and envelope spectrum of IMF3

图16 重构信号的波形及包络谱Fig.16 Waveform and envelope spectrum of reconstructed signal

5 结 论

(1)齿轮早期故障特征信号微弱，并受环境噪声及传递路径衰减影响，特征提取相对困难，传统包络谱方法可能诊断失效，利用MCKD算法分析早期故障信号，可有效增强微弱冲击特征。

(2)滤波器长度参数和解卷积周期参数的选择对MCKD算法处理结果起着决定性作用，利用粒子群算法自适应地搜索、优化MCKD算法的影响参数是行之有效的。

(3)利用稀疏编码收缩算法对解卷积信号进行处理，可实现非高斯故障冲击信号与高斯噪声干扰信号的有效隔离，达到降噪目的。

(4)实例分析及对比结果表明，最大相关峭度解卷积结合稀疏编码收缩的微弱故障特征提取方法可实现齿轮早期故障诊断，并具有一定可靠性和优势。

[1] 梅检民，肖云魁，贾继德，等.基于改进阶比的变速器微弱故障特征提取[J].振动工程学报，2012，25(3)：317—322.

MEI Jianmin，XIAO Yunkui，JIA Jide，et al. Weak fault characteristics extraction for automobile transmission based on improved order analysis[J].Journal of Vibration Engineering，2012，25(3)：317—322.

[2] 雷亚国,韩冬,林京，等.自适应随机共振新方法及其在故障诊断过中的应用[J].机械工程学报，2012，48(7)：62—67.

LEI Yaguo，HAN Dong，LIN Jing，et al. New adaptive stochastic resonance method and its application to fault diagnosis[J].Journal of Mechanical Engineering，2012，48(7)：62—67.

[3] 高立新, 吴丽娟, 张建宇.基于EMD解调方法的齿轮早期故障诊断[J].北京工业大学学报，2009，35(7)：876—881.

GAO Lixin, WU Lijuan, ZHANG Jianyu. An application in gear incipient failure diagnosis based on EMD demodulation method[J].Journal of Beijing University of Technology，2009,35(7):876—881.

[4] 杨洁明，熊诗波.小波包分析方法在齿轮早期故障特征提取中的应用[J].振动、测试与诊断，2000，20(4)：269—272.

YANG Jieming，XIONG Shibo. Extraction of early fault information of gears using wavelet packet[J].Journal of Vibration, Measurement & Diagnosis，2000，20(4)：269—272.

[5] 张文义，于德介，陈向民.基于信号共振稀疏分解与包络谱的齿轮故障诊断[J].中国机械工程，2013，24(24)：3 349—3 354.

ZHANG Wenyi，YU Dejie，CHEN Xiangmin. Fault diagnosis of gears based on resonance-based sparse signal decomposition and envelope spectrum[J]. China Mechanical Engineering，2013，24(24)：3 349—3 354.

[6] SARAVANAN N, RAMACHANDRAN K I. Incipient gear box 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)：4 168—4 181.

[7] 梅检民，肖云魁，沈虹，等. 多阶FRFT自适应滤波及齿轮微弱故障特征提取[J].振动工程学报，2013，26(5)：771—778.

MEI Jianmin，XIAO Yunkui，SHEN Hong，et al. Multi-order FRFT adaptive filter and its application to extract weak fault feature based on sparse signal decomposition [J].Journal of Vibration Engineering，2013，26(5)：771—778.

[8] MCDONALD G L，ZHAO Q，ZUO M J. Maximum correlated kurtosis deconvolution and application on gear tooth chip fault detection[J].Mechanical Systems and Signal Processing，2012，33：237—255.

[9] KERNNEDY J，EBERHART R C. Particle swarm optimization[A].Proceeding of IEEE International Conference on Neutral Networks[C]． Perth, Australia，1995：1 942—1 948.

[10]沈伋,韩丽川,沈益斌.基于粒子群算法的飞机总体参数优化[J].航空学报，2008，29(6)：1 538—1 541.

SHEN Ji，HAN Lichuan，SHEN Yibin. Optimization of airplane primary parameters based on particle swarm algorithm[J].ACTA AERONAUTICA ET ASTRONAUTICA SINICA，2008，29(6)：1 538—1 541.

[11]何水龙，訾艳阳，万志国，等.自适应提升多小波在螺旋伞齿轮故障诊断中的应用[J].仪器仪表学报，2014，35(1)：148—153.

HE Shuilong，ZI Yanyang，WAN Zhiguo, et al. Application of adaptive multiwavelets via lifting scheme in bevel gear fault diagnosis[J].Chinese Journal of Scientific Instrument，2014，35(1)：148—153.

[12]HYVARINEN A. Sparse code shrinkage: de-noising of nonGaussian data by maximum likelihood estimation[J]. Neural Computation,1999,11：1 739—1 768.

[13]LIN J,ZUO M J, FYFE K R. Mechanical fault detection based on the wavelet de-noising technique[J].Journal of Vibration and Acoustics,2004，126：9—16.

[14]RAFIEE J. Free gearbox vibration signals[EB/OL]. http://rafiee.us/index.html,2008.

Weak feature extraction of gear fault based on maximum correlated kurtosis deconvolution and sparse code shrinkage

TANGGui-ji，WANGXiao-long

(School of Energy，Power and Mechanical Engineering， North China Electric Power University， Baoding 071003， China)

Aiming at incipient fault diagnosis problem of gears under strong background noises, a feature extraction method for weak faults based on maximum correlated kurtosis deconvolution and sparse code shrinkage was proposed. As the processing result of maximum correlated kurtosis deconvolution algorithm was affected by filter length parameter and deconvolution period parameter, in order to adaptively achieve the best deconvolution result, the particle swarm optimization algorithm with excellent optimization characteristic was applied to search for the optimal combination of influencing parameters of maximum correlated kurtosis deconvolution algorithm. The impact characteristic of the original fault signal could be enhanced after processed by maximum correlated kurtosis deconvolution algorithm with optimized parameters. In order to eliminate the residual noise, the deconvolution signal was further processed by sparse code shrinkage de-noising algorithm. Then the fault characteristic frequency components could be identified by analyzing the envelope spectrum of the de-noising signal. The analysis results verified the effectiveness and reliability of this method.

fault diagnosis; gear； weak feature； particle swarm optimization； maximum correlated kurtosis deconvolution

2014-03-29；

2014-09-23

河北省自然科学基金资助项目(E2014502052)；中央高校基本科研业务费专项资金资助项目(2015XS120)

TH165+.3; TH132.417

A

1004-4523(2015)03-0478-09

10.16385/j.cnki.issn.1004-4523.2015.03.019

唐贵基(1962—)，男，教授,博士生导师。电话：18603123318；E-mail：tanggjlk@ncepubd.edu.cn

王晓龙(1989—)，男，博士研究生。电话：13400430727；E-mail：wangxiaolong0312@126.com