基于奇異值分解的信號消噪技術

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

24、有效秩階次和重構矩陣結構的確定3.1 有效秩階次的確定 由實測經(jīng)驗可知，當實際測量環(huán)境很好時，測量得到的信號一般為光滑曲線；而當受到外界的隨機噪聲干擾時，測量得到的信號中就會含有大量的“毛刺”。利用噪聲污染信號構造Hankel 矩陣進行奇異值分解降噪，就是對含噪信號進行逼近、剔除毛刺的過程。經(jīng)過對大量仿真結果的研究，發(fā)現(xiàn)去噪結果中剩余噪聲對信號的干擾影響能夠通過信號中“毛刺”的數(shù)量及大小來判斷，這與實測經(jīng)驗相吻合。當選取的奇異值數(shù)目較小時，大部分噪聲都被剔除掉了，但也損失了大量有用信號。隨著奇異值數(shù)目的逐漸增大，有用信號信息趨于完整，但噪聲的干擾也會逐漸增加。當奇異值數(shù)目達到某一值時，降噪后的

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

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

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

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

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

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

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

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

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

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

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

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

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

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

39、后信號波形與有用信號波形吻合較好，證明了這種消噪方法的實用性。實驗結果表明, 利用源信號主頻個數(shù)來確定有效秩的階次以及取信號數(shù)據(jù)長度的一半確定重構矩陣的行數(shù)的方法可以得到較好的降噪效果，該方法可以獲得較高的信噪比,同時也較好的保留了原信號的特征波形，證明該方法是十分有效的。5.2 Matlab程序clear;Fs=1;%采樣頻率n=400;%采樣個數(shù) L=n/2;%矩陣最佳行數(shù)，數(shù)據(jù)長度一半 A=1:L;%第一列 B=L:n;%最后一行 c=hankel(A,B);%構造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);%有用信號 y=wgn(size(c,1),size(c,2),10);%加噪聲 h=y+h;%源信號 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');%輸出源信號波形 grid; title('消噪前'); axis(0,n,-6,6);window=boxcar(length

41、(x1);%矩形窗nfft=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);%畫出源信號功率譜密度圖grid; U,S,V=svds(h,8);%奇異值分解，確定有效秩階次 h2=U*S*V'%重構有用信號的矩陣 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);%輸出消噪后信號波形 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');%輸出有用信號波形 title('消噪后'); axis(0,n,-6,6); 6 結束語 本文將Hankel矩陣與奇異值分解法相結合，根據(jù)

43、測量信號中有用奇異值是不變的，將測量信號構造的Hankel矩陣進行奇異值分解，得到一個對角矩陣，對角矩陣中包含的元素即測量信號中的奇異值（包括有用信號和噪聲信號的奇異值）。因為測量信號快速傅立葉變換結果中主頻率的個數(shù)的2倍即有效秩階次，可以據(jù)此剔除噪聲信號的奇異值。接著通過消噪信號的信噪比和均方差大小確定重構矩陣結構，利用奇異值分解的逆過程即可得到有用信號的矩陣，最后通過反對角線平均法得到消噪后的信號數(shù)據(jù)。仿真結果表明，對于不同頻率，不同噪聲強度的信號，該方法可以剔除大部分噪聲，獲得較高的信噪比,同時也能較好的保留原信號的特征波形，證明該方法是十分有效的。參考文獻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