劉景良 ， 鄭錦仰， 林友勤， 邱仁輝， 駱勇鵬

( 1. 福建農(nóng)林大學(xué) 交通與土木工程學(xué)院，福州 350002； 2. 福州大學(xué) 土木工程學(xué)院，福州 350108)

環(huán)境激勵下的土木工程結(jié)構(gòu)/構(gòu)件本質(zhì)上屬于時變和非線性系統(tǒng)，其響應(yīng)信號不但呈現(xiàn)非平穩(wěn)性，同時也具有模態(tài)密集性。通常情況下，為準(zhǔn)確識別時變結(jié)構(gòu)的特征參數(shù)需對多分量響應(yīng)信號進(jìn)行有效分解。目前，常見的信號分解方法有經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)[1]、局部均值分解(Local Mean Decomposition，LMD)[2]和集合經(jīng)驗?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)[3]等遞歸模式分解方法。其中，EMD與LMD都是根據(jù)信號時間尺度的局部特征對分量信號的包絡(luò)線進(jìn)行估計，但在估計過程中誤差不斷傳播，分解時有可能出現(xiàn)模態(tài)疊混現(xiàn)象，故無法有效分解頻率相近的信號分量。不僅如此，噪聲的干擾有可能使得原始信號的非極值點變?yōu)闃O值點，因而信號的包絡(luò)線估計會出現(xiàn)誤差[4]。EEMD通過加入白噪聲對EMD中出現(xiàn)的模態(tài)疊混現(xiàn)象進(jìn)行抑制，但該方法需要進(jìn)行成百上千次運算，且有可能分解出超出信號真實成分的分量個數(shù)[3]。此外，遞歸模式分解方法還存在端點效應(yīng)等問題[5]，這給實際應(yīng)用帶來了不少困難。

最近，Dragomiretskiy等[6]提出一種非遞歸模式分解方法，即變分模態(tài)分解(Variational Mode Decomposition，VMD)。該方法通過迭代搜尋變分模型的最優(yōu)解并以此確定所有分量的頻率中心和帶寬，從而直接從頻域中分離各分量信號。截至目前，VMD已在參數(shù)識別和損傷識別等領(lǐng)域得到廣泛應(yīng)用。如Wang等[7]利用VMD成功分解了結(jié)構(gòu)齒輪中的摩擦信號并準(zhǔn)確識別了齒輪故障位置。張蒙等[8]將基于VMD和EMD的多尺度排列熵進(jìn)行比較，結(jié)果表明VMD算法不僅能夠克服EMD算法不能分離模態(tài)疊混信號的缺點，且分解出的分量信號在各個尺度均包含更多細(xì)節(jié)。Chen等[9]基于VMD算法和信號調(diào)制方法提出能夠分解非窄帶信號的VVMD算法，進(jìn)一步拓寬了VMD算法的使用范圍。雖然VMD不具有遞歸模式分解方法中的一些缺點，但其算法本身也存在一般頻域分割算法無法避免的共同問題，即需要預(yù)先判斷原信號中分量信號的個數(shù)[6]。分量信號數(shù)目的判別可通過頻率相近原則[10-11]、粒子群優(yōu)化算法[12]或是本征信號的峭度值[13]等方法來實現(xiàn)，但上述方法均是通過對分量個數(shù)進(jìn)行試取從而得到最優(yōu)值，降低了算法的整體計算效率?？偠灾?，VMD作為一種新的多分量信號分解方法，其算法的改進(jìn)和在土木工程中的應(yīng)用還十分缺乏。

VMD分解得到的分量信號通常是調(diào)幅調(diào)頻信號，需要有效的時頻分析方法進(jìn)行分析處理。Daubechies等[14]提出的同步擠壓小波變換(Synchrosqueezing Wavelet Transform，SWT)將小波變換后的時頻圖進(jìn)行重組從而獲得較高頻率精度的時頻曲線。劉景良等[15]將同步擠壓小波變換引入土木工程領(lǐng)域，并研究了不同小波母函數(shù)對同步擠壓小波變換的影響。汪祥莉等[16]利用同步擠壓小波變換對混沌背景下的諧波信號進(jìn)行提取，驗證了同步擠壓小波變換具有較強(qiáng)的魯棒性。然而，同步擠壓小波變換并不能有效分離瞬時頻率相近的分量信號，且重構(gòu)的分量信號與理論值相差較遠(yuǎn)[17-18]。

基于此，本文將變分模態(tài)分解方法引入土木工程領(lǐng)域，并結(jié)合同步擠壓小波變換理論，提出一種識別時變結(jié)構(gòu)瞬時頻率的新方法。首先，對多分量信號進(jìn)行連續(xù)小波變換得到相應(yīng)的小波量圖，然后根據(jù)圖中的時頻分布情況準(zhǔn)確判斷分量信號的個數(shù)，從而大大提高了計算效率。其次，利用VMD對多分量信號進(jìn)行分解，彌補(bǔ)了SWT重構(gòu)單分量信號精度不足的問題。最后，采用同步擠壓小波變換識別分量信號的瞬時頻率，從而提高了其頻率精度。通過一個多分量信號數(shù)值算例、一個時變懸臂梁試驗以及一個時變拉索試驗對所提出的方法進(jìn)行驗證，并與基于希爾伯特-黃變換(Hilbert-Huang Transform，HHT)[19]和連續(xù)小波變換(Continuous Wavelet Transform，CWT)[20]的瞬時頻率識別結(jié)果進(jìn)行對比，從而驗證了該方法的有效性與準(zhǔn)確性。

1 基本理論

1.1 變分模態(tài)分解

VMD作為一種新型自適應(yīng)模態(tài)分解方法，其主要目的是將輸入信號分解為一系列具有稀疏特性的分量信號uk(t)。一般地，多分量信號x(t)可通過式(1)來表示。

x(t)=u1(t)+u2(t)+…+uk(t)

(1)

式中：分量信號uk(t)為具有中心頻率和一定帶寬的調(diào)幅調(diào)頻信號。

VMD分解方法主要是通過求解模態(tài)分量的變分問題來確定各分量信號的帶寬和中心頻率。在構(gòu)造變分問題時，其基本原理是對響應(yīng)信號使用希爾伯特變換、頻率混合等信號處理手段，在各階模態(tài)分量之和等于原信號的約束前提下，將有關(guān)模態(tài)分量的變分問題轉(zhuǎn)化為尋求估計帶寬之和最小的模態(tài)函數(shù)。其具體步驟如下：

首先對分量信號uk(t)進(jìn)行希爾伯特變換得到對應(yīng)的單邊頻譜；其次將得到的單邊頻譜與e(-iωk t)相乘，從而使每個分量信號的頻譜調(diào)整至以預(yù)估中心頻率ωk為中心的頻帶；最后計算頻率混合后信號梯度范數(shù)的平方并估計移頻后分量信號的帶寬，得到如式(2)所示的約束優(yōu)化問題。

(2)

為求解式(2)中的約束優(yōu)化問題，引入二次罰函數(shù)因子α以保證噪聲干擾下的信號重構(gòu)精度；同時，引入拉格朗日乘法算子λ以保證約束條件的嚴(yán)格性。最后采用拓展拉格朗日表達(dá)式L表示如式(3)所示的無約束優(yōu)化問題：

L({uk},{ωk},λ)=

(3)

(1) 初始化{u1k}，{ω1k}，λ1，n，并賦初始值為零。

(4)

式中：α建議取值為2 000。

(5)

(4) 根據(jù)式(6)迭代更新λn+1。

(6)

(5) 循環(huán)(2)～(5)步直至滿足式(7)所示的迭代收斂條件

(7)

1.2 同步擠壓小波變換

對于給定的母小波函數(shù)ψ(t)，對分量信號uk(t)進(jìn)行如式(8)所示的連續(xù)小波變換，其小波變換系數(shù)Wu(a,b)為

(8)

(9)

(10)

1.3 基于變分模態(tài)分解和同步擠壓小波變換識別時變結(jié)構(gòu)瞬時頻率

首先，基于變分模態(tài)分解和同步擠壓小波變換識別時變結(jié)構(gòu)瞬時頻率的方法根據(jù)小波量圖中的時頻曲線分布情況確定模態(tài)分量個數(shù)。其次，通過VMD將響應(yīng)信號分解成一系列單分量信號，在避免模態(tài)疊混現(xiàn)象的同時，也有利于后續(xù)的SWT算法。最后，采用SWT識別單分量信號的瞬時頻率。該算法的具體流程，如圖1所示。

圖1 基于變分模態(tài)分解和同步擠壓小波變換識別時變結(jié)構(gòu)瞬時頻率流程圖Fig.1 The flow chart of the instantaneousfrequency identification of time-varying structures based on VMD and SWT

2 數(shù)值算例驗證

考慮如下多分量信號x(t)：

x(t)=x1(t)+x2(t)=cos[2πt+sin(0.5πt)]+
[2+cos(2πt)]cos[4πt+ sin(πt)]

(11)

式中：信號采樣頻率為100 Hz，采樣時間為10 s。分量信號對應(yīng)的瞬時頻率理論值分別為f1=1+0.25cos(0.5πt) Hz和f2=2+0.5cos(πt) Hz。為考慮實際噪聲影響，對信號添加20%水平的高斯白噪聲，噪聲強(qiáng)度由信噪比定義。

(12)

添加20%噪聲水平后的信號如圖2所示。首先，對含噪信號進(jìn)行連續(xù)小波變換得到如圖3所示的小波量圖。由圖3可知，分量信號的瞬時頻率主要集中在[0.5 Hz，1.5 Hz]和[1.5 Hz，4 Hz]兩個頻率區(qū)間內(nèi)，因此可判斷分量信號個數(shù)為2。其次對信號x(t)進(jìn)行VMD分解，得到的分量信號如圖4所示。為證明VMD的有效性，采用EMD對原信號x(t)進(jìn)行分解直至殘差函數(shù)中的零點個數(shù)與極值點個數(shù)的差大于1，提取的前2階本征分量如圖4所示。由圖4可知，由于噪聲的影響，EMD分解出的分量信號與理論值有較大偏差，且端點效應(yīng)明顯大于VMD分解方法。而VMD分解得到的分量信號與理論值更加吻合，且沒有出現(xiàn)明顯的端點效應(yīng)。最后，對VMD分解后的分量信號進(jìn)行SWT并識別分量信號的瞬時頻率，如圖5所示。圖5中同時給出了HHT和CWT兩種方法的瞬時頻率識別結(jié)果。

圖2 添加20%水平噪聲的多分量信號Fig.2 The simulated multi-component signal with 20% Gaussian white noise

圖3 多分量信號的小波量圖Fig.3 The wavelet scalogram of the simulated multi-component signal

(a) x1

(b) x2圖4 提取的單分量信號Fig.4 The extracted mono-components

(a) x1

(b) x2圖5 瞬時頻率識別結(jié)果Fig.5 The instantaneous frequency extracted by the proposed method and HHT

從圖5可以看出，基于變分模態(tài)分解和同步擠壓小波變換方法(VMD+SWT)識別的瞬時頻率值與理論值更加吻合，且識別效果明顯優(yōu)于HHT和CWT。

為量化瞬時頻率的識別精度，采用瞬時頻率在整個時間歷程內(nèi)的均方根值作為精度指標(biāo)(Index of Accuracy, IA)：

(13)

式中：fd(t)為瞬時頻率識別值；fe(t)為瞬時頻率理論值。精度指標(biāo)IA值越小，說明識別值與理論值越接近，即精度越高。

表1給出了VMD+SWT、HHT和CWT三種方法所對應(yīng)的IA值。其中，IA1和IA2分別表示分量信號x1與x2的瞬時頻率識別精度指標(biāo)。由表1可知，針對多分量信號x(t)，VMD+SWT的瞬時頻率識別效果明顯優(yōu)于HHT和CWT兩種方法。

表1 添加20%高斯白噪聲多分量信號的瞬時頻率識別精度指標(biāo)Tab.1 IA of instantaneous frequency identification of the multi-componentsignal with 20% Gaussian white noise %

3 試驗驗證

為驗證變分模態(tài)分解和同步擠壓小波變換聯(lián)合方法識別時變結(jié)構(gòu)瞬時頻率的有效性和準(zhǔn)確性，設(shè)計一個質(zhì)量突變懸臂梁試驗與一個剛度時變拉索試驗對其進(jìn)行驗證。

3.1 質(zhì)量突變懸臂梁試驗

試驗中，鋁合金懸臂梁的截面尺寸為40 mm×15 mm，長度為500 mm，重量為0.81 kg，懸臂端的質(zhì)量塊重1 kg。質(zhì)量塊正上方為一塊方形磁鐵。整個試驗裝置如圖6所示。試驗前，測得帶重塊與不帶重塊時的懸臂梁的基頻分別為21.24 Hz和47.06 Hz。試驗過程中，先用力錘在懸臂梁末端施加激振力，并在2 s時利用磁鐵吸引質(zhì)量塊使之與懸臂梁分離，從而讓懸臂梁的質(zhì)量產(chǎn)生變化。在質(zhì)量塊與懸臂梁分離后，在2.3 s時再次對懸臂梁進(jìn)行力錘激勵。通過放置在懸臂梁跨中處的加速度傳感器(靈敏度：98.7 mV/m/s2)記錄懸臂梁在整個時間歷程的響應(yīng)信號，結(jié)果如圖7所示。

在VMD分解前，對響應(yīng)信號進(jìn)行連續(xù)小波變換得到如圖8所示的小波量圖。由圖8(a)可知，響應(yīng)信號的瞬時頻率主要集中在[20 Hz，50 Hz]和[250 Hz,280 Hz]頻率區(qū)間內(nèi),從而判斷分量信號的個數(shù)為2。而圖8(b)又表明了[20 Hz，50 Hz]頻率區(qū)間內(nèi)的瞬時頻率在2s附近發(fā)生突變，這與實際試驗情況相吻合。對原始響應(yīng)信號進(jìn)行VMD分解，得到的一階分量信號如圖9所示。最后，對分解后的分量信號進(jìn)行同步擠壓小波變換以識別懸臂梁結(jié)構(gòu)的瞬時頻率，結(jié)果如圖10所示。與HHT和CWT識別的瞬時頻率結(jié)果相比，采用VMD+SWT得到的瞬時頻率更接近理論值，且沒有明顯的端點效應(yīng)。表2中給出了三種方法對應(yīng)的瞬時頻率識別精度指標(biāo)IA。從表2可以看出，HHT和CWT的識別精度與VMD+SWT相比分別相差了9.1%和18.1%，這也再次證明了VMD+SWT方法的有效性與準(zhǔn)確性。

圖6 鋁合金懸臂梁試驗裝置圖(mm)Fig.6 The aluminum cantilever beam test rig(mm)

圖7 實測響應(yīng)信號Fig.7 The measured acceleration responses

(a) [20 Hz,50 Hz]和[250 Hz,280 Hz]頻率區(qū)間

(b) [20Hz,50Hz]頻率區(qū)間圖8 小波量圖Fig.8 The Wavelet scalogram

圖9 VMD和EMD分解后的分量信號Fig.9 The mono-component signal decomposed by VMD and EMD

圖10 質(zhì)量突變懸臂梁瞬時頻率識別值Fig.10 The identified instantaneous frequency of the cantilever beam with abrupt mass reduction

表2 質(zhì)量突變懸臂梁瞬時頻率識別精度指標(biāo)Tab.2 The IA of instantaneous frequency identification of the cantilever beam with abrupt mass reduction %

3.2 剛度時變拉索試驗

試驗所用拉索為一根7φs5的鋼絞線，一端用反力架錨固，另一端固定于電液伺服加載系統(tǒng)(MTS)加載系統(tǒng)。兩錨固點間的索長為4.55 m，將加速度傳感器豎向安裝在拉索中部。在使用MTS的作動器對索施加20 kN的預(yù)拉力后，通過調(diào)整作動器使索的拉力產(chǎn)生變化，從而模擬索的剛度隨時間變化的情況。改變索力的同時，用力錘敲擊拉索，采集索的豎向加速度沖擊響應(yīng)，采樣頻率為600 Hz，采樣時間為6 s。試驗裝置如圖11所示。

試驗過程中模擬了索拉力線性變化和正弦變化兩個工況。本文采用“凍結(jié)法[21]”，即假定在很小的時間間隔內(nèi)結(jié)構(gòu)參數(shù)保持不變，通過求解系統(tǒng)振動方程的特征值和特征向量近似得到拉索瞬時頻率的理論值。

首先，進(jìn)行拉力線性變化時索的瞬時頻率識別，試驗時拉力從20 kN開始以1.67 kN/s的速率線性增加。測得索的加速度響應(yīng)如圖12所示。同理，先對響應(yīng)信號進(jìn)行連續(xù)小波變換得到信號的小波量圖，如圖13所示。從圖13中可以看出，整個時頻面上共有兩個明顯的分量信號。因此，確定VMD中需要選取的模態(tài)分量個數(shù)為2。利用VMD得到如圖14所示的一階分量信號，然后對其進(jìn)行同步擠壓小波變換，得到的瞬時頻率識別結(jié)果如圖15所示。在圖15和表3中同時給出了采用VMD+SWT、HHT和CWT三種方法提取的瞬時頻率曲線以及對應(yīng)的精度指標(biāo)IA。由此可知：VMD+SWT方法具有較高的精度，且識別效果優(yōu)于HHT方法和CWT。值得注意的是，雖然CWT提取的瞬時頻率識別精度較高，其識別的瞬時頻率曲線實際上存在大量毛刺。

圖11 試驗裝置圖Fig.11 The cable testsetup

圖12 拉力線性變化時索的加速度響應(yīng)Fig.12 The measured cable acceleration responses with linearly varying cable tension force

圖13 拉力線性變化時實測信號小波量圖Fig.13 The wavelet scalogram of the measured cable acceleration responses with linearly varying cable tension force

圖14 分解后的加速度信號一階模態(tài)分量Fig.14 The first order modal component of the measured cable acceleration responses with linearly varying tension force

圖15 拉力線性變化時索瞬時頻率識別結(jié)果Fig.15 The identified instantaneous frequency of the cable with linearly varying cable tension force

表3 拉力線性變化時索瞬時頻率識別精度指標(biāo)Tab.3 IA of instantaneous frequency identification of the cable with linearly varying cable tension force %

其次，進(jìn)行拉力正弦變化時拉索的瞬時頻率識別。試驗時索的拉力成正弦變化，幅度為±10 kN。采集到的加速度響應(yīng)曲線如圖16所示。在采用連續(xù)小波變換得到圖17所示的小波量圖后，確定分量信號分布在[10 Hz，20 Hz]和[40 Hz,50 Hz]兩個頻率區(qū)間范圍內(nèi)，且數(shù)目為2。值得注意的是，圖17在[20 Hz，40 Hz]之間雖然有能量出現(xiàn)，但其并不連續(xù)且明顯小于[40 Hz,50 Hz]頻帶之間的能量，因此可判斷為噪聲干擾，在確定模態(tài)個數(shù)k時可忽略不計。然后對信號進(jìn)行VMD分解，得到如圖18所示的一階分量。最后，對分解后的一階分量信號進(jìn)行同步擠壓小波變換，得到的瞬時頻率識別結(jié)果如圖19所示。由圖19可知，VMD+SWT識別的瞬時頻率值雖然與理論值有一定的偏差，但其變換趨勢與理論值基本一致，且識別效果明顯好于HHT。雖然CWT也能夠識別出相應(yīng)的趨勢，但其識別的瞬時頻率曲線仍含有較多的毛刺。表4給出的瞬時頻率識別方法的精度指標(biāo)IA則再一次證明了VMD+SWT的準(zhǔn)確性。

圖16 拉力正弦變化時索的加速度響應(yīng)Fig.16 The measured cable acceleration responses with sinusoidally varying tension force

圖17 拉力正弦變化時實測信號小波量圖Fig.17 The wavelet scalogram of the measured cable acceleration responses with sinusoidally varying tension force

圖18 分解后的加速度信號一階模態(tài)分量Fig.18 The first order modal component of the measured cable acceleration responses with sinusoidally varying tension force

圖19 拉力正弦變化時索瞬時頻率識別結(jié)果Fig.19 The identified instantaneous frequency of the cable with sinusoidally varying cable tension force

表4 拉力正弦變化時索瞬時頻率識別精度指標(biāo)Tab.4IA of instantaneous frequency identification of the cable with sinusoidally varying cable tension force%

4 結(jié) 論

自然環(huán)境激勵下的土木工程結(jié)構(gòu)/構(gòu)件本質(zhì)上屬于時變和非線性系統(tǒng)，其響應(yīng)信號呈現(xiàn)非平穩(wěn)性。因此，本文基于變分模態(tài)分解和同步擠壓小波變換算法提出一種識別時變結(jié)構(gòu)響應(yīng)信號瞬時頻率的新方法。該方法不僅通過引入小波量圖快速便捷地確定響應(yīng)信號中的模態(tài)分量個數(shù)，解決了VMD這種非遞歸模式的頻率分割算法的典型問題，同時該方法引入的VMD算法解決了同步擠壓小波變換無法準(zhǔn)確分解響應(yīng)信號的問題，從而進(jìn)一步提高了同步擠壓小波變換識別瞬時頻率的精度。通過一個數(shù)值算例和兩個時變結(jié)構(gòu)試驗對提出的新方法進(jìn)行了驗證，研究結(jié)果表明：

(1)使用VMD時所需的模態(tài)分量個數(shù)可通過直接觀察小波量圖中頻率分量的個數(shù)直接判斷。但當(dāng)小波量圖上某一頻帶范圍內(nèi)的能量不連續(xù)且不明顯時，該分量可忽略不計。

(2)VMD算法能夠有效分離時變多分量響應(yīng)信號，一定程度上避免了EMD算法中的端點效應(yīng)和模態(tài)疊混等問題。

(3)基于變分模態(tài)分解和同步擠壓小波變換的聯(lián)合算法能夠有效識別時變結(jié)構(gòu)/構(gòu)件非平穩(wěn)響應(yīng)信號的瞬時頻率，且識別精度優(yōu)于HHT和CWT。