基于同步壓縮變換和局部替代數(shù)據(jù)的非平穩(wěn)振動信號分解方法

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

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

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

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

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

1 同步壓縮變換

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

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

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

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

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

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

Thakur等［10］利用同步壓縮變換的時頻譜提取了脊線，成功分解了多分量非平穩(wěn)信號，并應(yīng)用于分析古氣候?qū)W數(shù)據(jù)，證明了同步壓縮變換為時頻分析的有力工具，并提供了MATLAB工具箱［11］。

2 信號平穩(wěn)性檢測——時頻替代數(shù)據(jù)法

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

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

保持其傅里葉變換的幅度譜不變，而相位譜φf變?yōu)榫鶆蚍植荚冢郏?，π］的隨機(jī)相位。利用傅里葉反變換可以得到時域的替代數(shù)據(jù)，即:

由于替代數(shù)據(jù)的相位譜是均勻分布在［－π，π］的隨機(jī)相位，通過這種方法產(chǎn)生的替代數(shù)據(jù)，具備原數(shù)據(jù)相同的功率譜幅度以及一階矩二階矩的統(tǒng)計特性，是廣義平穩(wěn)的［3］。

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

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

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

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

設(shè)平穩(wěn)性零假設(shè)被拒絕，可定義平穩(wěn)度指標(biāo)為:

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

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

圖1 同步壓縮時頻替代數(shù)據(jù)平穩(wěn)性檢測Fig.1 Stationarity test with synchrosqueezing transformation and surrogate

3 局部替代數(shù)據(jù)

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

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

設(shè)定閾值

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

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

4 振動信號平穩(wěn)成分與非平穩(wěn)成分的分解

利用同步壓縮變換在時頻分析的優(yōu)勢，結(jié)合局部替代數(shù)據(jù)方法，提供了一種可以將原振動信號分解為平穩(wěn)成分和非平穩(wěn)成分的途徑。

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

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

求幅度譜的標(biāo)準(zhǔn)差:

設(shè)定平穩(wěn)性檢驗閾值:

其中λ為閾值調(diào)整參數(shù)。

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

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

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

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

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

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

5 數(shù)值算例

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

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

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

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

6 結(jié)論

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

［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，F(xiàn)landrin P.Testing stationarity with timefrequency surrogates［C］.15th European Signal Processing Conference［C］.Poznan，Poland，2007:2020－2024.

［3］ Xiao J，Borgnat P，F(xiàn)landrin 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，F(xiàn)landrin 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］任偉新，韓建剛，孫增壽.小波分析在土木工程結(jié)構(gòu)中的應(yīng)用［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，F(xiàn)ucˇ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.