鄭 堂 李世平 羅 鵬 鄔肖敏

(第二炮兵工程大學(xué)，西安 710025)

0 引言

在工程實(shí)際中，測(cè)試得到的信號(hào)往往含有大量的噪聲，這嚴(yán)重影響了下一步對(duì)信號(hào)的處理和分析，同時(shí)，對(duì)于一個(gè)含有多個(gè)頻率成分的信號(hào)，需要知道其中所包含的具體頻率成分，就要對(duì)該信號(hào)進(jìn)行分解。

在去除噪聲方面，奇異值分解(SVD)方法已經(jīng)得到了成功應(yīng)用并被證明是有效的[1-2]，SVD方法通過(guò)將反映噪聲的奇異值置零能夠有效的達(dá)到去除噪聲的目的。經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)方法根據(jù)不同的時(shí)間尺度將信號(hào)分解為多個(gè)滿足一定條件的單一震蕩模式的線性疊加，每一個(gè)震蕩模式稱為一個(gè)本征模態(tài)函數(shù)(IMF)，通過(guò)這一方法，可以將信號(hào)中所包含的有用成分提取出來(lái)。

本文提出一種基于SVD與EMD相結(jié)合的方法，對(duì)于采集的信號(hào)，首先用SVD方法濾除噪聲，使信號(hào)的信噪比滿足EMD分解的要求，然后用EMD對(duì)信號(hào)進(jìn)行分解，得到其中所包含的各頻率成分，最終達(dá)到檢測(cè)多頻小信號(hào)的目的。

1 SVD濾除噪聲原理

1.1 SVD相關(guān)理論[3]

設(shè)有m×n階矩陣：

(1)

假定其中n≤m，那么A的秩r≤n。存在n階正交陣V和m階正交陣U，使得：

UTAV=Λ

(2)

其中，Λ為m×n階非負(fù)對(duì)角陣，即：

(3)

式中，σ1,σ2,…,σr和σr+1=σr+2=…=σn=0稱為A的奇異值；U、V的列向量分別為A的左右奇異向量。

對(duì)式(2)作等效變換可得：

A=UΛVT

(4)

上式就是矩陣的奇異值分解形式。

1.2 SVD去噪原理[4-6]

設(shè)含有噪聲的信號(hào)為：

f(n)=s(n)+u(n)

(5)

其中，s(n)為信號(hào)成分；u(n)為噪聲成分。設(shè)f(n)的信號(hào)長(zhǎng)度為N，那么由f(n)可構(gòu)造Hankel矩陣：

(6)

對(duì)F作SVD，它的奇異值σ1,σ2,…,σr可反映信號(hào)和噪聲能力集中的情況，將其按遞減順序排列，即σ(1)≥σ(2)≥…≥σ(r)，那么前i個(gè)較大的奇異值主要反映信號(hào)，后面較小的奇異值主要反映噪聲，將反映噪聲的奇異值置零即可達(dá)到去噪的目的。

至此，該方法有兩個(gè)關(guān)鍵問(wèn)題需要解決，其一是重構(gòu)Hankel矩陣結(jié)構(gòu)的確定，即矩陣行列數(shù)的確定；其二是有效秩個(gè)數(shù)的確定，即判別反映信號(hào)的奇異值個(gè)數(shù)的問(wèn)題。

由式(6)可看出，對(duì)于同一組信號(hào)序列，可構(gòu)造出多種結(jié)構(gòu)的Hankel矩陣，而不同的矩陣結(jié)構(gòu)對(duì)信號(hào)處理的效果影響是很大的[7]，針對(duì)這一問(wèn)題，可做一個(gè)仿真實(shí)驗(yàn)，通過(guò)對(duì)比不同結(jié)構(gòu)Hankel矩陣的降噪效果來(lái)確定最佳的矩陣結(jié)構(gòu)。

降噪效果一般用均方誤差來(lái)衡量，即：

(7)

圖1 MSE與重構(gòu)矩陣行數(shù)L關(guān)系圖

從圖中MSE與重構(gòu)矩陣行數(shù)之間存在的關(guān)系可以看出，當(dāng)重構(gòu)矩陣行數(shù)選在數(shù)據(jù)序列長(zhǎng)度一半，也就是50左右時(shí)，MSE值最小，此時(shí)SVD降噪效果最好，行數(shù)多于或少于50，MSE都會(huì)逐漸增大，降噪效果逐漸變差。所以，對(duì)于一組含有噪聲的信號(hào)序列，為了確保SVD降噪的效果，構(gòu)造的Hankel矩陣行數(shù)應(yīng)選為數(shù)據(jù)長(zhǎng)度的一半為宜，若信號(hào)長(zhǎng)度N為奇數(shù)，則行數(shù)可取為(N-1)/2。

對(duì)于有效秩個(gè)數(shù)的確定問(wèn)題，文獻(xiàn)[6]指出，有效秩個(gè)數(shù)與原信號(hào)中主頻率個(gè)數(shù)存在一定關(guān)系，即有效秩個(gè)數(shù)為主頻率個(gè)數(shù)兩倍時(shí)，降噪效果最好。但是，該方法存在缺陷，當(dāng)有用信號(hào)淹沒(méi)于噪聲中時(shí)將無(wú)法判別出原信號(hào)中所包含的主頻率個(gè)數(shù)，這種方法失效，針對(duì)這一問(wèn)題，本文提出了如下的解決方案。由于自相關(guān)不會(huì)改變信號(hào)的頻率組成，對(duì)于主頻率無(wú)法分辨的情況，可以先對(duì)信號(hào)做自相關(guān)處理，突出主頻率分量，從而得到原信號(hào)中主頻率的個(gè)數(shù)。

為了檢驗(yàn)該方法的效果，我們用一組含有兩個(gè)正弦成分和噪聲的信號(hào)進(jìn)行仿真實(shí)驗(yàn)，結(jié)果如圖2和圖3。通過(guò)圖2可以看出，已經(jīng)無(wú)法分辨出原信號(hào)中有幾個(gè)主頻率成分，而經(jīng)過(guò)自相關(guān)后，從圖3中可以清晰地分辨出原信號(hào)中有兩個(gè)主頻率成分，這與仿真信號(hào)所設(shè)定的完全一致，證明了該方法的有效性。通過(guò)這一方法，就能有效得出原信號(hào)中所含的主頻率數(shù)，從而確定SVD運(yùn)算時(shí)有效秩的個(gè)數(shù)。

圖2 自相關(guān)前信號(hào)的FFT圖

圖3 自相關(guān)后信號(hào)的FFT圖

綜合以上所講的，運(yùn)用SVD降噪具體過(guò)程如下：1)根據(jù)信號(hào)序列的點(diǎn)數(shù)N，確定構(gòu)造矩陣的行數(shù)為N/2或(N-1)/2，根據(jù)式(6)構(gòu)造Hankel矩陣并對(duì)構(gòu)造的矩陣進(jìn)行奇異值分解；2)對(duì)信號(hào)序列進(jìn)行FFT，判斷主頻率個(gè)數(shù)，若無(wú)法分辨則先對(duì)信號(hào)進(jìn)行自相關(guān)，再對(duì)自相關(guān)函數(shù)進(jìn)行FFT，判斷主頻率個(gè)數(shù)，通過(guò)得到的主頻率個(gè)數(shù)M，選定2M為有效秩的個(gè)數(shù)；3)根據(jù)式(4)對(duì)矩陣A進(jìn)行重構(gòu)，對(duì)重構(gòu)矩陣中對(duì)應(yīng)位置的元素相加平均后得到重構(gòu)信號(hào)，這就是降噪后的信號(hào)。

2 EMD分解信號(hào)原理

EMD算法假設(shè)對(duì)于任何信號(hào)都是由若干有限的本征模態(tài)函數(shù)(IMF)組成的，IMF必須滿足兩個(gè)條件，其一是在所有數(shù)據(jù)序列中，極值點(diǎn)個(gè)數(shù)和過(guò)零點(diǎn)個(gè)數(shù)必須相等或最多相差一個(gè)；其二是連接信號(hào)的局部極大值和極小值的包絡(luò)線均值為零，每一個(gè)IMF的獲得方法如下：

找到原始信號(hào)x(t)的所有極大值點(diǎn)，擬合出其極大值點(diǎn)包絡(luò)線e+(t)(通常采用三次樣條插值法)，同理找到原始信號(hào)x(t)的極小值點(diǎn)并擬合出其極小值包絡(luò)線e-(t)，原信號(hào)的均值包絡(luò)m1(t)定義為上下包絡(luò)線的均值，即：

(8)

(9)

這個(gè)ε一般取為0.2～0.3之間。

令：

r1(t)=x(t)-c1(t)

那么r1(t)就是原始信號(hào)x(t)中去除了高頻分量的剩余分量。由于仍然含有低頻成分，將r1(t)當(dāng)作新的信號(hào)重復(fù)上述操作，這樣可得到第二個(gè)基本IMF分量c2(t)，如此反復(fù)進(jìn)行，可得到第n個(gè)基本IMF分量cn(t)，若剩余分量rn(t)滿足分解停止條件之一時(shí)，可以認(rèn)為無(wú)法再篩分出基本IMF分量，從而停止EMD分解過(guò)程。

分解停止條件有如下幾種：1)當(dāng)剩余分量變得比預(yù)定值小時(shí)停止；2)當(dāng)剩余分量單調(diào)或只含兩個(gè)以下極值點(diǎn)時(shí)停止；3)對(duì)于有趨勢(shì)項(xiàng)的數(shù)據(jù)，剩余分量為一個(gè)平均的趨勢(shì)時(shí)停止。具體的EMD分解流程圖如圖4所示：

經(jīng)過(guò)EMD分解后，原始信號(hào)x(t)可表示為：

(10)

圖4 EMD分解流程圖

3 SVD與EMD檢測(cè)多頻小信號(hào)具體步驟

對(duì)于一個(gè)含有噪聲的多頻信號(hào)，通過(guò)上述方法檢測(cè)出其中有用成分的具體步驟如下：

1)根據(jù)第2節(jié)詳細(xì)敘述的SVD降噪過(guò)程對(duì)待檢測(cè)信號(hào)進(jìn)行降噪，使其信噪比滿足EMD分解的要求；

2)對(duì)SVD降噪后的信號(hào)進(jìn)行EMD分解，得到各個(gè)本征模態(tài)函數(shù)IMF，從而檢測(cè)出原信號(hào)中包含的各有用成分。

4 仿真驗(yàn)證與分析

信號(hào)經(jīng)過(guò)SVD降噪滿足EMD對(duì)信噪比的要求后，再經(jīng)EMD分解，即可得到信號(hào)中所包含的多頻成分。

為了檢驗(yàn)上述方法在不同信噪比情況下的檢測(cè)效果，對(duì)如下信號(hào)進(jìn)行仿真實(shí)驗(yàn)：

f(t)=cos(2pf1t)+cos(2pf2t)+n(t)

其中，f1=10Hz，f2=25Hz，采樣頻率為fs=250Hz，n(t)為Matlab中用randn函數(shù)產(chǎn)生的高斯白噪聲，在仿真過(guò)程中取不同的幅值，從而產(chǎn)生不同信噪比的混合信號(hào)。

圖5是信噪比為0.34dB時(shí)的分析效果，從圖中可以看出，自相關(guān)前我們無(wú)法從信號(hào)的FFT圖中判斷出主頻率個(gè)數(shù)，自相關(guān)后可以清晰的判斷出主頻率的個(gè)數(shù)為2，所以有效秩個(gè)數(shù)選為4，通過(guò)SVD降噪和EMD分解，得到檢測(cè)出的信號(hào)如圖5中(c)所示，從圖中可以看出，該方法有效的檢測(cè)出了原信號(hào)中包含的兩個(gè)正弦成分，而且檢測(cè)出的兩正弦信號(hào)幅值和頻率都與原信號(hào)中兩正弦成分的幅值和頻率具有較好的相合度，所以，在信噪比為0.34dB以上時(shí)該方法能夠得到較為理想的檢測(cè)效果。

圖5 信噪比為0.34dB時(shí)信號(hào)分析效果圖

降低待檢測(cè)信號(hào)信噪比，依照上述步驟繼續(xù)做仿真，得到了如圖6和圖7的仿真結(jié)果：

圖6為信噪比為-2dB時(shí)的檢測(cè)效果圖，從圖中可以看出，前兩個(gè)IMF分量分別為檢測(cè)出的原信號(hào)中兩個(gè)正弦分量，但此時(shí)幅值的偏差明顯比0.34dB時(shí)檢測(cè)的幅值偏差大，效果已經(jīng)開(kāi)始下降。

當(dāng)信噪比繼續(xù)降低，到達(dá)-5dB時(shí)，檢測(cè)效果如圖7所示，從圖中可以看出，此時(shí)已經(jīng)完全無(wú)法從IMF分量中判別出原信號(hào)中的兩個(gè)正弦分量，該方法失效。

綜合上述的仿真結(jié)果可以看出，本文提出的基于SVD與EMD相結(jié)合檢測(cè)多頻信號(hào)的方法在0.34dB以上信噪比情況下能夠有效的檢測(cè)出信號(hào)中所含有的多個(gè)頻率成分，同時(shí)檢測(cè)出的信號(hào)與原信號(hào)具有較好的相合度，充分運(yùn)用了SVD在降噪方面的特性和EMD在信號(hào)分解方面的優(yōu)勢(shì)，獲得了較好的檢測(cè)效果。

圖6 信噪比為-2dB時(shí)檢測(cè)效果

圖7 信噪比為-5dB時(shí)檢測(cè)效果

5 結(jié)束語(yǔ)

本文研究了一種基于SVD和EMD的多頻信號(hào)檢測(cè)方法，對(duì)于一個(gè)含有噪聲的信號(hào)，首先用SVD對(duì)其降噪，使其滿足EMD分解的要求，而后用EMD對(duì)降噪后的信號(hào)進(jìn)行分解，得到其中所包含的各個(gè)頻率成分，文章還著重討論了SVD中的兩個(gè)核心問(wèn)題，即重構(gòu)矩陣結(jié)構(gòu)問(wèn)題和有效秩個(gè)數(shù)問(wèn)題，提出了獲取重構(gòu)矩陣行數(shù)和有效秩個(gè)數(shù)的方法并獲得了較好的效果。最后的仿真表明，在信噪比高于0.34dB時(shí)，本文提出的方法能夠有效的檢測(cè)出信號(hào)中所含有的多個(gè)頻率成分，效果十分明顯。

[1] 趙東升.最小二乘法直線度評(píng)定的不確定度分析[J].計(jì)量技術(shù)，2011(5)

[2] 孟宗.基于EMD與AR譜的軋機(jī)主傳動(dòng)系統(tǒng)故障診斷研究[J].計(jì)量學(xué)報(bào)，2011，32(4)

[3] 王益艷.基于特征均值的SVD信號(hào)去噪算法[J].計(jì)算機(jī)應(yīng)用與軟件，2012，29(5)

[4] 呂志民．基于奇異譜的降噪方法及其在故障診斷技術(shù)中的應(yīng)用[J]．機(jī)械工程學(xué)報(bào)，1999，35(3)

[5] 程云鵬．矩陣論[M]．西安：西北工業(yè)大學(xué)出版社，2006

[6] 王維．基于動(dòng)態(tài)聚類(lèi)的奇異值分解降噪方法研究[J]．振動(dòng)工程學(xué)報(bào)，2008，21(3)

[7] 段向陽(yáng)．基于奇異值分解的信號(hào)特征提取方法研究[J]．振動(dòng)與沖擊，2009，28(11)

[8] 錢(qián)征文．利用奇異值分解的信號(hào)降噪方法[J]．振動(dòng)、測(cè)試與診斷，2011，31(4)

[9] 趙學(xué)智．矩陣構(gòu)造對(duì)奇異值分解信號(hào)處理效果的影響[J]．華南理工大學(xué)學(xué)報(bào)，2008，36(9)

[10] 高占晉．微弱信號(hào)檢測(cè)[M]．北京：清華大學(xué)出版社，2004

[11] Norden E．Huang，Man-li C．Wu,Steven R．Long，etc．A confidence limit for the empirical mode decomposition and Hilbert spectral analysis.The Royal Society．Proc．R．lond．A(2003)459，2317-2345

[12] Huang N E Shen Zheng Long S R et al The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis[J] Proc.R.Soc.Lond.1998 A:903-995