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

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

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

引言

環(huán)境激勵(lì)下的在役土木工程結(jié)構(gòu)本質(zhì)上是時(shí)變和非線性結(jié)構(gòu)，其響應(yīng)信號呈現(xiàn)非平穩(wěn)性。時(shí)頻分析方法（time frequency analysis，簡稱TFA）能同時(shí)在時(shí)域和頻域提取信號瞬時(shí)特征，是分析非平穩(wěn)信號強(qiáng)有力的工具［1‐2］。其主要方法有Wigner‐Ville 分布［3］、短時(shí)傅里葉變換（short time Fourier trans‐form，簡稱STFT）［4］和小波變換［5］等。

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

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

1 基本理論

1.1 同步提取變換

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

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

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

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

設(shè)定改進(jìn)高斯窗φ(ω，t)=

則式（3）可改寫為

根據(jù)傅里葉變換性質(zhì)及Parseval 定理，式（4）可表示為

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

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

對式（7）求偏導(dǎo)可得

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

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

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

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

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

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

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

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

1.2 基于同步提取和時(shí)頻系數(shù)模極大值的聯(lián)合方法提取時(shí)頻脊線和瞬時(shí)頻率

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

1）根據(jù)式（7）對信號進(jìn)行STFT 得到時(shí)頻系數(shù)矩陣；

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

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

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

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

2 數(shù)值算例驗(yàn)證

2.1 單分量信號數(shù)值算例

考慮如式（14）所示的單分量調(diào)幅調(diào)頻信號

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

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

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

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

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

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

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

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

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

2.2 多分量信號數(shù)值算例

考慮如式（16）所示多分量調(diào)幅調(diào)頻信號

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

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

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

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

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

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

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

3 試驗(yàn)驗(yàn)證

筆者從剛度時(shí)變和質(zhì)量時(shí)變兩個(gè)角度來設(shè)計(jì)時(shí)變結(jié)構(gòu)動(dòng)力試驗(yàn)以進(jìn)一步驗(yàn)證所提CM 方法的有效性和準(zhǔn)確性，并將瞬時(shí)頻率識(shí)別結(jié)果與基于小波系數(shù)模極大值和同步擠壓小波變換的識(shí)別值進(jìn)行了對比。

3.1 剛度線性變化拉索試驗(yàn)

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

圖13 試驗(yàn)裝置圖Fig.13 The cable test setup

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

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

圖15 拉索響應(yīng)信號的時(shí)頻圖Fig.15 The time-frequency spectrogram of the response from the cable

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

圖17 拉索響應(yīng)信號的瞬時(shí)頻率識(shí)別結(jié)果對比Fig.17 Comparison of IF identification results of the re‐sponse from the cable

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

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

3.2 質(zhì)量突變懸臂梁試驗(yàn)

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

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

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

圖21 實(shí)測懸臂梁響應(yīng)信號Fig.21 The measured acceleration response from the cantile‐ver beam

圖22 懸臂梁響應(yīng)信號時(shí)頻圖Fig.22 The time-frequency spectrogram of the response from the cantilever beam

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

圖24 懸臂梁響應(yīng)信號的瞬時(shí)頻率識(shí)別結(jié)果對比Fig.24 Comparison of IF identification results of the re‐sponse from the cantilever beam

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

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

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

4 結(jié)論

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

2）CM 法能夠有效識(shí)別剛度和質(zhì)量時(shí)變結(jié)構(gòu)非平穩(wěn)響應(yīng)信號的瞬時(shí)頻率，且識(shí)別效果優(yōu)于小波系數(shù)模極大值和同步擠壓小波變換方法。

3）CM 法的識(shí)別結(jié)果不會(huì)隨著信號的減弱而精度下降，特別是在信號末端噪聲容易掩蓋真實(shí)信號的時(shí)候表現(xiàn)良好，因而更加適合實(shí)際響應(yīng)信號的瞬時(shí)特征參數(shù)提取。