CN114563181A - 基于改进变分模态提取的旋转机械故障特征提取方法 - Google Patents

Abstract

基于改进变分模态提取的旋转机械故障特征提取方法，包括：S1：采集振动信号；S2：时频变换处理得到时频谱；S3：从故障冲击特征所在的频率范围内选取一个频率值；S4：设置惩罚参数的变化范围和增大步长；S5：针对惩罚参数增长的每一步值执行变分模态提取，并计算对应期望模态的SDE指标值；S6：画出惩罚参数与SDE指标之间的关系曲线；S7：SDE指标最小值对应惩罚参数的最优值；S8：执行参数优化后的变分模态提取；S9：对最优期望模态进行平方包络谱分析，提取故障特征频率。本发明克服了现有的变分模态分解和变分模态提取中参数难以选取的问题，可以方便有效地从旋转机械振动信号中提取出故障特征，实现旋转机械故障诊断。

Description

基于改进变分模态提取的旋转机械故障特征提取方法

技术领域

本发明属于机械故障诊断技术领域，尤其涉及一种基于改进变分模态提取的旋转机械故障特征提取方法。

背景技术

齿轮与轴承是旋转机械设备中的重要组成零件，其健康状态对于设备的正常运行具有重要影响。但在实际工作中，这些零件经常要承受动态载荷甚至过载，容易产生各种类型的故障损伤，影响设备的使用效益，因此，对它们进行状态监测与故障诊断是非常有必要的。但由于设备中通常包含多个旋转运动的零件，不同零件激励的振动相互耦合叠加，而且故障损伤激发的冲击特征信号从损伤处传递到振动传感器的过程中，会有较大程度的衰减，再加上环境噪声的干扰，振动传感器采集的振动信号中，故障冲击特征通常比较微弱，尤其是在设备故障损伤早期。因此，从设备振动信号中检测或提取出故障相关的规律性冲击特征，进而识别故障特征频率，并不是一件容易的事情，这方面的研究也是旋转机械设备故障诊断研究中的一项重要内容。

目前，旋转机械设备故障特征提取的一个常用方法是变分模态分解(VMD)。VMD作为一种信号自适应分解方法，能够将信号分解成一系列具有不同中心频率的窄带本征模态函数，具有完备的数学理论支撑，且对信号采样要求不高、抗噪性好。但VMD也存在诸多不足之处，首先，VMD在执行过程中需要预先设定模态分解的层数，层数过多会导致过分解，产生虚假模态分量干扰；层数过少则会由于欠分解而导致模态混叠，难以提取目标模态；其次，VMD结果中有多个模态，有用的目标模态也不太好选取；此外，虽然可以采用各种优化算法对VMD的参数进行优化，但这些优化算法普遍计算量大、耗时长，严重降低VMD的使用效率。作为VMD的改进，变分模态提取(VME)最初提出用于ECG(心电图)信号的分析，虽然VME与VMD有着相同的理论基础，但VME只将信号分解成期望模态和残余信号两层，避免了VMD模态层数难以确定、目标模态难以选取、参数难以优化的问题。即便如此，VME依然存在期望模态中心频率初始值难以确定、惩罚参数难以优化的问题，尤其是在旋转机械振动信号分析中，这两个问题的解决是成功提取故障特征与实现故障诊断的关键。

发明内容

本发明要克服现有技术中变分模态分解和变分模态提取在旋转机械故障特征提取中存在的不足，提出一种基于改进变分模态提取的旋转机械故障特征提取方法。

本发明的基于改进变分模态提取的旋转机械故障特征提取方法，包括以下步骤：

S1：在旋转机械设备上，以f_s为采样频率，采集一段振动信号x(t)；

S2：对振动信号x(t)进行时频变换处理，得到时频谱TFP(t,f)，其中t表示时间，f表示频率；

S3：利用时频谱TFP(t,f)，确定故障冲击特征所在的频率范围，在该频率范围中间选取一个频率值f₀；

S4：设置惩罚参数α的变化范围为[α_min,α_max]，令α从α_min开始，以步长s_α增大，增大到第i步时的值为：

α_i＝α_min+(i-1)·s_α (I)

其中，i＝1,2,3,…,(α_max-α_min)/s_α；

S5：以2πf₀作为期望模态u_d(t)中心频率ω_d的初始值，以α_i作为惩罚参数，执行变分模态提取，得到一个期望模态u_{d_i}(t)，计算u_{d_i}(t)的SDE指标值；

S6：根据最终得到的(α_max-α_min)/s_α个SDE指标值，画出α与SDE指标之间的关系曲线；

S7：找出在惩罚参数α增大的过程中，SDE指标最小值所对应的α值，以该值作为惩罚参数α的最优值α_opt；

S8：以2πf₀作为最优期望模态u_{d_opt}(t)中心频率ω_d的初始值，以α_opt作为惩罚参数α的值，执行变分模态提取，得到最优期望模态u_{d_opt}(t)；

S9：对最优期望模态u_{d_opt}(t)进行平方包络谱分析，提取故障特征频率。

优选地，步骤S2中，时频变换处理方法包括短时傅里叶变换，计算公式为：

其中，t表示时间，f表示频率，g(τ-t)表示以时间t为中心的窗函数；针对短时傅里叶变换所得的时频矩阵，对其各个元素求模得到时频谱TFP(t,f)，表示为：

TFP(t,f)＝|STFT(t,f)| (III)

进一步优选地，步骤S5、步骤S8中，变分模态提取将振动信号x(t)分解成期望模态u_d(t)和残余信号x_r(t)两层，变分模态提取包括以下步骤：

S3-1：设定期望模态u_d(t)中心频率ω_d的初始值为

惩罚参数为α，期望模态u_d(t)的傅里叶变换初始化为

拉格朗日乘子λ(t)的傅里叶变换初始化为

更新步数n＝0；

S3-2：更新步数n＝n+1，整个算法开始循环执行；

S3-3：对所有的ω≥0，

按照公式(Ⅳ)进行更新：

其中，

为原始信号x(t)的傅里叶变换，

为第n步更新得到期望模态

的傅里叶变换，

为拉格朗日乘子λ(t)的傅里叶变换；ω_d按照公式(V)进行更新：

拉格朗日乘子

按照公式(VI)进行更新：

其中，σ为更新步长；

S3-4：给定精度值ε，判断公式(Ⅶ)表示的迭代更新终止条件是否满足：

如果不满足迭代更新终止条件，返回步骤S3-2，如果满足迭代更新终止条件，则结束整个循环，得到期望模态

S3-4：对频域的期望模态

进行傅里叶逆变换，得到时域的期望模态u_d(t)。

进一步优选地，步骤S5、步骤S6、步骤S7中，SDE指标值的计算包括以下步骤：

S4-1：利用Hilbert变换对期望模态u_d(t)进行解调，得到其上包络e_up(t)；

S4-2：在上包络e_up(t)中，找出大于某一阈值m_th的所有局部极大值所在位置，这些位置表示为p_k，k＝1,2,3,…,K，K为局部极大值的个数；

S4-3：计算位置序列p_k的差分序列，公式为：

d_j＝p_j+1-p_j (Ⅷ)

其中，j＝1,2,3,…,K-1；

S4-4：计算SDE指标值，公式为：

其中，

表示差分序列d_j的平均值，计算公式为：

进一步优选地，步骤S4-2中，阈值m_th设定为：

其中，max[e_up(t)]表示求上包络e_up(t)的最大值，min[e_up(t)]表示求上包络e_up(t)的最小值。

进一步优选地，步骤S9中，最优期望模态u_{d_opt}(t)的平方包络谱计算包括以下步骤：

S6-1：计算u_{d_opt}(t)的平方包络信号u_se(t)，公式为：

其中，j为虚数单位，

表示Hilbert变换，|·|表示求复数的模；

S6-2：计算u_{d_opt}(t)的平方包络谱U_se(ω)，公式为：

其中，

表示傅里叶变换。

本发明取得的积极效果在于：

(1)本发明利用变分模态提取方法，只将信号分解成期望模态和残余信号两层，在保持变分模态分解方法优点的同时，避免了变分模态分解的模态层数难以确定、目标模态难以选取的问题。

(2)本发明利用时频谱分析方法，能够有效展示旋转机械振动信号中故障冲击特征所在的频段，从而可以方便有效地选取变分模态提取的期望模态中心频率初始值。

(3)本发明所提出的SDE指标值可以有效地对变分模态提取的惩罚参数进行优化。

(4)本发明克服了现有的变分模态提取方法中期望模态中心频率初始值难以确定、惩罚参数难以优化的不足，利用改进变分模态提取方法得到优化的期望模态，再结合平方包络谱分析，可以成功地提取出故障特征频率，实现旋转机械故障诊断。

附图说明

图1为本发明的实施流程图；

图2为本发明实施例中，齿轮箱振动信号的时域波形图；

图3为本发明实施例中，齿轮箱振动信号的短时傅里叶变换时频谱图。

图4为本发明实施例中，变分模态提取中惩罚参数α与所提取期望模态SDE指标之间的关系曲线图。

图5为本发明实施例中，参数优化后变分模态提取所得到的最优期望模态时域波形图；

图6为本发明实施例中，最优期望模态的平方包络谱图。

具体实施方式

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

参照图1，基于改进变分模态提取的旋转机械故障特征提取方法，包括以下步骤：

S1：在旋转机械设备上，以f_s为采样频率，采集一段振动信号x(t)；

S2：对振动信号x(t)进行时频变换处理，得到时频谱TFP(t,f)，其中t表示时间，f表示频率；

S3：利用时频谱TFP(t,f)，确定故障冲击特征所在的频率范围，在该频率范围中间选取一个频率值f₀；

S4：设置惩罚参数α的变化范围为[α_min,α_max]，令α从α_min开始，以步长s_α增大，增大到第i步时的值为：

α_i＝α_min+(i-1)·s_α (I)

其中，i＝1,2,3,…,(α_max-α_min)/s_α；

S5：以2πf₀作为期望模态u_d(t)中心频率ω_d的初始值，以α_i作为惩罚参数，执行变分模态提取，得到一个期望模态u_{d_i}(t)，计算u_{d_i}(t)的SDE指标值；

S6：根据最终得到的(α_max-α_min)/s_α个SDE指标值，画出α与SDE指标之间的关系曲线；

S7：找出在惩罚参数α增大的过程中，SDE指标最小值所对应的α值，以该值作为惩罚参数α的最优值α_opt；

S8：以2πf₀作为最优期望模态u_{d_opt}(t)中心频率ω_d的初始值，以α_opt作为惩罚参数α的值，执行变分模态提取，得到最优期望模态u_{d_opt}(t)；

S9：对最优期望模态u_{d_opt}(t)进行平方包络谱分析，提取故障特征频率。

步骤S2中，时频变换处理方法包括短时傅里叶变换，计算公式为：

其中，t表示时间，f表示频率，g(τ-t)表示以时间t为中心的窗函数；针对短时傅里叶变换所得的时频矩阵，对其各个元素求模得到时频谱TFP(t,f)，表示为：

TFP(t,f)＝|STFT(t,f)| (III)

步骤S5、步骤S8中，变分模态提取将振动信号x(t)分解成期望模态u_d(t)和残余信号x_r(t)两层，变分模态提取包括以下步骤：

S3-1：设定期望模态u_d(t)中心频率ω_d的初始值为

惩罚参数为α，期望模态u_d(t)的傅里叶变换初始化为

拉格朗日乘子λ(t)的傅里叶变换初始化为

更新步数n＝0；

S3-2：更新步数n＝n+1，整个算法开始循环执行；

S3-3：对所有的ω≥0，

按照公式(Ⅳ)进行更新：

其中，

为原始信号x(t)的傅里叶变换，

为第n步更新得到期望模态

的傅里叶变换，

为拉格朗日乘子λ(t)的傅里叶变换；ω_d按照公式(V)进行更新：

拉格朗日乘子

按照公式(VI)进行更新：

其中，σ为更新步长；

S3-4：给定精度值ε，判断公式(Ⅶ)表示的迭代更新终止条件是否满足：

如果不满足迭代更新终止条件，返回步骤S3-2，如果满足迭代更新终止条件，则结束整个循环，得到期望模态

S3-4：对频域的期望模态

进行傅里叶逆变换，得到时域的期望模态u_d(t)。

步骤S5、步骤S6、步骤S7中，SDE指标值的计算包括以下步骤：

S4-1：利用Hilbert变换对期望模态u_d(t)进行解调，得到其上包络e_up(t)；

S4-2：在上包络e_up(t)中，找出大于某一阈值m_th的所有局部极大值所在位置，这些位置表示为p_k，k＝1,2,3,…,K，K为局部极大值的个数；

S4-3：计算位置序列p_k的差分序列，公式为：

d_j＝p_j+1-p_j (Ⅷ)

其中，j＝1,2,3,…,K-1；

S4-4：计算SDE指标值，公式为：

其中，

表示差分序列d_j的平均值，计算公式为：

步骤S4-2中，阈值m_th设定为：

其中，max[e_up(t)]表示求上包络e_up(t)的最大值，min[e_up(t)]表示求上包络e_up(t)的最小值。

步骤S9中，最优期望模态u_{d_opt}(t)的平方包络谱计算包括以下步骤：

S6-1：计算u_{d_opt}(t)的平方包络信号u_se(t)，公式为：

其中，j为虚数单位，

表示Hilbert变换，|·|表示求复数的模；

S6-2：计算u_{d_opt}(t)的平方包络谱U_se(ω)，公式为：

其中，

表示傅里叶变换。

将上述发明应用于一个齿轮箱试验台的振动信号处理中。试验齿轮箱为一个单级平行轴齿轮箱，工作在增速状态，传动比为2.45，设定试验台的驱动电机输出转速为900rpm，加载电机施加的载荷为3N·m，利用转速计测得试验齿轮箱的输出轴转速为1730r/min，其内部大齿轮的一个轮齿发生断裂损伤，通过计算可知该大齿轮的故障特征频率为11.77Hz，设定试验齿轮箱的振动信号采样频率f_s为5.12kHz。利用本发明提取故障特征的步骤如下。

S1：在试验齿轮箱上，以f_s为采样频率，采集到一段数据长度为4000的振动信号x(t)，如图2所示；

S2：对振动信号x(t)进行短时傅里叶变换处理，得到时频谱TFP(t,f)，如图3所示；

S3：利用时频谱TFP(t,f)，可以确定故障冲击特征所在的频率范围为1350～1450Hz，在该频率范围中间选取一个频率值f₀＝1400Hz；

S4：设置惩罚参数α的变化范围为[500,100000]，令α从500开始，以步长100增大，增大到第i步时的值为：

α_i＝500+100(i-1)

其中，i＝1,2,3,…,995；

S5：以2800πrad/s作为期望模态u_d(t)中心频率ω_d的初始值，以α_i作为惩罚参数，执行变分模态提取，得到一个期望模态u_{d_i}(t)，计算u_{d_i}(t)的SDE指标值；

S6：根据最终得到的995个SDE指标值，画出α与SDE指标之间的关系曲线，如图4所示；

S7：在惩罚参数α增大的过程中，SDE指标最小值为0.0092，其所对应的α值为61300，以该值作为惩罚参数α的最优值α_opt，即α_opt＝61300；

S8：以2800πrad/s作为最优期望模态u_{d_opt}(t)中心频率ω_d的初始值，以61300作为惩罚参数α的值，执行变分模态提取，得到最优期望模态u_{d_opt}(t)，如图5所示；

S9：对最优期望模态u_{d_opt}(t)进行平方包络谱分析，如图6所示，提取出的故障特征频率为11.52Hz，与试验齿轮箱中的大齿轮理论故障特征频率11.77Hz基本一致，说明该大齿轮发生了故障损伤。

本说明书实施例所述的内容仅仅是对发明构思的实现形式的列举，本发明的保护范围不应当被视为仅限于实施例所陈述的具体形式，本发明的保护范围也及于本领域技术人员根据本发明构思所能够想到的等同技术手段。

Claims (6)

1.基于改进变分模态提取的旋转机械故障特征提取方法，包括以下步骤：

S1：在旋转机械设备上，以f_s为采样频率，采集一段振动信号x(t)；

S2：对振动信号x(t)进行时频变换处理，得到时频谱TFP(t,f)，其中t表示时间，f表示频率；

S3：利用时频谱TFP(t,f)，确定故障冲击特征所在的频率范围，在该频率范围中间选取一个频率值f₀；

S4：设置惩罚参数α的变化范围为[α_min,α_max]，令α从α_min开始，以步长s_α增大，增大到第i步时的值为：

α_i＝α_min+(i-1)·s_α (I)

其中，i＝1,2,3,…,(α_max-α_min)/s_α；

S5：以2πf₀作为期望模态u_d(t)中心频率ω_d的初始值，以α_i作为惩罚参数，执行变分模态提取，得到一个期望模态u_{d_i}(t)，计算u_{d_i}(t)的SDE指标值；

S6：根据最终得到的(α_max-α_min)/s_α个SDE指标值，画出α与SDE指标之间的关系曲线；

S7：找出在惩罚参数α增大的过程中，SDE指标最小值所对应的α值，以该值作为惩罚参数α的最优值α_opt；

S8：以2πf₀作为最优期望模态u_{d_opt}(t)中心频率ω_d的初始值，以α_opt作为惩罚参数α的值，执行变分模态提取，得到最优期望模态u_{d_opt}(t)；

S9：对最优期望模态u_{d_opt}(t)进行平方包络谱分析，提取故障特征频率。

2.基于改进变分模态提取的旋转机械故障特征提取方法，其特征在于，步骤S2中，时频变换处理方法包括短时傅里叶变换，计算公式为：

其中，t表示时间，f表示频率，g(τ-t)表示以时间t为中心的窗函数；针对短时傅里叶变换所得的时频矩阵，对其各个元素求模得到时频谱TFP(t,f)，表示为：

TFP(t,f)＝|STFT(t,f)| (III)

3.基于改进变分模态提取的旋转机械故障特征提取方法，其特征在于，步骤S5、步骤S8中，变分模态提取将振动信号x(t)分解成期望模态u_d(t)和残余信号x_r(t)两层，变分模态提取包括以下步骤：

S3-1：设定期望模态u_d(t)中心频率ω_d的初始值为

惩罚参数为α，期望模态u_d(t)的傅里叶变换初始化为

拉格朗日乘子λ(t)的傅里叶变换初始化为

更新步数n＝0；

S3-2：更新步数n＝n+1，整个算法开始循环执行；

S3-3：对所有的ω≥0，

按照公式(Ⅳ)进行更新：

其中，

为原始信号x(t)的傅里叶变换，

为第n步更新得到期望模态

的傅里叶变换，

为拉格朗日乘子λ(t)的傅里叶变换；ω_d按照公式(V)进行更新：

拉格朗日乘子

按照公式(VI)进行更新：

其中，σ为更新步长；

S3-4：给定精度值ε，判断公式(Ⅶ)表示的迭代更新终止条件是否满足：

如果不满足迭代更新终止条件，返回步骤S3-2，如果满足迭代更新终止条件，则结束整个循环，得到期望模态

S3-4：对频域的期望模态

进行傅里叶逆变换，得到时域的期望模态u_d(t)。

4.基于改进变分模态提取的旋转机械故障特征提取方法，其特征在于，步骤S5、步骤S6、步骤S7中，SDE指标值的计算包括以下步骤：

S4-1：利用Hilbert变换对期望模态u_d(t)进行解调，得到其上包络e_up(t)；

S4-2：在上包络e_up(t)中，找出大于某一阈值m_th的所有局部极大值所在位置，这些位置表示为p_k，k＝1,2,3,…,K，K为局部极大值的个数；

S4-3：计算位置序列p_k的差分序列，公式为：

d_j＝p_j+1-p_j (Ⅷ)

其中，j＝1,2,3,…,K-1；

S4-4：计算SDE指标值，公式为：

其中，

表示差分序列d_j的平均值，计算公式为：

5.基于改进变分模态提取的旋转机械故障特征提取方法，其特征在于，步骤S4-2中，阈值m_th设定为：

其中，max[e_up(t)]表示求上包络e_up(t)的最大值，min[e_up(t)]表示求上包络e_up(t)的最小值。

6.基于改进变分模态提取的旋转机械故障特征提取方法，其特征在于，步骤S9中，最优期望模态u_{d_opt}(t)的平方包络谱计算包括以下步骤：

S6-1：计算u_{d_opt}(t)的平方包络信号u_se(t)，公式为：

其中，j为虚数单位，

表示Hilbert变换，|·|表示求复数的模；

S6-2：计算u_{d_opt}(t)的平方包络谱U_se(ω)，公式为：

其中，

表示傅里叶变换。

Images