基于向上取整的多重同步擠壓變換識別時變結構瞬時頻帶

2022-11-21 04:11:58李宇祖劉景良廖飛宇駱勇鵬

振動與沖擊 2022年21期

李宇祖, 劉景良，2, 廖飛宇，2, 駱勇鵬，2, 王芳

(1.福建農林大學交通與土木工程學院，福州 350108；2.“數字福”智能交通技術物聯網實驗室，福州 350108;3.福建省鍋爐壓力容器檢驗研究院，福州 350008)

在役土木工程結構所處的自然環境和受到的激勵對其動力學參數有顯著影響[1]。如車輛(行人)通過橋梁結構時，二者組成一個新的具有移動質量的系統，而且系統的質量分布和剛度分布隨車輛(行人)位置的移動而不斷發生變化[2]。環境激勵下的大跨度斜拉橋的索力變化將影響拉索的剛度并導致結構的模態參數隨時間變化，而且由于斜拉索的密布特點有可能造成各階模態密集或疊混[3]。近年來不斷涌現的新型組合橋梁的各個構件通過螺栓、銷釘等剪力連接件形成一個整體而共同工作，具有良好的力學性能[4]。然而，由于連接件的存在和組成材料的不同，組合結構橋梁本身表現出復雜的非線性力學行為。此外，在役組合結構橋梁在車輛、行人、風載、溫度等多重荷載作用下有較大概率發生橋面板腐蝕、螺栓連接件松動等損傷現象，其結果必然導致結構的剛度、質量等物理參數發生變化，而具有剛度或質量突變的結構系統的動力特征參數必將隨時間發生變化[5-6]。因此，環境激勵下的在役土木工程結構本質上屬于時變和非線性結構系統，其響應信號呈現非平穩特性。顯然，從非平穩信號處理的角度識別環境激勵下橋梁結構的時變模態參數更符合實際情況，對于深入理解結構損傷診斷、有限元模型修正、振動控制和安全評估具有重要的理論意義和工程應用價值。

時頻分析方法是一種將一維時域信號映射到二維時頻面的信號分析方法，能夠全面反映非平穩信號的時頻聯合特征，十分適用于非平穩信號的分析處理。然而，現有的時頻分析方法存在一些不足之處，如短時傅里葉變換[7](short time Fourier transform, STFT)雖然具有一定的局部分析能力，但是由于其窗函數的固定將導致STFT無法同時獲得較高的時間分辨率和頻率分辨率。Wigner-Ville分布[8-9]具有良好的時頻聚集性，但是當信號為頻率相近的多分量信號時會出現交叉項干擾，而時頻圖將出現模糊不清的現象。希爾伯特-黃變換[10-11]中的經驗模態分解算法雖然能夠成功地將原始信號分解成若干個本征函數(intrinsic mode function, IMF)，但是有可能出現同一頻率成分被分解到相鄰IMF的現象。小波變換[12]作為一種在低頻部分具有較高的頻率分辨率而在高頻部分具有較高的時間分辨率的時頻分析方法，可以更為客觀地反映結構的瞬時特征參數變化。然而，小波變換在分析土木工程結構響應信號中的長周期信號分量時，仍然缺少足夠的頻率精度，即時頻圖中低頻信號分量的時頻曲線分布較寬[13]。

Daubechies等[14]提出的同步擠壓小波變換通過重組小波變換后的時頻圖獲得了更高的頻率精度，然而該方法將所有尺度上的小波系數擠壓到估算的瞬時頻率位置，因而無法同時處理時間和頻率兩個維度方向的擴散。為解決這一問題，Yu等[15]提出了同步提取變換(synchroextracting transform, SET)。SET通過保留估算瞬時頻率位置處的時頻系數提高了時頻分辨率，然而SET的信號重構能力隨著信號非線性的增加而減弱。為此，Yu等[16]提出了多重同步擠壓變換(multi-synchrosqueezing transform，MSST)。該方法通過對估算的瞬時頻率進行多次迭代從而提高了時頻精度，然而它在估算瞬時頻率的過程中采用的四舍五入方法造成了時頻系數在同一時刻無法重排至一點，即未重排點問題。由于未重排點的存在，MSST識別的頻帶精度有所降低。于是Yu[17]提出改進的多重同步擠壓算法(improved multi-synchrosqueezing transform，IMSST)來改善未重排點問題，但是IMSST并未徹底解決這一問題。

基于此，本文提出基于向上取整的多重同步擠壓變換算法(ceiling based multi-synchrosqueezing transform, CBMSST)。它首先引入Burg算法對自回歸(auto regressive, AR)模型的參數進行估計[18]，然后通過AR模型對信號進行有效延拓，從而解決了MSST及IMSST中的端點效應問題。最后，通過對估算的瞬時頻率小數部分進行兩次向上取整，從而達到消除未重排點的目的。通過兩個數值算例以及一個拉索試驗對所提方法進行驗證，研究結果表明：CBMSST不但徹底排除了未重排點，而且提高了頻帶的識別精度。

1 基本理論

1.1 信號延拓

由于MSST和IMSST的端點效應，其時頻圖在端點處存在未重排點現象。為解決這一問題，擬引入AR模型對信號進行延拓。然而，采用AR模型進行信號延拓首要解決的一個問題就是參數φ的估計。截至目前，AR模型的參數估計方法主要有最小二乘法、最大似然法、Yule-Walker、Burg算法等[19]。其中，Burg算法不但避開了序列的自相關函數估計，而且具有良好的頻率分辨率，因此本文最終選用Burg算法對AR模型的參數φ進行估計，從而實現了對數據序列或者目標信號的延拓。

AR模型是利用前期若干時刻隨機變量的線性組合來預測后續某一時刻隨機變量的線性回歸模型,其基本理論如下。

設定Yn={y(1),y(2),…,y(n)}為一時間序列，其中y(n+1)是前n個序列的線性組合和誤差項函數之和，其表達式為

(1)

式中:φk為模型參數;en為服從均值為0、方差為δ的正態分布的白噪聲。

在建立AR模型的基礎上，Burg算法首先根據第k項的前向(第1～第k項)和后向(第k～第n項)預測誤差的平均功率最小準則估計出反射系數，然后再運用Levinson遞推公式求解AR模型的參數φ。Burg算法的具體步驟如下所示。

首先，定義前向和后向預測誤差fk(n)和gk(n)，其表達式分別為

(2)

(3)

然后，采用格型濾波器結構計算各階前向和后向誤差可得

(4)

式中，λk為格型濾波器反射系數。

而前向和后向預測誤差的平均功率如式(5)所示。

(5)

(6)

然后根據如式(7)所示的Levinson遞推公式求解AR模型的參數φk-1,φk-2直至φ0。

(7)

最后，將求得的參數(φk-1,φk-2直至φ0)代入式(1)可得y(n+1)。以此類推，可以逐步求解y(n+2)，y(n+3)，…，y(2n-1)，y(2n)，即實現了對序列Yn的前向延拓。在這里需要指出的是：Burg算法也可實現序列Yn的后向延拓，其算法與前向延拓過程類似，具體算法參見孫志強等的研究。

1.2 基于向上取整的多重同步擠壓變換

同步擠壓變換可以通過對時頻系數進行多次迭代以取得較高的時頻分辨率，但是對時頻系數進行多次迭代將導致計算效率較低。為解決這個問題，Yu等在同步擠壓變換的基礎上提出了MSST。MSST摒棄了對時頻系數進行多次迭代的思路，轉而對估算的瞬時頻率進行多次迭代從而在提高時頻精度的同時也減少了計算時間。然而，正是因為四舍五入取整算法的引入導致了未重排點的產生。為此，Yu在MSST的基礎上又提出了IMSST。該算法對擠壓后的瞬時頻率進行兩次四舍五入以改善未重排點問題，但是并未完全解決這一問題。基于此，本文提出一種基于向上取整的多重同步擠壓變換來徹底解決重排點問題，其具體算法如下所示。

設定如式(8)所示的調幅調頻信號。

s(t)=A(t)eiφ(t)

(8)

式中，A(t)和φ(t)分別為信號的幅值和相位。

若ε足夠小，且|A′(t)|≤ε和|φ?(t)|≤ε，根據二階泰勒展開公式，A(u)和φ(u)在t時刻可表示為A(u)=A(t)+O{A′(t)}，φ(u)=φ(t)+φ′(t)(u-t)+0.5φ″(t)(u-t)2+O{φ?(t)}，其中O{A′(t)}和O{φ?(t)}分別為A′(t)和φ?(t)的無窮小，可忽略不計，因此式(8)改寫為

s(u)=A(t)ei{φ(t)+φ′(t)(u-t)+0.5φ″(t)(u-t)2}

(9)

以高斯窗函數作為STFT的窗函數，對式(9)進行短時傅里葉變換，則

G(t,ω)=

(10)

式中:ω為頻率;τ=(u-t)。

對式(10)求關于t的偏導，可得如(11)所示的估算瞬時頻率表達式。

(11)

然后，取式(11)的實部作為瞬時頻率估算值，如式(12)所示。

(12)

(13)

(14)

(15)

由于

(16)

當迭代次數足夠多時，若存在足夠小的ε，則有

(17)

(18)

(19)

(20)

(21)

為解決這一問題，本文摒棄了四舍五入法，轉而提出基于向上取整的多重同步擠壓變換算法，如式(22)所示。

(22)

Ts{N}(t,η)=

(23)

由于基于向上取整的多重同步擠壓變換不但通過信號延拓改善了端點效應，而且通過兩次向上取整解決了MSST及IMSST存在的未重排點問題。同時，它也在一定程度上提高了頻帶的識別精度，因而不失為一種可靠的時頻分析新方法。

2 數值算例

2.1 單分量調頻信號

考慮如式(24)所示的單分量調頻信號。該信號的采樣頻率為100 Hz，采樣時間為6 s，理論瞬時頻率為f=0.1t2+10 Hz，其時域波形如圖1(a)所示。對信號分別進行MSST,IMSST及CBMSST變換，結果如圖1(b)～圖1(d)所示。為更好地對比MSST，IMSST和CBMSST，對圖1(b)～圖1(d)進行局部放大，結果如圖2所示。

(a) x(t)

(b) MSST

(d) CBMSST圖1 x(t)時域波形及其瞬時頻帶識別結果Fig.1 x(t) and its instantaneous frequency bands extracted by various methods

(a) MSST

(b) IMSST

(24)

由圖2(a)、圖2(b)可知：雖然IMSST未能完全消除未重排點，但是其未重排點的重疊長度明顯小于MSST，這說明IMSST的識別效果優于MSST。相對而言，圖2(c)不存在未重排點，這說明CBMSST已經完全解決了MSST及IMSST中的未重排點問題。同時，圖3也給出了上述3種方法在端點處的識別效果。由圖3可知：由于CBMSST進行了信號延拓，其在端點處并未產生未重排點。這證明CBMSST在端點處的識別效果最佳。因此，不論是端點處還是非端點處，CBMSST通過向上取整和引入AR模型排除了瞬時頻帶中的未重排點，其識別精度也明顯高于MSST及IMSST。

(a) MSST

(b) IMSST

為方便對比，對x(t)也進行了基于向下取整的多重同步擠壓變換，結果如圖4所示。雖然基于向下取整的多重同步擠壓變換也排除了未重排點，但是其瞬時頻率識別精度如何尚未可知。為此，引入動態規劃法對基于向上取整和向下取整兩種方式得到的瞬時頻帶進行脊線提取并與理論值進行對比，結果如圖5(a)所示。圖5(b)同時也給出了圖5(a)的局部放大圖。由圖5(b)可得：相對于向下取整而言，通過向上取整方式得到的瞬時頻率識別結果與理論值更加接近，因而識別精度更高，這也是本文采用向上取整而非向下取整方式來識別瞬時頻帶的主要原因。

圖4 向下取整的多重同步擠壓變化的瞬時頻帶識別結果Fig.4 The instaneous frequency band of x(t) by flooring based multi-synchrosqueezing transfrom methods

圖5 x(t)的瞬時頻率識別結果Fig.5 The instaneous frequency band of x(t) by various

2.2 多分量調頻信號

考慮如式(25)所示的多分量調頻信號。設定信號的采樣頻率為100 Hz，采樣時間為6 s，其時域波形如圖6(a)所示。對信號分別進行MSST,IMSST及CBMSST變換，結果如圖6(b)～圖6(d)所示。為更好地進行比對，對圖6(b)～圖6(d)進行局部放大，結果如圖7和8所示。

(a) y(t)

(b) MSST

(d) CBMSST圖6 y(t)及其瞬時頻帶識別結果Fig.6 y(t) and its instantaneous frequency bands extracted by various methods

(a) MSST

(b) IMSST

(a) MSST

(b) IMSST

y(t)=y1(t)+y2(t)

(25)

式中:y1(t)=sin{2π[17t+6cos(1.5t)]};y2(t)=sin{2π[35t+6cos(1.5t)]}。

由圖7(a)、圖7(b)可知：雖然IMSST提取的瞬時頻帶中仍然存在未重排點，但其個數已經明顯少于圖7(a)，這說明IMSST的識別效果優于MSST。與此對應的是，圖7(c)中不存在未重排點，因此CBMSST的頻帶識別效果好于IMSST。同理，通過對比圖8(a)～圖8(c)也可得出類似結論。圖9和圖10也給出了上述3種方法在端點處的頻帶識別效果。基于CBMSST的瞬時頻帶在端點處并未出現未重排點，這再次驗證了該方法在端點效應抑制上的優越性。總而言之，相比MSST及IMSST，CBMSST避免了未重排點，因而提高了瞬時頻帶識別精度。

(a) MSST

(b) IMSST

(a) MSST

(b) IMSST

3 試驗驗證

本章采用Wang等研究中提供的時變拉索結構試驗數據來驗證本文所提的方法的有效性及準確性。該試驗通過連續變更索的拉力來改變其剛度，從而實現該結構的時變特性。

試驗采用的拉索為一根7Φ5鋼絞線，彈性模量E=1.95×105MPa，截面積為1.374×10-4m2，線密度為1.1 kg/m。拉索的一端采用反力架錨固，而另一端則通過電液伺服加載系統(mechanical testing and simulation，MTS)的作動器施加拉力。兩個錨固點之間的索長為4.55 m。整個試驗裝置,如圖11所示。激勵設備采用LC-03沖擊力錘，試驗開始前，對拉索施加20 kN的預拉力，然后通過MTS作動器改變索的拉力，使其剛度隨時間發生變化，在改變索力的同時采用沖擊力錘敲擊拉索并通過德國HBM公司生產的動態信號采集系統和加速度傳感器采集拉索的豎向加速度響應。

圖11 拉索試驗裝置Fig.11 The cable test setup

設定采集儀的采樣頻率為600 Hz，采樣時間為6 s。拉索的拉力線性變化及信號采集儀采集的加速度響應數據，如圖12和圖13所示。對圖13所采集到的加速度響應分別進行MSST，IMSST和CBMSST分析，結果如圖14所示。為更好地對比MSST，IMSST，CBMSST識別結果，對圖14進行局部放大，結果如圖15所示。由圖15(a)和圖15(b)可知：雖然基于IMSST的時頻圖中的未重排點仍然存在，但是其重疊長度明顯較短，因此IMSST的識別效果相對MSST更佳。由于圖15(c)中不存在未重排點，這說明CBMSST的頻帶識別效果最佳。從圖15(a)～圖15(c)可以看出：由于信號延拓的作用，CBMSST在端點處(0附近)不存在未重排點，這再次驗證了CBMSST在端點處的識別效果優于既有的MSST和IMSST。由此可見，本文提出的CBMSST解決了MSST及IMSST存在的未重排點問題，提高了頻帶的時頻聚集性。

圖12 實測線性變化拉力Fig.12 The measured cable tension force with linear variation

圖13 拉索加速度響應Fig.13 The measured cable acceleration response

(a) MSST

(b) IMSST

(c) CBMSST圖14 采用3種方法識別拉索響應信號的瞬時頻帶Fig.14 The instantaneous frequency band of the response from the cable extracted by various methods

(a) MSST

(b) IMSST

(c) CBMSST圖15 位于0～1 s的瞬時頻帶識別結果Fig.15 The instantaneous frequency band of the response from the cable between 0 sec and 1 sec

此外，為明確CBMSST的頻率識別結果與理論頻率之間的誤差，通過動態規劃法對采用CBMSST得到的頻帶進行瞬時頻率曲線的提取并與理論值進行對比，結果如圖16所示。此外，在表1中給出了理論頻率值、識別值以及二者之間的相對誤差。其中，拉索的固有頻率的理論值通過“凍結法”求解。相比理論值，雖然通過CBMSST識別的瞬時頻率值存在一定的誤差，但是誤差基本在1%左右，最大值也不超過2.8%。因此，CBMSST具有很好的識別精度。

圖16 拉索響應信號瞬時頻率曲線識別結果Fig.16 The identified instantaneous frequency curve of the response from the cable

表1 拉力線性變化工況下不同瞬時頻率識別方法對應的相對誤差Tab.1 Relative error of instantaneous frequency identification of the cable with linearly varying cable tension force using various methods

4 結論

為解決MSST和IMSST算法存在未重排點的問題，本文提出一種新的基于向上取整的多重同步擠壓變換的瞬時頻帶識別方法。該方法首先通過Burg法對AR模型的參數進行估計，然后采用AR模型對目標信號進行延拓，最后對N次迭代后的估算瞬時頻率的小數部分進行兩次向上取整從而解決了端點處和非端點處的未重排點問題，同時也提高了頻帶的識別精度。通過兩組數值算例和一個時變拉索試驗對所提方法進行驗證，研究結果表明：

(1) 基于向上取整的多重同步擠壓變換算法采用Burg法對AR模型的參數進行有效估計，然后通過AR模型對目標信號進行延拓，從而有效避免了因端點效應而存在的未重排點問題。

(2) 基于向上取整的多重同步擠壓變換對估算的瞬時頻率進行兩次向上取整，明顯改善了MSST及IMSST中除端點外的未重排點問題，同時也提高了信號的時頻聚集性。

基于向上取整的多重同步擠壓變換識別時變結構瞬時頻帶

1 基本理論

1.1 信號延拓

1.2 基于向上取整的多重同步擠壓變換

2 數值算例

2.1 單分量調頻信號

2.2 多分量調頻信號

3 試驗驗證

4 結 論

4 結論