用SVD實現(xiàn)近頻諧波提取

張曉濤

（廣東機電職業(yè)技術(shù)學(xué)院，廣東 廣州 510515）

近頻諧波分析表征了一類帶有共性的信號處理需求，具有廣泛應(yīng)用背景。當(dāng)兩臺工作轉(zhuǎn)速相近的水泵、風(fēng)機等動力設(shè)備在相鄰基礎(chǔ)上運行時，二者的振動會相互傳遞而產(chǎn)生拍振[1]，導(dǎo)致所采集信號存在近頻諧波；生物信號[2-3]的頻率分布范圍與工頻范圍重疊，當(dāng)生物信號受工頻干擾時，生物信號與工頻形成近頻諧波；工頻干擾常與地震信號有效范圍重疊，所形成的近頻諧波影響地震資料疊加和偏移成像質(zhì)量[4]；許多動平衡試驗的工作轉(zhuǎn)速在10Hz甚至低至1Hz，所采集信號中包含成分復(fù)雜的干擾頻率，干擾頻率對動不平衡信號形成近頻干擾[5]。

SVD具有穩(wěn)定性好、結(jié)果無相移、波形失真小等優(yōu)點，在工頻干擾消除[6]、信號降噪[7]、特征提取[8]、過電壓辨識[9]等諸多領(lǐng)域獲得應(yīng)用，但在近頻諧波提取方面研究則鮮有報道。本文研究SVD在近頻諧波提取上的應(yīng)用，包括奇異值與諧波參數(shù)的關(guān)系、重構(gòu)誤差與諧波幅值的關(guān)系、重構(gòu)誤差與頻率差的關(guān)系等。將其用于油膜支承可傾瓦軸承振動信號近頻諧波提取，提取出完整的工頻基波，從而消除其對振動信號特征分析的影響。

1 SVD原理[10]

若實矩陣表示為A∈Rm×n，則存在正交矩陣U∈Rm×m和V∈Rn×n，使下式成立：

式中，U為m×m維正交矩陣，V為n×n維正交矩陣；∑=diag（σ1，σ2，…，σr）為對角矩陣，對角元素σ1，σ2，…，σr為按降序排列的奇異值，即σ1≥…≥σr＞0，r=min（m，n）是A的秩。

因∑為對角矩陣，可將秩為r的m×n階矩陣A表示為r個秩為1的m×n階子矩陣的和，每個子矩陣由兩個奇異向量和權(quán)值相乘得到，如式（2）所示：

式中，ui、vi分別為U、V的第i個列奇異向量，σi是A的第i個奇異值。應(yīng)用中通常根據(jù)式（1）的有效奇異值個數(shù)得到A的估計值，以實現(xiàn)特征提取或降噪。

2 近頻諧波提取的SVD方法分析

2.1 奇異值與諧波參數(shù)的關(guān)系

為提取諧波，需求出幅值、頻率、相位等諧波參數(shù)，諧波參數(shù)通常具有隨機性。SVD在諧波提取上具有獨特優(yōu)勢[11]，表現(xiàn)為奇異值只與諧波幅值有關(guān)而與其頻率、相位無關(guān)，在應(yīng)用中只需關(guān)注諧波幅值一項，從而降低諧波提取難度。

為說明SVD在諧波提取上的應(yīng)用，采用仿真信號進(jìn)行分析。仿真信號用x表示，設(shè)α=3.1、f=45Hz、Φ=0.3π，信號長度為1024。

將x構(gòu)造為512×513維的Hankel矩陣并SVD，奇異值序列中第1個奇異值σ1為796.24，第2個奇異值σ2為792.51，其余奇異值近似為零。

為說明奇異值與諧波幅值的關(guān)系，將α設(shè)為變量，而將f、Φ設(shè)為常量，取f=45Hz、Φ=0.3π。奇異值σ1同幅值α的關(guān)系如圖1所示，σ1與α呈線性關(guān)系，這表明奇異值同諧波幅值成正比。

圖1 奇異值與諧波幅值的關(guān)系

為說明奇異值與諧波頻率的關(guān)系，將f設(shè)為變量，而將α、Φ為常量，取α=3.1、Φ=0.3π。奇異值σ1同頻率f的關(guān)系如圖2所示，在[0，499]Hz范圍內(nèi)，改變頻率f的值，但奇異值σ1基本不變，這表明奇異值同頻率無關(guān)。

圖2 奇異值同諧波頻率的關(guān)系

為說明奇異值與諧波相位的關(guān)系，將Φ設(shè)為變量，將f、α為常量，即取f=45Hz，α=3.1。奇異值σ1同相位Φ的關(guān)系如圖3所示，改變相位Φ，但σ1基本不變，這表明奇異值同相位無關(guān)。

圖3 奇異值同諧波相位的關(guān)系

由上可知奇異值只與諧波幅值有關(guān)，而與諧波頻率和相位無關(guān)，奇異值與諧波幅值成線性關(guān)系，此結(jié)論對第2個奇異值σ2仍然成立。

一個諧波對應(yīng)于兩個連續(xù)分布的奇異值[12]，將這兩個奇異值置零，再作SVD逆變換就可以將該諧波消除。如果只保留這兩個奇異值而將其余奇異值置零，則可以將該諧波提取出來。如果信號包含多個諧波，當(dāng)各諧波的幅值相差較大時，與諧波相對應(yīng)的奇異值會與其余奇異值存在較大的區(qū)分度，而當(dāng)各諧波具有相同或相近的幅值時，顯然就無法通過奇異值將其區(qū)分開。

2.2 當(dāng)諧波幅值相近時提取結(jié)果分析

設(shè)仿真信號為x（t），其包含兩個諧波x1、x2，因奇異值同相位無關(guān)，故設(shè)Φ1、Φ2為0，同時設(shè)f1為49Hz，f2為152Hz。

設(shè)α1為定值1，α2取圖4中所示的若干值，采用SVD方法從x（t）分離出諧波x1，分離結(jié)果用表示。重構(gòu)誤差Δx1與α2的關(guān)系如圖4所示，重構(gòu)誤差定義為：

圖4 重構(gòu)誤差與幅值的關(guān)系

式中t∈[1，1024]，N=1024。

在圖4中，當(dāng)α2為1時Δx1最大，這表明當(dāng)x2和x1的幅值相同時，x2對x1的影響最大。當(dāng)α2值為1時，的頻譜如圖5所示，圖中既有諧波x1的頻譜（49Hz），同時也有諧波x2的頻譜（152Hz），這從另一個角度表明，如果重構(gòu)結(jié)果中存在多個諧波成分，則這些諧波的幅值可能相同或相近。

圖5 重構(gòu)諧波頻譜

2.3 用EMD分離諧波

由圖5知，當(dāng)諧波幅值相同或相近時，SVD無法將二者完全分離。此時可將與諧波相對的全部奇異值保留，再作SVD逆變換，以得到包含這些諧波的時域信號，該信號為降噪信號。經(jīng)驗?zāi)Ｊ椒纸鈁13]（Empirical Mode Decomposition，簡稱EMD）是一種自適應(yīng)的信號分析方法，可根據(jù)時間尺度實現(xiàn)信號成分分離，但若信號含噪聲，則分解結(jié)果會出現(xiàn)模態(tài)混疊。前述時域信號為降噪信號，可避免因噪聲而出現(xiàn)模態(tài)混疊，故用EMD方法分析該信號具有合理性。

設(shè)式（4）仿真信號x（t）的諧波幅值α1、α2都為1，其余參數(shù)不變。兩個諧波對應(yīng)4個奇異值，保留奇異值序列中與其相應(yīng)的4個奇異值，得到包含兩個諧波的時域信號，該信號的頻譜如圖6所示。應(yīng)用EMD方法將該信號分解為兩個諧波，分離的諧波頻譜如圖7所示，兩個諧波被有效分離出來。

圖6 重構(gòu)時域信號頻譜

圖7 EMD分解結(jié)果

2.4 諧波頻率差與重構(gòu)誤差的關(guān)系

當(dāng)兩諧波的頻率相近時，重構(gòu)結(jié)果可能會受到影響，因此有必要分析重構(gòu)誤差與頻率差的關(guān)系。

式（4）所示的仿真信號，設(shè)諧波參數(shù)為α1=2、f1=50Hz、α2=1、f2=f1+Δf（Hz），頻 率 差Δf取0.05～7.5范圍內(nèi)的若干值。x1的重構(gòu)誤差Δx1與頻率差Δf的關(guān)系如圖8所示，當(dāng)頻率差小于1.8Hz時，頻率差越小則重構(gòu)誤差越大，當(dāng)頻率差為零時重構(gòu)誤差最大；當(dāng)頻率差由1.8Hz逐漸增大時，重構(gòu)誤差逐漸減小，并隨著頻率差的增加，重構(gòu)誤差逐漸趨于穩(wěn)定。

圖8 重構(gòu)誤差與頻率差的關(guān)系

3 工程信號分析

油膜支承可傾瓦軸承[14]具有獨特雙層油膜結(jié)構(gòu)，同常規(guī)軸承相比具有優(yōu)異的減振特性。為開展試驗搭建了軸承－轉(zhuǎn)子試驗臺，采用電渦流位移傳感器測量振動，用西門子LMS系統(tǒng)采集振動信號。振動信號受工頻干擾影響明顯，研究采用SVD方法將其提取出來，以消除其對有效特征分析的影響。

所研究的油膜支承可傾瓦軸承具有特殊的用途，其在2730r/min（即45.5Hz）轉(zhuǎn)速下所表現(xiàn)出的減振特性是研究重點。以1024Hz采樣率采集1024點的信號，信號頻譜如圖9（a）所示。為進(jìn)一步觀察工頻諧波特征，將頻譜圖進(jìn)行放大，結(jié)果如圖9（b）所示，可見工頻基波較強。正常電網(wǎng)工頻基波波動范圍為50±0.2Hz，而信號基頻為45.5Hz，二者頻率差約為4.5Hz，構(gòu)成近頻諧波關(guān)系。

圖9 振動信號波形及其頻譜

將振動信號構(gòu)造為512×513維的Hankel矩陣，并作SVD。由2.1節(jié)有關(guān)結(jié)論，奇異值序列中第1個奇異值對應(yīng)于直流分量，其波形如圖10（a）所示；第2、3個奇異值對應(yīng)于基頻分量，其頻譜如圖10（b）所示；第4、5個奇異值對應(yīng)于二倍頻分量，其頻譜如圖10（c）所示；第6、7個奇異值對應(yīng)于工頻基波分量，如圖10（d）所示，圖中既有基波也有三次諧波。由2.2節(jié)有關(guān)結(jié)論，可斷定基波幅值和三次諧波幅值是相同或相近的，保留第6 ～ 9個奇異值，而將其他奇異值置零，再作SVD逆變換，得到包含基波和三次諧波的時域信號，其頻譜如圖11所示，圖中可見基波和三次諧波的幅值非常接近。對該時域信號作EMD分解，得到圖12所示的基波波形及頻譜，所提取的基波清晰完整。

圖10 信號各分量頻譜

圖11 時域信號頻譜

圖12 基波波形及頻譜

采用Prony算法[15]對基波波形作擬合，擬合結(jié)果如圖13所示，圖中只給出了0～0.1s范圍的波形，實線表示擬合前的波形，而虛線表示擬合后的波形，二者擬合度較高。由該算法求得基波幅值為0.0042mm，初相位為-1.4650rad，頻率為49.8562Hz，三次諧波幅值為0.0044mm，初相位為1.5409rad，頻率為149.4170Hz，基波幅值和三次諧波幅值幾乎相等。

圖13 基波擬合前、后波形

4 結(jié)論

（1）利用奇異值與諧波幅值成線性關(guān)系，而與頻率和相位無關(guān)的特性，將SVD用于油膜支承可傾瓦軸承振動信號工頻諧波提取，提取出的基波完整，基波頻率為49.8562Hz，該方法特點是無需精確求解諧波參數(shù)；

（2）諧波幅值越接近則重構(gòu)誤差越大，當(dāng)幅值相同時重構(gòu)誤差最大；

（3）諧波重構(gòu)誤差同頻率差有關(guān)，當(dāng)頻率差大于1.8Hz時重構(gòu)誤差較小并趨于穩(wěn)定。該頻率差是基于特定矩陣結(jié)構(gòu)下得到的，在其他矩陣結(jié)構(gòu)下該值可能不同。