基于改進多重同步擠壓廣義S變換的結構瞬時頻率識別研究

袁平平， 程雪莉， 王航航， 沈中祥， 任偉新， 張 健

(1. 江蘇科技大學 材料科學與工程學院，江蘇 鎮江 212100； 2. 江蘇新揚子造船有限公司，江蘇 靖江 214532； 3. 江蘇科技大學 土木工程與建筑學院，江蘇 鎮江 212100； 4. 江蘇科技大學 船舶與海洋工程學院，江蘇 鎮江 212100；5. 深圳大學 土木與交通工程學院，廣東 深圳 518061)

隨著社會的發展，土木工程逐漸趨于大型化、高層化、結構異形化，其在人們日常生活中扮演的角色也愈發重要[1]。土木工程結構的損傷必定會導致其動力特性發生變化，隨著損失的累積其一旦出現問題，輕則發生財產損失、影響人們的正常生活，重則造成人員傷亡[2-3]。Kijewski等[4]認為工程結構具有較長的周期運動，動力特性在其使用壽命內往往會隨著時間的推移而發生變化。因此，精確識別工程結構的瞬時頻率等時變動態參數至關重要，這將有助于工程結構的狀態監測、損傷識別及安全評估等。

時頻分析作為一種新興的信號處理方法，是研究非平穩信號的一個重要工具，其通過設計時間和頻率的聯合函數，把時域信號轉換到二維平面上，從而在時間域和頻率域上同時體現信號的能量和強度[5-6]。作為傳統時頻分析方法，短時傅里葉變換(short-time Fourier transform, STFT)、Wigner-Ville(Winger-Ville distribution, WWD)分布、S變換(S-transform, ST)、希爾伯特-黃變換(Hilbert-Huang Transform, HHT)等均存在一定的不足[7]。基于此，諸多學者進一步展開了相關研究。裴強等[8]研究了基于時頻分析ST的高層框架損傷識別。Zidelmal等[9]提出了基于緊湊型支持內核的ST。續秀忠等[10]提出了一種基于線性時頻表示和Hilbert變換的時變模態參數識別方法。Hou等[11]提出了一種基于連續小波的時變結構瞬時模態參數識別。連續小波變換與傅里葉變換相似，能夠實現信號變換與信號重構，且能保持原有信息的完整并具有可逆性。Daubechies等[12]提出了同步擠壓小波變換(synchrosqueezed wavelet transform, SST)，該方法以小波變換為基礎，通過對時頻周圍的能量進行擠壓，有效提高了時頻能量的集中度。張志禹等[13]將SST應用于地震信號瞬時屬性的提取。劉景良等[14-15]提出了基于SST的結構瞬時頻率識別并對其進行了改進。受同步擠壓變換算法的啟示，以實現理想的時頻分析為目標，Yu等[16]提出了一種新的時頻分析方法，即同步提取變換(synchroextracting transform, SET)。該方法是基于STFT的一種時頻分析方法，可提高瞬時頻率的識別精度。康佳星[17]研究了SET在地震信號分析中的應用。張文海等[18]提出了基于SET與Vold-Kalman濾波的燃機階次跟蹤方法。Li等[19]提出了用于地震時頻分析的時間同步提取廣義Chirplet變換。然而，SET無法實現信號的精確重構，因此，Yu等[20]進一步提出了多重同步壓縮變換(multi-synchrosqueezing transform, MSST)，該方法對STFT獲得的時頻圖執行多次同步壓縮處理，從而有效提高了時頻能量的聚集性。

本文對廣義S變換(generalized S-transform,GST)的窗函數進行了改進，并根據能量集中度(concentration measure, CM)[21]選取窗函數的參數，同時結合MSST，提出了一種新穎的改進多重同步擠壓廣義S變換(improved multi-synchrosqueezing generalized S-transform, IMSSGST)。在數值模擬方面，采用兩層剪切框架時變結構驗證該方法的準確性；試驗方面，采用某七層鋼筋混凝土(reinforced concrete, RC)剪力墻振動臺的試驗數據進行了時頻分析，用于驗證該方法在實際工程中的實用性和準確性。

1 理論基礎

1.1 改進廣義ST

ST是STFT和連續小波變換(continue wavelet transform,CWT)相結合聯合起來的一種時頻分析方法，其特點是引入了寬度和頻率成反比的高斯窗

(1)

式中：Sx(τ,f)為ST；x(t)為原始信號；g(t)為窗函數；τ為時移因子；f為頻率。g(t)的表達式為

(2)

將式(2)代入式(1)，可以得到ST的表達式，即

(3)

ST具有完全可逆性，其逆變換如式(4)所示

(4)

(5)

將式(5)代入式(1)即可得到GST，其計算結果為

(6)

(7)

該窗函數可以通過多項式來表達GST中的窗函數。將式(7)代入式(1)得到IGST(improred GST)的計算結果如式(8)所示

(8)

IGST的頻域表達式如式(9)所示

(9)

1.2 IGST的參數優化算法

本文采用能量CM的方法進行計算

(10)

對IGST進行歸一化處理，可得

(11)

則，IGST的優化問題可由式(12)表示

(12)

優化問題的約束條件與所分析窗口的寬度范圍有關，窗口不應太窄而改變時間分辨率，但也不能太寬而影響頻率分辨率，即

(13)

式中：Ts為采樣周期；f∈[fmin,fmax]；fmin=1 Hz；為滿足奈奎斯特采樣定理，fmax取采樣頻率的一半。式(13)可簡化為

(14)

在Zidelmal等和He等的研究中，相關學者令m∈(0,3]，p∈[0,3]，r=[0,1]，K=10，L=1 000。本文采用相同的邊界條件，因此，最終的優化問題如下所示

(15)

通過優化算法對上述問題進行求解即可得到IGST窗函數的參數值。

1.3 改進多重同步擠壓廣義S變換

原信號的瞬時頻率ωi(τ,f)可由式(16)得到

(16)

其中，

(17)

所以，對信號x(t)的IGST計算結果Gx(τ,f)進行單次壓縮，可得

(18)

進行多次擠壓，可得IMSSGST的計算結果如式(19)所示

…

(19)

其中，

(20)

2 數值模擬

2.1 框架結構模型

實際工程結構往往受到極限荷載或周期荷載的作用，結構的剛度也會隨之發生改變，其響應信號的瞬時頻率往往也會隨之改變。為了驗證該方法對結構瞬時頻率識別的可行性，建立兩層框架模型進行驗證，如圖1所示。該結構的剛度會隨著時間變化，結構的具體參數如表1所示。

圖1 兩層剪切框架結構Fig.1 Two-layer shear frame structure

表1 兩層剪切框架結構模型參數

框架結構相應的運動方程

(21)

結構的剛度按式(22)和式(23)變化，其時變剛度如圖2所示。

(22)

(23)

圖2 框架結構剛度Fig.2 Stiffness of frame structure

通過矩陣特征值計算，可得框架結構的理論瞬時頻率如圖3所示。

2.2 瞬時頻率識別

對該框架結構進行白噪聲激勵，用Runge-Kutta法求得結構的位移，速度及加速度響應。采樣頻率為50 Hz，采樣時間為30 s。其中：噪聲激勵如圖4所示；第一層結構加速度響應見如圖5所示。

圖4 白噪聲激勵Fig.4 White noise excitation

圖5 白噪聲激勵下第一層的加速度響應Fig.5 Acceleration response of first layer under white noise

為了確定IMSSGST的壓縮次數，本文采用瞬時頻率在整個時間歷程內的均方根值作為瞬時頻率識別精度的指標(index of accuracy，IA)。IA值越小，說明識別值與理論值越接近。隨著壓縮次數的增加，IA值會趨于平穩，以此來確定IMSSGST的壓縮次數。IA值的計算方法如下所示

(24)

式中：fd(t)為瞬時頻率識別值；fe(t)為瞬時頻率理論值。

MSST和IMSSGST識別的低頻的精度值基本保持一致，但在高頻處，IA值在第4次和第5次達到最小，識別的精度最高，如圖6所示。本文將采用一次、兩次和四次壓縮分別對信號進行頻率識別。

圖6 壓縮1～10次的IA值Fig.6 IA value compressed 1-10 times

對加速度響應信號進行IMSSGST及MSST時頻分析，壓縮次數為一次、兩次和四次，其中，通過優化算法得到IMSSGST的參數：m=3，p=2.155 1，r=0.645 4。時頻分析結果如圖7所示。從圖7可知，隨著壓縮次數的增加，IMSSGST、MSST的時頻能量更加集中。

提取四次壓縮后的頻率如圖8所示。從圖8可知，通過多次壓縮后，IMSSGST和MSST均能有效識別該框架結構的瞬時頻率，且識別值與理論值保持在較小的誤差范圍內，有較高的實用性。

圖7 時頻分析結果Fig.7 Time frequency analysis results

圖8 瞬時頻率識別結果Fig.8 Instantaneous frequency identification results

3 試驗驗證

為了驗證IMSSGST對實際工程結構瞬時頻率的識別效果，本文對加州大學圣地亞哥分校Panagiotou等[24]設計完成的七層RC剪力墻振動臺的試驗數據進行了時頻分析。施加于結構的地震波加速度及剪力墻的加速度響應，分別如圖9和圖10所示。

圖9 施加在剪力墻的地震波Fig.9 Seismic wave applied to the shear wall

圖10 剪力墻的加速度響應Fig.10 Acceleration response of shear wall

分別用IMSSGST、MSST對測得的加速度響應信號進行時頻分析，壓縮次數為一次、兩次和四次，時頻分析結果如圖10所示。在IMSSGST中，通過優化算法得到m=1.346 9，p=3，r=0.917 0。由圖11可知：隨著壓縮次數的增加，IMSSGST、MSST識別出的頻率脊線也更加準確、明顯；在相同壓縮次數下，與MSST相比，IMSSGST的曲線波動更小，這有利于提取結構的瞬時頻率。

圖11 剪力墻的多次壓縮時頻分析結果Fig.11 Time-frequency analysis results of shear wall

采用極值法提取IMSSGST、MSST四次壓縮后的瞬時頻率，頻率提取結果如圖12所示。從圖12可知，該RC剪力墻結構在地震激勵作用下的頻率范圍在0.8～2.1 Hz，IMSSGST和MSST都能有效地識別其頻率范圍，但是，IMSSGST識別的瞬時頻率的曲線更加光滑，曲線波動更小。

圖12 四次壓縮后剪力墻的頻率識別結果Fig.12 Frequency identification result of shear wall after four-synchrosqueezing

4 結 論

本文結合IGST和MSST，提出了新穎的IMSSGST，并將其應用到工程結構的瞬時頻率識別中。對兩層剪切框架和七層RC剪力墻結構進行了瞬時頻率識別，數值模擬和試驗結果表明：

(1) 將CM原理和參數優化算法結合起來，可以快速得到IGST中調節因子的參數，有效地縮短了IGST的計算時間。

(2) 通過多次壓縮后，IMSSGST識別的瞬時頻率更加準確，能量更加集中，有效地改善了IMSSGST的頻率識別效果。

(3) 與MSST算法相比，壓縮相同次數的情況下，IMSSGST識別的頻率曲線更加光滑，曲線波動更小，且與理論值更加接近，是一種準確、高效的時頻分析方法。