變分模態分解和同步擠壓小波變換識別時變結構瞬時頻率

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

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

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

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

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

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

1 基本理論

1.1 變分模態分解

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

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

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

(5)

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

(6)

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

(7)

1.2 同步擠壓小波變換

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

(8)

(9)

(10)

1.3 基于變分模態分解和同步擠壓小波變換識別時變結構瞬時頻率

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

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

2 數值算例驗證

考慮如下多分量信號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。分量信號對應的瞬時頻率理論值分別為f1=1+0.25cos(0.5πt) Hz和f2=2+0.5cos(πt) Hz。為考慮實際噪聲影響，對信號添加20%水平的高斯白噪聲，噪聲強度由信噪比定義。

(12)

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

圖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 瞬時頻率識別結果Fig.5 The instantaneous frequency extracted by the proposed method and HHT

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

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

(13)

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

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

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

3 試驗驗證

為驗證變分模態分解和同步擠壓小波變換聯合方法識別時變結構瞬時頻率的有效性和準確性，設計一個質量突變懸臂梁試驗與一個剛度時變拉索試驗對其進行驗證。

3.1 質量突變懸臂梁試驗

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

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

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

圖7 實測響應信號Fig.7 The measured acceleration responses

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

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

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

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

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

3.2 剛度時變拉索試驗

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

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

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

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

圖12 拉力線性變化時索的加速度響應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 分解后的加速度信號一階模態分量Fig.14 The first order modal component of the measured cable acceleration responses with linearly varying tension force

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

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

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

圖16 拉力正弦變化時索的加速度響應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 分解后的加速度信號一階模態分量Fig.18 The first order modal component of the measured cable acceleration responses with sinusoidally varying tension force

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

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

4 結 論

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

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

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

(3)基于變分模態分解和同步擠壓小波變換的聯合算法能夠有效識別時變結構/構件非平穩響應信號的瞬時頻率，且識別精度優于HHT和CWT。