基于CEEMD和小波包閾值的組合降噪及泄流結(jié)構(gòu)的模態(tài)識(shí)別方法

胡劍超， 練繼建， 馬 斌， 董霄峰

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

基于CEEMD和小波包閾值的組合降噪及泄流結(jié)構(gòu)的模態(tài)識(shí)別方法

胡劍超， 練繼建， 馬 斌， 董霄峰

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

水流荷載激勵(lì)下，泄流結(jié)構(gòu)模態(tài)識(shí)別的實(shí)測(cè)振動(dòng)信號(hào)常受噪聲干擾。以往EEMD算法有添加白噪聲造成非標(biāo)準(zhǔn)IMF導(dǎo)致的模態(tài)分裂問題及IMF有噪聲殘余不能完整重構(gòu)信號(hào)的完備性問題。而完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)法，通過在信號(hào)分解的每一層面添加特定高斯白噪聲，利用分解后第一階分量加總平均得到唯一余量計(jì)算固有模態(tài)函數(shù)，克服了EEMD的缺點(diǎn)；同時(shí)提出CEEMD與小波包閾值結(jié)合的組合降噪方法，運(yùn)用到向家壩水彈性模型實(shí)測(cè)振動(dòng)信號(hào)降噪中，驗(yàn)證了該組合方法降噪的有效性；為了提高帶噪振動(dòng)響應(yīng)模態(tài)識(shí)別的精度，基于數(shù)據(jù)相關(guān)技術(shù)，利用Markov參數(shù)構(gòu)造相關(guān)矩陣R，用該相關(guān)矩陣重構(gòu)Hankel矩陣后SVD分解得到系統(tǒng)最小實(shí)現(xiàn)，即數(shù)據(jù)相關(guān)特征系統(tǒng)實(shí)現(xiàn)法(ERA/DC)。將濾波降噪結(jié)合模態(tài)識(shí)別的整套方法，應(yīng)用到錦屏一級(jí)拱壩的泄流實(shí)測(cè)振動(dòng)響應(yīng)中，得到了較好的應(yīng)用效果。

泄流結(jié)構(gòu)；濾波降噪；完備總體經(jīng)驗(yàn)?zāi)B(tài)分解；小波包；模態(tài)識(shí)別；相關(guān)特征系統(tǒng)實(shí)現(xiàn)法

水工結(jié)構(gòu)的模態(tài)參數(shù)識(shí)別是確認(rèn)結(jié)構(gòu)運(yùn)行健康狀態(tài)的重要技術(shù)手段。由于數(shù)據(jù)采集環(huán)境復(fù)雜以及大型泄流結(jié)構(gòu)激勵(lì)難等狀況，必須消除或降低振動(dòng)響應(yīng)信號(hào)中的噪聲影響并在未知輸入激勵(lì)條件下進(jìn)行水工結(jié)構(gòu)的模態(tài)識(shí)別。

李成業(yè)等[1-3]研究的EMD/EEMD(Empirical Mode Decomposition With Ensemble EMD)+小波的濾波降噪方法具有諸多優(yōu)點(diǎn)：EEMD具有克服模態(tài)混疊、完全自適應(yīng)及IMFs都是單分量函數(shù)的特點(diǎn)，能彌補(bǔ)小波在處理非平穩(wěn)信號(hào)上的不足；同時(shí)小波變換具有多分辨率分析的優(yōu)點(diǎn)，利用小波閾值函數(shù)濾波，克服了EMD/EEMD濾除噪聲時(shí)，容易損失有用信號(hào)的缺點(diǎn)。

然而該組合方法仍有不足之處。EEMD通過原始信號(hào)中添加高斯白噪聲改進(jìn)EMD算法[4]。鄭近德等[5]指出這種方法使得分解后的信號(hào)包含殘余噪聲；每步構(gòu)造添加的白噪聲不同造成IMF分解數(shù)量不同，使得總體平均后的IMF不符合標(biāo)準(zhǔn)定義產(chǎn)生模態(tài)分裂；EEMD不能完整的重構(gòu)原始信號(hào)，不具完備性。基于多尺度分析的小波變換，通過尺度由大到小、時(shí)域-頻域變換的過程，由粗到精的展現(xiàn)信號(hào)的特征。然而由文獻(xiàn)[6]中小波算法可知，其只能對(duì)低頻小波細(xì)節(jié)子空間進(jìn)一步分解，而信號(hào)的高頻部分則沒有予以考慮。

由于方法自身的局限性，使得EEMD+小波閾值組合去噪方法在泄流結(jié)構(gòu)振動(dòng)響應(yīng)分析中受到一定的限制。為了克服EEMD的缺點(diǎn)，Torres等[7]提出完備總體模態(tài)經(jīng)驗(yàn)分解法(CEEMD)，而小波包變換[8](Wavelet Dacket Transform, WPT)既對(duì)信號(hào)低頻尺度又對(duì)高頻部分進(jìn)行精細(xì)分解。本文提出的基于CEEMD和小波包閾值的組合濾波降噪方法，不僅繼承了EEMD加小波方法的優(yōu)點(diǎn)，還改進(jìn)了該方法不足之處。文中通過向家壩水彈性模型實(shí)測(cè)振動(dòng)信號(hào)降噪，證明改進(jìn)組合降噪方法更加有效。

Lian等[9]運(yùn)用特征系統(tǒng)實(shí)現(xiàn)法(Eigensystem Realization Algolithm, ERA)，根據(jù)水工結(jié)構(gòu)工作狀態(tài)下只受環(huán)境激勵(lì)(如水流脈動(dòng)荷載)的特點(diǎn)，僅用拱壩泄流響應(yīng)成功進(jìn)行了模態(tài)識(shí)別。ERA實(shí)質(zhì)是利用脈沖或自由衰減響應(yīng)數(shù)據(jù)，構(gòu)造廣義Hankel矩陣進(jìn)行奇異值SVD分解，尋找系統(tǒng)的最小實(shí)現(xiàn)。Hankel矩陣的維數(shù)與模型階次相關(guān)聯(lián)，影響著系統(tǒng)矩陣A是否為系統(tǒng)最小實(shí)現(xiàn)，最終決定模態(tài)識(shí)別精度。若系統(tǒng)輸入響應(yīng)中混有噪聲時(shí)，對(duì)Hankel矩陣影響較大[10]。Juang[11-12]用相關(guān)響應(yīng)數(shù)據(jù)R代替脈沖響應(yīng)數(shù)據(jù)，構(gòu)造數(shù)據(jù)相關(guān)Hankel矩陣進(jìn)行SVD分解，即數(shù)據(jù)相關(guān)特征系統(tǒng)實(shí)現(xiàn)法(ERA/DC)。本文結(jié)合泄流激勵(lì)下水工結(jié)構(gòu)振動(dòng)特點(diǎn)，將ERA/DC引入泄流結(jié)構(gòu)模態(tài)參數(shù)識(shí)別中。并將基于組合降噪的ERA/DC應(yīng)用于錦屏一級(jí)拱壩實(shí)測(cè)振動(dòng)響應(yīng)分析中，得到了精確的識(shí)別結(jié)果。

1 完備總體模態(tài)經(jīng)驗(yàn)分解法CEEMD

以EEMD為基礎(chǔ)，CEEMD算法原理如下：

步驟1：x[n]為離散數(shù)n的原始信號(hào)；定義算子Ej(·)表示信號(hào)EMD分解后第j階IMF；wi[n](i=1,…,I)是單位方差零均值高斯白噪聲；εk系數(shù)是各階段信噪比；

(1)

步驟4：對(duì)于一階余量r1[n]中添加白噪聲信號(hào)，構(gòu)造信號(hào)r1[n]+ε1E1(wi[n])并進(jìn)行EMD分解，加總平均后得到第二階固有模態(tài)函數(shù)分量：

(2)

步驟6：對(duì)k階余量rk[n]中添加白噪聲信號(hào)，有構(gòu)造信號(hào)rk[n]+εkEk(wi[n])，進(jìn)行EMD分解，得到第k+1個(gè)固有模態(tài)函數(shù)分量：

(3)

步驟7：回到步驟(5)，計(jì)算下一個(gè)k。算法構(gòu)造循環(huán)取第一階IMFs分量加總平均，避免了EEMD同時(shí)分解的IMFs數(shù)量不同造成的模態(tài)分裂問題；

步驟8：重復(fù)步驟5～步驟7，直到余量不能再分解為止，最終得到：

(4)

式中：R是最終余量；K是固有模態(tài)分量總數(shù)。目標(biāo)信號(hào)x[n]可以表示為

(5)

式(5)為原始信號(hào)x[n]得到精確重構(gòu)，即該方法是具有完備性的分解方法。

下面結(jié)合一簡(jiǎn)單算例驗(yàn)證CEEMD方法的優(yōu)越性。設(shè)采樣頻率為200 Hz、采樣時(shí)間5 s的原始信號(hào)x，混入白噪聲k后得到含噪信號(hào)z，表達(dá)式如下：

z=x+k=3sin(πt)+5sin(2πt)+4.5sin(8πt)+5randn(1,1 000)

信號(hào)x，z如圖1，染噪信號(hào)z分別用EMD、EEMD、CEEMD方法濾波降噪(EEMD與CEEMD添加相同標(biāo)準(zhǔn)差與數(shù)量的白噪聲，一般取0.01倍～0.5倍原信號(hào)標(biāo)準(zhǔn)差，添加次數(shù)100～300次)，濾波后與x對(duì)比，如圖2所示。

(a)原始信號(hào)x(b)混入白噪聲后信號(hào)2

圖1 原始信號(hào)x與混入白噪聲后信號(hào)z

Fig.1 Original signalxand signalzmixed white noise

圖2 三種方法濾波降噪信號(hào)圖(局部)

從圖2(a)、圖2(b)可知，EMD或EEMD 方法濾波降噪都有不足之處，而圖(c)中CEEMD濾波降噪的重構(gòu)信號(hào)與原信號(hào)x波形最為接近。CEEMD方法改進(jìn)了EEMD不足之處，得到了顯著優(yōu)于EMD以及EEMD的濾波結(jié)果。

2 小波包(WPT)閾值濾波方法

小波包變換閾值濾波效果主要受閾值的影響。閾值函數(shù)可分為硬閾值函數(shù)和軟閾值函數(shù)。軟閾值降噪可實(shí)現(xiàn)最優(yōu)估計(jì)，保證降噪后的重構(gòu)信號(hào)與原信號(hào)有同樣的平滑，且不會(huì)產(chǎn)生附加震蕩，其表達(dá)式為

(6)

3 CEEMD與WPT閾值的組合濾波方法

壩體振動(dòng)主要受表中孔水流脈動(dòng)荷載激發(fā)，該脈動(dòng)荷載是具有一定頻帶寬度的平穩(wěn)隨機(jī)過程，表現(xiàn)為有色白噪聲譜的性質(zhì)。壩體動(dòng)力響應(yīng)通常是低頻或平穩(wěn)的振動(dòng)信號(hào)，而環(huán)境噪聲常為非平穩(wěn)的高頻信號(hào)。

CEEMD將實(shí)測(cè)振動(dòng)響應(yīng)信號(hào)的不同頻率模量由高到低自適應(yīng)清晰分離，分離出的帶有噪聲污染的信號(hào)固有模態(tài)分量，可利用小波包閾值的濾波方法有效去除噪聲。該方法適用于水工結(jié)構(gòu)的振動(dòng)信號(hào)濾波降噪處理，有利于下一步水工結(jié)構(gòu)模態(tài)精確識(shí)別。

基于CEEMD和WPT的組合濾波降噪方法步驟如下：

下面結(jié)合向家壩消力池內(nèi)導(dǎo)墻水彈性模型驗(yàn)證文中組合濾波方法的有效性，如圖3所示。

圖3 向家壩水彈模型隔振系統(tǒng)與右消力池內(nèi)導(dǎo)墻模型

Fig.3 Vibration isolation system of Xiang JiaBa hydroelastic model and Guide wall hydroelastic model in the right stilling basin

為了研究向家壩水電站泄洪誘發(fā)場(chǎng)地振動(dòng)規(guī)律，天津大學(xué)水利工程系建造了1∶80比尺的向家壩水彈性模型。為了減小外界振源對(duì)壩體振動(dòng)信號(hào)的干擾，采取在水彈模型混凝土基礎(chǔ)四周設(shè)置隔振溝、模型上下游連接段由止水橡膠軟連接等隔振措施，見圖3。試驗(yàn)證明，該減振隔振系統(tǒng)有效減小了隔振溝外振動(dòng)源的干擾[13]。

向家壩采用消力池消能，見圖3，右消力池內(nèi)設(shè)置一導(dǎo)墻水彈模型。導(dǎo)墻長(zhǎng)寬高為99 cm×4 cm×38 cm。材料模擬碾壓混凝土，彈模為105 MPa，泊松比0.167。布設(shè)4個(gè)動(dòng)位移測(cè)點(diǎn)。選擇測(cè)點(diǎn)4在上游水位370 m、表孔開度6 m工況實(shí)測(cè)導(dǎo)墻動(dòng)位移響應(yīng)數(shù)據(jù)進(jìn)行分析，采樣頻率200 Hz，取10 s數(shù)據(jù)分析。

圖4 實(shí)測(cè)導(dǎo)墻振動(dòng)信號(hào)與添加噪聲后振動(dòng)信號(hào)

由于設(shè)置減振隔振系統(tǒng)，該動(dòng)位移實(shí)測(cè)響應(yīng)數(shù)據(jù)x(n)的“毛刺”很少，信號(hào)較為光滑，可認(rèn)為其基本不含噪聲信號(hào)，如圖4(a)。在實(shí)測(cè)信號(hào)x(n)中加入白噪聲k，受噪聲污染信號(hào)z=x(n)+k(SNR=3.39 dB)，如圖4(b)所示，針對(duì)染噪信號(hào)z，分別用CEEMD、小波包(WPT)閾值、EEMD+小波、CEEMD+WPT四種方法進(jìn)行濾波降噪。

圖5為各方法濾波重構(gòu)波形與原信號(hào)x對(duì)比，表1為信噪比SNR、濾波信號(hào)與原始信號(hào)x互相關(guān)系數(shù)。

圖5 四種方法濾波降噪對(duì)比圖(局部)

信號(hào)類型降噪方法信噪比SNR/dB互相關(guān)系數(shù)構(gòu)造染噪信號(hào)z(SNR=3.39)CEEMD9.9070.9475小波包(WPT)11.2420.9626EEMD+小波11.8310.9666CEEMD+WPT13.0060.9747

由四種方法濾波降噪后的振動(dòng)信號(hào)與原始信號(hào)x(n)對(duì)比圖5看出：CEEMD分解加帶通濾波的做法，得到的SNR較低，重構(gòu)信號(hào)仍有較多毛刺，說明加入的噪聲有殘留且部分有用信號(hào)也被一并去除；小波包閾值濾波的信噪比大于CEEMD，但重構(gòu)波形仍存在毛刺；EEMD+小波的濾波方法，SNR較高，信號(hào)波形重構(gòu)較好、毛刺基本消除，說明噪聲基本被濾除，原始信號(hào)保留較好；而CEEMD+WPT的組合濾波降噪振動(dòng)信號(hào)，白噪聲“毛刺”不僅被很好的濾除干凈，從信噪比SNR=13.006 dB、相關(guān)系數(shù)=0.974 7來看，原始振動(dòng)信號(hào)x(n)保留最好、降噪效果最佳。

4 基于ERA/DC的泄流結(jié)構(gòu)模態(tài)參數(shù)識(shí)別方法

4.1數(shù)據(jù)相關(guān)特征系統(tǒng)實(shí)現(xiàn)法(ERA/DC)算法原理

根據(jù)振動(dòng)理論，線性系統(tǒng)的脈沖響應(yīng)函數(shù)矩陣，即馬可夫參數(shù)(Markov)為

Y0=D,Y1=CB,Y2=CAB,…,

Y(k)=CAk-1B

(7)

式中：Y(k)為l×m階矩陣，又是脈沖響應(yīng)函數(shù)矩陣。A為系統(tǒng)矩陣、B為輸入矩陣、C為輸出矩陣。由有限維線性時(shí)不變系統(tǒng)狀態(tài)方程知，系統(tǒng)模態(tài)識(shí)別問題即為系統(tǒng)的最小實(shí)現(xiàn)問題，即求得矩陣A，B，C并使其階次最小。由Markov參數(shù)Y(k)，構(gòu)造分塊Hankel矩陣：

(8)

式中：H(k-1)為lξ×mη階矩陣，ξ和η為任意整數(shù)。當(dāng)令k取1時(shí)，有Hankel矩陣H(0)。在ERA方法中，會(huì)針對(duì)H(0)進(jìn)行SVD分解。而在ERA/DC方法中，則構(gòu)建γ=lξ維方陣R(q)：

R(q)=H(q)HT(0)

(9)

式中：數(shù)據(jù)相關(guān)矩陣R(q)為lξ×lξ維，Hankel矩陣H(0)、H(q)為lξ×mη階矩陣。當(dāng)Hankel矩陣的列數(shù)很大時(shí)，R(q)的維數(shù)將?H(q)。若當(dāng)q=0時(shí)，有：

R(0)=H(0)HT(0)=

(10)

構(gòu)建(α+1)×(β+1)階塊相關(guān)矩陣U(q)：

U(q)=

(11)

式中：q為避免有噪聲存在時(shí)相關(guān)矩陣產(chǎn)生偏移的整數(shù)；r為防止相鄰R塊產(chǎn)生重疊的整數(shù)；α和β為分析中包含多少相關(guān)延遲。針對(duì)U(0)矩陣SVD分解：

U(0)=PDQT

(12)

式中：左奇異矩陣P和右奇異矩陣Q的維數(shù)為(α+1)γ×n和(β+1)γ×n，奇異值對(duì)角矩陣D為n×n階方陣。構(gòu)造U(0)偽逆矩陣U?，代入式(16)～式(18)推得：

A=D-1/2PTU(1)QD1/2

B=[EγTPD1/2]?H(0)EmC=EpTPD1/2

(13)

式中：EγT=[Iγ0γ]，EmT=[Im0m]，EpT=[Ip0p]，Iγ表示γ階單位矩陣，0γ表示γ階零矩陣。

針對(duì)系統(tǒng)矩陣A進(jìn)行奇異值分解，求得特征值矩陣Z=diag(Z1,Z2,…,Zn)，系統(tǒng)特征值與矩陣A特征值有關(guān)系：λi=ln(Zi)/T，T為采樣間隔。求得模態(tài)參數(shù)為

有阻尼頻率：ωdi=imag(λi)

阻尼比：ξi=-Re(λi)/|λi|

(14)

式中：Re(·)為取實(shí)部，imag(·)為取虛部，|·|表示取模。ERA/DC與ERA算法的主要區(qū)別在于首先進(jìn)行數(shù)據(jù)相關(guān)運(yùn)算后重構(gòu)Hankel矩陣U(0)，這樣可壓縮矩陣規(guī)模，減少奇異值分解的計(jì)算時(shí)間，并減少噪聲影響。相較于ERA，ERA/DC具有更好的識(shí)別精度。

4.2基于組合降噪的ERA/DC模態(tài)識(shí)別方法

ERA/DC通過Markov參數(shù)(即為脈沖響應(yīng)函數(shù))進(jìn)行模態(tài)參數(shù)識(shí)別。水工結(jié)構(gòu)在泄流激勵(lì)下振動(dòng)時(shí)，可假定為線性時(shí)不變系統(tǒng)，利用自然激勵(lì)技術(shù)(NExT)將參考點(diǎn)和響應(yīng)點(diǎn)數(shù)據(jù)變換為測(cè)點(diǎn)間的互相關(guān)函數(shù)，得到結(jié)構(gòu)的脈沖響應(yīng)函數(shù)。

在未知輸入激勵(lì)下，結(jié)構(gòu)的振動(dòng)模態(tài)階次通常是未知的。原始信號(hào)由延時(shí)嵌陷技術(shù)重構(gòu)的吸引子軌道矩陣D，與ERA方法中Hankel矩陣H(0)的構(gòu)造形式一致，基于奇異熵定階技術(shù)，H(0)的定階(即系統(tǒng)的定階)可根據(jù)奇異熵增量的變化進(jìn)行。觀察式11可知，ERA/DC用相關(guān)R矩陣構(gòu)造的Hankel矩陣，當(dāng)q=0時(shí)，U(0)的矩陣形式與H(0)相同，因此由ERA/DC進(jìn)行結(jié)構(gòu)模態(tài)識(shí)別，同樣可根據(jù)奇異熵增量譜確定系統(tǒng)階次。

為了進(jìn)一步減少實(shí)測(cè)信號(hào)中噪聲影響，模態(tài)識(shí)別之前須對(duì)響應(yīng)數(shù)據(jù)運(yùn)用CEEMD+WPT的組合濾波方法進(jìn)行降噪，這就是基于組合降噪方法的ERA/DC泄流結(jié)構(gòu)模態(tài)識(shí)別方法，其流程如圖6所示。

圖6 CEEMD加小波包組合濾波后的信號(hào)

為了說明ERA/DC方法針對(duì)帶噪信號(hào)識(shí)別的優(yōu)越性，建立混凝土懸臂梁有限元模型驗(yàn)證。梁長(zhǎng)20.0 m、截面尺寸1.4 m×0.6 m、動(dòng)彈模30 GPa。梁的一端全約束固定，如圖7。為了充分激發(fā)梁的各階模態(tài)，沿梁軸線選擇六個(gè)點(diǎn)P1～P6施加負(fù)Y向的典型水流脈動(dòng)荷載，荷載頻率100 Hz，時(shí)長(zhǎng)30 s，如圖8所示。

圖7 懸臂梁有限元模型與荷載施加方向圖

圖8 典型水流脈動(dòng)荷載時(shí)程線圖(局部)

計(jì)算懸臂梁動(dòng)力響應(yīng)，提取懸臂梁D1～D6測(cè)點(diǎn)振動(dòng)信號(hào)。以靠近懸臂梁固端的D6為參考點(diǎn)，其與梁中間D4測(cè)點(diǎn)的時(shí)程線如圖9，用NExT法計(jì)算D4-D6的脈沖響應(yīng)函數(shù)如圖10(a)所示。

得到脈沖響應(yīng)函數(shù)后計(jì)算奇異熵增量隨階次的變化如上圖10(b)所示，奇異熵階次為9階時(shí)，熵增量開始緩慢增長(zhǎng)并趨于平穩(wěn)，說明信號(hào)的有效特征信息量已經(jīng)趨于飽和，特征信息已經(jīng)基本完整。剔除特征值中的非模態(tài)項(xiàng)(非共軛根)和共軛項(xiàng)(重復(fù)項(xiàng))后，利用模態(tài)置信因子MAC準(zhǔn)則判別虛假模態(tài)[14]。得出結(jié)構(gòu)的前3階振動(dòng)階次進(jìn)行對(duì)比，如表2所示。

圖9 測(cè)點(diǎn)D4與D6振動(dòng)時(shí)程線圖(局部)

(a)脈沖響應(yīng)函數(shù)(b)脈沖響應(yīng)奇異熵增量

圖10 D4-D6點(diǎn)的脈沖響應(yīng)函數(shù)與奇異熵增量譜圖

由表2可知，ERA法識(shí)別頻率誤差<2%，而 ERA/DC方法頻率識(shí)別誤差<1%，都準(zhǔn)確識(shí)別出了懸臂梁前三階模態(tài)。

圖11 添加噪聲脈沖函數(shù)信號(hào)與濾波降噪信號(hào)

為了對(duì)比振動(dòng)信號(hào)受到噪聲污染情況下，ERA與ERA/DC方法的識(shí)別精度，將D4-D6脈沖響應(yīng)函數(shù)添加白噪聲(信噪比SNR=20 dB)，信號(hào)如圖11(a)所示。運(yùn)用兩種方法模態(tài)識(shí)別，結(jié)果如表3所示。

表3帶噪脈沖響應(yīng)信號(hào)的模態(tài)頻率識(shí)別對(duì)比(SNR=20)

Tab.3Modalfrequencyidentificationofpulse-responsesignalmixednoise(SNR=20dB)

階數(shù)數(shù)值模擬頻率ERA方法ERA/DC方法頻率/Hz誤差/%頻率/Hz誤差/%一階1.2121.2805.6111.2222.887二階2.8102.9946.5482.8261.281三階7.5418.1027.4397.5521.366

由表3看出，信號(hào)受噪聲污染后，ERA識(shí)別懸臂梁三階模態(tài)頻率，相較于無干擾情況誤差顯著增大。其一階、三階誤差都>5%，第三階誤差>7%。ERA/DC識(shí)別中，三階誤差都<3%，表明其模態(tài)識(shí)別抗噪聲干擾能力強(qiáng)于ERA方法。

將該染噪脈沖函數(shù)信號(hào)(D4點(diǎn))用CEEMD+小波包組合降噪，結(jié)果如圖11(b)所示。濾波后信號(hào)信噪比SNR=25.538、相關(guān)系數(shù)0.998 6。模態(tài)識(shí)別結(jié)果，如表4所示。

表4濾波脈沖函數(shù)響應(yīng)信號(hào)模態(tài)頻率識(shí)別對(duì)比

Tab.4Frequencyidentificationoffilteringpulse-responsesignal

階數(shù)數(shù)值模擬頻率ERA方法ERA/DC方法頻率/Hz誤差/%頻率/Hz誤差/%一階1.2121.2361.9801.2020.825二階2.8102.7472.2422.8441.210三階7.5417.7632.9447.5200.278

染噪信號(hào)濾波后識(shí)別結(jié)果由表4可知：ERA方法識(shí)別中，頻率誤差降至<3%，ERA/DC頻率識(shí)別誤差降至<1.5%，證明該組合降噪方法與ERA/DC方法結(jié)合的有效性與必要性。

5 錦屏一級(jí)拱壩實(shí)測(cè)振動(dòng)信號(hào)應(yīng)用分析

錦屏一級(jí)水電站是雅礱江干流下游河段的控制性梯級(jí)電站。其最大壩高305 m，是目前世界最高的混凝土雙曲拱壩。2014年天津大學(xué)水利工程系針對(duì)錦屏一級(jí)樞紐工程進(jìn)行了水力學(xué)原型觀測(cè)。沿壩頂1#～24#壩段，水平向順?biāo)鞑荚O(shè)H1～H7共7個(gè)壩頂動(dòng)位移傳感器。水庫汛期泄洪期間，實(shí)測(cè)了不同泄流工況下，泄流激勵(lì)下壩體振動(dòng)響應(yīng)。傳感器測(cè)點(diǎn)布置圖及錦屏一級(jí)拱壩有限元模型，如圖12所示。

圖12 傳感器測(cè)點(diǎn)布置及錦屏一級(jí)拱壩有限元模型圖

5.1濾波降噪

流激振動(dòng)造成的拱壩動(dòng)位移均方根通常為幾十個(gè)微米。由于振動(dòng)量級(jí)較小，混入的干擾信號(hào)和背景噪聲容易淹沒有用信號(hào)。分析各工況壩體振動(dòng)規(guī)律：拱冠梁處振動(dòng)最大但遠(yuǎn)離干擾源，壩肩振動(dòng)量小且易受噪聲干擾。因此選取工況7(3#表孔25%開度、泄流量329 m3/s、上游水位1 879.2 m、下游1 645.6 m)時(shí)，拱冠梁附近參考測(cè)點(diǎn)H3及右壩肩附近測(cè)點(diǎn)H6，采樣頻率200 Hz，時(shí)長(zhǎng)40 s，其時(shí)程線及歸一化功率譜，如圖13所示。

圖13 H3、H6測(cè)點(diǎn)振動(dòng)信號(hào)(局部)及歸一化功率譜圖

由錦屏有限元模態(tài)計(jì)算得拱壩前六階模態(tài)頻率，頻率范圍在1.0～3.5 Hz之間。測(cè)點(diǎn)H3位于拱冠梁附近，受到噪聲干擾小，由其功率譜圖13(b)可知其振動(dòng)能量集中在1.0～3.0 Hz頻帶部分，與有限元計(jì)算相符。

在現(xiàn)場(chǎng)泄流觀測(cè)期間，發(fā)現(xiàn)右壩肩附近存在少量施工作業(yè)。壩肩測(cè)點(diǎn)H6振動(dòng)信號(hào)“毛刺”較多，并出現(xiàn)漂移現(xiàn)象，運(yùn)用最小二乘擬合法消除趨勢(shì)項(xiàng)后如圖13(c)。功率譜如圖13(d)，<0.5 Hz的低頻部分有能量集中，且>4.5 Hz處有寬頻白噪聲存在，直接應(yīng)用H6實(shí)測(cè)振動(dòng)數(shù)據(jù)進(jìn)行模態(tài)識(shí)別會(huì)產(chǎn)生很大誤差。

為了得到拱壩真實(shí)振動(dòng)信號(hào)，運(yùn)用文中提出的組合濾波降噪方法濾除原始信號(hào)中的噪聲。將H6測(cè)點(diǎn)信號(hào)CEEMD分解，添加白噪聲標(biāo)準(zhǔn)差1.494 μm，添加數(shù)300次，最大迭代數(shù)5 000次。分解得IMF1～I(xiàn)MF13分量，IMFs分量從高頻到低頻依次排列如上圖14。

圖14 H6測(cè)點(diǎn)實(shí)測(cè)信號(hào)CEEMD分解圖

做出IMF1～I(xiàn)MF12分量歸一化功率譜圖如下圖15所示。前兩階IMF能量集中在高頻且頻帶很寬，為寬頻白噪聲分量，運(yùn)用低通濾波進(jìn)行濾除。IMF3～I(xiàn)MF6是振動(dòng)信號(hào)分量，但其中混入了大量噪聲信號(hào)，需要運(yùn)用WPT方法進(jìn)一步濾除噪聲。IMF7～I(xiàn)MF9分量振動(dòng)能量集中在1.0～3.0 Hz附近，判斷其屬于結(jié)構(gòu)真實(shí)振動(dòng)信號(hào)。而IMF10～13為低頻趨勢(shì)項(xiàng)干擾，予以高通濾除。

IMF3～I(xiàn)MF6混有白噪聲，運(yùn)用小波包閾值方法進(jìn)一步降噪。小波包分解采用dbN小波基，共N層分解，不同IMF分量采用不同的N值。采用軟閾值方法，默認(rèn)全局閾值。濾波后，由IMF3～I(xiàn)MF9重構(gòu)信號(hào)如圖16所示。

從圖16可知，原始信號(hào)毛刺被去除。從功率譜曲線可知，低頻干擾噪聲和高頻白噪聲被有效消除了，響應(yīng)中反映結(jié)構(gòu)模態(tài)信息的共振分量所占的比重大大增強(qiáng)。

圖15 H6信號(hào)IMF分量歸一化功率譜圖

圖16 H6組合濾波降噪振動(dòng)信號(hào)(局部)及功率譜圖

5.2模態(tài)識(shí)別

由圖17可知，H6測(cè)點(diǎn)實(shí)測(cè)振動(dòng)響應(yīng)濾除噪聲信號(hào)后，以拱冠梁附近的測(cè)點(diǎn)H3為參照點(diǎn)，運(yùn)用NExT方法計(jì)算H3與H6測(cè)點(diǎn)脈沖響應(yīng)函數(shù)；由H6-H3脈沖響應(yīng)函數(shù)，首先構(gòu)造4 000×4 000階Hankel矩陣H(0)以及2 000×2 000階Hankel矩陣U(0)，計(jì)算奇異熵增量隨階次變化。確定奇異熵階次為11階，采用文中ERA/DC方法進(jìn)行識(shí)別。剔除非模態(tài)項(xiàng)和共軛項(xiàng)后用MAC準(zhǔn)則判別虛假模態(tài)，得到結(jié)構(gòu)前6階振動(dòng)階次，識(shí)別結(jié)果列于表5。表5還列出了拱壩有限元模態(tài)計(jì)算結(jié)果，并將該振動(dòng)信號(hào)同時(shí)運(yùn)用ITD、STD、Prony(復(fù)指數(shù)法)、ARMA(時(shí)序分析法)、SSI(隨機(jī)子空間法)與ERA(特征系統(tǒng)實(shí)現(xiàn)法)進(jìn)行模態(tài)識(shí)別，識(shí)別結(jié)果一同列于表5中。

圖17 H3-H6測(cè)點(diǎn)脈沖響應(yīng)函數(shù)與奇異熵增量譜

表5 錦屏一級(jí)拱壩模態(tài)頻率識(shí)別對(duì)比

拱壩是一種典型的密頻結(jié)構(gòu)，綜合各方法的模態(tài)識(shí)別的結(jié)果，可得到錦屏一級(jí)拱壩的前六階工作頻率。可以看出，ERA/DC方法對(duì)于低頻、密頻具有較強(qiáng)的識(shí)別能力。

將各方法識(shí)別結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)識(shí)別結(jié)果都較為接近，都能識(shí)別出拱壩前六階頻率。而文中所提方法與拱壩有限元計(jì)算結(jié)果同樣相吻合。文中方法的計(jì)算結(jié)果不僅能夠滿足工程分析的需要，且具有相對(duì)更好的計(jì)算精度。

6 結(jié) 論

本文針對(duì)泄流結(jié)構(gòu)振動(dòng)響應(yīng)信號(hào)，提出一種組合濾波降噪結(jié)合模態(tài)識(shí)別的整套方法，應(yīng)用到實(shí)際泄流結(jié)構(gòu)的狀態(tài)參數(shù)識(shí)別中，得到以下結(jié)論：

(1) CEEMD解決了EEMD算法導(dǎo)致的模態(tài)分裂及完備性問題，小波包能獲得豐富的時(shí)頻局部化信息，更適合突變信號(hào)、非平穩(wěn)信號(hào)的分析與濾波。向家壩水彈性模型實(shí)測(cè)振動(dòng)信號(hào)降噪說明，CEEMD+小波包閾值的組合降噪方法較EEMD+小波方法具有更好的濾波降噪效果。

(2) 基于數(shù)據(jù)相關(guān)技術(shù)的ERA/DC方法克服了泄流結(jié)構(gòu)輸入激勵(lì)條件難獲得的問題，利用懸臂梁振動(dòng)模型進(jìn)行驗(yàn)證，表明其對(duì)含噪信號(hào)具有比ERA更好的模態(tài)識(shí)別精度，同時(shí)說明了其與組合降噪方法結(jié)合的必要性。

(3) 將CEEMD+小波包閾值組合降噪方法結(jié)合ERA/DC模態(tài)識(shí)別方法，應(yīng)用到錦屏一級(jí)拱壩泄流實(shí)測(cè)振動(dòng)響應(yīng)信號(hào)中，得到了很好的實(shí)際應(yīng)用效果。

[1] 李成業(yè)，練繼建，劉昉，等. EMD與小波閾值聯(lián)合濾波方法的改進(jìn)及其在泄流結(jié)構(gòu)振動(dòng)分析中的應(yīng)用[J]. 振動(dòng)與沖擊，2013，32(9)：63-70.

LI Chengye, LIAN Jijian, LIU Fang, et al. An improved filtering method based on EMD and wavelet-threshold and its application in vibration analysis for a flood discharge structure[J]. Journal of vibration and shock, 2013, 32(9): 63-70.

[2] ZHANG Yan, LIAN Jijian, LIU Fang. An improved filtering method based on EEMD and wavelet-threshold for modal parameter identification of hydraulic structure[J]. Mechanical Systems and Signal Processing, 2016, 68: 316-329.

[3] 邵忍平，曹精明，李永龍. 基于EMD小波閾值去噪和時(shí)頻分析的齒輪故障模式識(shí)別與診斷[J]. 振動(dòng)與沖擊，2012，31(8)：96-101.

SHAO Renping, CAO Jingming, LI Yonglong. Gear fault pattern identification and diagnosis using time-frequency analysis and wavelet threshold de-noising based on EMD[J]. Journal of vibration and shock, 2012, 31(8): 96-101.

[4] WU Z H, HUANG N E. Ensemble empirical mode decomposition: A noise assisted data analysis method[J]. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41.

[5] 鄭近德，程軍圣，楊宇. 改進(jìn)的EEMD算法及其應(yīng)用研究[J]. 振動(dòng)與沖擊，2013，32(21)：21-26.

ZHENG Jinde, CHENG Junsheng, YANG Yu. Modified EEMD algorithm and its applications[J]. Journal of Vibration and Shock, 2013, 32(21): 21-26.

[6] OUADFEUL S A, ALIOUANE L. Random seismic noise attenuation data using the discrete and continuous wavelet transforms[J]. Arabian Journal of Geosciences, 2014, 7(7): 2531-2537.

[7] TORRES M E, COLOMINAS M A, SCHLOTTHAUER G, et al.A complete ensemble mode decomposition with adaptive noise[C] // IEEE International Conference on Acoustics. Speech and signal Processing(ICASSP). USA, NJ: IEEE. 2011. 4144-4147.

[8] HEMMATI F, ORFALI W, GADALA M S. Roller bearing acoustic signature extraction by wavelet packet transform, applications in fault detection and size estimation[J]. Applied Acoustics, 2016, 104(1): 101-118.

[9] LIAN J J,LI H K,ZHANG J W. ERA modal identification method for hydraulic structures based on order determination and noise reduction of singular entropy[J]. Science in China series E-Technological Sciences, 2009, 52(2): 400-412.

[10] 李蕾紅，陸秋海，任革學(xué).特征系統(tǒng)實(shí)現(xiàn)法的識(shí)別特性研究及算法的推廣[J]. 工程力學(xué)，2002，19(1)：109-115.

LI Leihong, LU Qiuhai, REN Gexue. A study of eigensystem realization algorithm and its generalization[J]. Engineering Mechanics, 2002, 19(1): 109-115.

[11] JUANG J N, COOPER J E, WRIGHT J R. An eigensystem realization algorithm using data correlation (ERA/DC) for modal parameter identification[J]. Control Theory and Advanced Technology, 1988, 4(1): 5-14.

[12] CHIANG D Y, LIN C S. Identification of modal parameters from ambient vibration data using eigensystem realization algorithm with correlation technique[J]. Journal of Mechanical Science and Technology, 2010, 24(12): 2377-2382.

[13] 李淑君，練繼建，歐陽秋平. 高壩泄流誘發(fā)壩區(qū)及場(chǎng)地振動(dòng)的振源分析研究[J]. 水利水電技術(shù)，2014，45(9)：47-51.

LI Shujun, LIAN Jijian, OUYANG Qiupin. Study on source of dam discharge-induced vibration in dam region and site[J]. Water Resources and Hydropower Engineering, 2014, 45(9): 47-51.

[14] 董霄峰，練繼建，楊敏，等. 海上風(fēng)機(jī)結(jié)構(gòu)工作模態(tài)識(shí)別的組合降噪方法[J]. 天津大學(xué)學(xué)報(bào)(自然科學(xué)與工程技術(shù)版)，2015，48(3)：203-208.

DONG Xiaofeng, LIAN Jijian, YANG Min, et al. A compound de-noising method for operational modal identification of offshore wind turbine structure[J]. Journal of Tianjin University(Science and Technology),2015,48(3):203-208.

Ade-noisingandmodalidentificationcombinedmethodbasedonCEEMDandwaveletpacketthresholdforflooddischargestructures

HU Jianchao， LIAN Jijian， MA Bin， DONG Xiaofeng

(State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin, 300072 China)

Vibration signals measured for modal identification of a flood discharge structure under water load excitation are often interfered by noise. The EEMD algorithm has problems of modes splitting caused by nonstandard IMFs produced due to adding white noise and IMFs with residual noise not being able to completely reconstruct original signals. Here, the complete ensemble empirical mode decomposition (CEEMD) was proposed. Using it, the drawbacks of EEMD were overcome through adding a specific white noise to each stage of signal decomposition and using the unique residual obtained by adding the average to the first order component gained after decomposition to calculate IMFs. Furthermore, a de-noising method with CEEMD and the wavelet packet threshold combined was proposed to be applied in de-noising vibration signals measured from Xiangjia dam’s hydro-elastic model, the effectiveness of this method was verified. In order to improve the precision of modal identification with the noised vibration responses, based on the data correlation technology, Markov parameters were used to construct a correlation matrix R. The R matrix was employed to reconstruct a Hankel matrix. After its SVD, the minimum system realization was obtained. This algorithm was called the eigen-system realization algorithm with data correlation (ERA/DC). Finally, the proposed de-noising and modal identification combined method was applied in the measured vibration responses of Jinping arch dam to get a better application effect.

flood discharge structure; noise reduction; CEEMD; wavelet packet; modal identification; ERA/DC

國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)子課題“極端條件下高碾壓混凝土壩動(dòng)力響應(yīng)機(jī)制”(2013CB035905-1)

2016-06-06 修改稿收到日期：2016-08-08

胡劍超 男，碩士生，1990年生

練繼建 男，教授，博士生導(dǎo)師，1965年生

TV65

： A

10.13465/j.cnki.jvs.2017.17.001