基于ITD-SVD和MOMEDA的故障特征提取方法

2022-01-20杨静宗杨天晴吴丽玫

工业工程 2021年6期

杨静宗，杨天晴，吴丽玫

(保山学院大数据学院，云南保山 678000)

轴承是机械设备中的重要部件之一，目前已在国民经济中的各个领域得到普及应用，但轴承也是最容易发生故障的部件[1]。其一旦发生故障，不仅会导致设备无法安全运行，还将使得生产效率大大降低。因此，对轴承的运行状态展开监测并及时掌握设备的运行状况显得尤为重要。通常，从被检测设备的信号中提取反映信号特征的参数是识别其是否发生故障的基础。目前，基于时频分析的特征提取法已逐渐成为研究的热点方向。其中，最具代表性的方法之一是基于小波分析的特征提取方法，但其在使用的过程中往往存在频带选取不够灵活，基函数选取标准不一等缺陷[2]。另一种方法则是由Huang等[3]提出的经验模态分解法(empirical mode decomposition，EMD)，其能够将信号分解为众多固有模态分量，上述分量直观地反映了信号中的时频分布情况。因此，在高压断路器故障识别、风电齿轮故障识别、模拟电路故障诊断等领域[4-7]均得到了普及应用。但该方法在实际使用过程中常面临端点效应、模态混叠等问题，而且没有严格的理论证明[8-10]。固有时间尺度分解(intrinsic time scale decomposition，ITD)[11]是在EMD方法基础上提出的一种非平稳、非线性信号分析方法，ITD不但能够像EMD方法一样自适应地分解信号，而且计算复杂度更低，传统方法存在的模态混叠和端点效应问题也得到了优化。近年来，张朝林等[12]为实现对风力发电机机组运行情况的实时监测，提出了ITD和多尺度熵相结合的信号分析方法，并通过构建的极限学习机模型识别风电轴承的运行状态。结果显示，该方法取得了较好的识别效果。任东方等[13]针对通信辐射源信号的非线性、非平稳特性，提出了ITD和非线性分析相结合的识别方法。结果表明，该方法能够有效提取信号中的深层次信息，且识别准确率较高。赵雄鹏等[14]将ITD和分形模糊熵理论引入至自动机的故障诊断中，并通过提取分形维数和模糊熵完成故障识别模型的建模。结果表明，该方法成功地识别出了自动机的故障类型。

然而，对于复杂工况下的滚动轴承信号，经过时频分析法分解后的信号分量中仍会残留一定的噪声。如果直接进行故障分析，噪声干扰将会对故障诊断结果的准确性产生影响。基于此，本文提出了一种基于ITD−奇异值分解(singular value decomposition，SVD)和多点最优最小熵反褶积(multipoint optimal minimum entropy deconvolution adjusted，MOMEDA)相结合的故障特征提取方法。首先，采用ITD分解轴承故障信号，以获取信号中的固有旋转分量(proper rotation components，PRC)。然后，利用K-C组合权值筛选准则确定符合条件的PRC分量进行信号，并进行重构。其次，对重构后信号进行SVD滤波降噪，以降低噪声干扰。最后，利用MOMEDA算法增强信号中的冲击成分，并通过Hilbert包络谱分析得到诊断结果。通过滚动轴承实验数据分析，验证所提出方法的有效性。

1 基本理论介绍

1.1 固有时间尺度分解

对任意一段信号Xt(t表示信号Xt的其中某一段的序号值，t≥0)，固有时间尺度分解(ITD)有着类似EMD的功能，可以将其分解成若干个频率由高频至低频的旋转PRC分量和残余分量，该方法的分解步骤描述如下。

1) 定义一段低频基线提取算子L，使得原信号被分解为一段基线信号，同时也使剩余的信号组成PRC分量，即

其中，Lt=LXt代表基线信号，Ht=(1−L)Xt代表PRC旋转分量。

2) 从信号Xt中提取出局部极值Xk(k=1, 2, ···,M)和对应的时刻τk，M为极值点数。定义τ0=0，表示当τk中的k取值为0的时候对应的局部极值点的时间。

3) 在连续两个局部极值点的时间间隔内[ τk,τk+1]，定义基线提取算子

其中

α为控制提取固有旋转分量幅度的增益控制参数，0＜α＜1，通常取0.5。

4) 当第1个PRC分量被分离出来之后，把基线信号视为新的待分解信号，不断重复以上的过程，直到该信号成为单调函数之后，分解完成。信号Xt的分解过程为

1.2 SVD分解

SVD分解是一种特殊的矩阵变换，由于其可以很好地表征信号的本质特征，因而近年来在信号的降噪和周期性分量的提取中都有广泛的应用。设原始信号x={x1,x2,···,xN}，可以表示为xd=Sd+ωd,d=1,2,···,N。其中，Sd表示真实信号，ωd为噪声信号，下标d表示信号Sd其中某一段的序号值，N为数据长度。用上述数据构造如下m×n阶的Hankel矩阵。

式中，S为对应于真实信号在重构相空间内的轨迹矩阵，W为对应于噪声的轨迹矩阵。m=N+n−1，则 Hankel∈Rm×n，R代表全体实数。无论上述矩阵的行和列相关与否，恒存在正交矩阵U与V，使得

则所有的bi形成的序列B=[b1,b2,···,bq−1]称作奇异值的差分谱序列，其表示了相邻奇异值的变化过程。由于真实信号的能量较集中，根据奇异值的差分谱定义，一定含有一个相对较大的峰值bβ，将信号和噪声分隔开，即峰值后的奇异值主要由噪声主导，且相邻奇异值能量差值也必定较小。而峰值前的奇异值主要由真实信号主导，生成的谱峰会十分陡峭。因而可以选取符合上述条件的分界点∂，使得该点的能量峰值较大，在这之后的各点峰值相对较小。即前 ∂个奇异值对应的分量代表的就是真实信号，而噪声为∂ 点之后的奇异值分量。

1.3 多点最优最小熵反褶积

如果x为原始信号，y为冲击信号，h为系统频响函数，e为背景噪声，则采集到的信号可以表示为

MOMEDA方法的实质是利用一种非迭代方法寻找到一个最优的FIR滤波器，通过解卷积过程由x尽可能恢复振动冲击信号y，并且降低噪声的影响。解卷积过程为

式中，k=1,2,···，ς−J，f为滤波器系数，J为FIR滤波器的长度。针对周期性冲击信号的特点，该方法引入了多点D-范数。即

式中，g是用于确定目标冲击成分的位置和权重的常数矢量。通过求解多点D-范数的最大值得到最优滤波器，此时解卷积过程也得到最优解。求式(11)的极值，首先对滤波器系数f=(f1,f2,···,fJ)求导。

因此，式(12)可以写为

经查询相关参考文献，Gk及λ0仅表示一普通矩阵。令式(15)中的 λ0=[G1,G2,···,Gk]，则式(12)等价于

对式(16)进行整理，可得

由于y=λ0Tf，并假设( λ0λ0T)−1存在，则得到

式(18)的特解就是一组最优滤波器，即

将式(19)代入y=λ0Tf中，就可以重构出原冲击信号y。

2 有效分量的筛选

经过ITD方法自适应分解得到的分量信号虽然降低了模态混叠和端点效应，但可能存在一定的虚假分量或与原信号关联度较小的分量。为避免盲目选取和舍弃分量信号可能产生的故障特征信息提取不充分的问题，本文引入了基于K-C组合权重的分量信号筛选标准。

2.1 峭度

峭度可以衡量信号波形的尖峰程度，能够较好地感知信号中的冲击成分。对于冲击成分比重越大的信号，其峭度值也会越高。其计算方式为

其中，θi和x¯分别是原信号的具体值与平均值。ℓ为信号的标准差，µ为样本个数。

2.2 相关系数

根据互相关准则的定义可知，其取值越大表示与原信号的关联度越高。通过利用该准则，可以判断出采用时频分析法得到的分量信号与原信号存在的内在关系。其计算方式为

其中，δi和y¯ 分别是信号ϑ 的具体值与平均值。

2.3 K-C组合权重

虽然峭度对冲击成分较为敏感，但有用信号在噪声的干扰下，极容易受到影响。相关系数虽然可用作衡量两个信号间关联程度的指标，但其容易受到样本总数的影响。通过综合以上两个评价指标的优缺点，得到计算方法为

其中，ϕ 和φ分别对应于峭度和相关系数的权值，且满足 ϕ+φ=1。结合文献[15]的参数选取规则，本文经过多次实验，选取权值ϕ = 0.7，φ=0.3。

3 故障特征提取流程

本文提出一种基于ITD-SVD和MOMEDA相结合的故障特征提取方法。故障特征提取流程如图1所示。首先，求出经过ITD分解得到的各个信号分量的相关系数和峭度值，并计算出K-C组合权重，再以此为依据筛选符合条件的PRC分量进行信号重构。其次，通过SVD对重构信号进行降噪，以消除噪声干扰。为增强降噪后的轴承故障信号中的冲击成分，引入MOMEDA方法对信号进行滤波处理，最后求出滤波之后的信号的Hilbert包络谱，提取相应的故障频率，从而对故障进行准确诊断。其具体步骤如图1。

图1 基于ITD-SVD和MOMEDA的故障特征提取流程Figure1 Fault feature extraction process based on ITD-SVD and MOMEDA

1)选取轴承故障数据，并通过ITD将振动信号分解为若干个PRC分量。

2)根据峭度和互相关系数计算法构建K-C组合权重筛选体系，并进行有效分量的筛选。

3)将筛选得到的PRC分量进行重构，并对重构后的信号展开SVD降噪。

4)采用MOMEDA算法增强降噪后的信号中的周期性冲击成分。

5)对上一步得到的信号进行Hilbert包络解调，提取出故障频率并进行故障诊断。

4 实验分析

为验证所提出方法的有效性，采用凯斯西储大学(CWRU)的轴承数据展开实验[15]。故障检测实验台如图2所示。驱动端的轴承型号为6205-2RSJEMSKF，电机转速为1797r/min(转动频率为1797/60Hz=29.95 Hz)。采样频率为12000Hz，实验数据长度为2048。将损伤直径为0.007in(0.01778cm)的内圈振动信号用于实验分析。按文献[16]中的故障特征频率计算方法可得到对应的故障特征频率，结果如表1所示。

图2 故障检测实验台Figure2 Test bench of fault detection

表1 故障特征频率Table1 Fault characteristic frequency Hz

根据K-C组合权重的计算公式，计算得到轴承正常运行信号的K-C=1.752，详细结果如表2所示。基于该指标，对后续故障信号的分量进行筛选。

表2 正常信号K-C权重指标Table2 K-C weight index of normal signal

轴承内圈故障信号的时域和频域如图3和图4所示。虽然受到一定程度的噪声干扰，但仍旧可以看出部分周期性的脉冲。为检验所提出方法在强噪声干扰下的有效性，本文在原始内圈故障信号中添加信噪比为-10 dB的高斯白噪声，所生成的故障混合信号的时域和频域如图5和图6所示。由于噪声干扰的影响，从图5中已经无法分辨出周期性的脉冲特征，难以直接从时域和频域中提取出故障特征频率。

图3 内圈故障时域波形Figure3 Time domain waveform of inner ring fault

图4 内圈故障频域波形Figure4 Frequency domain waveform of inner ring fault

图5 添加噪声后的内圈故障时域波形Figure5 Time domain waveform of inner ring fault after adding noise

图6 添加噪声后的内圈故障频域波形Figure6 Frequency domain waveform of inner ring fault after adding noise

接下来，采用ITD算法分解添加过噪声的轴承内圈故障信号。信号经过分解后，可得到5个PRC分量和1个残余分量，结果如图7所示。为突出本文所提方法的优越性，实验同时采用传统的EMD算法对故障信号进行分析。经过EMD算法分解后，得到10个IMF分量和1个残余分量，结果如图8所示。通过比较图7和图8可以看出，经过EMD分解得到了更多的信号分量，且上述信号分量存在一定程度的模态混叠，同时产生了很多虚假分量。其次，分别计算基于ITD和EMD分解得到的各个信号分量的互相关系数和峭度，再根据K-C组合权重的计算公式求出对应指标值，计算结果如表3和表4所示。经过比较表3和表4中的数据可知，PRC2和PRC3信号分量，以及IMF1、IMF2和IMF4信号分量均大于所设定的阈值，说明其与原始信号的关联度较高，同时也保留了原始信号中较多的冲击特征。因此，选取上述信号分量进行重构。

表3 PRC分量的K-C组合权重系数Table3 K-C combination weight coefficient of PRC component

表4 IMF分量的K-C组合权重系数Table4 K-C combination weight coefficient of IMF component

图7 ITD分解结果图Figure7 ITD decomposition results

图8 EMD分解结果图Figure8 EMD decomposition results

通过对基于ITD算法和K-C组合权重筛选得到的分量信号进行重构后，再进一步构建Hankel矩阵，并进行SVD分解。所得到的奇异值谱图与奇异值能量差分谱图分别如图9和图10所示。为便于观察，实验所作的图只选取了前100个点。从图10中可看出信号第9个峰值点的能量相对较大，且在其之后的所有峰值点与之相比相对较小，所以选择此峰值对应的重构阶数(此处为18)。

图9 信号奇异值分布曲线Figure9 Distribution curve of signal singular value

图10 奇异值能量差分谱Figure10 Singular value energy difference spectrum

图11和图12分别是基于ITD-SVD方法和基于EMD-SVD方法降噪后的信号时域波形图。经过SVD降噪，噪声干扰均得到不同程度的降低。为定量比较两种方法的降噪效果，选取均方误差(MSE)、平均绝对误差(MAE)、均方根误差(RMSE)作为评价指标，计算结果如表5所示。从中可知，基于ITDSVD方法得到的均方误差、平均绝对误差、均方根误差指标均优于基于EMD-SVD方法得到的结果。为进一步作对比分析，本文对未经过降噪的内圈故障混合信号进行包络谱分析，结果如图13所示。将基于上述两种方法处理得到的信号分别进行包络谱分析，结果如图14和图15所示。

图11 基于ITD-SVD降噪后的信号时域波形Figure11 Time domain waveform of signal after noise reduction based on ITD-SVD

图12 基于EMD-SVD降噪后的信号时域波形Figure12 Time domain waveform of signal denoising based on EMD-SVD

表5 ITD和EMD分解法的性能比较Table5 Performance comparison of ITD and EMD decomposition method

图13 原始含噪内圈故障信号包络谱图Figure13 Envelope spectrum of fault signal of inner ring with noise

通过分析图13可知，由于受到复杂噪声的干扰，轴承内圈故障信号的包络谱中无关频率成分对故障诊断的干扰较大，难以分辨出故障频率。通过分析图14和图15可知，经过基于ITD-SVD方法和基于EMD-SVD方法降噪，噪声对有效信号的影响均有着不同程度的降低。同时，基于ITD-SVD方法降噪后的信号包络谱图中已经可以提取到轴承内圈的故障基频(164.1 Hz)和2倍频(328.2 Hz)，且上述频率很接近于内圈故障特征频率的理论计算值162.185 2 Hz以及基频的2倍频。但是，从基于EMD-SVD方法降噪后的包络谱中仍难以提取轴承内圈的故障基频。此外，从图14中也可看出，故障信号的包络谱中出现了部分无关频率成分，这对故障诊断产生了一定程度的干扰。因此，还需通过MOMEDA滤波方法进一步增强降噪后的信号中的周期性冲击成分。

图14 基于ITD-SVD降噪后的信号包络谱图Figure14 Signal envelope spectrum after noise reduction based on ITD-SVD

图15 基于EMD-SVD降噪后的信号包络谱图Figure15 Signal envelope spectrum after noise reduction based on EMD-SVD

通过MOMEDA滤波和Hilbert包络解调后，得到的结果如图16所示。从图16中的包络谱分析结果来看，经过MOMEDA滤波后，已经有效地将降噪之后的信号中的周期性故障冲击成分提取出来，出现了更大幅度的峰值，并弱化了与故障冲击无关的信号成分。图中的故障特征频率更加明显，不但可以明显分辨出轴承内圈故障的基频，而且二倍频(328.1 Hz)、三倍频(492.2 Hz)、四倍频(656.3 Hz)等多倍频均清晰可见，且与理论计算得到的故障特征频率及倍频很接近。因此，根据特征频率可以准确诊断出轴承内圈发生了故障。