基于VMD和廣義Morse小波的結構瞬時頻率識別*

王 超， 毛 羚

(1.湖北工業大學土木建筑與環境學院 武漢， 430034) (2. 武漢理工大學土木工程與建筑學院 武漢，430070)

引 言

對于土木工程結構損傷診斷、健康監測與安全評估來說，模態參數識別是一個非常重要的研究課題。最近幾十年，各國學者提出了許多模態參數識別技術，如頻域分解法、時域分解法、隨機子空間識別、基于響應傳遞比的參數識別法[1]和基于貝葉斯理論的模態參數識別法[2]等。這些參數識別方法都是針對時不變系統。然而，實際土木工程結構在運營過程中由于環境或者荷載變化可能導致結構發生時變，如不同時刻構件不均勻溫度變換會引起結構模態頻率隨時間改變；橋上車輛荷載的變化引起結構拉索索力的變化，導致拉索剛度時變；車輛過橋時由于車輛位置不斷改變引起車橋系統質量時變。因此許多學者提出了時變參數識別方法進行系統識別，如狀態空間小波識別法[3]、S變換法[4]、HHT(Hilbert-Huang transform, 簡稱HHT)法[5]、自適應稀疏時頻分析法[6]、同步擠壓小波變換法[7]和解析模態分解法(analytical modal decomposition簡稱AMD)[8]等。

變分模態分解(variational mode decomposition,簡稱VMD)是最近提出的一種新的信號自適應分解方法[9],與HHT類似，它可以將非平穩的時變信號自適應的分解成多個具有有限帶寬的模態分量信號之和，但是它沒有HHT分解時遇到的模態混疊問題。由VMD分解得到模態分量需再采取合適方法提取信號的瞬時頻率。HHT是采用Hilbert變換的方法進行分析得到Hilbert譜，但是它會產生負頻率和端點效應問題。利用小波變換技術也可以提取信號瞬時頻率，通常采用的復Morlet小波并不是完全解析的，當中心頻率較小時可能產生負頻率泄露，導致時頻分析結果產生偏差。

因此，筆者采用VMD和廣義Morse小波(generalized Morse wavelets，簡稱GMW)相結合的技術來識別時變結構瞬時頻率。首先,利用信號在小波時頻域具有稀疏性的特性，引入快速閾值迭代技術對信號進行去噪處理;然后,采用VMD將去噪信號分解為多個模態分量;最后,采用GMW對每個模態分量進行時頻分析提取對應的瞬時頻率。通過數值算例和一個移動小車試驗驗證了提出的方法。

1 快速閾值迭代降噪

現場實測信號通常會受到噪聲污染,這里引入快速閾值迭代算法(fast iterative shrinkage-thresholding algorithm, 簡稱FISTA)進行降噪。含噪信號通過小波變換可以轉化到時間-尺度域，噪聲能量將分布到非常廣泛的時間-尺度空間，而有效信號能量將集中在一個局部的時間-尺度域空間，具有稀疏特性，因此，通過最優稀疏性可以實現將噪聲和信號盡量分離。假設實測的含噪信號y表示為

y=x+n=Iw(C)+n

(1)

其中:x為真實無噪信號;C為x小波變換系數矩陣;Iw(·)為小波逆變換;n為噪音。

稀疏性對應于l0范數最小化問題，對于非凸復雜性優化問題，可以采用l1范數來替代求解[10],因此，上述問題可以轉換為求解優化問題

(2)

其中：s.t.為約束條件;ε為小量。

式(2)有約束優化問題可以通過施加拉格朗日乘子正則化方法轉變為無約束優化問題來求解[11]

(3)

其中:λ為拉格朗日乘子。

這里采用Beck等[12]提出的快速閾值迭代算法(IST)通過下式快速迭代求解式(3)的解

(4)

其中：WT(·)為小波變換;Tλ(·)為軟閾值函數;k為迭代數;t為迭代參數。

其中：c為小波系數矩陣C中的小波系數。

為了提高收斂效率，λ按指數衰減變化取值

(7)

其中：N為迭代求解次數;λmax=0.99max|C0|;λmin=max|C0|/q，q為系數C0的數量。

通過上述算法對響應數據進行分析，可以有效降低噪音干擾。

2 變分模態分解

VMD是一種新的信號自適應分解技術，將實信號分解為具有最優特定稀疏特性的多個分量信號(模態)之和。這里的稀疏特性采用分解信號的頻譜帶寬來衡量， VMD要求分解的每個模態分量具有圍繞各自中心頻率最緊的帶寬。

為了估計每個模態分量的帶寬，VMD利用希爾伯特變換計算各模態對應的解析信號, 得到其單邊頻譜, 然后通過頻率平移將每個模態頻譜平移到各自估計的中心頻率位置，最后計算頻移后信號梯度范數的平方來估計各模態分量的帶寬。這樣最優模態分解轉化為如下有約束的變分優化問題來實現[13]

(8)

其中：{si}和{ωi}分別為分解的模態分量和對應的中心頻率;δ(t)為狄拉克函數;*為卷積;f為分解的信號。

通過引入二次罰函數項和拉格朗日乘子ζ來處理上式的重構約束條件，拉格朗日乘子ζ用來施加嚴格的約束條件，而二次罰函數項可以提高收斂性能，將式(8)的約束優化問題轉變為式(9)無約束的函數優化問題

(9)

其中：α為帶寬約束參數。

α增加可能導致分解的模態分量帶寬減小，中心頻率誤差估計增大；α減小可能導致分解的模態帶寬增加，包含更多的噪音成分。α建議取值區間為100～10 000, 文中取2 000。

通過交替方向乘子法(alternating direction method of multipliers,簡稱ADMM)可迭代求取原變分優化問題的解。步驟如下：

(10)

其中：F-1(·)表示傅里葉逆變換。

(11)

4) 按式(12)進行迭代更新

(12)

循環迭代步驟2～4直到滿足收斂條件

(13)

其中:ζ為迭代收斂容差，文中取1e-7。

3 廣義Morse小波瞬時頻率識別

對于N自由度的時變結構系統，其振動響應z(t)可分解為N階振型模態響應的疊加和

(14)

其中:zm(t)為第m階振型模態響應。

每階模態響應可表示為時變幅值Am(t)和相位φm(t)的乘積

xm(t)=Am(t)eiφm(t)

(15)

通過VMD對結構響應信號z(t)進行分析，可以將結構不同階振型模態分離成各模態分量，然后分別對各模態分量采用連續小波變換進行時頻分析，提取對應的小波脊線，由小波尺度與信號的頻率關系可以得到結構的瞬時頻率[7]。

連續小波時頻分析通常采用復Morlet小波，然而，復Morlet小波只是當中心頻率較大時近似解析，當中心頻率較小時可能產生負頻率泄露，導致分析結果產生一定偏差。文中引入GMW[13]進行時頻分析, 該小波是完全解析的，不存在負頻率泄露問題。GMW頻域定義為

Ψβ, γ(ω)=H(ω)2(eγ/β)β/γ)ωβe-ω γ

(16)

其中：H(ω)為階躍函數;β和γ為控制小波時域和頻域衰減的參數。

當γ=3時，GMW接近高斯函數，具有很好的時頻分辨率。文中采用的GMW參數為γ=3和β=27，對應的時域和頻域波形如圖1所示(示意圖，幅值無量綱)。

圖1 GMW時域和頻域波形Fig.1 Wave form of GMW in time domainand frequency domain

4 數值算例

為驗證提出的算法, 首先采用如下具有不同時變頻率成分組合成的振動信號進行數值驗證

(17)

信號采樣頻率為100 Hz，并添加10%的高斯白噪聲。這里噪聲強度10%定義為噪聲的2-范數與信號的2-范數之比。

首先,利用FISTA對含噪信號進行降噪處理, 降噪前后信號的時頻分布圖(又稱小波量圖，尺度坐標為指數坐標，無量綱)如圖2所示。

圖2 信號小波量圖Fig.2 Signal wavelet scalogram

圖中顏色亮度表示小波系數幅值大小，越亮表示越大，對應能量就越大。可以看出，信號成分在小波域分布具有稀疏性，比較集中，能量較大，而噪聲分布較廣，通過FISTA處理后噪聲明顯減少。

隨后,采用VMD對降噪后的信號進行分析，分解得到如圖3所示的4個模態分量，原信號各成分也顯示在圖中用于比較。

圖3 VMD分解的各模態分量(黑線)和對應的原信號xi(藍線)Fig.3 The decomposed mode components by VMD (dark line) and corresponding signal xi (blue line)

從圖可以看出,VMD將原信號不同成分完好的分解開來，其中第2個成分x2具有突變特性，VMD將其突變前后分解成兩個部分，分解的各部分與原信號對應各成分吻合的非常好。

采用GMW對VMD分解的各模態分量進行連續小波分析提取對應的小波脊線，識別其時變頻率。另外，也采用HHT的方法的識別結果用于對比，原信號的瞬時頻率的理論值用來驗證識別結果的準確性，分析結果如圖4所示。

圖4 信號各成分時變頻率識別結果Fig.4 The identified time-varying frequency of each component

由圖4可知，HHT方法識別的3個信號成分瞬時頻率受噪音干擾波動較大，對于信號突變成分突變時刻存在一定的偏差，而且HHT方法的端點效應影響很大。提出的方法3個信號成分瞬時頻率識別結果與理論值吻合的很好，只在端部由于小波端點效應稍有偏離，但是遠小于HHT的影響。

5 移動荷載試驗

5.1 試驗模型

為進一步驗證提出的方法，設計了一個移動小車試驗，采用馬達牽引小車通過一個簡支鋼板梁來模擬移動荷載作用下的結構時變特性。板梁尺寸為1 000 mm×150 mm×9 mm，質量為10.555 kg，采用圓鋼桿支撐模擬的簡支支撐。為避免小車進入和離開主梁時產生大的振動與沖擊，在板梁首尾分別連接了一段過渡的引梁和尾梁，但它們與主梁斷開。小車設有兩軸，軸距為0.117 m，質量為2.308 kg，小車內通過砝碼調整重量, 本試驗時小車總重為5.296 kg。主梁上設置了一根與主梁中線平行的穿過小車的的鋼絲，使小車能沿板梁中線移動。試驗模型如圖5所示。

圖5 試驗模型圖Fig. 5 The experimental mode

在梁底板起點和終點，1/4,1/2,3/4跨處布設了應變片，在節點3,5,7,9處(見圖6)布設了加速度計，測點布置如圖6所示。

圖6 測點布置圖Fig. 6 The layout of measure points

5.2 模態試驗

首先進行模態試驗。在節點4處用力錘錘擊主梁，采集結構的加速度響應，采樣頻率為1 000 Hz。進行模態分析并對結構的有限元模型進行校正。隨后將小車當做移動質量塊，利用校正后的有限元模型計算小車重心通過不同節點位置時結構的第1階模態頻率，結果如表1所示。

表1 小車重心通過不同節點時結構第1階頻率(FEA)Tab.1 The first frequency of structure when the vehicle cross different node (FEA)

5.3 移動小車試驗

通過馬達牽引小車勻速通過試驗梁，小車移動時導致主梁上質量分布發生改變, 引起結構的時變。采集響應的應變和加速度響應，采樣頻率為200 Hz。節點1和11處測得的應變響應如圖7所示。

圖7 節點1和節點11處測點應變響應Fig.7 Response of strain by sensor at node 1 and node11

通過首尾節點應變測點響應、主梁長度、小車軸距可以計算小車通過主梁時速度。應變峰值對應小車前后輪正好經過測點位置，由圖7測得前輪通過節點1的時刻為4.925 s，后輪通過節點11的時刻為7.895 s。由此估算得到小車速度為0.376 m/s。

采用提出的方法對節點9處加速度測點測得的響應進行分析識別結構的第一階瞬時頻率。為進行比較，也分別采用HHT和復Morlet連續小波變換方法(continuous wavelet transform,簡稱CWT)進行時變頻率識別，結果如圖8所示。

圖8 節點9處測點識別的瞬時頻率Fig.8 Identified instantaneous frequency by sensor at node 9

從識別的結果可以看出，采用提出了VMD和GMW方法識別的結果比較光滑，性能較好，通常的CWT識別結構波動稍大，主要存在不同模態間干擾和噪音干擾，HHT識別結果最差，端點效應影響很大。

采用GMW方法提取小波脊線時，對識別精度的影響主要是信號小波脊線對應的各階成分是否被充分激起，另外信號各階成分頻率間距太小，間距過密可能會造成能量的相互干擾。

隨后對所有測點(節點3,5,7,9)測得的加速度響應進行分析，首先用FISTA算法進行降噪，然后采用VMD和GMW相結合的方法識別結構的第1階瞬時頻率。由應變估計的小車移動速度和通過主梁的起止時間，可以推導出小車在不同時刻在主梁上的位置。結合表1中有限元計算的結果，可以得到小車通過主梁不同時刻的第1階瞬時頻率有限元計算結果，與提出的方法進行對比分析，如圖9所示。

圖9 所有測點識別的結構瞬時頻率Fig.9 Identified instantaneous frequency by all sensors

由圖可見，所有測點數據的識別結果與有限元計算結果吻合的較好，主要在兩端位置稍有偏差，分析偏差的原因:a.尾梁跟主梁連接位置還是存在一定沖擊振動干擾;b.小波端點效應有一定影響。另外，所有測點識別結果相差很小，表明提出的方法可以有效識別結構時變參數。

6 結束語

土木工程結構在運營過程中通常表現出時變特性，筆者提出了基于VMD和GMW的方法識別時變結構瞬時頻率。對于實測信號的噪聲污染問題，提出采用快速閾值迭代算法(FISTA)進行降噪處理。數值模擬結果表明，提出的FISTA算法通過將信號轉換到小波域利用信號稀疏特性可以較好地對信號進行降噪，VMD和GMW法可以有效準確地識別時變結構瞬時頻率。移動小車試驗進一步驗證了方法的有效性和精度，提出的方法比HHT和CWT具有更好的效果。