基于拟合误差最小化原则的奇异值分解降噪有效秩阶次确定方法

崔伟成， 许爱强， 李 伟， 孟凡磊

(1. 海军航空工程学院 飞行器工程系，山东 烟台 264001； 2. 海军航空工程学院 飞行器检测与应用研究所，山东 烟台 264001)

基于拟合误差最小化原则的奇异值分解降噪有效秩阶次确定方法

崔伟成1， 许爱强2， 李 伟1， 孟凡磊1

(1. 海军航空工程学院 飞行器工程系，山东 烟台 264001； 2. 海军航空工程学院 飞行器检测与应用研究所，山东 烟台 264001)

为了最大限度地提高旋转机械设备故障振动信号的信噪比，研究了奇异值分解降噪的原理，提出了一种新的奇异值分解降噪有效秩阶次的确定方法。首先，对振动信号进行相空间重构，对吸引子轨迹矩阵进行奇异值分解；然后，按不同的阶数，将奇异值分成信号组和噪声组，对每次分组的结果，以阶数为自变量、以奇异值为因变量，拟合成信号特征奇异值曲线和噪声特征奇异值曲线，并求拟合误差；最后，将拟合误差最小值对应的奇异值阶数确定为有效秩阶次，并进行奇异值分解降噪。通过数值仿真和实际齿轮故障数据分析，表明该方法可以有效地提高信号的信噪比，为后期的故障特征提取创造有利条件。

奇异值分解；降噪；有效秩阶次；拟合误差最小化

旋转机械设备发生故障时，其振动信号很容易受到周围设备的噪声干扰，尤其是早期故障，由于调制源弱，振动信号的故障特征微弱，往往淹没于系统噪声中，难于识别。因此，研究测量信号的预处理技术，最大限度地提高信噪比，是故障特征提取的重要环节。

奇异值分解(Singular Value Decomposition，SVD)降噪是一种非线性滤波方法，可以有效消除信号中的随机噪声，得到相对纯净的振动信号[1]。由于该算法具有良好的稳定性和不变性，在工程上已经得到了广泛应用。

在SVD降噪技术中奇异值有效秩阶次(有效奇异值的个数)的确定问题是当前研究的热点和难点。工程中应用较多的有效秩阶次确定方法是试凑法和寻找奇异值突变点法，这两种方法对操作者的经验要求相对较高，不易掌握。因此，研究人员倾向于寻找操作性强、有量化判据的确定方法。文献[2]提出了奇异值均值法，选取奇异值均值点对应的阶数作为奇异值有效秩阶次，该方法操作简单、贴近工程实用，但估计的阶次往往偏大，导致欠降噪；文献[3]提出了奇异值差分谱法，根据差分谱的最大峰值位置确定有效秩阶次，在文献[3-4]的基础上提出了单边极大值方法，在奇异值差分谱中从右至左，选择第一个至少单边与其相邻峰值差距绝对值最大的极大峰值的对应点位置来确定有效秩阶次，这两种方法在信噪比较高的情况下降噪效果明显，但未考虑有用信号分量的奇异值差分谱也可能出现全谱最大值或单边极大值的情况，可能会得到武断的结论，使估计阶次偏小，导致过降噪；文献[5]采用奇异值相对变化率的最大值确定有效秩阶次，该方法与文献[3]方法类似，相当于对奇异值差分值增加了权值，其效果及存在的问题与文献[3]方法基本一致；文献[6]根据信号快速傅里叶变换结果中主频个数的2倍来确定奇异值有效秩阶次，该方法在仿真信号分析中取得了较好的效果，但是实际旋转机械设备的故障振动信号一般具有调幅-调频特征，频带会发生大幅度迁移，加之噪声的影响，信号频谱往往杂乱，主频个数难以区分，致使方法失效；文献[7]采用结构风险最小化原则确定奇异值有效秩阶次，由于概念抽象，且文中未给出具体实现过程，应用较少。

为此，本文提出了一种新的奇异值有效秩阶次的确定方法，该方法以奇异值及奇异值阶数的拟合误差为判据，选择拟合误差最小值对应的阶数作为有效秩阶次。研究的目的在于探讨一种计算简单、原理清晰、结果准确的振动信号降噪方法，为旋转机械故障诊断的工程应用奠定基础。

1 SVD降噪原理

SVD降噪的原理是利用信号与噪声的能量可分性，对含噪信号构成的矩阵进行分解，仅保留信号特征奇异值，达到去除噪声目的。该方法具体步骤为：

1.1 相空间重构

设原始信号为X=[x1,x2,…,xN]，利用相空间重构理论重构吸引子轨迹矩阵：

(1)

式中：n=N-(m-1)×τ；A为m×n阶的Hankel矩阵；τ为延迟步长；m为嵌入维数[8]。

1.2 奇异值分解

对式(1)所示的Hankel矩阵进行其奇异值分解：

(2)

式中：U是m×m阶正交矩阵,V是n×n阶正交矩阵,Σ是按降序排列的对角矩阵，其对角元素为矩阵A的奇异值σi,ui和vi为列向量,d为非零奇异值的个数,Ai=ui·vi。

1.3 确定奇异值有效秩阶次

根据奇异值分解理论和Frobeious范数意义下矩阵最佳逼近定理可知：前r(r

1.4 信号重构

根据奇异值有效秩阶次，去掉代表噪声信号的较小奇异值(置零)，再进行奇异值分解的逆运算即可实现信号的重构。也就是将上述分组得到的信号分量矩阵转化成长度为N的序列。

SVD降噪的关键问题是Hankel矩阵的构造及奇异值有效秩阶次的确定。前者的理论研究已相对成熟：延迟步长一般取τ=1[9]；在嵌入维数的选取方面，文献[10-11]推荐当N为奇数时，m取中值；当N为偶数时，m=N/2。本文只研究奇异值有效秩阶次的确定问题。

2 基于拟合误差最小原则的有效秩阶次确定

2.1 算法原理

奇异值分解降噪的本质是根据奇异值的大小关系对奇异值分组，即寻找奇异值的突变点。当信号信噪比较高时，奇异值的突变点即为奇异值分布曲线的“肘部”。当信号信噪比较低时，信号特征奇异值与噪声特征奇异值相差不大，整个曲线存在较大的“肘部”，此时，奇异值的突变点已不易寻找。

本文从奇异值分布曲线来寻找信号、噪声分量矩阵的最佳分组。图1给出了低信噪比情况下的原始信号吸引子轨迹矩阵奇异值分布曲线。从曲线上可以看出，奇异值分布曲线可分成三部分：“肘部”、左侧及右侧。左侧的奇异值(反映信号)迅速减小，呈阶梯状分布；右侧的奇异值(反映噪声)缓慢减小，呈直线分布；“肘部”区域的曲线则比较弯曲，该区域属争议区，区内有一点对应有效秩阶次。在“肘部”左侧、右侧分别拟合两条直线，可得到奇异值分布的趋势线，则“肘部”某一点与两条趋势线的接近程度决定了本身的归属。

图1 奇异值分布曲线Fig.1 The distribution of singular values

因此，本文提出的奇异值有效秩阶次确定方法为：按不同的奇异值阶数将奇异值分成信号特征组和噪声特征组；对每次分组的结果，分别采用最小二乘法，以阶数为自变量、以奇异值为因变量，线性拟合成信号特征奇异值曲线和噪声特征奇异值曲线，并求拟合误差；将拟合误差最小所对应的奇异值阶数确定为有效秩阶次。

2.2 实现流程

(1)生成非零奇异值序列Λ及阶数序列n：

Λ=[σ1,σ2,…,σd]

n=[1,2,…,d]

(3)

(2)确定初始阶数r=3；

(3)依据阶数r，将Λ,n分成两组：

Λ1=[σ1,σ2,…,σr]

n1=[1,2,…,r]

(4)

Λ2=[σr+1,σr+2,…,σd]

n2=[r+1,r+2,…,d]

(5)

(4)以阶数为自变量、以奇异值为因变量，分别采用最小二乘法拟合线性函数f1、f2；

(5)计算拟合误差(Fit Square Error，FSE)

(6)

(6)增加阶数r=r+1；

(7)重复步骤(3)～(6)d-5次；

(8)寻找拟合误差FSE的最小值，对应的阶数即为有效秩阶次。

3 仿真数据分析

考察下式所示的仿真信号：

(7)

式中：仿真信号x(t)由调幅-调频信号x1(t)、正弦信号x2(t)及高斯白噪声x3(t)合成,x1(t)+x2(t)为净信号,噪声x3(t)的均值为0、方差为1。采样频率为fs=1 024 Hz，仿真时间t∈[0,1]。仿真信号的时域波形见图2。

按照式(8)计算输入信号信噪比(Signal to Noise Ratio，SNR)：

(8)

式中：s(i)为净信号的第i个数据，n(i)为噪声的第i个数据。得到输入信噪比为-2.618 0 dB。

图2 仿真信号时域波形Fig.2 Simulated signal in time domain waveform

由于文献[4]已分析其方法优于文献[3]方法，文献[6]方法在本文所给例子中均失效，文献[7]方法应用较少，本文只与文献[2,4-5]所给方法进行比较。

文献[2,4-5]及本文方法确定有效秩阶次结果见图3。文献[2]方法确定的阶次为212，文献[4]方法确定的阶次为6，文献[5]方法确定的阶次为2，本文方法确定的阶次为15。

图3 确定仿真信号SVD降噪有效秩阶次Fig.3 Determining effective order rank of SVD denoising

根据有效秩阶次，分别对原始信号进行SVD降噪，降噪结果在图4中给出。可以看出：

(1)文献[2]方法降噪信号(图4(a))与净信号差别较大，信号凌乱、高频噪声明显，降噪不够彻底；

(2)文献[4]方法降噪信号(图4(b))波形形状与净信号类似，但幅值明显偏小，且净信号的调幅特征已被滤掉，出现了过降噪现象；

(3)文献[5]方法降噪信号(图4(c))波形已严重失真，明显降噪过大；

(4)本文方法降噪信号(图4(d))与净信号趋势一致，幅值基本相当，直观上看降噪效果较好。

图4 仿真信号SVD降噪结果时域波形Fig.4 The SVD denoising simulated signal in time domain waveform

为了从数值上比较降噪效果，采用信号均方误差(Mean Square Error，MSE)和信噪比(Signal to Noise Ratio，SNR)来评价4种方法，其定义分别如下：

(9)

(10)

表1是4种方法降噪后的MSE和SNR，可以看出，基于本文方法的SVD降噪效果最好。

表1 四种方法降噪后的信号MSE和SNR

需要说明的是，虽然文献[4-5]方法降噪的MSE和SNR优于文献[2]方法，但这两种方法降噪过大，已将原始信号的有用分量滤掉，造成信息流失。而采用文献[2]这种偏保守的阶次估计方法，虽然降噪不足，信噪比提高不多，但可通过采用其他方法进行二次降噪。因此，有效秩阶次的估计时要尽量避免阶次过小。

4 试验数据分析

齿轮故障试验在QPZZ-Ⅱ旋转机械振动分析及故障诊断试验系统上进行，试验系统的齿轮箱传动结构如图5所示。变频调速电机通过联轴节驱动小齿轮，大齿轮与小齿轮直接啮合，大齿轮通过联轴节带动负载；大、小齿轮均为圆柱齿轮，大齿轮齿数为75，小齿轮齿数为55。

图5 试验系统齿轮箱传动结构图Fig.5 Equipment assembly of gearing box

人为对小齿轮做断齿处理，模拟小齿轮断齿故障。试验中，设置电机轴转速为880 r/min，实测转速866 r/min，则小齿轮的转频f1=14.3 Hz，大齿轮的转频f2=10.6 Hz，齿轮啮合频率为793.8Hz。采用加速度传感器采集振动信号，传感器放置在输出轴电机侧，信号采样频率fs=5 120 Hz，计算采样点N=1 024。

图6给出了原始信号的时域波形图及包络谱。从时域波形图上隐约可以看出原始信号具有调幅-调频特性，周期性冲击信号不显著。在包络谱上可以看出30 Hz处具有谱峰，但不明显，高频干扰较大。

图6 原始信号时域波形及包络谱Fig.6 The original signal in time domain waveform and envelope spectrum

图7 确定试验信号SVD降噪有效秩阶次Fig.7 Determining effective order rank of singular value decomposition denoising

分别采用文献[2,4-5]及本文方法确定SVD降噪有效秩阶次，图7分别给出了相应的结果。文献[2]、文献[4]、文献[5]及本文方法确定的降噪阶次分别为283、5、3及32。

根据有效秩阶次，分别对原始信号进行SVD降噪，结果如图8所示。可以看出：

(1)文献[2]方法降噪信号(图8(a))与原始信号差别不大，原始信号的趋势得以保留，但信号仍然比较凌乱；

(2)文献[4-5]方法降噪信号(图8(b)、图8(c))的脉冲量基本被消除，且幅值明显偏小，很难体现原始信号的时域特征，已将原始信号的有用分量滤掉，明显降噪过大，没有进一步分析的价值；

(3)本文方法降噪信号(图8(d))也保证信号趋势，脉冲信号的周期性有所增强，信号凌乱程度有所减小。

图8 SVD降噪信号时域波形Fig.8 The SVD denoising signal in time domain waveform

为进一步分析降噪效果，对文献[2]及本文方法的降噪信号分别求取Hilbert包络谱，相应结果见图9。可以看出：

(1)文献[2]方法降噪信号包络谱(图9(a))上可看出30 Hz处存在谱峰，但不够明显，且高频噪声依然很大，信噪比较原始信号改善不大；

(2)本文方法降噪信号包络谱(图9(b))上则可以清楚看到小齿轮转频2倍频(30 Hz)处存在全谱最大峰值，且高频噪声明显减小，可直接给出小齿轮故障的正确结论。

图9 SVD降噪信号包络谱Fig.9 The SVD denoising signal envelope spectrum

5 结 论

从仿真信号及试验数据的分析过程可以看出，在奇异值分解降噪有效秩阶次的确定过程中，若采用奇异值差分值(或奇异值相对变化率、奇异熵)的最大值(或单边极大值)作为依据，由于净信号的形式多种多样，难免出现净信号有用分量被误判、进而降噪过大的情况；而采用奇异值均值法确定有效秩阶次又容易出现降噪不足。本文综合考虑信号特征奇异值和噪声特征奇异值的分布曲线，提出了一种新的有效秩阶次确定方法，以奇异值及奇异值阶数的曲线拟合误差为判据，选择拟合误差最小值对应的阶数作为有效秩阶次。该方法便于实施，物理意义清晰，确定的阶次准确，降噪效果好，无过、欠降噪现象，能在保留信号特征的前提下，有效地消除噪声。其应用不止于旋转机械振动信号的降噪，具有一定的通用性。

[ 1 ] GOLYANDINA N, ZHIGLJAVSKY A. Singular spectrum analysis for time series [M]. Springer, 2013.

[ 2 ] 王益艳. 基于特征均值的SVD信号去噪算法[J]. 计算机应用与软件, 2012, 29(5): 121-123. WANG Yiyan. Mean value of evenval-based SVD signal denoising algorithm[J]. Computer Applications and Software，2012，29(5): 121-123.

[ 3 ] 赵学智, 叶邦彦, 陈统坚. 奇异值差分谱理论及其在车床主轴箱故障诊断中的应用[J]. 机械工程学报, 2010, 46(1): 100-108.

ZHAO Xuezhi, YE Bangyan, CHEN Tongjian. Difference spectrum theory of singular value and its application to the fault diagnosis of headstock of lathe[J]. Journal of Mechanical Engineering, 2010, 46(1):100-108.

[ 4 ] 王建国, 李健, 刘颖源. 一种确定奇异值分解降噪有效秩阶次的改进方法[J]. 振动与冲击, 2014, 33(12): 176-180. WANG Jianguo, LI Jian, LIU Yingyuan. An improved method for determining effective order rank of SVD denoising [J]. Journal of Vibration and Shock, 2014, 33(12): 176-180.

[ 5 ] 王树青, 林裕裕, 孟元栋,等. 一种基于奇异值分解技术的模型定阶方法[J]. 振动与冲击, 2012, 31(15): 87-91. WANG Shuqing, LIN Yuyu, MENG Yuandong, et al. Model order determination based on singular valued decomposition[J]. Journal of Vibration and Shock, 2012, 31(15):87-91.

[ 6 ] 钱征文, 程礼, 李应红. 利用奇异值分解的信号降噪方法[J].振动、测试与诊断, 2011, 31(4): 459-463. QIAN Zhengwen, CHENG Li, LI Yinghong. Signal denoising method by means of SVD[J]. Journal of Vibration Measurement & Diagnosis, 2011, 31(4): 459-463.

[ 7 ] 朱启兵, 刘杰, 李允公. 基于结构风险最小化原则的奇异值分解降噪研究[J]. 振动工程学报, 2005, 18(2): 204-207. ZHU Qibing, LIU Jie, LI Yungong. Study on noise reduction in singular value decomposition based on structural risk minimization[J]. Journal of Vibration Engineering, 2005, 18(2): 204-207.

[ 8 ] 赵学智, 陈统坚, 叶邦彦. 基于奇异值分解的铣削力信号处理与铣床状态信息分离[J]. 机械工程学报, 2007, 43(6): 169-174. ZHAO Xuezhi, CHEN Tongjian, YE Bangyan. Processingof milling force signal and isolation of state information of milling machine based on singular value decomposition[J]. Chinese Journal of Mechanical Engineering, 2007, 43(6): 169-174.

[ 9 ] 孟智慧, 王晶. 基于SVD降噪和谱峭度的滚动轴承故障诊断[J]. 轴承, 2013(10): 52-55. MENG Zhihui, WANG Jing. Faultdiagnosis for rolling bearings based on SVD denoising and spectral kurtosis[J]. Bearing, 2013(10): 52-55.

[10] GOLYANDINA N. On the choice of parameters in singular spectrum analysis and related subspace-based methods [DB/OL]. http://arxiv.org/abs/1005.4374, 2010-07-20.

[11] MAHMOUDVAND R, ZOKAEI M. On the singular values of the Hankel matrix with application in singular spectrum analysis[J]. Chilean Journal of Statistics, 2012, 3(1): 43-56.

A new method for determining effective rank order of singularvalue decomposition denoising based on fitting error minimum principle

CUI Weicheng1, XU Aiqiang2, LI Wei1, MENG Fanlei1

(1. Department of Aircraft Engineering, Naval Aeronautical and Astronautical University ,Yantai 264001, China;2. Institute of Aircraft Detection and Application, Naval Aeronautical and Astronautical University, Yantai 264001, China)

In order to maximize the signal-to-noise ratio of a rotating mechanical equipment’s fault vibration signals, the singular value decomposition (SVD) de-noising method was studied, and a new method for determining its effective rank order was proposed. Firstly, a vibration signal was reconstructed in phase space, and the singular value decomposition of the attractor trajectory matrix was performed. Secondly, singular values were divided into a signal group and a noise group. For the results of each grouping, rank and singular value were taken as independent variable and dependent one, respectively. The feature singular value curve of signal and the feature singular value curve of noise were fitted, then the fitting errors were solved. At last, the singular value order corresponding to the minimum fitting error was taken as the effective rank order, and the SVD de-noising was performed. The results of numerical simulation and actual gear fault data analysis showed that the proposed method can effectively improve signal-to-noise ratios of signals, and create a beneficial condition for the subsequent fault feature extraction.

singular value decomposition (SVD); noise reduction; fitting error minimum principle

国家部委预研基金资助(9140A27020214JB1446)

2015-09-28 修改稿收到日期：2016-01-05

崔伟成 男，博士生，讲师，1981年6月生

TN911.7；TP206.3

A

10.13465/j.cnki.jvs.2017.03.021