一种基于奇异谱分解的旋转机械故障诊断方法与流程

本发明属于机械工程技术领域，具体涉及一种旋转机械故障诊断方法。

背景技术：

工业现场中常见的零部件，如轴承、齿轮、转子等，是旋转机械设备中的重要组成构件，由于在长期高速运转、交变载荷等恶劣工况下，这些关键部件极易产生局部损伤并演化成晚期故障，极大地影响着整个旋转机械传递系统的工作性能。因此，针对旋转机械设备中的关键部件进行有效故障检测具有重要的现实意义。由于实际工程中旋转机械振动信号表现为非平稳特性，其受到部件间多传递路径耦合作用的影响，振动信号中掺杂着强烈背景噪声和多干扰源信号，致使旋转机械的故障特征提取极其复杂和困难。因此，如何从振动信号中有效地提取有用故障特征信息，是判断旋转机械设备是否发生故障需要解决的关键难题。

现有的旋转机械故障诊断方法主要包括：小波变换、经验模式分解和奇异谱分析。小波变换受海森堡测不准原理的约束，需选择基函数且缺乏自适应性，对非平稳性信号分析存在局限性。经验模式分解具备自适应性，能够对非平稳信号进行有效分析，但其采用插值拟合包络原理进行分解具备较大拟合误差，同时存在过包络、欠包络、端点效应和模态混叠等缺陷。奇异谱分析是一种基于主成分分析的非参数谱估计方法，包括嵌入维数的选取、奇异值分解及分量序列的重构等过程，能有效用于旋转机械测量信号的降噪处理，但其主要依赖经验性选取嵌入维数长度构建轨迹矩阵，并基于轨迹矩阵的主成分进行奇异值分解、特征重组及信号重构，很大程度上影响其诊断结果。

技术实现要素：

针对上述技术存在的问题，本发明提供了一种基于奇异谱分解的旋转机械故障诊断方法，旨在提高旋转机械故障诊断的精确度。

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

一种基于奇异谱分解的旋转机械故障诊断方法，包括以下步骤：

步骤1，在旋转机械的关键部件附近安装传感器进行测量，采集测量信号作为源信号；

步骤2，利用奇异谱分解对源信号进行分解，得到若干个瞬时频率具有物理意义的奇异谱分量；

步骤3，根据特征能量因子最大准则挑选蕴含丰富故障特征信息的分解分量作为主奇异谱分量；

步骤4，对主奇异谱分量进行希尔伯特解调，获得相应的包络谱；

步骤5，从包络谱中观察故障特征频率处是否存在明显峰值，从而实现旋转机械故障类型的准确判别。

步骤1中，所述的旋转机械的关键部件为轴承、齿轮和转子。

所述步骤2中利用奇异谱分解对源信号进行分解的步骤为

步骤2.1，对源信号构建一个新的轨迹矩阵表达式：以源信号x(n)＝{1,2,3,4,5}为例，其嵌入维数M为3，则对应的矩阵X为

式中，矩阵X的虚线部分对应于奇异谱分析中的新轨迹矩阵，且新矩阵中各对角线元素相同且数量相等；

步骤2.2，自适应选择第j次迭代所需的嵌入维数M：首先计算第j次迭代时残余分量的功率谱密度PSD，即

式中，j＝2,3,4...，υj(n)表示第j次迭代时获得的残余分量，υk(n)表示第k次迭代时获得的残余分量；

随后估计功率谱密度PSD中最大峰值对应的频率fmax；在第1次迭代中，如果归一化频率fmax/Fs小于给定阈值10^-3，则残余量被视为一个大趋势项，此时M设置为N/3，其中Fs是采样频率；否则，当迭代次数j＞1，嵌入维数被设置为M＝1.2×(Fs/fmax)；

步骤2.3，按照高频到低频的顺序依次重构第j个奇异谱分量：在第1次迭代中，如果检测到一个大趋势项，只使用第1个左右特征向量去获取g⁽¹⁾(n)，使得且g⁽¹⁾(n)从X1的对角平均中获得；否则，当迭代次数j＞1，必须获得一个分量序列g⁽¹⁾(n)去描述一个物理意义明确的时间尺度；在这个意义上，它的频率成分主要集中在频带[fmax-δf,fmax+δf]之间，其中δf表示残余项功率谱密度中主峰的半带宽；因此，根据左特征向量在频谱[fmax-δf,fmax+δf]范围中具有突出主频率的所有特征组和对选取模态分量的主峰能量贡献最大的一个特征组，创建一个子集Ij(Ij＝{i1,i2,…,i3})；随后通过矩阵XIj＝Xi1+…+Xip的对角平均法重构这对应的奇异谱分量；

步骤2.4，设置迭代终止条件：将迭代估计出的奇异谱分量从源信号中分离出，得到一个残余项计算所得残余项和源信号之间的归一化均方差NMSE，即

式中，表示第j次迭代时获得的奇异谱分量，υ^(j+1)(n)表示第j次迭代时获得的残余项，υ^(j)(n)表示第j次迭代时设定的源信号；

当归一化均方差NMSE小于给定阈值th＝1％时，整个分解过程终止，否则，将残余项作为输入信号重复上述迭代过程获得最终分解结果，即

式中m是获得的奇异谱分量个数，在每次迭代后残余项υ^(j+1)(n)的能量都会减小。

所述步骤3中根据特征能量因子最大准则选取蕴含最丰富故障特征信息的主奇异谱分量的步骤为：

步骤3.1，首先将源信号通过奇异谱分解得到若干个瞬时频率具有物理意义的奇异谱分量；

步骤3.2，然后计算若干个奇异谱分量的特征能量因子λ，其计算公式为

λ＝E/E^*

式中，E是旋转机械故障特征频率前m倍频处的幅值能量和，E*是源信号的包络谱幅值总能量；

步骤3.3，最后将最大特征能量因子λmax对应的奇异谱分量视为蕴含最丰富故障特征信息的主奇异谱分量。

步骤4的具体步骤为：

步骤4.1，对选取的主奇异谱分量进行希尔伯特变换处理得到希尔伯特变换对，即

式中x(t)为选取的主奇异谱分量；

步骤4.2，以主奇异谱分量为实部、希尔伯特变换对为虚部构成解析信号，即

z(t)＝x(t)+jx(t)

步骤4.3，将解析信号进行求模运算，得到主奇异谱分量的包络信号为

步骤4.4，对包络信号作傅里叶变换得到主奇异谱分量相对应的包络谱。

步骤5中，如果存在突出故障特征频率，则诊断旋转机械发生故障，否则旋转机械运行正常。

有益效果：本发明的技术效果在于：

1)本发明运用奇异谱分解对旋转机械的测量信号进行分解，能自适应地获得若干个瞬时频率具有物理意义的奇异谱分量，克服了奇异谱分析按经验性选取嵌入维数长度的缺陷，无需小波变换基函数，避开了经验模式分解采用插值拟合包络带来的拟合误差，对测量信号中的背景噪声及干扰信号具有较好的分离能力，可初步提取测量信号的冲击波形。

2)通过计算奇异谱分量的特征能量因子选取蕴含丰富故障特征信息的主奇异谱分量进行解调分析，能有效避免工程实际中因人为主观因素选取奇异谱分量引起的错误诊断或无法诊断故障的现象，从而提高故障诊断的精确度。

3)对挑选的主奇异谱分量进行希尔伯特解调，获得相应的包络谱，能够提取明显的故障特征频率，从而实现旋转机械故障类型的准确判别。此外，本发明简单易行，适用于旋转机械的在线监测与故障诊断。

附图说明

图1是本发明的流程图；

图2是轴承外圈故障的振动加速度信号的时域波形图；

图3是经奇异谱分解获得的若干个奇异谱分量的时域波形图；

图4是奇异谱分量的特征能量因子的分布图；

图5是被选取的蕴含丰富故障特征信息的主奇异谱分量的包络谱图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

图1为本发明的流程图。下面结合流程图对本发明的步骤作详细说明。

1)首先将传感器安装在旋转机械的关键部件附近进行测量，采集测量信号作为源信号。其中，关键部件为轴承、齿轮和转子等。

2)对源信号进行奇异谱分解得到若干个瞬时频率具有物理意义的奇异谱分量。其具体过程为：

2.1)首先针对源信号构建一个新的轨迹矩阵。设源信号为x(n)，其数据长度和嵌入维数分别为N和M，将其构建为一个M行N列的矩阵X，矩阵X的第i行为xi＝(x(i),…,x(N),x(1),…,x(i-1))且i＝1,…,M，即矩阵为了更好地理解矩阵构建的过程，以源信号x(n)＝{1,2,3,4,5}为例，其嵌入维数M为3，则对应的矩阵X为

式中，矩阵X的虚线部分左3行左3列对应于奇异谱分析中的轨迹矩阵。为了增强原待分析信号中的振荡成分，并使得迭代后残余分量的能量呈现出递减规律，将矩阵X中的右下角3个元素搬移至矩阵X的左上角位置，构建了一个新的矩阵表达式，即

式中，矩阵X的虚线部分对应于奇异谱分析中的新轨迹矩阵，且新矩阵中各对角线元素相同且数量相等。

2.2)自适应选择嵌入维数大小M。考虑到奇异谱分析按经验选取嵌入维数的缺陷，这里采用一种自适应法则选取第j次迭代使用的嵌入维数大小M。首先计算第j次迭代时残余分量的功率谱密度PSD，即

式中，j＝2,3,4…，υj(n)表示第j次迭代时获得的残余分量，υk(n)表示第k次迭代时获得的残余分量；

随后估计PSD中最大峰值对应的频率fmax。在第1次迭代中，如果归一化频率fmax/Fs小于给定阈值10-³，则残余量被视为一个大趋势项，此时M设置为N/3，其中Fs是采样频率。否则，当迭代次数j＞1，嵌入维数被设置为M＝1.2×(Fs/fmax)。

2.3)按照高频到低频的顺序依次重构第j个奇异谱分量。在第1次迭代中，如果检测到一个大趋势项，只使用第1个左右特征向量去获取g⁽¹⁾(n)，使得且g⁽¹⁾(n)可以从X1的对角平均中获得。否则，当迭代次数j＞1，必须获得一个分量序列g⁽¹⁾(n)去描述一个物理意义明确的时间尺度。在这个意义上，它的频率成分主要集中在频带[fmax-δf,fmax+δf]之间，其中δf表示残余项功率谱密度中主峰的半带宽。因此，根据左特征向量在频谱[fmax-δf,fmax+δf]范围中具有突出主频率的所有特征组和对选取模态分量的主峰能量贡献最大的一个特征组，创建一个子集Ij(Ij＝{i1,i2,…,i3})。随后通过矩阵XIj＝Xi1+…+Xip的对角平均法重构这对应的奇异谱分量。

在这一过程中，为了更好地估计主峰带宽δf，构建了一个具有叠加高斯函数的谱模型用于拟合功率谱密度的轮廓，该模型定义为三个高斯函数之和，且每个高斯函数表示一个谱峰，即

式中Ai是第i个高斯函数的幅值，ui是它的位置，σi是它的带宽和θ＝[Aσ]^T是参数矢量，其中A＝[A1,A2,A3]，σ＝[σ1,σ2,σ3]。第1个高斯函数记录主谱峰对应的频率fmax，第2个高斯函数记录次谱峰对应的频率f2，第3个高斯函数记录前两谱峰间任意峰值对应的频率。因此，依据此模型可得

u1＝fmax，u2＝f2，

通过对模型的加权最小二乘拟合可获得模型参数Ai。首先给定模型拟合的初始参数值，即

随后采用Levenberg-Marquardt法确定模型参数Ai最优值。给定σ1的估计值，则主峰带宽δf的值可由δf＝2.5σ1获得。此外，为了恢复第j个分量序列开始第二次迭代，采用尺度因子调整g^(j)(n)与残余时间序列υ^(j)(n)的差值，即

式中满足

2.4)设置迭代停止条件。将迭代估计出的奇异谱分量从源信号中分离出，得到一个残余项将其作为输入信号重复上述迭代过程，计算所得残余项和原信号之间的归一化均方差NMSE，即

式中，表示第j次迭代时获得的奇异谱分量，υ^(j+1)(n)表示第j次迭代时获得的残余项，υ^(j)(n)表示第j次迭代时设定的源信号；

当归一化均方差NMSE小于给定阈值th＝1％时，整个分解过程终止，获得最终分解结果，即

式中m是获得的奇异谱分量个数，在每次迭代后残余项υ^(j+1)(n)的能量都会减小。

3)根据特征能量因子最大准则筛选出蕴含最丰富故障特征信息的奇异谱分量作为主奇异谱分量。其具体过程为

3.1)计算各奇异谱分量包络解调后的频谱序列为Y，其故障特征频率及其倍频处对应的幅值为Y(i)，其中i＝1,2,…,m，分别是故障特征频率的第i倍频，则故障特征前m倍频处的幅值能量和E为

3.2)计算故障特征前m倍频处的幅值能量和与各奇异谱分量信号包络谱幅值总能量的比值λ，即得到各奇异谱分量的特征能量因子为

λ＝E/E^*，E^*＝∑Y²

式中，E^*是各奇异谱分量的包络谱幅值总能量。

3.3)将最大特征能量因子λmax对应的奇异谱分量视为主奇异谱分量进行特性提取。

4)对选取的主奇异谱分量进行希尔伯特解调，获得相应的包络谱。其具体过程为

4.1)对选取的主奇异谱分量进行希尔伯特变换处理得到希尔伯特变换对，即

式中x(t)为选取的主奇异谱分量。

4.2)以主奇异谱分量为实部、希尔伯特变换对为虚部构成解析信号，即

z(t)＝x(t)+jx(t)

4.3)将解析信号进行求模运算，得到主奇异谱分量的包络信号为

4.4)对包络信号作傅里叶变换得到主奇异谱分量相对应的包络谱。

5)分析包络谱中是否存在具有明显峰值的故障特征频率，如果存在突出故障特征频率，则诊断旋转机械发生故障，否则旋转机械运行正常。

图2为设置有外圈故障的SKF6205型滚动轴承振动信号的时域波形。该振动信号通过安装在轴承座上方的压电式振动传感器进行采集。实验测试过程中，电机转速为1 470r/min，转频fr为24.5Hz，采样频率为12 800Hz，采样点数为12 800。其中外圈故障是通过电火花加工出沟槽来模拟设置的，沟槽宽为0.2mm，沟槽深为1.5mm。经理论公式计算得到轴承外圈故障特征频率fa为88Hz。

图3为轴承外圈故障振动信号经奇异谱分解后得到的奇异谱分量的时域波形图。

图4为奇异谱分量的特征能量因子的分布图，图中可清晰辨识最大特征能量因子对应于奇异谱分量1，这表明奇异谱分量1蕴含最丰富的故障特征信息，从而选取奇异谱分量1作为主奇异谱分量进行后续包络谱分析。

图5为主奇异谱分量的包络谱图，也就是奇异谱分量1的包络谱图，从图中能明显发现88Hz处的外圈故障特征频率fa及多倍频成分，进而判断轴承携带外圈故障，由此说明本文方法在旋转机械滚动轴承故障诊断中的实用性和精确性。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。