CN103810394A - 旋转设备故障信号奇异值分解降噪的设计方法 - Google Patents

Abstract

本发明是一种旋转设备故障信号奇异值分解降噪的设计方法，其特点是，它包括SVD降噪、重构矩阵的维数和有效秩阶数的确定等内容，通过选择第一个至少单边与其相邻峰值比较，差距绝对值最大的极大峰值的对应点位置，来确定重构信号的有效秩阶数，从而完成对有用信号的重构和对噪声的有效消除。具有既能够直观有效地确定奇异值有效秩降噪阶数，又能够降低算法的复杂程度，降噪效果好，信号的信噪比高等优点。

Description

旋转设备故障信号奇异值分解降噪的设计方法

技术领域

本发明涉及旋转设备的故障诊断研究领域，是一种旋转设备故障信号奇异值分解降噪的设计方法。

背景技术

旋转设备是电厂的重要机械设备，包括汽轮机、风机、水泵、电动机等，电厂的正常生产运行离不开旋转设备，其安全稳定的运行是电厂高效、经济运行的关键。旋转机械发生故障时总是会伴随一定的状态的变化，异常状态有很多，其中振动和噪声是最为常见的。而生产设备由于在系统中相互关联，一旦某一部件发生故障就会导致连锁反应，不仅使设备本身受到损坏，失去了预定的功能，而且会严重的影响电厂的生产过程，不但使电厂在综合经济效益上蒙受巨大的损失，设备受到不可修复的损坏，甚至会发生威胁现场工作人员人身安全的重大安全事故。

有效消除噪声影响一直是旋转设备故障诊断研究的重要内容之一，特别是故障早期，由于调制源弱，早期故障信号微弱，且受周围设备的噪声干扰，导致故障特征难以识别。因此，最大限度地提高振动测量信号的信噪比，是为故障特征信号的提取做好前期工作的重要环节。振动信号Hankel矩阵奇异值分解，作为一种非线性滤波方法，可以用于消除信号中的随机噪声成分，提取信号中的周期成分，得到相对纯净的故障信号。

如何确定奇异值有效秩阶次是该方法的关键技术之一，目前采用方法相对较多的是试凑法和均阀法，然而这两种方法对操作者的经验要求相对较高，且不易掌握。本发明将基于Hankel矩阵和奇异值差分谱单边极大值分解技术，研究奇异值分解降噪特性消除系统信号中混合噪声的改进设计方法。

发明内容

本发明的目的是针对旋转设备故障信号的噪声干扰问题，克服现有技术奇异值分解降噪有效秩单纯的选取差分谱中最大值峰值所对应点的缺陷，既能够直观有效地确定奇异值有效秩降噪阶数，而且又能够降低算法的复杂程度，降噪效果好，信号的信噪比高的旋转设备故障信号奇异值分解降噪的设计方法。

实现本发明目的所采用的技术方案是：一种旋转设备故障信号奇异值分解降噪的设计方法，其特征是，它包括以下内容：

1）SVD降噪：

设从滚动轴承测得的含有噪声的数据信号为y=[y₁,y₂，…，y_N]，基于相空间重构理论，将上述数据构造成p×q阶Hankel矩阵:

$H_m=\left[\begin{array}{cccc}{y}_1& y_2& ...& y_q\\ y_2& y_3& ...& y_{q+1}\\ ...& ...& ...& ...\\ y_p& y_{p+1}& ...& y_N\end{array}\right]---\left(1\right)$

H_m为p×q阶矩阵，其中N为信号长度，N=p+q-1并且p≥q；

矩阵H_m通过重构吸引子的特征揭示了它在重构空间的动态特性，因此，将H_m表示为H_m=D+W的形式，D表示光滑信号在重构空间的p×q矩阵，W表示噪声干扰信号的p×q矩阵。所以如何对原始信号进行降噪，就是怎样寻找到D的最佳逼近矩阵；

对H_m进行奇异值分解可以得到：

H_m=USV^T      (2)

其中U和V^T分别为p×p和q×q矩阵，S为p×q的对角矩阵，主对角线元素为λ_i(1,2,…,k)，k=min(p,q)即：

S=diag(λ₁,λ₂,…,λ_k)     (3)

式中λ₁,λ₂,…,λ_k为矩阵H_m的奇异值，且λ₁≥λ₂≥…≥λ_k≥0，U和V^T表示左右奇异阵；根据奇异值分解理论和Frobeious范数意义下矩阵最佳逼近定理得到：有用的信号主要由前r个较大的奇异值反映，噪声信号由后面较小的奇异值反映，去掉代表噪声信号的较小奇异值，则源信号中的噪声被去除，再进行奇异值分解的逆过程演算，最终得到矩阵

那么矩阵

就是H_m的秩为r的最佳逼近矩。

相对于H_m而言其噪声被压缩，将矩阵

中的反对角线元素相加平均，就得到降噪后的信号；

2）重构矩阵的维数和有效秩阶数的确定

（1）确定重构矩阵的维数

Hankel重构矩阵维数的确定，直接影响降噪效果的好坏，通过对多组不同长度及频率的源信号进行分析后发现，最佳维数基本在p=N/2处的一个邻域内产生，并且在此邻域所取的维数值的降噪效果最佳，能满足工程要求，因此，重构矩阵的结构根据振动信号长度N来确定，工程应用中取p=N/2；

（2）确定有效秩阶数

奇异值的差分谱序列b_i=λ_i-λ_i+1i=1,2,…,q-1能够描述奇异值序列的具体变化，差分后形成序列B=(b₁,b₂,…,b_q-1)充分反应了相邻两个奇异值的变化，在奇异值差分谱中从右至左，选择第一个至少单边与其相邻峰值比较，差距绝对值最大的极大峰值的对应点位置，来确定重构信号的有效秩阶数，从而完成对有用信号的重构和对噪声的有效消除。

本发明的旋转设备故障信号奇异值分解降噪的设计方法，提供了一个有效确定奇异值秩阶次的选取原则，使得操作者能够有效地、快速地确定选取阶次，而不用在单纯依靠操作者的经验来进行选取，提高了可操作性和实用性，避免了传统差分方法只取最大值点的局限性，采用单边极大值原则使得选取方法有效避免了过降噪现象发生，为后续故障诊断提供了有利的保证。具有既能够直观有效地确定奇异值有效秩降噪阶数，而且又能够降低算法的复杂程度，降噪效果好，信号的信噪比高等优点。

附图说明

图1为旋转设备中的一种送风机滚动轴承结构示意图；

图2为x(t)₁加噪后的信号示意图；

图3为x(t)₁奇异值分布曲线示意图；

图4为x(t)₁奇异值差分谱示意图；

图5为选取有效秩阶数为2的x(t)₁降噪图；

图6为选取有效秩阶数为4的x(t)₁降噪图；

图7为选取有效秩阶数为6的x(t)₁降噪图；

图8为x(t)₂加噪后的信号示意图；

图9为x(t)₂奇异值分布曲线示意图；

图10为x(t)₂奇异值差分谱示意图；

图11为选取有效秩阶数为4的x(t)₂降噪图；

图12为选取有效秩阶数为8的x(t)₂降噪图；

图13为选取有效秩阶数为10的x(t)₂降噪图；

图14为图1的轴承外圈故障信号奇异值分布曲线示意图；

图15为图1的轴承外圈故障信号奇异值差分谱示意图；

图16为图1的轴承外圈故障信号原始信号图；

图17为图1的轴承外圈故障信号选取有效秩阶数为44的降噪图；

图18为图1的轴承外圈故障信号选取有效秩阶数为4的降噪图；

图19为图1的轴承外圈故障信号选取有效秩为44的降噪后的包络谱示意图；

图20为图1的轴承外圈故障信号选取有效秩为4的降噪后的包络谱示意图；

图中的标号、符号和线条等说明如下：

图1中，β为接触角；D为滚道节径；d为滚动体直径。

图2-图7为信号x(t)₁见仿真分析的降噪结果，横坐标表示采样点数，图2纵坐标表示幅值，单位为mv；图3纵坐标表示奇异值；图4纵坐标表示奇异值差分值，图5、图6、图7纵坐标表示幅值，单位为mv。

图8-图13为信号x(t)₂见仿真分析的降噪结果，横坐标表示采样点数，图8纵坐标表示幅值，单位为mv；图9纵坐标表示奇异值；图10纵坐标表示奇异值差分值；图11、图12、图13纵坐标表示幅值，单位为mv。

图14-图18为轴承外圈故障信号降噪结果，图14横坐标表示奇异值阶数，纵坐标表示奇异值；图15横坐标表示奇异值阶数，纵坐标表示奇异值差分值；图16、图17、图18横坐标表示时间，纵坐标表示幅值，单位为m·s^-2。

图19、图20为轴承外圈故障信号降噪后包络谱，横坐标表示采样点数，纵坐标表示幅值，单位为m·s^-2。

具体实施方式

下面利用附图和实施例对发明作进一步说明。

本发明的旋转设备故障信号奇异值分解降噪的设计方法，以旋转设备之一的送风机为例：送风机滚动轴承结构如图1所示，其设计的详细过程如下：

1）SVD降噪

设从滚动轴承测得的含有噪声的数据信号为y=[y₁,y₂，…，y_N]，基于相空间重构理论，将上述数据构造成p×q阶Hankel矩阵:

$H_m=\left[\begin{array}{cccc}{y}_1& y_2& ...& y_q\\ y_2& y_3& ...& y_{q+1}\\ ...& ...& ...& ...\\ y_p& y_{p+1}& ...& y_N\end{array}\right]---\left(4\right)$

H_m为p×q阶矩阵，其中N为信号长度，N=p+q-1并且p≥q。

矩阵H_m通过重构吸引子的特征揭示了它在重构空间的动态特性，因此，将H_m表示为H_m=D+W的形式，D表示光滑信号在重构空间的p×q矩阵，W表示噪声干扰信号的p×q矩阵。所以如何对原始信号进行降噪，就是怎样寻找到D的最佳逼近矩阵。

对H_m进行奇异值分解可以得到：

H_m=USV^T      (5)

其中U和V^T分别为p×p和q×q矩阵，S为p×q的对角矩阵，主对角线元素为λ_i(1,2,…,k)，k=min(p,q)即：

S=diag(λ₁,λ₂,…,λ_k)      (6)

式中λ₁,λ₂,…,λ_k为矩阵H_m的奇异值，且λ₁≥λ₂≥…≥λ_k≥0，U和V^T表示左右奇异阵。根据奇异值分解理论和Frobeious范数意义下矩阵最佳逼近定理得到：有用的信号主要由前r个较大的奇异值反映，噪声信号由后面较小的奇异值反映，去掉代表噪声信号的较小奇异值，则源信号中的噪声被去除，再进行奇异值分解的逆过程演算，最终得到矩阵

那么矩阵

就是H_m的秩为r的最佳逼近矩。

相对于H_m而言其噪声已被大大压缩。将矩阵

中的反对角线元素相加平均，就可以得到降噪后的信号。

2）重构矩阵的维数和有效秩阶数的确定

（1）确定重构矩阵的维数

Hankel重构矩阵维数的确定，直接影响降噪效果的好坏。通过对多组不同长度及频率的源信号进行分析后发现，最佳维数基本在p=N/2处的一个邻域内产生，并且在此邻域所取的维数值的降噪效果比较理想，而且能满足工程要求。因此，重构矩阵的结构可以根据振动信号长度N来确定，工程应用中可以取p=N/2。经验证该方法可行且效果较好。

（2）确定有效秩阶数

奇异值的差分谱序列b_i=λ_i-λ_i+1i=1,2,…,q-1描述了奇异值序列的具体变化情况。差分后形成序列B=(b₁,b₂,…,b_q-1)充分的反应了相邻的两个奇异值的变化。根据定义的差分谱序列可知，两个相邻的奇异值差别越大，那么他们在整个差分谱中所表现出的特征也越明显，也就是说它们之间产生的峰值也就相对越大。而这些较大的峰值通常携带有非常重要的状态信息，因此较大峰值尤其应受到关注，这其中也包含着最大峰值，但并不意味着最大峰值就是最值得关注的信息。而这些较大峰值的出现就是由于有用信号和噪声信号的不相关而导致的。因此，根据所得峰值的不同，选则较恰当的奇异值有效秩阶次就能够去除噪声信号从而得到有用信号。采用单边极大值原则，即在奇异值差分谱中从右至左，选择第一个至少单边与其相邻峰值比较，差距绝对值最大的极大峰值的对应点位置，来确定重构信号的有效秩阶数，从而完成对有用信号的重构和对噪声的有效消除。

具体实施例：

1、仿真验证及分析

为验证本发明旋转设备故障信号奇异值分解降噪的设计方法的有效性，分别用不同频率的信号源加载不同能量分贝的噪声对该方法进行验证。现选取x(t)₁和x(t)₂两个信号进行分析，并通过对比降噪后信号和加噪前信号的图谱接近程度和信噪比的大小来衡量该方法的高效实用性。

x(t)₁=(1+0.6t)sin(80πt)+2cos(30πt)+ζ(n)

x(t)₂==5sin(3πt)+4sin(8πt)cos(2πt+3π)+8sin(8πt)+ζ(n)

其中：ζ(n)分别表示强度为5dB和15dB的高斯白噪声。

以上两个信号源采样长度均取1024，Hankel矩阵行数取512。跟据奇异值差分谱我们能够看到x(t)₁较大峰值处大奇异值分别为2和4，其中最大峰值处为4，x(t)₂较大峰值处大奇异值分别是2、4和8，其中最大峰值处为2，根据本发明提到的奇异值差分谱单边极大值原则，所以两个信号源有效秩阶数分别确定为4和8，降噪后分别对应于图5和图11。为了进行比较，对信号x(t)₁的奇异值2、6和x(t)₂的奇异值6、10也做了相应处理，对应于图6和图13，降噪前后的SNR如表1所示。

表1信号x(t)₁降噪前后SNR

表2信号x(t)₂降噪前后SNR

从表1与表2中可以看出，噪声信号通过奇异值差分谱单边极大值原则降噪，SNR得到大大提高，而图5和图11是选取有效秩阶数为4和8的信噪比，为最大。x(t)₁信号图2降噪前的SNR之所以出现负值，是由于源信号的能量小于噪声能量的原因，即源信号被噪声信号所覆盖所造成。为了能更清晰的观察，将信号x(t)₁降噪图和x(t)₂降噪图中奇异值差分谱值分别取前100个奇异值进行放大，从中可以看到差分谱有效地反映了所包含大信息量的奇异值个数，并且根据单边极大值原则选择了有效秩阶数分别为4和8进行降噪，而不是选取了最大峰值所对应的有效秩阶数。

x(t)₁选择有效秩阶数为4，是因为4的峰值与其相邻峰值的2比较差距最大，而x(t)₂之所以选择8是因为遵循本发明所提原则，从右至左8与4所对应峰值的差距最大，虽然比较后为负值但是原则提到绝对值最大，且从右至左8点为第一个符合要求的点，所以选择8为x(t)₂有效秩阶数。降噪后的信号图谱比较光滑并且与加噪之前的信号图谱对比基本吻合。

图5为有效秩阶数为2的降噪图形，可以看出降噪后的图形明显失真，主要原因是因为有效秩阶数的选取偏低，造成过降噪后果，图7为有效秩阶数为6的降噪图形，由于有效秩阶数选取偏高导致降噪后的效果存在大量毛刺，说明噪声没有完全去掉，仍有保留。图11有效秩阶数为4、图13有效秩阶数为10，也是由于有效秩的选取不恰当出现了欠降噪和过降噪的结果。

2、实验验证及分析

如示意图1所示，本发明研究对象为送风机滚动轴型号N205的轴承外圈故障信号，滚动体个数Z=12，轴承节圆直径D=39mm，滚动体直径d=7.5mm，压力角β=0。

信号采集频率为f=12800Hz，转速为1440r/min，采样点个数取3000，实验结果如图14-图20所示。

从图14-图18中可以看出，根据本发明提出的单边极大值原则选择有效秩阶数为44，而没有选择最大峰值阶数4。从图17和图16的比较可以看出，降噪后的信号不仅有效的保留了源信号的周期部分，而且大大降低了源信号的幅值，有效地去除了噪声干扰信号，为后期的故障特征信号提取做了充分准备。本发明同时作了经验法经常容易选取的有效秩为4的降噪效果图，如图18所示，幅值同样降低很多，但是为了进一步表征轴承故障的物理意义，对图17和图18两组信号做了如下处理。

根据轴承外圈故障特征频率公式可得：

$f_o=0.5\mathrm{Zf}(1-\frac{d}{D}\cos \mathrm{\β})=0.5\·12\·24\·(1-\frac{7.5}{39})=116.2\mathrm{Hz}---\left(7\right)$

为了提取外圈故障特征频率，进一步对降噪后的信号做Hilbert包络并进行频率谱分析，结果如图19-图20所示。

从图19功率谱的分析可以发现存在一个明显的峰值116.4Hz，对应轴承外圈故障特征频率，同时存在232.8Hz、348.5Hz和464.9Hz频率成分，这刚好与故障特征频率的二倍频、三倍频、四倍频相对应，跟轴承外圈故障特征相吻合。而图20中不但没有显示出明确的故障频率信息，而且功率谱图像严重失真，缺少大量信息。故此可以分析，该谱是由于降噪的有效秩阶数过高，从而导致严重过降噪，失去有用信号所造成。所以，通过对比试验，进一步的验证了本发明旋转设备故障信号奇异值分解降噪的设计方法的有效性。

具体实施方式仅是对本发明的说明，并不构成对权利要求保护范围的限制，本领域技术人员不经过创造性劳动的等同替代，均在本发明保护范围内。

Claims (1)

1.一种旋转设备故障信号奇异值分解降噪的设计方法，其特征是，包括以下内容：

1）SVD降噪：

设从滚动轴承测得的含有噪声的数据信号为y=[y₁,y₂，…，y_N]，基于相空间重构理论，将上述数据构造成p×q阶Hankel矩阵:

$H_m=\left[\begin{array}{cccc}{y}_1& y_2& ...& y_q\\ y_2& y_3& ...& y_{q+1}\\ ...& ...& ...& ...\\ y_p& y_{p+1}& ...& y_N\end{array}\right]---\left(1\right)$ H_m为p×q阶矩阵，其中N为信号长度，N=p+q-1并且p≥q；

矩阵H_m通过重构吸引子的特征揭示了它在重构空间的动态特性，因此，将H_m表示为H_m=D+W的形式，D表示光滑信号在重构空间的p×q矩阵，W表示噪声干扰信号的p×q矩阵。所以如何对原始信号进行降噪，就是怎样寻找到D的最佳逼近矩阵；

对H_m进行奇异值分解可以得到：

H_m=USV^T      (2)

其中U和V^T分别为p×p和q×q矩阵，S为p×q的对角矩阵，主对角线元素为λ_i(1,2,…,k)，k=min(p,q)即：

S=diag(λ₁,λ₂,…,λ_k)     (3)

式中λ₁,λ₂,…,λ_k为矩阵H_m的奇异值，且λ₁≥λ₂≥…≥λ_k≥0，U和V^T表示左右奇异阵；根据奇异值分解理论和Frobeious范数意义下矩阵最佳逼近定理得到：有用的信号主要由前r个较大的奇异值反映，噪声信号由后面较小的奇异值反映，去掉代表噪声信号的较小奇异值，则源信号中的噪声被去除，再进行奇异值分解的逆过程演算，最终得到矩阵

那么矩阵

就是H_m的秩为r的最佳逼近矩。

相对于H_m而言其噪声被压缩，将矩阵中的反对角线元素相加平均，就得到降噪后的信号；

2）重构矩阵的维数和有效秩阶数的确定

（1）确定重构矩阵的维数

Hankel重构矩阵维数的确定，直接影响降噪效果的好坏，通过对多组不同长度及频率的源信号进行分析后发现，最佳维数基本在p=N/2处的一个邻域内产生，并且在此邻域所取的维数值的降噪效果最佳，能满足工程要求，因此，重构矩阵的结构根据振动信号长度N来确定，工程应用中取p=N/2；

（2）确定有效秩阶数

奇异值的差分谱序列b_i=λ_i-λ_i+1i=1,2,…,q-1能够描述奇异值序列的具体变化，差分后形成序列B=(b₁,b₂,…,b_q-1)充分反应了相邻两个奇异值的变化，在奇异值差分谱中从右至左，选择第一个至少单边与其相邻峰值比较，差距绝对值最大的极大峰值的对应点位置，来确定重构信号的有效秩阶数，从而完成对有用信号的重构和对噪声的有效消除。

Images