基于變分模態分解和同步提取變換識別時變結構瞬時頻率

唐 蕾， 黃天立， 萬 熹

(中南大學 土木工程學院，湖南 長沙 410075)

在實際運營過程中，由于環境和荷載的變化會導致結構出現時變特性，如斜拉橋拉索的索力變化、列車過橋引起的橋梁振動以及環境溫度等都會導致結構動力特性隨時間而改變。這使得結構具有時變特征，結構響應表現出非平穩特性。因此，開展時變結構的瞬時頻率等特征參數識別方法研究，對于結構健康監測具有重要意義。目前國內外對時變結構瞬時頻率識別的方法可以分為兩類。第一，基于時頻分析的方法，即從信號時頻譜圖中獲取時變結構瞬時頻率等參數，并進一步通過時頻重排、能量重組等方法提高識別精度；第二，基于信號自適應分解的方法，即將多分量信號分解成多個單分量信號，再對各個單分量信號進行瞬時頻率等參數識別。

時頻分析方法能同時在時域和頻域內分析信號的時頻特征，是分析非平穩信號的有力工具[1]。傳統的時頻分析方法，如Wigner-Ville分布、Gabor變換、短時Fourier變換(short-time Fourier transform，STFT)[2]、小波變換[3-4]等，已在時變模態參數識別和非線性系統識別等領域得到了廣泛的應用。為提高時頻分析方法的精度，發展了針對時頻譜進行能量重組的方法。Daubechies等[5]提出了同步壓縮變換(synchrosqueezing transform，SST)，該方法將小波變換后的小波時頻譜圖進行重組，將小波系數壓縮到中心頻率附近獲取更高精度的時頻表示。劉景良等[6-7]采用同步壓縮小波變換對Duffing系統自由振動響應和兩層剪切型框架模型在地震激勵下的響應進行了結構瞬時頻率識別。王超等[8]采用同步壓縮小波變換識別了移動車輛荷載作用下的橋梁時變參數。Oberlin等[9]在同步壓縮小波變換的基礎上，提出了同步壓縮Fourier變換(fourier-based synchrosqueezing transform，FSST)。徐曉迪等[10]利用同步壓縮短時Fourier變換，提取了高速列車軸箱振動加速度信號的瞬時頻率。FSST提高了時頻分辨率，但噪聲同樣也會被壓縮，使得該方法的抗噪性能差。針對FSST存在的問題，于剛等[11]提出同步提取變換(synchroextrcting transform，SET)，SET提取STFT譜圖在瞬時頻率位置上的時頻系數，具有較強的時頻聚焦性。使用SET識別時變結構瞬時頻率可以較好地剔除環境或溫度等引起的噪聲。應該指出，SET方法是針對STFT的改進，在處理多分量信號時，各分量成分信號的頻譜之間須具有一定的距離，才能獲得精度較高的時頻譜圖。因此，尋找適合的信號分解方法，有效分解近距離頻率成分分量對于SET的應用至關重要。

信號的自適應分解方法中，Huang等[12]提出的經驗模式分解(empirical mode decomposition，EMD)方法影響較大，通過EMD分解得到的固有模式函數(intrinsic mode function，IMF)，結合Hilbert變換可提取多分量信號的瞬時頻率等參數。Shi等[13-14]基于EMD和Hilbert變換識別了剛度慢變、周期改變和突變3種類型單自由度時變系統在自由振動情況下和多自由度時變系統在受迫振動情況下的瞬時頻率和阻尼系數。EMD通過包絡線擬合的方法將多分量信號分解成多個固有模式函數之和，但是在擬合的過程中容易出現欠包絡或過包絡，從而導致模態混疊。許多學者利用EMD分解的思路發展了許多多分量信號分解方法，如局部均值分解[15]、局部特征尺度分解[16]等。相比EMD分解方法這些方法在一定程度上改善了分解效果，但是由于這些方法仍然都是基于包絡線擬合的方法，不能完全解決模態混疊的問題。變分模態分解(variational mode decomposition，VMD)[17]是一種新近提出的信號自適應分解方法，它舍棄了包絡線擬合，將基于極值點模態獲取問題變為變分模型的構造，將分解方式變為變分模型的求解，解決了模態混疊的問題，同時VMD相當于自適應的維納濾波器組，具有較強的抗噪性能。結合信號分解方法與單分量信號參數識別方法可以實現結構瞬時頻率等參數的識別。趙亞軍等[18]利用VMD和經驗包絡法對密集模態和瞬態系統的瞬時頻率和瞬時阻尼比進行了識別。王超等[19]提出基于VMD和廣義Morse小波相結合的時變結構瞬時頻率識別方法，用具有時變特性的移動小車試驗驗證了所提方法的有效性。

針對SET不能分離頻率成分間隔相近的多分量信號的問題，本文結合VMD和SET，提出了高能量聚集性的時變結構瞬時頻率識別方法。首先，通過傅里葉變換確定預設模態數量，利用VMD對多分量信號進行分解處理得到多個模態分量；然后，采用SET對每個模態分量進行時頻分析獲取瞬時頻率；最后，將各個模態分量的時頻譜圖疊加得到完整的多分量信號時頻譜圖。該方法結合了VMD和SET的優勢，解決了SET處理具有近距離頻率成分多分量信號的不足，準確識別了時變結構的瞬時頻率。多分量時變信號和兩自由度時變結構自由振動響應信號數值算例和時變拉索試驗驗證了方法的正確性和適用性。

1 基本原理

1.1 變分模態分解

VMD是一種新近提出的信號自適應分解方法，它假設信號的每個模態分量具有不同中心頻率且具有圍繞各自中心頻率最緊的帶寬。VMD分解主要基于維納濾波、希爾伯特變換、頻率混合這3個概念，為了使得估計帶寬之和最小，將待分解信號引入變分模型，求解變分模型的最優解獲取信號的最佳分解。

1.1.1 變分模型構造

原始信號s(t)由K個分量信號構成，通過傅里葉變換等方法確定預設模態數量，VMD假定原始信號s(t)可被分解為K個IMF分量，定義各IMF分量為調頻調幅信號uk(t)=Ak(t)cos(φk(t))，構造變分模型。

對uk(t)進行Hilbert變換，得到解析信號

(1)

式中：δ(t)為Dirac函數；j為虛數單位；*為卷積運算。

利用解析信號的頻移特性，對其混合一個中心頻率e-jωkt，將uk(t)的頻譜移動到相應的基頻帶

(2)

計算頻移后信號梯度范數的平方估計各模態分量的帶寬，為使各IMF分量的帶寬之和最小，建立約束變分模型如下

(3)

式中：{uk}和{ωk}為VMD分解得到的模態分量和對應的中心頻率；s(t)為被分解信號。

1.1.2 變分模型求解

步驟1為了求解約束變分模型，引入拉格朗日乘子λ對約束嚴格控制，再引入二次懲罰因子α提高該變分問題的收斂性，將有約束變分問題轉化為無約束變分問題

(4)

(5)

(6)

(7)

式中，τ為迭代步長。

步驟4設定判別精度e>0，若滿足收斂條件式(8)，則分解結束，否則循環迭代步驟2和步驟3直到滿足收斂條件

(8)

1.2 同步提取變換

SET是一種新近提出的高分辨率時頻分析方法，其僅提取STFT譜圖在瞬時頻率位置上的時頻系數，大大減少了噪聲的影響，具有優越的抗噪性能；同時SET克服了同步壓縮Fourier變換在能量壓縮過程中將噪聲分量壓縮到時頻譜圖的缺陷。 SET基本步驟如下。

步驟1由STFT計算信號s(t)的時頻譜圖

(9)

令gω(u)=g(u-t)·ejωu，由于窗函數g通常采用實函數，故窗函數的復共軛等于它本身，根據Parseval定理，式(9)可以寫成

(10)

(11)

(12)

(13)

步驟3提取STFT在瞬時頻率位置的時頻系數。為盡可能減輕噪聲影響，采用δ函數僅提取ω=ω0處的時頻系數，得到一個具有高分辨率的時頻譜Te(t,ω)，

Te(t,ω)=Ge(t,ω)δ(ω-ω0(t,ω))

(14)

式中，δ(ω-ω0(t,ω))為同步提取算子(synchroextracting operator,SEO)。

考慮到計算誤差，同時實際應用中要用到SEO的實部，建議用式(15)計算SEO

(15)

式中，Δω為離散頻率之間的間隔。

對于多分量信號，STFT時頻譜圖為

(16)

根據相位信息估計的瞬時頻率為

(17)

此時，SET的表達式為

Te(t,ω)=Ge(t,ω)δ(ω-φ′(t,ω))

(18)

1.3 時變結構瞬時頻率識別流程

VMD可以有效地分解多分量信號，對于多自由時變結構也適用，解決了SET無法識別頻率間隔較近多分量信號的不足，同時VMD+SET方法利用了SET高能量聚集性的優點，可以識別時變結構的瞬時頻率。圖1給出了基于VMD+SET方法識別時變結構瞬時頻率的流程圖，其基本步驟如下。

圖1 VMD+SET方法識別時變結構瞬時頻率流程圖

步驟1基于待分析時變結構響應信號的傅里葉幅值譜圖確定預設模態數量K。

步驟2采用VMD對多分量信號進行分解得到K個模態分量。

步驟3采用SET對每個模態分量進行時頻分析獲取瞬時頻率。

步驟4將各個模態分量的時頻譜圖疊加得到完整的多分量信號時頻譜圖，從而揭示時變結構的時變特性。

2 數值算例

2.1 多分量時變信號

考慮一模擬多分量時變信號s(t)，由3個分量信號s1(t)、s2(t)和s3(t)組成，各分量瞬時頻率分別隨時間產生線性、突變和二次變化。取信號時長為60 s，采樣頻率為50 Hz，考慮到實際工作測試中噪聲的影響，對信號添加10 dB的高斯白噪聲，其時程曲線如圖2所示。

圖2 模擬的含噪多分量時變信號s(t)

(19)

采用高斯窗函數g(t)=e-πt2/0.322，時寬為2 s，對含噪信號s(t)分別進行STFT、FSST、SET處理得到時頻譜圖，如圖3所示，其局部放大圖如圖4所示。從圖3(a)和圖4(a)可以看出，信號經STFT變換后，其時頻譜圖上的頻率軌跡線比較模糊，雖然可以大致看出信號瞬時頻率的變化趨勢，但是不能準確地確定每個時刻各分量信號對應的瞬時頻率。從圖3(b)和圖3(c)可以看出，信號經過FSST和SET變換后，其時頻譜圖上的頻率軌跡線更為清晰，噪聲影響明顯減少。對比局部放大圖4(b)和圖4(c)可以發現SET的時頻譜圖能量較FSST更為集中，時頻分辨率更高，FSST的譜圖中的噪點較少，在能量壓縮過程中噪聲也不可避免地一同壓入。由此可見，SET能較好地識別頻率隨時間線性、突變和二次變化的信號，且比STFT和FSST的能量聚集性強，SET適用于識別時變結構的瞬時頻率。

圖3 含噪聲信號s(t)的時頻譜圖

圖4 含噪聲信號s(t)的時頻譜圖(局部放大圖)

進一步采用VMD結合SET的方法識別該多分量信號的瞬時頻率。根據傅里葉頻譜圖確定該多分量信號的模態數量K=4，用VMD對信號進行分解得到4個分量信號，如圖5所示。從圖5(b)和圖5(c)可以看出，VMD分解出了突變信號，但是單獨使用VMD無法判斷各模態分量的具體成分，結合SET對VMD分解的各模態分量進行瞬時頻率的識別，再將各個模態分量的時頻譜圖疊加得到完整的多分量時變信號的時頻譜圖，如圖6所示。從圖6可以看出，對比理論值(點線)，VMD結合SET的方法能準確識別頻率隨時間產生線性、突變和二次變化的3種時變信號的瞬時頻率。

圖5 VMD分解的各模態分量

圖6 含噪聲信號s(t)經VMD+SET處理的時頻譜圖

采用Rényi熵定量評價STFT、FSST、SET和VMD+SET 4種時頻分析方法的時頻能量聚集特性。Rényi熵是描述系統混亂度的一種衡量指標，熵值越低則代表時頻譜圖的能量聚集性更好，其計算公式為

(20)

式中：W(t,ω)為時頻系數；α為Rényi熵的次序，本文采用3次Rényi熵進行評價。

對信號s(t)分別施加信噪比(signal-noise ratio ,SNR)從0～20 dB的噪聲，計算分別采用STFT、FSST、SET和VMD+SET 4種時頻分析方法獲取的信號時頻圖的Rényi熵值，如圖7所示。從圖7可以看出，基于SET時頻譜圖計算得到的Rényi熵值較STFT和FSST時頻譜圖計算得到的Rényi熵值更低，驗證了SET具有較好的時頻聚焦能力。將信號先用VMD分解再采用SET處理，相當于對信號進行了降噪處理。因此，基于VMD+SET時頻譜圖計算得到的Rényi熵值最低。由此可見，在評價的4種方法中，VMD+SET方法具有最高的時頻能量聚集特性和最優的抗噪性能，可以用于高精度的時變結構瞬時頻率識別。

圖7 噪聲信號s(t)(信噪比0～20 dB)的STFT、FSST、SET和VMD+SET時頻譜圖的Rényi熵

2.2 兩自由度時變結構

圖8所示為一兩自由度時變結構，調整其剛度使其兩階自振頻率接近，不滿足SET處理多分量信號的限制條件。在一定的初位移和初速度條件下，結構發生自由振動，振動過程中假設結構質量和阻尼不隨時間變化，振動過程中由于損傷結構剛度隨時間變化，結構的運動方程為

圖8 兩自由度時變結構

(21)

式中，m1=m2=1 kg，c1=c2=0.15 N·s/m，k2=300 N/m，剛度k1按式(21)隨時間線性變化

(22)

設定初始條件x1(0)=0.3 m，x2(0)=-0.06 m，由四階Runge-Kutta法求解結構的自由振動響應，其采樣頻率為100 Hz，信號時長為15 s。考慮測量噪聲的影響，對質量塊m1的位移響應信號添加10 dB的高斯白噪聲，如圖9所示。

圖9 質量塊m1的含噪聲位移響應信號

根據“凍結法”的原理，假定短時間間隔內時變結構為時不變結構，通過結構動力學方法求解得到各時刻時不變結構的固有頻率，作為此時變結構的瞬時頻率理論值。圖10(a)～圖10(c)給出了對圖9所示含噪聲位移響應分別進行STFT、FSST、SET處理后的時頻譜圖，采用高斯窗函數g(t)=e-πt2/0.322，時寬為1 s。從圖10可以看出，STFT時頻譜圖中能量主要集中在4～6 Hz內，兩階頻率離得比較近，發散部分產生了交叉重疊，不能滿足間隔距離條件，導致兩階模態頻率沒有清晰地分離開；而FSST和SET作為對STFT 時頻譜圖的后處理手段，也不能將兩階模態分量清晰地分離開，由此使得低能量的一階模態頻率無法識別。采用VMD+SET方法對該兩自由度時變結構的自由振動位移響應進行分析，首先由VMD分解信號得到振幅逐漸衰減的兩階模態分量信號，如圖11所示，再利用SET獲得時頻譜圖，如圖12所示(圖12中實線為理論頻率值)。從圖12可以看出，VMD+SET方法很好地分離了時變兩自由度結構的兩階模態分量，且識別的瞬時頻率值與理論頻率值基本一致。

圖10 位移響應的時頻譜圖

圖11 VMD分解的兩階模態分量

圖12 兩自由度時變結構自由振動位移響應的VMD+SET時頻譜圖

3 時變拉索試驗驗證

為進一步驗證VMD+SET方法識別時變結構瞬時頻率的有效性，設計一個時變拉索裝置，通過改變拉索在不同時刻的張拉力，模擬拉索的剛度時變特性。

3.1 拉索試驗裝置

圖13給出了拉索試驗裝置的示意圖和現場布置圖[21]。試驗拉索為7φS5一根的鋼絞線，水平放置，彈性模量為E=1.95×105MPa，截面積為1.374×10-4m2，線質量密度為1.1 kg/m，兩錨固點間拉索總長度為4.55 m。拉索一端錨固于反力架上，一端用電液伺服加載系統(electro hydraulic servo,MTS)作動器對其施加拉力，如圖13所示。試驗開始時，首先對拉索施加20 kN的預拉力，然后通過MTS作動器調整索力的大小，使得拉索剛度隨時間發生變化，從而導致拉索固有頻率隨之改變。在索力變化過程中，通過力錘敲擊拉索，并利用在拉索中部安裝的一個加速度傳感器，采集拉索的豎向加速度響應。

圖13 拉索試驗裝置

試驗研究索力線性和正弦變化兩種工況，對比基于MTS實測索力換算得到的拉索理論頻率，驗證基于VMD+SET的時變結構瞬時頻率識別方法識別拉索瞬時頻率的正確性。

3.2 拉力線性變化時拉索瞬時頻率識別

試驗開始時，拉索預拉力為20 kN，拉索拉力通過MTS作動器以1.67 kN/s的速率線性增長，MTS作動器同步記錄拉力變化曲線，如圖14所示。索力變化過程中，采用力錘敲擊拉索后，同時利用拉索中部安裝的加速度傳感器采集拉索豎向沖擊加速度響應，如圖15所示。響應時長6 s，采樣頻率為600 Hz。

圖14 MTS實測拉力(線性)

圖15 索力線性變化時拉索的加速度響應

對實測加速度響應信號采用VMD+SET方法識別拉索的瞬時頻率，采用時寬為0.5 s的高斯窗函數g(t)=e-πt2/0.322，圖16給出了識別得到的拉索第1階模態瞬時頻率以及線性拉力作用下，基于MTS實測拉力數據換算得到的拉索第1階模態瞬時頻率理論值(虛線)。從圖16中可以看出，拉索的瞬時頻率線性增加，識別的瞬時頻率與理論瞬時頻率非常接近，識別準確度高。應該指出，由于振動衰減，振幅減小，噪聲影響較大，加速度響應信號尾部的VMD+SET時頻譜圖不太明顯，因此識別的瞬時頻率值稍有偏差。此外，信號兩端由于“端點效應”影響，識別結果與理論結果之間也存在一定的誤差。

圖16 理論頻率與識別頻率結果(線性)

3.3 拉力正弦變化時拉索瞬時頻率識別

試驗開始時，拉索預拉力為20 kN，拉索拉力通過MTS作動器按正弦規律變化，變化幅度為±4 kN，MTS作動器同步記錄拉力變化曲線，如圖17所示。索力變化過程中，采用力錘敲擊拉索后，同時利用拉索中部安裝的加速度傳感器采集拉索豎向沖擊加速度響應，如圖18所示。響應時長6 s，采樣頻率為600 Hz。

圖17 MTS實測拉力(正弦)

圖18 索力正弦變化時拉索的加速度響應

對實測加速度響應信號采用VMD+SET方法識別拉索的瞬時頻率，采用時寬為0.5 s的高斯窗函數g(t)=e-πt2/0.322，圖19給出了識別得到的拉索第1階模態瞬時頻率以及正弦變化拉力作用下，基于MTS實測拉力數據換算得到的拉索瞬時頻率理論值(虛線)。從圖19中可以看出，識別的拉索瞬時頻率略低于理論頻率值，與理論值稍有偏差，能明顯看出具有正弦變化規律。同樣由于振動衰減，以及“端點效應”影響，識別結果與理論結果之間也存在一定的誤差。

圖19 理論頻率與識別頻率結果(正弦)

4 結 論

本文提出了一種結合VMD和SET識別時變結構瞬時頻率的方法，并進行了數值模擬和試驗驗證，結果表明：

(1)多分量時變信號數值模擬結果表明，VMD+SET方法能準確識別多分量信號中頻率隨時間線性、突變和二次變化的時變分量信號的瞬時頻率，具有很好的時頻能量聚集特性和抗噪性能，可用于高精度的時變結構瞬時頻率識別。

(2)兩自由度時變結構數值模擬結果表明，VMD+SET方法能很好地分離了具有近距離頻率成分的結構振動信號，解決了SET處理具有近距離頻率成分多分量信號的不足，能準確識別時變結構的瞬時頻率。

(3)時變拉索試驗表明，VMD+SET方法可準確識別索力線性工況下的瞬時頻率，識別精度高，正弦工況下的識別結果與理論結果雖然略有偏差，但是可以反映瞬時頻率正弦變化規律，驗證了該方法識別時變結構瞬時頻率的適用性。