基于同步提取和時頻系數模極值的瞬時頻率識別?

劉景良，王新宇，鄭錦仰，2，盛 葉，駱勇鵬

（1.福建農林大學交通與土木工程學院 福州，350002）（2.福建省建筑科學研究院 福州，350025）

引言

環境激勵下的在役土木工程結構本質上是時變和非線性結構，其響應信號呈現非平穩性。時頻分析方法（time frequency analysis，簡稱TFA）能同時在時域和頻域提取信號瞬時特征，是分析非平穩信號強有力的工具［1‐2］。其主要方法有Wigner‐Ville 分布［3］、短時傅里葉變換（short time Fourier trans‐form，簡稱STFT）［4］和小波變換［5］等。

采用時頻分析方法識別時變結構響應信號瞬時頻率需要解決的一個關鍵問題是脊線提取，而最常見的脊線提取方法是時頻系數模極大值法。如Liu等［6］提出基于小波系數模局部極大值的小波脊線提取方法并將之應用于干涉圖分析，該方法的主要優點是運算簡單、計算量小，但易受噪聲和端點效應的干擾［7］。為此，Carmona 等［8］提出基于隨機走動的瘋狂爬坡算法來提取信號各分量的小波脊線，但該方法的計算效率不高。Wang 等［9］提出基于動態規劃的非平穩信號瞬時頻率提取方法在一定程度上消除了端點效應，但其算法較為復雜。近年來，一些后處理方法被引入以解決噪聲等外部因素對信號時頻脊線提取的影響，如重新分配法［10‐11］、同步擠壓小波變換［12‐14］、參數化TFA［15‐17］、解調TFA［18］和同步提取變換（synchroextracting transform，簡稱SET）［19］等。其中，SET 可以實現時頻系數能量的高度聚集和響應信號的完整重構，但是目前尚未運用于土木工程領域。

通過對時變結構響應信號進行同步提取變換，瞬時頻率被鎖定在一個范圍內，有效避免了隨機噪聲對瞬時頻率的干擾。在此基礎上，筆者將同步提取與時頻系數模極大值脊線提取算法相結合，提出一種新的聯合方法（combined method，簡稱CM）來精確識別響應信號的瞬時頻率。首先，對響應信號進行STFT 得到時頻圖；其次，通過同步提取算法從時頻圖中提取一條或多條狹窄的頻帶，在降低噪聲干擾的同時也減小了時頻系數模極大值算法的計算量；最后，在限定的頻帶范圍內逐點搜索時頻系數模極大值點，并按時間順序連接起來從而得到精確的時頻脊線和瞬時頻率曲線。

1 基本理論

1.1 同步提取變換

由于諧波信號經過時頻變換得到的時頻系數實部的頻率與原信號的瞬時頻率相等，因此同步提取變換通過計算各時頻系數實部的頻率來確定信號的瞬時頻率（instantaneous frequency，簡稱IF）［19］。同步提取變換作為一種后處理技術，可以和STFT 等多種TFA 方法結合。以STFT 為例，原信號s(t)的STFT 可表示為

其中：τ為時移因子；ω為頻率；t為時間；g(τ-t)為可移動的緊支高斯窗函數。

其中：σ為決定窗寬的尺度因子。

將式（2）代入式（1）可得

設定改進高斯窗φ(ω，t)=

則式（3）可改寫為

根據傅里葉變換性質及Parseval 定理，式（4）可表示為

若以諧波信號s(t)=為例，其傅里葉變換可表示為

將式（6）代入式（5）得

對式（7）求偏導可得

至此，原信號s(t)的瞬時頻率，即式（8）中的ω0可表示為

由于同步提取變換僅保留時頻面中瞬時頻率位置的時頻系數而令其余位置的時頻系數歸零，因此可通過對時頻系數矩陣乘以平移到ω0(τ，ω)的δ(ω)函數來實現同步提取，如式（10）所示

Te(τ，ω)=SSTFT(τ，ω)δ[ω-ω0(τ，ω)] （10）

其中：δ[ω-ω0(τ，ω)]為克羅內克函數。

考慮到計算誤差的存在，將式（11）改寫為

其中：Δω為采樣頻率間隔。

將式（12）代入式（10）得到同步提取的最后結果，如式（13）所示

由于同步提取變換算法涉及到求導、相除等運算，將造成一定的計算誤差，而定義頻帶范圍就是為了更好地考慮在瞬時頻率提取過程中出現的識別誤差。

1.2 基于同步提取和時頻系數模極大值的聯合方法提取時頻脊線和瞬時頻率

響應信號經時頻變換后其時頻圖上的某些位置會出現亮點，這表明該處的時頻系數模值較大。因此，可通過搜索時頻系數模極大值位置來獲取時頻脊點，然后將脊點按照時間順序連接起來即為時頻脊線。然而，時頻系數模極大值算法由于抗噪性較差，其提取的瞬時頻率曲線含有大量毛刺，影響了識別精度?；诖?，筆者引入同步提取算法并與時頻系數模極大值法相結合來解決上述問題。該聯合算法的具體步驟如下：

1）根據式（7）對信號進行STFT 得到時頻系數矩陣；

2）按照式（9）初步估算出信號的瞬時頻率；

3）根據式（13）將時頻面上非瞬時頻率范圍位置上的時頻系數歸零；

4）令i，j分別表示時頻系數矩陣的列和行，初始時令i=1，搜索第i列在限定頻帶范圍內的時頻系數最大值，最后令j(i)=j，i=i+1，并開始搜索下一列在限定頻帶范圍內的時頻系數最大值；

5）重復步驟4，直至求出時頻系數矩陣所有列中限定頻帶范圍內的時頻系數模極大值，并將極大值點連接成線，即為待求時頻脊線和瞬時頻率曲線。

2 數值算例驗證

2.1 單分量信號數值算例

考慮如式（14）所示的單分量調幅調頻信號

根據瞬時頻率定義，x(t)的理論頻率為f=25+5cos(πt)。設定信號采樣頻率為600 Hz，采樣時長為10 s。為模擬真實情況中噪聲對信號的影響，對上述信號添加50%水平的高斯白噪聲，其噪聲強度由信噪比（signal noise ratio，簡稱SNR）定義

其中：Asignal和Anoise分別為信號和噪聲的均方根值；噪聲水平是指A2noise與A2signal之間的比值。

添加50%水平高斯白噪聲后的信號如圖1 所示。首先，對信號進行STFT 得到如圖2 所示的時頻圖；其次，對信號x(t)進行同步提取變換，得到的時頻曲帶如圖3 所示；最后，采用時頻系數模極大值法在鎖定的瞬時頻率曲帶范圍內搜索時頻系數模極大值。通過比較圖2 和圖3 可以看出，同步提取變換將瞬時頻率鎖定在一個曲帶范圍內，同時也去除了圖2 中所含的毛刺。為方便比較，圖4 同時給出了基于CM 法、小波系數模極大值法［6］和同步擠壓小波變換方法［5］的瞬時頻率識別值以及理論值。圖5 為上述3 種方法瞬時頻率識別結果的局部放大效果圖。圖6 給出了3 種方法的相對誤差。由圖4和圖5 可知：筆者提出的CM 方法識別結果與理論值吻合較好，且具有較高的抗噪性；同步擠壓小波變換的識別效果次之，但是同步擠壓小波變換是在每一中心頻率附近區間固定一個頻率值來細化時頻曲線，同時最終提取的瞬時頻率曲線也是時頻圖中能量最高或曲率最小的部分；該方法雖然提高了頻率分辨率，但也降低了時間分辨率［20］；小波系數模極大值法識別的瞬時頻率曲線因含有大量毛刺而與理論頻率值存在較大的偏差，整體識別效果較差。

圖1 添加50%水平高斯白噪聲的x(t)Fig.1 x(t)with 50% Gaussian white noise

圖2 x(t)的時頻圖Fig.2 The time-frequency spectrogram of x(t)

圖3 同步提取變換識別x(t)的瞬時頻率曲帶Fig.3 The IF curve band of x(t)extracted by SET

圖4 x(t)的瞬時頻率識別結果對比Fig.4 Comparison of IF identification results of x(t)

圖5 x(t)的瞬時頻率識別結果局部放大效果Fig.5 Local magnification of IF identification results of x(t)

圖6 x(t)的瞬時頻率識別相對誤差Fig.6 Relative errors of IF identification of x(t)

2.2 多分量信號數值算例

考慮如式（16）所示多分量調幅調頻信號

設定采樣頻率為600 Hz，采樣時間為10 s。瞬時頻率理論值分別為f1=10+2.5 cos(0.5πt)和f2=30+2.8 cos(1.4πt)。對信號添加50%水平的高斯白噪聲，其含噪信號如圖7 所示。對信號進行STFT 得到如圖8 所示的時頻圖。利用同步提取變換去除噪聲干擾并鎖定瞬時頻率的限定范圍，結果如圖9 所示。圖10 同時給出了基于CM、小波系數模極大值法和同步擠壓小波變換提取的瞬時頻率識別結果以及理論瞬時頻率曲線。為更清楚地比較3種方法的瞬時頻率識別效果，圖11 給出了局部放大結果，圖12 則給出了3 種方法的相對誤差。通過對比圖11（a）和圖11（b）可以看出：在識別f1時，由于f1變化緩慢，該信號為慢變信號，CM 和同步擠壓小波變換識別的瞬時頻率結果均與理論值十分吻合；在識別f2時即瞬時頻率相對快變的情況下，CM 法仍具有很好的穩定性且識別結果與理論值吻合較好，但是同步擠壓小波變換識別結果由于時間分辨率下降的原因其識別效果相對較差；小波系數模極大值法識別的瞬時頻率曲線則含有大量毛刺，識別效果不佳。總之，CM 識別的瞬時頻率曲線不僅最平滑而且基本沒有毛刺，其識別值也最接近理論值，在3 種瞬時頻率識別方法中表現最優。

圖7 添加50%水平高斯白噪聲的y(t)Fig.7 y(t)with 50% Gaussian white noise

圖8 y(t)的時頻圖Fig.8 The time-frequency spectrogram of y(t)

圖9 同步提取變換識別y(t)的瞬時頻率曲帶Fig.9 The IF curve band of y(t)extracted by SET

圖10 y(t)的瞬時頻率識別結果對比Fig.10 Comparison of IF identification results of y(t)

圖11 y(t)的瞬時頻率識別結果局部放大效果Fig.11 Local magnification of IF identification results of y(t)

圖12 y(t)的瞬時頻率識別相對誤差Fig.12 Relative errors of IF identification of y(t)

3 試驗驗證

筆者從剛度時變和質量時變兩個角度來設計時變結構動力試驗以進一步驗證所提CM 方法的有效性和準確性，并將瞬時頻率識別結果與基于小波系數模極大值和同步擠壓小波變換的識別值進行了對比。

3.1 剛度線性變化拉索試驗

剛度線性變化拉索試驗裝置如圖13 所示。試驗所用拉索為1 根4.55 m，7?5 的預應力鋼絞線拉索，一端用反力架固定，另一端固定于電液伺服加載系 統（mechanical testing &simulation，簡稱MTS）上。在拉索中部豎向安裝加速度傳感器并采集加速度沖擊響應。在初始階段對拉索施加20 kN 的預拉力后，通過調整MTS 的作動器使索的拉力以1.67 kN/s 的速率線性增加，拉力的變化將引起拉索剛度的改變。在改變拉力的同時，用力錘敲擊拉索以采集拉索的加速度沖擊響應，采樣頻率設為600 Hz，采樣時間為6 s，所得響應信號如圖14 所示。

圖13 試驗裝置圖Fig.13 The cable test setup

圖14 拉索線性變化時的加速度響應Fig.14 The measured cable acceleration responses with lin‐early varying cable tension force

首先，采用解析模態分解方法［21］對實測響應信號進行分解，得到1 階分量信號并對其進行STFT得，到如圖15 所示的時頻圖；其次，對信號進行同步提取變換以去除噪聲影響并將瞬時頻率鎖定在一定范圍內，結果如圖16 所示；最后，在限定范圍內提取信號的時頻系數模極大值并識別其瞬時頻率。為方便比較，圖17 給出了CM 法、小波系數模極大值法、同步擠壓小波變換的瞬時頻率識別結果以及基于“凍結法”［22］求解的瞬時頻率理論值。圖18 和圖19則分別給出了相對誤差曲線圖和相對誤差曲線局部放大圖。由圖17 可知，CM 法識別的拉索瞬時頻率變化趨勢為線性，而且最接近理論值。在信號的末端由于隨機噪聲掩蓋了真實響應信號，此時同步擠壓小波變換容易出現端點效應，其瞬時頻率識別結果偏離理論值較大［5］；而CM 法在信號末端并沒有出現毛刺現象和明顯的端點效應，其針對實測信號的識別效果較佳。

圖15 拉索響應信號的時頻圖Fig.15 The time-frequency spectrogram of the response from the cable

圖16 同步提取變換識別拉索響應信號的瞬時頻率曲帶Fig.16 The IF curve band of the response from the cable ex‐tracted by SET

圖17 拉索響應信號的瞬時頻率識別結果對比Fig.17 Comparison of IF identification results of the re‐sponse from the cable

圖18 拉索信號的瞬時頻率識別相對誤差Fig.18 Relative errors of IF identification of the response from the cable

圖19 拉索響應信號的瞬時頻率相對誤差局部放大效果Fig.19 Local magnification of IF relative errors of the re‐sponse from the cable

3.2 質量突變懸臂梁試驗

質量突變懸臂梁試驗裝置如圖20 所示。鋁合金懸臂梁的質量為0.81 kg，截面尺寸為40 mm×15 mm，長度為500 mm。懸臂端的質量塊為1 kg，在質量塊上方放置一塊懸掛的磁鐵。在梁上每隔100 mm 設置一錘擊點，然后通過力錘敲擊懸臂梁的自由端，并在某個時刻放下細繩使永磁鐵垂直靠近質量塊并將其吸起以實現結構質量的改變。設定采樣頻率為2 kHz，在梁的跨中位置安裝加速度傳感器以采集加速度沖擊響應。試驗前，基于“凍結法”測得帶質量塊與不帶質量塊時懸臂梁的基頻分別為21.24 和47.06 Hz。試驗過程中，先在懸臂梁末端用力錘施加激振力，并在2 s 附近利用磁鐵吸引質量塊以改變懸臂梁質量，最后在2.3 s 再次錘擊懸臂梁以避免響應信號衰減過快。

圖20 鋁合金懸臂梁試驗裝置圖(單位：mm)Fig.20 The aluminum cantilever beam test rig (unit:mm)

加速度傳感器記錄的響應信號如圖21 所示。首先，采用解析模態分解定理對實測響應信號進行分解得到1 階分量信號，再對1 階分量信號進行ST‐FT 得到如圖22 所示的時頻圖，由圖可知頻率在2s附近發生突變，與實際情況相符；其次，對響應信號進行同步提取變換，得到如圖23 所示的瞬時頻率曲帶；最后，采用時頻系數模極大值法提取實測響應信號的瞬時頻率。為方便比較，圖24 同時給出了CM法、小波系數模極大值法和同步擠壓小波變換法提取的瞬時頻率曲線以及理論值。圖25 給出了瞬時頻率識別結果的局部放大圖。圖26 和27 分別給出了相對誤差曲線和局部相對誤差曲線。由圖24 和圖25 可知：小波系數模極大值法識別的瞬時頻率曲線存在大量毛刺；雖然同步擠壓小波變換識別的曲線比小波系數模極大值法平滑，但其識別效果仍然不如CM 法，且在2.0 s 附近出現了明顯波動。

圖21 實測懸臂梁響應信號Fig.21 The measured acceleration response from the cantile‐ver beam

圖22 懸臂梁響應信號時頻圖Fig.22 The time-frequency spectrogram of the response from the cantilever beam

圖23 同步提取變換識別懸臂梁響應信號的瞬時頻率曲帶Fig.23 The IF curve band of the response from the cantile‐ver beam extracted by SET

圖24 懸臂梁響應信號的瞬時頻率識別結果對比Fig.24 Comparison of IF identification results of the re‐sponse from the cantilever beam

圖25 懸臂梁響應信號的瞬時頻率識別結果局部放大效果Fig.25 Local magnification of IF identification results of the response from the cantilever beam

圖26 懸臂梁響應信號的瞬時頻率識別相對誤差Fig.26 Relative errors of IF identification of the response from the cantilever beam

圖27 懸臂梁響應信號的瞬時頻率識別相對誤差局部放大效果Fig.27 Local magnification of IF relative errors of the re‐sponse from the cantilever beam

4 結論

1）通過同步提取算法鎖定了瞬時頻率范圍，避免了時頻系數模極大值易受噪聲干擾的缺點，同時也減小了時頻系數模極大值法的搜索范圍。

2）CM 法能夠有效識別剛度和質量時變結構非平穩響應信號的瞬時頻率，且識別效果優于小波系數模極大值和同步擠壓小波變換方法。

3）CM 法的識別結果不會隨著信號的減弱而精度下降，特別是在信號末端噪聲容易掩蓋真實信號的時候表現良好，因而更加適合實際響應信號的瞬時特征參數提取。