基于變參數域和短時高斯線性調頻基的自適應信號分解算法

郭劍峰，劉金朝，王衛東

（1．中國鐵道科學研究院研究生部，北京 100081；2．中國鐵道科學研究院基礎設施檢測研究所，北京 100081）

基于變參數域和短時高斯線性調頻基的自適應信號分解算法

郭劍峰1，2，劉金朝2，王衛東2

（1．中國鐵道科學研究院研究生部，北京 100081；2．中國鐵道科學研究院基礎設施檢測研究所，北京 100081）

基于高斯線性調頻基的參數化時頻分析方法由于具有很高的時域和頻域分辨能力，而被廣泛應用于非線性非穩態信號的分解和特征提取中，但其巨大的計算量常常讓工程人員望而生畏。因此結合變參數域和短時傅里葉變換的方法提出了一種改進的短時高斯線性調頻基自適應信號分解算法，將四參數優化問題轉化成窄帶范圍的兩參數優化問題，提高了參數化時頻分析的時效性。利用改進算法對四原子組合的非線性解析信號和動檢列車軸箱振動加速度信號進行分解，結果表明該方法能有效消除交叉項干擾，時頻分辨率高，而且具有計算量小，速度快的優點，對分析動檢列車軸箱振動與輪軌短波沖擊有實際意義。

變參數域；短時傅里葉變換；高斯線性調頻基；自適應分解；軸箱振動

時頻分析法［1］能同時展示信號的時間和頻率特性，隨著時頻分析方法的出現、發展和成熟，人們不斷地利用時頻分析方法，如短時傅里葉變換［2］、Wigner-Ville分布［3］、小波變換［4］、Hilbert Huang Transform（HHT）［5］等分析非平穩非線性信號的頻率和幅值特性。但上述時頻分析方法均有一定局限性。短時傅里葉變換受到窗長限制，不可能得到任意時間的時間分辨率和頻率分辨率，如：長窗短時傅里葉變換能很好地分辨振動信號中的穩態成分，但不能很好地分辨瞬態高頻成分；Wigner-Ville分布的交叉項一直苦惱著工程人員，雖然有些改進方法，如Choi-Williams分布［6］，但這些方法在降低交叉項的同時也降低了分辨率；小波變換受所謂的Q品質數的制約。直到參數化時頻分析誕生后，才進入到高分辨率、無交叉項自適應時頻分析階段。

Qian等［7］提出的基于自適應高斯基的信號時頻聯合表示開啟了參數化時頻分析的大門。通過調整高斯基函數的三個參數：時間尺度因子、時域中心和頻域中心達到基與待分解信號的最佳匹配，Mallat等［8］提出的時頻字典匹配追蹤算法的延伸，但采用全局搜索法進行參數估計，計算量大。

Mann等［9］提出了高斯包絡線性調頻基變換并稱之為Chirplet變換，這種變換的基函數與小波基相比，其時寬、帶寬和時頻中心可通過改變基函數的四個參數進行調整，并在時頻面上構成一條調頻直線，能更好的逼近信號。

Yin等［10］提出了一種高斯線性調頻基的自適應信號分解的改進快速算法，將高斯線性調頻基函數的四個參數估計問題轉換成求解超越方程組的代數根。與全局搜索算法相比，其計算速度有很大改善，但仍然難以滿足工程應用的要求。

呂貴州等［11］提出了一種基于優化初值的高斯線性調頻基自適應信號分解算法；尚朝軒等［12］提出了一種基于短時高斯線性調頻基自適應信號分解算法，這兩種算法對信號進行STFT變換，用網格劃分時頻面找到能量峰值，并以峰值點作為時頻中心，通過3 db帶寬限定時間尺度因子搜索范圍，建立高斯線性調頻基分解信號，取得了更好的效果，但調頻率的搜索范圍仍不夠明確。

彭志科等［13］提出了一種新的估計多項式相位信號瞬時頻率的參數化時頻分析方法，通過多項式非線性核函數取代線性調頻小波變換中的線性核函數，實現了單分量多項式相位信號瞬時頻率和參量的精確估計，但對多分量多項式相位信號的分析有待于進一步研究。

針對上述問題，本文提出基于變參數域和短時高斯線性調頻基的自適應信號分解算法（CPST-AGCD），通過轉變參數域的方法，把高斯基函數的兩個參數：時間尺度因子和調頻率轉變到時頻域中估計。通過對信號加長度逐漸變化的矩形窗觀察被截取的信號片段在時頻域中時寬的變化，得到穩定的時寬和帶寬，并估計時間尺度因子和調頻率的搜索范圍，再通過短時傅里葉變換準確定位時頻中心，提高了參數估計的準確性和分解計算速度。

1 基于高斯線性調頻基的自適應信號分解

信號分解可以概括為用已知類型的基函數的線性組合對未知類型的信號進行逼近的過程。自適應分解是指用來逼近待分解信號的基函數可以根據信號局部時頻特點，通過改變某些參數實現自適應調整，達到與待分解信號的最佳匹配。常用的基函數有傅里葉基、小波基和高斯包絡線性調頻基（chirple基），其中chirplet基由于本身具有最好的時頻分辨率而被廣泛用于信號的自適應分解中。

1.1 高斯線性調頻基

高斯線性調頻基函數的形式如下：

式中：αm為時間尺度因子；βm為調頻率；（tm，ωm）為基函數的時間、頻率中心高斯線性調頻基函數雖然在時域、頻域均不構成緊致框架，但卻有最好的時頻分辨率，即最小的時間－頻率帶寬積1／（4π），因此可以準確捕捉到信號的局部時頻特點。

1.2 正交匹配投影法

展開系數Bp定義為sp（t）和hp（t）的內積，其大小反映了信號和基函數的相似程度。sp（t）是sp－1（t）到hp－1（t）正交投影的信號殘量。

該過程見圖1。

圖1 正交匹配投影法示意圖Fig．1 Schematic diagram of orthogonal projection matchesmethod

通過式（4），式（5）可知，當每次分解內積最大時，信號殘量的能量最小，這樣能使分解過程逐漸收斂，當信號殘量的能量小到滿足一定要求時分解結束，可得到一組信號的近似表示。

1.3 已有參數估計方法存在的問題

信號分解的每一步要找滿足式（5）的一個高斯線

式（3）中sp（t）為經過p－1步分解后信號的殘量。因為高斯線性調頻基具有單位能量，所以經過p－1步分解后，有：性調頻基，從式（1）之中可以看出基的搜索問題實際上是估計四參數的組合優化問題。不同的參數估計方法中，各參數的取值范圍不同，即優化問題約束條件不同，計算優化問題所采用的方法也不相同，所以，分解的結果、計算速度也都不相同。已有參數估計算法存在如下四個問題：

（1）全局搜索計算量大，迭代次數多，搜索速度很慢；

（2）曲線擬合算法需要求解超越方程，計算較復雜，搜索速度較慢；

（3）基于優化初值的OI-AGCD分解算法雖然提高了分解速度，但分解存在偽能量峰值信號效果尚待改善；

（4）短時高斯包絡線性調頻基的ST-AGCD分解算法對調頻率的估計方法仍有待于研究。

2 基于變參數域和短時高斯線性調頻基的自適應信號分解算法

針對上述算法存在的問題，本文從變參數域的角度研究各參數取值范圍，使分解更準確，計算更快速。

2.1 分解算法

基于變參數域和短時高斯線性調頻基的自適應信號分解算法的步驟可描述如下：

步驟（1） m＝0，將信號s（t）賦給sm（t）準備分解。

步驟（2） 對信號sm（t）做STFT，得到能量峰值點（tm，ωm）。此時可能搜索到多個能量相等的峰值點，任選一個峰值點作為信號的能量峰值點。

步驟（3） 以該峰值點作為時頻中心，用小的初始矩形窗在時頻面截取STFT片段，考慮到STFT時頻分析方法特點，初始矩形窗的長度一般不小于33個數據采集點。

步驟（4） 仍以該峰值點作為時頻中心，加寬矩形窗，重新截取STFT片段，每次加寬的寬度可以取初始矩形窗寬度或初始矩形窗寬度的一半。當矩形窗截取的信號到達信號長度邊界時停止（見圖2）。

圖2 漸變矩形窗截取信號Fig．2 Changing rectangular window gradually to intercept signal

步驟（5） 對第j個被截取的片段，用式（6），式（7）求其時間帶寬Tj和頻率帶寬Bj，j＝0，1，2…；

步驟（6） 觀察時間帶寬Tj的變化，當時間帶寬在某個區間內穩定時，得到Tstable和Bstable。通過穩定的時間帶寬和頻率帶寬構成的Heisenberg盒計算出穩定的αstable和βstable的值如下：

式中：G為穩定的時間帶寬Tstable和時間尺度因子αstable間的函數，H為穩定的時間帶寬Tstable和穩定的頻率帶寬Bstable與調頻率βstable間的函數，求解方法將在下一節中詳細推導。

步驟（7） 以該次搜索的能量峰值點（tm，ωm）作為時頻中心，在式（8）求出的αstable和βstable附近5%的搜索范圍內生成K個基函數hmk（αmk，tm，ωm，βmk），k∈1，2…，K，計算內積，將使內積最大的一組參數作為該步分解最佳匹配的高斯線性調頻基hm（αm，tm，ωm，βm）的參數。

步驟（8） 從信號sm（t）中除去sm（t）在hm（t）的投影，得到殘余信號sm+1（t），計算殘余信號的能量，滿足要求時停止分解，不滿足時重復步驟（2）～步驟（8）繼續分解。

2.2 時間尺度因子和調頻率確定方法

通過變參數域的方法估計參數，關鍵是找到穩定的時間帶寬Tstable和基函數的時間尺度因子astable之間的函數關系。考慮高斯函數的幅值：

在時間tm取最大值，時域波形見圖3。

圖3 高斯基函數時域波形圖Fig．3 Gaussian function's domain waveform

對于被截取的信號，其時域峰值出現在時間中心點上，設在該點的函數值為emax。因為候選高斯線性調頻基的時頻中心等于信號的時頻中心，令基函數的時間帶寬與信號的時間帶寬相等，從而建立起astable與Tstable間的函數關系：αstable＝G（Tstable）。解析上，二者間的函數關系難以描述，無法直接求出G的解析形式。但根據高斯基函數的性質，其95%的能量集中在時間帶寬內，令e＝0．05*emax，t－tm＝Tstable，式（9）變為：

確定了穩定的時間帶寬 Tstable和時間尺度因子αstable間的函數G后，繼續求穩定的時間帶寬Tstable和穩定的頻率帶寬Bstable與調頻率βstable間的函數H。考慮高斯線性調頻基的瞬時頻率如下：

瞬時頻率在時頻面上構成一條斜率為βm的調頻直線（見圖4）。

圖4 瞬時頻率與調頻直線Fig．4 Instantaneous frequency and linear frequencymodulation rate

因此瞬時頻率取值與頻率帶寬和時間帶寬的比值有關。穩定的時間帶寬Tstable和穩定的頻率帶寬Bstable與調頻率βstable間的函數H為：

3 算例分析

3.1 非平穩信號分解實例

通過如下四個高斯線性調頻信號疊加成的非平穩信號的分解驗證算法性能：

s（t）＝s1（t）+s2（t）+s3（t）+s4（t）

首先，信號的Wigner-Ville分布見圖5，由圖5可知，四組chirplet信號線性疊加后的Wigner-Ville分布存在著明顯的交叉項。

圖5 原始信號的Wigner-Ville分布Fig．5Wigner-Ville distribution of original signal

對該信號做STFT變換，信號的時域、頻域、STFT時頻分布分別如圖6～圖8所示。

圖6 原始信號的時域波形圖Fig．6 Time domain waveform of original signal

圖7 原始信號的頻域波形圖Fig．7 Frequency domain waveform of original signal

圖8 原始信號的STFT時頻分布圖Fig．8 STFT waveform of original signal

進行第一步分解時信號的時頻中心為：（384，1．256 6），以該點作為時頻中心，對信號進行局部時頻分析，取初始矩形窗為33點，每次將矩形窗加寬16點截取信號，所得被截取片段的時間帶寬變化趨勢（見圖9）。

圖9 信號局部時頻特點之帶寬變化趨勢Fig．9 Local time-frequency characteristics of the signal

從圖9可知，第4次加寬矩形窗后（增加至97點后），信號時間帶寬趨于穩定，根據式（14），在第4次～第14次加寬矩形窗時選擇時間尺度因子αm和調頻率βm的搜索范圍。確定了搜索范圍后，以STFT找到的時頻中心作為基函數的時頻中心，在該范圍內搜索基函數，并滿足式（5）的最大內積條件，則可確定每一步分解的高斯線性調頻基函數。經過15步分解后，原始信號和分解后恢復出的信號見圖10。

圖10 自適應分解恢復的信號與原始信號對比Fig．10 Contrast of original signal and recovery signal

圖11 CPST-AGCD算法恢復信號的時頻分布Fig．11 Time-frequency distribution of recovery signal

結果表明，經15步分解后，原始信號和恢復出的信號之間的相關系數達到0．987 9，且信號殘量的能量不足原始信號的2．44%，整個分解過程僅需3 s即可完成，由此可見，基于信號變參數域和短時高斯線性調頻基信號分解算法有效。使用本算法恢復的信號的時頻分布見圖11。

將圖11與圖5的時頻分布比較，可知用CPSTAGCD算法對信號自適應分解后的時頻分布克服了Wigner-Ville分布中的交叉項，更清晰的揭示了信號的時頻特性。

3.2 算法復雜度分析

目前的幾種高斯線性調頻基信號分解算法中，短時高斯包絡線性調頻基自適應信號分解算法（STAGCD）效果較好。

該方法首先用短時矩形窗截取信號，并移動矩形窗分段分解。在每個矩形窗內對被截取信號做STFT，劃分時頻網格找到能量峰值，提高了時頻中心定位的準確性。此外，通過控制采樣基時寬（3db帶寬范圍）獲取有效的時間尺度因子的取值范圍，分解結果較好。本節將對兩種算法進行比較分析。

使用ST-AGCD算法對給定算例分解15步結果見圖12，與圖10比較可以看出在細節部分ST-AGCD分解效果欠佳。

圖12 ST-AGCD算法分解后恢復出的信號Fig．12 CPST-AGCD compared with ST-AGCD algorithm

對給定算例，用本文CPST-AGCD算法分解15步得到的中間過程參數見表1。

表1 CPST-AGCD分解15步中間結果Tab.1 CPST-AGCD 15 steps decomposition detailed results

比較兩種算法，計算時間和恢復信號與原始信號相似性比較結果見圖13、圖14。

圖13 兩種算法時間復雜度比較Fig．13 Comparison of two algorithms'calculating time

圖14 兩種算法分解恢復結果相關系數比較Fig．14 Comparison of two algorithms'correlation coefficient

圖13、圖14中實線為CPST-AGCD算法結果，藍色虛線為ST-AGCD算法結果。可以看出，經過15步分解后，用ST-AGCD算法恢復出的信號與原始信號之間的相似程度為0．943 2，略小于CPST-AGCD算法的結果0．987 9。CPST-AGCD自適應信號分解算法可以有效的對信號進行稀疏分解。

此外，ST-AGCD算法在時頻面上劃分網格以便準確找到能量峰值點的位置，因此在計算中需要更多的時間。

4 動檢列車軸箱振動信號分析

高速綜合檢測列車（動檢列車）是鐵路基礎設施綜合檢測的重要技術裝備，為高速鐵路運營安全評估和指導各鐵路局的養護維修提供技術支撐。動檢列車的軸箱直接與車軸相連，其振動狀態可以反映出軌道的短波不平順對車輪的沖擊作用，所以在動檢列車軸箱上安裝了加速度傳感器監測軸箱的垂向和橫向振動，安裝位置見圖15。

2013年6月，盤營高速鐵路鐵聯調聯試時，動檢列車行駛至盤營下行盤錦站里程K29+500至K29+600道岔處的速度里程圖及1車軸箱垂向振動加速度時域波形見圖16。

圖15 動檢列車軸箱加速度傳感器安裝圖Fig．15 Acceleration sensor installed on axle box

圖16 動檢列車軸箱加速度和速度里程波形圖Fig．16Waveforms of axle box acceleration and speed-mileage

圖17 STFT振動信號分析結果Fig．17 Analysis result by STFT

使用短時傅里葉變換對列車軸箱振動信號進行時頻分析的結果見圖17，可以看出STFT的時頻分辨率較低，沒有很好的時頻聚集性。

使用本文提出的CPST-AGCD對列車軸箱振動信號分解50步后的結果見圖18，可以觀察出列車在行駛到K29+540至K29+550受到道岔沖擊作用時產生頻率為378．4 Hz至391．8 Hz的持續振動。

圖18 CPST-AGCD振動信號分析結果Fig．18 Analysis result by CPST-AGCD

5 結 論

基于變參數域和短時高斯線性調頻基自適應信號分解算法根據待分解信號局部時頻特性，從信號的時間帶寬、頻率帶寬的變化規律中得到時間尺度因子、調頻率的搜索范圍，為chirplet自適應信號分解技術提出了一種新的參數估計方法，該方法有效提高了分解的快速性和自適應性，并可用于分析動檢列車的軸箱振動與輪軌短波沖擊作用等工程應用中。

［1］Abed M，Belouchrani A，Cheriet M，et al．Time-frequency distributions based on compact support kernels：properties and performance evaluation［J］．Signal Processing，IEEE Transactions on，2012，60（6）：2814－2827．

［2］Allen J．Short time spectral analysis，synthesis and modification by discrete fourier transform［J］．IEEE Transactions on Acoustics，Speech，and Signal Processing，1977，25（3）：235－238．

［3］Boashash B，O'shea P．Use of the cross Wigner-Ville distribution for estimation of instantaneous frequency［J］．Signal Processing，IEEE Transactions on，1993，41（3）：1439－1445．

［4］彭志科，何永勇，褚福磊，等．小波尺度譜在振動信號分析中的應用研究［J］．機械工程學報，2002，38（3）：122－126．

PENG Zhi-ke，HE Yong-yong，CHU Fu-lei，et al．Using wavelet scalogram for vibration signals anaylsis［J］．Chinese Journal of Mechanical Engineering，2002，38（3）：122－126．

［5］Huang N E，Wu Z．A review on Hilbert-Huang transform：Method and its applications to geophysical studies［J］．Reviews of Geophysics，2008，46（2）：RG2006．

［6］Melia U S P，VallverdúM，Claria F，et al．Choi-williams distribution to describe coding and non-coding regions in primary transcript pre-mRNA［J］．Journal of Medical and Biological Engineering，2013，33（5）：504－512．

［7］Qian S，Chen D．Signal representation using adaptive normalized gaussian functions［J］．Signal Processing，1994，36（1）：1－11．

［8］Mallat SG，Zhang Z．Matching pursuits with time-frequency dictionaries［J］．IEEE Transactions on Signal Processing，1993，41（12）：3397－3415．

［9］Mann S，Haykin S．The chirplet transform：Physical considerations［J］．IEEE Trans Signal Processing．1995，43（11）：2745－2761．

［10］Yin Q，Qian S，Feng A．A fast refinement for adaptive gaussian chirplet decomposition［J］．Signal Processing，IEEE Transactions on，2002，50（6）：1298－1306．

［11］呂貴洲，何強，魏震生，等．基于優化初值選擇的自適應高斯包絡線性調頻基信號分解［J］．信號處理，2006，22（4）：506－510．

LüGui-zhou，HE Qiang，WEI Zhen-sheng．Optimized initial-adaptive gaussian chirplet signal decomposition［J］．Signal Processing，2006，4：506－510．

［12］尚朝軒，羅賢全，何強，等．短時高斯包絡線性調頻基自適應信號分解算法［J］．信號處理，2008，24（6）：917－922．SHANG Chao-xuan，LUO Xian-quan，HE Qiang，et al．Short time-adaptive gaussian chirplet signal decomposition［J］．Signal Processing，2008，24（6）：917－922．

［13］方楊，彭志科，孟光，等．一種新的估計多項式相位信號瞬時頻率的參數化時頻分析方法［J］．噪聲與振動控制，2012，32（3）：7－11．

FANG Yang，PENG Zhi-ke，MENG Guang，et al．A new parametric time-frequency analysis method for instantaneous frequency estimation of polynomial phase signal［J］．Noise and Vibraiton Control，2012，32（3）：7－11．

Variable parameters domain and short time adaptive gaussian chirplet signal decomposition algorithm

GUO Jian-feng1，2，LIU Jin-zhao2，WANGWei-dong2
（1．Postgraduate Department，China Academy of Railway Sciences，Beijing 100081，China；2．Infrastructure Inspection Research Institute，China Academy of Railway Sciences，Beijing 100081，China）

The parametric time-frequency analysismethod based on Gaussian chirplet function has the best timefrequency resolution，so，it is widely used in non-linear and non-stationary signal decomposition and feature extraction．But it needs a large amount of computation．A reformed short time Gaussian chirplet signal decomposition algorithm based on variable parameters domain method and short time Fourier transform（STFT）was proposed．Taking as an example，it tranfers a four parameters optimization problem to two parameters one in a narrow range and improve the efficiency of computation．The reformed algorithm was used to decompose a four atoms non-linear analytic signal and the vibration accelation signal of a high speed comprehensive inspection train's axle box．The results show the algorithm can avoid the cross-term's interferer and achieve fast computation．It can be applied to analyze the vibration of axle box and the wheelrail shortwave shock．

variable parameters domain；short time Fourier transform；Gaussian chirplet function；adaptive decomposition；axle box vibration

TN911．72

A

10．13465／j．cnki．jvs．2015．12．023

國家973重點基礎研究發展計劃項目（2013CB329406）；國家自然科學基金資助項目（51178464）；中國鐵道科學研究院基金項目（2013YJ068，2013YJ069）

2014－04－21 修改稿收到日期：2014－06－19

郭劍峰男，博士生，1987年8月生

王衛東 男，研究員，博士生導師，1963年生