基于参数优化变分模态分解和多尺度熵偏均值的行星变速箱故障特征提取

杨大为， 赵永东， 冯辅周， 江鹏程， 丁闯

(陆军装甲兵学院 机械工程系， 北京 100072)

0 引言

某型装甲车辆行星变速箱长期工作于高温重载的恶劣环境，齿轮故障常有发生，而其负载常在大范围瞬时波动，弱化了故障产生的异常进而掩盖了故障，很难及时发现并进行有效处理，故障往往进一步发展恶化，进而严重影响变速箱正常运转，造成车辆机动性能下降，减弱整车的战斗力。而某型行星变速箱结构复杂，含有多个定轴轮系和行星轮系，工作时多对齿轮啮合相互影响，信号分解难度较大。振动传感器采集到的信号存在大量噪声干扰且受复杂的传递路径影响衰减严重，属于典型非线性非平稳信号，给故障特征提取工作带来困难。行星变速箱的3个行星轮系含有3个或6个行星轮。单一行星轮故障对轮系影响本身就较太阳轮小，加之行星轮运动为复合运动，而振动传感器安装位置固定，不能随其旋转，相对位置改变会对信号进行二次调制，因此行星轮故障更难于区分[1]。目前，国内外研究人员针对行星变速箱故障诊断在信号分析处理方面做了大量工作，并取得了一定成果，如小波包分解、包络分析、倒频谱分析、经验模式分解等[2]。但对复合行星变速箱的故障诊断研究工作较少，且其故障模拟实验台搭建经验不足，实验条件仍不够成熟[3]。

变分模态分解(VMD)算法由Dragomiretskiy等[4]于2014年提出，能够将信号自适应地分解为一系列分量，摆脱了传统信号分解的递归筛分剥离模式约束。Ram等[5]验证了VMD算法可以缓解模态混叠和边界效应，同时具有强鲁棒性和高运算效率。Tang等[6]将VMD算法应用于滚动轴承复合故障诊断，并验证了VMD算法在分解强噪声信号方面的优越性。由于使用VMD算法时，需选定预设尺度数k和二次惩罚因子α，两个参数的取值都会对分解结果产生影响，故使用VMD算法的难点和关键在于如何选定合适的k值和α值。对此，姜万录等[7]采用中心频率观察法选取k值，凭借经验选取α值，处理信号取得一定效果。但如果分别考虑单一参数对分解的影响，忽略了两个参数间的相互作用，得到的结果并不能保证达到最优，并且实际信号复杂多变，仅凭借经验选取参数难以保证合理性。唐贵基等[8]引入粒子群优化(PSO)算法对VMD算法参数进行优化选取，避免人为因素干预，同时考虑两个参数的影响，优化结果更加合理。

样本熵由Richman在近似熵基础上提出，能计算时间序列复杂度，适合处理非线性非平稳信号，具有无需自我匹配、计算快、精度高的优点[9]。由于时间序列的复杂度与尺度有着密切关系，Costa等在样本熵的基础上进行多尺度分析，即多尺度熵[10]。Pan等[11]、Zheng等[12]将多尺度熵作为特征向量输入支持向量机，较单一的样本熵取得更好分类效果，但支持向量机建立知识库需大量样本，并只能做出对故障进行定性识别，不能定量化描述故障程度。张龙等[13]提出了多尺度熵偏均值，综合多尺度熵值和变化趋势全面反映信号的复杂度，能很好地衡量故障程度和追踪故障发展，区分不同故障效果很好，无需再结合分类器进一步分类，简化了故障诊断流程。

综合考虑VMD算法和样本熵的优越性，本文采用PSO算法和人工观察算法对VMD算法参数进行优化，基于互信息对分量进行筛选重构。在综合样本熵和多尺度熵优点的基础上，使用多尺度熵偏均值来综合衡量重构信号在不同尺度下的特征，更加完整和全面地反映信号的复杂程度，对行星变速箱的故障进行有效特征提取。

1 VMD算法介绍

1.1 VMD算法参数优化

作为一种新的非递归算法，VMD算法能自适应地将信号分解为k个中心频率为ωk的模态函数uk. VMD算法的实质为构造和求解变分问题的过程。

VMD算法构造的约束变分问题可表示为

(1)

(2)

式中：uk：={u1,u2,…,uk}为各模态函数；ωk：={ω1,ω2,…,ωk}为各模态函数的中心频率。

引入拉格朗日乘子λ(t)和二次惩罚因子α，将约束变分问题变为非约束变分问题，得到增广拉格朗日表达式：

L({uk},{ωk},λ):=

(3)

采用乘法算子交替方向法，通过迭代更新，求得扩展拉格朗日表达式的鞍点。算法具体过程较为复杂，在此不作赘述，具体过程参考文献[4-6].

VMD算法中，预设尺度数k决定分解所得分量个数，进而影响分解结果。二次惩罚因子α影响分解得到的各分量带宽，α越小，分解得到的各分量带宽越大，α越大则带宽越小。PSO算法能够根据某个适应度函数的最大值、最小值作为优化准则，来查找所需的最优参数组合。使用PSO算法搜寻VMD算法最优参数组合，需要确定一个适应度函数。

齿轮发生故障时，会产生以啮合频率及其倍频为载波，故障齿轮转频为调制频率的调制现象。齿轮故障特征提取常通过包络分析提取齿轮信号中的调制信息，通过识别由齿轮故障产生的特征频带对齿轮状态进行判别[14]。唐贵基等[8]提出了包络熵Ep，即将信号经Hilbert解调后所得的包络信号视为概率分布序列，并计算其熵值。已知具有N个数据点的时间信号x(j)，其Ep定义为

(4)

(5)

式中：pj为a(j)的归一化形式，j=1,2,…,N；a(j)为信号x(j)经Hilbert解调后得到的包络信号。

Ep可定量衡量原始信号的稀疏性，若齿轮振动信号含有明显由故障引起的调制现象和周期性冲击，则信号稀疏性越强，Ep值越小；反之，信号若含有大量噪声干扰，由故障产生的调制现象和周期性冲击越不明显，则信号稀疏性越弱，Ep值越大。

当VMD算法取某个[k,α]组合时，即粒子位于某一位置时，求得VMD算法的k个全部分量，计算它们的Ep值，k个Ep值中最小的一个称为局部极小熵值。选择局部Ep值作为整个搜索过程中的适应度值，即选取Ep作为适应度函数，将局部极小值作为搜索目标。

PSO算法优化流程如图1所示。

1.2 多尺度熵偏均值

样本熵可衡量行星变速箱振动信号复杂程度，系统不同的运行状态对应着不同的样本熵，可用于定量表征系统的运行状态。如果振动信号成分单一，周期性越明显，信号噪声干扰越少，信号复杂程度越低，样本熵值越小;反之，信号噪声干扰越多，信号复杂程度越高，样本熵值越大。样本熵具体计算过程可参考文献[9].

多尺度熵是在数据的多个时间尺度上分别计算样本熵值，不同尺度时间序列通过粗粒化过程获得[10-11]：

1)原始数据X={x(n),n=1,2,…,N}，数据长度为N，建立新的粗粒序列：

(6)

式中：τ为尺度因子，本文取τ=20.

2)分别计算每一个粗粒序列的样本熵，得到τ个粗粒序列的样本熵值，即多尺度熵。

同一信号在不同尺度因子下的样本熵值不同，并且不同信号的样本熵值随尺度因子变化趋势也不同。为了综合分析多尺度熵值，使用平均值观察多尺度熵的集中趋势，使用偏斜度来比较不同信号的多尺度熵偏斜程度[13]。在结合多尺度熵和偏斜度的基础上，多尺度熵偏均值的计算步骤如下：

1)计算时间序列的多尺度熵值，并将其记为MSE(τ)=[MSE(1),MSE(2),…,MSE(20)].

2)计算多尺度熵的偏斜度

Ske=3(MSEa-MSEm)/MSEs,

(7)

式中：MSEa、MSEm和MSEs为多尺度熵的平均值、中位数和标准差。偏斜度Ske范围为[-3,+3]，0为完全对称，+3为极右偏态，-3为极左偏态，绝对值越大，偏斜程度越大。

3)计算多尺度熵偏均值

PMMSE=(1+|Ske|/3)MSEa.

(8)

1.3 信号分析流程

结合VMD算法参数优化和多尺度熵偏均值的分析方法基本步骤如下：

1)采用基于Ep的PSO算法，把Ep作为适应度函数，搜寻局部Ep的最小值，从而优化VMD算法参数，并结合人工观察法最终确定参数，采用具有最优参数的VMD算法对振动信号进行分解。

2)计算分解所得的各分量与原信号的互信息并确定互信息阈值，选取大于阈值的有效分量进行叠加，从而完成对信号的重构。

3)计算重构信号的多尺度熵，进而计算多尺度熵偏均值，并将其作为故障特征参数对行星变速箱进行故障诊断。

2 实例分析

2.1 实验介绍

实验对象为某型三轴式离合器换挡3自由度行星变速箱，其主动轴齿轮和中间轴齿轮为定轴齿轮，主轴上有3个行星排，K1排为外啮合双行星排，其他2排为简单行星排。实验工况设定输入转速为1 500 r/min，负载为900 N·m，此时齿轮振动最为剧烈，便于提取故障特征。行星变速箱挡位为4挡，此时只有K3行星排以行星轮系方式工作。故障设置在K3行星排齿数为30，太阳轮和齿数为15行星轮的某轮齿齿面上。单向振动传感器安装在箱体内部K3行星排齿数为60的内齿圈正上方行星框架上，该测点距振源近，受传递路径干扰小，能更好地采集振动信号。实验台具体情况如图2所示。

行星轮系的运动方式不同于定轴轮系为复合运动，行星轮既自转又公转，采集的数据应包括同一行星齿轮回到初始位置的完整周期，以保证数据的可用性。如果采样所取完整周期数不同，测得信号对应的物理过程不同，故障产生的冲击数也不同，这样将多个振动信号的样本熵做横向比较则失去意义。依据齿轮相关参数和理论计算，可得行星轮旋转到初始位置的周期历时0.056 6 s. 本次实验采样频率20 kHz，每次采样取30个周期，即每次采样1.7 s，连续采集50组数据用于处理。

2.2 数据分析

2.2.1 VMD算法参数优化

本文以太阳轮断齿故障信号为例，详细介绍数据处理过程。使用PSO算法选取VMD算法参数，分别处理50组太阳轮断齿信号，得到50组数据的优化结果如图3所示。

由图3可知，k值很稳定，除3组数据外均为8，而α值在较小范围波动，基本保持稳定在2 510.

由于PSO算法可能存在局部最优情况，仅根据迭代次数终止迭代确定优化结果不尽合理，本文采用一种基于人工观察的算法进一步验证PSO算法结果，避免PSO算法优化结果陷入局部最优。

由VMD算法分解得到的本征模函数(IMF)分量的中心频率是由低频至高频分布的，如预设尺度数k从小到大取值，最后一层IMF分量的中心频率首次达到最大值时，将不会出现分解不足问题。若随着k值增大，最大中心频率仍然保持相对稳定，则可认为此时k值为最佳。

当k取不同值时，对太阳轮故障信号进行VMD算法分解，得到的IMF分量中心频率最小值和最大值如图4所示。

当k=5时，取到中心频率最小值，并随着k值增大，IMF分量中心频率最小值趋于稳定。当k=8时，取到中心频率最大值，并随k值增大中心频率最大值趋于稳定，保证VMD算法分解不会遗漏更高或更低的中心频率。当k=8时，各分量中心频率如图5所示。由图5可知：中心频率为1 kHz时，对应行星变速箱主动轴定轴齿轮的啮合频率，是频谱低频段的主要特征频率；中心频率为8 kHz时，对应着信号频谱最后一个明显峰值，符合信号实际情况，并且由中心频率迭代过程可知没有发生中心频率混叠。

行星变速箱正常工作状态时，齿轮振动幅度较小，振动信号相对复杂，并无明显规律性，故其样本熵值较大。而如果齿轮发生故障，出现与故障相关的周期性冲击，振动信号呈现出较强的规律性，样本熵值较正常状态变小。但是，如果信号干扰和噪声较强，与故障相关的周期性冲击脉冲不凸显，则故障信号的样本熵值仍较大，与正常状态区不明显。故选取最佳二次惩罚因子α更好地去除背景噪声的同时凸显了冲击成分，使降噪后信号的样本熵区分故障状态能力最大，即降噪后信号的样本熵对故障状态最为敏感。

为量化样本熵对行星变速箱状态的分类能力，引入双样本Z值(敏感度)检验法。双样本Z值检验法能有效地对两类样本在统计上的差异进行评估，计算结果Z值可作为特征参数分类能力的评价指标。在相同情况下，Z值越大，两类样本分类距离越大，说明特征参数分类能力越强。

X1{x11,x12,…,x1j}和X2{x21,x22,…,x2j}分别为两类样本，定义特征参量的Z值为

(9)

本文在预设尺度k=8条件下采用不同二次惩罚因子处理实验50组太阳轮故障和正常工况数据，得到降噪后信号的样本熵区分故障的敏感度均值如图6所示。

由图6可知，当α=2 500时，降噪后信号的样本熵对故障最敏感。

PSO算法与人工观察方法所得结果相近，证明了PSO算法的可用性，综合两种方法，避免陷入局部最优，又能同时考虑两种参数相互影响，保证VMD算法的分解性能。

2.2.2 信号重构

为进一步有效抑制信号中噪声和干扰成分，提高特征参量对齿轮故障状态的敏感度，以便更好地提取故障特征，需选取分解结果中的有效分量进行重构。对有效分量的选取常采用相关系数法，而相关系数仅能反映分量与原信号的线性关系，不能反映分量与原信号非线性关系。互信息法可衡量一个分量包含关于另一个分量的信息量，能够定量地表示两个变量之间的相互依赖程度，比相关系数法更加精确[15]，所以本文采用互信息法选取有效分量。变量X、Y之间的互信息定义如下：

(10)

式中：pXY(x,y)为变量X和变量Y的联合概率密度函数；pX(x)和pY(y)分别为变量X和变量Y的边缘概率密度函数。

各分量与原信号的互信息值越小者，表明它与原信号的依赖程度越低，包含原信号的可用信息越少，可能属于噪声干扰成分。反之，说明分量与原信号依赖程度越高，包含原信号的信息越多，为有效分量[16]。所以计算各分量与原始信号的互信息，可以用于筛选有效分量。

计算各分量与原信号的互信息MIi，并做归一化处理：

(11)

(12)

式中：λi为归一化后MIi中第i个分量与原信号的互信息；λm为归一化的互信息阈值。若λi>λm，则认为该分量为有效分量，否则予以剔除。

计算所得各分量和原信号的互信息如表1所示。

表1 各分量与原信号的互信息

由阈值确定方法计算处理太阳轮故障信号的阈值为0.515 6，选取大于阈值的前5阶有效分量进行重构，得重构前后信号时域波形如图7所示。

由图7可知，重构后的信号周期性更加明显，噪声干扰得到有效抑制，证明了VMD算法重构方法的有效性。

使用集合经验模态分解(EEMD)方法重构信号，添加的白噪声幅值系数取0.2，总体平均次数取300，重构阈值选取方法与VMD算法一致，得到对比如表2所示。

表2 重构方法对比

由表2可知：VMD算法较EEMD重构信噪比更高，耗时更短；且VMD算法重构后信号样本熵的均值和标准差更小，敏感度也更大，说明VMD算法较EEMD方法可以更好地去除背景噪声，凸显冲击成分，稳定性更好。

采用同样方法对正常状态，太阳轮断齿，行星轮断齿的信号分别进行参数优化VMD算法重构，得到重构后的信号频谱局部放大如图8所示。

由图8可知，行星变速箱在不同状态下的振动信号频率主要成分均为定轴部分特征频率的2倍频(900 Hz)和3倍频(1 350 Hz)。太阳轮断齿时，出现了明显的特征频率(529 Hz)及其2倍频(1 058 Hz)、3倍频(1 587 Hz)，并且定轴部分特征频率的2倍频(900 Hz)处边频带更加突出，调制现象更加明显。行星轮断齿时，只出现了微弱3倍频(1 587 Hz)与正常状况频谱差异不大，特征频率微弱。原因在于行星变速箱的行星轮系振动相对于定轴轮系微弱很多，加上行星排上离合器摩擦片的衰减作用，反映在频谱峰值就更加微弱。而K3行星排有6个行星轮，设置多个行星轮的目的就是为了使动力传递更加平稳，一个行星轮断齿对变速箱的影响不如太阳轮断齿大，由故障产生的振动也相对微弱，加上振动传感器位置固定不能有效采集行星轮的振动信号，这也是行星轮故障特征难于提取的原因。为更好地区分不同故障类型，还需结合多尺度熵偏均值进一步分析处理。

2.2.3 多尺度熵偏均值

对重构后的50组信号计算多尺度熵，得到50组信号在20个尺度下的多尺度熵平均值如图9所示。

同一信号在不同尺度因子下的复杂度是不同的[17]，并且不同信号的复杂度随尺度因子变化趋势也是不同的，说明不同类型信号的复杂度与尺度因子有着不同的相关性，考虑多尺度熵值变化趋势是非常有必要的。多尺度熵对太阳齿轮故障的区分度较好，但在区分行星轮故障和正常状况时出现了部分混叠，不能做到准确区分。为了更好地区分行星轮故障，引入能同时综合反映多尺度熵值大小和变化趋势的多尺度熵偏均值来衡量信号复杂度。

正常状态时，太阳轮故障、行星轮故障的50组数据分别使用样本熵、VMD算法重构样本熵、VMD算法重构多尺度熵偏均值方法处理结果如图10所示。

由图10可知，经过VMD算法重构后的样本熵减小，区分不同故障能力显著增强，由此表明VMD算法在保留不同状态信号特征差异的基础上，有效降低了信号中的随机成分，在更好地去除背景噪声的同时凸显了冲击成分，证明了重构方法的有效性。断齿故障状态的多尺度熵偏均值和样本熵值均小于正常状态，因为发生齿轮断齿故障时，振动信号中会产生与故障相关的有规律冲击[18]，导致信号复杂程度降低，多尺度熵偏均值和样本熵值减小。太阳轮故障较行星轮故障产生的冲击更明显，调制现象也更明显，信号复杂程度降低更多，熵值也更小，样本熵和多尺度熵偏均值都能有效区分太阳轮故障。样本熵只能反映信号在单一尺度下的复杂度，而多尺度熵偏均值可以反映信号在多尺度下的熵值和变化趋势，从而能够更全面反映信号的复杂度。行星轮故障信号更加复杂，故障特征更加微弱，更难于提取。单一尺度的样本熵反映的信号特征信息较少，不能有效区分行星轮故障，存在一定局限性。而多尺度熵偏均值对行星轮故障有较好区分能力，证明了多尺度熵偏均值处理非线性非平稳复杂信号相对于样本熵的优势。

使用本文方法处理不同工况下的行星轮断齿信号，结果如表3所示。

表3 不同工况对比

由表3可知，在行星变速箱不同挡位、不同转速、不同负载工况下，多尺度熵偏均值具有较高敏感度。当变速箱挡位设置在4挡，只有K3行星排传动方式为行星传动，K1行星排和K2行星排整体旋转，可以排除其他行星排的干扰，而在其他挡位时，多个挡位同时工作，振动信号成分更为复杂，尤其在3挡时3个行星排一起工作，样本熵基本不能区分故障，存在自身局限性。而在同一挡位下，转速越高，负载越大，齿轮断齿产生的冲击越明显，采集到的信号越规律，样本熵区分齿轮状态的敏感度越高。与样本熵相比，多尺度熵偏均值能从更深层次反映信号特性，受转速和负载影响较小，从而进一步证明了本文方法可以更加有效应用于行星变速箱故障特征提取。

3 结论

本文采用参数优化VMD算法和多尺度熵偏均值方法处理行星变速箱齿轮故障实验数据，结果表明：

1)使用PSO算法和人工观察方法优化VMD算法参数可行，摒除凭借经验选取的弊端，保证了VMD算法分解性能最优。同时使用互信息法对VMD算法分解后的分量进行筛选，有效去除了信号的噪声干扰成分，使与故障有关的规律性冲击信号更为明显，说明了重构方法的有效性。

2)多尺度熵偏均值将多尺度熵与偏均值结合，可以展现信号在不同尺度下的特征，能够更加完整和全面地反映信号的复杂性，从而反映行星变速箱的工作状态。通过对不同实验工况采集信号进行分析，多尺度熵偏均值衡量信号复杂程度能力强于样本熵，有效区分了行星轮故障和太阳轮故障，值得深入研究。