何周杰，涂晓彤，王 凯，李富才，包文杰，包 隽

（1.上海交通大学机械系统与振动国家重点实验室 上海，200240）

（2.博世华域转向系统有限公司 上海，201821）

引言

机械系统的故障诊断中，传感器直接获取的一般是时域信号，时域信号能准确地反映机械系统在不同时刻下的运行状态。对时域信号进行傅里叶变换（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］。在此之后，为了进一步提高信号能量的集中程度，同步压缩出现了不同的发展方向：通过将匹配解调与同步压缩相结合，发展出了广义同步压缩变换［5］；将二阶或高阶多项式引入瞬时频率模型之中，得到二阶或高阶同步压缩变换［6-8］；将同步压缩过程和脊线提取过程相结合，得到同步提取变换［1］。上述方法均是建立在SST的基础之上提出的，由于它们的压缩重排是沿着频率轴方向进行的，因此很难处理信号脊线是平行于频率轴的时频压缩，例如在轴承冲击信号的时频表示中，常常出现平行于频率轴的周期性脊线特征。为此人们提出了时间同步压缩，它的压缩重排是沿着时间轴方向进行的，因此其具备处理强时变信号的能力［9］。

笔者以SST和TSST为基础，通过对比两者所采用的STFT来说明两种同步压缩的区别以及各自的应用场合，然后按照同步压缩的流程推导出它们的实现算法，最后通过转轴的碰摩故障和轴承外圈冲击故障的实验数据来验证两种算法的有效性。

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的 区 别

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的函数，故的相位沿着时间轴和频率轴方向都将发生相位振荡，因此对沿着时间轴方向压缩的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分别表示图中两个椭圆的内部区域。其中区域1代表信号频率缓慢变化的部分，区域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值在区域1中基本相同而在区域2中差异较大，故SST仅能在区域1内取得很好的压缩效果（如图3（d）所示）。相反，由图3（c）可以看出，在水平方向上所有的点的GDO值在区域1中差异较大而在区域2在中基本相同，故TSST仅在区域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），发 现 当L1=L2时，的计算量大于的计算量，但在实际情况下，信号数据量N会很大，为了取得较好的时频表示结果，L2也会取得很大；而由奈奎斯特频率和最小频率分辨率的限制，L1相对较小且不会随着信号数据量N增大而增大。因此在实际情况中，采用作为STFT时计算效率更高。

在1.3节中已经论述，TSST采用TSTFT是为了避免出现相位振荡，因此为了提高TSST计算效率，希望MSTFT和TSTFT之间能够相互转化。仔细观察可以发现

两种STFT公式只差一个相位ejωt，故它们的谱图相同，这解释了两种STFT都能表征信号能量的原因。因此在计算TSTFT时，可以先按照MSTFT计算，然后在所有的计算点[b，k]上乘以相位修正值e-jbk即可得到TSTFT，这个过程可以表示为

2.2 IFO和GDO的 估 计

在两种同步压缩中，估计IFO和GDO是至关重要的一步，它们分别由式（2）和式（7）给出，通过观察可以发现它们都涉及对STFT结果求偏导，如果直接由差分计算，会放大噪声误差。因此在知道窗函数解析式的条件下，可以利用STFT本身的性质来精确计算两个估计值。

2.2.1 IFO的估计

对于MSTFT，有

结合式（2）有

对式（34）进行离散化，可得

2.2.2 GDO的估计

对于TSTFT，有

结合式（7）有

同样，将式（37）进行离散化，可以得到

值得指出的是，为了保证算法的稳定性，需要给出一个控制阈值γ0，只有时频图上或的点才能计算IFO或GDO，并计入后续的重排阶段。

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）中能看出比较规则的周期性时域信号，但很难发现有明显的故障特征；同样从图7（b）中也仅能看出其一阶转频（1X）及其高阶倍频成分，却得不到任何的故障信息；图7（c）为信号的STFT谱图，由于机组为稳定匀速工作，因此振动信号在时频图上表示为一组平行于时间轴的水平直线，其中最亮的代表1X分量。

为了更加清楚地描述1X分量的细节部分，将分析窗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）可以发现周期性的脉冲振荡，但由于时域峰值的不规则性，很难通过时域波形的测量计算出轴承外圈的故障周期。图11（b）展示了信号的频带，其能量主要集中在2 500～4 000 Hz，因此在时频分析时，重点关注此频带的细节。

图12 （a）～（c）分别展示了轴承冲击信号在STFT，SST和TSST下的时频表示结果，可以看出虽然STFT也能观察到周期性冲击规律，但因其时间分辨率不够，难以测量故障周期；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最大的区别是它们有着不同的压缩重排方向。对于SST来说，由于其采用MSTFT进行计算，避免了STFT系数在频率轴方向的相位振荡，从而可以实现频率轴上的压缩重排；类似地，对于TSST来说，由于采用了TSTFT进行计算，避免STFT系数在时间轴方向的相位振荡，从而可以完成时间轴上的压缩重排。虽然两种同步压缩采用了不同的STFT公式，但是由于MSTFT和TSTFT可以实现相互转换，因此为算法的实现带来了方便。在工程应用中，SST因其在频率轴方向压缩的特点适合处理类似于复正弦函数的缓变信号，如识别定转速下转轴的碰摩故障；相反，TSST因其在时间轴方向压缩的特点更加适合处理类似脉冲函数的快变信号，如计算轴承外圈的故障频率。因此将SST和TSST两种方法相结合，就可以处理工程中常见的故障信号，具有广泛的应用价值。