时频分析是处理非平稳随机信号强有力的工具，可以提供时变信号的时间-频率联合域分布，有利于了解信号瞬时频率随时间变化的规律特点^[1]。目前时频分析方法已经广泛应用于通信、地质、航空、勘探、工业生产、目标检测等各个领域^[2-3]。在人造金刚石合成加工过程中，顶压机顶锤破裂是一种常见的生产事故，若操作过程中没有及时预警发现，极易带来巨大的经济损失和安全隐患^[3]。金属破裂信号属于非平稳随机信号，应用时频分析方法在线监测加工过程，提供破裂判决依据，对于生产安全有着重要的意义。

20世纪中叶以来，随着各种时频分析方法的提出、发展与演化，逐渐形成较为完善的算法体系^[4]。目前广泛应用的方法有短时傅里叶变换(short-time Fourier transform, STFT)^[5]、连续小波变换(continuous wavelet transform, CWT)^[6]、Wigner-Ville分布(Wigner-Ville distribution, WVD)^[7]、S变换^[8]、希尔伯特-黄变换(Hilbert-Huang transformation, HHT)^[9]以及同步压缩小波变换(synchrosqueezing transform, SST)^[10]等。1946年，Gabor首先提出了STFT，但是STFT受到固定窗函数的影响，在分析过程中无法兼顾高频分量和低频分量^[5-11]。1948年，物理学家Ville将Wigner分布应用于信号处理领域，形成了WVD。不过WVD对于多分量信号的时频分析而言，存在严重的交叉干扰项^[7]。20世纪80年代，Morlet和Grossman等学者共同发展了小波变换，这是一种多频率分辨率的分析方法^[6]。小波变换有效地解决了频率分辨率和时间分辨率的矛盾，然而由于每种小波基影响的不确定性，增加了实际应用中的难度。1996年，Stockwell等发现S变换的逆变换易于实现，提出了S变换^[8]。1998年，Huang等人在经验模态分解方法的基础上引入希尔伯特谱的概念和分析方法，出现了希尔伯特-黄变换(HHT)^[9]。基于EMD的HHT方法能够较好地反映信号的真实物理过程，具有完全自适应性。HHT的数学支撑理论尚不完善，并且存在端点效应和模态混淆等问题^[9]。受不确定性原理或窗函数的影响，这些传统方法易产生时频模糊或交叉干扰等现象，限制了它们在信号处理领域的发展^[12]。随着研究的深入出现了一些新兴时频分析方法，2011年，Daubechies提出了一种时频分析后处理方法，即同步压缩变换(SST)。SST是一种时频重排技术，其变换核可以基于小波变换或者基于STFT变换，SST支持信号重构，其时频输出具有高锐化聚集度的特点^[10]。Lin等人提出将倒谱与非线性时频分析方法结合的De-shape SST算法^[13]，该算法可以减轻波形函数对时频分布的干扰影响。随着信号理论的发展，为了更灵活、自适应和简地的表示信号，Mallat等人于1993年提出在过完备原子库上分解信号的思想^[14]。常用的稀疏分解算法有匹配追踪(matching pursuit, MP)^[14]及其各种变异优化的改进算法。

De-shape SST算法在时间-频率平面上有良好的时频聚集度，但其抗噪性能及鲁棒性较差。文献[15]提到，对于信噪比低的信号，传统降噪算法不能获得令人满意的降噪效果，稀疏分解理论根据是否为稀疏成分，将信号和噪声区分开来。对高强度噪声的去噪效果表现更为优良。文献[16]也提到稀疏分解在图像领域的良好的降噪效果。文献[17]中提出的NMP(nonlinear matching pursuit)算法，与传统MP算法相比计算复杂度降低，算法适用性更广，稀疏分解效果更加优异。本文分析的金属破裂信号在采集过程中存在着大量的背景噪声，鉴于此，本文选择文献[17]中提出的NMP稀疏降噪重构算法作为本文去噪算法。利用NMP良好的稀疏能力，结合De-shape SST算法优越的时频锐化性能，实现基于非线性匹配追踪的De-shape SST算法(nonlinear matching pursuit De-shape synchrosqueezing transform, NDSST)。通过数值仿真实验说明，改进后的算法可以得到信号准确的时频分布，抗噪声性能良好。金属破裂信号在采集过程中不可避免会掺杂环境噪声，尝试利用NDSST算法分析金属破裂样本数据，通过稀疏重构抑噪后，得到信号的时间-频率分布特点。

1 De-shape SST算法

对于某类实际信号，通常会选择相应的数学模型来模拟该信号，本文中采用自适应非谐波(adaptive non-harmonic，ANH)模型量化实际待处理的动态系统^[11]。针对这一系统的信号，De-shape SST算法可以提供清晰的频率随时间即时变化分布，倒谱的应用可以消除由非平凡波形函数在信号的时频表示中产生的影响^{[13, 18-19]}。

1.1 短时倒谱变换

倒谱的传统概念是将信号的频谱取对数，再将结果进行傅里叶变换，得到信号在倒频域的倒谱。传统倒谱思想存在一定的局限性，例如，复倒谱需要相位解缠过程，实倒谱受对数运算的限制等。De-shape SST算法中采用广义化根倒谱，广义化倒谱和根倒谱的定义如式(1)和式(2)所示^{[13, 20-21]}:

$ g_{\gamma}(x)=|x|^{\gamma} $

(1)

$ \tilde{f}_{\gamma}(q)=\int g_{\gamma}(\hat{f}(\xi)) \mathrm{e}^{2 \mathtt{π} i q \xi} \mathrm{d} \xi $

(2)

式中，参数γ需要根据不同的实际信号选择适合该信号的值。当γ=2时，式(1)和式(2)等效于信号的自相关函数，不适用于多尺度信号^[13]。为解决上述问题，在本文中取γ < 2。参考倒谱测量“每赫兹的谐波峰值的速率”的思想，求解信号的周期，即倒谱中的突出峰值指示信号的周期和它们的倍数^[13]。将倒谱变换与短时傅里叶变换结合便可以得到短时倒谱变换(short-time cepstral transform，STCT)。STFT的定义为：

$ V_{f}^{(h)}(t, \xi)=\int f(\tau) h(\tau-t) \mathrm{e}^{-i 2 \mathtt{π} \xi(\tau-t)} \mathrm{d} \tau $

(3)

式中，ξ为信号的瞬时频率；f(τ)为待分析信号；h为窗函数。STCT的定义为^[13]：

$ C_{f}^{(h, \gamma)}(t, q)=\int g_{\gamma}\left(V_{f}^{(h)}(t, \xi)\right) \mathrm{e}^{-i 2 \mathtt{π}q \xi} \mathrm{d} \xi $

(4)

式中，C_f^{(h, γ)}(t, q)中的q称作倒频，它以秒或者其他时域中的单位作为计量单位。显然C_f^{(h, γ)}(t, ·)是信号f(·)h(·-t)的倒频谱，且通常C_f^{(h, γ)}(t, q)不总是正值。

C_f^{(h, γ)}(t, q)是倒频轴上的连续函数，倒谱能够提供待分析信号的周期信息，STCT可以实现信号的时间-周期表示^[13]。基于周期和频率的关系，在STCT中的基频周期的倒数是该ANH函数的基本频率，即在STCT变换中将信号的瞬时频率信息镶嵌于瞬时周期信息之中^[13]。

根据对STCT的讨论和谐波约束原则^[13]的引导，得到逆短时倒谱变换(inverse short-time cepstral transform，iSTCT)，定义为^[13]：

$ \left\{\begin{array}{l}{\psi : q \rightarrow 1 / q} \\ {U_{f}^{(h, \gamma)}(t, \xi)=C_{f}^{(h, \gamma)}(t, \psi \xi)}\end{array}\right. $

(5)

式中，算子ψ操作翻转时间-周期表示的倒频轴，当q>0时，通过频率和周期的倒数关系，可以在iSTCT过程中获得信号的频率信息。

1.2 De-shape SST

De-shape STFT算法定义为iSTCT与STFT结合，用来提取信号的基本频率信息^[13]，如式(6)所示：

$ W_{f}^{(h, \gamma)}(t, \xi)=V_{f}^{(h)}(t, \xi) U_{f}^{(h, \gamma)}(t, \xi) $

(6)

De-shape STFT可以减轻波形函数的影响，但依然受限于不确定性原理。经过上述过程得到的时频平面存在时频模糊、频率带过宽的局限性。针对此问题，应用SST算法锐化De-shape STFT的输出矩阵，时频重排规则^{[10, 18, 22]}如下：

$ \mathit{\Omega}_{f}^{(h, \nu)}(t, \xi)=\\ \left\{\begin{array}{ll}{-\mathfrak{I} \frac{V_{f}^{(Dh)}(t, \xi)}{2 \mathtt{π} V_{f}^{(h)}(t, \xi)}, } & {\left|V_{f}^{(h)}(t, \xi)\right|>v} \\ {-\infty, } & {\left|V_{f}^{(h)}(t, \xi)\right| \leqslant v}\end{array}\right. $

(7)

式中，Dh代表窗函数的导数；$\mathfrak{I}$代表虚数部分，将v>0作为阈值，可避免因|V_f^(h)(t, ξ)|太小而造成的数值不稳定性。由此可得De-shape SST的定义^[13]如下：

$ \begin{gathered} SW_f^{\left( {h,\gamma ,v} \right)}\left( {t,\xi } \right) = \int {_{{\mu _{v(t)}}}W_f^{\left( {h,\gamma } \right)}\left( {t,\eta } \right)\frac{1}{\alpha }g \cdot } \hfill \\ \left( {\frac{{\left| {\xi - \mathit{\Omega} _f^{\left( {h,v} \right)}\left( {t,\eta } \right)} \right|}}{\alpha }} \right){\text{d}}\eta \hfill \\ \end{gathered} $

(8)

De-shape SST算法的主要优点是利用时频重排技术，抑制波形函数对信号的时间和频率轴的影响，得到准确清晰的时频表示^[13]。

2 NDSST算法

De-shape SST在算法原理上可以消除波形函数运算过程中带来的高次谐波影响，同时也能消除部分噪声干扰，但De-shape SST算法本身没有除噪过程，其时频分布结果受噪声影响的大小存在不确定性。当信号的信噪比(signal noise ratio, SNR)较低时，De-shape SST无法得到准确清晰的信号的时频表示。针对De-shape SST的上述局限性，本文提出NMP与De-shape SST相结合的强鲁棒性和时频高锐化算法。NMP算法的主要作用是对信号进行稀疏分解，去除信号中的噪声，得到信号的稀疏表达式，用最少的基函数逼近原始信号本质信息^[17]。

NMP算法的基本思想是建立过完备字典并实现逐点迭代。具体步骤为：①首先建立过完备的原子字典，比较残差信号和原子之间的相似度，然后找出表示最小误差的目标原子；②去掉旧残差信号的投影以获得新的残差信号；③重复上述过程，直到分解结果的误差足够小或者满足其他收敛条件，便停止迭代，完成整个稀疏分解过程^[13]。这里的迭代次数根据信号分解的需求设置，迭代次数越多，分解的精度就越高，但是耗时也就越长。NMP迭代原子过程如式(9)至式(12)所示^{[17, 23]}。

首先建立一个过完备字典D：

$ \begin{array}{c}{D=\{a(t) \cos \varphi(t) | a(t) \in V(\varphi)} ,\\ {a(t)>0, \frac{\mathrm{d} \varphi(t)}{\mathrm{d} t}>0 \}}\end{array} $

(9)

式中, V(φ)是所有比cosφ(t)更平滑的函数集合，对于一个给定的φ(t)，V(φ)是一组基于傅里叶变换的过完备线性空间：

$ V(\varphi)=\operatorname{span}\left\{1, \cos \left(\frac{k \varphi}{2 L_{\varphi}}\right), \sin \left(\frac{k \varphi}{2 L_{\varphi}}\right), k=1 \cdots 2 \lambda L_{\varphi}\right\} $

(10)

式中，$L_{\varphi}=\left[\frac{\varphi(1)-\varphi(0)}{2 \mathtt{π}}\right] $；λ是控制函数V(φ)平滑度的参数，那么：

$ D=\left\{a(t) \cos \varphi(t) : \varphi^{\prime}(t) \geqslant 0, a(t), \varphi^{\prime}(t) \in V(\varphi)\right\} $

(11)

然后进行初始化r⁰=f(t)，解决正则非线性最小二乘问题：

$ P : \gamma\left\|\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a}_{k}\right\|_{l^{1}}+\left\|r^{k-1}-a_{k} \cos \varphi_{k}\right\|_{l^{1}}^{2}(最小化) $

(12)

式中，γ是正则化参数；$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over a}_{k} $和a_k是V(φ)的空间坐标。

最后更新信号残差，并根据预先设置的迭代停止条件，不断重复上述过程，直到满足条件为止。

NMP的稀疏分解能力强于传统MP，但是对于大数据来说，解决正则化最小二乘问题依然非常耗时。本文采用快速傅里叶变换(fast Fourier transform，FFT)的方法提高其计算速度^[17]，该方法的稀疏表示可以减少存储的数据量，提高压缩效率，且可抑制信号中的噪声。NDSST算法流程如图 1所示。

首先将待分析信号通过NMP算法得到代替原始信号的稀疏表示，该稀疏表示作为De-shape SST算法的输入，通过De-shape SST过程得到信号的时频分布表示。该算法结合二者的优点，适合处理含噪声的非平稳信号，尤其对低SNR信号的去噪效果和时频高锐化描述有良好的表现。

3 数值实验 3.1 理想信号时频分析实验

选择4种不同类型的仿真信号来检验NDSST算法的时频描述效果和抗噪性能，如图 2所示。

图 2(a)为单频信号x=sin(10πt)；图 2(b)为余弦频率信号x=cos(2π(4t+cos(t/2)))；图 2(c)为双分量线性调频信号x=sin(8πt+2.4t²)+sin(16πt+2.4t²)；图 2(d)为交叉频率信号x=sin(t²+t+sint)+sin(8t)。实验的信号类型多样，有恒频、变频、线性调频和组合频率信号，NDSST算法的时频分析表现较好，时频分布的正确性和聚集度的效果达到预期。

3.2 NDSST算法鲁棒性实验

将NDSST算法与De-shape SST、基于SST的STFT算法(STFT-SST)和SWT就算法鲁棒性方面进行对比测试。首先测试4种方法在不同SNR下的时频分布直观效果，再用时频聚集度(concentration measure, CM)对算法性能进行定量评价，得到这4种方法处理同一仿真信号CM值的大小关系，以及CM值随SNR的变化规律。

对测试信号x=cos(2π(4t+cos(t/2)))添加加性噪声，SNR分别设置成-5 dB、-2 dB、0 dB、5 dB。随着SNR的改变，即使在-2 dB和-5 dB的低信噪比条件下，NDSST算法也可得到比较理想的时频表示。

在信号中添加-5 dB、-2 dB、0 dB、5 dB的高斯白噪声，用SWT、STFT-SST、De-shape SST和NDSST 4种算法分别处理测试信号，得到如图 3-6所示的时频分布图。据图 3可知，相同条件下，NDSST的分析结果更加清晰准确，时频分布图中噪声波纹与3种对比方法相比明显减少，其他信噪比条件下的分析结果与上述结论相同。纵向观察同一方法不同信噪比情况下，SWT在低信噪比时难以得到正确的信号时频分布。随着信噪比的增加，STFT-SST和De-shape SST算法能够得到相对准确的时频表示，同时De-shape SST算法的去噪性能优于STFT-SST，时频分布效果有一些改善，但仍不尽人意。低信噪比时，STFT-SST和De-shape SST很难从背景噪声中观察出真实的时频表示。NDSST在低信噪比时，去噪性能良好，能够得到准确易于直观识别的时频分布，随着信噪比的增加，时频分布越来越清晰。

3.3 时频分析聚集度评价

为了能够更好地说明NDSST算法的优越性，采用时频聚集度(CM)对算法性能进行定量分析。CM值定义为^[24]：

$ \mathrm{CM}=\frac{\iint\left|W_{s}(\omega, t)\right|^{4} \mathrm{d} \omega \mathrm{d} t}{\left(\iint |W_{s}\left.(\omega, t)\right|^{2} \mathrm{d} \omega \mathrm{d} t\right)^{2}} $

(13)

式中，W_s(ω, t)是信号最终得到的时频分布；CM值代表了频率分量的峰值能量与总能量之间的归一化比率。CM值越高代表频率聚集度越好，算法的性能越优越。

对测试信号x=cos(2π(4t+cos(t/2)))进行仿真，如图 7所示，是应用§3.2中4种方法分别处理信号，计算得出每种方法随SNR增加CM的变化趋势。实验表明，4种方法的CM值随SNR增加均呈先增加后趋于稳定的变化规律。图 7中红色曲线代表NDSST在不同SNR下的CM值，可以看出在SNR从-15 dB增加至25 dB的过程中，其CM值均大于其他3种方法。在-15 ~0 dB范围内，本文算法的时频聚集度明显高于3种对比算法，说明NDSST的抗噪声性能优良，在低信噪比情况下，时频聚集度依然保持良好，本算法与其他算法相比更加适合于分析含大量背景噪声的实际信号，此结论与§3.2中的实验结果分析相符。

为了进一步说明NDSST与其他算法相比的性能提升效果，本文采用性能提升百分比(percentage increase in performance，PIP)更加直观说明NDSST的强鲁棒性。其计算公式为：

$ \left\{ {\begin{array}{*{20}{c}} {{\rm{PI}}{{\rm{P}}_{{\rm{max}}}} = \frac{{\left( {{\rm{C}}{{\rm{M}}_{{\rm{this}}}} - {\rm{C}}{{\rm{M}}_{\min }}} \right)}}{{{\rm{C}}{{\rm{M}}_{{\rm{min}}}}}} \times 100\% }\\ {{\rm{PI}}{{\rm{P}}_{\min }} = \frac{{\left( {{\rm{C}}{{\rm{M}}_{{\rm{this}}}} - {\rm{C}}{{\rm{M}}_{{\rm{max}}}}} \right)}}{{{\rm{C}}{{\rm{M}}_{{\rm{max}}}}}} \times 100\% } \end{array}} \right. $

(14)

式中，PIP_max表示本文所比较的几种方法中，本文方法相对于其他算法性能提升的最高百分比；PIP_min表示最低百分比；CM_this代表本文方法在各信噪比下的CM值；CM_min表示除本文方法外其他算法在相应信噪比下的最小CM值，CM_max表示相应信噪比下的最大CM值。根据式(14)计算得性能提升百分比如图 8所示。

从图 8中可以看出，在低信噪比(SNR≤0 dB)情况下，本文算法性能提升最低为0 dB时，提升56.25%。信噪比大于0 dB时，通过图 7可以看出本文算法性能已经趋于稳定，并且优于其他3种方法，从图 8中也可以看出性能提升趋于稳定。随着其他方法在SNR提高的同时算法性能不断提升，本文方法的性能提升百分数有所下降，最低提升百分比基本稳定在15%~20%之间，最高提升百分比稳定在60%左右。

4 金属破裂信号时频分析应用

对于非平稳随机信号，不局限于仅从频域或者时域分析，更需要将时域和频域联合起来研究信号的瞬时频率变化特点。金属破裂信号是声发射波的一种，属于宽频带非平稳随机的多模式波^[4]。本文选择4组传感器接收的顶压机顶锤破裂前后的样本数据，应用De-shape SST算法和NDSST算法进行处理，求解时域幅度突增时频率变化的范围。通过对比发现，NDSST算法能够清晰得到信号瞬时频率随时间变化的细节，可为工程实践中设置监测传感器的阈值提供判断依据，如图 9-12所示。

金属破裂信号幅度非常小，即使幅度突增，增大后的信号幅度一般也不超过0.05个单位，因此在生产过程中，很难准确检测出来，如图 9-12所显示的时域信号波形图。以第一组样本信号为例，由图 9(a)可知，信号在采样数据点的第4 000个数据前后开始出现幅度突增，表示金属开始出现破裂。通过De-shape SST算法分析信号得到的时频表示如图 9(b)所示，可以大致看出发生金属破裂的时间，但是在该时刻信号的频率段分布于0~1 000 kHz整个频率轴上并且存在大量杂乱波纹，时频表示能量分散，时频聚集度低。由NDSST算法得到的时频分布图如9(c)所示，可以清晰地显示，在该疑似破裂时刻，时频分布图上出现600~1 000 kHz的大幅度频率分量，此特点为仪器监测示警条件提供判断参考。在图 10-12中，金属发生破裂时间极短，异常信号能量相比第一组样本信号更小，用De-shape SST算法进行处理得到的时频表示不仅能量分散，更是很难观察到突变信号的特点。与De-shape SST算法的分析结果相比，通过NDSST得到的时频分布非常集中，抑制了噪声和非突变信息对破裂信息的影响，便于观察和捕捉金属破裂出现的时间和频率窗口。图 10-12的实验结果均表现出类似规律，说明该改进后的算法对非平稳信号具有良好的时频高锐化表示和抗噪性能。

5 结语

本文提出了一种基于NMP的De-shape SST算法，该算法通过NMP对信号进行稀疏分解以达到去除噪声的目的，通过De-shape SST变换可得到聚集度良好的时频分布。数值实验结果表明，NDSST算法针对非平稳信号有较高的锐化CM值和强鲁棒性。将本文算法应用于金属破裂信号的时频分析，得到较为清晰的信号样本的时间-频率分布，在金属破裂期间出现600~1 000 kHz幅值较大的频率集束，可作为监测仪器示警条件之一。