高成文， 孫云鵬

(1 92853部隊(duì)4分隊(duì) 興城 125106 2 92853部隊(duì)裝備部 興城 125106)

基于Hankel矩陣奇異值分解的速變參數(shù)去噪方法實(shí)現(xiàn)及仿真分析

高成文1， 孫云鵬2

(1 92853部隊(duì)4分隊(duì) 興城 125106 2 92853部隊(duì)裝備部 興城 125106)

針對(duì)現(xiàn)有降噪方法的不足，提出基于Hankel矩陣奇異值分解(SVD)的速變參數(shù)去噪方法。介紹奇異值分解及Hankel矩陣的相關(guān)理論，描述信號(hào)處理流程，并應(yīng)用Matlab對(duì)方法進(jìn)行實(shí)現(xiàn)。通過(guò)實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證表明，該方法對(duì)平穩(wěn)、非平穩(wěn)信號(hào)都具有較好的去噪效果，適用于速變參數(shù)的降噪處理。

遙測(cè)參數(shù)； Hankel矩陣； 奇異值分解； 去噪

引 言

速變參數(shù)通常指信號(hào)變化頻率高于10Hz的遙測(cè)參數(shù)(GJB 727—89)。在各類飛行器的飛行試驗(yàn)中，對(duì)該類隨機(jī)信息需應(yīng)用統(tǒng)計(jì)的方法進(jìn)行分析，描述它們?cè)跁r(shí)域、頻域的基本特性，為飛行器的改進(jìn)及性能分析提供可靠的依據(jù)。

由于飛行器上搭載的遙測(cè)分系統(tǒng)工作環(huán)境惡劣，測(cè)量及變換環(huán)節(jié)多，遙測(cè)信號(hào)在采集、轉(zhuǎn)換和傳輸過(guò)程中經(jīng)常受到設(shè)備、環(huán)境等因素的影響，致使實(shí)測(cè)信號(hào)中混有噪聲[1]。對(duì)染噪的遙測(cè)速變參數(shù)進(jìn)行降噪是數(shù)據(jù)預(yù)處理中的重要環(huán)節(jié)?；跀?shù)據(jù)平滑處理的降噪方法[2]，如五點(diǎn)三次平滑法，其優(yōu)點(diǎn)是對(duì)于時(shí)域數(shù)據(jù)能減少混入的高頻噪聲，對(duì)于頻域數(shù)據(jù)能使譜線變得光滑，但其也存在如下缺陷：頻域數(shù)據(jù)的多次平滑處理會(huì)使譜線的峰值降低、體態(tài)變寬，可造成識(shí)別參數(shù)誤差的增大。基于小波變換的降噪方法[3～5]，可將信號(hào)和噪聲正交分解于不同的頻率范圍中，它具有較強(qiáng)的濾除隨機(jī)信號(hào)噪聲的能力，不足之處在于：降噪效果與小波基及閾值的選取策略密切相關(guān)。本文提出一種基于Hankel矩陣奇異值分解(SVD)的遙測(cè)速變參數(shù)降噪方法，并應(yīng)用Matlab進(jìn)行編程實(shí)現(xiàn)。實(shí)驗(yàn)數(shù)據(jù)仿真分析表明，該方法對(duì)速變參數(shù)具有較強(qiáng)的去噪能力，實(shí)用性較好。

1 矩陣奇異值分解(SVD)的相關(guān)理論及Hankel矩陣

1.1 矩陣的奇異值及奇異值分解定理[6]

設(shè)A是m×n階實(shí)矩陣，稱n階方陣ATA的非零特征值的算術(shù)平方根為矩陣A的奇異值。

由于(ATA)T＝ATA，所以方陣ATA是實(shí)對(duì)稱矩陣，任取x∈Rn，記y＝Ax，則y∈Rm，且：xTATAx＝由此得ATA是半正定的，它的非零特征值都是正數(shù)。記它的全部非零特征值為λ1，λ2，…，λr，則矩陣A的奇異值為σ1＝它們皆取正值。

若記Im、In分別為m階和n階單位矩陣，由于：

所以n階方陣ATA和m階方陣AAT具有相同的非零特征值，因此，A的奇異值也可以定義為方陣AAT的非零特征值的算術(shù)平方根。顯然，若A對(duì)稱，則其奇異值就是它的非零特征值的絕對(duì)值。

矩陣奇異值分解定理：設(shè)A∈Rm×n，秩為r(r≤min(m，n))，則存在m階正交矩陣U和n階正交矩陣V，使得

這里，Σ＝diag(σ1，σ2，…，σr)，σ1，σ2，…，σr為A的奇異值。

1.2 Hankel矩陣[7]

Hankel矩陣是一類構(gòu)成特殊的矩陣。該型矩陣在數(shù)值分析、優(yōu)化理論、系統(tǒng)辯識(shí)等眾多領(lǐng)域具有廣泛的應(yīng)用。設(shè)離散數(shù)字序列x(i)的長(zhǎng)度為N，即i＝1，2，…，N，由x(i)構(gòu)造的Hankel矩陣具有如下形式：

式(3)中，m、n和N具有m+n-1＝N的關(guān)系。矩陣中的各行是由序列x(i)左移一個(gè)元素而形成的，每一行和前一行有n-1個(gè)元素是相同的，各行之間存在很大的相關(guān)性。

2 基于Hankel矩陣奇異值分解(SVD)的信號(hào)去噪方法及處理流程

文獻(xiàn)[8]從理論上研究了采用Hankel矩陣時(shí)SVD的信號(hào)分解本質(zhì)，從空間基的角度與小波分解過(guò)程進(jìn)行了剖析對(duì)比，分析了兩者在信號(hào)處理機(jī)理上的相似性。矩陣奇異值分解定理中A的左右正交矩陣U、V分別構(gòu)成了m維和n維的規(guī)范正交基。SVD分解結(jié)果是利用這兩個(gè)規(guī)范正交基實(shí)現(xiàn)的，結(jié)果反映的是原信號(hào)在規(guī)范正交基上投影的大小。基于奇異值分解(SVD)的信號(hào)去噪方法的基本思路是：根據(jù)真實(shí)信號(hào)和噪聲的不同特征，將含噪的測(cè)量信號(hào)構(gòu)成的Hankel矩陣分解成兩個(gè)互不相關(guān)的空間——真實(shí)信號(hào)空間和噪聲空間，利用相應(yīng)的方法對(duì)兩個(gè)空間的關(guān)聯(lián)矩陣——奇異值矩陣處理后，再重構(gòu)信號(hào)，實(shí)現(xiàn)去除測(cè)量信號(hào)噪聲的目的。

設(shè)實(shí)際測(cè)量信號(hào)序列x(i)具有如下形式：

其中，s(i)為真實(shí)信號(hào)，u(i)為噪聲。基于Hankel矩陣奇異值分解(SVD)的信號(hào)去噪方法的處理流程如下：

①利用實(shí)際測(cè)量信號(hào)x(i)構(gòu)造Hankel矩陣A，矩陣行數(shù)m、列數(shù)n和信號(hào)長(zhǎng)度N須滿足m+n-1＝N的關(guān)系。

②對(duì)矩陣A進(jìn)行奇異值分解，得到U?Rm×m、S?Rm×n和V?Rn×n，對(duì)應(yīng)關(guān)系為

③確定奇異值閾值(本文采用差分方法確定奇異值突變點(diǎn)作為閾值)，對(duì)奇異值矩陣S進(jìn)行處理，使小于閾值的奇異值為0，處理后的矩陣為Snew。

④利用U、V和Snew構(gòu)造新的矩陣Anew(Anew?Rm×n)。

⑤由矩陣Anew第一行的所有元素和第二行第n列至第m行第n列的m-1個(gè)元素，共同構(gòu)造長(zhǎng)度為N的新測(cè)量數(shù)據(jù)序列xnew(n)。

步驟③中提及的閾值處理方法有多種，如特征均值方法、奇異值中值方法及奇異熵增量方法等。文獻(xiàn)[7]驗(yàn)證：當(dāng)利用測(cè)量數(shù)據(jù)序列構(gòu)造的Hankel矩陣接近方陣時(shí)，去噪效果最為理想。

3 去噪效果檢驗(yàn)及分析

本文方法的核心是構(gòu)造的Hankel矩陣的形式及有效奇異值的確定方法。為檢驗(yàn)基于Hankel矩陣奇異值分解(SVD)的信號(hào)去噪方法的有效性，本文利用Matlab對(duì)該方法進(jìn)行了實(shí)現(xiàn)，并分別利用生成的平穩(wěn)和非平穩(wěn)仿真信號(hào)進(jìn)行了驗(yàn)證。構(gòu)造的Hankel矩陣A為512×513階。

①平穩(wěn)信號(hào)形式

本例中仿真信號(hào)由確定性信號(hào)x(n)＝sin(0.03·n)+sin(0.02·n)與服從N(0，σ2)分布的噪聲疊加而成，噪聲標(biāo)準(zhǔn)差σ＝0.8，采樣點(diǎn)數(shù)N＝1024，信噪比約為1.59dB。

圖1 信號(hào)、噪聲及含噪信號(hào)SVD分解奇異值分布

為分析奇異值的分布規(guī)律，分別對(duì)由純凈信號(hào)、噪聲及含噪信號(hào)構(gòu)成的Hankel矩陣進(jìn)行奇異值分解，結(jié)果如圖1所示。圖1(a)為疊加噪聲前信號(hào)奇異值分布，得到的奇異值矩陣的秩為4，奇異值幅值較大；圖1(b)為噪聲的奇異值分布，奇異值幅值較小，變化平緩；圖1(c)為含噪信號(hào)的奇異值分布。對(duì)比圖1中的三種奇異值分布可以發(fā)現(xiàn)，對(duì)于該類含噪信號(hào)，由反映信號(hào)的奇異值過(guò)渡到反映噪聲的奇異值時(shí)，奇異值在幅度上有突變。

圖2為仿真數(shù)據(jù)及去噪后信號(hào)重構(gòu)結(jié)果。其中，圖2(c)為利用本文方法確定的第5個(gè)奇異值作為閾值對(duì)奇異值矩陣進(jìn)行去噪處理后的信號(hào)重構(gòu)結(jié)果，圖2(d)為利用特征均值方法確定的第179個(gè)奇異值作為閾值時(shí)的處理結(jié)果。從圖中可以直觀地看出，圖2(c)和圖2(a)具有高度的相似性，表明本文方法對(duì)該類染噪信號(hào)具有很好的去噪效果。對(duì)比圖2(c)和圖2(d)可見(jiàn)，本文方法要優(yōu)于特征均值方法。

圖2 仿真數(shù)據(jù)及去噪后信號(hào)重構(gòu)結(jié)果

②非平穩(wěn)信號(hào)形式

圖3 加噪非平穩(wěn)信號(hào)去噪處理后重構(gòu)結(jié)果

圖3為加噪非平穩(wěn)信號(hào)去噪處理后重構(gòu)結(jié)果。其中，圖3(c)為前100個(gè)奇異值幅度分布情況，圖3 (d)為利用本文方法確定的第21個(gè)奇異值作為閾值對(duì)奇異值矩陣進(jìn)行去噪處理后的信號(hào)重構(gòu)結(jié)果，圖3 (e)為利用特征均值方法確定的第98個(gè)奇異值作為閾值時(shí)的處理結(jié)果。從圖中可以直觀地看出，圖3(d)的去噪效果較好，但個(gè)別波形在幅度上有失真；圖3(e)的去噪效果不如圖3(d)。剔除變頻信號(hào)中的噪聲比較困難，采用一般的降噪方法會(huì)造成真實(shí)信號(hào)的較大損失，本例說(shuō)明對(duì)于該類信號(hào)本文方法仍然適用。

通過(guò)數(shù)據(jù)仿真發(fā)現(xiàn)，對(duì)于含噪平穩(wěn)信號(hào)采用本文方法確定奇異值閾值重構(gòu)信號(hào)濾除噪聲的效果很好，而文獻(xiàn)[9]提出的基于特征均值的SVD信號(hào)去噪算法的去噪效果略遜一籌。對(duì)于含噪非平穩(wěn)信號(hào)，當(dāng)噪聲標(biāo)準(zhǔn)差σ＞0.18時(shí)，采用本文方法和特征均值方法處理奇異值，重構(gòu)信號(hào)波形都出現(xiàn)了失真。

下面對(duì)本文方法、傳統(tǒng)濾波方法及小波閾值去噪方法的性能進(jìn)行比較。從噪聲尤其是白噪聲頻率分布規(guī)律的角度分析，應(yīng)用低通或高通濾波都不能有效根本濾除噪聲。從信號(hào)和噪聲能量分布的角度分析，真實(shí)信號(hào)能量比較集中而噪聲能量比較分散，小波分解正是根據(jù)信號(hào)和噪聲的能量在小波系數(shù)上分布的差異，通過(guò)對(duì)分解后各層小波空間的細(xì)節(jié)系數(shù)進(jìn)行閾值處理來(lái)實(shí)現(xiàn)信號(hào)去噪，閾值的確定方法直接影響去噪效果。利用SVD降噪，同樣是基于信號(hào)和噪聲能量考慮，對(duì)由含噪測(cè)量數(shù)據(jù)構(gòu)成的Hankel矩陣進(jìn)行奇異值分解時(shí)，得到的奇異值序列中反映信號(hào)的奇異值較大，反映噪聲的奇異值較小，通過(guò)將反映噪聲的奇異值置零，再重構(gòu)矩陣，實(shí)現(xiàn)信號(hào)降噪。對(duì)比能量分布規(guī)律和處理策略，本文方法較優(yōu)。

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

經(jīng)Matlab編程實(shí)現(xiàn)及數(shù)據(jù)仿真驗(yàn)證表明，基于Hankel矩陣奇異值分解(SVD)的遙測(cè)速變參數(shù)去噪方法對(duì)平穩(wěn)、非平穩(wěn)信號(hào)都具有較強(qiáng)的去噪能力，具體體現(xiàn)在：①對(duì)于含噪平穩(wěn)信號(hào)，本文方法去噪效果很好(當(dāng)信號(hào)僅含有單頻確定性信號(hào)和高斯白噪聲時(shí)，信噪比為0dB時(shí)的去噪效果仍然很理想)；②對(duì)于含噪非平穩(wěn)信號(hào)，本文方法的去噪效果較好，但當(dāng)信噪比低于12.90dB時(shí)，重構(gòu)信號(hào)在幅度上有些許失真。本文也進(jìn)一步驗(yàn)證了以奇異值突變點(diǎn)作為閾值更具實(shí)用性。本文提出的基于Hankel矩陣奇異值分解(SVD)的去噪方法便于軟件實(shí)現(xiàn)，采用差分方式確定奇異值閾值方法簡(jiǎn)單易行，適用于遙測(cè)速變參數(shù)的降噪處理。

[1] 中國(guó)人民解放軍總裝備部.遙測(cè)數(shù)據(jù)處理[M].北京：國(guó)防工業(yè)出版社，2002.

[2] 王 濟(jì)，胡 曉.Matlab在振動(dòng)信號(hào)處理中的應(yīng)用[M].北京：中國(guó)水利水電出版社，2006.

[3] Azzalini A，F(xiàn)arge M，Schneider K.Nonlinear Wavelet Thresh-olding：a Recursive Method to Determine the Optimal Denoising Threshold[J].App.1 Comput.Harmon.Ana，2005，18：177～185.

[4] 張吉先，鐘秋海.小波門(mén)限消噪法應(yīng)用中分解層數(shù)及閾值的確定[J].中國(guó)電機(jī)工程學(xué)報(bào)，2004，2(24)：118～122.

[5] 程發(fā)斌，湯寶平，鐘佑明.基于最優(yōu)Morlet小波和SVD的濾波消噪方法及故障診斷的應(yīng)用[J].振動(dòng)與沖擊，2008，27(2)：91～94.

[6] 蔣長(zhǎng)錦.線性代數(shù)計(jì)算方法[M].合肥：中國(guó)科技大學(xué)出版社，2003.

[7] 張 波，李健君.基于Hankel矩陣與奇異值分解(SVD)的濾波方法及在飛機(jī)顫振試驗(yàn)數(shù)據(jù)預(yù)處理中的應(yīng)用[J].振動(dòng)與沖擊，2009，28(2)：162～166.

[8] 趙學(xué)智，葉邦寬.SVD和小波變換的信號(hào)處理效果相似性及其機(jī)理分析[J].電子學(xué)報(bào)，2008，36(8)：1582～1588.

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

Simulation and Implementation of Denoising Method for Fast Changing Parameter Data Based on Hankel Matrix and SVD

Gao Chengwen， Sun Yunpeng

Aiming at the weakness of existing denoising methods，this paper puts forward a denoising method for the fast changing parameter data based on the Hankel matrix and SVD.The correlative theories of the Hankel matrix and SVD are introduced，the flow of signal processing is described，and themethod is realized in the Matlab software.The data of experiments prove that this method has good effects for the stationary signal and non-stationary signal，and is applicable to the denoising processing of fast changing parameter data.

Telemetry data； Hankelmatrix； SVD； Denoising

TP393

A

CN11-1780(2014)01-0060-05

2013-05-08 收修改稿日期：2013-06-27

高成文 1966年生，學(xué)士，高級(jí)工程師，主要從事飛行器試驗(yàn)遙測(cè)數(shù)據(jù)處理工作。

孫云鵬 1982年生，碩士，工程師，主要從事測(cè)試裝備總體及管理工作。