1、基于奇異值分解的信號(hào)消噪技術(shù)摘 要模態(tài)參數(shù)識(shí)別是從結(jié)構(gòu)不同位置的動(dòng)力響應(yīng)信號(hào)中提取出結(jié)構(gòu)的模態(tài)參數(shù)，即：從動(dòng)力測(cè)試響應(yīng)信號(hào)數(shù)據(jù)中確定結(jié)構(gòu)的模態(tài)參數(shù)(模態(tài)振型、固有頻率和阻尼比)。每一個(gè)結(jié)構(gòu)都有其固有的模態(tài)參數(shù)，并且如果結(jié)構(gòu)動(dòng)力特性發(fā)生變化了，那么結(jié)構(gòu)的模態(tài)參數(shù)也將發(fā)生相應(yīng)的變化。顯見，結(jié)構(gòu)的模態(tài)參數(shù)識(shí)別是非常重要的，為診斷結(jié)構(gòu)健康狀況提供了依據(jù)。基于輸出的模態(tài)參數(shù)識(shí)別方法利用的信息主要是系統(tǒng)的自由振動(dòng)信號(hào)，要獲得自由振動(dòng)信號(hào)首先需獲得結(jié)構(gòu)的響應(yīng)信號(hào)。由于環(huán)境激勵(lì)的不充分和噪聲等干擾因素的存在，導(dǎo)致信號(hào)測(cè)試信號(hào)不能直接用于參數(shù)辨識(shí)，需要對(duì)信號(hào)進(jìn)行消噪處理。即從大量背景噪聲中提取出可用于模態(tài)參數(shù)辨

2、識(shí)的有用信號(hào)成分,剔除干擾因素，提取有用信息。此時(shí)，信號(hào)消噪技術(shù)研究變得尤為重要。 本文采用了一種將Hankel矩陣和奇異值分解相結(jié)合的消噪方法。該方法首先對(duì)測(cè)量信號(hào)構(gòu)造的Hankel矩陣進(jìn)行奇異值分解，再利用測(cè)量信號(hào)快速傅立葉變換結(jié)果中主頻率的個(gè)數(shù)來確定有效秩階次，接著通過消噪信號(hào)的信噪比和均方差大小確定重構(gòu)矩陣結(jié)構(gòu)，最后通過反對(duì)角線平均法得到消噪后的信號(hào)數(shù)據(jù)。通過數(shù)值仿真，對(duì)不同信號(hào)進(jìn)行定秩和消噪，從結(jié)果可以知道這種方法具有較好的消噪效果。關(guān)鍵詞：信號(hào)消噪；奇異值分解；快速傅立葉變換；信噪比；均方差 A Method for Noise Reduction Based on Singula

3、r Value DecompositionAbstract Accurate estimate of the modal parameters of an offshore structure is crucial to many practical engineering issues, such as finite element (FE) model updating and validation, damage detection, etc. Modal parameter identification method uses the the response signal of st

4、ructure ,but actual response signal often contains a lot of noise, which will affect the accuracy of signal recognition. The test signal de-noising processing is an important step in signal processing. Using Singular Value Decomposition（SVD）of constructed Hankel matrix by measured signal is an effec

5、tive method for eliminating the random noise. The key is to choose the rank of the Hankel matrix and determine the structure of the reconstruction matrix. In this paper, it is using the number of the main frequency in the result of using signal fast Fourier transform to determine the rank of the Han

6、kel matrix, and through SNR(Signal to Noise Ratio) and MSE(Mean Square Error) to determine reconstruction matrix structure.Simulation and experiment validated this method. The results shows that the number of rank is double of the main frequency, and the best lines of reconstruction matrix is half o

7、f the length of the signal data. You can easy to choose the rank of the matrix and get a better noise elimination result.Keywords:Signal de-noising; Singular value decomposition; Fast Fourier transform; Signal to noise ratio; Mean square error目 錄1 引言12 SVD分解消噪理論52.1 Hankel矩陣52.2 SVD分解的基本理論62.3 對(duì)測(cè)量信號(hào)

8、進(jìn)行SVD分解63 有效秩階次和重構(gòu)矩陣結(jié)構(gòu)的確定83.1 有效秩階次的確定83.2 重構(gòu)矩陣結(jié)構(gòu)的確定144 消噪后的信號(hào)重構(gòu)165 數(shù)值仿真175.1 Matlab仿真結(jié)果分析175.2 Matlab程序236 結(jié)束語26參考文獻(xiàn)27基于奇異值分解的信號(hào)消噪技術(shù)1 引言 隨著社會(huì)的發(fā)展，人類社會(huì)對(duì)石油的需求日益提高，海上采油區(qū)域不斷擴(kuò)大，有越來越多的海洋平臺(tái)建造并投入使用，而這些海洋平臺(tái)結(jié)構(gòu)在復(fù)雜的服役環(huán)境中將受到設(shè)計(jì)載荷的作用以及各種突發(fā)性外在因素的影響而面臨結(jié)構(gòu)的損傷積累的問題，從而使結(jié)構(gòu)的安全受到威脅。大型工程結(jié)構(gòu)一旦出現(xiàn)事故，所帶來的命和財(cái)產(chǎn)損失將是巨大的，對(duì)社會(huì)的影響更是深遠(yuǎn)和難

9、以估量的。特別對(duì)于海洋平臺(tái)而言，其結(jié)構(gòu)復(fù)雜，造價(jià)昂貴，一旦發(fā)生事故，不僅會(huì)對(duì)海洋環(huán)境造成很大的污染，還會(huì)帶來不可估量的經(jīng)濟(jì)損失和人員傷亡，造成不好的社會(huì)政治影響。眾所周知，海洋平臺(tái)結(jié)構(gòu)長(zhǎng)期服役在惡劣的海洋環(huán)境中，并受到各種載荷的交互作用，如風(fēng)載荷、海流、波浪載荷、冰載荷等，有時(shí)還要遭到地震、臺(tái)風(fēng)、海嘯、船碰撞等意外打擊，結(jié)構(gòu)本身還要遭受環(huán)境腐蝕、海洋生物附著、海底沖刷等影響的作用。在這些惡劣的環(huán)境載荷長(zhǎng)期作用下，再加上設(shè)計(jì)或使用的不當(dāng)，結(jié)構(gòu)容易產(chǎn)生各種形式的損傷，使結(jié)構(gòu)的承載能力下降，嚴(yán)重的還會(huì)導(dǎo)致平臺(tái)失效嘲。在國(guó)內(nèi)外海洋開發(fā)工程中，曾發(fā)生過多起災(zāi)難性海洋平臺(tái)事故，造成了巨大的人員傷亡、經(jīng)濟(jì)損

10、失以及不良的社會(huì)影響。 隨著石油開采向海洋發(fā)展，海洋平臺(tái)的數(shù)量成倍增加，合適的設(shè)計(jì)方法確保結(jié)構(gòu)能夠抵抗住不可預(yù)測(cè)的載荷造成的損傷，但是損傷在海洋平臺(tái)結(jié)構(gòu)的服役期間是不可避免的，確保人的生命安全和減少財(cái)產(chǎn)損失的唯一方法是診斷出結(jié)構(gòu)的損傷，并能及時(shí)進(jìn)行修復(fù)。由此可見，提高海洋平臺(tái)結(jié)構(gòu)及設(shè)備的可靠性，確保海洋作業(yè)安全的問題日益突出，新平臺(tái)的質(zhì)量評(píng)價(jià)、舊平臺(tái)的殘余壽命估計(jì)和在役平臺(tái)的結(jié)構(gòu)安全保證將成為日益突出的問題，海洋平臺(tái)結(jié)構(gòu)的健康監(jiān)測(cè)與損傷診斷已成為刻不容緩的重要課題。 當(dāng)前的結(jié)構(gòu)損傷檢測(cè)的方法很多，除了人工目測(cè)外，還有超聲波、磁場(chǎng)法、放射法、熱力場(chǎng)等局部檢測(cè)方法。然而，較弱的視覺觀測(cè)條件以及損傷

11、部位有可能被生長(zhǎng)的海洋生物覆蓋著，所以利用這類局部損傷檢測(cè)技術(shù)方法對(duì)海洋平臺(tái)結(jié)構(gòu)進(jìn)行損傷診斷是不可靠的。此外，這些技術(shù)要求結(jié)構(gòu)的損傷區(qū)域是己知作為先決條件，要求配備特殊額外的測(cè)試設(shè)備和專業(yè)人員，因此，這些方法的檢測(cè)成本較昂貴。與上述方法相比較，基于振動(dòng)測(cè)試的結(jié)構(gòu)健康監(jiān)測(cè)技術(shù)是相對(duì)簡(jiǎn)單、成本較低的，被公認(rèn)為是較有發(fā)展前景的全局性損傷診斷方法。這種方法的基本原理是：損傷將導(dǎo)致結(jié)構(gòu)的系統(tǒng)剛度和阻尼矩陣發(fā)生改變，因而導(dǎo)致結(jié)構(gòu)的動(dòng)力特性參數(shù)(如結(jié)構(gòu)的頻響函數(shù)，模態(tài)參數(shù)等)的變化。換言之，結(jié)構(gòu)動(dòng)力特性參數(shù)能夠作為結(jié)構(gòu)損傷診斷的指標(biāo)。這類方法最突出的優(yōu)點(diǎn)是整個(gè)損傷診斷操作過程不會(huì)影響結(jié)構(gòu)的正常工作。 在土木

12、工程領(lǐng)域，海洋平臺(tái)結(jié)構(gòu)、橋梁和大壩等工程結(jié)構(gòu)被視為“系統(tǒng)”，而“識(shí)別”則意味著從振動(dòng)測(cè)試數(shù)據(jù)中識(shí)別出結(jié)構(gòu)的動(dòng)力特性參數(shù)(模態(tài)參數(shù))。對(duì)結(jié)構(gòu)物而言，模態(tài)參數(shù)是結(jié)構(gòu)的“指紋”，它是一系列獨(dú)特的數(shù)據(jù)，能夠反映結(jié)構(gòu)本身的固有動(dòng)力特性。每一個(gè)結(jié)構(gòu)都有其固有的模態(tài)參數(shù)，并且如果結(jié)構(gòu)動(dòng)力特性發(fā)生變化了，那么結(jié)構(gòu)的“指紋”也將發(fā)生相應(yīng)的變化。因此，模態(tài)分析是結(jié)構(gòu)動(dòng)態(tài)設(shè)計(jì)以及設(shè)備故障診斷的重要方法。基于輸出的模態(tài)參數(shù)識(shí)別方法利用的信息主要是系統(tǒng)的自由振動(dòng)信號(hào)，要獲得自由振動(dòng)信號(hào)首先需獲得結(jié)構(gòu)的響應(yīng)信號(hào)。由于環(huán)境激勵(lì)的不充分和噪聲等干擾因素的存在，響應(yīng)信號(hào)中包含的有用信息十分微弱，尤其對(duì)海洋平臺(tái)結(jié)構(gòu)而言，其服役環(huán)

13、境惡劣而復(fù)雜，現(xiàn)場(chǎng)測(cè)試的信號(hào)中包含了較多的噪聲成分，因此不能直接用于參數(shù)辨識(shí)，這是導(dǎo)致參數(shù)識(shí)別方法難以奏效的主要原因。因此，需要對(duì)信號(hào)進(jìn)行消噪處理，從大量背景噪聲中提取出可用于模態(tài)參數(shù)辨識(shí)的有用信號(hào)成分,剔除響應(yīng)信號(hào)中的干擾因素，提取有用信息成為關(guān)鍵。此時(shí)，信號(hào)消噪技術(shù)研究變得尤為重要。 目前，有大量文獻(xiàn)對(duì)信號(hào)降噪技術(shù)進(jìn)行了研究，提出了多種降噪方法，如時(shí)域平均法1、小波降噪技術(shù)2、頻域特征抽取技術(shù)3、自適應(yīng)濾波技術(shù)4等。然而，各種方法在實(shí)際應(yīng)用中都有各自的局限性，時(shí)域平均法需要有足夠的數(shù)據(jù)量，并且在使用過程中必須有時(shí)標(biāo)信息的支持；小波降噪和自適應(yīng)濾波技術(shù)很大程度上依賴于濾波器性能；頻域特征抽取

14、技術(shù)過度依賴于信號(hào)的幅值、頻率、相位信息，計(jì)算起來很不方便，而且對(duì)多譜勒等變頻信號(hào)而言，無法成功降噪。 近年也發(fā)展了一些其它方法，其中基于奇異值分解(SVD)的降噪技術(shù)因其計(jì)算方法簡(jiǎn)單易用引起了國(guó)內(nèi)外相當(dāng)一部分專家學(xué)者的重視。奇異值分解技術(shù)在聲學(xué)、智能控制、電子學(xué)、信號(hào)處理等領(lǐng)域得到了廣泛的應(yīng)用。Fort等人5利用矩陣SVD方法和一個(gè)新的標(biāo)準(zhǔn)(dynamic mean evaluation，DME)確定模型階次，再利用AR模型進(jìn)行譜估計(jì)，將其應(yīng)用到benchmark和Doppler信號(hào)分析中；Sanliturk等人6將Hankel矩陣和SVD算法結(jié)合，從復(fù)雜的噪聲信號(hào)中獲得較高精確度的頻響函數(shù)

15、；Vrabie等人7在SVD算法的基礎(chǔ)上引入獨(dú)立分量分析的概念，可以在傳感器相互干擾很大的情況下，在低通子空間很好的分離原波形，并應(yīng)用于垂直地震剖面的分析；馬寨璞8等人在卡爾曼濾波的基礎(chǔ)上，提出了利用矩陣的奇異值將數(shù)據(jù)矩陣進(jìn)行SVD分解的新的簡(jiǎn)化方法，得到能夠描述原狀態(tài)向量的新的較少維數(shù)向量的有效秩，并和HAMSON模式結(jié)合，利用渤海區(qū)域的SST進(jìn)行實(shí)驗(yàn)得到驗(yàn)證。Eckart和Young9提出了截?cái)嗥娈愔捣纸饧夹g(shù)，即只保留系統(tǒng)階數(shù)范圍內(nèi)的幾個(gè)最大奇異值,范圍外的奇異值均設(shè)為零。該理論成為低秩逼近問題的核心,為后來的降噪技術(shù)發(fā)展提供了新的思路。近年來,低秩逼近問題己發(fā)展到結(jié)構(gòu)低秩逼近,矩陣從一般

16、形式發(fā)展到特殊形式，如Hankel、Toeplitz矩陣等。已有很多學(xué)者提出結(jié)構(gòu)完全最小二乘法。Aoki和Yue最早做了對(duì)結(jié)構(gòu)完全最小二乘法問題的研究工作，在多年后的文獻(xiàn)中才首次出現(xiàn)但這一問題的提法。Cadzowt,Bresler和Macovski提出了另一種求解方法，最終這些方法被證明只滿足次優(yōu)的L2優(yōu)化準(zhǔn)則,盡管如此,它們由于簡(jiǎn)單易用而被廣泛采用。基于SVD分解的消噪技術(shù),它是從矩陣的角度出發(fā),將包含信號(hào)特征的矩陣分解到一系列奇異值和奇異值矢量對(duì)應(yīng)的子空間中，近些年的研究對(duì)SVD算法進(jìn)行了有益的應(yīng)用和改進(jìn)。該方法有兩個(gè)關(guān)鍵: 1：如何確定分解后重構(gòu)的有效秩階次; 2：如何確定重構(gòu)矩陣的行列

17、數(shù)。 針對(duì)有效秩階次的選擇， 常用的方法是試湊法和閾值法，均依賴于經(jīng)驗(yàn)， 缺乏依據(jù)。目前,對(duì)此問題,已經(jīng)出現(xiàn)了一些研究方法。如穩(wěn)定圖法，該方法通過在頻譜圖上標(biāo)示出滿足一定條件的穩(wěn)定極點(diǎn)，并被認(rèn)為是系統(tǒng)的真實(shí)極點(diǎn)。但這種方法不能完全排除噪聲模態(tài)，特別是隨著模型階數(shù)的升高，一些擬合模型的擬合模態(tài)往往容易趨于穩(wěn)定，用穩(wěn)定圖很難完全并正確確定模型階次。有人提出了基于奇異值分解的模型定階與降噪技術(shù)，該類方法一般將奇異值由大到小按降序排列，并將奇異值以最大值歸一化，通過畫出奇異值歸一化曲線，在曲線上找到突降的位置，該處對(duì)應(yīng)的奇異值個(gè)數(shù)即為模型階次，也即信號(hào)中包含模態(tài)數(shù)目的兩倍。對(duì)受噪聲影響的數(shù)據(jù)，奇異值曲

18、線突降不明顯，而是趨向于一條水平漸近線。一般認(rèn)為奇異值曲線開始變?yōu)樗降狞c(diǎn)對(duì)應(yīng)模型的階次。此外還有朱啟兵10提出了基于結(jié)構(gòu)風(fēng)險(xiǎn)最小化原則的奇異值分解降噪方法，該方法依據(jù)統(tǒng)計(jì)學(xué)習(xí)理論， 把有效秩階次的選擇看作是一個(gè)學(xué)習(xí)過程，利用結(jié)構(gòu)風(fēng)險(xiǎn)最小化原則來代替?zhèn)鹘y(tǒng)的經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化，從而自動(dòng)得到奇異值分解降噪中矩陣的有效秩。王維11提出了基于非監(jiān)督動(dòng)態(tài)聚類算法來確定矩陣有效重構(gòu)階次的方法，該方法利用含噪聲信號(hào)的奇異譜圖中表征噪聲的噪聲平臺(tái)平緩和集中的特性，通過向譜圖縱軸投影，應(yīng)用動(dòng)態(tài)聚類合理確定噪聲平臺(tái)的邊界，進(jìn)而有效地確定奇異值分解降噪中矩陣的有效重構(gòu)階次。康春玉12采用主分量分析的方法，根據(jù)奇異值的大

19、小來確定有效秩的階次。孫鑫暉13提出了通過奇異熵增量確定降噪階次的方法，信號(hào)的奇異熵是信息熵的一種改進(jìn)形式，能夠反映信號(hào)包含信息量的多少。當(dāng)信號(hào)受到寬頻帶噪聲干擾，信號(hào)的奇異熵隨著階次升高一直增大。在較低階次時(shí)，在信號(hào)與噪聲的共同作用下，奇異熵增長(zhǎng)速度較快。當(dāng)達(dá)到一定階次后，奇異熵增長(zhǎng)速度放緩。這時(shí)信號(hào)的有效特征信息量已經(jīng)趨于飽和，之后的增量是由于噪聲所致。因此,選擇奇異熵達(dá)到飽和階次作為重構(gòu)階次，能夠保留信號(hào)信息同時(shí)去除噪聲。 針對(duì)重構(gòu)矩陣的行列數(shù)的選擇問題，普遍采取的方法是根據(jù)具體信號(hào)選擇不同行列數(shù)進(jìn)行試湊，這種試湊法需要大量的計(jì)算，并且嚴(yán)重依賴使用者的信號(hào)分析經(jīng)驗(yàn)。Kanjilal14提

20、出了通過奇異值比譜來確定行列數(shù)的方法，但是該方法在信號(hào)由多個(gè)周期分量組成或噪聲較嚴(yán)重的情況下效果并不明顯。趙學(xué)智15提出通過分析所有分解分量信號(hào)所含信息量的變化趨勢(shì)來確定合理的矩陣結(jié)構(gòu)。上述方法在實(shí)際應(yīng)用中取得了較好的效果，但也存在著一定的局限性。本文將Hankel矩陣與SVD分解相結(jié)合，首先對(duì)測(cè)量信號(hào)構(gòu)造的Hankel矩陣進(jìn)行奇異值分解，得到測(cè)量信號(hào)的奇異值（包括有用信號(hào)和噪聲信號(hào)的奇異值），再利用測(cè)量信號(hào)快速傅立葉變換結(jié)果中主頻率的個(gè)數(shù)來確定有效秩階次，剔除噪聲信號(hào)的奇異值，接著通過消噪信號(hào)的信噪比和均方差大小確定重構(gòu)矩陣結(jié)構(gòu)，得到有用信號(hào)的矩陣，最后通過反對(duì)角線平均法得到消噪后的信號(hào)數(shù)據(jù)

21、。2 SVD分解消噪理論2.1 Hankel矩陣 在信息與信號(hào)處理領(lǐng)域，要想將有用信號(hào)不失真地變換和處理幾乎是不可能的。誤差的來源有：信息傳輸處理時(shí)，信道或設(shè)備不理想造成的；傳輸處理過程中串入的一些其它信號(hào)也能引起誤差。噪聲是使信號(hào)產(chǎn)生失真的誤差源。來自外部的噪聲也稱為干擾。從功率譜的角度看，如果一個(gè)隨機(jī)過程的功率譜密度是常數(shù)，無論是什么分布，都稱為白噪聲。白噪聲的頻率成分非常豐富。假設(shè)有一個(gè)載有信息的、以時(shí)間t為變量的信號(hào)，由于在傳輸過程中受到了某種加性噪聲或干擾的污染，致使觀察到的信號(hào)發(fā)生變化，其可以表示為： 對(duì)于一個(gè)測(cè)得的信號(hào)，其中為有用信號(hào)，為噪聲信號(hào)，基于相空間重構(gòu)理論16，可以由其

22、構(gòu)造Hankel矩陣 消噪的目的就是從A矩陣中求出的最佳逼近矩陣，再求出。 2.2 SVD分解的基本理論 設(shè)是一個(gè)秩為r的維矩陣，則的奇異值分解是指，存在矩陣和以及，使得： 其中，U，V分別為，維正交矩陣，0為零元素矩陣。 為維對(duì)角陣，其對(duì)角線元素為矩陣H的非零奇異值，且以非增順序排列，即 ，且有。 矩陣的秩為r，從式（2）中除去H的零奇異值，得到奇異值分解的精簡(jiǎn)形式 其中：，分別為U，V的第個(gè)行向量。2.3 對(duì)測(cè)量信號(hào)進(jìn)行SVD分解 根據(jù)式（2）對(duì)構(gòu)造的Hankel矩陣（即式（1）進(jìn)行SVD分解： 都為維矩陣，為維正交矩陣，為維正交矩陣，為維對(duì)角矩陣，其對(duì)角元素分別為的奇異值，同時(shí)還可看出。

23、上述奇異值分解式根據(jù)式（4）還可以進(jìn)一步寫成 由上式可以看出，如果信號(hào)中沒有噪聲或信噪比特別高，則矩陣是奇異的， 即的數(shù)目，的大小和系統(tǒng)有關(guān)。如果有噪聲或信噪比不高，則矩陣是非奇異的，即的數(shù)目。因源信號(hào)是由有用信號(hào)和噪聲信號(hào)共同組成，則矩陣也是由有用信號(hào)和噪聲信號(hào)共同組成的矩陣，那么矩陣的奇異值可以反映信號(hào)和噪聲信號(hào)能量集中的情況。前r個(gè)較大的奇異值將主要反映有用信號(hào)，較小的奇異值則主要反映噪聲信號(hào)，把這部分反映噪聲的奇異值置零就可以去除信號(hào)中的噪聲。再利用奇異值分解的逆過程得到矩陣，即矩陣的秩為r的最佳逼近矩陣。相對(duì)于其噪聲已被大大壓縮。將中對(duì)應(yīng)的元素相加平均，就可以得到降噪后的信號(hào) 。3

24、有效秩階次和重構(gòu)矩陣結(jié)構(gòu)的確定3.1 有效秩階次的確定 由實(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)行逼近、剔除毛刺的過程。經(jīng)過對(duì)大量仿真結(jié)果的研究，發(fā)現(xiàn)去噪結(jié)果中剩余噪聲對(duì)信號(hào)的干擾影響能夠通過信號(hào)中“毛刺”的數(shù)量及大小來判斷，這與實(shí)測(cè)經(jīng)驗(yàn)相吻合。當(dāng)選取的奇異值數(shù)目較小時(shí)，大部分噪聲都被剔除掉了，但也損失了大量有用信號(hào)。隨著奇異值數(shù)目的逐漸增大，有用信號(hào)信息趨于完整，但噪聲的干擾也會(huì)逐漸增加。當(dāng)奇異值數(shù)目達(dá)到某一值時(shí)，降噪后的

25、信號(hào)既保留了大部分有用信息，也剔除掉了大部分噪聲，這就是要選取的最佳奇異值數(shù)目。由第個(gè)非零奇異值重構(gòu)得到重構(gòu)信號(hào)分量， 分別對(duì)這些分量信號(hào)進(jìn)行傅里葉變換后發(fā)現(xiàn)，分量的頻率成分均為源信號(hào)的頻率成分組成，而由較大的奇異值重構(gòu)得到的分量信號(hào)其頻率成分與源信號(hào)中主頻率相對(duì)應(yīng)。顯然， 當(dāng)有用信號(hào)未被噪聲完全淹沒時(shí)，源信號(hào)中主頻率是有用信號(hào)的頻率。因此可以通過奇異值與有用信號(hào)頻率之間的某種對(duì)應(yīng)關(guān)系來確定有效秩的階次，從而達(dá)到降噪的效果。 給定一個(gè)源信號(hào)，其中：為有用信號(hào)，為強(qiáng)度為1的高斯白噪聲。將和快速傅立葉變換FFT的結(jié)果如圖1和圖2所示，比較兩張圖，發(fā)現(xiàn)它們都包含兩個(gè)主頻率成分，圖2不含其他頻率，由此

26、可知這兩個(gè)主頻率為有用信號(hào)的頻率成分，其他均為噪聲的頻率成分。圖1 源信號(hào)的時(shí)域圖和功率譜密度圖圖2 有用信號(hào)的時(shí)域圖和功率譜密度圖分別對(duì)源信號(hào)和有用信號(hào)利用式(1)構(gòu)造相應(yīng)的Hankel矩陣，并對(duì)其進(jìn)行奇異值分解。隨著矩陣的行數(shù)變化，其奇異值的個(gè)數(shù)分布如圖3和圖4所示：圖3 源信號(hào)奇異值個(gè)數(shù)隨矩陣行數(shù)變化情況圖4 有用信號(hào)奇異值個(gè)數(shù)隨矩陣行數(shù)變化情況 可以看到：經(jīng)過對(duì)無噪聲信號(hào)和受噪聲污染信號(hào)的研究，發(fā)現(xiàn)一般由無噪聲理想信號(hào)構(gòu)造的Hankel矩陣的大部分奇異值為零。根據(jù)非零奇異值的數(shù)目，可以很容易地判斷矩陣的有效階次。因噪聲具有隨機(jī)和不相關(guān)的特點(diǎn)，因此，由受隨機(jī)噪聲污染信號(hào)構(gòu)造的Hankel

27、 矩陣呈列滿秩或行滿秩狀態(tài)（取決于行和列哪個(gè)維數(shù)更小）。具體分析如下：1：當(dāng)行數(shù)或列數(shù)小于主頻個(gè)數(shù)的兩倍時(shí)，無論是有用信號(hào)還是源信號(hào)，其非零奇異值的個(gè)數(shù)等于矩陣行數(shù)和列數(shù)之中的最小值。2：當(dāng)行數(shù)大于主頻個(gè)數(shù)的兩倍時(shí)，有用信號(hào)的非零奇異值一直是4個(gè)，且不隨矩陣行數(shù)的變化而變化；當(dāng)行數(shù)小于列數(shù)時(shí)，源信號(hào)的非零奇異值等于行數(shù)的值，當(dāng)行數(shù)大于列數(shù)時(shí)，源信號(hào)的非零奇異值等于列數(shù)的值。由此可知，源信號(hào)經(jīng)過奇異值分解等到的對(duì)角矩陣的對(duì)角元素都不為零。 再對(duì)源信號(hào)和有用信號(hào)的矩陣經(jīng)過奇異值分解之后得到的對(duì)角矩陣的對(duì)角元素大小進(jìn)行分析，如圖5和圖6所示：圖5 源信號(hào)奇異值大小隨矩陣行數(shù)變化情況因上圖點(diǎn)數(shù)較多，不

28、能確定大奇異值個(gè)數(shù)，隨選擇較短行數(shù)以便觀察。圖51 矩陣行數(shù)在1-20之間時(shí)，奇異值大小情況圖52 矩陣行數(shù)在70-120之間時(shí)，奇異值大小情況圖6 有用信號(hào)奇異值大小隨矩陣行數(shù)變化情況對(duì)比圖5和圖6中的奇異值曲線可知，每個(gè)奇異值受到噪聲干擾的程度不同：較大的奇異值受噪聲干擾的影響較小，而較小的奇異值受噪聲干擾的影響則較大。由于由受噪聲污染信號(hào)構(gòu)造的Hankel 矩陣的每個(gè)奇異值都是由“有用信號(hào)”和“噪聲信號(hào)”兩部分組成，實(shí)際選取奇異值數(shù)目就是一個(gè)對(duì)有用信號(hào)和噪聲信號(hào)進(jìn)行取舍的過程。從圖6中可以看出，隨著行數(shù)的增加，奇異值的個(gè)數(shù)恒為4，是源信號(hào)快速傅立葉變換中主頻個(gè)數(shù)的2倍；而圖5中較大奇異值

29、的個(gè)數(shù)也一直為4，說明其主要反映有用信號(hào)的信息，將這些奇異值稱為大奇異值，而圖5中其他奇異值相對(duì)較小且分布比較集中，說明其反映出了噪聲的特點(diǎn)。所以，在消噪過程中，選取的奇異值個(gè)數(shù)越接近有用信號(hào)的奇異值個(gè)數(shù)，就越是接近不受噪聲污染的原始信號(hào)，當(dāng)奇異值數(shù)目過小時(shí)，就會(huì)漏掉有用信號(hào)，奇異值過大時(shí)，降噪后的信號(hào)中就會(huì)含有大量的噪聲。由此可知：有效秩的階次與源信號(hào)的主頻個(gè)數(shù)存在一個(gè)確定的關(guān)系， 即有效秩的階次為源信號(hào)快速傅里葉變換后主頻的個(gè)數(shù)的兩倍。利用這種倍數(shù)關(guān)系可以確定有效秩的階次，取得最佳的奇異值數(shù)目，不僅可以過濾掉大量的噪聲，而且還能最大程度地保留原始信號(hào)的信息。3.2 重構(gòu)矩陣結(jié)構(gòu)的確定確定有

30、效秩的階次之后，就可以提取出信號(hào)對(duì)角矩陣中的有用奇異值，即剔除了測(cè)試信號(hào)中的噪聲奇異值，接著要重構(gòu)有用信號(hào)的矩陣，這需要確定重構(gòu)矩陣的結(jié)構(gòu)，即矩陣的行列數(shù)，因?yàn)椴煌辛袛?shù)會(huì)導(dǎo)致不同的消噪效果。本文根據(jù)降噪效果的好壞來確定最佳的矩陣結(jié)構(gòu)。 降噪效果一般用信號(hào)的均方誤差( mean square error，簡(jiǎn)稱MSE)和信噪比(signal to noise ratio，簡(jiǎn)稱SNR) 來衡量: MSE 越小，SNR 越大，降噪效果越好。MSE 與SNR 的定義形式如下 其中：為含噪聲信號(hào)的第k個(gè)數(shù)據(jù)點(diǎn)；為無噪聲信號(hào)的第k個(gè)數(shù)據(jù)點(diǎn)；N為信號(hào)長(zhǎng)度。 選取不同行數(shù)L(L>5)對(duì)信號(hào)(N = 20

31、0)，噪聲強(qiáng)度為1dB，進(jìn)行奇異值分解與重構(gòu)，降噪信號(hào)的SNR與MSE隨L變化如圖7所示。可以看出，當(dāng)L增加到一定程度時(shí)，降噪信號(hào)的SNR基本穩(wěn)定在22dB左右，僅有小的波動(dòng)，此時(shí)MSE的變化也呈現(xiàn)出這種規(guī)律。由結(jié)果可知，當(dāng)L=102時(shí)，降噪信號(hào)的MSE取得最小值，SNR取得最大值，此時(shí)降噪效果最好。圖7 均方誤差和信噪比隨行數(shù)變化 進(jìn)行分析后發(fā)現(xiàn)，SNR與MSE分別取最大和最小值時(shí)的L值不一定相等， 但最佳L值基本出現(xiàn)在處的一個(gè)領(lǐng)域內(nèi)，并且L在該鄰域內(nèi)取值時(shí)，降噪效果較好且差異較小，均能滿足要求。因此，重構(gòu)矩陣的結(jié)構(gòu)可以根據(jù)N來確定，實(shí)驗(yàn)應(yīng)用中不妨取( 當(dāng)N不是偶數(shù)時(shí)，舍棄最后一個(gè)數(shù)據(jù)點(diǎn)，不

32、影響最終結(jié)果)。4 消噪后的信號(hào)重構(gòu) 基于SVD分解的消噪技術(shù),它是從矩陣的角度出發(fā),將包含信號(hào)特征的矩陣分解到一系列奇異值和奇異值矢量對(duì)應(yīng)的子空間中，通快速過傅立葉變換結(jié)果中主頻率個(gè)數(shù)確定有用信號(hào)的奇異值個(gè)數(shù)以保留矩陣的有用信號(hào)奇異值，并將其余奇異值置零以剔除信號(hào)噪聲的奇異值。再根據(jù)有用信號(hào)奇異值重構(gòu)有用信號(hào)的矩陣，因?yàn)椴煌辛袛?shù)會(huì)導(dǎo)致不同的消噪效果，根據(jù)均方誤差和信噪比可以知道當(dāng)行數(shù)為信號(hào)數(shù)據(jù)長(zhǎng)度的一半時(shí)（即矩陣為方陣），降噪效果最好。奇異值個(gè)數(shù)和重構(gòu)矩陣行數(shù)確定后，就可以依據(jù)此得到消噪后的信號(hào)矩陣，但矩陣并不等于由真實(shí)信號(hào)構(gòu)成的Hankel矩陣A，不是嚴(yán)格的Hankel矩陣。通過觀察式（

33、1），可以看出真實(shí)信號(hào)與矩陣A的各元素之間存在如下關(guān)系： 其中，。也就是對(duì)矩陣A的反對(duì)角線求平均值即可得到真實(shí)信號(hào)在每一時(shí)刻的值，。根據(jù)這一思路，同樣對(duì)矩陣的反對(duì)角線求平均值，從而得到經(jīng)過消噪后，信號(hào)在每一時(shí)刻的估計(jì)值，即： 其中，。5 數(shù)值仿真在信息與信號(hào)處理領(lǐng)域，將有用信號(hào)不失真地變換和處理是不可能的，因?yàn)樵谛畔鬏斕幚頃r(shí)，信道或設(shè)備的不理想會(huì)造成誤差，或者在傳輸處理過程中會(huì)串入一些其它信號(hào)即噪聲。本文將Hankel矩陣與SVD分解相結(jié)合，首先對(duì)測(cè)量信號(hào)構(gòu)造的Hankel矩陣進(jìn)行奇異值分解，因源信號(hào)是由有用信號(hào)和噪聲信號(hào)共同組成，則矩陣也是由有用信號(hào)和噪聲信號(hào)共同組成的矩陣，那么矩陣的奇異

34、值可以反映信號(hào)和噪聲信號(hào)能量集中的情況。前r個(gè)較大的奇異值將主要反映有用信號(hào)，較小的奇異值則主要反映噪聲信號(hào)，把這部分反映噪聲的奇異值置零就可以去除信號(hào)中的噪聲。本文利用測(cè)量信號(hào)快速傅立葉變換結(jié)果中主頻率的個(gè)數(shù)來確定有效秩階次，剔除噪聲信號(hào)的奇異值。因?yàn)橹貥?gòu)矩陣行數(shù)的不同會(huì)影響降噪效果，通過消噪信號(hào)的信噪比和均方差大小比較可發(fā)現(xiàn)，SNR與MSE分別取最大和最小值時(shí)的矩陣行數(shù)L值不一定相等，但最佳L值基本出現(xiàn)在處的一個(gè)領(lǐng)域內(nèi)，N為信號(hào)數(shù)據(jù)長(zhǎng)度，并且L在該鄰域內(nèi)取值時(shí)，降噪效果較好且差異較小。最后根據(jù)重構(gòu)矩陣和信號(hào)方程之間的關(guān)系，通過反對(duì)角線平均法得到消噪后的信號(hào)數(shù)據(jù)。5.1 Matlab仿真結(jié)果

35、分析 基于上述分析，對(duì)于一個(gè)含噪聲的測(cè)試信號(hào)，其降噪的基本步驟如下:(1) 取信號(hào)數(shù)據(jù)長(zhǎng)度的一半作為重構(gòu)矩陣的行數(shù)，根據(jù)式(1)構(gòu)造Hankel矩陣并進(jìn)行奇異值分解;(2) 對(duì)信號(hào)進(jìn)行快速傅里葉變換，確定主頻個(gè)數(shù)n，以2n作為有效秩的階次;(3) 用前2n個(gè)奇異值根據(jù)式(4)進(jìn)行重構(gòu)，得到重構(gòu)矩陣，將中對(duì)應(yīng)的元素根據(jù)式（10）相加后平均就可得到降噪后的信號(hào)。 分別用不同頻率成分的信號(hào)對(duì)該方法進(jìn)行驗(yàn)證：信號(hào)1：，高斯白噪聲強(qiáng)度為1dB，10dB，15dB；信號(hào)2：,高斯白噪聲強(qiáng)度為10dB；信號(hào)3：，高斯白噪聲強(qiáng)度為10dB；信號(hào)4：，高斯白噪聲強(qiáng)度為10dB。信號(hào)1，2，3取數(shù)據(jù)長(zhǎng)度為200，

36、可將其構(gòu)成矩陣的行數(shù)設(shè)為100；信號(hào)4取數(shù)據(jù)長(zhǎng)度為400，可將其構(gòu)成矩陣的行數(shù)設(shè)為200。表1 信號(hào)1消噪前后SNR和MSE的變化 為了考察這種降噪方法在不同噪聲水平下的表現(xiàn),分別進(jìn)行了對(duì)含1dB, 10dB, 15dB高斯白噪聲信號(hào)1的消噪處理, 發(fā)現(xiàn)在不同信噪比下, 通過這種方法得到的消噪波形都能較好地保留目標(biāo)信號(hào)的波形特征,分別對(duì)各個(gè)信號(hào)進(jìn)行噪聲消除，得到結(jié)果如下：信號(hào)1噪聲強(qiáng)度/dB源信號(hào)信噪比SNR消噪后信號(hào)信噪比SNR源信號(hào)均方誤差MSE消噪后信號(hào)均方誤差MSE136.63237.4870.0558960.0490791011.43114.4680.866270.65448157.

37、88239.69091.8171.6677圖8 信號(hào)1-噪聲1dB圖9 信號(hào)1-噪聲10dB圖10 信號(hào)1-15dB 從表1的可以看出，降噪后信號(hào)的信噪比SNR都提高了，均方誤差MSE都明顯降低了，這說明本文給出的降噪方法是有效的。再看圖8到圖10中的波形，可以看出，隨著噪聲強(qiáng)度的增加，源信號(hào)波形失真越嚴(yán)重，功率譜密度圖中噪聲頻率的個(gè)數(shù)也逐漸增多，并且逐漸淹沒主頻率。當(dāng)噪聲強(qiáng)度不至于淹沒有用信號(hào)時(shí)，降噪后的波形與原波形吻合較好，證明了這種方法的降噪效果。圖11 信號(hào)2-10dB圖12 信號(hào)3-10dB圖13 信號(hào)4-10dB 從傅立葉變換結(jié)果可以看出信號(hào)2，3，4的主頻個(gè)數(shù)分別為2，3，4，那

38、么可以確定重構(gòu)矩陣的有效秩階次分別為4，6，8，由此可以剔除信號(hào)中的噪聲信號(hào)。因?yàn)樾盘?hào)2和3的數(shù)據(jù)長(zhǎng)度為200，信號(hào)4的數(shù)據(jù)長(zhǎng)度為400，根據(jù)上文分析得出的當(dāng)重構(gòu)矩陣的行數(shù)為數(shù)據(jù)長(zhǎng)度的一半時(shí)信噪比最大的理論，信號(hào)2和3重構(gòu)矩陣行數(shù)為100，信號(hào)4的重構(gòu)矩陣的行數(shù)為200。因?yàn)楫?dāng)信號(hào)受到噪聲干擾時(shí)，測(cè)量得到的信號(hào)中含有大量的“毛刺”，利用噪聲污染信號(hào)構(gòu)造Hankel 矩陣進(jìn)行奇異值分解降噪，就是對(duì)含噪信號(hào)進(jìn)行逼近、剔除毛刺的過程，所以可以根據(jù)去噪結(jié)果中“毛刺”的數(shù)量及大小來判斷消噪效果。觀察各圖波形，發(fā)現(xiàn)消噪后信號(hào)波形基本都是光滑曲線，幾乎沒有“毛刺”，這表明大部分噪聲都被剔除掉了。并且發(fā)現(xiàn)消噪

39、后信號(hào)波形與有用信號(hào)波形吻合較好，證明了這種消噪方法的實(shí)用性。實(shí)驗(yàn)結(jié)果表明, 利用源信號(hào)主頻個(gè)數(shù)來確定有效秩的階次以及取信號(hào)數(shù)據(jù)長(zhǎng)度的一半確定重構(gòu)矩陣的行數(shù)的方法可以得到較好的降噪效果，該方法可以獲得較高的信噪比,同時(shí)也較好的保留了原信號(hào)的特征波形，證明該方法是十分有效的。5.2 Matlab程序clear;Fs=1;%采樣頻率n=400;%采樣個(gè)數(shù) L=n/2;%矩陣最佳行數(shù)，數(shù)據(jù)長(zhǎng)度一半 A=1:L;%第一列 B=L:n;%最后一行 c=hankel(A,B);%構(gòu)造hankel矩陣 h=sin(0.01*pi*c)+2*cos(0.03*pi*c)+cos(0.05*pi*c)-2*si

40、n(0.07*pi*c);%有用信號(hào) y=wgn(size(c,1),size(c,2),10);%加噪聲 h=y+h;%源信號(hào) for i=1:n p=max(1,i-L+1); q=min(n-L+1,i); he=0; for j=p:q he=he+(h(i-j+1,j); end ca=q-p+1; x1(i)=he/ca; end figure(2); subplot(311); i=1:n; plot(i,x1,'r');%輸出源信號(hào)波形 grid; title('消噪前'); axis(0,n,-6,6);window=boxcar(length

41、(x1);%矩形窗nfft=n;%采樣點(diǎn)數(shù)Pxx,f=periodogram(x1,window,nfft,Fs);%直接法求功率譜密度subplot(312);plot(f,10*log10(Pxx);title('功率譜密度圖');xlabel('頻率Hz');ylabel('功率譜密度');axis(0,0.5,0,30);%畫出源信號(hào)功率譜密度圖grid; U,S,V=svds(h,8);%奇異值分解，確定有效秩階次 h2=U*S*V'%重構(gòu)有用信號(hào)的矩陣 for i=1:n p=max(1,i-L+1); q=min(n-L+1

42、,i); he=0; for j=p:q he=he+(h2(i-j+1,j); end ca=q-p+1; x2(i)=he/ca; end subplot(313); i=1:n; plot(i,x2);%輸出消噪后信號(hào)波形 hold on; grid; x0=sin(0.01*pi*i)+2*cos(0.03*pi*i)+cos(0.05*pi*i)-2*sin(0.07*pi*i); plot(i,x0,'-r');%輸出有用信號(hào)波形 title('消噪后'); axis(0,n,-6,6); 6 結(jié)束語 本文將Hankel矩陣與奇異值分解法相結(jié)合，根據(jù)

43、測(cè)量信號(hào)中有用奇異值是不變的，將測(cè)量信號(hào)構(gòu)造的Hankel矩陣進(jìn)行奇異值分解，得到一個(gè)對(duì)角矩陣，對(duì)角矩陣中包含的元素即測(cè)量信號(hào)中的奇異值（包括有用信號(hào)和噪聲信號(hào)的奇異值）。因?yàn)闇y(cè)量信號(hào)快速傅立葉變換結(jié)果中主頻率的個(gè)數(shù)的2倍即有效秩階次，可以據(jù)此剔除噪聲信號(hào)的奇異值。接著通過消噪信號(hào)的信噪比和均方差大小確定重構(gòu)矩陣結(jié)構(gòu)，利用奇異值分解的逆過程即可得到有用信號(hào)的矩陣，最后通過反對(duì)角線平均法得到消噪后的信號(hào)數(shù)據(jù)。仿真結(jié)果表明，對(duì)于不同頻率，不同噪聲強(qiáng)度的信號(hào)，該方法可以剔除大部分噪聲，獲得較高的信噪比,同時(shí)也能較好的保留原信號(hào)的特征波形，證明該方法是十分有效的。參考文獻(xiàn)1G.J.JanacekPractical time seriesLondon：Arnold，20012YY Kim，JC Hong，NYLeeFrequency response function estimation via a robustwavelet de-noising methodJ Sound Vib2001，244(4)：635-6493LLPresti，GOlm