广义瞬时速度同步化分步解调变换及其对旋转机械振动信号分析

石娟娟，花泽晖，沈长青,，江星星，冯毅雄，朱忠奎，孔 林

(1.苏州大学 轨道交通学院，江苏 苏州 215131；2.浙江大学 流体动力与机电系统国家重点实验室，杭州 310027；3.长光卫星技术有限公司，长春 130102)

旋转机械如汽轮机、风机、发动机被广泛用于电力和航空航天等领域。齿轮和轴承作为最常见的旋转机械零部件，在机械系统中发挥着连接和支撑机械结构的重要作用。随着实际需求的发展，旋转机械系统常被要求服役于高速重载等恶劣的复杂工况下，而这些关键零部件一旦发生意外故障，会直接影响到整个系统的平稳运行，严重时甚至还会造成经济损失和安全事故。同时，机械系统设计结构和工艺的复杂性以及可能存在的安装误差也会引起机械系统的不良振动并导致故障。因此，实时监测机械系统中关键零部件的健康状态可更加快速和精确地定位故障并做出相应的诊断和维护措施。通常旋转机械的故障诊断方法可通过振动、噪声、温度、声发射和油液等方法来进行，其中，振动检测涉及信号处理方面，而旋转机械振动信号常包含揭示系统健康状况的有用信息，是最常用的诊断方法。当旋转机械发生故障时，常表现为振动频率的变化，可通过检测振动频率、加速度、相位等参数进行分析。当旋转机械运行于变速工况下，信号的故障特征频率也是随时间而变化的，因此需要采用时频分析方法来揭示信号中频率随时间的变化。加上信号采集过程中不可避免地存在噪声干扰，还有加速度计有时无法安装来获取瞬时频率的信息等问题都给变转速工况下振动信号处理带来了新的问题和挑战。因此，需要提出更加先进和有效的时频分析方法来提高旋转机械振动信号中特征提取的准确性以及时获知系统关键零部件的健康状态，从而保障整个机械系统的安全运行和维护。

时频分析方法因能同时揭示信号的时间和频率信息，广泛用于非平稳信号的处理[1]。在时频分析方法中，最常见的是线性分析方法，如短时傅里叶变换(short-time Fourier transform, STFT)和小波变换(wavelet transform, WT)等[2]。之所以称为线性变换方法，是因为这类方法都利用了加窗信号和基函数的内积。这类方法受有限时频分辨率影响存在严重的时频模糊问题。类似的，还有一类双线性变换方法，如维格维纳分布、Cohen经典分布等。这类方法虽能提供较线性变换更高的时频分辨率，但存在严重的交叉项干扰问题。对线性变换得到的时频表示进行时频重映射这一类方法称为时频后处理，可进一步解决时频分析中因固定的时频分辨率导致的时频模糊问题。如Auger等[3]提出沿着时间和频率方向重分配时频分布的RM(reassignment method)方法来提升时频表示的可读性；Daubechies等[4]提出了同步压缩变换(synchrosqueezing transform, SST)结合小波变换实现信号分解。与RM在时间和频率同时压缩不同，SST仅沿频率方向重新分配时频表示，保留了信号的重构能力。虽然SST将时频分布集中在目标时频脊线周围，但是当信号的频率变化较快时，由时频分辨率本身存在的模糊问题仍未解决。并且当待分析信号包含较大的噪声时，噪声引起的时频模糊也被挤压到了时频脊线上，导致时频表示中目标频率分量处的能量聚集性降低，从而使时频脊线不清晰。类似的，Yu等[5]提出了同步提取变换(synchroextracting transform, SET)，通过将得到的时频表示与一个脉冲函数相乘，达到只保留目标频率处时频分布的效果从而提升时频表示可读性。这些时频重分配方法都是在不改变时频分辨率的情况下通过后处理从而进一步提升时频表示。除此之外，还有一类信号分解方法，如经验模态分解(empirical mode decomposition, EMD)、局部均值分解(local mean decomposition, LMD)和改进的变分模态分解(variational mode decomposition, VMD)等，这些方法旨在从原始的多分量信号中分离出目标频率分量(如旋转机械振动信号中的与故障特征相关的频率分量)，然后使用时频分析方法再对单个分量进行分析。这一类方法更适用于定转速工况下的轴承故障诊断，因为定转速工况下故障特征频率是一个定值，不随时间而变化，因此可根据频带划分而分离出来特征频率分量。因此这类方法还常与滤波和谱峭度结合，谱峭度指导选取含有最多相关特征所在的频带范围，然后经过带通滤波后再次分析可得到更加清晰的结果[6-8]。

为了解决时频分析中时间和频率上存在的时频模糊问题，广义解调(generalized demodulation, GD)因其能将变化的时频曲线解调到固定频率从而提升能量聚集性这一特性可有效提升时频表示的可读性。为了解决SST在分析频率快速变化信号时存在的频率模糊问题，提出了GD和SST结合的新方法[9-10]。该组合的优势在于广义解调可增强时频表示的能量聚集性，SST用于进一步锐化时频脊线，从而可进一步提升时频图的可读性。Wang等[11]提出迭代使用匹配解调变换方法(matching demodulation transform, MDT)来改善提取的瞬时频率的准确性；Shi等[12-13]提出了基于分形维数和广义解调的变速轴承故障诊断方法，避免了处理变速信号时的重采样过程。上述方法都已成功用于仿真分析和试验验证。但是，前述方法依赖瞬时频率相关先验知识，需要有迭代过程来改善瞬时频率估计的准确性和处理多分量信号从而提升时频表示[14]。

为了不依赖瞬时频率等先验知识并成功地运用广义解调以增强时频表示的能量聚集性，本文提出了广义瞬时频率同步化分步解调变换。该变换不需对瞬时频率进行预估计，并且避免了在处理多分量信号时所需的迭代操作。在足够小的时间窗内信号的瞬时频率脊线可近似为直线，该直线的斜率以及与横坐标轴的夹角也是唯一确定的。然后根据使用真实解调因子进行解调变换可将信号解调至水平直线并得到最高的能量聚集性这一特点，提出了最大峭度指导的自适应角度选取方法。为有效处理多分量信号，对原变换基函数进行扩展，得到与瞬时频率同步的线性变换基函数，由此可同步增强多个频率分量的时频表示，适合处理旋转机械所产生的多分量振动信号并且避免了算法中的迭代过程。此外，考虑到窗长的变化也会对时频图可读性产生影响，因此同时变化角度和窗长，利用峭度最大准则，提出在每个时刻选取合适参数组合的方法，使所提方法可在一定的噪声下自适应选择合适的参数并避免对先验知识的依赖。根据提升的时频表示，可更加准确地从非平稳振动信号中提取出相关时频特征并用于变转速旋转机械的健康状态监测。

1 广义瞬时速度同步化分步解调变换

1.1 广义分步解调变换

广义解调可处理非平稳信号并增强时频聚集性，在于该方法可将时变的频率曲线解调到固定频率上。通过这一步解调操作，还可有效地把多分量信号分解成多个代表不同分量的单分量信号。给定信号x(t)，广义解调变换可表示成

(1)

式中：d(t)为解调因子；e-j2πd(t)为变换基函数。可以发现,x(t)的广义解调变换等同于对x(t)e-j2πd(t)的标准傅里叶变换。当d(t)=0，式(1)就等价于标准的傅里叶变换。类似地，根据傅里叶变换的性质，可定义逆广义傅里叶变换为

(2)

观察式(2)，可以发现如果XG(f)≡δ(f-f0)，x(t)可以被写成

x(t)=ej2π[f0t+d(t)]=ejφ(t)

(3)

式中,φ(t)为信号的瞬时相位。由此可发现，信号原先的瞬时频率f(t)=f0+d′(t)=dφ(t)/2πdt被投影到一固定频率f=f0，其中d′(t)表示d(t)的一阶导。根据式(3)可知，如需将时变频率曲线映射到一固定频率，只需计算相关的解调因子并对信号进行解调。根据上述分析，可总结出广义解调的两个主要特性：①时变的瞬时频率曲线可被解调到某一恒定频率；②解调后信号能量集中在频率f=f0，较解调之前的信号具有更高的能量集中性。

据此，Shi等提出广义分步解调来研究振动信号时频特性，针对每一段截取信号将所提方法分为两步：①前向解调增强时频聚集性；②后向解调将时频能量聚集性得到增强的瞬时频率脊线返回至真实频率。具体介绍如下：短时窗内的信号频率可写成

f(t)=f0+d′(t)

(4)

前向解调因子可计算为

xd f(t)=x(t)·df(t)=ej2π[f0t+d(t)]e-j2πd(t)=ej2πf0t

(5)

式中：df(t)为前向解调因子；xd f为前向解调后的信号。式(5)的物理含义即为信号在乘上前向解调因子后，原先变化的频率f(t)被变换到了f0处。

为更清晰说明上述过程，广义分步解调变换的示意图，如图1所示。从图1中可以看出，在时刻t=τ原先的频率曲线f0+d′(t)被前向解调到了频率f0处，比较图1(b)中的幅值和频率的关系可发现解调后的时频脊线能量更加集中，即时频分布的能量相比解调之前更加聚集在频率f0周围，但在时间段[τ-Δt,τ+Δt]内，信号前向解调后的频率并非信号原先的真实频率，因此还需要后续操作(如后向解调)来恢复分析信号的真实频率，如图1(c)所示。在后向解调中，时频聚集性与前向解调之后一致，仅发生频移，由此后向解调可定义为

FD为前向解调；BD为后向解调。

xd(t)=xd f(t)db(t)=ej2πf0tej2πd′(τ)t=ej2π[f0t+d′(τ)t]

(6)

式中：db(t)为后向解调因子；d′(τ)>0为脊线上移；d′(τ)<0为下移。通过上述分步解调过程，可得到能量分布更加集中的时频表示。但是，上述的解调过程依赖瞬时频率的先验知识，即需首先对瞬时频率进行预估计，据此构建正确的解调因子来对信号进行解调变换。在不知道瞬时频率的前提下，目前已有文献大多是从得到的时频表示中先粗略地估计瞬时频率来对信号进行解调，然后利用迭代操作来提升瞬时频率估计的准确性。一方面，解调因子的准确性极大地依赖瞬时频率的预估计；另一方面，涉及的迭代操作使算法本身更复杂，并且在处理多分量谐波信号时，还需要再次涉及迭代操作来提升多分量的时频表示。

1.2 广义瞬时速度同步化分步解调变换

为避免对瞬时频率预估计的依赖并简化算法对多分量信号所需的迭代过程，提出了本文的研究方法。

与广义分步解调变换相同，在短时窗截取的信号，该段时间窗内的瞬时频率可近似看成线性变化，根据泰勒展开多项式，该瞬时频率可写成

f(t)≈f′(τ)(t-τ)+υ

(7)

式中：f(t)为信号的瞬时频率；f′(τ)为瞬时频率在t=τ处的斜率；υ为该时刻t=τ的频率。该瞬时频率还可写成

(8)

引入参数α，定义为

(9)

将式(9)代入式(8)，瞬时频率可表示成

f(t)=υ[tanα(τ)(t-τ)+1]

(10)

式(10)表明，引入随时间变化的参数α来对斜率进行描述，由此可用角度有限的变化范围(-π/2, π/2)来映射斜率的变化范围。根据式(4)，解调因子满足

(11)

将式(10)代入式(11)，即可计算得到前向解调因子的表达式，写成

(12)

式中：f0为解调后的频率；υ为时刻t=τ的频率。可看出，式(12)中前向解调因子是与tanα有关的函数。只要确定tanα(τ)，便可对信号进行前向解调，前向解调后的信号可表示成

(13)

式(13)表明，只要角度α(τ)选取正确，能正确匹配待分析信号的真实频率变化，该加窗信号即可被解调到固定频率f=f0处。

然而，前向解调后信号所处频率并不是待分析信号在该时刻的真实频率，因此还需对前向解调后的信号进行后向解调以恢复其真实频率。根据傅里叶变换性质，频域内发生频移，等同于在时域内乘上复指数函数，可表示为

ej2πυ(t-τ)

(14)

式中,db(t)为后向解调因子。通过观察式(14)可发现，两次解调后的真实频率得以保持，而解调后信号频谱在当前时刻的能量聚集性得到了提高，即原先线性的频率直线先被解调至固定频率以聚集能量，然后再通过第二次解调恢复信号的真实频率。整合前、后向解调因子，可得完整的解调因子为

d(t)=df(t)·db(t)=

(15)

根据式(15)定义的解调因子，可完成对信号的前、后向解调。原始信号的加窗时频表示为

XSTFT(τ,υ)=〈x(t)g(t-τ),ej2πυ(t-τ)〉=

(16)

式中：g(t-τ)为窗函数；υ为时刻t=τ的频率。考虑式(15)定义的解调因子，则解调后信号时频表示可写成

XG[τ,υ,α(τ)]=

(17)

式(17)即为引入了倾斜角α的广义分布解调变换。

所提变换还可满足对信号中目标频率分量的时域重构，具体证明为

2πg(0)x(τ)

(18)

式中，δ为脉冲函数。

根据式(15)～式(17)可发现，所提变换关键在于不依赖先验知识而找到正确的角度参数α以成功解调待分析信号并提升其时频聚集性。同时，为对解调之后信号的时频聚集性进行量化描述，引入峭度指标。通过离散化角度α，分别对所截取的信号进行分析，并计算出每个角度所对应的峭度值，其中，正确的角度因能匹配真实频率(与真实频率最接近)，可将该分段信号解调至一固定频率并具有最高的时频聚集性，即对应峭度的最大值。因此，可根据峭度最大值来选取准确的角度α对原始信号进行解调，从而获得更佳的时频聚集性。相反，若所选取角度α不能反映该截取信号真实的瞬时频率，则解调之后信号的频率不为定值(解调后的时频脊线仍无法平行于时间轴)。根据广义解调的性质，此时该时频表示中信号能量相对分散，对应的峭度值也更小。

为对上述分析进行说明，定义仿真信号x(t)，写成

(19)

式中,f(t)为瞬时频率，写成

f(t)=40sin(0.36t+0.6)+10sin(2t-2.5)+

5.4sin(3.6t-7.6)

(20)

所定义信号原先的时频表示及广义解调后信号的时频表示如图2(a)所示，时变的瞬时频率首先被变换到了固定频率f0。图2(b)中，解调因子由真实角度计算得到，而在图2(c)中，解调因子是由有别于真实值的角度计算所得，对比两个时频表示的谱图，如图2(d)所示，用真实解调因子解调所得信号的时频能量更集中(实线代表采用真实角度计算得到的解调因子对信号进行解调所得时频幅值；虚线为与真实值不同的角度得到的结果)，同时将两种解调结果与真实值对比(原先真实的时频能量如图2(d)中点划线所示)。频谱分析结果也说明用真实的解调因子可将频率曲线解调至水平线并获得较高的能量，即解调后能量的聚集性可用于指导角度的选取，从而实现在无需预先知晓瞬时频率(瞬时速度)的前提下对信号进行广义解调变换以增强时频聚集性。

图2 采用不同角度对信号进行解调分析的结果

对于含N个点的频谱S，峭度K可被计算成

(21)

为了正确选取式(9)定义的角度，对角度进行离散化，写为

(22)

式中，I决定了角度的个数。将式(22)代入式(12)，前向解调因子可写成

(23)

式(23)表明，每个角度α都可计算得到一个特定的解调因子。角度离散化的个数增多会提高解调因子估计的准确性，但同时也会导致计算量的增大。角度和最大峭度的对应关系可表示为

(24)

式中，α*(τ)为时刻τ处由频谱计算得到的最大峭度值所对应的角度。

为验证最大峭度指导的角度是否正确，对式(20)提出的仿真信号进行详细分析，得到的结果如图3所示。该单分量信号的短时傅里叶分析结果如图3(a)所示，矩形代表分步解调的时间区间，虚线代表t=1.155 s时刻的时频表示切片。根据式(21)和式(23)，图3(a)中虚线处的时频分布在经不同角度构建的解调因子作用后所得时频表示的峭度变化规律如图3(b)所示，在图3(b)中三角标记出了一系列角度对应的峭度的最大值。通过滑动时间窗，同步估计每个时刻点最大峭度值对应的角度，结果如图3(c)所示。另外，图3(c)还给出了由式(20)直接计算得到的真实角度来验证峭度指导角度选取的准确性，经对比可以发现，由峭度指导选取的角度变化规律与真实情况基本吻合，说明峭度可用于判断所选的角度是否合适。在数值上，分析了所估计的角度与真实值的平均相对误差，经计算误差为0.98%，说明了所提峭度指导的角度估计策略的有效性。最终，所提方法得到的时频表示如图3(d)所示。通过与STFT结果(见图3(a))对比可以发现所提方法可得到时频能量更加集中的时频脊线和更高的时频聚集性，尤其是在待分析信号的频率变化较快时，如2～4 s。

图3 单分量信号的角度估计

除了常见的峭度指标，平滑因子(smoothness index, SI)也常被用于衡量信号的冲击程度。为进行比较，将平滑因子也用于指导角度的选取。选用式(21)中的频谱S，平滑因子可定义为

(25)

式(25)表明，平滑因子SSI可看成是几何均值与算术均值的比值，范围在0～1。当频谱中存在明显的脉冲聚集时，SSI更接近于0。考虑到时频表示中信号更聚集的时频能量在频谱中表现为幅值更高的脉冲成分，如图2(b)～图2(d)所示。因此，所估计角度和平滑因子的对应关系为

(26)

式中，α*(τ)为时刻τ处由频谱计算得到的最小平滑因子所对应的角度。

再次对图3中所用信号进行分析，峭度和平滑因子指导的角度结果，如图4所示。计算平滑因子和峭度所估计角度的平均相对误差，分别为4.86%和0.99%。经对比可发现平滑因子虽能揭示角度大致变化趋势，但在信号首尾两端估计精度没有峭度高且同时存在振荡现象。因此可得出结论：峭度指标可更准确地指导角度的选取。

图4 平滑因子和峭度所估计角度

该信号的重构结果如图5所示，经计算，重构信号与原始信号的平均相对误差为1.07%，说明了重构的准确性。

图5 单分量信号时域重构结果

所提方法的简要流程如图6所示。在每个时刻依据最大峭度准则估计得到合适的角度，通过构建合适的解调因子匹配待分析信号的频率变化趋势，从而提升时频表示。

图6 最大峭度指导下的角度估计流程图

1.3 改进的线性变换基函数

上述分析是以假设信号中仅有单个频率分量为前提进行的，但是在实际应用中，信号常常包含多个分量，如：旋转机械的振动信号中就包含多个分量，有轴转速及其谐波和故障特征频率及其谐波等多个分量。因此在分析多分量信号时，传统的分析方法(如可将STFT看成用角度α=0来解调信号)还可能会产生额外的时频模糊问题。为了能同步增强多个频率分量的时频表示，研究扩展后的线性变换基函数。先定义一个含有多个频率分量的信号，写成

(27)

式中,fk(t)=Ckf(t)，Ck=2.0,Ck=3.2,Ck=3.5,Ck=3.8,Ck=5.0，f(t)为瞬时频率并在式(20)定义。式(17)中定义的线性变换基函数，还满足如下关系

(28)

对该多分量信号进行时频分析，STFT结果如图7(a)所示，图7(a)中矩形区域经放大后的结果如图7(b)所示。结合这两个图可以看出，STFT的时频表示受较严重的交叉项干扰，图中的三个相对较近的频率分量难以被区分开，并且得到的时频聚集性欠缺。采用所提方法对该信号进行分析，得到的分析结果如图7(c)所示。图7(c)中矩形区域对应的时频局部放大结果如图7(d)所示。与图7(a)和图7(b)中STFT的结果相比，可以发现，所提方法可得到能量更加集中的时频表示和更容易分辨的时频脊线，并且在一定程度上解决了STFT中因基函数频率与待分析信号时频特征不匹配而引起的交叉项干扰。

广义瞬时频率同步化分步解调变换在同步增强多分量时频表示时，峭度指导下的角度选取方法与单分量一致，估计的角度如图7(e)和图7(f)所示。在时刻t=1.155 s处峭度随角度的变化如图7(e)所示，图7(e)中，三角表示最大峭度，表明该最大峭度值对应的角度为该时刻最终选取的角度。同步估计所有时刻的角度，结果如图7(f)所示，其中，实线代表估计所得角度，虚线代表真实角度，可见所估计的角度与真实值基本吻合。经计算，角度估计的平均相对误差为0.94%，说明在最大峭度指导下所提方法能同时提升多个频率分量的时频聚集性。对所得时频表示进行瞬时频率估计，然后根据估计的目标频率分量进行时域信号重构，结果如图8所示。经计算，重构的信号与原始信号的平均相对误差为1.28%，说明所提方法拥有信号的重构性能[16]。

图7 多分量信号角度估计

图8 多分量信号的重构

为比较所提方法对噪声的抗干扰能力，对图7中所分析信号添加适量噪声以进一步分析。含噪信号的信噪比从-5 dB变化至20 dB，估计角度的平均相对误差如图9所示。从图9中可以看出，随着信噪比降低，误差也逐渐增大，说明了噪声对角度估计的准确性有一定的影响。从图9可以看出，在分析信噪比低为-5 dB的信号时，所提方法仍能比较准确地估计角度，计算的误差均小于5%。

图9 所估计角度的平均误差随输入信噪比变化

1.4 窗长的确定

上述分析中，单分量和多分量信号分析结果均是采用固定长度的窗函数分析所得到的，考虑到窗长的变化对时频表示的分辨率也会有较大影响，故本节进一步分析窗长变化对最终时频表示的影响。同样地，根据峭度最大准则来选取每个时刻合适的分析窗长。同时考虑窗长和角度变化的时频表示记为XG[τ,ω,g(τ),α(τ)]。根据式(24)，在同时考虑窗长与角度变化情况下最大峭度对应的参数选择可表示为

[g*(τ),α*(τ)]=

(29)

式中，g*(τ)和α*(τ)分别为τ时刻峭度指导选取的窗长和角度。再次分析1.3节中的多分量信号，结果如图10所示。图10(a)给出了所提方法在t=1.155 s处角度、窗长与峭度的对应关系，并用三角标出了计算得到的峭度最大值，通过该点的坐标可知与最大峭度所对应的角度与分析窗长。依据式(29)，所分析信号在整个时段选择的角度和窗长如图10(b)和图10(c)所示。从图中可以看出所选取角度基本与真实值吻合，经计算平均相对误差为1.24%。在图10(c)中，分析窗长在前2 s内经历了较多的波动，但仍可以看出，在频率变化相对缓慢时，可选用相对较长的窗来分析信号，比如0～2 s左右的窗相对3 s更长(3 s左右信号频率变化相对0～2 s更快)。最终结果如图10(d)所示，与图7(c)获得的时频表示(固定窗长下的多分量分析结果)相比，可以发现同时考虑角度和窗长可得到更加集中的时频脊线，也说明了选取合适的窗长可进一步增强时频表示。

图10 所提方法中角度和窗长的选择

2 轴承仿真信号分析

实际采集的信号不可避免会包含噪声，并且噪声也会对峭度指导的参数选择产生影响。因此本章中在添加噪声后进一步验证所提方法的有效性并将结果与常见时频分析方法所得结果进行对比。

考虑旋转机械振动信号常包含多个分量，定义的含局部故障的轴承仿真多分量信号为

x(t)=

(30)

式中：fISRF为仿真信号的轴转频；n(t)为噪声。

该仿真信号的采样频率设为8 000 Hz，共振频率设为3 500 Hz，信号持续1 s。设定的故障特征系数为2.7,即fIFCF=2.7fISRF。

添加适量噪声，将信噪比设定为3 dB进行分析。含噪信号的时域波形如图11(a)所示，用所提方法对信号进行分析，在t=0.5 s角度、窗长与峭度值的对应关系，如图11(b)所示，图11(b)中三角代表该时刻计算得到的最大峭度值及其对应的角度和窗长。同样，对整个信号分析角度和窗长变化的影响，得到的结果分别如图11(c)和图11(d)所示。图11(c)中，所估计角度的平均相对误差为1.62%。在图11(d)中可看出，前0.5 s选取的窗长比在1.5～2.0 s左右更长，这与时频分析方法对窗长的要求基本吻合。

图11 含噪多分量信号分析

再次采用不同时频分析方法对信号进行分析，包括STFT，广义线性调频变换(general linear chirplet transform, GLCT)[17]，多项式调频变换(polynomial chirplet transform, PCT)[18], SST,SET和所提方法，得到的时频表示如图12所示。在分析该仿真信号时，所有时频方法的窗长统一设为0.4 s以方便对比。图12(a)中所示STFT结果仅在第1阶转频fISRF以及t=1.7 s左右有较好的时频聚集性。GLCT和PCT方法的分析结果如图12(b)和图12(c)所示, GLCT算法中所用调频率在结果中导致了较为严重的时频交叉干扰，从而难以获得准确的瞬时频率估计且难以观测得到2倍转频对应的频率分量存在；而PCT因采用固定的线性变换基函数，适合处理具有相同变化趋势的信号。图12(c)中，基函数设置与fISRF一致，因此fIFCF及其谐波分量的能量相对分散。此外,PCT还需依赖转频相关的先验知识。

SST和SET时频后处理方法的分析结果如图12(d)和图12(e)所示。SST因仅沿频率方向压缩，受快变频率和噪声的影响，得到的结果在0.5～0.8 s遭受了较为严重的时频模糊问题，难以观测得到清晰的瞬时频率脊线，如2fIFCF；而SET通过类似脉冲函数的同步提取因子，其结果为提取因子与STFT时频表示的乘积，即在STFT所得时频表示的基础上提取目标频率分量(瞬时频率估计因子)及其周围的时频分布，因此不可避免地受到STFT原始时频图中时频模糊问题的影响，只能提供相比图12(a)更加锐利的时频脊线，而无法提升能量的聚集性。所提方法的结果如图12(f)所示，从图12(f)中可清晰观测得到fISRF和fIFCF的倍频及其谐波存在。虽然得到的结果中在某些时刻的时频脊线不如后处理方法(见图12(d)和图12(e))集中，但时频表示的能量聚集性相对更高，且不需要瞬时频率相关的先验知识。因此，所提方法与其他方法相比具有一定的优越性，可有效用于提升时频表示可读性并具有一定的抗噪性能。

图12 不同时频分析方法所得时频表示(含噪多分量信号)

3 旋转机械故障诊断的应用

变转速工况下旋转机械振动信号中同样含有多个频率分量，如旋转频率、故障特征频率及噪声干扰等。相比传统的广义解调方法在分析时需要考虑迭代，所提方法因与转速同步，可同时增强多个频率分量的能量聚集性从而简化算法。在获得时频表示后，根据时频图中的瞬时频率脊线可诊断出该旋转机械系统中是否存在故障及故障类型。

3.1 行星齿轮信号分析

将所提方法用于分析试验信号，该行星齿轮信号的故障试验台如图13所示。具体配置如表1所示，根据齿轮箱的配置可计算出各自的故障特征频率系数。

表1 齿轮箱的参数配置

1.电机;2.测速计;3.定轴齿轮变速箱;4.行星齿轮箱一级传动轴;5.行星齿轮箱二级传动轴;6.加速度计;7.制动器。

该旋转机械系统的转速fd是由交流电机驱动控制的，转速先增加然后降低。采集信号的时域波形如图14(a)所示，该信号的采样频率设为20 kHz，信号采集时间持续30 s。转频fd的变化如图14(b)所示，先从40 Hz变化到60 Hz，再变化到40 Hz。测速计采集的转速信息主要用来验证所提角度选择方法是否准确。对该信号进行频谱分析，0～400 Hz内的频谱如图14(c)所示，由于该信号是在变转速工况下采集得到，因此在频谱中观测不到定转速工况下特别突出的特征频率等峰值。通常试验信号包含故障所引起的共振频段，为避免失真，采样频率设置较大；而当试验信号采集时间也较长时，分析中若考虑窗长变化会极大地增加计算量。因此，采用固定窗长来比较所获得时频表示的能量聚集性。

取时刻t=20 s分析峭度随所选取的角度的变化规律，得到的结果如图14(d)所示，图14(d)中用三角标记了最大峭度及其对应角度。同步估计所有时刻的角度，结果如图14(e)所示，其中，虚线代表由测速计采集信号直接计算得到的真实角度，实线为峭度指导选取的角度。通过对比可发现，所选取角度的变化趋势基本与真实情况吻合。STFT,GLCT和所提方法的时频分析结果分别如图14(f)～图14(h)所示。对比时频分析结果表明：采用相同窗长分析时，所提方法由于找到了合适的角度来同步解调信号，可得到更高的能量聚集性和更集中的时频脊线，如时频表示中转速二倍频2fd处。GLCT也是通过匹配信号调频率来增强时频表示，但得到的时频分析结果同时受窗长和不准确的调频率影响，从而导致时频模糊问题并降低时频表示可读性。从图14(h)中，可清晰地观测到转频fd及其前5阶谐波分量和第一级齿轮啮合的频率分量fmesh1(fmesh1=100/27fd)，而并未发现与该行星齿轮故障相关的边频带存在，这些频率分量可说明该行星齿轮是健康的。

图14 行星齿轮箱振动信号分析

3.2 轴承内圈故障诊断

同样将所提方法用于实际轴承振动信号分析[19]。该轴承故障模拟试验台的示意图，如图15所示。该系统中包括两个轴承，右边的轴承存在内圈故障。信号采样频率设为20 kHz，持续10 s。

图15 轴承内圈故障试验台

采集得到的振动信号如图16(a)所示，该信号的瞬时旋转频率(电机转频)fISRF变化如图16(b)所示，这里采集的fISRF主要用于计算真实角度来验证所提方法根据最大峭度选取角度的准确性。该轴承的内圈故障特征系数为5.43，即故障特征频率fIFCF=5.43fISRF。选取时刻t=4 s分析峭度随角度的变化规律，结果如图16(c)所示，图16(c)中用三角标记出了最大峭度值和对应的角度。对整个信号分析的结果如图16(d)所示。可以发现，根据最大峭度选择的角度基本是与真实值一致的，经计算所估计角度的平均相对误差为1.19%。采用STFT以及SST作对比，分析结果如图17(a)和图17(b)所示。所提方法结果如图17(c)所示。经对比可以发现所提方法能够得到能量更加集中的时频脊线并表征出信号中微弱的频率特征成分如fIFCF+fISRF,2fIFCF+fISRF和3fIFCF-2fISRF。因所提方法同时增强多个频率分量的时频聚集性，可得到更精确的时频脊线，如图17(d)所示。根据搜寻得到的时频脊线以及存在的转频调制现象均可表明此试验轴承存在内圈故障[20]。因此所提方法可成功用于诊断旋转机械系统潜在的健康状态。

图16 轴承振动信号的角度估计

图17 不同方法对轴承振动信号的分析结果

4 结 论

本文针对变转速工况下旋转机械非平稳信号处理存在时频模糊和脊线提取困难等问题，基于广义解调提出了广义瞬时速度同步化分步解调变换，用以对非平稳信号展开分析，并成功应用于变转速下旋转机械的健康状态监测。首先，引入角度变量并计算了广义分步解调变换中的前向和后向解调因子；然后，为减少对瞬时频率预估计的依赖、实现瞬时速度同步化，提出了峭度指导的角度(瞬时频率)自适应选择策略，可避免人为干预和对先验知识的依赖；最后，通过新的解调因子扩展原先的线性变换基函数，使所提方法能同时对多个分量进行解调并同步增强多个频率分量在时频表示中的能量聚集性。

针对信号重构，展开了证明，并对仿真信号进行了重构分析，结果表明：所提时频分析方法可根据提取的目标时频脊线重构出原始时域信号。此外，对窗长变化对时频分析的影响也展开了探讨，说明动态窗长可进一步提升信号的时频表示。仿真分析和试验验证说明该方法能得到更加集中的时频脊线和更高的时频能量聚集性从而提升时频表示可读性。对增强后的时频表示进行分析，可更准确地提取出变转速工况下旋转机械振动信号中与故障相关的特征，从而实现对变转速旋转机械的故障诊断。