基于矩陣秩最小化和統計修正的信號降噪方法研究

基于矩陣秩最小化和統計修正的信號降噪方法研究

李文峰1,2，許愛強1，戴豪民1，王豐3

(1.海軍航空工程學院科研部,山東煙臺264001；2.92635部隊山東青島266041； 3.海軍航空工程學院91206部隊,山東青島266108)

摘要：針對奇異值分解在信號降噪時有效秩的選擇問題，提出一種基于矩陣秩最小化和統計修正的信號降噪方法。首先，利用矩陣秩最小化理論將奇異值有效秩選擇問題轉化為秩的約束優化問題；然后，通過凸優化求解，得到干凈信號的Hankel矩陣，實現一次降噪；最后，根據奇異值子集標準差對干凈信號Hankel矩陣進行統計修正，進一步優化降噪效果。模擬信號和真實信號的實驗結果表明：該方法可以有效的消除脈沖干擾和高斯噪聲，能夠最大限度的降低信號均方誤差，提高信噪比，增強了奇異值分解在信號降噪中的通用性。

關鍵詞：奇異值分解；奇異值擾動理論；矩陣秩最小化；奇異值子集標準差

中圖分類號:TN911.7文獻標志碼：A

Signal denoising method based on matrix rank minimization and statistical modification

LIWen-feng1,2,XUAi-qiang1,DAIHao-min1,WANGFeng3(1. Research Department, Naval Aeronautical and Astronautical University, Yantai 264001, China;2. No. 92635 Unit of PLA,Qingdao 266041, China；3. No. 91206 Unit of PLA, Naval Aeronautical and Astronautical University, Qingdao 266108, China)

Abstract:Aiming at the selection problem of effective ranks in singular value decomposition for signal denoising, a signal denoising method based on matrix rank minimization and statistical modification was proposed. Firstly, the effective rank selection problem of singular value decomposition was transformed into a constrained optimization problem of rank by using the matrix rank minimization theory. Secondly, Hankel matrix of a clean signal was obtained with a convex optimization to realize the first noise reduction. At last, the statistical correction was performed for the clean signal’s Hankel matrix with subset standard deviation of a singular value to further optimize the noise reduction effect. The simulated signal and real signal test results showed that the method can effectively eliminate pulse noise and Gaussian noise; at the same time, the method can reduce the maximum signal mean square error and improve the signal-to-noise ratio; so the method can enhance the universality of singular value decomposition in signal denoising.

Key words:singular value decomposition (SVD); matrix perturbation theory; matrix rank minimization; subset standard deviation of singular value

信號在產生、采集和傳遞時往往伴有噪聲的干擾，因此信號預處理中的信號降噪是保證后續信號分析的重要前提。傳統的降噪方法如小波降噪和基于譜分析的濾波方法，其原理是利用信號與噪聲的頻譜不同達到降噪的效果。但是，當有用信號和噪聲頻帶混疊嚴重時，傳統的降噪方法則具有較大的局限性。近年來，奇異值分解(Singular Value Decomposition，SVD)在信號降噪問題上已經得到了廣泛應用，其中奇異值有效秩選擇問題一直是當前研究的熱點和難點。Donoho等[1]提出了奇異值優化硬閾值法，給出了硬閾值的優化選擇標準；趙學智等[2]提出了奇異值差分譜法，根據差分譜峰值位置可實現對有用分量個數的確定，王建國等[3]在奇異值差分譜的基礎上提出單邊極大值選擇方法，即在奇異值差分譜中從右至左，選擇第一個至少單邊與其相鄰峰值比較，差距絕對值最大的極大峰值的對應點位置，來確定重構信號的有效秩階次，錢征文等[4]采用信號快速傅里葉變換結果中主頻個數來確定奇異值數目,王樹青等[5]采用奇異值相對變化率確定信號降噪時奇異值的有效秩。

以上方法在高斯噪聲的背景下均取得了較好的降噪效果。但是也存在兩個問題：一是噪聲中如果包含稀疏大噪聲，比如脈沖信號，則上述方法效果欠佳；二是在奇異值個數選擇上，均是一次性截斷奇異值，容易造成了奇異值有效秩選擇的武斷性。針對上述問題，本文首先從矩陣奇異值擾動理論角度分析截斷奇異值方法的局限性；其次，利用矩陣秩最小化理論將有效秩選擇問題轉化為秩的約束優化問題，并通過凸優化求解得到干凈信號的Hankel矩陣。最后根據奇異值子集標準差對干凈信號Hankel矩陣進行統計修正，進一步優化降噪效果。實驗結果表明：相比于其他方法，該算法可以消除混合噪聲對信號的影響，增強了SVD在信號降噪中的可用性和有效性。

1奇異值分解降噪原理

(1)

式中，K=N-L+1，表示矩陣的嵌入維數，L是滑動窗口長度。矩陣Z是一個Hankel矩陣，即矩陣反對角線上元素相等。對矩陣Z進行奇異值分解，得：

(2)

式中，ui和vi分別是Z的左右奇異向量，d表示矩陣Z的秩，即非零奇異值個數，且d≤min(L,K)，λ1,λ2,…,λd是按降序排列的Z的奇異值。SVD降噪的本質就是通過選擇合適奇異值的個數，實現干凈信號矩陣X和噪聲矩陣E的分離，即Z=X+E。最后，對分離出的干凈信號矩陣X進行對角平均化，即可將該矩陣重構為長度為N的一維序列，重構序列的計算公式如下：

(3)

在利用SVD進行信號降噪時有兩個問題需要考慮，一是構造軌跡矩陣時的窗口長度L的選擇問題。文獻[6]從重構序列各成分的分離性和重構誤差的角度上推薦L取序列長度N的中值(如果N為偶數，取L=N/2)；文獻[7]研究了軌跡矩陣的奇異值隨窗口長度L的變化規律，從理論上證明了當L取N的中值時，對應的奇異值會取得最大值，同時該文獻也推薦在大多數情況下L取N的中值是最佳的選擇。二是分離干凈信號矩陣所需的奇異值個數，即有效秩的選擇問題。傳統方法大多從奇異值分布曲線的變化突點，即曲線“肘部”處截斷奇異值來確定干凈信號矩陣的秩，此類方法操作簡單，但是對脈沖噪聲或脈沖和高斯混合噪聲就會失效，接下來從矩陣奇異值擾動理論來說明截斷奇異值方法的失效原因。

2矩陣奇異值擾動理論

(4)

(5)

(6)

(7)

從式(7)中可以看出，在高信噪比的情況下，通過選擇前r個較大的奇異值逼近干凈信號矩陣是合理的。但是在低信噪比情況下或者信號受到稀疏大噪聲污染，雖然通過一次性截斷奇異值能夠移除部分噪聲，仍會有部分噪聲對前r個奇異值產生擾動，所以重構信號也將產生較大誤差。

3基于矩陣秩最小化的降噪方法

3.1矩陣秩最小化理論

Wright等[9]提出了魯棒主成分分析(Robust Principal Component Analysis，RPCA)算法，證明了通過解優化問題能夠從被污染矩陣中精確恢復出所需的低秩矩陣。因此，將奇異值有效秩選擇問題轉化為矩陣秩的約束優化問題，可以客觀地解決奇異值有效秩選擇問題。

當數據矩陣具有低秩性的結構[10-11]，并且只有少部分元素被噪聲嚴重污染，即噪聲是稀疏的，那么矩陣降噪問題可用如下模型來描述：

minrank(S)+γ‖W‖0

s.t.S+W=Z

(8)

min‖S‖*+γ‖W‖1

s.t.S+W=Z

(9)

這樣就將求解矩陣秩的非凸問題轉化為凸優化問題，使得該問題具有唯一解。

更為一般的情況，假定低秩矩陣S既被輕微高斯噪聲E所污染，也被稀疏大噪聲W所污染，并且E中元素滿足Ei,j～N(0,σ2)，‖E‖F≤δ，δ>0，那么矩陣降噪就可以用如下優化問題來描述[13]：

min‖S‖*+γ‖W‖1

s.t.‖Z-S-W‖F≤δ

(10)

目前，求解式(9)和式(10)最優化問題有多種算法，最為常用的是非精確增廣拉格朗日乘子法(Inexact Augmented Lagrange Multiplier，IALM)，算法的詳細描述及收斂性證明參見文獻[14]。

3.2仿真實驗

綜合上述分析，為去除信號中的噪聲干擾，將秩最小化理論應用于SVD中實現降噪的步驟如下：

(1)通過嵌入操作，將長度N的含噪序列ZN=[z1,z2,…,zN]T映射為軌跡矩陣Z∈RL×K，K=N-L+1，L=N/2。

(2)根據式(9)或式(10)，利用IALM算法恢復出低秩矩陣S。

(3)將低秩矩陣S對角平均化，轉化成所需的長度為N的干凈信號序列，最終實現信噪分離。

為了驗證降噪效果，采用信號均方誤差(Mean Square Error，MSE)和信噪比(Signal to Noise Ratio，SNR)來評價各種方法的優劣，其定義形式如下：

(11)

(12)

實驗1 考察以下仿真信號：

z(t)=s(t)+w(t)=

式中，w(t)是隨機添加的21個脈沖噪聲。

圖1　三種方法降噪效果 Fig.1 Denoising effect of three methods

文獻[1-5]采用的噪聲實例均是高斯噪聲，筆者經仿真實驗發現這些方法不能有效地去除脈沖噪聲。由于中值濾波(Median Filter，MF)具有抑制脈沖干擾的作用，所以采用中值濾波以及常見的小波閾值降噪與本文方法對比。各種方法的降噪效果如圖1所示。

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

從圖1可以看出，小波閾值降噪方法無法消除脈沖噪聲，中值濾波方法可以較好的抑制脈沖干擾，但是中值濾波受窗口長度的影響，抑制脈沖干擾效果并不穩定，過長的窗口長度容易使信號更加平滑，誤差較大，過短的窗口長度則不能有效抑制過于“密集”的脈沖干擾(如圖1a放大窗口所示)。相比較前兩種方法，本文方法在抑制脈沖干擾時效果最佳。表1是各種方法降噪后的MES和SNR，從表1也可以看出，基于矩陣秩最小化的SVD降噪方法表現特別優秀，幾乎完全恢復干凈信號，其誤差精度可以忽略不計。

實驗2考察以下仿真信號：

x(t)=s(t)+e(t)+w(t)=

式中，e(t)為信噪比是5 dB的高斯白噪聲，w(t)是隨機添加的17個脈沖噪聲。首先采用中值濾波方法濾除仿真信號的脈沖干擾，然后結合文獻[1,3-5]方法與矩陣秩最小化方法對剩余混合噪聲進行降噪效果對比。

圖2　文獻[3]、[4]和[5]的有效秩選擇方法 Fig.2 Effective rank selection method of literature [3],[4] and [5]

文獻[1]提出的奇異值優化硬閾值方法，通過計算優化硬閾值為24.669 6；文獻[3]提出的奇異值差分譜單邊極大值方法，從圖2(a)中計算得到所需的奇異值有效秩為5；文獻[4]提出用于恢復有用信號的奇異值個數等于信號功率譜主頻個數的兩倍，由于圖2(b)中功率譜主頻個數并不清晰可見，所以近似選擇奇異值有效秩為10；文獻[5]提出的奇異值相對變化率方法，從圖2(c)中得到所需的奇異值有效秩為15。

圖3　五種方法降噪效果 Fig.3 Denoising effect of five methods

從圖4和表2可以看出，在高斯噪聲和脈沖噪聲組成的混合噪聲背景下，矩陣秩最小化方法降噪效果優于其他四種降噪方法，可以完全濾除脈沖信號的干擾，但是由于高斯噪聲的影響，恢復出的信號與原始信號相比，還存在部分失真。針對上述問題，接下來采用奇異值子集標準差方法對矩陣秩最小化降噪后的信號進行統計修正，進一步優化降噪效果。

4基于奇異值子集標準差的統計修正降噪

4.1奇異值子集標準差統計修正方法

矩陣秩最小化降噪方法雖然可以濾除大部分混合噪聲，但仍存在部分高斯噪聲造成恢復出的信號失真。需要注意的是，此時恢復出的信號具有高信噪比的特點，因此可以根據奇異值分布曲線的“肘部”來選擇奇異值的個數。為了能夠更加清晰的觀察奇異值“肘部”處的變化情況，本文采用文獻[15]提出的奇異值子集標準差方法來確定奇異值個數，實現二次降噪，具體步驟如下：

(1)對秩最小化恢復出的干凈信號矩陣進行奇異值分解；

(2)將得到的奇異值按升序排列，即τ1≤τ2≤...≤τn；

(3)構建奇異值子集Γ1={τ1},Γ2={τ1,τ2},…,Γn={τ1,τ2,…,τn}；

(4)計算每個子集的標準差{σ1,σ2,…,σn}，并向下取整，即iσi=?σi,i=1,2,…,n，若kσk≠0，則計算終止。存儲從第k個一直到最后一個奇異值，然后將存儲的n-k+1個奇異值進行信號重構，得到最終的干凈信號。

采用奇異值子集標準差進行統計修正的詳細流程見圖4。

圖4　奇異值子集標準差統計修正方法流程 Fig.4 The process based on subset standard deviation of the singular value for statistical correction

4.2仿真分析

圖5是采用子集標準差方法對矩陣秩最小化方法恢復出的信號進行統計修正后的奇異值分布曲線。從圖5可以清晰的看出奇異值“肘部”的變化情況，從第467個奇異值一直到最后一個奇異值，其標準差均不為0，所以選擇467到512個奇異值進行信號重構。從圖6的重構信號可以看出，信號的失真程度已經變小。從表2也可以看出，在對混合噪聲進行降噪時，基于矩陣秩最小化結合統計修正的信號降噪方法要明顯優于其他四種降噪方法。

圖5　奇異值子集標準差選擇奇異值有效秩 Fig.5 The selection problem of effective rank

圖6　干凈信號S 0及矩陣秩最小化 和統計修正降噪信號S 6 Fig.6 The clean signal S 0 and statistical correction denoising signal S 6

評價指標MF+文獻[1]方法MF+文獻[3]方法MF+文獻[4]方法MF+文獻[5]方法矩陣秩最小化方法統計修正后方法MSE0.20650.44090.33340.39310.11370.0338SNR6.15012.85564.06973.35368.740914.0123

5實例分析

圖7a為電機軸承正常狀態下的振動加速度信號，信號中含有脈沖干擾。當電機軸承正常運行時，其振動信號應是平穩的隨機波形，所以可以認為圖中的尖峰是由異常運行情況導致，如果不消除這些脈沖干擾，將會對后續的信號分析及軸承故障診斷時的狀態分類產生較大影響。此外，本文提出的信號降噪方法同時也適用于解決圖像分析與處理、計算機視覺等其他領域。采用本文提出的基于矩陣秩最小化和統計修正方法對振動信號進行降噪后的信號如圖8(a)所示。圖7(b)、(c)和圖8(b)、(c)分別是對實測信號和降噪后信號進行經驗模態分解(Empirical Mode Decomposition，EMD)后得到的固有模態函數(Intrinsic Mode Function，IMF)，從圖中可以看出，原實測信號經EMD分解后產生11個imf分量和1個殘余分量，各模態之間產生混疊現象，而降噪后信號經EMD分解產生7個imf分量和1個殘余分量，能夠較好反映真實信號的成分。

圖7　電機軸承實測信號及其imf分量 Fig.7 The vibration signal of motor bearing and its imf components

圖8　本文方法降噪信號及其imf分量 Fig.8 The denoising signal by this paper method and its imf components

6結論

(1)從矩陣奇異值擾動理論角度分析了傳統奇異值分解方法在信號降噪時的局限性，提出利用矩陣秩最小化理論將奇異值有效秩選擇問題轉化為秩的約束優化問題，并通過凸優化求解，實現一次降噪。

(2)提出利用奇異值子集標準差方法將矩陣秩最小化恢復出的干凈信號矩陣進行統計修正，實現二次降噪。

(3)通過模擬信號和實測振動信號進行降噪仿真分析，驗證了本文方法的有效性，增強了奇異值分解在信號降噪中通用性。

參考文獻

[2]趙學智,葉彥邦,陳統堅.奇異值差分譜理論及其在車床主軸箱故障診斷中的應用[J]. 機械工程學報, 2010, 46(1):100-108.

ZHAO Xue-zhi, YE Yan-bang, CHEN Tong-jian. 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.

[3]王建國,李健,劉穎源.一種確定奇異值分解降噪有效秩階次的改進方法[J]. 振動與沖擊, 2014, 33(12):176-180.

WANG Jian-guo, LI Jian, LIU Yian-yuan. Animproved method for determining effective order rank of SVD denoising [J]. Journal of Vibration and Shock, 2014, 33(12):176-180.

[4]錢征文,程禮,李應紅.利用奇異值分解的信號降噪方法[J].振動,測試與診斷,2011, 31(4): 459-463.

QIAN Zhen-wen, CHENG Li, LI Ying-hong. Signal denoising method by means of SVD[J]. Journal of Vibration Measurement & Diagnosis, 2011, 31(4):459-463.

[5]王樹青,林裕裕,孟元棟,等. 一種基于奇異值分解技術的模型定階方法[J]. 振動與沖擊, 2012,31(15):87-91.

WANG Shu-qing, LIN Yu-yu, MENG Yuan-dong, et al. Model order determination based on singular valued decomposition[J]. Journal of Vibration and Shock, 2012,31(15):87-91.

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

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

[8]Hassani H, Xu Z Y, Zhigljavsky A. Singular spectrum analysis based on the perturbation theory[J]. Nonlinear Analysis : Real World Applications,2011,12:2752-2766.

[9]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, 2009: 2080-2088.

[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]彭義剛,索津莉,戴瓊海,等.從壓縮傳感到低秩矩陣恢復: 理論與應用[J].自動化學報, 2013, 39(7): 981-994.

PENG Yi-gang, SUO Jin-li, DAI Qiong-hai, et al. From compressed sensing to low-rank matrix recovery: theory and applications[J]. Acta Automatica Sinica, 2013, 39(7): 981-994.

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

[14]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.org/abs/1009.5055v3.

[15]Ashtiani M B, Shahrtash S M. Partial discharge de-noising employing adaptive singular value decomposition[J]. IEEE Transactions on Dielectrics and Electrical Insulation, 2014, 21(2): 775-782.