基于改进VMD-SVD降噪的齿轮箱故障特征提取*

2021-03-26蒋丽英潘宗博刘佳鑫

组合机床与自动化加工技术 2021年3期

蒋丽英，潘宗博，刘佳鑫

(沈阳航空航天大学自动化学院，沈阳 110136)

0 引言

齿轮传动在机械设备运行中是重要的环节之一，当齿轮发生故障时会使机械设备运行不正常，甚至造成重大的财产损失和人员伤亡。因此，对齿轮故障特征提取及诊断进行研究有着重要的意义。齿轮箱的振动信号常常表现出很强的非平稳非线性，故障特征频谱通常以调频调幅的形式存在于啮合频率及其倍频两边[1]，由于信号传播路径等因素的影响使得这个故障特征被噪声干扰，造成提取故障特征以及辨识故障类型困难的问题。

信号降噪处理的方法有很多，每种方法都有各自的优势和不足。经验模态分解(EMD)[2]能将待处理信号分解为一系列固有模态函数(IMF)，能对待处理信号进行不同尺度的描述；然而，EMD存在边界效应及模态混叠现象。固有时间尺度分解(ITD)[3]能将待处理信号自适应地分解为一系列频率由高到低的旋转分量(PR)和一个单调趋势的残余量，克服了EMD的模态混叠现象；但由于ITD的基线信号是基于相邻极值点间每一段信号的线性变换，PR分量有可能出现明显的波形失真，且仍存在端点效应[4]。小波分解(WT)[5]通过将原始信号分解成一系列高频与低频分量，在以一定的准则选取单个或多个分量重构信号实现降噪过程；若故障特征频率为低频，小波降噪可以去除大量噪声成分；若故障特征频率为高频，小波降噪在保留故障特征的同时也会保留较宽频率段的噪声。

本文基于变分模态分解(VMD)[6]和奇异值分解(SVD)[7]对齿轮箱振动信号的降噪问题进行研究。VMD将模态的估计转变为变分问题的求解，能够准确的将信号中各频率成分分解为不同中心频率和带宽的分量，具有很强的自适应性，较好地克服了EMD和ITD在信号分解中的不足。SVD是线性代数中最重要的工具之一, 在统计分析、信号处理和故障检测等方面都有重要应用[8]；在信号处理方面，SVD主要用于周期成分提取和信号降噪，具有极好的不变性和稳定性[9]；分解后选取合适的奇异值阈值重构信号，可以有效地保留有用信息、消除原始信号中的噪声干扰；相比于小波降噪，SVD没有频率限制，适用于全频率段的降噪处理。

周洋等[10]和马增强等[11]使用了传统的VMD-SVD降噪方法应用于滚动轴承振动信号，并成功提取出故障特征；但是该方法在强噪声背景下难以判定VMD处理后的各IMF分量是否保留，且由于齿轮箱啮合频率处的幅值远高于其他频率的特点难以找到合适的奇异值阈值用于SVD降噪处理，因此该方法并不适用于强背景噪声下齿轮箱振动信号的降噪处理。为有效提取出强背景噪声干扰下的齿轮箱振动信号的故障特征，本文提出了改进VMD-SVD的降噪方法。以齿轮箱断齿故障数据为试验样本，对降噪后的信号进一步做希尔伯特解调处理[12]，准确的提取出故障特征。

1 变分模态分解基本原理

变分模态分解是一种信号分解方法，通过迭代可以将原始信号分解为拥有不同中心频率和带宽的模态分量(IMF分量)。其以各模态分量的估计带宽之和最小、各模态分量之和等于输入信号为约束。受约束的变分模型如式(1)所示[13]：

(1)

式中，∂t为t的偏导数，δ(t)为脉冲函数，μk为第k个模态函数，ωk每个模态的中心频率，f为原始信号。

在此基础上，引入二次惩罚因子α和Lagrange乘法算子λ(t)，将上述变分模型转化为非约束变分模型，扩展的Lagrange表达式如式(2)所示：

(2)

式中，α参数用于保证重构信号的精度，λ(t)使约束条件更加严格。

(3)

过程如下：

②设置循环条件n=n+1；

(4)

④更新ωk：

(5)

⑤更新λ：

(6)

式中，τ为噪声容限参数，设τ=0时可以达到良好的降噪效果。

⑥循环步骤②～⑤，满足判定表达式时终止迭代。判定表达式如式(7)所示。

(7)

式中，e为判别精度。

最终，变分模态分解将原始信号f分解为K个窄带模态分量μk，其中，k∈1,2,…,K。

2 奇异值分解原理

2.1 奇异值分解基本原理

奇异值分解是一种分解矩阵的方法，可以将维数为m×n的矩阵X分解得到如下形式：

X=U×Σ×VT

(8)

式中，U为m×m阶酉矩阵，Σ为半正定m×n阶对角矩阵，V为n×n阶酉矩阵。Σ的对角线元素σ1,σ2,σ3,…,σi，(i=min(m,n))即为矩阵X的奇异值，且满足σ1≥σ2≥σ3≥…≥σi。

2.2 奇异值降噪原理

奇异值包含了分解前矩阵的全部信息，通过对奇异值的处理对原始数据进行分析，可以实现预期的需求。针对一维时域数据，奇异值降噪原理如下：

(1)由式(8)可以看出，奇异值分解是对矩阵进行处理，首先要将一维时域数据转化为矩阵形式。可以通过构造Hankel矩阵实现上述要求，将时域数据构造成Hankel矩阵的过程如下：

①假设时域信号数据中有L个采样点数据，对其编号为x(1),x(2),x(3),…,x(L)。

②按Hankel矩阵形式排列L个采样点数据，如式(9)所示。

(9)

式中，矩阵的行数满足1

(2)对构造的Hankel矩阵进行奇异值分解处理，得到U矩阵，Σ矩阵和V矩阵。

(3)选取Σ矩阵中的某一个奇异值σi为阈值，将矩阵Σ中大于等于阈值的奇异值保留，其余的奇异值置零，如式(10)所示：

Σ′=diag(σ1,σ2,…,σi,0,0,…,0)

(10)

(4)使用新构造的奇异值矩阵Σ′和奇异值分解生成的酉矩阵U和酉矩阵V做奇异值分解反运算，如式(11)所示：

H′=U×Σ′×VT

(11)

(5)新生成的H′矩阵负对角元素并不相同，对H′矩阵的负对角元素求均值得到降噪后的信号x′(1),x′(2),x′(3),…,x′(L)，如式(12)所示[14]：

(12)

式中，α=max(1,i-m+1)，β=min(n,i)。

2.3 奇异值降噪阈值选取原理

常用的降噪阈值选取方法有基于奇异值差分谱的方法、基于奇异值均值的方法和基于奇异值中值的方法。基于奇异值差分谱的方法有两种阈值确定方式：采用差分谱最大值对应的奇异值为降噪阈值，该阈值确定方式可以自适应的确定阈值，但不适用于频率调制信号，且噪声比较强时使用该阈值进行奇异值降噪处理时无法保留全部有用信息；从右向左观察奇异值差分谱，选取第一个较大的峰值对应的奇异值为降噪阈值，该阈值确定方式主观性强，选取的阈值可能不是最佳阈值，影响降噪效果。基于奇异值均值的方法将阈值定为所有奇异值的均值，基于奇异值中值的方法将阈值定为所有奇异值的中值，使用这两种方法选取的阈值做降噪处理时，信号的降噪效果不理想。

基于以上不足，本文提出了一种新的奇异值降噪阈值选取方法，可以自适应的确定适合的奇异值阈值，方法原理如下：

(1)根据信号特点确定信号保留比例γ。

(2)计算所有奇异值的和，如式(13)所示：

(13)

(3)选取第a阶奇异值为降噪阈值，其中a满足式(14)：

(14)

3 改进VMD-SVD降噪方法

本文针对强噪声背景下齿轮箱故障振动信号的特征提取问题，提出了改进VMD-SVD的降噪方法。首先使用VMD将原始信号分解为多个IMF分量，然后使用SVD对每个IMF分量单独进行降噪处理，最后将处理后的IMF分量线性叠加。具体实现步骤如下：

(1)获取试验样本数据。

在故障振动信号数据中随机截取长度为L的时间序列XL=[x(1),x(2),…x(L)]作为试验样本。

(2)对样本数据XL进行VMD分析。

首先利用中心频率观察法确定参数K；然后使用VMD分解XL，得到K个IMF分量，记作IMFi(i=1,2,…,K)。

(3)构造Hankel矩阵。

按式(9)方式分别构造各IMF的Hankel矩阵Hi(i=1,2,…,K)。

(4)使用SVD分解Hi(i=1,2,…,K)。

对Hi(i=1,2,…,K)进行奇异值分解，分别得到Hi(i=1,2,…,K)对应的奇异值对角矩阵和两个酉矩阵：Σi、Ui和Vi(i=1,2,…,K)。

(5)选取奇异值降噪阈值。

首先确定信号保留比例γ；然后按式(13)方式求取所有的奇异值之和sum；最后从小到大迭代a，满足式(14)的a对应的奇异值σa即为选取的奇异值降噪阈值。

(6)重构奇异值对角矩阵。

(7)反奇异值分解

(8)重构信号

通过线性叠加降噪后的各IMF分量重构信号，即为改进VMD-SVD降噪方法处理后的故障信号。

改进VMD-SVD降噪方法的示意图如图1所示。

图1 基于VMD和SVD的降噪方法示意图

4 仿真信号分析

为了验证改进VMD-SVD降噪方法的降噪性能，本文构造了一个拥有较强噪声的仿真信号，通过仿真验证了该方法的有效性。

构造的仿真信号为：

x(t)=x1(t)+x2(t)+δ(t)
x1(t)=0.2sin(1200πt)
x2(t)=0.3sin(3000πt)

(15)

式中，x1(t)为频率为600 Hz的正弦信号，x2(t)为频率为1500 Hz的正弦信号，δ(t)为高斯白噪声信号。

得到的仿真信号x(t)的时域波形和频率谱密度如图2所示，结合频率谱密度可以看出信号噪声很大，几乎要掩盖600 Hz的正弦信号。

图2 仿真信号的时域波形和频率谱密度

使用变分模态分解对x(t)进行分解，通过中心频率观察法确定最佳分解参数K=2，分解得到模态分量IMF1和IMF2。分别对IMF1和IMF2分量构造Hankel矩阵并进行奇异值分解。选择信号保留比例λ=10%，确定IMF1分量的奇异值降噪阈值阶数为12，IMF2分量的奇异值降噪阈值阶数为13。将反奇异值分解生成矩阵转化为时域信号即完成了奇异值分解对各模态分量的降噪处理。将降噪后的IMF1分量和IMF2分量线性叠加在一起，得到降噪处理后的故障信号，信号的时域波形和频率谱密度如图3所示。与图2相比，有用信号的幅值被最大程度地保留了下来，而噪声信号被滤掉大部分，验证了本文方法有效性。

图3 降噪后故障信号的时域波形和频率谱密度

5 实测信号分析

为验证本文所提方法的实用性，使用转频为1450 r/min、采样频率5120 Hz、主动轮齿数为55、从动轮齿数为75的齿轮箱断齿数据进行改进VMD-SVD降噪处理。经理论计算，断齿齿轮转频为17.5 Hz，啮合频率为1328 Hz。原始故障数据的时域波形、频率谱密度和解调频率谱密度如图4所示。

图4 原始故障数据的时域波形、频率谱密度和解调频率谱密度

可以看出，原始故障数据的噪声成分很少，啮合频率、故障齿轮转频及其2倍频和3倍频的幅值都很突出。为了突出本文所提方法的效果，在原始故障数据中加入了加性高斯白噪声，加入噪声的故障信号的时域波形、频率谱密度和解调频率谱密度如图5所示，信号的频率谱密度中除了啮合频率其他信号成分几乎淹没在噪声中，解调频率谱密度中故障转频也没有显现出来。

图5 加噪数据的时域波形、频率谱密度和解调频率谱密度

使用变分模态分解对加噪数据进行分解，通过中心频率观察法确定最佳分解参数K=3，分解得到模态分量IMF1、IMF2和IMF3。分别对IMF1、IMF2和IMF3分量构造Hankel矩阵并进行奇异值分解。选择信号保留比例λ=55%，确定IMF1分量的奇异值降噪阈值阶数为122，IMF2分量的奇异值降噪阈值阶数为111，IMF3分量的奇异值降噪阈值阶数为129。将反奇异值分解生成矩阵转化为时域信号即完成了奇异值分解对各模态分量的降噪处理。线性叠加降噪后的各模态分量得到降噪后的故障信号，时域波形、频率谱密度和解调频率谱密度如图6所示。与图5相比较，啮合频率、故障齿轮转频及其2倍频的幅值重新凸显出来，故障齿轮转频的3倍频有些许偏差。