奇異值分解降噪的改進(jìn)方法

張 磊 彭偉才 原春暉 劉 彥

中國(guó)艦船研究設(shè)計(jì)中心船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430064

0 引 言

在船舶振動(dòng)噪聲信號(hào)的測(cè)試過(guò)程中，總存在著由現(xiàn)場(chǎng)測(cè)試環(huán)境及儀器本身產(chǎn)生的誤差源［1-4］，其中一個(gè)不可避免的誤差源就是隨機(jī)噪聲。測(cè)量信號(hào)中的噪聲水平直接決定了分析的可靠性，為此，必須消除或者最小化信號(hào)中的噪聲干擾。

作為信號(hào)降噪的有效方法之一，奇異值分解（SVD）已被廣泛應(yīng)用于許多工程領(lǐng)域，如信號(hào)濾波和矩陣秩的估計(jì)。將受到噪聲污染的信號(hào)構(gòu)造成Hankel矩陣，然后對(duì)其進(jìn)行奇異值分解，將包含信號(hào)特征的矩陣分解到一系列奇異值和奇異值矢量對(duì)應(yīng)的子空間中，是對(duì)Hankel矩陣進(jìn)行噪聲過(guò)濾的一種非線性濾波方法。從SVD降噪的基本原理來(lái)看，其關(guān)鍵就是確定Hankel矩陣的有效奇異值數(shù)目，奇異值數(shù)目選擇不當(dāng)，會(huì)極大地影響降噪效果。而目前用于確定奇異值數(shù)目的方法，如奇異值曲線、奇異熵增量［5］、信噪比經(jīng)驗(yàn)［6］等，都不能明確地給出有效階次，往往只能依靠經(jīng)驗(yàn)選取。

本文將通過(guò)仿真分析與實(shí)驗(yàn)測(cè)試來(lái)研究奇異值數(shù)目變化時(shí)噪聲對(duì)信號(hào)的干擾影響。

1 奇異值分解降噪基本原理

對(duì)于一個(gè)受噪聲污染的離散信號(hào) X={x1，x2，…，xL}，構(gòu)造成m×n(m≤n)維的Hankel矩陣為：

式中，A為Hankel矩陣；m為嵌入維數(shù)，并且滿(mǎn)足m+n-1=L。

對(duì)Hankel矩陣進(jìn)行奇異值分解可得到：

式中，U為m×m維正交矩陣；V為n×n維正交矩陣；Σ為m×n維矩陣。主對(duì)角元素為矩陣的奇異值，并且從大到小排列。

矩陣A為由受噪聲污染的信號(hào)構(gòu)成的Hankel矩陣，可以表示為未受噪聲污染的信號(hào)子空間和噪聲子空間之和：

保留對(duì)角矩陣的前k個(gè)有效奇異值，將其他的奇異值設(shè)置為零，利用奇異值分解的逆過(guò)程得到重構(gòu)矩陣。一般來(lái)說(shuō)，這時(shí)的重構(gòu)矩陣不再是Hankel矩陣的形式。為了得到降噪后的信號(hào)，需要對(duì)重構(gòu)矩陣中的反對(duì)角元素采用下式進(jìn)行平均：

2 改進(jìn)方法

由實(shí)測(cè)經(jīng)驗(yàn)可知，當(dāng)實(shí)際測(cè)量環(huán)境很好時(shí)，測(cè)量得到的信號(hào)一般為光滑曲線；而當(dāng)受到外界的隨機(jī)噪聲干擾時(shí)，測(cè)量得到的信號(hào)中就會(huì)含有大量的“毛刺”。

利用噪聲污染信號(hào)構(gòu)造Hankel矩陣進(jìn)行奇異值分解降噪，就是對(duì)含噪信號(hào)進(jìn)行逼近、剔除毛刺的過(guò)程。經(jīng)過(guò)對(duì)大量仿真結(jié)果的研究，發(fā)現(xiàn)去噪結(jié)果中剩余噪聲對(duì)信號(hào)的干擾影響能夠通過(guò)信號(hào)中“毛刺”的數(shù)量及大小來(lái)判斷，這與實(shí)測(cè)經(jīng)驗(yàn)相吻合。當(dāng)選取的奇異值數(shù)目較小時(shí)，大部分噪聲都被剔除掉了，但也損失了大量有用信號(hào)。隨著奇異值數(shù)目的逐漸增大，有用信號(hào)信息趨于完整，但噪聲的干擾也會(huì)逐漸增加。當(dāng)奇異值數(shù)目達(dá)到某一值時(shí)，降噪后的信號(hào)既保留了大部分有用信息，也剔除掉了大部分噪聲，這就是要選取的最佳奇異值數(shù)目。當(dāng)奇異值數(shù)目再增加時(shí)，基于隨機(jī)噪聲對(duì)信號(hào)干擾的特點(diǎn)，降噪后的信號(hào)中會(huì)出現(xiàn)大量“毛刺”，波動(dòng)明顯，在數(shù)學(xué)意義上可表示為有大量極值點(diǎn)出現(xiàn)。本文根據(jù)降噪后信號(hào)極值點(diǎn)（均指極大值點(diǎn)）數(shù)量的變化，找到了突變點(diǎn)，可清晰地判斷應(yīng)該選取的有效奇異值數(shù)目。

本文信噪比計(jì)算所用的公式為：

式中，SNR為信噪比，dB；n(i)為含有噪聲的信號(hào)；x(i)為真實(shí)信號(hào)；N為信號(hào)長(zhǎng)度。

3 仿真分析

3.1 噪聲干擾對(duì)奇異值的影響

經(jīng)過(guò)對(duì)無(wú)噪聲理想信號(hào)和受噪聲污染信號(hào)的研究，發(fā)現(xiàn)一般由無(wú)噪聲理想信號(hào)構(gòu)造的Hankel矩陣的大部分奇異值為零。根據(jù)非零奇異值的數(shù)目，可以很容易地判斷矩陣的有效階次。因噪聲具有隨機(jī)和不相關(guān)的特點(diǎn)，因此，由受隨機(jī)噪聲污染信號(hào)構(gòu)造的Hankel矩陣呈列滿(mǎn)秩或行滿(mǎn)秩狀態(tài)（取決于行和列哪個(gè)維數(shù)更?。?。

現(xiàn)取如圖1所示的無(wú)噪聲理想信號(hào)和噪聲污染信號(hào)進(jìn)行仿真分析。無(wú)噪聲理想信號(hào)clean signal＝sin（3×t）＋sin（5×t）＋sin（8×t），單位為 V，其中 t為時(shí)間間隔為0.01的從0～2變化的時(shí)間序列。在無(wú)噪聲理想信號(hào)中加入信噪比約為6dB（通過(guò)式（5）計(jì)算得到）的高斯白噪聲，得到了圖1中的噪聲污染信號(hào)，噪聲干擾明顯。

圖1 無(wú)噪聲理想信號(hào)與噪聲污染信號(hào)Fig.1 Clean signal and noisy signal

對(duì)圖1中的兩個(gè)信號(hào)取維數(shù)m＝80，n＝122構(gòu)造Hankel矩陣（本文在此不討論如何選取的維數(shù)，具體方法可參考文獻(xiàn)［8］），進(jìn)行奇異值分解，得到奇異值序列圖如圖2所示（只顯示了80個(gè)奇異值中的前30個(gè)）。

圖2 無(wú)噪聲理想信號(hào)和噪聲污染信號(hào)的奇異值序列圖Fig.2 Singular value curves of clean signal and noisy signal

對(duì)比圖2中的兩條奇異值曲線可知，每個(gè)奇異值受到噪聲干擾的程度不同：較大的奇異值受噪聲干擾的影響較小，而較小的奇異值受噪聲干擾的影響則較大。由于由受噪聲污染信號(hào)構(gòu)造的Hankel矩陣的每個(gè)奇異值都是由“有用信號(hào)”和“噪聲信號(hào)”兩部分組成，實(shí)際選取奇異值數(shù)目就是一個(gè)對(duì)有用信號(hào)和噪聲信號(hào)進(jìn)行取舍的過(guò)程，所以，就要考慮在所選取的數(shù)目對(duì)應(yīng)的那個(gè)奇異值中，有用信號(hào)和噪聲信號(hào)的比例關(guān)系。

3.2 方法對(duì)比

3.2.1 奇異值曲線和奇異熵增量

奇異值曲線和奇異熵增量判斷有效奇異值數(shù)目的方法，分別是以奇異值擬合曲線和奇異熵增量曲線的拐點(diǎn)來(lái)作為閾值設(shè)置的依據(jù)。實(shí)際計(jì)算結(jié)果表明，兩種曲線除幅值有所變化外，形狀并沒(méi)有多大區(qū)別。

對(duì)圖1中受噪聲污染的信號(hào)取維數(shù)m＝80，n＝122構(gòu)造Hankel矩陣，進(jìn)行奇異值分解，得到奇異值序列圖（圖3）和奇異熵增量隨奇異值數(shù)目變化的曲線圖（圖4）。通過(guò)對(duì)比，證實(shí)上述結(jié)論是正確的。

由于降噪效果往往對(duì)奇異值的數(shù)目很敏感，因此，奇異值數(shù)目的微小差別都將導(dǎo)致信號(hào)去噪效果有較大差異。上述方法存在的主要問(wèn)題就是拐點(diǎn)不夠清晰，即使存在拐點(diǎn)，選取對(duì)應(yīng)的奇異值數(shù)目一般也不能得到最優(yōu)的降噪結(jié)果。采用上述兩種方法，并通過(guò)觀察圖3和圖4，可以選取奇異值的數(shù)目為5。

圖3 奇異值序列圖（噪聲污染仿真信號(hào)）Fig.3 Singular value curve（noisy simulation signal）

圖4 奇異熵增量隨奇異值數(shù)目變化的曲線圖（噪聲污染仿真信號(hào)）Fig.4 Increment of the singularity entropy versus the number of the singular values（noisy simulation signal）

當(dāng)選取的奇異值數(shù)目為5時(shí)，通過(guò)計(jì)算，得到降噪后的信噪比為18.5 dB。去噪后的信號(hào)與無(wú)噪聲理想信號(hào)的對(duì)比如圖5所示。由圖可見(jiàn)，信號(hào)中含有很多“毛刺”，尤其是后半段受噪聲影響較嚴(yán)重。

圖5 奇異值數(shù)目為5對(duì)應(yīng)的降噪信號(hào)和無(wú)噪聲理想信號(hào)Fig.5 Noise elimination signal with 5 singular values and clean signal

3.2.2 改進(jìn)方法——極值點(diǎn)數(shù)量

現(xiàn)利用改進(jìn)的方法，采用與第3.2.1小節(jié)中相同的維數(shù)構(gòu)造Hankel矩陣，進(jìn)行奇異值分解。通過(guò)計(jì)算得到的信號(hào)降噪后極值點(diǎn)數(shù)量隨選取奇異值數(shù)目變化的關(guān)系如圖6所示。極值點(diǎn)數(shù)量突變非常明顯，要選取的有效奇異值數(shù)目為4。

圖6 極值點(diǎn)數(shù)量隨奇異值數(shù)目變化圖（噪聲污染仿真信號(hào)）Fig.6 Number of the extremum points versus the number of singular values（noisy simulation signal）

當(dāng)奇異值數(shù)目為4時(shí)，通過(guò)計(jì)算得到降噪后的信噪比為 22.3 dB，相比第 3.2.1小節(jié)中用兩種方法選取奇異值數(shù)目為5時(shí)的信噪比（18.5 dB），有明顯的提高。降噪信號(hào)和無(wú)噪聲理想信號(hào)的對(duì)比如圖7所示。由圖可見(jiàn)，其明顯優(yōu)于圖5中的降噪效果。

圖7 奇異值數(shù)目為4對(duì)應(yīng)的降噪信號(hào)和無(wú)噪聲理想信號(hào)Fig.7 Noise elimination signal with 4 singular values and clean signal

當(dāng)奇異值數(shù)目小于最佳的選項(xiàng)4時(shí)，為使效果更明顯，例如，選取1時(shí)，其降噪后的信噪比為6.3 dB，相比奇異值數(shù)目為 4時(shí)的信噪比 22.3 dB，有較大的衰減。去噪后的信號(hào)與無(wú)噪聲理想信號(hào)的對(duì)比如圖8所示。由圖可見(jiàn)，降噪信號(hào)中沒(méi)有“毛刺”，但形狀有明顯的改變，第2和第3個(gè)峰值已基本被過(guò)濾掉。

圖8 奇異值數(shù)目為1對(duì)應(yīng)的降噪信號(hào)和無(wú)噪聲理想信號(hào)Fig.8 Noise elimination signal with 1 singular value and clean signal

圖9所示為奇異值數(shù)目從1～8變化時(shí)信號(hào)降噪后的信噪比曲線圖。由圖可見(jiàn)，奇異值數(shù)目為4時(shí)，信噪比最高；當(dāng)奇異值數(shù)目大于或者小于4時(shí)，信噪比都有所下降，故4為最優(yōu)選擇。

圖9 降噪信號(hào)信噪比隨選取奇異值數(shù)目變化的曲線圖（噪聲污染仿真信號(hào)）Fig.9 Signal to noise ratio versus the number of the singular values（noisy simulation signal）

綜上所述，Hankel矩陣有效奇異值數(shù)目的確定至關(guān)重要，選取不當(dāng)，就會(huì)造成結(jié)果偏差較大。當(dāng)奇異值數(shù)目選取過(guò)小時(shí)，因奇異值的組成中主要是有用信號(hào)，噪聲的干擾可以忽略不計(jì)，所以降噪后的信號(hào)中沒(méi)有噪聲引起的明顯波動(dòng)。但由于閾值選取過(guò)大，很多包含在奇異值中的有用信號(hào)也被剔除掉了，表現(xiàn)為有用信號(hào)的損失。當(dāng)選取最佳的奇異值數(shù)目時(shí)，不僅過(guò)濾掉了大量的噪聲，而且還最大程度地保留了原始信號(hào)的信息，與原始信號(hào)吻合很好。所以，在降噪后剩余噪聲對(duì)信號(hào)沒(méi)有明顯影響的情況下，選取的奇異值數(shù)目越大，信噪比越高，就越是接近不受噪聲污染的原始信號(hào)。當(dāng)奇異值數(shù)目選取過(guò)大時(shí)，降噪后的信號(hào)中仍然含有大量的噪聲。例如，奇異值數(shù)目選取為5時(shí)，與最佳的奇異值數(shù)目4相比，雖然僅僅只大1，但信號(hào)中卻出現(xiàn)了大量的“毛刺”，噪聲干擾很明顯。使用本文的方法，根據(jù)信號(hào)降噪后極值點(diǎn)數(shù)量的變化，能夠清晰地反映剩余噪聲對(duì)信號(hào)的影響，從而確定應(yīng)該選取的最佳奇異值數(shù)目。

4 實(shí)驗(yàn)驗(yàn)證

4.1 實(shí)測(cè)頻響函數(shù)加高斯白噪聲

文獻(xiàn)［9］中提到，經(jīng)過(guò)對(duì)分別使用頻響函數(shù)（FRF）和脈沖響應(yīng)函數(shù)（IRF）構(gòu)造的Hankel矩陣進(jìn)行奇異值分解降噪的大量研究發(fā)現(xiàn)，使用FRF得到的奇異值曲線沒(méi)有明顯的拐點(diǎn)，不容易找到合理的閾值，在判斷奇異值分解的有效秩時(shí)，使用IRF更具優(yōu)越性。但是，使用IRF對(duì)奇異值進(jìn)行分解降噪，就需要把FRF先進(jìn)行逆傅立葉變換，待奇異值分解降噪處理后，再通過(guò)傅立葉變換將IRF變換成FRF。

根據(jù)降噪后信號(hào)極值點(diǎn)數(shù)量的突變來(lái)判斷矩陣的有效奇異值數(shù)目，彌補(bǔ)了上述方法的不足，即利用FRF構(gòu)造Hankel矩陣進(jìn)行奇異值分解降噪仍能判斷出最佳的奇異值數(shù)目，避免了正、逆傅立葉變換增加的計(jì)算量。

本文通過(guò)實(shí)驗(yàn)測(cè)得了頻率范圍為0～500 Hz的FRF信號(hào)，分辨率為1 Hz。由于實(shí)驗(yàn)時(shí)隨機(jī)噪聲干擾不大，為了更明確地展示改進(jìn)方法的有效性，對(duì)該實(shí)測(cè)信號(hào)加入了信噪比為15 dB的高斯白噪聲。加入了高噪聲污染的實(shí)測(cè)信號(hào)如圖10所示。

圖10 實(shí)測(cè)頻響函數(shù)信號(hào)與實(shí)測(cè)頻響函數(shù)信號(hào)加噪聲Fig.10 Measured FRF signal and measured FRF signal with additive noise

取m＝210，n＝292構(gòu)造Hankel矩陣并進(jìn)行奇異值分解，依照第3.2.1小節(jié)中的方法，得到奇異值序列圖如圖11所示。通過(guò)觀察圖11可發(fā)現(xiàn)，奇異值序列圖的拐點(diǎn)處呈“弧狀過(guò)渡區(qū)”，有效奇異值數(shù)目的選取不明顯，無(wú)法作出判斷。而依照第3.2.2小節(jié)中的改進(jìn)方法，計(jì)算得到的極值點(diǎn)數(shù)量隨選取奇異值數(shù)目變化的關(guān)系如圖12所示，發(fā)現(xiàn)29為最優(yōu)選擇。圖13所示為降噪后信號(hào)與實(shí)測(cè)信號(hào)的對(duì)比圖，其降噪后的信噪比為20.9 dB。

圖11 奇異值序列圖（頻響函數(shù)）Fig.11 Curve of the singular values（FRF）

圖12 極值點(diǎn)數(shù)量隨奇異值數(shù)目變化的關(guān)系圖（頻響函數(shù)）Fig.12 Number of the extremum points versus the number of singular values（FRF）

圖13 實(shí)測(cè)信號(hào)與奇異值數(shù)目為29對(duì)應(yīng)的降噪信號(hào)Fig.13 Measured signal and noise elimination signal with 29 singular values

圖14所示為奇異值數(shù)目從21～38變化時(shí)信號(hào)降噪后的信噪比曲線圖。由圖可見(jiàn)，奇異值數(shù)目為29時(shí)，信噪比最高；當(dāng)奇異值數(shù)目大于或者小于29時(shí)，信噪比都有所下降，故29為最優(yōu)選擇。

圖14 降噪信號(hào)信噪比隨選取奇異值數(shù)目變化的曲線圖（頻響函數(shù)）Fig.14 Signal to noise ratio versus the number of the singular values（FRF）

4.2 實(shí)測(cè)受噪聲干擾自功率譜

圖15為受噪聲干擾明顯的實(shí)測(cè)加速度自功率譜，頻率范圍為 0～499 Hz，分辨率為 1 Hz，同樣采用第3.2.2中的改進(jìn)方法進(jìn)行處理。取 m＝200，n＝301構(gòu)造Hankel矩陣并進(jìn)行奇異值分解，通過(guò)計(jì)算得到的降噪后信號(hào)極值點(diǎn)數(shù)量隨奇異值數(shù)目變化的關(guān)系如圖16所示。由圖可見(jiàn)，奇異值數(shù)目為11時(shí)為最優(yōu)選擇，其降噪后的信號(hào)如圖17所示。

圖15 實(shí)測(cè)受噪聲干擾的自功率譜信號(hào)Fig.15 Measured noisy autospectrum signal

圖16 極值點(diǎn)數(shù)量隨奇異值數(shù)目變化的關(guān)系圖（自功率譜）Fig.16 Number of the extremum points versus the number of singular values（autospectrum）

圖17 奇異值數(shù)目為11對(duì)應(yīng)的自功率譜降噪信號(hào)Fig.17 Noise elimination autospectrum signal with 11 singular values

5 結(jié) 語(yǔ)

目前用于確定奇異值數(shù)目的方法，如奇異值曲線、奇異熵增量、信噪比經(jīng)驗(yàn)等，都不能明確給出有效的階次，往往只能依賴(lài)經(jīng)驗(yàn)進(jìn)行選取。本文提出了根據(jù)降噪后信號(hào)中極值點(diǎn)數(shù)量的變化來(lái)確定與最優(yōu)降噪結(jié)果對(duì)應(yīng)的奇異值數(shù)目的方法。仿真及實(shí)驗(yàn)結(jié)果表明，該方法準(zhǔn)確、有效。在機(jī)械設(shè)備、船體結(jié)構(gòu)等多種噪聲振動(dòng)工程測(cè)試領(lǐng)域，尤其是在面臨特殊的測(cè)試環(huán)境（例如，海浪海風(fēng)的影響）時(shí)，利用該方法對(duì)諸如FRF、自功率譜等測(cè)試到的信號(hào)進(jìn)行降噪處理能夠得到最優(yōu)的降噪效果，使測(cè)試信號(hào)具有更高的可靠性，可為正確分析后續(xù)信號(hào)提供保障。

［1］MARUDACHALAM K，WICKS A L.An attempt to quantify the errors in the experimental modal analysis［C］//The 9th International Modal Analysis Conference.Florence，1991：1522-1527.

［2］JUNG H，EWINS D J.On the use of simulated“experimental”data for evaluation of modal analysis methods［C］//Proceedings of the 10th International Modal Analysis Conference，1992：421-429.

［3］MITCHELL L D.Modal testing methods：quality，quantity and the unobtainable［J］.Sound and Vibration，1994，28（11）：10-17.

［4］FAHEY S O，WICKS A L.Noise sources in mechanical measurements［J］.Experimental Techniques ，2000，24（2）：40-43.

［5］楊文獻(xiàn)，姜節(jié)勝.機(jī)械］［J］.機(jī)械工程學(xué)報(bào)，2000，36（12）：9-13.YANG W X，JIANG J S.Study on the singular entropy of mechanical signal［J］.Journal of Mechanical Engineering，2000，36（12）：9-13.

［6］呂永樂(lè)，郎榮玲，梁家誠(chéng).基于信噪比經(jīng)驗(yàn)值的奇異值分解濾波門(mén)限確定［J］.計(jì)算機(jī)應(yīng)用研究，2009，26（9）：3253-3255.LV Y L，LANG R L，LIANG J C.Decision of threshold for singular value decomposition filter based on SNR’s empirical value［J］.Application Research of Computers，2009，26（9）：3253-3255.

［7］GOLUB G H，VAN C F，LOAN V V.Matrix computations［M］.3rd Edition.Maryland：John Hopkins Studues in Mathematical Sciences，1996.

［8］修春波，劉向東，張宇河.相空間重構(gòu)延遲時(shí)間與嵌入維數(shù)的選擇［J］.北京理工大學(xué)學(xué)報(bào)：自然科學(xué)版，2003，23（2）：219-224.XIU C B，LIU X D，ZHANG Y H.Selection of embedding dimension and delay time in the phase space reconstruction［J］.Journal of Beijing Institute of Technology（Natural Science Edition），2003，23（2）：219-224.

［9］SANLITURK K Y，CAKAR O.Noise elimination from measured frequency response functions［J］.Mechanical Systems and Signal Processing，2005，19（3）：615-631.