ISEGST在非線性RC結構瞬時頻率識別中的應用

沈中祥，袁平平，劉 義

(1.江蘇科技大學 土木工程與建筑學院，江蘇 鎮江 212003；2.河海大學 港口海岸與近海工程學院，南京 210098；3.江蘇揚子江船業集團公司，江蘇 靖江 214532)

在RC結構中，由于鋼筋、混凝土的材料非線性特性，結構的振動響應一般為非平穩信號，因此其頻率具有時變特性[1-9]。準確地識別出非線性RC結構在振動荷載作用下的瞬時頻率等特征，對結構的健康監測、狀態評估和維修加固等起著至關重要的作用[10]。近年來，時頻分析技術在信號處理領域顯示出了巨大的潛力，已成為一個新的研究熱點。然而，以往的時頻分析方法都存在一定的不足，例如短時傅里葉變換[11](short-time Fourier transform，STFT)的窗函數寬度固定不變，使其時頻分析能力較弱。在處理多分量復雜信號時，Wigner-Ville分布[12](Wigner-Ville distribution，WVD)會受到交叉項的影響而產生虛假模態。連續小波變換[13](continuous wavelet transform，CWT)在處理非平穩信號時也可能會產生模糊的時頻脊線，這將不利于瞬時頻率的提取。

為了使時頻脊線的能量更集中，瞬時頻率的識別效果更精準，Stockwell等[14]將STFT和CWT相結合，提出了S變換(S-transform，ST)。ST克服了STFT窗口時寬不變的缺陷，能根據頻率的變化調整分析時寬和提供直觀的時間頻率特征。隨后，一系列針對窗函數的改進算法，即廣義S變換(generalized S-transform,GST)被提出。其中，Zidelmal等[15-16]對GST的窗函數進行了改進，并通過能量集中度(concentration measure, CM)[17]給出了GST中窗函數的參數優化算法。為了提高聚集性，Daubechies等[18]提出了同步擠壓變換(synchrosqueezed transform，SST)，該方法以小波變換為基礎，通過對時頻周圍的能量進行擠壓，有效提高了時頻能量的集中度。趙朋[19]提出了基于改進S變換時頻譜重排的時變結構瞬時識別方法。Yu等[20]以SST為基礎，進一步提出了多重壓縮變換(multi-synchrosqueezing transform, MSST)，該方法通過多次壓縮變換，有效地提高了時頻能量的聚集性。同時，Yu等[21]還借鑒SST算法，提出了同步提取變換(synchroextracting transform，SET)。隨后相關學者對其不斷進行改進，使SET的時頻分辨率顯著提高[22-23]。Chen等[24]提出了基于GST的高精度時頻分析同步擠壓算法，即同步擠壓廣義S變換(synchrosqueezing generalized S-transform，SSGST)。康佳星[25]將SET和GST聯合起來，提出了同步提取廣義S變換(synchroextracting generalized S-transform，SEGST)，并將其成功應用到地震信號分析中。

本文采用了一種改進的廣義S變換(improved generalized S-transformation, IGST)，并結合能量CM原理推導了窗函數參數的優化算法，最后聯合SET引入了ISEGST。為了探究ISEGST對非線性RC結構瞬時頻率的識別能力，在數值模擬方面，本文基于OpenSees結構分析軟件，建立了一個三層RC框架結構，并計算其在地震作用下的動力響應，通過ISEGST對其加速度響應信號進行時頻分析獲取瞬時頻率，從而驗證所提方法的可行性。試驗方面，應用ISEGST對某七層剪力墻振動臺試驗數據進行時頻分析和瞬時頻率提取，進一步驗證該方法的實用性。

1 理論知識

1.1 改進廣義S變換

ST是地球物理學家Stockwell于1996年在STFT和CWT基礎上提出的一種時頻分析方法，其特點是引入了寬度和頻率成反比的高斯窗，其表達式為

(1)

ST具有完全可逆性，其逆變換為

(2)

對ST的窗函數進行改進，得到GST為

(3)

(4)

將改進后的窗函數代入到ST中，得到IGST為

特別地，當p=0時，IGST退化為GST；當m=1，p=0，r=1時，IGST則退化為傳統的ST。

通過轉換，得到IGST在頻域的表達式為

(6)

1.2 參數優化算法

GST中的關鍵問題是如何選擇窗函數的參數，IGST同樣面臨著這個問題。本文采用能量CM的方法進行計算

(7)

將IGST數據進行歸一化處理

(8)

優化問題可以表示為

(9)

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

(10)

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

(11)

本文中，參考Zidelmal等和Moukadem等的研究，令m∈(0，3]，p∈[0,3]，r∈[0,1]，K=10,L=1 000。所以，最終設置的優化問題為

(12)

式(12)可以看作是多元函數的優化問題，即在一系列約束條件下，尋找一組參數值，使某個或某一組函數的目標值達到最優。通過優化算法對上述問題進行求解即可得到IGST窗函數的參數值。

1.3 同步提取算法

借鑒SET原理，可以得到ISEGST為

SEIGST(τ,f)=IGST(τ,f)·δ(f-ωi(t,f))

(13)

式中：δ為同步提取算子，即

(14)

通過IGST(τ,f)來計算每一個時頻系數對應的瞬時頻率ωi(t,f)

由式(15)可得

(16)

其中，

實際應用時，考慮到計算誤差，可以采用式(18)近似計算同步提取算子δ，即

因此，ISEGST具體為

ISEGST(τ,f)=

(19)

式中，Δf為頻率間隔，一般可取為1。

2 數值模擬

2.1 模型計算

本文設計了一個三層非線性鋼筋混凝土框架結構，其尺寸如圖1所示。混凝土板厚度為120 mm，柱構件截面尺寸為500 mm×500 mm，梁構件截面尺寸為500 mm×250 mm，如圖2所示。混凝土類型為C30，鋼筋型號為HRB400，其材料本構具體如圖3所示。樓面恒載為3.5 kN/m2，活載為2 kN/m2，采用1940年EL-Centro地震波對結構進行X軸方向的激勵，地震波最大峰值加速度為2.2 m/s2。采用OpenSees軟件對結構進行分析，從而得到結構的動力響應。在OpenSees模擬中，混凝土本構關系采用Concrete02，鋼筋本構采用Steel01，鋼筋和混凝土的本構關系見圖3。箍筋能夠約束混凝土，提高混凝土的抗壓強度。考慮到箍筋對混凝土強度的影響，本文將梁和柱的核心混凝土的強度提高了40%。

圖1 框架結構的平面及立面布置圖(mm)

圖2 梁、柱截面配筋(mm)

圖3 材料本構圖

2.2 瞬時頻率提取

通過OpenSees軟件計算分析，得到框架結構頂層中間節點的加速度響應曲線，如圖4所示，選取該加速度響應信號為研究對象，通過STFT、SST、SET、IGST、ISEGST等五種方法對其進行時頻分析研究。在IGST中，窗函數的參數m=1.346 9，p=3，r=0.917。

圖4 框架結構的加速度響應

時頻分析的結果如圖5所示。對比圖5(a)和圖5(d)可以看出IGST的能量聚集性更高，這更利于ISEGST的頻率提取。從圖5(b)、圖5(c)和圖5(e)中可以看出SST的時頻分析結果較差，SET和ISEGST的時頻分析結果較好，而且ISEGST的結果更加平滑，時頻能量的聚集性更高。

圖5 框架結構的時頻分析結果

框架結構的頻率識別結果如圖6所示。可以看出，三種方法都能有效識別出結構的瞬時頻率。SET的識別結果波動較大，SST的識別結果與ISEGST的識別結果幾乎相等，但是，ISEGST的頻率識別曲線更加平滑，波動較小，結果更加準確。

圖6 框架結構頻率識別結果

3 試驗驗證

為了進一步研究ISEGST對實際結構運動響應信號的瞬時頻率識別效果，本文引用加州大學圣地亞哥分校(USCD-NEES)Panagiotou等[26]設計完成的七層RC剪力墻振動臺試驗數據進行了時頻分析，其中地震波加速度響應如圖7所示，實測的結構加速度響應如圖8所示。

圖7 施加于剪力墻的地震波

圖8 剪力墻結構的加速度響應

對所得加速度響應信號進行時頻分析，結果如圖9所示。其中，IGST的窗函數參數m=3，p=2.137 9，r=0.575 9。通過圖9(b)、圖9(c)、圖9(e)可以看出IGST的時頻分析結果更加平滑，效果更為理想。

圖9 剪力墻的時頻分析結果

通過極值法對SST、SET、ISEGST時頻分析結果的頻率進行進一步提取，提取結果如圖10所示。可以看出SST的時頻分析結果較差，SEI和ISEGST的時頻分析結果更為理想。但是，相對于SET, ISEGST的結果波動較小，平滑性更高。

圖10 剪力墻結構的頻率識別結果

4 結 論

本文將ISEGST應用到非線性RC結構瞬時頻率識別中。對三層RC框架結構和七層剪力墻結構的加速度響應進行時頻分析，從而提取結構的瞬時頻率。數值模擬和試驗結果表明ISEGST具有以下優點：

(1)結合CM原理和優化算法，可以計算得到IGST中窗函數的調節因子參數，有效地提高了IGST的計算效率和精度。

(2)延續了IGSE和SET的優點，ISEGST識別出的瞬時頻率曲線更加平滑、準確，是一種準確性較高的非線性RC結構瞬時頻率識別方法。