基于VMD和廣義Morse小波的結(jié)構(gòu)瞬時(shí)頻率識(shí)別*

王 超， 毛 羚

(1.湖北工業(yè)大學(xué)土木建筑與環(huán)境學(xué)院 武漢， 430034) (2. 武漢理工大學(xué)土木工程與建筑學(xué)院 武漢，430070)

引 言

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

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

因此，筆者采用VMD和廣義Morse小波(generalized Morse wavelets，簡(jiǎn)稱GMW)相結(jié)合的技術(shù)來(lái)識(shí)別時(shí)變結(jié)構(gòu)瞬時(shí)頻率。首先,利用信號(hào)在小波時(shí)頻域具有稀疏性的特性，引入快速閾值迭代技術(shù)對(duì)信號(hào)進(jìn)行去噪處理;然后,采用VMD將去噪信號(hào)分解為多個(gè)模態(tài)分量;最后,采用GMW對(duì)每個(gè)模態(tài)分量進(jìn)行時(shí)頻分析提取對(duì)應(yīng)的瞬時(shí)頻率。通過(guò)數(shù)值算例和一個(gè)移動(dòng)小車試驗(yàn)驗(yàn)證了提出的方法。

1 快速閾值迭代降噪

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

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

(1)

其中:x為真實(shí)無(wú)噪信號(hào);C為x小波變換系數(shù)矩陣;Iw(·)為小波逆變換;n為噪音。

稀疏性對(duì)應(yīng)于l0范數(shù)最小化問(wèn)題，對(duì)于非凸復(fù)雜性優(yōu)化問(wèn)題，可以采用l1范數(shù)來(lái)替代求解[10],因此，上述問(wèn)題可以轉(zhuǎn)換為求解優(yōu)化問(wèn)題

(2)

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

式(2)有約束優(yōu)化問(wèn)題可以通過(guò)施加拉格朗日乘子正則化方法轉(zhuǎn)變?yōu)闊o(wú)約束優(yōu)化問(wèn)題來(lái)求解[11]

(3)

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

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

(4)

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

其中：c為小波系數(shù)矩陣C中的小波系數(shù)。

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

(7)

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

通過(guò)上述算法對(duì)響應(yīng)數(shù)據(jù)進(jìn)行分析，可以有效降低噪音干擾。

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

VMD是一種新的信號(hào)自適應(yīng)分解技術(shù)，將實(shí)信號(hào)分解為具有最優(yōu)特定稀疏特性的多個(gè)分量信號(hào)(模態(tài))之和。這里的稀疏特性采用分解信號(hào)的頻譜帶寬來(lái)衡量， VMD要求分解的每個(gè)模態(tài)分量具有圍繞各自中心頻率最緊的帶寬。

為了估計(jì)每個(gè)模態(tài)分量的帶寬，VMD利用希爾伯特變換計(jì)算各模態(tài)對(duì)應(yīng)的解析信號(hào), 得到其單邊頻譜, 然后通過(guò)頻率平移將每個(gè)模態(tài)頻譜平移到各自估計(jì)的中心頻率位置，最后計(jì)算頻移后信號(hào)梯度范數(shù)的平方來(lái)估計(jì)各模態(tài)分量的帶寬。這樣最優(yōu)模態(tài)分解轉(zhuǎn)化為如下有約束的變分優(yōu)化問(wèn)題來(lái)實(shí)現(xiàn)[13]

(8)

其中：{si}和{ωi}分別為分解的模態(tài)分量和對(duì)應(yīng)的中心頻率;δ(t)為狄拉克函數(shù);*為卷積;f為分解的信號(hào)。

通過(guò)引入二次罰函數(shù)項(xiàng)和拉格朗日乘子ζ來(lái)處理上式的重構(gòu)約束條件，拉格朗日乘子ζ用來(lái)施加嚴(yán)格的約束條件，而二次罰函數(shù)項(xiàng)可以提高收斂性能，將式(8)的約束優(yōu)化問(wèn)題轉(zhuǎn)變?yōu)槭?9)無(wú)約束的函數(shù)優(yōu)化問(wèn)題

(9)

其中：α為帶寬約束參數(shù)。

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

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

(10)

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

(11)

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

(12)

循環(huán)迭代步驟2～4直到滿足收斂條件

(13)

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

3 廣義Morse小波瞬時(shí)頻率識(shí)別

對(duì)于N自由度的時(shí)變結(jié)構(gòu)系統(tǒng)，其振動(dòng)響應(yīng)z(t)可分解為N階振型模態(tài)響應(yīng)的疊加和

(14)

其中:zm(t)為第m階振型模態(tài)響應(yīng)。

每階模態(tài)響應(yīng)可表示為時(shí)變幅值A(chǔ)m(t)和相位φm(t)的乘積

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

(15)

通過(guò)VMD對(duì)結(jié)構(gòu)響應(yīng)信號(hào)z(t)進(jìn)行分析，可以將結(jié)構(gòu)不同階振型模態(tài)分離成各模態(tài)分量，然后分別對(duì)各模態(tài)分量采用連續(xù)小波變換進(jìn)行時(shí)頻分析，提取對(duì)應(yīng)的小波脊線，由小波尺度與信號(hào)的頻率關(guān)系可以得到結(jié)構(gòu)的瞬時(shí)頻率[7]。

連續(xù)小波時(shí)頻分析通常采用復(fù)Morlet小波，然而，復(fù)Morlet小波只是當(dāng)中心頻率較大時(shí)近似解析，當(dāng)中心頻率較小時(shí)可能產(chǎn)生負(fù)頻率泄露，導(dǎo)致分析結(jié)果產(chǎn)生一定偏差。文中引入GMW[13]進(jìn)行時(shí)頻分析, 該小波是完全解析的，不存在負(fù)頻率泄露問(wèn)題。GMW頻域定義為

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

(16)

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

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

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

4 數(shù)值算例

為驗(yàn)證提出的算法, 首先采用如下具有不同時(shí)變頻率成分組合成的振動(dòng)信號(hào)進(jìn)行數(shù)值驗(yàn)證

(17)

信號(hào)采樣頻率為100 Hz，并添加10%的高斯白噪聲。這里噪聲強(qiáng)度10%定義為噪聲的2-范數(shù)與信號(hào)的2-范數(shù)之比。

首先,利用FISTA對(duì)含噪信號(hào)進(jìn)行降噪處理, 降噪前后信號(hào)的時(shí)頻分布圖(又稱小波量圖，尺度坐標(biāo)為指數(shù)坐標(biāo)，無(wú)量綱)如圖2所示。

圖2 信號(hào)小波量圖Fig.2 Signal wavelet scalogram

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

隨后,采用VMD對(duì)降噪后的信號(hào)進(jìn)行分析，分解得到如圖3所示的4個(gè)模態(tài)分量，原信號(hào)各成分也顯示在圖中用于比較。

圖3 VMD分解的各模態(tài)分量(黑線)和對(duì)應(yīng)的原信號(hào)xi(藍(lán)線)Fig.3 The decomposed mode components by VMD (dark line) and corresponding signal xi (blue line)

從圖可以看出,VMD將原信號(hào)不同成分完好的分解開來(lái)，其中第2個(gè)成分x2具有突變特性，VMD將其突變前后分解成兩個(gè)部分，分解的各部分與原信號(hào)對(duì)應(yīng)各成分吻合的非常好。

采用GMW對(duì)VMD分解的各模態(tài)分量進(jìn)行連續(xù)小波分析提取對(duì)應(yīng)的小波脊線，識(shí)別其時(shí)變頻率。另外，也采用HHT的方法的識(shí)別結(jié)果用于對(duì)比，原信號(hào)的瞬時(shí)頻率的理論值用來(lái)驗(yàn)證識(shí)別結(jié)果的準(zhǔn)確性，分析結(jié)果如圖4所示。

圖4 信號(hào)各成分時(shí)變頻率識(shí)別結(jié)果Fig.4 The identified time-varying frequency of each component

由圖4可知，HHT方法識(shí)別的3個(gè)信號(hào)成分瞬時(shí)頻率受噪音干擾波動(dòng)較大，對(duì)于信號(hào)突變成分突變時(shí)刻存在一定的偏差，而且HHT方法的端點(diǎn)效應(yīng)影響很大。提出的方法3個(gè)信號(hào)成分瞬時(shí)頻率識(shí)別結(jié)果與理論值吻合的很好，只在端部由于小波端點(diǎn)效應(yīng)稍有偏離，但是遠(yuǎn)小于HHT的影響。

5 移動(dòng)荷載試驗(yàn)

5.1 試驗(yàn)?zāi)Ｐ?/h3>

為進(jìn)一步驗(yàn)證提出的方法，設(shè)計(jì)了一個(gè)移動(dòng)小車試驗(yàn)，采用馬達(dá)牽引小車通過(guò)一個(gè)簡(jiǎn)支鋼板梁來(lái)模擬移動(dòng)荷載作用下的結(jié)構(gòu)時(shí)變特性。板梁尺寸為1 000 mm×150 mm×9 mm，質(zhì)量為10.555 kg，采用圓鋼桿支撐模擬的簡(jiǎn)支支撐。為避免小車進(jìn)入和離開主梁時(shí)產(chǎn)生大的振動(dòng)與沖擊，在板梁首尾分別連接了一段過(guò)渡的引梁和尾梁，但它們與主梁斷開。小車設(shè)有兩軸，軸距為0.117 m，質(zhì)量為2.308 kg，小車內(nèi)通過(guò)砝碼調(diào)整重量, 本試驗(yàn)時(shí)小車總重為5.296 kg。主梁上設(shè)置了一根與主梁中線平行的穿過(guò)小車的的鋼絲，使小車能沿板梁中線移動(dòng)。試驗(yàn)?zāi)Ｐ腿鐖D5所示。

圖5 試驗(yàn)?zāi)Ｐ蛨DFig. 5 The experimental mode

在梁底板起點(diǎn)和終點(diǎn)，1/4,1/2,3/4跨處布設(shè)了應(yīng)變片，在節(jié)點(diǎn)3,5,7,9處(見圖6)布設(shè)了加速度計(jì)，測(cè)點(diǎn)布置如圖6所示。

圖6 測(cè)點(diǎn)布置圖Fig. 6 The layout of measure points

5.2 模態(tài)試驗(yàn)

首先進(jìn)行模態(tài)試驗(yàn)。在節(jié)點(diǎn)4處用力錘錘擊主梁，采集結(jié)構(gòu)的加速度響應(yīng)，采樣頻率為1 000 Hz。進(jìn)行模態(tài)分析并對(duì)結(jié)構(gòu)的有限元模型進(jìn)行校正。隨后將小車當(dāng)做移動(dòng)質(zhì)量塊，利用校正后的有限元模型計(jì)算小車重心通過(guò)不同節(jié)點(diǎn)位置時(shí)結(jié)構(gòu)的第1階模態(tài)頻率，結(jié)果如表1所示。

表1 小車重心通過(guò)不同節(jié)點(diǎn)時(shí)結(jié)構(gòu)第1階頻率(FEA)Tab.1 The first frequency of structure when the vehicle cross different node (FEA)

5.3 移動(dòng)小車試驗(yàn)

通過(guò)馬達(dá)牽引小車勻速通過(guò)試驗(yàn)梁，小車移動(dòng)時(shí)導(dǎo)致主梁上質(zhì)量分布發(fā)生改變, 引起結(jié)構(gòu)的時(shí)變。采集響應(yīng)的應(yīng)變和加速度響應(yīng)，采樣頻率為200 Hz。節(jié)點(diǎn)1和11處測(cè)得的應(yīng)變響應(yīng)如圖7所示。

圖7 節(jié)點(diǎn)1和節(jié)點(diǎn)11處測(cè)點(diǎn)應(yīng)變響應(yīng)Fig.7 Response of strain by sensor at node 1 and node11

通過(guò)首尾節(jié)點(diǎn)應(yīng)變測(cè)點(diǎn)響應(yīng)、主梁長(zhǎng)度、小車軸距可以計(jì)算小車通過(guò)主梁時(shí)速度。應(yīng)變峰值對(duì)應(yīng)小車前后輪正好經(jīng)過(guò)測(cè)點(diǎn)位置，由圖7測(cè)得前輪通過(guò)節(jié)點(diǎn)1的時(shí)刻為4.925 s，后輪通過(guò)節(jié)點(diǎn)11的時(shí)刻為7.895 s。由此估算得到小車速度為0.376 m/s。

采用提出的方法對(duì)節(jié)點(diǎn)9處加速度測(cè)點(diǎn)測(cè)得的響應(yīng)進(jìn)行分析識(shí)別結(jié)構(gòu)的第一階瞬時(shí)頻率。為進(jìn)行比較，也分別采用HHT和復(fù)Morlet連續(xù)小波變換方法(continuous wavelet transform,簡(jiǎn)稱CWT)進(jìn)行時(shí)變頻率識(shí)別，結(jié)果如圖8所示。

圖8 節(jié)點(diǎn)9處測(cè)點(diǎn)識(shí)別的瞬時(shí)頻率Fig.8 Identified instantaneous frequency by sensor at node 9

從識(shí)別的結(jié)果可以看出，采用提出了VMD和GMW方法識(shí)別的結(jié)果比較光滑，性能較好，通常的CWT識(shí)別結(jié)構(gòu)波動(dòng)稍大，主要存在不同模態(tài)間干擾和噪音干擾，HHT識(shí)別結(jié)果最差，端點(diǎn)效應(yīng)影響很大。

采用GMW方法提取小波脊線時(shí)，對(duì)識(shí)別精度的影響主要是信號(hào)小波脊線對(duì)應(yīng)的各階成分是否被充分激起，另外信號(hào)各階成分頻率間距太小，間距過(guò)密可能會(huì)造成能量的相互干擾。

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

圖9 所有測(cè)點(diǎn)識(shí)別的結(jié)構(gòu)瞬時(shí)頻率Fig.9 Identified instantaneous frequency by all sensors

由圖可見，所有測(cè)點(diǎn)數(shù)據(jù)的識(shí)別結(jié)果與有限元計(jì)算結(jié)果吻合的較好，主要在兩端位置稍有偏差，分析偏差的原因:a.尾梁跟主梁連接位置還是存在一定沖擊振動(dòng)干擾;b.小波端點(diǎn)效應(yīng)有一定影響。另外，所有測(cè)點(diǎn)識(shí)別結(jié)果相差很小，表明提出的方法可以有效識(shí)別結(jié)構(gòu)時(shí)變參數(shù)。

6 結(jié)束語(yǔ)

土木工程結(jié)構(gòu)在運(yùn)營(yíng)過(guò)程中通常表現(xiàn)出時(shí)變特性，筆者提出了基于VMD和GMW的方法識(shí)別時(shí)變結(jié)構(gòu)瞬時(shí)頻率。對(duì)于實(shí)測(cè)信號(hào)的噪聲污染問(wèn)題，提出采用快速閾值迭代算法(FISTA)進(jìn)行降噪處理。數(shù)值模擬結(jié)果表明，提出的FISTA算法通過(guò)將信號(hào)轉(zhuǎn)換到小波域利用信號(hào)稀疏特性可以較好地對(duì)信號(hào)進(jìn)行降噪，VMD和GMW法可以有效準(zhǔn)確地識(shí)別時(shí)變結(jié)構(gòu)瞬時(shí)頻率。移動(dòng)小車試驗(yàn)進(jìn)一步驗(yàn)證了方法的有效性和精度，提出的方法比HHT和CWT具有更好的效果。