基于二维二阶同步压缩变换的非平稳信号处理方法

技术领域

本发明涉及信号处理技术领域，具体涉及一种基于二维二阶同步压缩变换的非平稳信号处理方法。

背景技术

针对于音频信号、医疗信号、振动信号等非平稳多分量信号，同步压缩变换是一种增强时频表征的后处理算法，能够将时频系数从初始位置转移到信号能量沿频率轴分布的重心，不仅能够提高信号时频图的分辨率进而提高可读性，还可以通过时频结果重构信号。由于同步压缩变换受到信号分量频率调制的限制，该方法仅适用于纯谐波信号，然而生产生活中的信号往往是复杂且掺杂高噪声的，同步压缩变换应用在复杂的实际信号上往往不能产生集中的时频表示。

因此，如何通过改进算法提高同步压缩变换应用在多类型信号中的时频聚集性，以及提高同步压缩变换抗干扰性是一项挑战。

改进同步压缩变换的关键问题是为真实的瞬时频率提供无偏估计，原始的同步压缩变换方法仅通过频率重分配算子实现，一方面无法应用于如类脉冲信号的其他复杂类型信号，另一方面不精确的重分配算子导致了瞬时频率的有偏估计。因此本发明提出了二阶近似的瞬时频率估计方法，并将该估计方法分别应用于频率和时间两个维度，并与重分配算法相结合，提出二维二阶同步压缩变换。

发明内容

本发明为解决现有同步压缩变换方法通过频率重分配算子实现，存在无未能应用与复杂类型信号，并导致存在瞬时频率的有偏估计等问题，提供一种基于二维二阶同步压缩变换的非平稳信号处理方法。

基于二维二阶同步压缩变换的非平稳信号处理方法，该方法由以下步骤实现：

步骤一、计算信号的短时傅里叶变换得到信号时频特征；

步骤二、根据步骤一获得的信号时频特征，构造二阶时间重分配算子及二阶频率重分配算子估计二维瞬时频率，获得瞬时频率估计值；

步骤三、根据步骤二的获得的瞬时频率估计值，对短时傅里叶时频特征进行重分配；分别获得时间维度上的瞬时频率估计值和频率维度上的瞬时频率估计值，完成对非平稳信号的处理。

本发明的有益效果：本发明所述的非平稳信号处理方法，以非平稳信号的短时傅里叶时频表示为基础，通过构造重分配算子计算时域中幅度和相位的二阶近似得到精确的瞬时频率估计值，并在时间和频率双维度上执行重分配得到清晰的时频表示，更适用于描述强时变和变频信号。该方法可以准确反映瞬时频率随时间变化的趋势，在低信噪比环境下抗干扰能力强，不仅能够对复杂的多分量信号进行高分辨率分解，而且能够在高噪声情况下中提取主要特征。

本发明所述的非平稳信号处理方法，能够更好地消除噪声对信号的影响，得到更精确的多分量信号分解结果。引入定义的二阶重分配算子能够近似得到瞬时频率的无偏估计，增强了信号的分析精度，能够将多分量信号准确分解成不同尺度下的中频成分，从而更清晰的揭示信号不同模态的时频特征，并且该发明对于噪声具有很好的鲁棒性，更容易从受噪声影响的信号中准确提取模态，以便于更好地预测、检测、调控信号。

本发明所述的非平稳信号处理方法，对信号的普适性强，实际生产生活中的信号成分复杂，且受噪声影响严重，目前已存在的重分配变换方法只能处理单一类型信号，如同步压缩变换仅针对类谐波信号有效，相较于已有的方法，本发明可以用于处理多种类型的复杂信号，可应用的领域更加广泛。

附图说明

图1为本发明所述的基于二维二阶同步压缩变换的非平稳信号处理方法的总体流程图；

图2为信号f

图3为信号F

图4为信号s(t)的时频分布效果图；(a)为短时傅里叶变换的时频效果图，(b)为(a)的局部放大图，(c)为同步压缩变换的时频效果图，(d)为(c)的局部放大图，(e)为短时傅里叶变换的时频效果图，(f)为(e)的局部放大图；

图5为同步压缩变换与二维二阶同步压缩变换信号重构的结果对比图；其中，(a)为同步压缩变换的瞬时频率检测结果图，(b)为二维二阶同步压缩变换的瞬时频率检测结果图，(c)为应用同步压缩变换后信号分量一的重构结果图，(d)为应用二维二阶同步压缩变换后信号分量一的重构结果图，(e)为应用同步压缩变换后信号分量二的重构结果图，(f)为应用二维二阶同步压缩变换后信号分量二的重构结果图，(g)为应用同步压缩变换后整体信号的重构结果图，(h)为应用二维二阶同步压缩变换后整体信号的重构结果图；

图6为二维二阶同步压缩变换与另外四种信号处理算法去噪性能的对比图；

图7为振动信号采集环境示意图；

图8为采集信号x

图9为采集信号x

图10为采集信号x

图11为采集信号x

图12为采用本发明所述的获得二维二阶同步压缩变换信号重构的效果图。

具体实施方式

具体实施方式一、结合图1说明本实施方式，基于二维二阶同步压缩变换的非平稳信号处理方法，该方法基于同步压缩变换提出二维二阶同步压缩变换改进算法，以短时傅里叶变换为框架，提出了二阶近似的瞬时频率估计方法，并将该估计方法分别应用于频率和时间两个维度，并与重分配算法相结合，提出二维二阶同步压缩变换。该方法由以下步骤实现：

步骤一、计算信号的短时傅里叶变换得到时频特征；

构造时间变化信号为：

当窗函数为

其中t和τ代表时间及其变化量，ω代表频率，f(t)代表随时间变化的信号。由式(1)可得信号f(t)进行短时傅里叶变换得到的时频特征

其中g(t)和g(ξ)为窗函数，

步骤二、构造二阶时间、频率重分配算子估计二维瞬时频率；具体过程为：

对于非平稳信号，已知瞬时频率的估计值为相位的一阶导数，由此可推到出初始瞬时频率估计值为：

其中

根据短时傅里叶时频特征的能量分布定义时间重分配算子，由Parseval定律可得能量分布

其中f(t)为随时间变化的信号、，

定义

其中

构造二阶时间重分配算子

其中

由二阶时间重分配算子类比推导可得二阶频率重分配算子，针对一个随时间变化的信号f(t)，由式(1)可得信号f(t)进行短时傅里叶变换得到的时频特征

其中，

构造二阶频率重分配算子

其中

引入二阶频率重分配算子后的瞬时频率的无偏估计如下：

步骤3：获取二维二阶同步压缩变换的时频特征；

同步压缩变换(Synchro squeezing Transform,SST)中定义根据瞬时频率估计值对短时傅里叶时频特征进行重分配的结果如下：

其中τ为时间，ξ和ω代表频率及其变化量，

将式(3)和(8)代入式(12)中可以得到时间维度上二阶同步压缩变换表达式：

其中t和τ代表时间及其变化量，ω代表频率，

将式(2)和(11)代入式(12)中可以得到频率维度上二阶同步压缩变换表达式：

其中τ为时间，ξ和ω代表频率及其变化量，

将时间和频率的一维同步压缩变换融合并拓展到二维，给出二维二阶同步压缩变换的定义：

其中t和τ代表时间及其变化量，ξ和ω代表频率及其变化量，

本实施方式所述的方法可以准确反映瞬时频率随时间变化的趋势，在低信噪比环境下抗干扰能力强，不仅能够对复杂的多分量信号进行高分辨率分解，而且能够在高噪声情况下中提取主要特征。本发明可以用于处理多种类型的复杂信号，可应用的领域更加广泛。

具体实施方式二、结合图2至图12说明本实施方式，本实施方式为具体实施方式一所述的基于二维二阶同步压缩变换的非平稳信号处理方法的实施例：本实例首先分别构造了类谐波信号、类脉冲信号、加噪多分量信号等模拟信号来验证本实施方式所提出的方法在提取信号主要成分、去噪等方面的性能，之后为了证明普适性与实用性，本实例将实施方式应用在采集到的振动信号中，并在低信噪比环境下提取振动信号主要成分。

步骤1：定义类谐波信号如下式所示：

f

其中，采样频率设置为100Hz，窗函数采用

由于海森堡不确定原理，短时傅里叶变换的结果非常模糊，如图2中的(a)-(b)。同步压缩变换虽然一定程度上改善了时频能量的聚集，但在处理时频垂直于时间轴的信号部分会产生严重的拖尾现象，如图2中的(c)-(d)。二维二阶同步压缩变换结果如图2中的(e)-(f)所示。与上述方法相比，本实施方式提出的方法通过定义重分配时频系数在时域和频域进行二维压缩，提供了最集中的时频聚集结果，极大地增强了信号全范围内的时频聚合能力。

为了定量描述各种时频处理方法的时频聚集能力，本实施方式采用Rényi熵进行评价，Rényi熵与时频聚集程度呈负相关，时频分布越集中，Rényi熵越小，不同时频处理方法计算得到的Rényi熵结果如表1所示，从表1中可以看出二维二阶同步压缩变换得到的Rényi熵值最小，间接证明了二维二阶同步压缩变换的时频能量聚集性最好。

表1

步骤2：本实施方式中定义类脉冲信号，如下式所示：

F

其中，采样频率设置为600Hz，窗函数采用

从图3的(a)-(d)中可以看出短时傅里叶变换和同步压缩变换很难给出聚集的时频表示结果，特别是针对信号中近似垂直于频率轴的部分，二维二阶同步压缩变换相较于同步压缩变换更适用于类脉冲信号，如图3的(e)-(f)，通过引入重分配算子的方式弥补了时频聚集上的不足，一方面提高了信号模态集中和能量聚集能力，另一方面提高了重分配算法应用于信号的适用性。

类脉冲信号应用不同时频处理方法计算得到的Rényi熵结果如表2所示，从表2中可以看出二维二阶同步压缩变换得到的Rényi熵值最小，进一步证明了二维二阶同步压缩变换的时频能量聚集性最好。

表2

步骤3：本实例定义的加噪多分量信号由以下信号加入白噪声构造：

s(t)＝sin(2π(18(t-0.5)+4sin(2.5(t-0.5))))+sin(2π(40t+sin(3t))) (18)

其中，采样频率设置为100Hz，窗函数采用

在图4的(a)-(b)中，黄绿色的轨迹是信号的时频能量分布轨迹，信号周围分布的是背景噪声，从图4(a)中可以看出，短时傅里叶变换受到噪声的严重影响。加入白噪声后只能粗略描述能量集中的轨迹。同步压缩变换的结果不可避免地受到噪声的影响，导致时频映射非常模糊，如图的4(c)-(d)。二维二阶同步压缩变换可以显示更清晰的时频信息，如图4的(e)-(f)，因为它通过改进中频估计器提高了近似瞬时频率的能力。此外，经过二维二阶同步压缩变换处理后，能量沿着有用信号的时频能量分布轨迹聚集，且受到的背景噪声影响较小。间接证明了本文算法的抗干扰能力。

短时傅里叶变换不能够从时频图中恢复信号，本文的方法可以从时频结果中恢复信号，提取信号的瞬时频率，根据瞬时频率对信号进行重构。

为了从时频结果中分解和重建两种模式的信号，首先检测瞬时频率轨迹。然后根据瞬时频率轨迹对两个模态信号进行分解。最后，对分解后的模态进行叠加，完成对原始信号的重构。整个过程如图5所示。本方法的分解重建结果与同步压缩变换进行了比较。实例表明，两种方法都能对多分量信号进行分解和重构。然而，二维二阶同步压缩变换在提取瞬时频率和重建噪声环境下的信号方面具有更好的精度。

为了评价二维二阶同步压缩变换的去噪能力，本实施方式构造了基于等式(13)的信噪比范围为-10dB到10dB的加噪信号，并将本文方法的去噪性能与目前流行的几种信号处理方法进行了比较，结果如图6所示。显然，该方法在提高信噪比方面比传统的信号处理方法更加突出，特别是在低信噪比(-10dB～-5dB)的情况下。

由此，综合考虑现有时频分析方法的时频聚集性与去噪性能，认为本发明所提出的二维二阶同步压缩变换具有较好的表现。

步骤4：将本实施方式所述的方法应用在振动信号上，该信号是供水管道发生泄漏时所产生的振动信号。当供水管道发生泄漏时会产生流固耦合现象，进而带动管道外壁发生强烈振动，采集到的振动信号可以用来对泄漏点进行异常诊断和定位，然而采集信号由泄露引起的振动和大量背景噪声两部分组成，因此如何从采集信号中提取并恢复由泄露引起的振动信号是关键。

图7中S

传感器S

将短时傅里叶变换、同步压缩变换以及本实例中的实施方法分别应用x

通过计算Rényi熵来定量反映三种方法的时频聚集能力，如表3所示，从表3中可以看出，本实例的实施方式Rényi熵最小，证明了时频聚集能力最佳。

表3

假设采集信号x

以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。