相关系数峭度提取行星齿轮箱故障特征【机械制造与检测技术】
相关系数峭度提取行星齿轮箱故障特征摘要：针对行星齿轮箱振动信号源多，信号间耦合调制强烈，信号成分复杂，信噪比低导致其故障特征提取困难等问题，提出了一种新的基于VMD相关系数峭度的行星齿轮箱故障特征提取方法。即通过行星齿轮箱齿轮故障模拟试验采集振动信号，并对其进行变分模态分解(VMD)，得到本征模态函数分量(IMF)，计算IMF分量的相关系数，选择最优IMF分量作为特征分量，提取其峭度作为齿轮故障特征。实例分析结果表明，采用本文提出的方法较EEMD相关系数峭度方法、VMD相关系数样本熵方法更能有效地提取行星齿轮箱齿轮故障特征。
行星齿轮传动系统具有质量轻、体积小、传动比大、承载能力强、传动效率高等诸多优点，因此被广泛应用于作战飞机、舰船、装甲车辆、自行火炮、工程机械等大型复杂机械装备中[1--3]。由于经常承受大量级的交变载荷，导致其太阳轮、行星轮、内齿圈、行星架等关键部件出现严重磨损、疲劳裂纹、断齿和其他故障，尽早发现故障征兆并快速定位隔离，能有效避免故障恶化引发二次故障，降低维修成本和难度，提高设备可靠性和安全性。
从故障机理、检测与诊断方法等方面考虑，振动信号常被用来诊断行星齿轮传动系统故障，行星齿轮箱传动机理复杂，测得的振动信号调制强烈，频率成分复杂，行星轮和太阳轮故障难于区分，而现在很多研究还是套用定轴齿轮箱故障诊断方法，不能针对行星齿轮箱特点进行有效的分析[4]。某型行星齿轮箱为复合行星排行星齿轮箱，齿轮数量多，自由度维数高，同时存在多个定轴轮系和行星轮系，结构复杂，传感器测点选择十分困难，这些原因给该型行星齿轮箱故障诊断工作带来更大的困难。
Dragomiretskiy[5]等在2014年提出一种全新的非递归、自适应的信号处理方法变分模态分解(variational mode decomposition，VMD)，通过构造、迭代搜索约束变分模型的最优解以确定每个模态分量的中心频率和带宽，实现信号频谱的划分及各模态分量的有效分解[6]。不同于经验模态分解(empirical mode decomposition，EMD)、集合经验模态分解(ensemble empirical mode decomposition，EEMD)等传统递归算法，VMD摆脱了传统分解算法对信号分解的约束，有着坚实的数学基础，可以抑制或避免模态混叠现象，并且具有更高的运算效率及更好的鲁棒性[7-12]。采用VMD进行行星齿轮箱振动信号分解后，如何选取特征分量是特征提取前需要解决的问题[13-17]。选取的特征分量既要包含故障特征频率，又不能跟原始信号太相近，否则失去分解的意义，通过相关系数计算特征分量与原始信号的相关性来选择特征分量是目前常用的方法。行星齿轮箱齿轮产生磨损、疲劳裂纹、断齿等故障时，信号中会有周期性冲击成分，峭度与旋转部件的转速、尺寸和载荷无关，特别适用于旋转机械损伤类故障特征提取。
本研究针对传统以峭度为特征参数不能有效提取行星齿轮箱故障特征的问题，提出基于VMD相关系数峭度的故障特征提取新方法，并将其应用于某型行星齿轮箱。首先将信号进行VMD分解，得到若干IMF分量，然后依据相关系数理论选取与原始信号相关性好的分量作为特征分量，提取特征分量的峭度特征进行故障特征提取。新的组合方法克服了仅仅利用峭度进行故障特征提取噪声大、干扰多、效果不明显的问题，试验结果表明，该方法能有效提取故障特征。
相关系数峭度的基本原理行星齿轮箱振动信号经VMD分解后，得到K个中心频率为ωi的IMF分量，分解的实质就是变分问题构建与求解。
1) 变分问题构建[18]
1) 将实际输入信号S分解成K个IMF(t)分量，对其进行Hilbert变换，得到其解析函数和单边谱。
2) 在解析函数中加入预估指数e-jwkt，将IMF(t)频移到其中心频率。
3) 计算频移后IMF(t)的L2范数，解调估计其带宽。VMD变分问题可表示为：
(1)
式中: IMFk为所得模态分量；S为原始信号。
2) 变分问题求解
为求解该式(1)，引入增广拉格朗日函数，如式(2)所示。
(2)
式中：α为二次惩罚因子；λ为拉格朗日乘子；{IMFk}表示分解所得的K个IMF分量的集合；{ωk}={ω1，…，ωk}为中心频率集合；σ(t)为脉冲函数。采用乘法算子交替方向法求解式(2)，具体步骤[18]如下：
1) 初始化
2) 迭代次数n=n+1，循环开始。
3) 执行第1次内循环，并依据下式更新模态IMFk。
(3)
4) k=k+1，重复步骤3)，直到k=n，第一次内循环结束。
5) 执行第二次内循环，根据下式更新中心频率ωk。
(4)
6) k=k+1，重复步骤5)直到k=n，第二次内循环结束。
7) 根据式(5)更新拉格朗日乘子：
(5)
重复步骤2)至步骤7)，直到满足：
(6)
循环结束，得到K个IMF分量。
相关系数是衡量变量间相关程度的统计量，两信号相关系数越大，相关性越强，反之越弱。两信号x(t)、y(t)相关系数rxy计算为[19]：
(7)
峭度用来描述IMF分量的尖峰度，只与IMF分量中冲击信号有关，与复合行星齿轮尺寸、旋转速度以及载荷分布等无关，其计算方法如下：
(8)
式中： N为信号点数；σ为标准差；
为期望值。
行星齿轮箱齿轮故障模拟实验台及其结构框图如图1、图2所示，主要由驱动电机、转速转矩仪、传动箱、离合器、行星齿轮箱、负载电机、液压站等组成。行星齿轮箱行星传动部分由K1、K2、K3行星排组成，除了空档外，齿轮箱处于不同变速比时，均有≥2个行星排处于承载状态。本文以K3行星排太阳轮断齿故障模拟实验为例，由于故障设置在K3的某太阳轮轮齿面上，考虑到振动信号传递过程中存在衰减，因此振动传感器尽量安装在离振源最近的地方，如图2所示。为使实验更符合实际工况及负载情况，实验中行星变速箱输入转速设置为常用转速1 200 r/min，负载设置为接近实际负载的900 N·m。变速箱档位设置为1档，行星部分由行星排K2和行星排K3承载，行星排K1空载，K3行星排太阳轮固定。
图1 行星齿轮箱故障模拟实验台
图2 行星齿轮箱故障模拟实验台结构框图
图3 行星齿轮箱故障模拟实验设置(a)和K3太阳轮断齿(b)
依据齿轮相关参数和传动比理论计算可得行星排参数如表1所示。由表1可知，行星架转1圈历时t=1/3.245 436≈0.308 1250 s，采样频率20 kHz，为保证信号的周期性和可重复性，选取10个周期的信号长度进行分析，即3.081 25 s，即61 625点。
表1 K3行星排参数
参数数值太阳轮齿数30行星轮齿数15内齿圈齿数60行星轮个数6行星架转频/Hz3.245 436
对于太阳轮，如果发生局部故障，则其轮齿在旋转周期内与所有行星轮啮合，产生冲击，太阳轮局部故障特征频率为[20]：
(9)
式中： zs为太阳轮齿数； N为行星轮数量； fm为啮合频率。通过计算可知太阳轮局部故障特征频率为40 Hz。
实测正常和太阳轮断齿信号时域波形和频谱图如图4所示，从时域上看，故障信号无明显冲击脉冲，且周期不明显；从频域上看，实测信号频带成分复杂，干扰较多，找不到故障特征频率及其倍频，仅从时域和频域都不能区分正常和故障信号。
图4 实测信号时域及频谱图
对振动信号进行VMD分解过程中，模态数K的选择决定了信号的分解精度。K值过低VMD分解不充分，K值过高会导致信号过分解。本文通过中心频率观察法来判断是否发生过分解，进而确定K的取值。当K取不同值时，取10组太阳轮断齿故障信号(每组采样点数为61 625点)进行VMD分解，得到IMF分量的中心频率平均值如表2所示，二次惩罚因子α取默认值为2 000。
由表2可知，当K=5时，IMF5取到中心频率最小值，并且当K≥5时，IMF分量的中心频率最小值趋于稳定；当K=8时，IMF8取到中心频率最大值，并且当K≥8时，IMF分量的中心频率最大值趋于稳定，所以预设模态数K=8。
表2 取不同K值时VMD分解后各IMF分量中心频率
KIMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF8IMF911 08021 0602 26031 0602 2656 33141 0502 3203 5906 46056482 1603 3404 4306 68065802 1392 5403 6006 1207 24075782 1602 5203 5604 4306 3317 42085401 0602 2003 4704 3406 1417 0608 10095601 1602 3203 4354 3015 0406 2807 1608 166
分别取K3行星排正常和太阳轮故障信号各一组(采样点数为61 625)，VMD分解K取8，α取默认值2 000，将采集到的太阳轮故障信号和正常信号进行VMD分解，求解VMD分解后的每个IMF分量的相关系数，得到结果如表3所示。
表3 K3行星排正常和太阳轮故障IMF分量U3～U5的相关系数
IMF分量U1U2U3U4U5U6U7U8相关系数(正常)0.175 50.411 80.631 60.480 50.429 00.214 20.191 20.141 1相关系数(故障)0.203 60.417 90.671 00.488 70.351 50.242 60.227 20.189 5
由表3可知，K3行星排太阳轮正常和故障时，U3分量的相关系数均为最大，正常时，其次为U4和U5分量，为便于故障与正常比较，本文选取U3、U4、U5分量作为特征分量，计算其峭度，得到结果如表4所示。
表4 K3行星排正常和太阳轮故障IMF分量U3～U5的峭度
IMF分量U3U4U5峭度值(正常)2.216 52.897 02.510 7峭度值(故障)3.063 43.123 93.524 9
由表4中各IMF分量的峭度可知，K3行星排太阳轮故障时，U3、U4、U5分量的峭度值明显高于太阳轮正常时的峭度值。由峭度理论可知，U3、U4、U5分量中有较强的冲击信号成分，提取U3、U4、U5分量的峭度可以有效表征K3太阳轮断齿故障。
为验证VMD相关系数峭度方法的有效性，本研究首先通过与EEMD相关系数峭度方法对比验证VMD相关系数峭度方法中VMD分解的有效性，然后通过与VMD相关系数样本熵方法对比验证VMD相关系数峭度方法中峭度的有效性，从而验证本研究提出的基于VMD相关系数峭度的复合行星齿轮箱故障特征提取的有效性。
3.3.1 与EEMD相关系数峭度方法对比
基于EEMD相关系数峭度的故障特征提取方法是将信号进行EEMD分解，通过相关系数选取IMF分量作为特征分量，提取其峭度特征。复合行星齿轮箱振动信号经EEMD分解，得到的K3行星排齿轮正常和太阳轮故障时的IMF分量与原始信号的相关系数如表5所示。
表5 IMF分量的相关系数
IMF分量U1U2U3U4U5U6U7U8相关系数(正常)0.730 80.742 60.342 70.101 30.051 70.016 20.030 60.038 8相关系数(故障)0.767 60.717 20.364 20.118 70.041 20.021 20.026 20.021 3
由表5可知，K3行星排正常和太阳轮故障信号的IMF分量中，U1、U2、U3分量相关系数明显高于其他分量的相关系数，因此选取U1、U2、U3分量作为特征分量，计算其峭度，得到的结果如表6所示。
表6 IMF分量的峭度
IMF分量U1U2U3相关系数(正常)2.797 62.598 12.998 5相关系数(故障)3.821 93.372 63.000 6
由表6可知，K3行星排太阳轮故障时，U1、U2的峭度明显高于K3行星排正常时的U1、U2峭度特征。但是U3分量峭度特征则只是略高于与K3行星排正常时的峭度特征。再选取正常信号、故障信号各9组，分别计算经VMD分解和EEMD分解后K3太阳轮故障时IMF特征分量的峭度特征与K3正常时峭度特征的欧式距离[21]，得到10组故障信号、10组正常信号经2种方法分析得到的10组峭度欧式距离如图5所示。
图5 2种方法提取的峭度特征的距离
由图5可知，经VMD相关系数峭度方法计算得到的峭度欧式距离明显要高于经EEMD相关系数峭度方法计算得到的峭度欧式距离。即VMD相关系数峭度方法得到的K3正常和K3太阳轮故障时的特征差异更明显，间接验证了VMD相关系数峭度方法中VMD分解的有效性。
3.3.2 与VMD相关系数样本熵方法的对比
基于VMD相关系数样本熵的故障特征提取方法[22]是将信号进行VMD分解，通过相关系数选取IMF分量作为特征分量，提取其样本熵特征。基于3.2节VMD分解后信号，提取IMF特征分量的样本熵，得到的样本熵结果如表7所示。由表7可知，K3太阳轮故障时U3、U4要低于K3行星轮正常时的样本熵，K3太阳轮故障时U5样本熵明显高于K3行星轮正常时U5的样本熵。计算3.3.1节另外9组组故障信号、9组正常信号U3、U4、U5分量的样本熵，也得到了同样的结果。10组样本熵欧式距离与峭度欧式距离如图6所示。
表7 IMF分量的样本熵
IMF分量U3U4U5样本熵(正常)0.174 10.290 30.032 1样本熵(故障)0.162 10.230 50.137 2
图6 样本熵欧式距离与峭度欧式距离曲线
由图6可知，经VMD分解后的IMF特征分量的峭度欧式距离明显优于IMF特征分量的样本熵欧式距离，即正常和故障信号的峭度特征差异更明显。即VMD相关系数峭度方法中提取峭度特征有效。此外，计算过程中，由于样本熵重构的样本空间大，在内存16G的Windows 10 64位操作系统上计算耗时长达几分钟，相比峭度计算1 s内出结果，样本熵的计算效率太低。
将VMD、相关系数和峭度理论相结合，提出了基于VMD相关系数峭度的复合行星齿轮故障特征提取方法，并通过与EEMD相关系数峭度方法、VMD相关系数样本熵方法对比，验证了VMD相关系数峭度组合方法的有效性。相比单行星排齿轮箱，复合行星齿轮箱振动信号源更多，信号源间耦合调制更强烈，信号成分更复杂，信噪比更低，融合各种诊断理论与方法的优点。相应的组合方法可有效提取复杂信号的故障特征，也是故障诊断领域重要的研究内容之一。
[1] LI M，XIE L Y，DING L J.Load sharing analysis and reliability prediction for planetary gear train of helicopter[J].Mechanism and Machine Theory，2017，115：97-113.
[2] 周建星，刘更，马尚君.船用两级行星减速器箱体振动噪声分析[J].船舶力学，2014，18(1/2)：201-208.
[3] 孙灿飞，王友仁．直升机行星传动轮系故障诊断研究进展[J]．航空学报，2017，38(7)：106-119.
[4] 雷亚国，何正嘉，林京，等.行星齿轮故障诊断技术的研究进展[J].机械工程学报，2011，47[19]:59-67.
[5] DRAGOMIRETSKIY K，ZOSSO D.Variational mode decomposition[J].IEEE Transactions on Signal Processing，2014，62(3):531-54.
[6] ZHENG J D，PAN H Y，YANG S B，et al.Generalized variational mode decomposition and its applications to gearbox fault diagnosis under variable conditions[J].Journal of Vibration Engineering，2017，30(3):502-509.
[7] 钱林，康敏，傅秀清，等.基于VMD的自适应形态学在轴承故障诊断中的应用[J].振动与冲击，2017，36(3)：227-233.
[8] MIAO Y H，ZHAO M，LIN J.Identification of mechanical compound-fault based on the improved parameter-adaptive variational mode decomposition[J].ISA Transactions，2019，84：82-95.
[9] JIN S，KIM J S，LEE S K.Sensitive method for detecting tooth faults in gearboxes based on wavelet denoising and empirical mode decomposition[J].Journal of Mechanical Science & Technology，2015，29(8)：3165-3173.
[10] SMITH J S.The local mean decomposition and its application to EEG perception data[J].Journal of the Royal Society Interface，2005，2 (5)：443-454.
[11] FENG Z P，ZUO M J，QU J，et al.Joint amplitude and frequency demodulation analysis based on local mean decomposition for fault diagnosis of planetary gearboxes[J].Mechanical Systems & Signal Processing，2013，40(1):56-75.
[12] RAM R，MOHANTY M N.Comparative analysis of EMD and VMD algorithm in speech enhancement[J].International Journal of Natural Computing Research，2017，6(1):17-35.
[13] WANG Z J，WANG J Y，DU W H，et al.Application of parameter optimized variational mode decomposition method in fault diagnosis of gearbox[J].IEEE Access，2019(7)：44871-44882.
[14] WANG Z J，WANG J Y，DU W H.Research on fault diagnosis of gearbox with improved variational mode decomposition[J].Sensors，2018，18(10)：3510-3524.
[15] WANG C G，LI H K，HUANG G J，et al.Early Fault Diagnosis for Planetary Gearbox Base on Adaptive Parameter Optimized VMD and Singular Kurtosis Difference Spectrum[J].IEEE Access，2019(7)：31501-31516.
[16] ZHANG X，MIAO Q，ZHANG H，et al.A parameter-adaptive VMD method based on grasshopper optimization algorithm to analyze vibration signals from rotating machinery[J].Mechanical Systems and Signal Processing，2018，108:58-72.
[17] WANG Z J，WANG J Y，DU W H.Research on fault diagnosis of gearbox with improved variational mode decomposition[J].Sensors，2018，18(10)：3510-3524.
[18] 苏航.基于VMD的风机传动系统故障诊断研究[D].北京:华北电力大学，2016.
[19] 郭伟杰.基于电机电流的转子系统典型故障识别方法[D].太原：太原理工大学，2017.
[20] 门兰城.基于电机电流的行星齿轮箱故障诊断方法研究[D].太原：太原理工大学，2019.
[21] DAVID M J.Tax，Robert Duin，Dick De Ridder.An Engineering Approach Using MATLAB[J].Parameter Estimation and State Estimation，2014，0-470-09013-8.
[22] 杨大为，冯辅周，赵永东，等.VMD样本熵特征提取方法及其在行星变速箱故障诊断中的应用[J].振动与冲击，2018，37(16)：198-205.
本文引用格式：李玉豪，刘远宏.基于VMD相关系数峭度提取行星齿轮箱故障特征[J].兵器装备工程学报，2021,42(06):262-267.
Citation format:LI Yuhao， LIU Yuanhong.Fault Feature Extraction of Planetary Gearbox Based on VMD Correlation Coefficient Kurtosis[J].Journal of Ordnance Equipment Engineering,2021,42(06):262-267.