基于擬合誤差最小化原則的奇異值分解降噪有效秩階次確定方法

崔偉成， 許愛強， 李 偉， 孟凡磊

(1. 海軍航空工程學院 飛行器工程系，山東 煙臺 264001； 2. 海軍航空工程學院 飛行器檢測與應用研究所，山東 煙臺 264001)

基于擬合誤差最小化原則的奇異值分解降噪有效秩階次確定方法

崔偉成1， 許愛強2， 李 偉1， 孟凡磊1

(1. 海軍航空工程學院 飛行器工程系，山東 煙臺 264001； 2. 海軍航空工程學院 飛行器檢測與應用研究所，山東 煙臺 264001)

為了最大限度地提高旋轉機械設備故障振動信號的信噪比，研究了奇異值分解降噪的原理，提出了一種新的奇異值分解降噪有效秩階次的確定方法。首先，對振動信號進行相空間重構，對吸引子軌跡矩陣進行奇異值分解；然后，按不同的階數(shù)，將奇異值分成信號組和噪聲組，對每次分組的結果，以階數(shù)為自變量、以奇異值為因變量，擬合成信號特征奇異值曲線和噪聲特征奇異值曲線，并求擬合誤差；最后，將擬合誤差最小值對應的奇異值階數(shù)確定為有效秩階次，并進行奇異值分解降噪。通過數(shù)值仿真和實際齒輪故障數(shù)據(jù)分析，表明該方法可以有效地提高信號的信噪比，為后期的故障特征提取創(chuàng)造有利條件。

奇異值分解；降噪；有效秩階次；擬合誤差最小化

旋轉機械設備發(fā)生故障時，其振動信號很容易受到周圍設備的噪聲干擾，尤其是早期故障，由于調(diào)制源弱，振動信號的故障特征微弱，往往淹沒于系統(tǒng)噪聲中，難于識別。因此，研究測量信號的預處理技術，最大限度地提高信噪比，是故障特征提取的重要環(huán)節(jié)。

奇異值分解(Singular Value Decomposition，SVD)降噪是一種非線性濾波方法，可以有效消除信號中的隨機噪聲，得到相對純凈的振動信號[1]。由于該算法具有良好的穩(wěn)定性和不變性，在工程上已經(jīng)得到了廣泛應用。

在SVD降噪技術中奇異值有效秩階次(有效奇異值的個數(shù))的確定問題是當前研究的熱點和難點。工程中應用較多的有效秩階次確定方法是試湊法和尋找奇異值突變點法，這兩種方法對操作者的經(jīng)驗要求相對較高，不易掌握。因此，研究人員傾向于尋找操作性強、有量化判據(jù)的確定方法。文獻[2]提出了奇異值均值法，選取奇異值均值點對應的階數(shù)作為奇異值有效秩階次，該方法操作簡單、貼近工程實用，但估計的階次往往偏大，導致欠降噪；文獻[3]提出了奇異值差分譜法，根據(jù)差分譜的最大峰值位置確定有效秩階次，在文獻[3-4]的基礎上提出了單邊極大值方法，在奇異值差分譜中從右至左，選擇第一個至少單邊與其相鄰峰值差距絕對值最大的極大峰值的對應點位置來確定有效秩階次，這兩種方法在信噪比較高的情況下降噪效果明顯，但未考慮有用信號分量的奇異值差分譜也可能出現(xiàn)全譜最大值或單邊極大值的情況，可能會得到武斷的結論，使估計階次偏小，導致過降噪；文獻[5]采用奇異值相對變化率的最大值確定有效秩階次，該方法與文獻[3]方法類似，相當于對奇異值差分值增加了權值，其效果及存在的問題與文獻[3]方法基本一致；文獻[6]根據(jù)信號快速傅里葉變換結果中主頻個數(shù)的2倍來確定奇異值有效秩階次，該方法在仿真信號分析中取得了較好的效果，但是實際旋轉機械設備的故障振動信號一般具有調(diào)幅-調(diào)頻特征，頻帶會發(fā)生大幅度遷移，加之噪聲的影響，信號頻譜往往雜亂，主頻個數(shù)難以區(qū)分，致使方法失效；文獻[7]采用結構風險最小化原則確定奇異值有效秩階次，由于概念抽象，且文中未給出具體實現(xiàn)過程，應用較少。

為此，本文提出了一種新的奇異值有效秩階次的確定方法，該方法以奇異值及奇異值階數(shù)的擬合誤差為判據(jù)，選擇擬合誤差最小值對應的階數(shù)作為有效秩階次。研究的目的在于探討一種計算簡單、原理清晰、結果準確的振動信號降噪方法，為旋轉機械故障診斷的工程應用奠定基礎。

1 SVD降噪原理

SVD降噪的原理是利用信號與噪聲的能量可分性，對含噪信號構成的矩陣進行分解，僅保留信號特征奇異值，達到去除噪聲目的。該方法具體步驟為：

1.1 相空間重構

設原始信號為X=[x1,x2,…,xN]，利用相空間重構理論重構吸引子軌跡矩陣：

(1)

式中：n=N-(m-1)×τ；A為m×n階的Hankel矩陣；τ為延遲步長；m為嵌入維數(shù)[8]。

1.2 奇異值分解

對式(1)所示的Hankel矩陣進行其奇異值分解：

(2)

式中：U是m×m階正交矩陣,V是n×n階正交矩陣,Σ是按降序排列的對角矩陣，其對角元素為矩陣A的奇異值σi,ui和vi為列向量,d為非零奇異值的個數(shù),Ai=ui·vi。

1.3 確定奇異值有效秩階次

根據(jù)奇異值分解理論和Frobeious范數(shù)意義下矩陣最佳逼近定理可知：前r(r

1.4 信號重構

根據(jù)奇異值有效秩階次，去掉代表噪聲信號的較小奇異值(置零)，再進行奇異值分解的逆運算即可實現(xiàn)信號的重構。也就是將上述分組得到的信號分量矩陣轉化成長度為N的序列。

SVD降噪的關鍵問題是Hankel矩陣的構造及奇異值有效秩階次的確定。前者的理論研究已相對成熟：延遲步長一般取τ=1[9]；在嵌入維數(shù)的選取方面，文獻[10-11]推薦當N為奇數(shù)時，m取中值；當N為偶數(shù)時，m=N/2。本文只研究奇異值有效秩階次的確定問題。

2 基于擬合誤差最小原則的有效秩階次確定

2.1 算法原理

奇異值分解降噪的本質(zhì)是根據(jù)奇異值的大小關系對奇異值分組，即尋找奇異值的突變點。當信號信噪比較高時，奇異值的突變點即為奇異值分布曲線的“肘部”。當信號信噪比較低時，信號特征奇異值與噪聲特征奇異值相差不大，整個曲線存在較大的“肘部”，此時，奇異值的突變點已不易尋找。

本文從奇異值分布曲線來尋找信號、噪聲分量矩陣的最佳分組。圖1給出了低信噪比情況下的原始信號吸引子軌跡矩陣奇異值分布曲線。從曲線上可以看出，奇異值分布曲線可分成三部分：“肘部”、左側及右側。左側的奇異值(反映信號)迅速減小，呈階梯狀分布；右側的奇異值(反映噪聲)緩慢減小，呈直線分布；“肘部”區(qū)域的曲線則比較彎曲，該區(qū)域屬爭議區(qū)，區(qū)內(nèi)有一點對應有效秩階次。在“肘部”左側、右側分別擬合兩條直線，可得到奇異值分布的趨勢線，則“肘部”某一點與兩條趨勢線的接近程度決定了本身的歸屬。

圖1 奇異值分布曲線Fig.1 The distribution of singular values

因此，本文提出的奇異值有效秩階次確定方法為：按不同的奇異值階數(shù)將奇異值分成信號特征組和噪聲特征組；對每次分組的結果，分別采用最小二乘法，以階數(shù)為自變量、以奇異值為因變量，線性擬合成信號特征奇異值曲線和噪聲特征奇異值曲線，并求擬合誤差；將擬合誤差最小所對應的奇異值階數(shù)確定為有效秩階次。

2.2 實現(xiàn)流程

(1)生成非零奇異值序列Λ及階數(shù)序列n：

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

n=[1,2,…,d]

(3)

(2)確定初始階數(shù)r=3；

(3)依據(jù)階數(shù)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)以階數(shù)為自變量、以奇異值為因變量，分別采用最小二乘法擬合線性函數(shù)f1、f2；

(5)計算擬合誤差(Fit Square Error，F(xiàn)SE)

(6)

(6)增加階數(shù)r=r+1；

(7)重復步驟(3)～(6)d-5次；

(8)尋找擬合誤差FSE的最小值，對應的階數(shù)即為有效秩階次。

3 仿真數(shù)據(jù)分析

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

(7)

式中：仿真信號x(t)由調(diào)幅-調(diào)頻信號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個數(shù)據(jù)，n(i)為噪聲的第i個數(shù)據(jù)。得到輸入信噪比為-2.618 0 dB。

圖2 仿真信號時域波形Fig.2 Simulated signal in time domain waveform

由于文獻[4]已分析其方法優(yōu)于文獻[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

根據(jù)有效秩階次，分別對原始信號進行SVD降噪，降噪結果在圖4中給出。可以看出：

(1)文獻[2]方法降噪信號(圖4(a))與凈信號差別較大，信號凌亂、高頻噪聲明顯，降噪不夠徹底；

(2)文獻[4]方法降噪信號(圖4(b))波形形狀與凈信號類似，但幅值明顯偏小，且凈信號的調(diào)幅特征已被濾掉，出現(xiàn)了過降噪現(xiàn)象；

(3)文獻[5]方法降噪信號(圖4(c))波形已嚴重失真，明顯降噪過大；

(4)本文方法降噪信號(圖4(d))與凈信號趨勢一致，幅值基本相當，直觀上看降噪效果較好。

圖4 仿真信號SVD降噪結果時域波形Fig.4 The SVD denoising simulated signal in time domain waveform

為了從數(shù)值上比較降噪效果，采用信號均方誤差(Mean Square Error，MSE)和信噪比(Signal to Noise Ratio，SNR)來評價4種方法，其定義分別如下：

(9)

(10)

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

表1 四種方法降噪后的信號MSE和SNR

需要說明的是，雖然文獻[4-5]方法降噪的MSE和SNR優(yōu)于文獻[2]方法，但這兩種方法降噪過大，已將原始信號的有用分量濾掉，造成信息流失。而采用文獻[2]這種偏保守的階次估計方法，雖然降噪不足，信噪比提高不多，但可通過采用其他方法進行二次降噪。因此，有效秩階次的估計時要盡量避免階次過小。

4 試驗數(shù)據(jù)分析

齒輪故障試驗在QPZZ-Ⅱ旋轉機械振動分析及故障診斷試驗系統(tǒng)上進行，試驗系統(tǒng)的齒輪箱傳動結構如圖5所示。變頻調(diào)速電機通過聯(lián)軸節(jié)驅動小齒輪，大齒輪與小齒輪直接嚙合，大齒輪通過聯(lián)軸節(jié)帶動負載；大、小齒輪均為圓柱齒輪，大齒輪齒數(shù)為75，小齒輪齒數(shù)為55。

圖5 試驗系統(tǒng)齒輪箱傳動結構圖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給出了原始信號的時域波形圖及包絡譜。從時域波形圖上隱約可以看出原始信號具有調(diào)幅-調(diào)頻特性，周期性沖擊信號不顯著。在包絡譜上可以看出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。

根據(jù)有效秩階次，分別對原始信號進行SVD降噪，結果如圖8所示?？梢钥闯觯?/p>

(1)文獻[2]方法降噪信號(圖8(a))與原始信號差別不大，原始信號的趨勢得以保留，但信號仍然比較凌亂；

(2)文獻[4-5]方法降噪信號(圖8(b)、圖8(c))的脈沖量基本被消除，且幅值明顯偏小，很難體現(xiàn)原始信號的時域特征，已將原始信號的有用分量濾掉，明顯降噪過大，沒有進一步分析的價值；

(3)本文方法降噪信號(圖8(d))也保證信號趨勢，脈沖信號的周期性有所增強，信號凌亂程度有所減小。

圖8 SVD降噪信號時域波形Fig.8 The SVD denoising signal in time domain waveform

為進一步分析降噪效果，對文獻[2]及本文方法的降噪信號分別求取Hilbert包絡譜，相應結果見圖9?？梢钥闯觯?/p>

(1)文獻[2]方法降噪信號包絡譜(圖9(a))上可看出30 Hz處存在譜峰，但不夠明顯，且高頻噪聲依然很大，信噪比較原始信號改善不大；

(2)本文方法降噪信號包絡譜(圖9(b))上則可以清楚看到小齒輪轉頻2倍頻(30 Hz)處存在全譜最大峰值，且高頻噪聲明顯減小，可直接給出小齒輪故障的正確結論。

圖9 SVD降噪信號包絡譜Fig.9 The SVD denoising signal envelope spectrum

5 結 論

從仿真信號及試驗數(shù)據(jù)的分析過程可以看出，在奇異值分解降噪有效秩階次的確定過程中，若采用奇異值差分值(或奇異值相對變化率、奇異熵)的最大值(或單邊極大值)作為依據(jù)，由于凈信號的形式多種多樣，難免出現(xiàn)凈信號有用分量被誤判、進而降噪過大的情況；而采用奇異值均值法確定有效秩階次又容易出現(xiàn)降噪不足。本文綜合考慮信號特征奇異值和噪聲特征奇異值的分布曲線，提出了一種新的有效秩階次確定方法，以奇異值及奇異值階數(shù)的曲線擬合誤差為判據(jù)，選擇擬合誤差最小值對應的階數(shù)作為有效秩階次。該方法便于實施，物理意義清晰，確定的階次準確，降噪效果好，無過、欠降噪現(xiàn)象，能在保留信號特征的前提下，有效地消除噪聲。其應用不止于旋轉機械振動信號的降噪，具有一定的通用性。

[ 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 ] 趙學智, 葉邦彥, 陳統(tǒng)堅. 奇異值差分譜理論及其在車床主軸箱故障診斷中的應用[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 ] 趙學智, 陳統(tǒng)堅, 葉邦彥. 基于奇異值分解的銑削力信號處理與銑床狀態(tài)信息分離[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