基于矩陣秩最小化和統(tǒng)計(jì)修正的信號降噪方法研究

基于矩陣秩最小化和統(tǒng)計(jì)修正的信號降噪方法研究

李文峰1,2，許愛強(qiáng)1，戴豪民1，王豐3

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

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

關(guān)鍵詞：奇異值分解；奇異值擾動理論；矩陣秩最小化；奇異值子集標(biāo)準(zhǔn)差

中圖分類號:TN911.7文獻(xiàn)標(biāo)志碼：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

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

以上方法在高斯噪聲的背景下均取得了較好的降噪效果。但是也存在兩個(gè)問題：一是噪聲中如果包含稀疏大噪聲，比如脈沖信號，則上述方法效果欠佳；二是在奇異值個(gè)數(shù)選擇上，均是一次性截?cái)嗥娈愔?，容易造成了奇異值有效秩選擇的武斷性。針對上述問題，本文首先從矩陣奇異值擾動理論角度分析截?cái)嗥娈愔捣椒ǖ木窒扌?；其次，利用矩陣秩最小化理論將有效秩選擇問題轉(zhuǎn)化為秩的約束優(yōu)化問題，并通過凸優(yōu)化求解得到干凈信號的Hankel矩陣。最后根據(jù)奇異值子集標(biāo)準(zhǔn)差對干凈信號Hankel矩陣進(jìn)行統(tǒng)計(jì)修正，進(jìn)一步優(yōu)化降噪效果。實(shí)驗(yàn)結(jié)果表明：相比于其他方法，該算法可以消除混合噪聲對信號的影響，增強(qiáng)了SVD在信號降噪中的可用性和有效性。

1奇異值分解降噪原理

(1)

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

(2)

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

(3)

在利用SVD進(jìn)行信號降噪時(shí)有兩個(gè)問題需要考慮，一是構(gòu)造軌跡矩陣時(shí)的窗口長度L的選擇問題。文獻(xiàn)[6]從重構(gòu)序列各成分的分離性和重構(gòu)誤差的角度上推薦L取序列長度N的中值(如果N為偶數(shù)，取L=N/2)；文獻(xiàn)[7]研究了軌跡矩陣的奇異值隨窗口長度L的變化規(guī)律，從理論上證明了當(dāng)L取N的中值時(shí)，對應(yīng)的奇異值會取得最大值，同時(shí)該文獻(xiàn)也推薦在大多數(shù)情況下L取N的中值是最佳的選擇。二是分離干凈信號矩陣所需的奇異值個(gè)數(shù)，即有效秩的選擇問題。傳統(tǒng)方法大多從奇異值分布曲線的變化突點(diǎn)，即曲線“肘部”處截?cái)嗥娈愔祦泶_定干凈信號矩陣的秩，此類方法操作簡單，但是對脈沖噪聲或脈沖和高斯混合噪聲就會失效，接下來從矩陣奇異值擾動理論來說明截?cái)嗥娈愔捣椒ǖ氖г颉?/p>

2矩陣奇異值擾動理論

(4)

(5)

(6)

(7)

從式(7)中可以看出，在高信噪比的情況下，通過選擇前r個(gè)較大的奇異值逼近干凈信號矩陣是合理的。但是在低信噪比情況下或者信號受到稀疏大噪聲污染，雖然通過一次性截?cái)嗥娈愔的軌蛞瞥糠衷肼暎詴胁糠衷肼晫η皉個(gè)奇異值產(chǎn)生擾動，所以重構(gòu)信號也將產(chǎn)生較大誤差。

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

3.1矩陣秩最小化理論

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

當(dāng)數(shù)據(jù)矩陣具有低秩性的結(jié)構(gòu)[10-11]，并且只有少部分元素被噪聲嚴(yán)重污染，即噪聲是稀疏的，那么矩陣降噪問題可用如下模型來描述：

minrank(S)+γ‖W‖0

s.t.S+W=Z

(8)

min‖S‖*+γ‖W‖1

s.t.S+W=Z

(9)

這樣就將求解矩陣秩的非凸問題轉(zhuǎn)化為凸優(yōu)化問題，使得該問題具有唯一解。

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

min‖S‖*+γ‖W‖1

s.t.‖Z-S-W‖F(xiàn)≤δ

(10)

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

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

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

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

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

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

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

(11)

(12)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4基于奇異值子集標(biāo)準(zhǔn)差的統(tǒng)計(jì)修正降噪

4.1奇異值子集標(biāo)準(zhǔn)差統(tǒng)計(jì)修正方法

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

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

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

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

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

采用奇異值子集標(biāo)準(zhǔn)差進(jìn)行統(tǒng)計(jì)修正的詳細(xì)流程見圖4。

圖4　奇異值子集標(biāo)準(zhǔn)差統(tǒng)計(jì)修正方法流程 Fig.4 The process based on subset standard deviation of the singular value for statistical correction

4.2仿真分析

圖5是采用子集標(biāo)準(zhǔn)差方法對矩陣秩最小化方法恢復(fù)出的信號進(jìn)行統(tǒng)計(jì)修正后的奇異值分布曲線。從圖5可以清晰的看出奇異值“肘部”的變化情況，從第467個(gè)奇異值一直到最后一個(gè)奇異值，其標(biāo)準(zhǔn)差均不為0，所以選擇467到512個(gè)奇異值進(jìn)行信號重構(gòu)。從圖6的重構(gòu)信號可以看出，信號的失真程度已經(jīng)變小。從表2也可以看出，在對混合噪聲進(jìn)行降噪時(shí)，基于矩陣秩最小化結(jié)合統(tǒng)計(jì)修正的信號降噪方法要明顯優(yōu)于其他四種降噪方法。

圖5　奇異值子集標(biāo)準(zhǔn)差選擇奇異值有效秩 Fig.5 The selection problem of effective rank

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

評價(jià)指標(biāo)MF+文獻(xiàn)[1]方法MF+文獻(xiàn)[3]方法MF+文獻(xiàn)[4]方法MF+文獻(xiàn)[5]方法矩陣秩最小化方法統(tǒng)計(jì)修正后方法MSE0.20650.44090.33340.39310.11370.0338SNR6.15012.85564.06973.35368.740914.0123

5實(shí)例分析

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

圖7　電機(jī)軸承實(shí)測信號及其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結(jié)論

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

(2)提出利用奇異值子集標(biāo)準(zhǔn)差方法將矩陣秩最小化恢復(fù)出的干凈信號矩陣進(jìn)行統(tǒng)計(jì)修正，實(shí)現(xiàn)二次降噪。

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

參考文獻(xiàn)

[2]趙學(xué)智,葉彥邦,陳統(tǒng)堅(jiān).奇異值差分譜理論及其在車床主軸箱故障診斷中的應(yīng)用[J]. 機(jī)械工程學(xué)報(bào), 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ìn)方法[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]錢征文,程禮,李應(yīng)紅.利用奇異值分解的信號降噪方法[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]王樹青,林裕裕,孟元棟,等. 一種基于奇異值分解技術(shù)的模型定階方法[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]彭義剛,索津莉,戴瓊海,等.從壓縮傳感到低秩矩陣恢復(fù): 理論與應(yīng)用[J].自動化學(xué)報(bào), 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.