一种基于差分非线性模式分解的谐波信号检测方法

摘要

本发明公开了一种基于差分非线性模式分解的谐波信号检测方法，这种新的谐波检测方法先将差分算法应用于原始信号，然后利用非线性模式分解算法获得一系列有意义的非线性模式分量；对每个非线性模式分量积分后再次进行非线性模式分解，并提取出第一个非线性模式分量作为原始信号的非线性模式分量；最后，用希尔伯特边际谱提取信号的各种谐波信号。基于差分非线性模式分解的谐波信号检测方法继承了非线性模式分解方法的噪声鲁棒特性，并且能够抑制混沌和噪声信号的干扰，在提取小幅度的高次谐波分量方面有优势。

说明书

技术领域

本发明涉及谱估计和信号处理应用领域，尤其涉及一种基于差分非线性模式分解的谐波信号检测方法。

背景技术

自适应时频分析方法被广泛应用于诸多领域，例如语音信号分析、声呐信号处理和机械故障诊断等。自然现象中存在许多可观测的混沌信号，如海洋坏境噪声、由电磁脉冲在海洋表面产生的杂波等。在混沌背景下对信号进行检测和分析具有一定的难度。

傅里叶变换具有简单和计算高效的优点、是分析谐波最常用的算法。但是，它存在三个主要缺陷：混叠、栅栏效应和频谱泄漏，不适合分析非平稳信号。基于径向基函数神经网络的方法也可用于检测待测信号的谐波幅度。然而，神经网络算法的主要缺点是训练过程复杂、收敛速度慢和容易陷入局部极小值。

希尔伯特-黄变换(Hilbert-Huang Transform，HHT)是一种新的非平稳信号分析方法，采用经验模态分解(Empirical Mode Decomposition，EMD)将复杂信号分解为一系列本征模函数(Intrinsic Mode Function，IMF)，进一步对每一个IMF作Hilbert变换，得到关于原始信号的时间频率联合分布。EMD广泛应用于非线性和非平稳信号分析。但是，在实际应用中EMD会出现模式混叠现象，其表现在两个方面：一是同一个IMF分量中包含了尺度差异较大的信号分量；二是同一个尺度的信号分量出现在了不同的IMF中。为了解决这个问题，给信号加上高斯白噪声并应用EMD，于是，集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)应运而生。差分经验模态分解(Differential Empirical Mode Decomposition，DEMD)首先对原始信号差分运算然后应用EMD，可以提取出EMD和EEMD无法分离的小幅度高次谐波分量。然而，以上三种方法都存在对噪声敏感的缺陷。

非线性模式分解方法(NMD)可将信号分解为一系列物理上有意义的非线性模式分量，同时消除噪声。相对于过去的方法，NMD因为可以自适应地选择参数，所以具有极好的噪声鲁棒特性。但是，NMD在小幅度谐波信号检测时的表现并不是很好，特别是在混沌背景存在的情况下，其表现一般。

发明内容

本发明的目的在于提供一种基于差分非线性模式分解的谐波信号检测方法,该方法综合了现有方法的优势，克服了经验模态分解方法、集合经验模态分解方法和差分经验模态分解方法对混沌和噪声信号敏感以及非线性模式分解方法在小信号谐波检测时性能不好的缺点。

为实现上述目的，本发明采用的技术方案为：一种基于差分非线性模式分解的谐波信号检测方法，包括以下步骤：

步骤A：准备待检测处理的原始信号s(t)，其采样率为f_s，数据长度为N；

步骤B：对原始信号s(t)进行差分运算获得新的信号s′(t)；

步骤C：对信号s′(t)进行非线性模式分解获得非线性模式分量c_i′(t)；

步骤C-1：计算信号s′(t)的小波变换W_s′(ω,t)，小波变换定义为：

其中，是s′(t)的傅立叶变换；s′⁺(t)是信号s′(t)的正频率部分，其表达式为：ψ(t)为小波变换的小波函数，为ψ(t)的傅立叶变换，满足条件上标*表示共轭运算；是小波的峰值频率；小波函数采用对数正态分布小波，其和ω_ψ表示如下：

其中，f₀是权衡变换过程中时间和频率分辨率的分辨率参数，通常默认f₀＝1；

步骤C-2：检查小波变换是否为最佳时频表示，若不是，则采用加窗傅里叶变换G_s′(ω,t)，加窗傅里叶变换定义如下：

其中，g(t)是加窗傅里叶变换的窗函数，为g(t)的傅立叶变换，满足条件：选择高斯窗作为加窗傅里叶变换的窗函数，其表达式为：

步骤C-3：找出信号s′(t)的时频表示中所有的脊曲线这里定义是h次谐波的脊曲线；

在某个时刻t_n，利用下式算法，可找出h个极大值点；

上式中n＝1,2,…,N；N为数据长度；H_s′(ω,t)是经过上述傅里叶变换或加窗傅里叶变换后的时频表示，即为W_s′(ω,t)或G_s′(ω,t)；

将H_s′(ω,t)中找出的所有时刻的脊点连线，可构成h条脊曲线

步骤C-4：利用脊曲线重构h次谐波分量其中A^(h)(t)、和分别是其幅度、相位以及频率；计算过程如下：

若信号s′(t)的时频表示采用傅里叶变换，则第h次谐波分量x^(h)(t)由公式(1)得到：

若信号s′(t)的时频表示采用加窗傅里叶变换，则第h次谐波分量x^(h)(t)由公式(2)得到：

其中，和分别为傅里叶变换和加窗傅里叶变换由抛物线插值所产生的离散化影响的修正；

步骤C-5：用抗噪性替代检验方法确定有效的谐波分量；该方法利用替代数据鉴别提取出的谐波分量的真伪，筛选出所有的真实的谐波分量，并且当连续三个谐波分量被判断为假时停止分解过程；具体步骤如下：

(1)计算提取出的某一谐波分量的辨识统计量D₀(α_A,α_ν)；

提取出的每一个谐波分量的幅值A^(h)(t)和频率ν^(h)(t)的有序度可以用其谱熵和来定量地衡量，其中，和分别为A^(h)(t)和ν^(h)(t)的傅立叶变换，辨识统计量D定义如下：

其中，α_A和α_v分别为和的权值系数；

(3)为信号s′(t)创建N_s个傅立叶变换替代数据，其产生方法为：

其中，φ_ξ服从[0,2π)上的均匀分布，每一个φ_ξ对应一个傅立叶变换替代数据；

(3)计算与每一个替代数据对应的时频表示，并分别从中提取出各次谐波分量，计算出各个替代数据的辨识统计量这里定义显著性水平指标为：

式中，为满足D_s>D₀的替代数据的个数；假设创建N_s个替代数据并且将显著性水平指标设置为p，即至少有N_s×p个替代数据满足D_s>D₀才认为该分量不是噪声，从而继续分解过程；这个检测使用三组不同值的参数(α_A,α_ν)，即分别计算出D(1,1)，D(0,1)和D(1,0)的值，只要其中至少有一个值不符合无效假设，则认为满足D_s>D₀；

(4)计算谐波之间的相关度的综合度量值

其中，

式中，w_A,w_φ,w_v代表的权值；默认使用ρ^(h)≡ρ^(h)(1,1,0)为幅度和相位一致性分配相等的权值，且对频率一致性不分配权值；

(5)为了减少对真实谐波分量的错误判断，定义综合度量值的阈值为：

(6)当一个谐波分量的综合度量值指标满足ρ^(h)≥ρ_min且显著性水平指标significance-level≥p＝95％时，则认为此谐波分量通过检验，为真实的谐波分量；如果不能通过抗噪性替代检验，则停止非线性模式分解；

步骤C-6：将所有的真实的谐波分量相加构成一个非线性模式分量c₁′(t)；

步骤C-7：从信号s′(t)中减去该非线性模式分量，重复步骤C-1到步骤C-6，得到所有的非线性模式分量c_i′(t)；

步骤D：对得到的非线性模式分量c_i′(t)积分获得b_i(t)；

步骤E：对b_i(t)重新用非线性模式分解方法进行非线性模式分解；

步骤E-1：按步骤C-1到步骤C-6的过程对每个b_i(t)进行非线性模式分解；

步骤E-2：提取出每个b_i(t)的第一个非线性模式分量作为原始信号s(t)的非线性模式分量c_i(t)；最后得到原信号的非线性模式分解结果：

步骤F：对得到的s(t)进行谱分析，提取谐波信号；

步骤F-1：对c_i(t)做希尔伯特变换产生正交分量

步骤F-2：构造复信号z_i(t)，其表达式为：

步骤F-3：将复信号z_i(t)转化为极坐标形式，求得瞬时包络a_i(t)和瞬时频率ω_i(t)；z_i(t)的极坐标形式为：

其中，瞬时包络a_i(t)和瞬时相位φ_i(t)表示如下：

其瞬时频率ω_i(t)为瞬时相位φ_i(t)对时间的导数，表示为：

步骤F-4：求所有非线性模式分量c_i(t)的希尔伯特谱H(ω,t)；

希尔伯特谱定义为：

步骤F-5：求希尔伯特边际谱，提取信号谐波分量的谱峰；

希尔伯特边际谱的表达式为：

其中，T表示采样时间，其计算公式为：T＝N/f_s。

本发明与现有技术相比，其显著优点是：

1.差分非线性模式分解方法能够更有效地抑制噪声和混沌信号所产生的干扰；

2.对于具有多个谐波分量的信号，差分非线性模式分解方法能够提取出其中幅度较小的谐波分量。

附图说明

图1为混沌信号d(t)和原始信号s(t)的时域波形；

图2为SNR＝25dB时采用差分非线性模式分解方法得到的希尔伯特边际谱；

图3为SNR＝25dB时采用非线性模式分解方法得到的希尔伯特边际谱；

图4为SNR＝25dB时采用经验模态分解方法得到的希尔伯特边际谱；

图5为SNR＝25dB时采用集合经验模态分解方法得到的希尔伯特边际谱；

图6为SNR＝25dB时采用差分经验模态分解方法得到的希尔伯特边际谱。

具体实施方式

以下以加有高斯白噪声和混沌背景的原始信号为例，结合附图，详细说明本发明的实施方式。

本发明首先对原始信号s(t)进行差分运算得到新的信号，而后对其进行非线性模式分解得到一组非线性模式分量c_i′(t)；对得到的非线性模式分量进行积分得到b_i(t)；再对其进行非线性模式分解，提取出b_i(t)的第一个非线性模式分量作为原始信号的非线性模式分量，记为c_i(t)；最后，利用谱分析技术在信号处理领域的优点，对最终得到的非线性模式分量进行希尔伯特边际谱分析，实现谐波信号的检测。

假设原始信号s(t)为：

s(t)＝cos(2πft)+0.3cos(2π·3ft)+0.07cos(2π·5ft)+n(t)+d(t)

其中，n(t)为高斯白噪声，混沌信号d(t)＝0.05x(t)。这里，设基波频率f＝30Hz，采样频率f_s＝600Hz，数据长度为N＝650。

x(t)由典型的Lorenz混沌系统产生，其系统模型描述如下：

其中，σ＝10，r＝28，b＝8/3，初始值x₀＝y₀＝z₀＝0.1。d(t)和s(t)时域波形如图1所示。

本发明是一种基于DNMD(Differential Nonlinear Mode Decomposition，差分非线性模式分解)的谐波信号检测方法，包括如下步骤：

步骤A：准备待检测处理的原始信号s(t)，其采样率为f_s＝600Hz，数据长度为N＝650；

步骤B：对原始信号s(t)进行差分运算获得新的信号s′(t)；

步骤C：对信号s′(t)进行非线性模式分解获得非线性模式分量c_i′(t)；

步骤C-1：计算信号s′(t)的小波变换W_s′(ω,t)，小波变换定义为：

其中，是s′(t)的傅立叶变换；s′⁺(t)是信号s′(t)的正频率部分，其表达式为：ψ(t)为小波变换的小波函数，为ψ(t)的傅立叶变换，满足条件上标*表示共轭运算；是小波的峰值频率；小波函数采用对数正态分布小波，其和ω_ψ表示如下：

其中，f₀是权衡变换过程中时间和频率分辨率的分辨率参数，通常默认f₀＝1；

步骤C-2：检查小波变换是否为最佳时频表示，若不是，则采用加窗傅里叶变换G_s′(ω,t)，加窗傅里叶变换定义如下：

其中，g(t)是加窗傅里叶变换的窗函数，为g(t)的傅立叶变换，满足条件：选择高斯窗作为加窗傅里叶变换的窗函数，其表达式为：

步骤C-3：找出信号s′(t)的时频表示中所有的脊曲线所谓的脊曲线就是时频图上一些局部极大值点所连接而成的曲线，各个极大值点就叫脊点；这里定义是h次谐波的脊曲线；

在某个时刻t_n，利用下式算法，可以找出h个极大值点；

上式中n＝1,2,…,N；N为数据长度；H_s′(ω,t)是经过上述傅里叶变换或加窗傅里叶变换后的时频表示，即为W_s′(ω,t)或G_s′(ω,t)；

将H_s′(ω,t)中找出的所有时刻的脊点连线，可以构成h条脊曲线

步骤C-4：利用脊曲线重构h次谐波分量其中A^(h)(t)、和分别是其幅度、相位以及频率；计算过程如下：

若信号s′(t)的时频表示采用傅里叶变换，则第h次谐波分量x^(h)(t)由公式(1)得到：

若信号s′(t)的时频表示采用加窗傅里叶变换，则第h次谐波分量x^(h)(t)由公式(2)得到：

其中，和分别为傅里叶变换和加窗傅里叶变换由抛物线插值所产生的离散化影响的修正；

步骤C-5：用抗噪性替代检验方法确定有效的谐波分量；该方法利用替代数据鉴别提取出的谐波分量的真伪，筛选出所有的真实的谐波分量，并且当连续三个谐波分量被判断为假时停止分解过程；具体步骤如下：

(1)计算提取出的某一谐波分量的辨识统计量D₀(α_A,α_ν)；

提取出的每一个谐波分量的幅值A^(h)(t)和频率ν^(h)(t)的有序度可以用其谱熵和来定量地衡量，其中，和分别为A^(h)(t)和ν^(h)(t)的傅立叶变换，辨识统计量D定义如下：

其中，α_A和α_v分别为和的权值系数；

(4)为信号s′(t)创建N_s个傅立叶变换替代数据，其产生方法为：

其中，φ_ξ服从[0,2π)上的均匀分布，每一个φ_ξ对应一个傅立叶变换替代数据；

(3)计算与每一个替代数据对应的时频表示，并分别从中提取出各次谐波分量，计算出各个替代数据的辨识统计量这里定义显著性水平指标为：

式中，为满足D_s>D₀的替代数据的个数；假设创建N_s个替代数据并且将显著性水平指标设置为p，即至少有N_s×p个替代数据满足D_s>D₀才认为该分量不是噪声，从而继续分解过程；这个检测使用三组不同值的参数(α_A,α_ν)，即分别计算出D(1,1)，D(0,1)和D(1,0)的值，只要其中至少有一个值不符合无效假设，则认为满足D_s>D₀；

(4)计算谐波之间的相关度的综合度量值

其中，

式中，w_A,w_φ,w_v代表的权值；在这里，默认使用ρ^(h)≡ρ^(h)(1,1,0)为幅度和相位一致性分配相等的权值，且对频率一致性不分配权值；

(5)为了减少对真实谐波分量的错误判断，定义综合度量值的阈值为：

(6)当一个谐波分量的综合度量值指标满足ρ^(h)≥ρ_min且显著性水平指标significance-level≥p＝95％时，则认为此谐波分量通过检验，为真实的谐波分量；如果不能通过抗噪性替代检验，则停止非线性模式分解；

步骤C-6：将所有的真实的谐波分量相加构成一个非线性模式分量c₁′(t)；

步骤C-7：从信号s′(t)中减去该非线性模式分量，重复步骤C-1到步骤C-6，得到所有的非线性模式分量c_i′(t)；

步骤D：对得到的非线性模式分量c_i′(t)积分获得b_i(t)；

步骤E：对b_i(t)重新用非线性模式分解方法进行非线性模式分解；

步骤E-1：按步骤C-1到步骤C-6的过程对每个b_i(t)进行非线性模式分解；

步骤E-2：提取出每个b_i(t)的第一个非线性模式分量作为原始信号s(t)的非线性模式分量c_i(t)；最后得到原信号的非线性模式分解结果：

步骤F：对得到的s(t)进行谱分析，提取谐波信号；

步骤F-1：对c_i(t)做希尔伯特变换产生正交分量

步骤F-2：构造复信号z_i(t)，其表达式为：

步骤F-3：将复信号z_i(t)转化为极坐标形式，求得瞬时包络a_i(t)和瞬时频率ω_i(t)；z_i(t)的极坐标形式为：

其中，瞬时包络a_i(t)和瞬时相位φ_i(t)表示如下：

其瞬时频率ω_i(t)为瞬时相位φ_i(t)对时间的导数，表示为：

步骤F-4：求所有非线性模式分量c_i(t)的希尔伯特谱H(ω,t)；

希尔伯特谱定义为：

步骤F-5：求希尔伯特边际谱，提取信号谐波分量的谱峰；

希尔伯特边际谱的表达式为：

其中，T表示采样时间，其计算公式为：T＝N/f_s。

图2到图6依次为当信噪比SNR＝25dB时，分别采用差分非线性模式分解方法、非线性模式分解方法、经验模态分解方法、集合经验模态分解方法和差分经验模态分解方法得到的希尔伯特边际谱。从图2中可以看出，差分非线性模式分解方法能够清楚地提取出原始信号的基波频率(30Hz)、以及三次和五次谐波频率(90Hz、150Hz)，并且每个频率分量的谱峰十分尖锐；而非线性模式分解方法(图3所示)仅提取出了基波分量(30Hz)和三次谐波分量(90Hz)，未能提取出五次谐波(150Hz)分量；从图4-图6中看出，经验模态分解方法、集合经验模态分解方法和差分经验模态分解方法受噪声和混沌的干扰影响严重，仅仅能提取出信号的基波分量(30Hz)，无法准确提取出信号的其它谐波分量，仅在三次谐波频率(90Hz)附近有一定程度的突起。

为定量地比较差分非线性模式分解方法与其它方法的性能，计算指标R_k，定义为其中，k为谐波的次数，A_k是对应谐波频率的幅度，是被分解分量的谱的平均值。每种方法在不同信噪比下的R_k值如表1所示。

表1每种方法在不同信噪比下的R_k值

表1中DNMD为差分非线性模式分解方法，NMD为非线性模式分解方法，EMD为经验模态分解方法、EEMD集合经验模态分解方法和DEMD为差分经验模态分解方法，SNR为信噪比。

从表1所列数据可以看出：当SNR＝20dB时，差分非线性模式分解方法能够提取出原始信号的基波频率(30Hz)、三次和五次谐波频率(90Hz，150Hz)，非线性模式分解方法仅提取出了基波分量(30Hz)和三次谐波分量(90Hz)，而经验模态分解方法、集合经验模态分解方法和差分经验模态分解方法只提取出了基波分量(30Hz)。虽然在SNR＝15dB时，差分非线性模式分解方法也失去了提取五次谐波(150Hz)的能力，但其R₁＝52.74dB和R₃＝28.30dB均大于非线性模式分解方法的R₁＝51.90dB和R₃＝26.28dB。

综上所述，以上实施例仅用以说明本发明的技术方案，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，其均应涵盖在本发明的权利要求范围当中。