李 穎， 盧洪超， 周 琳， 陳文文， 齊聰山， 劉福順

(1. 浙江科技學(xué)院 中德工程師學(xué)院， 杭州 310023； 2. 中國海洋大學(xué) 工程學(xué)院， 青島 266100；3. 海洋石油工程(青島)有限公司， 山東 青島 266520)

對結(jié)構(gòu)實測振動響應(yīng)信號進行分析是獲取結(jié)構(gòu)動力特性的重要手段。然而，在土木工程、海洋工程領(lǐng)域，測試設(shè)備、環(huán)境等因素導(dǎo)致測試的結(jié)構(gòu)振動信號中不可避免的包含各種噪聲。信號中的噪聲不僅使結(jié)構(gòu)模態(tài)識別的精度降低，還會產(chǎn)生“噪聲模態(tài)”[1]。當(dāng)實測信號中包含可能淹沒結(jié)構(gòu)固有模態(tài)的高背景噪聲時，結(jié)構(gòu)模態(tài)參數(shù)識別十分困難。為此，在對工程結(jié)構(gòu)進行模態(tài)參數(shù)識別、結(jié)構(gòu)損傷診斷和健康檢測時，為了保證分析結(jié)果的正確性，往往需要采取合理的技術(shù)手段降低測試噪聲的影響[2]。

目前，已經(jīng)有許多信號消噪技術(shù)應(yīng)用于工程實踐。Cadzow[3]提出了基于奇異值分解(Singular Value Decomposition, SVD)的低秩逼近算法，該方法運用結(jié)構(gòu)響應(yīng)數(shù)據(jù)構(gòu)建Hankel矩陣，進行SVD分解，認為較小的奇異值對應(yīng)噪聲，將其置為零并重構(gòu)Hankel矩陣，從而達到消噪的目的。該方法對信噪比(Signal-Noise Ratio, SNR)較小的信號消噪效果較好，當(dāng)噪聲能量較大時，噪聲成分對應(yīng)的奇異值可能比結(jié)構(gòu)固有成分對應(yīng)的奇異值更大，從而表現(xiàn)出消噪局限性；另外，當(dāng)響應(yīng)數(shù)據(jù)較大時，導(dǎo)致Hankel矩陣的維數(shù)較大，進行SVD分解會消耗大量計算機資源[4]，甚至無法完成計算。Shumway等[5]對時域平均法作了詳細闡述，時域平均法用一段數(shù)據(jù)的平均值代替中心點的數(shù)值，進行平均的方式包括直接平均、多項式平均和周期回歸平均等方法，并將各種方法應(yīng)用于時間序列，對時間序列進行光滑。時域平均法具有一定的消噪效果，但各種平均方法的選取缺乏參照，沒有對消噪效果進行定量討論，另外運用該方法還需要時標(biāo)信息，而且要有足夠的數(shù)據(jù)量。Kim等[6]在頻域內(nèi)對頻響函數(shù)(Frequency Response Function, FRF)運用小波分解進行消噪，并將該方法運用于2自由度質(zhì)量-彈簧-阻尼系統(tǒng)，能夠有效消除信號中添加的高斯白噪聲，但小波消噪算法對小波核函數(shù)的選擇有較大依賴性，在實際運用中受到一定局限。Wu等[7-8]對經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decompositon, EMD)的研究表明，白噪聲經(jīng)EMD分解得到的固有模態(tài)函數(shù)(Intrinsic Mode Function, IMF)分量的統(tǒng)計特征滿足正態(tài)分布，IMF的能量密度與對應(yīng)的平均周期的積為常量，這些研究結(jié)果為EMD消噪提供了理論基礎(chǔ)。對含噪聲的信號進行EMD分解，得到頻率從高到低的IMF，高頻部分受噪聲影響較大，低頻部分受噪聲影響較小，而信號的能量也主要集中在低頻段，使用低頻段的IMF進行重構(gòu)完成消噪過程[9]。基于EMD分解的消噪方法對模擬的高信噪比信號具有較好的消噪效果，但該方法可能舍棄部分有用信號的能量，而且無法準(zhǔn)確判定選取有效IMF的個數(shù)。另外，近年來還有學(xué)者運用EMD分解方法對非平穩(wěn)信號的消噪進行消噪研究，并應(yīng)用于轉(zhuǎn)子振動信號的去噪，取得了較好效果[10]。

本文將從信號分解的角度出發(fā)，提出一種基于極值-留數(shù)的信號分解消噪方法，該方法提取代表信號中固有信息的極值和留數(shù)，并根據(jù)這些信息進行重構(gòu)，從而達到消噪目的。與傳統(tǒng)方法相比，本文方法避免反復(fù)進行SVD分解運算，具有較高的計算效率。另外，進行信號分解時，無需積分變換的核函數(shù)，便于實施。該方法能夠解決實際工程中在遭遇高背景噪聲時，其結(jié)構(gòu)固有信息往往被當(dāng)做環(huán)境噪聲而剔除的問題。文中將首先應(yīng)用一質(zhì)量-彈簧-阻尼模型進行驗證，然后應(yīng)用實測海洋平臺加速度響應(yīng)數(shù)據(jù)進一步研究新方法在實際工程中的應(yīng)用效果。

1 信號分解方式

1.1 傅里葉序列

連續(xù)函數(shù)y(t)的周期為T，其傅里葉分解為

(1)

(2)

對于實測的離散數(shù)字信號yk(k=0,2,…,N-1)，其傅里葉分解為

(3)

其中xn=ei(2πn/N)。

1.2 Prony序列

連續(xù)信號y(t)分解成復(fù)指數(shù)的形式為[11]

(4)

式中:p為分解階次;λn為極值，由于y(t)為實數(shù)，λn必須為實數(shù)或為共軛復(fù)數(shù)對。令λn≡-αn+iωn，αn為阻尼系數(shù)，ωn為角頻率。系數(shù)γn為留數(shù)，它與λn成對出現(xiàn)，或為實數(shù)，或為共軛復(fù)數(shù)。令γn≡Aneiθn，其中An為各個成分的幅值，θn為對應(yīng)相位。

等間隔采樣離散信號yk，其復(fù)指數(shù)分解為

(5)

式中:zn=eλnΔt。

2 基于極值-留數(shù)的信號降噪方法

對信號進行Prony分解，得到一系列極值和留數(shù)，每個極值和對應(yīng)的留數(shù)表示信號中的一個成分。通過式(4)進行疊加，即可用極值和留數(shù)表示實測信號。當(dāng)信號中包含噪聲成分時，噪聲成分也有對應(yīng)的極值和留數(shù)，本文方法就是通過剔除噪聲成分的極值和留數(shù)，從而達到消噪的目的。因此本文消噪方法的關(guān)鍵在于準(zhǔn)確求解信號所對應(yīng)的極值和留數(shù)。

為了求解式(4)中的極值和留數(shù)，引入高階微分方程

(6)

(7)

式(4)為式(6)的通解，其中λn為式(6)對應(yīng)特征方程式(7)的根。

對于實測的等間隔采樣離散數(shù)據(jù)yk(k=0,1,…,N-1)，式(6)高階微分方程變?yōu)椴罘址匠?/p>

(8)

式中:bn為常系數(shù)。令bn=1，式(9)對應(yīng)的特征方程為

(9)

通過以下步驟求解極值和留數(shù)

步驟1將式(8)寫成矩陣形式，由于Y和y′已知，可求出向量b，即系數(shù)bk(k=0,1,…,p-1)

Yb=-y′

(10)

其中

(11)

(12)

(13)

步驟2將bk代入式(9)中，求得特征值zk(k=1,2,…,p)

步驟3將特征值zk代入式(5)中，即可求得對應(yīng)的留數(shù)γk。

2.1 一階差分方程求解方法

由于求解方程式(10)時，存在病態(tài)問題，很小的誤差會使結(jié)果產(chǎn)生較大的偏差。為了解決病態(tài)問題，本文將高階差分方程轉(zhuǎn)化為一階差分方程，即令

x1,k=yk,x2,k=yk+1,…,xp,k=yk+p-1

(14)

則

(15)

式(8)可以寫成

xk+1=Fxk

(16)

其中

(17)

式(16)表明，其特征值即是方程式(7)的根。故可由響應(yīng)數(shù)據(jù)yk構(gòu)建Hankel矩陣

(18)

式中:ξ、η分別為Hankel矩陣的行數(shù)和列數(shù)，為了提高計算效率，ξ和η的值都取最接近于數(shù)據(jù)序列長度一半的整數(shù)。當(dāng)k分別取0和1時，對H(0)進行奇異值分解，得

(19)

(20)

由式(21)得到系統(tǒng)矩陣A

(21)

對A進行特征值分析，得到特征值zn，則響應(yīng)信號的極值λn=ln(zn)/Δt，將λn代入式(5)即可求得對應(yīng)的留數(shù)γn。

2.2 噪聲信號分離與剔除

海洋結(jié)構(gòu)實測振動數(shù)據(jù)y(n)可以看作由兩部分組成[12]

y(n)=yreal(n)+ynoise(n)

(22)

式中:yreal(n)為結(jié)構(gòu)的真實響應(yīng)信號;ynoise(n)為實測信號中的噪聲信號。

對實測信號進行極值-留數(shù)分解，得到極值和留數(shù)

(23)

(24)

去除噪聲對應(yīng)的極值和留數(shù)后，得到實真響應(yīng)信號的極值和留數(shù)

(25)

(26)

根據(jù)式(5)重構(gòu)響應(yīng)信號，得到真實信號

(27)

極值-留數(shù)分解消噪算法步驟如下：

步驟1對海洋結(jié)構(gòu)實測振動信號y(n)進行傅里葉變換，得到響應(yīng)信號的頻譜；

步驟2根據(jù)頻譜找到峰值對應(yīng)的頻率，并以此頻率為中心頻率設(shè)置頻率窗口；

步驟3對實測信號y(n)進行極值-留數(shù)分解，取極值的絕對值，將復(fù)頻率轉(zhuǎn)化為實頻率，選出步驟二頻率窗口中的實頻率對應(yīng)的極值和留數(shù)；

步驟4通過式(5)，將步驟三得到極值和留數(shù)重構(gòu)成時域信號，完成消噪過程。

3 數(shù)值算例

為了驗證本文方法的有效性，構(gòu)建一個6自由度質(zhì)量-彈簧-阻尼模型，如圖1所示。系統(tǒng)中各質(zhì)量塊的質(zhì)量為mi=100 kg，剛度為ki=5×107N/m，阻尼系數(shù)為ci=50 N·s/m，各質(zhì)量塊的位移分別為xi，其中i=1,2,…,6。

圖1 6自由度質(zhì)量-彈簧-阻尼模型

根據(jù)上述模型的系統(tǒng)狀態(tài)空間方程，利用MATLAB軟件的impulse函數(shù)生成離散脈沖響應(yīng)數(shù)據(jù)。脈沖響應(yīng)的采樣頻率設(shè)為500 Hz，信號序列包含2 048個點。對于不同的激勵點，整個6自由度系統(tǒng)會產(chǎn)生36個響應(yīng)序列，為了便于分析，選取在質(zhì)量m1處輸入、m1處輸出的脈沖響應(yīng)數(shù)據(jù)進行研究，其時程如圖2所示。

圖2 分析數(shù)據(jù)的時程圖

為了研究消噪效果，在信號中添加白噪聲模擬實測信號中的噪聲。使用SNR來描述噪聲水平，SNR為信號功率Ps與噪聲功率Pn之比，即SNR=10lg(Ps/Pn)。對分析數(shù)據(jù)分別添加不同信噪比的白噪聲，模擬實測信號中的不同強度噪聲。

對含信噪比為40 dB噪聲的信號和不含噪聲的信號進行傅里葉分析，得到信號的頻譜，它們的時程和頻譜對比如圖3所示。

圖3 不含噪聲信號與含信噪比40 dB噪聲的信號時程與頻譜對比

Fig.3 The time domain and frequency spectrum comparison between clear signal and signal containing SNR 40 dB noise

根據(jù)頻譜的峰值，分別設(shè)置頻率范圍(26.5，27.5)、(79，80)、(127，128)、(168，169)、(198.5，199.5)和(218，219)，6個峰值分別包含于6個頻率范圍之內(nèi)。將信噪比為40 dB的信號代入式(18)構(gòu)建Hankel矩陣，根據(jù)式(21)得到系統(tǒng)矩陣A，進而計算極值和留數(shù)。對極值取模獲得實頻率，根據(jù)實頻率選取各頻率范圍內(nèi)對應(yīng)的極值和留數(shù)，代入式(5)中進行重構(gòu)，得到消噪后信號。消噪后信號與不含噪聲信號的時程和頻譜對比,如圖4所示。結(jié)果表明本文方法能有效消除信噪比為40 dB的噪聲。

為了進一步驗證本文方法消噪效果，將該方法應(yīng)用于噪聲水平較高的信噪比為10 dB的信號，消噪前后的時程和頻譜如圖5和圖6所示。消噪前第六階頻率幾乎淹沒于噪聲中，消噪后得到的頻譜與不含噪聲的頻譜重合較好，表明本文方法對高噪聲水平的信號也有較好的消噪效果。

圖4 不含噪聲信號與含信噪比40 dB噪聲的信號消噪后的時程與頻譜對比

Fig.4 The time domain and frequency spectrum comparison between clear signal and signals eliminated noise

圖5 不含噪聲信號與含信噪比10 dB噪聲的信號時程與頻譜對比

Fig.5 The time domain and frequency spectrum comparison between clear signal and signal containing SNR 10 dB noise

圖6 不含噪聲信號與含信噪比10 dB噪聲的信號消噪后的時程與頻譜對比

Fig.6 The time domain and frequency spectrum comparison between clear signal and signals eliminated noise

為了定量分析本文方法的消噪效果，引入噪聲抑制率(Noise Suppression Ration, NSR)、信號失真率(Signal Distortion Ratio, SDR)[13]和決定系數(shù)(Coefficient of Determination)三個指標(biāo)。

NSR=

(28)

(29)

式中：s(n)為不含噪聲信號；f(n)為添加噪聲后f′(n)為消噪后信號。

(30)

對于信噪比為40 dB的信號，本文消噪方法消噪后的噪聲抑制率為0.797 8，信號失真率為0.002 008 4，決定系數(shù)為0.999 96。運用本文方法分別對含信噪比為30 dB、20 dB和10 dB噪聲的信號進行消噪，其噪聲抑制率和信號失真率列于表1中，消噪前后的決定系數(shù)如圖7所示。對比文獻[13]，其提出的消噪方法的噪聲抑制率為0.872 7，信號失真率為0.174 4。本文方法消噪后得到不同噪聲水平的噪聲抑制率都約為0.8，信號失真率隨著噪聲的增大而逐漸增大，但總都維持在比較低的水平；消噪后的決定系數(shù)R2都大于0.995，隨著噪聲水平的提高，決定系數(shù)R2逐漸減小。三個指標(biāo)的分析結(jié)果表明本文的消噪方法消噪效果較好。

表1 不同信噪比信號的噪聲抑制率(NSR)和信號失真率(SDR)

4 實測數(shù)據(jù)應(yīng)用

海洋結(jié)構(gòu)物實測振動響應(yīng)信號往往含有較高水平的背景噪聲，尤其環(huán)境激勵下其結(jié)構(gòu)固有模態(tài)往往淹沒于噪聲之中，導(dǎo)致結(jié)構(gòu)動力特性識別困難。為研究方法的有效性，選用冰擊激勵作用下的某海洋石油平臺(見圖8)實測數(shù)據(jù)進行研究。結(jié)構(gòu)振動測試時，使用三向加速度傳感器采集其加速度響應(yīng)信號,傳感器安裝位置如圖9所示。采樣頻率為200 Hz。選取主振動方向300～320 s的實測數(shù)據(jù)進行分析，其加速度響應(yīng)時程如圖10所示。對加速度響應(yīng)數(shù)據(jù)進行傅里葉分析，得到頻譜如圖11所示。由于海洋平臺實際環(huán)境激勵情況比較復(fù)雜，激勵能量較弱，得到頻譜非常復(fù)雜，包含大量噪聲成分。這些噪聲使響應(yīng)信號的進一步分析面臨很大困難。

圖7 決定系數(shù)隨信噪比變化

圖8 海洋石油平臺

圖9 加速度傳感器安裝位置

1994年，該平臺所有單位對其進行了振動檢測，應(yīng)用傳統(tǒng)的頻譜分析技術(shù)，得到該石油平臺的前2階頻率，分別為0.9 Hz、1.15 Hz與1.175 Hz?？紤]到該平臺至今仍在服役，其結(jié)構(gòu)固有特性相對穩(wěn)定。故將應(yīng)用本文方法對圖10所示信號進行消噪處理。設(shè)置頻率窗口為0.8～1.3 Hz，并采用獲取圖5同樣的技術(shù)流程，可以得到消噪后的信號，其頻譜如圖12 所示。從圖12可知，實測原始信號中小于0.8 Hz、大于1.3 Hz的頻率成分被成功剔除，其0.8～1.3 Hz范圍內(nèi)有2個明顯的峰值，峰值處對應(yīng)的頻率分別約為0.85 Hz和0.95 Hz。為進一步驗證消噪信號的正確性，對其應(yīng)用ERA方法進行模態(tài)參數(shù)識別，其結(jié)果為0.916 19 Hz和1.036 7 Hz，相對1994年結(jié)果計算精度得到大幅提高。

圖10 加速度響應(yīng)信號

圖11 頻譜圖

圖12 消噪后頻譜圖

5 結(jié) 論

本文從信號分解的角度提出了一種基于極值-留數(shù)的信號分解消噪方法，解決了實際工程中在遭遇高背景噪聲其結(jié)構(gòu)固有信息往往被當(dāng)作環(huán)境噪聲而剔除的問題。選用的6自由度質(zhì)量-彈簧-阻尼系統(tǒng)已經(jīng)證實：當(dāng)信噪比(SNR)分別為40 dB和10 dB時，消噪后信號與不含噪聲信號的時程圖和頻譜圖幾乎重合，表明該方法能有效消除信號中強噪聲的干擾。計算獲得的不同信噪比下的噪聲抑制率(NSR)和信號失真率(SDR)，NSR約為0.8，SDR維持在較低水平，當(dāng)噪聲能量變化較大時，使用本文方法消噪得到的NSR變化較小，SDR較大。不同信噪比下，消噪后信號與不含噪聲信號的決定系數(shù)都大于0.995，表明消噪效果較好。為研究方法的有效性，選用了冰擊激勵作用下的某海洋石油平臺實測數(shù)據(jù)進行研究，獲得了預(yù)期頻率區(qū)間的頻率成分，基于此進行模態(tài)參數(shù)識別，可以大大提高模態(tài)參數(shù)識別精度，具有潛在的工程應(yīng)用價值。