0 引　言

船用螺杆式滑油泵是船舶常见的辅机设备，对于主轴推力轴承、主机轴承的润滑必不可少，准确高效开展该类设备故障诊断有着积极的意义。其振动故障常见的原因有泵管路进气、泵进口管路真空度过大引起气蚀、滑油压力波动大、轴承故障等，其中气蚀、轴承故障等都产生冲击信号，开展故障诊断时不容易分辨，导致不能精准排除故障。

冲击信号常用的分析方法有基于Hilbert变换的包络谱分析方法^[1]，以及在此基础上发展起来的结合去噪前处理的包络谱分析方法^[2-5]。但部分去噪方法运用了ICEEMD和VMD等模态分解时频分析方法，而模态分解时频分析方法实际上隐含信号具有慢时变特性的假设，因而对于信号中的冲击成分较为敏感^[6]，可能得不到较准确的分析结果。

文献[7]提出了一个应用频率变化模型描述类似冲击信号的TSST（time-reassigned SST）方法。文献[8]在该方法基础上提出了一种基于不动点迭代算法的高能量集中TSST方法，能够更好地用于分析具有快速变化特性的复杂信号，该方法称为时间重新分配多次同步压缩变换（time-reassigned multisynchrosqueezing transform, TMSST）。TMSST方法针对TSST方法因应用狄拉克δ函数仅能处理弱频变信号方面的局限性，采用迭代法解决了模糊时频分析问题，然后提出了一种提取冲击特征用于信号重构的算法并编程实现^[9]。总之，TMSST方法能够较好处理冲击信号，有望在船用螺杆式滑油泵故障诊断中得到应用。

1 方法介绍 1.1 时间重分配同步压缩变换（TSST）

频率变化的单分量信号可以描述为：

$ \hat s(\omega ) = A(\omega ){e^{i\varphi (\omega )}}，$

(1)

其中 $A(\omega )$ 和 $\varphi (\omega )$ 表示频域中的信号幅度和相位， $ - {\varphi ^\prime }(\omega )$ 表示群时延(group delay，GD)。当该信号出现在传递函数域时，理想时频分析可以表示为：

$ {ITFR} (t,\omega ) = A(\omega ){e^{i\varphi (\omega )}}\delta \left( {t + {\varphi ^\prime }(\omega )} \right)，$

(2)

其中 $\delta ()$ 是狄拉克函数。根据式(2)，在理想情况下，信号的时频特征应该只出现在群时延（GD）轨迹上。短时傅里叶分析（short time Fourier analysis, STFT）可用于将式(1)给出的信号扩展到时频域。在频域中，使用移动窗函数 $\hat g(\xi )$ 信号的STFT结果可以写为

$ G(t,\omega ) = {(2\text{π} )^{ - 1}}\int_{ - \infty }^{ + \infty } {\hat s} (\xi )\hat g(\xi - \omega ){e^{i(\xi - \omega )t}}{\rm{d}}\xi 。$

(3)

为了进一步研究由式(1)给出信号的STFT结果，假设分析的信号是弱频变信号，即 $\exists \varepsilon $ 足够小， $\left| {{A^\prime }(\omega )} \right| \leqslant $ $ \varepsilon $ 和 $\left| {{\varphi ^{\prime \prime }}(\omega )} \right| \leqslant \varepsilon $ ，则得到信号的一阶展开式：

$ \hat s(\xi ) = A(\omega ){e^{i\left( {\varphi (\omega ) + {\varphi ^\prime }(\omega )(\xi - \omega )} \right)}}，$

(4)

把式(4)代入式(3)得

$ \begin{split} G(t,\omega ) =& {(2\text{π} )^{ - 1}}\int_{ - \infty }^{ + \infty} A (\omega ){e^{i\left( {\varphi (\omega ) + {\varphi ' }(\omega )(\xi - \omega )} \right)}}\hat g(\xi - \omega ){e^{i(\xi - \omega )t}}{\rm{d}}\xi =\\ & {(2\text{π} )^{ - 1}}A(\omega ){e^{i\varphi (\omega )}}\int_{ - \infty }^{ + \infty } {{e^{i\left( {t + {\varphi ' }(\omega )} \right)(\xi - \omega )}}} \hat g(\xi - \omega ){\rm{d}}\xi =\\ & A(\omega ){e^{i\varphi (\omega )}}g\left( {t + {\varphi ^\prime }(\omega )} \right)。\\[-12pt] \end{split} $

(5)

式中， $g(t)$ 为时域窗函数，其中 ${supp} (g) \in [ - \Delta ,\Delta ]$ 。由式(5)可知时频能量围绕GD轨迹进行扩散。为了提高式(5)的能量浓度，二维GD轨迹估算为：

$ \hat t(t,\omega ) = {Re} \left( {\frac{{i{\partial _\omega }G(t,\omega )}}{{G(t,\omega )}}} \right) ，$

(6)

把式(5)代入式(6)得

$ \hat t(t,\omega ) = - {\varphi ^\prime }(\omega )。$

(7)

沿时间方向进行一维积分，将模糊的TF能量压缩到GD轨迹中：

$ Ts(u,\omega ) = \int_{ - \infty }^{ + \infty } G (t,\omega )\delta (u - \hat t(t,\omega )){\rm{d}}t，$

(8)

式(7)和式(8)相结合得

$ \begin{split} Ts(u,\omega ) =& {(2\text{π} )^{ - 1}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {\hat s} } (\xi )\hat g(\xi - \omega ){e^{i(\xi - \omega )t}} {\rm{d}}\xi {\rm{d}}t\\ &\delta \left( {u + {\varphi ^\prime }(\omega )} \right)=\\ &\int_{ - \infty }^{ + \infty } {\hat s} (\xi )\hat g(\xi - \omega )\delta (\omega - \xi ){\rm{d}}\xi \delta \left( {u + {\varphi ^\prime }(\omega )} \right) = \\ &\hat s(\omega )\hat g(0)\delta \left( {u + {\varphi ^\prime }(\omega )} \right) 。\end{split} $

(9)

说明对于式(4)给出的微弱的频率变化信号，TSST可以通过将模糊的TF能量压缩到GD轨迹中来产生理想的信号表示。

1.2 时间重分配多同步压缩变换（TMSST）

1）TMSST方法的建立

一个强频变信号，即 $\exists \varepsilon $ 足够小， $\left| {{A^\prime }(\omega )} \right| \leqslant \varepsilon $ 和 $\left| {{\varphi ^{\prime \prime \prime }}(\omega )} \right| \leqslant \varepsilon $ 。对于这样的信号，由式(1)给出的信号可以扩展为：

$ \hat s(\xi ) = A(\omega ){e^{i\left( {\varphi (\omega ) + {\varphi ^\prime }(\omega )(\xi - \omega ) + 0.5{\varphi ^{\prime \prime }}(\omega ){{(\xi - \omega )}^2}} \right)}}，$

(10)

STFT中使用的高斯窗函数的傅里叶变换可以表示为：

$ \hat g(\omega ) = \sqrt {2\sigma \text{π} } {e^{ - 0.5\sigma {\omega ^2}}} ，$

(11)

把式(10)和式(11)代入式(3)得

$ \begin{split} G(t,\omega ) & = {(2\text{π} )^{ - 1}}\int_{ - \infty }^{ + \infty } A (\omega ){e^{i\left( {\varphi (\omega ) + {\varphi ^\prime }(\omega )(\xi - \omega ) + 0.5{\varphi ^\varphi }(\omega ){{(\xi - \omega )}^2}} \right)}}\times\\ &\sqrt {2\sigma \text{π} } {e^{ - \frac{{\sigma {{(\xi - \omega )}^2}}}{2}}}{e^{i(\xi - \omega )t}}{\rm{d}}\xi =\sqrt {2\sigma \text{π} } A(\omega )\times{e^{i\varphi (\omega )}}{(2\text{π} )^{ - 1}}\\ &\int_{ - \infty }^{ + \infty } {{e^{0.5\left( {i{\varphi ^0}(\omega ) - \sigma } \right){{(\xi - \omega )}^2}}}} {e^{i\left( {t + {\varphi ^\prime }(\omega )} \right)(\xi - \omega )}}{\rm{d}}\xi = \\ &A(\omega ){e^{i\varphi (\omega )}}\sqrt {\frac{\sigma }{{\sigma - i{\varphi ^{\prime \prime }}(\omega )}}} {e^{ - \frac{{{{\left( {t + {\varphi ^\prime }(\omega )} \right)}^2}}}{{2\sigma - 2i{\varphi ^n}(\omega )}}}} 。\\[-15pt] \end{split} $

(12)

根据式(6)可以得到式(12)的二维GD轨迹估计值为：

$ \hat t(t,\omega ) = - {\varphi ^\prime }(\omega ) + \frac{{{\varphi ^{\prime \prime }}{{(\omega )}^2}}}{{{\sigma ^2} + {\varphi ^{\prime \prime }}{{(\omega )}^2}}}\left( {t + {\varphi ^\prime }(\omega )} \right)。$

(13)

根据式(13)，对于强频率变化的信号，式(6)的表达式不能准确估计信号的真实GD轨迹。为解决该问题，把 $t = - {\varphi ^\prime }(\omega )$ 代入(13)得

$ \hat t\left( { - {\varphi ^\prime }(\omega ),\omega } \right) = - {\varphi ^\prime }(\omega )。$

(14)

式(14)表示群时延 $ - {\varphi ^\prime }(\omega )$ 为 $\hat t(t,\omega )$ 的固定点。意味着可以采用定点迭代算法来减少 $ - {\varphi ^\prime }(\omega )$ 和 $\hat t(t,\omega )$ 之间的误差。第一个迭代可以写为：

$ \hat t(\hat t(t,\omega ),\omega ) = - {\varphi ^\prime }(\omega ) + {\left( {\frac{{{\varphi ^{\prime \prime }}{{(\omega )}^2}}}{{{\sigma ^2} + {\varphi ^{\prime \prime }}{{(\omega )}^2}}}} \right)^2}\left( {t + {\varphi ^\prime }(\omega )} \right) ，$

(15)

由式(15)可知，固定点迭代算法实际上可以构造一个新的二维GD轨迹估计 $\hat t(\hat t(t,\omega ),\omega )$ ，推导出

$ \hat t(\hat t(t,\omega ),\omega ) + {\varphi ^\prime }(\omega ) < \hat t(t,\omega ) + {\varphi' }(\omega ) 。$

(16)

式(16)意味着通过一次迭代，新的二维GD轨迹估计值 $t(t(t,\infty )$ 已经比 $\hat t(t,\omega )$ 更接近 $ - {\varphi ^\prime }(\omega )$ 。执行第二次迭代，则可以得出

$ \hat t(\hat t(\hat t(t,\omega ),\omega ),\omega ) = - {\varphi ^\prime }(\omega ) + {\left( {\frac{{{\varphi ^{\prime \prime }}{{(\omega )}^2}}}{{{\sigma ^2} + {\varphi ^{\prime \prime }}{{(\omega )}^2}}}} \right)^3}\left( {t + {\varphi ^\prime }(\omega )} \right), $

(17)

由式(17)可以得出

$ \hat t(\hat t(\hat t(t,\omega ),\omega ),\omega ) + {\varphi' }(\omega ) < \hat t(\hat t(t,\omega ),\omega ) + {\varphi' }(\omega )，$

(18)

将式(18)与式(16)进行比较，结果表明，对于每一次迭代，新构造的二维GD轨迹估计都更接近真实的 $ - {\varphi ^\prime }(\omega )$ 。令 ${\hat t^{[N]}}(t,\omega )$ 表示N次迭代后新构造的二维GD轨迹估计，则有：

$ {\hat t^{[N]}}(t,\omega ) = - {\varphi ^\prime }(\omega ) + {\left( {\frac{{{\varphi ^{\prime \prime }}{{(\omega )}^2}}}{{{\sigma ^2} + {\varphi ^{\prime \prime }}{{(\omega )}^2}}}} \right)^N}\left( {t + {\varphi ^\prime }(\omega )} \right)。$

(19)

当迭代次数足够大时， ${\hat t^{[N]}}(t,\omega )$ 将无限接近 $ - {\varphi ^\prime }(\omega )$ ，即

$ \mathop {\lim }\limits_{N \to \infty } {\hat t^{[N]}}(t,\omega ) = - {\varphi ^\prime }(\omega )，$

(20)

将 $\hat t(t,\omega )$ 和 ${\hat t^{[N]}}(t,\omega )$ 替换到式(8)中有

$ T_S^{[N]}(u,\omega ) = \int_{ - \infty }^{ + \infty } G (t,\omega )S(u - {t^{\rm{T}}}(t,\omega )){\rm{d}}t，$

(21)

经过足够多的迭代，可以得出

$ \mathop {\lim }\limits_{N \to \infty } Ts(u,\omega ) = \hat s(\omega )\hat g(0)\delta \left( {u + {\varphi ^\prime }(\omega )} \right) 。$

(22)

由式(22)可知，在经过足够次数的迭代后，即使是较强的频变信号，式(21)的TF能量也可以有效地压缩到GD轨迹中。由于该算法采用多个不动点迭代，因此该技术被称为时间重分配多同步压缩变换(TMSST)。

2）TMSST的信号重建能力

计算TMSST结果沿时间方向的积分，即

$ \begin{split} \int_{ - \infty }^{ + \infty } T {s^{[N]}}(u,\omega ){\rm{d}}u = & \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } G } (t,\omega )\delta \left( {u - {{\hat t}^{[N]}}(t,\omega )} \right){\rm{d}}t{\rm{d}}u =\\ & \int_{ - \infty }^{ + \infty } G (t,\omega ){\rm{d}}t = \hat g(0)\hat s(\omega )。\\[-15pt] \end{split} $

(23)

因此，信号在频域上

$ \hat s(\omega ) = \hat g{(0)^{ - 1}}\int_{ - \infty }^{ + \infty } T {s^{[N]}}(u,\omega ){\rm{d}}u，$

(24)

在时域内的恢复可以使用

$ \begin{split} s(t) = &{(2\text{π} )^{ - 1}}\int_{ - \infty }^{ + \infty } {\hat s} (\omega ){e^{i\omega t}}{\rm{d}}\omega =\\ &{(2\text{π} \hat g(0))^{ - 1}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } T } {s^{[N]}}(u,\omega ){e^{i\omega t}}{\rm{d}}u{\rm{d}}\omega 。\end{split} $

(25)

对于迭代算法来说，终止条件非常重要。考虑到随着迭代次数的增加，2个连续迭代的二维GD轨迹估计将越来越接近，即

$ \left| {{{\hat t}^{[N]}}(t,\omega ) - {{\hat t}^{[N - 1]}}(t,\omega )} \right| = \left| {\frac{{{\sigma ^2}\left( {t + {\varphi ^\prime }(\omega )} \right)}}{{{\sigma ^2} + {\varphi ^{\prime \prime }}{{(\omega )}^2}}}{{\left( {\frac{{{\varphi ^{\prime \prime }}{{(\omega )}^2}}}{{{\sigma ^2} + {\varphi ^{\prime \prime }}{{(\omega )}^2}}}} \right)}^{N - 1}}} \right| ，$

(26)

因此，当满足下列不等式时，迭代可以终止:

$ \iint {\left| {{{\hat t}^{[N]}}(t,\omega ) - {{\hat t}^{[N - 1]}}(t,\omega )} \right|}{\rm{d}}t{\rm{d}}\omega /{L^2} < \lambda。$

(27)

其中 $\lambda $ 为小阈值， $L$ 为离散样本的长度。

TMSST技术的离散实现可以分别通过让 $G[h,m]$ 和 $\hat t[h,m]$ 表示离散STFT结果 $G(t,\omega )$ 和二维GD估计值 $\hat t(t,\omega )$ 来实现。假设所分析的信号的数据长度为偶数 $L$ ，当使用长度为 $M$ 的窗函数来分析信号时，离散 $G[h,m]$ 和 $\hat t[h,m]$ 的大小应该是一个 $L \times M$ 二维矩阵，其中 $M \leqslant L$ 。

3）冲击特征的提取

考虑到冲击往往具有宽频带性质，时频域应存在一个幅值最显著的频率点。该频率点可以用来描述冲击之间的间隔。基于此，TMSST表示冲击频率点的包络谱的最大值为：

$ TFES(\omega )=\max\left|{\int }_{-\infty }^{+\infty }\left(\left|T{s}^{[N]}(t,\omega )\right|-\varphi (\omega )\right){e}^{-i\xi t}{\rm{d}}t\right|，$

(28)

式中 $\phi (\omega )$ 为频率点 $\omega $ 中TMSST结果的平均值。

式(28)可以表示TMSST结果具有最显著冲击特征的频率点。因此，式(28)给出的表达式称为时频包络谱(TFES)。可以用TMSST结果中TFES值最大的频率点来表示冲击类型故障的冲击特征，即 $\left| {T{s^{[N]}}(t,\mathop {\arg \max }\limits_\omega ({TFES} (\omega )))} \right|$ 。

4）时频掩蔽的信号重建

由式(24)和式(25)可知，利用傅里叶反变换可以精确地重构式(1)中给出的频率变化信号，并恢复时域波形。然而，当信号受到背景噪声和其他意外来源的污染时，也有必要将与机械故障密切相关的冲击分量从其他成分中分离出来。从TFES分析可知，振幅最高的频段包含最显著的冲击特征。这就导致了基于TFES频谱在频率方向上设计二值化时频掩蔽(binarized time-frequence masking，TFM)。在频率方向上的TFM (TFM_F)以下式计算：

$ TF{M_ - }F(t,\omega ) = \left\{ {\begin{array}{*{20}{l}} {1,\quad {\text{ if }}\quad {TFES} (\omega ) > \alpha } ，\\ {0,\quad {\text{ otherwise }}} 。\end{array}} \right. $

(29)

其中， $\alpha $ 为TFES的平均值。从式(29)可以得到一个与TFM_F绑定的时频分析，即 $TF{M_ - }F(t,\omega )\left| {T{s^{[N]}}(t,\omega )} \right|$ 。2个连续冲击之间的信号视为背景噪声也应去除。为此，以式(30)获取时间方向的TFM(TFM_T)：

$ \begin{split}&TF{M}_-T(t,\omega )=\\ &\left\{\begin{array}{*{20}{l}}1,&{\rm{if}}\;{\displaystyle \int T}F{M}_-F(t,\omega )\left|{T}^{[N]}(t,\omega )\right|{\rm{d}}\omega > ，\beta\\ 0,&{\rm{otherwise}} 。\end{array}\right.\end{split} $

(30)

其中， $\ \beta $ 为 $\int T F{M_ - }F(t,\omega )\left| {T{s^{[N]}}(t,\omega )} \right|{\rm{d}}\omega $ 的平均值。

以下式计算TFM：

$ {TFM} (t,\omega ) = TF{M_ - }T(t,\omega ) \cdot TF{M_ - }F(t,\omega ) ，$

(31)

去除噪声后的时间序列冲击信号为：

$ s(t) = {(2\text{π} \hat g(0))^{ - 1}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } T } {s^{[N]}}(u,\omega )TFM(t,\omega ){e^{i\omega t}}{\rm{d}}u{\rm{d}}\omega。$

(32)

从式(28)～式(31)得知TFM的构造仅基于TMSST结果，意味着该方法可以在不进行任何人工干预的情况下实现信号的重构过程，因此适合于冲击故障诊断。

2 试验分析

某船用螺杆滑油泵振动异常，该设备为三螺杆泵，额度转速1 400 r/min，额定功率1.1 kW。设备的时域波形和STFT分析结果见图1。从图中可见时域波形中存在密集出现的冲击信号，由于分析时间较短，STFT分析能够观察到冲击产生的谱线，但未能对冲击的发生时刻进行定位。

TMSST方法的分析结果如图2所示，放大约6倍的重见结果如图3所示。图中可见一组清晰的时频谱线。为分析该组谱线是否与设备部件运动产生的频率有关，对TMSST分析结果进行频域、时域切片，使分析结果频率峰值时频分布特性呈现于图4。图4（a）分图为包络谱，显示谱线沿频率的分布情况，可见TMSST方法能够较好得到冲击信号能量的频域分布特性。观察到幅值峰值对应的频率，243 Hz和725 Hz等，均不是螺杆泵部件的故障特征频率。图4（b）是各个冲击在时间轴上的定位，图4（c）左右图是图4（b）的放大，可以在图上测量各冲击的时间间隔。显然时间间隔非定值，说明冲击不是设备部件运动产生。考虑到该设备是流体设备，判断该故障应是流体运动产生，即气蚀。

对螺杆滑油泵及其附件进行拆检，发现滑油进口管路滤器堵塞。判断是该故障引起泵进口真空度增加，压力低于滑油气化压力而引起了气蚀，并激发冲击信号。清洗滤器后故障消失。

为对比TMSST方法在低信噪比情况下的应用效果，对信号添加高斯噪声，噪声水平为SNR= 0 dB，实际上若开展信号直接处理通常要求具备SNR>6 dB。添加噪声后的时域波形和STFT分析结果见图5，可见时域波形中冲击成分已不甚明显。

运用TMSST方法处理加噪信号的时频结果、冲击信号重建、频率峰值时频分布特性分别如图6～图8所示。观察各图并对比图2～图4可知，低信噪比情况下，运用TMSST方法仍能得到冲击信号的清晰谱线，频率峰值时频分布特性仍能反映冲击能量的特点，冲击信号仍较好重建。

分别利用基于Hilbert变换的包络谱分析方法处理添加噪声前后的信号，结果如图9所示。显然，添加噪声前，方法得到了较为清晰的包络谱，1阶频率为243 Hz，与TMSST分析结果一致，但添加噪声后方法已失效。说明该方法难以胜任低信噪比冲击信号的处理需求。而TMSST方法基本不受影响。

3 结　语

TMSST方法采用多个不动点迭代进行群时延轨迹估计，使得方法适应于强频变信号。针对冲击信号的抽取和重建需求进行了频点筛选和二值化时频掩蔽，充分考虑了冲击信号特点和分析需求。对船用螺杆滑油泵气蚀故障的实验分析结果表明，该方法在低信噪比的情况下，仍能清晰提取冲击信号的有关特征支撑开展冲击信号的故障诊断，具有较好的工程实际应用前景。