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

（1.上海交通大學(xué)機械系統(tǒng)與振動國家重點實驗室 上海，200240）

（2.博世華域轉(zhuǎn)向系統(tǒng)有限公司 上海，201821）

引言

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

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

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

1 兩種同步壓縮原理對比

1.1 SST簡述

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

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

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

對應(yīng)于MSTFT的逆變換公式為

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

因此，只需把移動后聚集在脊線附近的系數(shù)再進(jìn)行一次頻率積分就可以得到原始信號，故SST的重構(gòu)公式表達(dá)為

1.2 TSST簡述

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

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

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

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

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

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

1.3 SST與TSST的 區(qū) 別

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 兩種同步壓縮算法實現(xiàn)

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

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

2.1 MSTFT和TSTFT的 算 法 實 現(xiàn)

2.1.1 MSTFT的算法實現(xiàn)

對于SST，首先假設(shè)則MSTFT可以表示為

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

結(jié)合式（23）和式（24），得到MSTFT計算公式

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

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

2.1.2 TSTFT的算法實現(xiàn)

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

其中：Fτ[·]表示對變量τ進(jìn)行FT。

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

考慮式（29）計算復(fù)雜度，其復(fù)數(shù)乘法次數(shù)為

2.1.3 兩種短時傅里葉變換公式的轉(zhuǎn)換

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

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

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

2.2 IFO和GDO的 估 計

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

2.2.1 IFO的估計

對于MSTFT，有

結(jié)合式（2）有

對式（34）進(jìn)行離散化，可得

2.2.2 GDO的估計

對于TSTFT，有

結(jié)合式（7）有

同樣，將式（37）進(jìn)行離散化，可以得到

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

2.3 頻率重排和時間重排的實現(xiàn)

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

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

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

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

2.4 SST信號重構(gòu)和TSST信號重構(gòu)的實現(xiàn)

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

同樣，對于TSST離散化重構(gòu)公式，由式（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 兩種同步壓縮在故障診斷中的應(yīng)用

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

3.1 轉(zhuǎn)軸的碰摩故障識別

轉(zhuǎn)軸碰摩信號采集自一個重油催化機組［1，10］，其結(jié)構(gòu)如圖6所示，它是由燃?xì)廨啓C、壓縮機、變速箱以及電動機構(gòu)成，測試軸承（1#～4#）用于支撐相應(yīng)的轉(zhuǎn)軸，轉(zhuǎn)軸轉(zhuǎn)速為5 381 r/min，振動信號的采樣率設(shè)置為2 000 Hz。

圖6 重油催化機組結(jié)構(gòu)示意圖Fig.6 The structure of a heavy oil catalytic machine set

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

為了更加清楚地描述1X分量的細(xì)節(jié)部分，將分析窗g(t)的時域?qū)挾热≌?，并只計?0～140 Hz的時頻表示，結(jié)果如圖8所示。

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

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

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

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

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

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

圖10 測試臺結(jié)構(gòu)Fig.10 The structure of the test rig

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

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

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

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

圖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 結(jié)束語

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