頻率同步壓縮與時間同步壓縮的對比和應用?

何周杰，涂曉彤，王 凱，李富才，包文杰，包 雋

（1.上海交通大學機械系統與振動國家重點實驗室 上海，200240）

（2.博世華域轉向系統有限公司 上海，201821）

引言

機械系統的故障診斷中，傳感器直接獲取的一般是時域信號，時域信號能準確地反映機械系統在不同時刻下的運行狀態(tài)。對時域信號進行傅里葉變換（Fourier transform，簡稱FT）可以得到信號的頻域表示，頻域表示刻畫了整個分析時長下信號所包含的頻率成分，但在實際情況下，單一地從時域或者頻域中分析故障特征是較困難的。為解決這一難題，人們提出一系列的時頻聯合分析方法。

時頻聯合分析方法的本質是將信號的一維時域信息轉化為二維時頻信息，從時間和頻率兩方面同時表征信號。常用的時頻聯合分析方法包括短時傅里葉變換和小波變換（continuous wavelet transform，簡稱CWT），它們都屬于線性變換，但由于測不準原理，時間和頻率不可能同時獲得較高的分辨率［1］。因此，為了提高時頻表示的可讀性，時頻重排法（reassignment method，簡稱RM）應運而生［2］，該方法通過瞬時頻率算子（instantaneous frequency operator，簡稱IFO）和群延遲算子（group delay operator，簡稱GDO）的估計，將彌散的信號能量重新映射到脊線的中心，從而達到提高時頻表示可讀性的目的。但由于RM是沿著頻率軸和時間軸方向同時壓縮信號的能量，且忽略了相位信息，因此它不具備信號重構的特性。為了解決信號的重構難題，人們提出了頻率同步壓縮（frequency-reassigned synchrosqueezing transform，簡稱SST），因其只考慮頻率軸上的壓縮，從而具備信號重構的能力［3］。最初的SST是在CWT的基礎上推導得出的，后來經過改進，得到了在STFT下的同步壓縮變換［4］。在此之后，為了進一步提高信號能量的集中程度，同步壓縮出現了不同的發(fā)展方向：通過將匹配解調與同步壓縮相結合，發(fā)展出了廣義同步壓縮變換［5］；將二階或高階多項式引入瞬時頻率模型之中，得到二階或高階同步壓縮變換［6-8］；將同步壓縮過程和脊線提取過程相結合，得到同步提取變換［1］。上述方法均是建立在SST的基礎之上提出的，由于它們的壓縮重排是沿著頻率軸方向進行的，因此很難處理信號脊線是平行于頻率軸的時頻壓縮，例如在軸承沖擊信號的時頻表示中，常常出現平行于頻率軸的周期性脊線特征。為此人們提出了時間同步壓縮，它的壓縮重排是沿著時間軸方向進行的，因此其具備處理強時變信號的能力［9］。

筆者以SST和TSST為基礎，通過對比兩者所采用的STFT來說明兩種同步壓縮的區(qū)別以及各自的應用場合，然后按照同步壓縮的流程推導出它們的實現算法，最后通過轉軸的碰摩故障和軸承外圈沖擊故障的實驗數據來驗證兩種算法的有效性。

1 兩種同步壓縮原理對比

1.1 SST簡述

為了方便，將SST所采用的STFT稱為改進短時傅里葉變換（modified short-time Fourier transform，簡稱MSTFT），它可以表示為其中：X(ω)為信號x(t)的頻域表示；G(ω)為窗函數g(t)的頻域表示；VMgx(t，ω)表示MSTFT計算結果；文中*為共軛符號。

SST是根據IFO估計值來完成MSTFT系數在頻率軸上的壓縮重排，因此IFO的計算是SST中至關重要的一步，IFO一階表達式ω?(t，ω)可以寫為

其中：ω?(t，ω)表示IFO；?(·)表示對復數取實部。

對應于MSTFT的逆變換公式為

由式（3）可以看出，MSTFT的逆變換只在頻率方向積分，因此如果MSTFT系數只在頻率方向重排，是不會影響到信號的重構，SST重排公式可寫為

因此，只需把移動后聚集在脊線附近的系數再進行一次頻率積分就可以得到原始信號，故SST的重構公式表達為

1.2 TSST簡述

與SST不同，TSST采用的STFT公式是其最初被提出時所采用的形式，在文中被稱為傳統短時傅里葉變換（traditional short-time Fourier transform，簡稱TSTFT），其可以被表示為其中TSTFT計算結果。

TSST是根據GDO估計值來完成TSTFT系數在時間軸上的壓縮重排，因此GDO的計算是TSST中最重要的一步，GDO一階表達式τ?(t，ω)為

其中：τ?(t，ω)為GDO；?(·)為復數取虛部。對應于TSTFT的逆變換公式為

其中：F-1ω(·)為變量ω的傅里葉逆變換。

由式（8）可知道，TSTFT的反變換只在時間方向積分，因此如果TSTFT系數只在時間軸方向重排，不會影響到信號的重構，TSST的重排公式可寫為

因此，只需把移動后聚集在脊線附近的系數再進行一次時間積分就可以得到原始信號的頻域表示，故SST的重構公式可以表達為

1.3 SST與TSST的 區(qū) 別

1.3.1 兩種不同的短時傅里葉變換

信號處理中復正弦信號和脈沖信號是兩種典型的信號。復正弦信號的脊線為平行于時間軸的直線，圖1為復正弦信號的STFT，SST和TSST結果，其時域表達式和頻域表達式分別為x1(t)和X1(ω)

圖1 復正弦信號的時頻表示結果Fig.1 The time-frequency representations of a complex sinusoidal signal

脈沖信號的脊線為平行于頻率軸的直線，圖2為脈沖信號的STFT、SST和TSST結果，而其時域和頻域表達式分別由x2(t)和X2(ω)表示

考察復正弦信號x1(t)，其MSTFT可以表示為

如圖1（a）所示，因為窗函數在頻域上為緊支撐函數，所以x1(t)的譜圖集中于ω=ω0的水平帶狀分布。在式（15）中ejω0t為t的函數，故的相位僅會沿著時間軸方向產生振蕩，因此對沿著頻率軸方向壓縮的SST來說，不會因為相位振蕩而出現正負STFT系數抵消的現象，如圖1（b）所示。

而x1(t)的TSTFT表示為式（16）中，ej(ω0-ω)t為ω，t的函數，故的相位沿著時間軸和頻率軸方向都將發(fā)生相位振蕩，因此對沿著時間軸方向壓縮的TSST來說，會因為相位振蕩而出現正負STFT系數抵消的現象，如圖1（c）所示，其亮帶寬度和圖1（a）的相同但能量反而降低，故不能達到壓縮脊線的目的。

接下來分析脈沖信號x2(t)，其MSTFT表示為

如圖2（a）所示，因為窗函數在時域上為緊支撐函數，所以x2(t)的譜圖是集中于t=t0的垂直帶狀分布。在式（17）中由于e-jω(t0-t)為ω，t的函數，故VMgx2(t，ω)的相位沿著時間軸和頻率軸都將產生相位振蕩，因此對沿著頻率軸方向壓縮的SST來說，會因為相位振蕩而出現正負STFT系數抵消的現象，如圖2（b）所示，其亮帶寬度和圖2（a）的相同但能量反而降低，故不能達到壓縮脊線的目的。

圖2 脈沖信號的時頻表示結果Fig.2 The time-frequency representations of a pulse signal

而x2(t)的TSTFT可以表示為式（18）中e-jωt0為ω的函數，故的相位僅會沿著頻率軸方向產生振蕩，因此對沿著時間軸方向壓縮的TSST來說，是不會因為相位振蕩而出現正負STFT系數抵消的現象，如圖2（c）所示。

1.3.2 兩種不同的壓縮重排過程

考察SST和TSST的壓縮重排過程，定義仿真測試信號x3(t)

其STFT譜圖如圖3（a）所示，用1，2分別表示圖中兩個橢圓的內部區(qū)域。其中區(qū)域1代表信號頻率緩慢變化的部分，區(qū)域2表示信號頻率急劇變化的部分；圖3（b）和圖3（c）分別是IFO平面和GDO平面，其中紅色小箭頭的指向表示兩種同步壓縮的重排方向；圖3（d）和圖3（e）分別是SST和TSST壓縮后的時頻平面。

SST壓縮重排是將圖3（a）中STFT系數，按照圖3（b）中對應點所計算出來的IFO，沿著頻率軸方向移動到新的位置，從而完成從圖3（a）～（d）的轉換。TSST壓縮重排是將圖3（a）中的STFT系數，按照圖3（c）中對應點所計算出來的GDO，沿著時間軸方向移動到新的位置，完成從圖3（a）～（e）的轉換。

由圖3（b）可以看出，在垂直方向上所有的點的IFO值在區(qū)域1中基本相同而在區(qū)域2中差異較大，故SST僅能在區(qū)域1內取得很好的壓縮效果（如圖3（d）所示）。相反，由圖3（c）可以看出，在水平方向上所有的點的GDO值在區(qū)域1中差異較大而在區(qū)域2在中基本相同，故TSST僅在區(qū)域2內可以取得很好的壓縮效果（如圖3（e）所示）。

綜上所述，SST由于采用了MSTFT公式，從而避免其相位沿著頻率軸方向的振蕩，因此在頻率方向上壓縮緩變信號能取得很好的效果，適合診斷定轉速下轉軸的碰摩類的故障；而TSST由于采用了TSTFT公式，從而避免其相位沿著時間軸方向上的振蕩，因此在時間方向上壓縮快變信號有更好的效果，適合診斷軸承沖擊類的故障。

圖3 SST和TSST對測試信號x3(t)的壓縮重排過程Fig.3 The SST and the TSST compression reassignment processes for the test signal x3(t)

2 兩種同步壓縮算法實現

不論是SST還是TSST，同步壓縮的實現都包含以下4個主要步驟：信號STFT的計算、IFO或GDO的估計、STFT系數的壓縮重排及信號的重構。

首先對信號作離散化處理，假設信號x(t)的采樣率為Fs，采樣時長為T，則采樣點數為：N=T×Fs，故 離 散 化 信 號 可 表 示 為N-1}，對x[b]作快速傅里葉變換，可得信號的N點離散頻譜：{X[m]|m∈Z+0，m≤N-1}。

2.1 MSTFT和TSTFT的 算 法 實 現

2.1.1 MSTFT的算法實現

對于SST，首先假設則MSTFT可以表示為

為了使式（21）可以寫成卷積形式，令m(t)為m(t)的共軛取反

結合式（23）和式（24），得到MSTFT計算公式

將式（25）中的頻移ω離散化成長度為L1的序列：{k|k∈Z+0，k≤L1-1}，得離散化的MSTFT

計算式（26）的復雜度，由于G(v)可通過窗函數公式得到解析表達式，故不必把求解窗函數的FT的計算復雜度計入，因此式（26）復數乘法次數有

2.1.2 TSTFT的算法實現

接下來討論TSTFT的計算，它可以看成窗函數時移t以后再取共軛，然后與原始信號對應相乘，最后做一次FT得到，故TSTFT可以寫成

其中：Fτ[·]表示對變量τ進行FT。

現將式（28）中的時移t離散化成長度為L2的序列：{n|n∈Z+0，n≤L2-1}，得到離散化的TSTFT

考慮式（29）計算復雜度，其復數乘法次數為

2.1.3 兩種短時傅里葉變換公式的轉換

對 比 式（27）和 式（30），發(fā) 現 當L1=L2時，的計算量大于的計算量，但在實際情況下，信號數據量N會很大，為了取得較好的時頻表示結果，L2也會取得很大；而由奈奎斯特頻率和最小頻率分辨率的限制，L1相對較小且不會隨著信號數據量N增大而增大。因此在實際情況中，采用作為STFT時計算效率更高。

在1.3節(jié)中已經論述，TSST采用TSTFT是為了避免出現相位振蕩，因此為了提高TSST計算效率，希望MSTFT和TSTFT之間能夠相互轉化。仔細觀察可以發(fā)現

兩種STFT公式只差一個相位ejωt，故它們的譜圖相同，這解釋了兩種STFT都能表征信號能量的原因。因此在計算TSTFT時，可以先按照MSTFT計算，然后在所有的計算點[b，k]上乘以相位修正值e-jbk即可得到TSTFT，這個過程可以表示為

2.2 IFO和GDO的 估 計

在兩種同步壓縮中，估計IFO和GDO是至關重要的一步，它們分別由式（2）和式（7）給出，通過觀察可以發(fā)現它們都涉及對STFT結果求偏導，如果直接由差分計算，會放大噪聲誤差。因此在知道窗函數解析式的條件下，可以利用STFT本身的性質來精確計算兩個估計值。

2.2.1 IFO的估計

對于MSTFT，有

結合式（2）有

對式（34）進行離散化，可得

2.2.2 GDO的估計

對于TSTFT，有

結合式（7）有

同樣，將式（37）進行離散化，可以得到

值得指出的是，為了保證算法的穩(wěn)定性，需要給出一個控制閾值γ0，只有時頻圖上或的點才能計算IFO或GDO，并計入后續(xù)的重排階段。

2.3 頻率重排和時間重排的實現

第3個步驟是按照估計算子重新排布STFT系數，對于SST重排過程，將式（4）離散化可得

其中：ω[p]=Δω×p，Δω為頻率軸的劃分間隔；p為頻率軸離散序號。

同樣的方式，對于TSST重排過程，將式（9）離散化可得

其中：τ[q]=Δτ×q，Δτ為時間軸的劃分間隔；q為時間軸離散序號。

2.4 SST信號重構和TSST信號重構的實現

接下來考慮兩種同步壓縮的重構實現。對于SST的重構離散化公式，可以由式（5）得到

同樣，對于TSST離散化重構公式，由式（10）可得

為了更清楚地表示SST和TSST的算法流程，整個計算過程以流程圖的方式給出，圖4為SST的計算流程，圖5為TSST的計算流程。

圖4 SST計算流程圖Fig.4 The flow diagram of the calculation for the SST

圖5 TSST計算流程圖Fig.5 The flow diagram of the calculation for the TSST

3 兩種同步壓縮在故障診斷中的應用

為了說明兩種同步壓縮各自的應用場合，文中采用定轉速下轉軸的碰摩故障數據和軸承外圈的沖擊故障數據分別對兩種同步壓縮算法進行驗證。

3.1 轉軸的碰摩故障識別

轉軸碰摩信號采集自一個重油催化機組［1，10］，其結構如圖6所示，它是由燃氣輪機、壓縮機、變速箱以及電動機構成，測試軸承（1#～4#）用于支撐相應的轉軸，轉軸轉速為5 381 r/min，振動信號的采樣率設置為2 000 Hz。

圖6 重油催化機組結構示意圖Fig.6 The structure of a heavy oil catalytic machine set

在測試過程中，來自2#測試軸承處的振動最大并且超過了報警閾值，因此重點分析2#測試軸承的數據。圖7展示了振動信號的波形及其頻譜圖，從圖7（a）中能看出比較規(guī)則的周期性時域信號，但很難發(fā)現有明顯的故障特征；同樣從圖7（b）中也僅能看出其一階轉頻（1X）及其高階倍頻成分，卻得不到任何的故障信息；圖7（c）為信號的STFT譜圖，由于機組為穩(wěn)定勻速工作，因此振動信號在時頻圖上表示為一組平行于時間軸的水平直線，其中最亮的代表1X分量。

為了更加清楚地描述1X分量的細節(jié)部分，將分析窗g(t)的時域寬度取窄，并只計算40～140 Hz的時頻表示，結果如圖8所示。

從圖8（a）和圖8（c）中可以看出，信號的STFT和TSST的能量在頻率方向完全彌散開來，而在圖8（b）中可以觀察到瞬時頻率微小波動的故障現象，因此將SST時頻表示結果中的瞬時頻率軌跡提取出來（見圖9（a）），并對其進行FT分析。

通過圖9（b）可以看出，瞬時頻率軌跡的頻率成分仍然為1X轉頻，這說明每當轉軸旋轉一圈，由于轉軸和機組固定件間的摩擦，轉軸就會產生一次局部的升速和降速，形成如圖8（b）所示的微小波動。因此在定速條件下，SST比TSST更容易識別轉軸碰摩故障。

圖7 碰摩故障信號Fig.7 The rub-impact fault signal

圖8 碰摩故障的時頻表示Fig.8 Time-frequency representations of the rub-impact fault

圖9 提取的瞬時頻率軌跡結果Fig.9 The extracted instantaneous frequency trajectory results

3.2 軸承外圈的沖擊故障檢測

圖10 測試臺結構Fig.10 The structure of the test rig

軸承外圈沖擊故障數據采用由凱斯西儲大學公布的故障軸承數據集［11-13］，其測試臺結構如圖10所示，它由電機、扭矩傳感器和測力計等部分組成。電機轉軸工作轉速為1 797 r/min，并由測試軸承所支撐，其中測試軸承外圈通過電火花加工引入了單點缺陷，加速度計安裝在電機一端并將采樣率設置為12 kHz，軸承的波形及其頻譜如圖11所示。

圖11 軸承沖擊信號Fig.11 Bearing shock signal

通過圖11（a）可以發(fā)現周期性的脈沖振蕩，但由于時域峰值的不規(guī)則性，很難通過時域波形的測量計算出軸承外圈的故障周期。圖11（b）展示了信號的頻帶，其能量主要集中在2 500～4 000 Hz，因此在時頻分析時，重點關注此頻帶的細節(jié)。

圖12 （a）～（c）分別展示了軸承沖擊信號在STFT，SST和TSST下的時頻表示結果，可以看出雖然STFT也能觀察到周期性沖擊規(guī)律，但因其時間分辨率不夠，難以測量故障周期；SST由于其壓縮方向是沿著頻率軸方向的原因，因此也不能提高其時間分辨率，最終難以獲得沖擊的故障特征；而TSST能清晰地表示故障特征。通過圖12（c）中的測量，可以知道軸承外圈沖擊故障周期為9.25 ms，因此軸承外圈的沖擊故障頻率為108.11 Hz，這與通過軸承尺寸（表1）所計算得到的理論外圈故障頻率（表2）107.36 Hz接近。故TSST能夠實現對軸承外圈沖擊故障頻率的檢測。

圖12 軸承沖擊故障的時頻表示Fig.12 Time-frequency representations of bearing shock fault

表1 軸承尺寸Tab.1 Drive end bearing size

表2 軸承缺陷頻率Tab.2 Bearing defect frequencies Hz

4 結束語

SST和TSST均可作為STFT的后處理工具，它們通過壓縮STFT系數達到細化脊線、提高時頻表示可讀性的目的，并且兩者都具備信號重構能力。SST和TSST最大的區(qū)別是它們有著不同的壓縮重排方向。對于SST來說，由于其采用MSTFT進行計算，避免了STFT系數在頻率軸方向的相位振蕩，從而可以實現頻率軸上的壓縮重排；類似地，對于TSST來說，由于采用了TSTFT進行計算，避免STFT系數在時間軸方向的相位振蕩，從而可以完成時間軸上的壓縮重排。雖然兩種同步壓縮采用了不同的STFT公式，但是由于MSTFT和TSTFT可以實現相互轉換，因此為算法的實現帶來了方便。在工程應用中，SST因其在頻率軸方向壓縮的特點適合處理類似于復正弦函數的緩變信號，如識別定轉速下轉軸的碰摩故障；相反，TSST因其在時間軸方向壓縮的特點更加適合處理類似脈沖函數的快變信號，如計算軸承外圈的故障頻率。因此將SST和TSST兩種方法相結合，就可以處理工程中常見的故障信號，具有廣泛的應用價值。