CN109029977A - 一种基于vmd-amckd的行星齿轮箱早期故障诊断方法 - Google Patents

Abstract

本发明涉及一种基于VMD‑AMCKD的行星齿轮箱早期故障诊断方法，包括步骤S1获取行星齿轮箱的振动加速度信号，并确定VMD参数；步骤S2基于VMD参数，对原振动信号进行VMD分解，并根据相关系数最大准则，选取最优分量；步骤S3:根据GOA寻优算法对所述最优分量进行MCKD参数的寻优；步骤S4:根据步骤S3得到的MCKD参数，对最优分量信号进行MCKD分析，再对解卷积后信号进行包络解调；步骤S5:将轴承理论故障特征频率值与包络谱中峰值明显的谱线进行对照，从而诊断出故障类型，确定故障部位。本发明采用GOA算法克服了MCKD算法中参数难以确定问题，采用VMD实现了信号降噪及各频率带分离，采用MCKD进一步增强故障冲击成分，对于行星齿轮箱早期故障诊断更为精准。

Description

一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法

技术领域

本发明属于旋转类机械故障诊断领域，具体涉及一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法。

背景技术

行星齿轮箱体积小，承载能力强，传动比大等优点，故其在工业、运输业、军事领域都得到了广泛的运用。如果没有及时发现行星齿轮箱早期故障，那么其可能导致的后果不堪设想。行星齿轮箱早期故障振动信号十分微弱，极易被噪声掩盖；再加上实际工况的复杂性，使得行星齿轮箱早期故障信号难以被检测。

行星齿轮箱故障信号具有非线性、非平稳性，针对这一类信号，学者们将许多处理方法引入到故障诊断领域。常见信号处理的方法有：小波变换，EMD分解，EEMD分解，LMD分解等。尽管这些方法在轴承故障诊断领域取得了不错的效果，但是还存在如下问题：(1)小波分解难以根据不同实际信号选择小波基和分解层数；(2)EMD分解、EEMD分解、LMD分解均属于递归式模态分解，都存在模态混叠等问题，且没有严格的数学推导。不同于上述信号处理方法，变分模态分解(VMD)能够有效避免模态混叠、端点效应等问题。VMD算法中参数惩罚因子α及模态分解个数K的选择对算法结果影响很大，有学者采用智能优化算法实现参数的自适应选择，但增加了算法的复杂性。

MCKD通过构造一系列FIR滤波器，保留信号连续冲击成分，能够在低信噪比的情况下提取出故障特征信号。MCKD中滤波器长度参数L及解卷积周期T的选择对结果影响非常大。

发明内容

有鉴于此，本发明的目的在于提供一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，针对行星齿轮箱早期故障产生的冲击十分微弱，易被系统噪声干扰造成行星齿轮箱早期故障诊断困难问题，且VMD算法、MCKD算法参数难以确定问题。

为实现上述目的，本发明采用如下技术方案：

一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：包括如下步骤：

步骤S1：获取行星齿轮箱的振动加速度信号，并确定VMD参数；

步骤S2:基于VMD参数，对原振动信号进行VMD分解，并根据相关系数最大准则，选取最优分量；

步骤S3:根据GOA寻优算法对所述最优分量进行MCKD参数的寻优；

步骤S4:根据步骤S3得到的MCKD参数，对最优分量信号进行MCKD分析，再对解卷积后信号进行包络解调；

步骤S5:将轴承理论故障特征频率值与包络谱中峰值明显的谱线进行对照，从而诊断出故障类型，确定故障部位。

进一步的，步骤S1所述确定VMD参数具体为：设定VMD中的惩罚因子α₀＝2000，基于中心频率观察法确定模态分解个数K₀。

进一步的，所述中心频率观察法具体为：固定惩罚因子α₀，观察K＝K₀+1与K＝K₀时分解后的分量的频谱图，判断最后一个分量的中心频率是否相近，若相近，则确定模态分解个数为K₀。

进一步的，所述AMCKD算法具体包含GOA优化算法和MCKD算法，即利用GOA优化算法实现MCKD算法中两个参数L、T的自适应选择。其中，MCKD算法能够通过保留信号中的连续脉冲，提取出故障特征。

进一步的，所述相关系数最大准则，即计算模态分量信号序列与原信号序列的相关系数，相关系数越大，则该分量越优。其中，相关系数的计算公式为：

式中：C表示分量u_k与原信号的相关系数；E(·)表示数学期望；表示模态分量序列的平均值；表示原信号序列的平均值。

进一步的，所述GOA寻优算法为一种自然启发式算法，模拟自然界中的蚱蜢群寻找食物的行为，具体的模拟蚱蜢群集行为的数学模型如下所示：

P_i＝S_i+G_i+W_i

式中：P_i为第i只蚱蜢的位置；S_i为第i只蚱蜢受到的社会干预；G_i是第i只蚱蜢受到的地心引力；W_i为风向，设定为指向最优的方向。具体的，不考虑重力的影响，对社会干预进行进一步优化，有下面公式：

式中：P_i ^d表示蚱蜢i在第d维的位置；N为蚱蜢群规模；c为递减系数；ub_d和lb_d为第d维的上下位置边界；s(·)用于定义社会干预(吸引和排斥)；(p_i-p_j)/d_ij为第i只蚱蜢和第j只蚱蜢之间的单位向量；为指向优化目标的向量(该迭代为止最优的值)；为使GOA取得更好的优化效果，这里递减系数c为按线性变化的动态值，即：

式中：h为迭代次数；H为最大迭代数；c_max为最大递减系数；c_min为最小递减系数；参数c保证了GOA算法不能过快的收敛于目标，从而避免局部最优，并能在最后的几次迭代中加快收敛速度。

进一步的，所述GOA寻优算法参数设置为蚱蜢群规模N＝30，最大迭代次数H＝10，最大递减系数c_max＝1，最小递减系数c_min＝0.0004。

进一步的，所述GOA寻优算法的适应度函数为包络谱幅值比ESR的负值，ESR定义为：

式中：f_t为包络谱中最大幅值对应的频率，f_t∈[f_r,f_o]；f_r为输入电机转速，fo大于行星齿轮箱系统中的齿轮故障最大频率；E(f_t)到E(5f_t)为f_t到5f_t对应的包络谱幅值；E为包络谱有效幅值的总和，其中有效幅值为频率范围(0,f_a)的包络谱幅值，f_a应至少大于5f_o。

进一步的，所述步骤S5具体为：根据MCKD分析后的信号进行Hilbert包络解调并得到包络谱图，谱图中的特征频率与行星齿轮箱理论故障频率的对比，最终判定行星齿轮箱系统的故障状态，包括运行正常和存在故障两种状态。

本发明与现有技术相比具有以下有益效果：

(1)、本发明既能够发挥VMD在降噪方面的优越性，又能彰显MCKD能突出被噪声所掩盖的连续脉冲的优点，改善传统方法在行星齿轮箱微弱故障信号特征提取上的困难；能够克服单一使用MCKD或者单一使用VMD方法难以实现微弱故障诊断的缺点

(2)、本发明采用中心频率观察法确定VMD中的参数，克服了VMD算法参数难以选择问题，相比寻优算法具有运算时间短，效率高的特点。

(3)、本发明引入新型智能优化算法——GOA算法实现了MCKD参数的自适应选择，避免人为选择参数导致的错误诊断结果，GOA算法比一些其他寻优算法如：GA算法、PSO算法，寻优结果更加准确且具有更高的效率。

(4)、本发明在GOA算法中，提出新的适应度指标——ESR能够较好的反映故障特征频率对应幅值的强度，并限制适应度函数求解范围，避免了随机突触频率对两种寻优算法的影响。

附图说明

图1是本发明诊断方法流程图；

图2是本发明一实施例中行星齿轮箱太阳轮裂纹故障实验信号时域波形；

图3是本发明一实施例中太阳轮裂纹故障实验信号的频谱图；

图4是本发明一实施例中太阳轮裂纹故障实验信号的包络谱图；

图5为本发明一实施例中VMD对信号分解的最终结果图；

图6为本发明一实施例中各模态分量与故障信号的相关系数的直方图；

图7为本发明一实施例中GOA算法对MCKD算法中参数寻优的流程图；

图8为本发明一实施例中GOA算法对MCKD寻优，ESR随迭代次数变化曲线图；

图9为本发明一实施例中MCKD算法分析后的时域波形；

图10为本发明一实施例中MCKD算法分析后的包络谱。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

本发明提供一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，流程图如图1所示，其主要步骤如下：

步骤S1：获取行星齿轮箱的振动加速度信号，并确定VMD参数；

步骤S2:基于VMD参数，对原振动信号进行VMD分解，并根据相关系数最大准则，选取最优分量；

步骤S3:根据GOA寻优算法对所述最优分量进行MCKD参数的寻优；

步骤S4:根据步骤S3得到的MCKD参数，对最优分量信号进行MCKD分析，再对解卷积后信号进行包络解调；

步骤S5:将轴承理论故障特征频率值与包络谱中峰值明显的谱线进行对照，从而诊断出故障类型，确定故障部位。

实施例1：

搭建动力传动故障模拟试验台(DDS)，在太阳轮齿根位置处，人为加工微弱裂纹模拟行星齿轮箱早期裂纹故障。驱动电机实际输入转速为39.80Hz。使用加速度传感器采集行星齿轮箱箱体的振动信号，采样频率f_s为12800Hz，采样总点数为12800点。行星齿轮箱的齿数参数见表1，由齿数参数和输入转速可计算表2所示齿轮箱特征频率。

表1单级行星齿轮箱齿数参数

齿轮	太阳轮	行星轮	内齿圈
				齿数	28	36	100

表2行星齿轮箱相应特征频率(单位/Hz)

啮合频率fm	故障频率fp	行星架转频fc	太阳轮转频fsr
				870.63	124.38	8.71	39.80

本发明的基于VMD-AMCKD方法的行星齿轮箱早期故障诊断方法具体流程如图1所示，具体包括以下步骤。

步骤一：利用加速度传感器对所述旋转机械设备(DDS)传动轴的径向进行测量，以获得振动加速度信号。图2、图3、图4分别为实验信号的时域波形、频谱及其包络谱。图2中故障轴承冲击成分被噪声及其他成分掩盖，并不能从中观察到太阳轮裂纹故障信息。图3可以明显的观察到频率40Hz，该频率成分为输入电机的转频，由于频率分辨率，造成其与实际输入转速39.80Hz有点偏差；在啮频附近很难找到调制边频带。图4的包络谱在低频部分的频率成分十分冗杂，在故障频率124.38Hz附近更是杂乱，无法寻找有效的故障信息。故应进一步分析，先要确定VMD算法中的两个参数。本发明采用中心频率观察法确定模态分解个数K，即设定惩罚因子α₀＝2000，根据不同K值分解产生的最后一个模态分量的中心频率，判断出K的取值。所述VMD算法能够将一个实际信号x分解K个模态分量u_k，并且确定每个模态分量的中心频率ω_k和带宽。该算法模型构造及求解步骤如下：

1)由希尔伯特变换计算每个模态u_k(t)的解析信号，从而获得模态的单边频谱：

式中：δ(t)为脉冲函数；t为时间；j为虚数单位；“*”表示卷积。

2)向单边频谱添加一个指数项来进行频率混合，将各模态的频谱调制到相应的基频带：

3)计算解调信号的梯度，用其二范数平方，来估计各模态信号的带宽；所有分量相加等于原信号作为约束条件，约束变分模型描述如下：

式中：{u_k}＝{u₁,…,u_K}，{ω_k}＝{ω₁,…,ω_K}；表示函数对t求偏导。

4)为求解上式的变分模型，引入拉格朗日乘子λ(t)及二次惩罚因子α，将约束性变分问题变为非约束性变分问题。其中，α能够有效地降低高斯噪声干扰，λ(t)能够增强约束严格性。扩展的拉格朗日表达式为：

5)运用交替方向乘子算法(Alternate Direction Method of Multiplies，ADMM)，不断迭代更新{u_k ⁿ⁺¹}、{ω_k ⁿ⁺¹}、λⁿ⁺¹寻求上式的“鞍点”。迭代方法为：

式中：n为迭代次数；g∈[1,K]；Г为更新因子；ε为一个大于0的正数，代表精度。为更新{u_k ⁿ⁺¹}、{ω_k ⁿ⁺¹}、λⁿ⁺¹的收敛条件。

6)利用Parseval/Plancheral傅里叶等距变换；对第5步中三式求解可得：

式中：分别表示x(t)、λⁿ所对应的傅里叶变换。VMD算法是将各模态分量频域不断更新，然后通过傅里叶逆变换到时域。

7)综上所述，可制定VMD的具体实施流程如下：

7.1)初始化{u_k ¹}、{ω_k ¹}、λ¹，n＝0；

7.2)n＝n+1，开始整个算法的循环；

7.3)k＝k+1，直到k＝K，更新u_k、ω_k；

7.4)更新λ；

7.5)判断是否满足收敛条件，如果满足收敛条件则停止迭代，否则返回第7.2步。

根据不同K值对应的最后一个模态分量的中心频率，确定算法模态分解个数K0，不同K值对应的各分量如表3所示。

表3不同K值分量的中心频率

由表3可知，当K＝5的最后一个分量，与K＝6时的最后一个分量的中心频率，均为4312Hz；因此取K₀＝5。

步骤二：基于步骤一确定VMD算法中的参数，对故障信号进行VMD分解。如图5所示为VMD分解结果，从中观察右边模态分量的频率部分，不存在模态混叠现象。计算每个分量u_k与故障信号x的相关系数C，C越大，则判断该分量越优；C的计算公式如下：

式中：C表示分量u_k与原信号x的相关系数；E(·)表示数学期望；表示模态分量序列的平均值；表示原信号序列的平均值。

如图6所示为各分量与原信号相关系数的值，从中可知，分量2为最优分量。

步骤三：对分量2进行下一步的分析，基于GOA优化算法，确定MCKD参数中滤波器长度参数L和解卷积周期T。MCKD算法实质上是通过寻找一系列FIR滤波器f_l，使得原始周期性冲击序列的相关峭度最大。MCKD的具体算法步骤如下：1)周期冲击信号y_i(i＝1,2,…,N)的相关峭度定义为：

式中：T为冲击信号周期；M为移位数；m∈[0,M]；增加M，会增加算法序列脉冲数，考虑到M过大会影响精度，本实施例选取M＝7。

2)对实际信号进行滤波，使得相关峭度最大：

式中：y、x分别是周期信号y_i、实际信号x_i的向量形式；f＝[f₁,…,f_L]^T；L为滤波器长度参数；l∈[0,L]。

3)求解上式等价于求解(l＝1,2,…,L)，最终将求得的滤波器系数的结果使用矩阵的形式表示：

式中：r＝[0,T,…,mT]；

4)MCKD算法的具体实施过程为：

4.1)初始化解卷积周期T、位移数M、及滤波器长度L等参数；

4.2)计算输入信号x的X_T，X₀ ^T，(X₀X₀ ^T)^-1；

4.3)计算滤波后的输出信号y；

4.4)根据y计算ψ和β；

4.5)更新滤波器f的系数；

4.6)若滤波前与滤波后信号的相关峭度值ΔCK_M(T)小于阈值，结束迭代；否则，返回第4.3步。

图7所示为GOA算法对MCKD参数寻优的步骤流程图，GOA寻优算法为一种自然启发式算法，它模拟自然界中的蚱蜢群寻找食物的行为。模拟蚱蜢群集行为的数学模型如下所示：

P_i＝S_i+G_i+W_i

式中：P_i为第i只蚱蜢的位置；S_i为第i只蚱蜢受到的社会干预；G_i是第i只蚱蜢受到的地心引力；W_i为风向，设定为指向最优的方向。具体的，不考虑重力的影响，对社会干预进行进一步优化，有下面公式：

式中：P_i ^d表示蚱蜢i在第d维的位置；N为蚱蜢群规模；c为递减系数；ub_d和lb_d为第d维的上下位置边界；s(·)用于定义社会干预(吸引和排斥)；(p_i-p_j)/d_ij为第i只蚱蜢和第j只蚱蜢之间的单位向量；为指向优化目标(该迭代为止最优的值)的向量。为使GOA取得更好的优化效果，这里递减系数c为按线性变化的动态值，即：

式中：h为迭代次数；H为最大迭代数；c_max为最大递减系数；c_min为最小递减系数。参数c保证了GOA算法不能过快的收敛于目标，从而避免局部最优，并能在最后的几次迭代中加快收敛速度。

本发明中的GOA寻优算法参数设置为蚱蜢群规模N＝30，最大迭代次数H＝10，最大递减系数c_max＝1，最小递减系数c_min＝0.0004。

GOA寻优算法的适应度函数为包络谱幅值比ESR的负值，ESR定义为：

式中：f_t∈[f_r,f_o]，为包络谱中最大幅值对应的频率，f_r为输入电机转速39.80Hz，f_o为130Hz，大于行星齿轮箱系统中的齿轮故障最大频率；E(f_t)到E(5f_t)为f_t到5f_t对应的包络谱幅值；E为包络谱有效幅值的总和，其中有效幅值为频率范围(0,1000)的包络谱幅值。

图8所示为GOA算法优化MCKD参数迭代曲线图，GOA算法在第4代收敛，得到的最佳影响参数组合[L₀,T₀]为[720,103]。

步骤四：基于步骤三确定的参数，对最优分量进行MCKD分析。如图9，图10所示为MCKD分析后的时域波形及其包络谱。从图9可以发现7个冲击成分的总时间为0.0532，因此可以计算出特征频率为124.29Hz。进一步的，在包络谱中可以明显的发现124Hz及其倍频。

步骤五：时域波形的冲击间隔计算出的特征频率及所述包络谱的突出谱线，与太阳轮故障特征频率124.38Hz相近，从而判断该行星齿轮箱发生齿轮故障，故障类型为太阳轮裂纹。需要说明的是，图中突出谱线及其倍频与裂纹故障频率的差值可以认为是频率分辨率及在齿轮箱的传动误差造成的。

综上，本发明基于VMD-AMCKD的滚动轴承微弱故障诊断方法，能够成功提取行星齿轮箱中的早期微弱故障特征。VMD与MCKD结合能够更加精准的实现行星齿轮箱微弱故障诊断；GOA算法能够较好的实现MCKD参数的选择

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。

Claims (8)

1.一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：包括如下步骤：

步骤S1：获取行星齿轮箱的振动加速度信号，并确定VMD参数；

步骤S2:基于VMD参数，对原振动信号进行VMD分解，并根据相关系数最大准则，选取最优分量；

步骤S3:根据GOA寻优算法对所述最优分量进行MCKD参数的寻优；

步骤S4:根据步骤S3得到的MCKD参数，对最优分量信号进行MCKD分析，再对解卷积后信号进行包络解调；

步骤S5:将轴承理论故障特征频率值与包络谱中峰值明显的谱线进行对照，从而诊断出故障类型，确定故障部位。

2.根据权利要求1所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：步骤S1所述确定VMD参数具体为：设定VMD中的惩罚因子α₀＝2000，基于中心频率观察法确定模态分解个数K₀。

3.根据权利要求2所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：所述中心频率观察法具体为：固定惩罚因子α₀，观察K＝K₀+1与K＝K₀时分解后的分量的频谱图，判断最后一个分量的中心频率是否相近，若相近，则确定模态分解个数为K₀。

4.根据权利要求1所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：步骤S2所述相关系数最大准则，即计算模态分量信号序列与原信号序列的相关系数，相关系数越大，则该分量越优。其中，相关系数的计算公式为：

式中：C表示分量u_k与原信号的相关系数；E(·)表示数学期望；表示模态分量序列的平均值；表示原信号序列的平均值。

5.根据权利要求1所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：步骤S3所述GOA寻优算法为一种自然启发式算法，模拟自然界中的蚱蜢群寻找食物的行为，具体的模拟蚱蜢群集行为的数学模型如下所示：

P_i＝S_i+G_i+W_i

式中：P_i为第i只蚱蜢的位置；S_i为第i只蚱蜢受到的社会干预；G_i是第i只蚱蜢受到的地心引力；W_i为风向，设定为指向最优的方向。具体的，不考虑重力的影响，对社会干预进行进一步优化，有下面公式：

式中：P_i ^d表示蚱蜢i在第d维的位置；N为蚱蜢群规模；c为递减系数；ub_d和lb_d为第d维的上下位置边界；s(·)用于定义社会干预；(p_i-p_j)/d_ij为第i只蚱蜢和第j只蚱蜢之间的单位向量；为指向优化目标的向量；为使GOA取得更好的优化效果，这里递减系数c为按线性变化的动态值，即：

式中：h为迭代次数；H为最大迭代数；c_max为最大递减系数；c_min为最小递减系数；参数c保证了GOA算法不能过快的收敛于目标，从而避免局部最优，并能在最后的几次迭代中加快收敛速度。

6.根据权利要求6所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：所述GOA寻优算法参数设置为蚱蜢群规模N＝30，最大迭代次数H＝10，最大递减系数c_max＝1，最小递减系数c_min＝0.0004。

7.根据权利要求7所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：所述GOA寻优算法的适应度函数为包络谱幅值比ESR的负值，ESR定义为：

式中：f_t为包络谱中最大幅值对应的频率，f_t∈[f_r,f_o]；f_r为输入电机转速，f_o大于行星齿轮箱系统中的齿轮故障最大频率；E(f_t)到E(5f_t)为f_t到5f_t对应的包络谱幅值；E为包络谱有效幅值的总和，其中有效幅值为频率范围(0,f_a)的包络谱幅值，f_a应至少大于5f_o。

8.根据权利要求1所述的一种基于VMD-AMCKD的行星齿轮箱早期故障诊断方法，其特征在于：所述步骤S5具体为：根据MCKD分析后的信号进行Hilbert包络解调并得到包络谱图，谱图中的特征频率与行星齿轮箱理论故障频率的对比，最终判定行星齿轮箱系统的故障状态，包括运行正常和存在故障两种状态。