基于同步壓縮變換和局部替代數據的非平穩振動信號分解方法

胡異丁，任偉新，楊 棟，4

振動工程不可避免地涉及到響應信號的處理，通常假定信號是平穩的。然而，結構的時變和非線性以及外界激勵的非平穩特征都有可能使得響應信號呈現出非平穩性。響應中可能包含了激勵的某些非平穩特征和結構的某些固有特征，從而信號從整體看呈現非平穩性，但其中仍可能包含平穩成分信息［1］。因此數據的平穩性檢驗尤為必要。此外，找到方法從響應信號中分解出平穩成分和非平穩成分，則可為結構系統辨識更為細致的分類提供一種可能。

Xiao 和 Borgnat等［2－5］建立了在工作(Operational)意義下平穩性的研究理論框架，指出所謂平穩性應該相對于一定的持續時間內，及給定的觀測尺度而言的，且可以由時頻譜是否存在進化(Evolutionary)特性來決定，其概念的適用范圍包括隨機信號以及確定性的信號。換句話說，若在所有不同時間間隔內的局部的頻譜統計特性與整個觀測尺度內的全局頻譜特性基本一致，則可認為平穩。基于這樣的理解，利用替代數據(Surrogate)和多窗譜時頻分析，比較信號全局和局部頻譜特征來評價信號平穩性，并采用了零假設檢驗判定信號的平穩性。

小波變換由于其可變窗口，使其在低頻部分獲得較高的頻率分辨率和較低的時間分辨率，而在高頻部分獲得較高的時間分辨率和較低的頻率分辨率［6］。McCullough等［7］采用小波時間－尺度平面和替代數據評價信號的平穩性，并且提供了一種稱為局部替代數據(Local Surrogate)的方法，在時頻譜上非平穩信息存在的時間和頻率位置反映出來。由于小波在時頻上仍然比較模糊，降低了時頻可讀性。最近，Daubechies等［8－9］提出了一種基于小波變換的同步壓縮變換方法，該方法在小波尺度方向上利用一種特別的重分配方法，將時間－尺度平面轉化為時間－頻率平面，同時提高時頻聚集性，消除干擾項，從而進一步能洞察非平穩信號內部組成成分。

本文基于同步壓縮變換和替代數據法實現了信號平穩性的檢測，并且結合局部替代數據法提出了一種信號的平穩成分和非平穩成分分解的方法。數值算例說明該方法的有效性。

1 同步壓縮變換

高質量的時頻表示對于信號的可靠分析至關重要。信號可以看作由多個調制(調幅和調頻)分量的疊加，小波變換作為一種時頻工具廣泛應用于提取脊線來獲得這些分量的瞬時特征(如振動頻率)。然而小波變換仍然比較模糊，降低了脊線的可讀性。Daubechies提出的同步壓縮變換大大增強了時頻分辨率，使得在時頻譜上獲取脊線及其分量成為了可能。

信號的同步壓縮變換是以小波變換為基礎。首先給定小波母函數ψ(t)，對信號x(t)進行連續小波變換(CWT):

其中:a為尺度因子，b為平移因子，ψ(t )為 ψ (t)的共軛。

對任意信號x(t)，在不至于混淆的情況下記ωx(a，b)為 ω(a，b)，這樣建立了(a，b)→(ω(a，b)，b)的映射關系，同步壓縮變換再對時間－尺度平面的能量進行重分配將其轉化為時間－頻率平面，其離散計算式為:

其中:ωl為第l個離散的角頻率，ak為第k個離散的尺度點，Δω = ωl－ ωl－1，(Δak)=ak－ ak－1。最后瞬時頻率可由f=ω/2π得到。

同步壓縮變換是可逆的，對于離散計算，可以通過式(4)重構原始信號:

Thakur等［10］利用同步壓縮變換的時頻譜提取了脊線，成功分解了多分量非平穩信號，并應用于分析古氣候學數據，證明了同步壓縮變換為時頻分析的有力工具，并提供了MATLAB工具箱［11］。

2 信號平穩性檢測——時頻替代數據法

輪次檢驗和逆排序檢驗是時間序列平穩性的檢驗常用的方法［12］，但是對于具有時變頻率的非平穩信號效果往往不佳［13］。信號平穩性利用替代數據法在可以在時頻上來檢驗，其基本思想是對于給定的觀測數據樣本，在時頻譜上比較該觀測數據的全局譜和局部譜特征，若頻譜特征選在整個觀測時間內有較大變化，表明信號具有非平穩性，反之則是平穩的［5］。

替代數據最初是用于檢驗時間序列非線性的一種檢驗統計學方法［14］。根據Wiener-Khintchin定理，平穩隨機過程的功率譜密度就是信號自相關函數的傅里葉變換［12］。由于研究對象是過程的一個觀測樣本，功率譜可看作該樣本的傅里葉變換模平方除以觀測時間長度。設觀測數據x(t)的傅里葉變換為:

保持其傅里葉變換的幅度譜不變，而相位譜φf變為均勻分布在［－π，π］的隨機相位。利用傅里葉反變換可以得到時域的替代數據，即:

由于替代數據的相位譜是均勻分布在［－π，π］的隨機相位，通過這種方法產生的替代數據，具備原數據相同的功率譜幅度以及一階矩二階矩的統計特性，是廣義平穩的［3］。

給定觀測信號x(n)，通過上述方法產生一組替代數據xj'(n)(j=1～J)，對每個替代數據進行時頻變換得到TFj(f，i)(j=1～J)，其中 f 為頻率，i為時間。替代數據的全局譜分布可通過其時頻譜邊際化來獲得，局部譜為時頻譜上某個時刻的頻率分布。選擇合適的距離測度來描述和全局譜分布和某個時刻局部譜分布的差別。則對于第j個替代數據，局部譜和全局譜的距離隨時間變化的波動可用其經驗方差來描述:

其中:Di為全局譜和第i個時刻局部譜的距離，為所有時刻該距離測度的均值，N為時間點數，θ0(j)表示第j個替代數據的距離分布的方差。

同樣可求出原數據局部譜和全局譜的距離隨時間變化的波動方差θ1:

利用伽馬分布擬合θ0(j)，置信度設為95%，即可確定閾值γ，從而實現平穩性的零假設檢驗:

設平穩性零假設被拒絕，可定義平穩度指標為:

若信號是平穩的，則INS值在1附近，信號不平穩程度越大，INS的值也越大。

在同步壓縮時頻上應用替代數據法檢測數據平穩性。選擇歐氏距離作為評價局部譜分布和全局譜分布的距離測度。分析數據為汶川地震波，幅度縮小了100倍，采樣頻率為50 Hz，分析時長100 s。生成替代數據J=50次。分析結果如圖1所示。其中(a)為原始數據和一個替代數據，(b)為原始數據和替代數據經同步壓縮變換的時頻圖。(c)為檢測結果，圖中垂直虛線表示閾值 γ，垂直實線為 θ1，θ1＞γ表明 θ1落在 θ0為平穩的概率密度區域之外，因此為檢測結果等于1，為非平穩，非平穩指標INS=1.333 9。

圖1 同步壓縮時頻替代數據平穩性檢測Fig.1 Stationarity test with synchrosqueezing transformation and surrogate

3 局部替代數據

McCullough等［7］采用小波時頻和替代數據評價信號的平穩性，并且提供了一種稱為局部替代數據的方法。首先對一非平穩信號x(t)生成一批替代數據，并分別做小波變換得到其小波系數，由于替代數據是平穩的，因此不包含原信號中有意義的非平穩信息，但是可能包含了由噪聲引起的非平穩。對這些小波系數做平均:

其中:N為生成的替代數據的個數，Wi(a，b)為第i個替代數據的小波系數。在每個尺度及時間位置的標準差為:

設定閾值

其中λ為一調整參數，可依據噪聲情況而設定。將原信號x(t)的時頻Wx(a，b)與該閾值比較:

這樣，在時頻圖上原信號的平穩成分將被濾除，同時由噪聲或者泄露引起的非平穩信息也將被濾除，而保留部分就是原信號中有意義的非平穩成分，并且從時頻圖上可以清楚的反映出非平穩時間和頻率信息。但是小波時頻特性清晰度不如同步壓縮變換，且方法未能恢復出時域信號。

4 振動信號平穩成分與非平穩成分的分解

利用同步壓縮變換在時頻分析的優勢，結合局部替代數據方法，提供了一種可以將原振動信號分解為平穩成分和非平穩成分的途徑。

首先將原信號x(t)利用式(6)生成N個替代數據si(t)，這些替代數據具備原信號相同的邊際譜，并且是平穩過程，因此可利用生成的N個替代數據在時頻上的統計規律作為平穩與非平穩成分的判據，具體操作類似于前述的局部替代數據法，這里的時頻分析則采用了更為精細的同步壓縮變換。具體步驟如下。

對每個替代數據進行同步壓縮變換Ti(f，b)，對其幅度譜求均值得:

求幅度譜的標準差:

設定平穩性檢驗閾值:

其中λ為閾值調整參數。

對原信號x(t)進行同步壓縮變換并取幅度譜為Tx(f，b )，為了能濾除非平穩成分的時頻信息而保留平穩成分時頻信息，因此采用了與公式(14)相反的平穩性檢驗閾值判定條件，如式(18):

需注意，采用式(18)的判定條件，則平穩性檢驗閾值調整參數λ越小則濾除非平穩越多，反之則濾除非平穩越少。

由于得到的只是幅度譜，為了恢復時域信號，還必須保證留下的平穩成分的時頻譜上相位與原來信號一致，因此計算原信號同步壓縮變換的相位:

將保留的平穩成分的幅度譜和求出的相位譜合成為平穩成分的同步壓縮譜:

最后利用式(4)同步壓縮反變換恢復出時域的平穩成分的估計x^stationary(t)，采用同步壓縮替代數據檢測該成分是否平穩，如果不平穩則重新調整參數λ設定閾值，再重復上述過程，直至最后的平穩成分的估計x^stationary(t)滿足平穩性檢測。將原信號減去該平穩成分即得非平穩成分x^nonstationary(t)=x(t)－x^stationary(t)。算法流程如圖2所示。

圖2 算法流程圖Fig.2 Algorithm flow chart

5 數值算例

算例1 分析對象為時長40 s的信號，采樣頻率50 Hz，仿真信號包含的平穩成分為40 s幅度為1、頻率3 Hz正弦信號，非平穩成分為一段10 s幅度為1、頻率為10 Hz的正弦信號，疊加一個一段10 s幅度為1的調頻(FM)信號。分解結果如圖3所示，其中計算出λ=2.55。其中(a)、(b)、(c)分別為原始信號、分解出的平穩成分、分解出的非平穩成分的同步壓縮變換的時頻圖，(d)為它們時域波形比較，從上到下依次為原信號，平穩信號成分，非平穩信號成分。可見，采用本文提出的方法可將信號平穩及非平穩兩部分較好的分解出來，證明了該方法對調頻以及幅度變換明顯的非平穩信號分解的有效性。

圖3 正弦信號與含有明顯調幅、調頻信號的分解Fig.3 Decomposition of a sine signal plus a signal with apparent AM and FM

算例2 前面已經檢測汶川地震波是非平穩信號，因此用該信號疊加兩個幅度為1、頻率分別為3 Hz和10 Hz的正弦信號，時長100 s，采樣50 Hz。分解結果如圖4所示，取λ=1.7。其中(a)為原始信號同步壓縮變換的時頻圖，(b)圖為分解出的平穩成分時頻圖，可以看出其時頻譜基本不隨時間變化，(c)為分解出的非平穩成分時頻圖，可看出頻率隨時間變化較大，(d)為它們時域波形比較，從上到下依次為原信號，平穩信號成分，非平穩信號成分。

圖4 兩個諧波成分與汶川地震波的分解Fig.4 Decomposition of two harmonic components plus Wenchuan seismic wave

6 結論

本文提出了基于同步壓縮變換所獲得的清晰時頻圖和替代數據法相結合的信號平穩性檢測方法。提出了基于同步壓縮變換時頻圖和局部替代數據法相結合的非平穩信號分解方法。該方法通過在時頻面上分析時頻幅度譜的統計特性，利用平穩性檢測方法確定的分解閾值，在時頻圖上對信號的平穩和非平穩成分進行分解;在保持原信號相位不變的基礎上，采用同步壓縮反變換重構出時域信號的平穩和非平穩成分。數值算例驗證了所提出的信號平穩性檢測方法和非平穩信號分解方法的有效性。

［1］ Bünau P，Meinecke F C，Király F C，et al.Finding stationary subspaces in multivariate time series［J］.Physical Review Letters，2009，103(21):21401.

［2］ Xiao J，Borgnat P，Flandrin P.Testing stationarity with timefrequency surrogates［C］.15th European Signal Processing Conference［C］.Poznan，Poland，2007:2020－2024.

［3］ Xiao J，Borgnat P，Flandrin P，et al.Testing stationarity with surrogates-A one-class SVM approach［A］.IEEE Stat.Sig.Proc.Workshop SSP07［C］.Madison，WI，2007.

［4］Borgnat P， Flandrin P. Time-frequency surrogates for nonstationary signal analysis［A］.8th IMA International Conference on Mathematics in Signal Processing［C］.Cirencester，UK，2008.

［5］Borgnat P，Flandrin P，Paul H，et al.Testing stationarity with surrogates:A time-frequency approach［J］.IEEE Transactions on Signal Processing，2010，58(7):3459－3470.

［6］任偉新，韓建剛，孫增壽.小波分析在土木工程結構中的應用［M］.北京:中國鐵道出版社，2006.

［7］McCullough M，Kareem A.Testing stationarity with wavelet based surrogates［J］.Journal of Engineering Mechanics，2012，139(2):200 －209.

［8］ Daubechies I， Maes S. A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models.in:A.Aldroubi，M.Unser(Eds.)，Wavelets in Medicine and Biology［M］.CRC Press，1996:527 546.

［9］ Daubechies I，Lu J F，Wu H T.Synchrosqueezed wavelet transforms:An empirical mode decomposition-like tool［J］.Applied and Computational Harmonic Analysis， 2011，30(2):243－261.

［10］ Thakur G，Brevdo E，Fucˇkar N S，et al.The synchrosqueezing algorithm for time-varying spectral analysis:robustness properties and new paleoclimate applications［J］.Signal Processing，2013，93(5):1079－1094.

［11］ Brevdo E，Wu H T.The synchrosqueezing toolbox［EB/OL］.https://web.math.princeton.edu/～ ebrevdo/synsq/，2012，6.

［12］ Bendat J S，Piersol A G.Random data:analysis and measurement procedures，4 Ed.［M］.John Wiley and Sons，Hoboken，NJ.，2010.

［13］ Beck T W，Housh T J，Weir J P，et al.An examination of the runs test，reverse arrangements test，and modified reverse arrangements test for assessing surface EMGsignal stationarity［J］.Journal of Neuroscience Methods，2006，156(1 －2):242－248.

［14］ Theiler J， Eubank S， Longtin A， et al. Testing for nonlinearity in time series:the method of surrogate data［J］.Physica D，1992，58:77－94.