基于改進(jìn)奇異譜分析的信號(hào)去噪方法

戴豪民，許愛強(qiáng)，孫偉超

(海軍航空工程學(xué)院 飛行器檢測(cè)與應(yīng)用研究所，山東，煙臺(tái) 264001)

?

基于改進(jìn)奇異譜分析的信號(hào)去噪方法

戴豪民，許愛強(qiáng)，孫偉超

(海軍航空工程學(xué)院 飛行器檢測(cè)與應(yīng)用研究所，山東，煙臺(tái) 264001)

傳統(tǒng)的去噪方法，比如小波閾值去噪，它只對(duì)高斯噪聲有效，對(duì)于脈沖噪聲卻無能為力. 近年來發(fā)展起來的奇異譜分析方法可以在高信噪比的條件下很好地濾除上述兩類噪聲，但該方法降噪過程涉及了一定的主觀因素，并且受矩陣擾動(dòng)理論的限制，該方法隨著信噪比的降低，去噪能力也隨之下降. 針對(duì)上述情況，提出一種改進(jìn)算法，將矩陣秩最小化理論應(yīng)用于奇異譜分析方法中. 仿真結(jié)果表明，改進(jìn)算法去噪效果明顯，能夠最大限度降低信號(hào)均方誤差，提高信噪比，增強(qiáng)奇異譜分析方法的通用性.

奇異譜分析；奇異值分解；矩陣擾動(dòng)理論；秩最小化理論

信號(hào)的產(chǎn)生、采集和傳遞都不可避免地會(huì)受到噪聲的干擾，因此在信號(hào)預(yù)處理時(shí)進(jìn)行去噪是十分必要的. 傳統(tǒng)的去噪方法，比如基于譜分析的濾波器和小波去噪，先將原始信號(hào)從時(shí)域變換到頻域，通過去除信號(hào)的某些頻段成分的方法實(shí)現(xiàn)降噪. 通常，噪聲表現(xiàn)為高頻小幅值信號(hào)，信號(hào)通過低通濾波或小波閾值處理即可去除高頻成分. 然而，有些情況下噪聲的頻帶會(huì)分布在整個(gè)頻率軸上，這樣傳統(tǒng)的方法就很難將有用信號(hào)和噪聲的頻譜區(qū)分開. 因此，傳統(tǒng)方法具有一定的局限性. 近年來，在奇異值分解的基礎(chǔ)上發(fā)展起來的奇異譜分析(singular spectrum analysis，SSA)技術(shù)已經(jīng)應(yīng)用于信號(hào)的去噪處理中[1]，由于該技術(shù)獨(dú)立于信號(hào)模型，并且具有不受噪聲頻譜分布影響的特點(diǎn)，所以能夠達(dá)到自適應(yīng)降噪的效果[2]. 但SSA亦有其不足之處，降噪過程涉及了一定的主觀因素. 當(dāng)原始信號(hào)受到輕微隨機(jī)高斯噪聲污染時(shí)，可以通過人為設(shè)定閾值，提取原信號(hào)軌跡矩陣的較大奇異值集合構(gòu)成信號(hào)子空間，較小奇異值集合構(gòu)成噪聲子空間，實(shí)現(xiàn)信噪分離. 當(dāng)原始信號(hào)被噪聲嚴(yán)重污染，即使是極少數(shù)的大噪聲，也會(huì)使該方法分析失效，尤其在復(fù)雜背景噪聲環(huán)境中，噪聲可能既包括輕微高斯噪聲，又包括非平穩(wěn)沖擊噪聲. 為此，針對(duì)SSA固有缺陷問題，本文提出一種改進(jìn)SSA算法，將矩陣秩最小化理論應(yīng)用于SSA中，實(shí)現(xiàn)信噪分離，增強(qiáng)SSA的通用性.

1 奇異譜分析

SSA是一種全局分析法，利用相空間重構(gòu)的基本思想，通過奇異值分解以達(dá)到識(shí)別原始信號(hào)成份(如趨勢(shì)、周期或準(zhǔn)周期、噪聲等)的目的[3]. 該方法包括分解和重構(gòu)兩個(gè)階段，分解階段包括嵌入操作和奇異值分解，重構(gòu)階段包括分組和對(duì)角平均化，其具體步驟如下：

① 嵌入. 嵌入操作可以視為一維序列向多維序列的映射. 設(shè)長度為N的一維實(shí)序列FN=(f1，f2，…，fN)，正整數(shù)L為滑動(dòng)窗口長度，1

Xj=[fjfj+1… fj+L-1]T∈RL，

其中K=N-L+1，j=1，2，…，K. 映射的結(jié)果形成軌跡矩陣

X=[X1X2…XK]∈RL×K，

該軌跡矩陣是Hankel矩陣，可以表示為

(1)

② 奇異值分解. 對(duì)軌跡矩陣X進(jìn)行奇異值分解，得到

(2)

其中d為X的非零奇異值個(gè)數(shù)，顯然d=rank(X)≤min(L，M)，λ1，λ2，…，λd為按降序排列的X的奇異值，Ui和Vi分別是X的左右奇異向量.

③ 分組. 分組的目的是分離信號(hào)中的加性成分. 如果對(duì)原始信號(hào)進(jìn)行降噪處理，那么分組操作就是將原始序列FN構(gòu)造的軌跡矩陣X表示成有用信號(hào)S和噪聲E之和，即X=S+E. 在信號(hào)處理領(lǐng)域，通常認(rèn)為前r(r

(3)

2 SSA的參數(shù)選擇

2.1 窗口長度L的選擇

上述SSA過程中，涉及到兩個(gè)待定參數(shù)，即窗口長度L和重構(gòu)信號(hào)的奇異值數(shù)目r. 由于窗口長度分別為L和K=N-L+1的軌跡矩陣奇異值分解是對(duì)稱的，所以窗口長度L超過N/2是沒有必要的. 同時(shí)，窗口長度L越大，序列的分解就會(huì)越精細(xì)，過小的L往往會(huì)造成序列中可解釋的成分(如趨勢(shì)、周期或準(zhǔn)周期、噪聲等)混雜在一起，不利于這些成分的提取[4]. 如果軌跡矩陣的秩遠(yuǎn)小于L，那么過大的L也是無意義，反而會(huì)增加奇異值分解的運(yùn)算時(shí)間. 現(xiàn)實(shí)中不同信號(hào)軌跡矩陣的秩往往不同，但有用信號(hào)矩陣S往往是低秩的，即rank(S)

對(duì)于完全由噪聲序列構(gòu)成的矩陣E，往往是滿秩矩陣，并且所有的奇異值相差較小，奇異值衰減較慢. 所以選擇大的L往往可以挖掘出隱藏在原始信號(hào)中的噪聲成分. 文獻(xiàn)[5]從重構(gòu)序列各成分的分離性和重構(gòu)誤差的角度上推薦L取序列長度N的中值(如果N為偶數(shù)，取L=N/2). 文獻(xiàn)[6]研究了軌跡矩陣的奇異值隨窗口長度L的變化規(guī)律，從理論上證明了當(dāng)L取N的中值時(shí)，對(duì)應(yīng)的奇異值會(huì)取得最大值，同時(shí)該文獻(xiàn)也推薦在大多數(shù)情況下L取N的中值是最佳的選擇.

2.2 奇異值數(shù)目r的選擇

分組過程中，如何客觀選取奇異值數(shù)目r進(jìn)行重構(gòu)信號(hào)，目前還沒有通用方法. 相關(guān)文獻(xiàn)的選擇方法往往具有一定的主觀性. 文獻(xiàn)[7]提出奇異值均值法，即將所有低于奇異值平均值的奇異值置0，該方法只對(duì)原始信號(hào)受到輕微隨機(jī)高斯噪聲污染時(shí)有效，此時(shí)原始信號(hào)的軌跡矩陣前10%甚至1%的奇異值的和就占了全部奇異值之和的99%以上了，去掉一些較小的奇異值，就可以最大程度逼近有用信號(hào)矩陣，這也與主成分分析的思想一致. 文獻(xiàn)[8]利用信號(hào)快速傅里葉變換(fastFouriertransform，F(xiàn)FT)結(jié)果中主頻率個(gè)數(shù)確定奇異值數(shù)目，該方法只對(duì)平穩(wěn)信號(hào)有效，對(duì)于非平穩(wěn)信號(hào)，如果信號(hào)經(jīng)過快速傅里葉變換后主頻特征不明顯或者被噪聲頻率淹沒，那么就無法正確判定主頻率個(gè)數(shù). 尤其在復(fù)雜噪聲背景下，上述這些方法就會(huì)失效，接下來本文從矩陣擾動(dòng)理論角度闡述失效的原因.

3 矩陣擾動(dòng)理論

設(shè)含噪序列FN=(f1，f2，…，fN)，其加法模型可以表示為FN=SN+EN，SN為有用信號(hào)序列，EN為噪聲序列，N為序列長度. 設(shè)正整數(shù)L為滑動(dòng)窗口長度，1

(4)

假設(shè)噪聲與有用信號(hào)統(tǒng)計(jì)不相關(guān)，則

XXT=SST+EET.

(5)

(6)

其中I為單位矩陣. 令rank(SST)=r≤L，如果SST的特征值分解為SST=UΛUT，U為特征向量，那么矩陣XXT的特征值分解可以表示為

XXT=UΛUT+σ2eI=U(Λ+σ2eI)UT=UΣUT.

(7)

(8)

(9)

所以選擇前r個(gè)奇異值逼近有用信號(hào)矩陣是合理的. 但是，在低信噪比情況下或者信號(hào)受到稀疏大噪聲污染，雖然可以通過截?cái)嗥娈愔狄瞥糠衷肼?，但仍有部分噪聲?duì)含噪矩陣前r個(gè)奇異值產(chǎn)生擾動(dòng)，重構(gòu)信號(hào)也將產(chǎn)生較大誤差.

4 基于秩最小化理論的去噪方法

相關(guān)文獻(xiàn)[10-11]指出高維的冗余數(shù)據(jù)本質(zhì)上位于低維的子空間，這樣的高維數(shù)據(jù)可以認(rèn)為是稀疏的，其稀疏性表現(xiàn)在低秩性上，即可以通過計(jì)算數(shù)據(jù)矩陣X的最佳低秩矩陣S獲得其最佳逼近. 在工程領(lǐng)域，潛在的低秩矩陣S可能被輕微噪聲E所污染，也可能被稀疏大噪聲W所污染. 當(dāng)矩陣元素被嚴(yán)重破壞，傳統(tǒng)方法已經(jīng)不能有效從污染的數(shù)據(jù)中恢復(fù)出“干凈”的數(shù)據(jù). 受到壓縮感知理論的啟發(fā)，John Wright提出了魯棒主成分分析(robust principal component analysis，RPCA)算法，并且證明了通過解優(yōu)化問題就能夠從被嚴(yán)重污染矩陣中精確恢復(fù)出所需的低秩矩陣[12]. 當(dāng)數(shù)據(jù)矩陣有非常良好的結(jié)構(gòu)，即低秩的，并且只有很少一部分元素被嚴(yán)重破壞，即噪聲是稀疏的但大小可以任意，矩陣去噪(或矩陣恢復(fù))可用如下優(yōu)化問題來描述：

(10)

(11)

這樣就將求解矩陣秩的非凸問題轉(zhuǎn)換到凸優(yōu)化問題，使得這個(gè)問題具有唯一解.

(12)

目前，求解上述最優(yōu)化問題有多種算法，如迭代閾值法[16](iterativethresholding，IT)、奇異值閾值法[17](singularvaluethresholding，SVT)、加速近端梯度法[18](acceleratedproximalgradient，APG)、精確增廣拉格朗日乘子法[19](exactaugmentedLagrangemultiplier，EALM)以及非精確增廣拉格朗日乘子法(inexactaugmentedLagrangemultiplier，IALM)等. 從各算法的時(shí)間開銷和誤差精度角度上比較，IALM法最佳[20]. 在仿真實(shí)驗(yàn)中，本文就是利用IALM法解決式(12)和(13)中的優(yōu)化問題，算法的詳細(xì)描述及收斂性證明請(qǐng)參見文獻(xiàn)[19].

5 仿真實(shí)驗(yàn)

綜合上述分析，為去除信號(hào)中的噪聲干擾，將秩最小化理論應(yīng)用于SSA中實(shí)現(xiàn)去噪的步驟如下：

① 通過嵌入操作，將長度N的含噪序列FN=(f1，f2，…，fN)映射為軌跡矩陣

X=[X1X2…XK]∈RL×K，

K=N-L+1，L=N/2.

② 根據(jù)模型(13)，利用IALM算法恢復(fù)出低秩矩陣S.

③ 將低秩矩陣S對(duì)角平均化，轉(zhuǎn)化成所需的長度為N的有用信號(hào)序列，最終實(shí)現(xiàn)信噪分離.

為了驗(yàn)證改進(jìn)SSA去噪效果，本文分別采用原始SSA和小波閾值去噪方法與之進(jìn)行對(duì)比. 此外，為了考察本文方法的適用性，以文獻(xiàn)[20]的調(diào)幅調(diào)頻信號(hào)為例，分別通過3個(gè)實(shí)驗(yàn)考慮3種背景噪聲環(huán)境下的去噪效果，同時(shí)，采用信號(hào)均方誤差(meansquareerror，MSE)和信噪比(signaltonoiseratio，SNR)來衡量各種方法的優(yōu)劣，其定義形式如下：

(13)

(14)

實(shí)驗(yàn)1 考察以下仿真信號(hào)：

0.2sin10πt)+sin80πt+e(t).

式中，e(t)分別為信噪比是-5，3和10 dB的高斯白噪聲，原始SSA算法采用文獻(xiàn)[7]提出的奇異值均值法確定奇異值有效階次. 信號(hào)采樣頻率為1 024 Hz，采樣時(shí)間1 s，表1和表2是各種方法去噪后的MSE和SNR. 從這兩個(gè)表中可以看出，在低信噪比情況下，小波閾值去噪和文獻(xiàn)[7]方法都不能有效的濾除高斯白噪聲；在中等信噪比情況下，這兩種方法的去噪效果均有所改善，小波閾值去噪方法更好；在高信噪比情況下，3種方法都可以很好地去除掉高斯白噪聲. 但上述3種噪聲背景環(huán)境，本文方法表現(xiàn)俱佳.

表1 3種方法去噪后的信號(hào)MSE

實(shí)驗(yàn)2 考察以下仿真信號(hào)：

0.2sin10πt)+sin80πt+w(t).

表2 3種方法去噪后的信號(hào)SNR

其中，w(t)為隨機(jī)添加的17個(gè)脈沖噪聲. 原始SSA算法采用文獻(xiàn)[8]提出的方法，該文獻(xiàn)指出重構(gòu)信號(hào)所需的奇異值個(gè)數(shù)是仿真信號(hào)功率譜主頻數(shù)目的2倍，從圖1可以看出，仿真信號(hào)功率譜的主頻數(shù)目是2，所以本文采用前4個(gè)奇異值恢復(fù)原始干凈信號(hào).

表3是各種方法去噪后的信號(hào)MSE和SNR. 在稀疏大噪聲背景環(huán)境下，本文方法表現(xiàn)特別優(yōu)秀，幾乎完全恢復(fù)干凈信號(hào)，其誤差精度可以忽略不計(jì). 文獻(xiàn)[8]方法也能濾除稀疏大噪聲的干擾，其MSE和SNR與本文方法比均有明顯劣勢(shì). 小波閾值去噪方法則不能去除掉這種脈沖噪聲. 從圖2也可以看出，文獻(xiàn)[8]方法恢復(fù)出的信號(hào)已經(jīng)出現(xiàn)失真，而小波閾值去噪方法中脈沖干擾依然存在.

表3 3種方法去噪后的信號(hào)MES和SNR

實(shí)驗(yàn)3 考察以下仿真信號(hào)：

(1+sin5πt)cos(20πt+0.2sin10πt)+

式中：e(t)為信噪比是-5 dB的高斯白噪聲；w(t)為隨機(jī)添加的17個(gè)脈沖噪聲. 文獻(xiàn)[20]提出對(duì)信號(hào)進(jìn)行中值濾波可以有效抑制脈沖干擾，原始SSA算法采用文獻(xiàn)[20]和文獻(xiàn)[8]相結(jié)合的方法，首先對(duì)信號(hào)進(jìn)行中值濾波，然后采用FFT確定信號(hào)主頻個(gè)數(shù)以及奇異值數(shù)目. 從圖3可以看出，中值濾波后信號(hào)功率譜的主頻數(shù)目是2，所以本文采用前4個(gè)奇異值恢復(fù)干凈信號(hào).

表4是各種方法去噪后的MES和SNR. 在這種復(fù)雜噪聲背景環(huán)境中，中值濾波和FFT相結(jié)合的方法去噪效果較好，中值濾波起到了抑制脈沖干擾的作用，但中值濾波也會(huì)使信號(hào)變得更為光滑，從圖4可以看出，該方法恢復(fù)出的信號(hào)與原始干凈信號(hào)相比失真明顯. 小波閾值去噪方法只對(duì)高信噪比的高斯白噪聲有效，對(duì)于脈沖干擾幾乎沒有效果.

表4 4種方法去噪后的信號(hào)MES和SNR

6 結(jié) 論

復(fù)雜的噪聲背景環(huán)境中，信號(hào)往往既被高斯噪聲污染，又被脈沖噪聲污染，傳統(tǒng)的去噪方法比如小波閾值去噪只對(duì)高斯噪聲有效，對(duì)于脈沖干擾卻無能為力. 原始SSA方法結(jié)合中值濾波能夠?qū)γ}沖噪聲起到一定的抑制作用，但同時(shí)也會(huì)使原始信號(hào)嚴(yán)重失真. 本文將矩陣秩最小化理論應(yīng)用于SSA中，從仿真結(jié)果上看，去噪效果明顯，能夠最大限度的降低信號(hào)均方誤差，提高信噪比尤其對(duì)脈沖噪聲，幾乎可以完全實(shí)現(xiàn)信噪分離. 此外，由于仿真實(shí)驗(yàn)僅僅以高斯白噪聲為參考，未涉及其他噪聲類型的去噪效果，所以這也是本文下一步研究方向.

[1] Hassani H，Mahmoudvand R，Zokaei M，et al. On the separability between signal and noise in singular spectrum analysis[J]. Fluctuation and Noise Letters，2012，11(2):1-11.

[2] Zhigljavsky A. Singular spectrum analysis for time series: introduction to this special issue[J]. Statistics and Its Interface，2010，3(3):255-258.

[3] Golyandina N，Zhigljavsky A. Singular spectrum analysis for time series[M]. [S.l.]: Springer，2013.

[4] Hassani H，Mahmoudvand R，Zokaei M. Separability and window length in singular spectrum analysis[J]. Comptes Rendus Mathematique，2011，349(17):987-990.

[5] Golyandina N. On the choice of parameters in singular spectrum analysis and related subspace-based methods[DB/OL]. [2010-07-20].http://arxiv.org/abs/1005. 4374.

[6] 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.

[7] 王益艷.基于特征均值的SVD信號(hào)去噪算法[J].計(jì)算機(jī)應(yīng)用與軟件，2012，29(5):121-123.

Wang Y Y. Mean value of eigenvalue-based SVD signal denoising algorithm[J]. Computer Applications and Software，2012，29 (5) :121-123. (in Chinese)

[8] 錢征文，程禮，李應(yīng)紅.利用奇異值分解的信號(hào)降噪方法[J].振動(dòng)，測(cè)試與診斷，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. (in Chinese)

[9] Stewart G W. Perturbation theory for the singular value decomposition[R]. College Park: Computer Science Department，University of Maryland，1998.

[10] Candès E J，Tao T. The power of convex relaxation: Near-optimal matrix completion[J]. IEEE Transactions on Information Theory，2010，56(5):2053-2080.

[11] Candes E J，Plan Y. Matrix completion with noise[J]. Proceedings of the IEEE，2010，98(6):925-936.

[12] Wright J，Ganesh A，Rao S，et al. Robust principal component analysis: exact recovery of corrupted low-rank matrices via convex optimization[C]∥The Twenty-Fourth Annual Conference on Neural Information Processing Systems. Vancouver:[s.n.]，2009:2080-2088.[13] Candès E J，Li X，Ma Y，et al. Robust principal component analysis？[J]. Journal of the ACM ，2011，58(3):1101-1137.

[14] 彭義剛，索津莉，戴瓊海，等.從壓縮傳感到低秩矩陣恢復(fù):理論與應(yīng)用[J].自動(dòng)化學(xué)報(bào)，2013，39(7):981-994.

Peng Y G，Suo J L，Dai Q H，et al. From compressed sensing to low-rank matrix recovery: theory and applications[J]. Acta Automatica Sinica，2013，39(7):981-994. (in Chinese)

[15] Zhou Z，Li X，Wright J，et al. Stable principal component pursuit[C]∥Proceedings of the IEEE International Symposium on Information Theory. Austin: IEEE，2010:1518-1522.

[16] Ganesh A，Lin Z，Wright J，et al. Fast algorithms for recovering a corrupted low rank matrix[C]∥The third IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing. Radisson, Aruba: IEEE，2009:213-216.

[17] Cai J F，Candès E J，Shen Z. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization，2010，20(4):1956-1982.

[18] Toh K C，Yun S. An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[J]. Pacific Journal of Optimization，2010，15(6):615-640.

[19] Lin Z，Chen M，Ma Y. The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[DB/OL].2010-09-26. http:// arXiv preprint arXiv: 1009.5055.

[20] 徐鋒，劉云飛.基于中值濾波-奇異值分解的膠合板拉伸聲發(fā)射信號(hào)降噪方法研究[J].振動(dòng)與沖擊，2011，30(12):135-140.

Xu Feng，Liu Yunfei. Noise reduction of acoustic emission signals in a plywood based on median filtering-singular value decomposition[J]. Journal of Vibration and Shock，2011，30(12):135-140. (in Chinese)

(責(zé)任編輯：劉芳)

Signal Denoising Method Based on Improve Singular Spectrum Analysis

DAI Hao-min，XU Ai-qiang，SUN Wei-chao

(Institute of Aircraft Detection and Application，College of Naval Aeronautical and Engineering，Yantai，Shandong 264001，China)

The traditional denoising method，such as wavelet threshold，is only valid for Gaussian noise，but is powerless for pulse jamming. Singular spectrum analysis developed in recent years can be a good filter at high SNR conditions of these two types of noises，but the process of noise reduction involves a certain subjective factors，and is subject to restrictions of matrix perturbation theory. Moreover，the ability of denoising will decrease with the lower SNR. For the above situation，an improved algorithm was proposed，which applied rank minimization theory to singular spectrum analysis. Simulation results show that denoising effect of the improved algorithm is obvious，which can maximize the reduction of the mean square error of the signals and improve signal to noise ratio; enhance versatility of singular spectrum analysis.

singular spectrum analysis; singular value decomposition; matrix perturbation theory; rank minimization theory

2014-04-01

國家部委預(yù)研基金資助項(xiàng)目(9140A27020212JB14311)

戴豪民(1982—)，男，博士生，E-mail：daihaomin521@163.com.

許愛強(qiáng)(1963—)，男，教授，博士生導(dǎo)師，E-mail：1397245858@qq.com.

TH 911

A

1001-0645(2016)07-0727-07

10.15918/j.tbit1001-0645.2016.07.013