一种基于非线性模式分解和自适应最优核的时频分析方法与流程

本发明涉及非平稳信号的时频分析领域，尤其是一种基于非线性模式分解和自适应最优核的时频分析方法。

背景技术：

一个好的时频分析方法在非平稳信号的估计分析中的作用不言而喻。信号本身具有很多属性，而对于信号估计，频域特性是很重要的属性。

目前有很多信号处理和估计的方法，最经典的当然是傅立叶变换，它可以反映信号的频率特性，傅立叶变换对于平稳信号拥有良好的频域分辨能力，变换后可以得到信号的频谱，但是傅立叶变换丢失了时间信息，即各个频谱出现的时间并不知道，所以傅立叶变换不适合处理非平稳信号。对于实际环境里的非平稳信号，常用的线性时频表示方法有短时傅里叶变换、小波变换、S变换等，但是它们的时频精度有待提高；而Wigner-Ville分布、Cohen类时频分布等双线性时频分布虽然具有较高的时频精度，但却不可避免地在计算过程中引入了交叉项。

经验模态分解(Empirical Mode Decomposition，EMD)可将复杂信号分解为一系列本征模函数(Intrinsic Mode Function，IMF)。尽管EMD在抑制多分量信号的交叉项上有着不俗的表现，但是其缺点也很明显，抗噪性能差，并且会因为不连续信号或者噪声的存在而带来模态混叠现象。为了克服EMD的缺点，出现了一种基于噪声辅助的分析方法，叫做集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)。虽然EEMD在抗噪性能上要比EMD提升了一些，但是仍然达不到期望值。

时频表示(Time-Frequency Representation，TFR)是一种十分有效的瞬时频率估计方法，尤其是对于多分量信号。一般的TFR方法只有固定的窗函数或核函数，因此它们只对某一类信号的处理效果比较好。自适应最优核时频表示(AOK TFR)是一种基于信号的瞬时频率估计方法，采用可以随着信号变化而变化的径向高斯核函数，所以，AOK具有较好时频聚焦性能和抑制交叉项能力。

EMD+AOK可以在分析多分量信号的时候减小交叉项的影响，但是对噪声过于敏感。EEMD+AOK可以一定程度上提高抗噪性能，但是还远远不够。

非线性模式分解(Nonlinear Mode Decomposition，NMD)可以将信号分解为一系列的有物理意义的非线性模式分量，具有很强的噪声鲁棒性。它是将时频分析、数据替代检验以及谐波鉴别等多种方法融合起来的一种新型算法。然而，通过NMD分解得到的非线性模式分量并不是像IMF一样的单分量，需用AOK方法的核函数能自适应随信号的变化而变化的特点来进行优化。

基于NMD和AOK的新的时频表示方法，不仅利用了NMD对多分量信号有效的分解性能，同时还继承了AOK的优秀时频聚焦性能和有效抑制交叉项的能力，解决了EMD+AOK和EEMD+AOK的抗噪性能差以及仍含有一定量的交叉项的问题。

技术实现要素：

发明目的：为解决上述技术问题，本发明提供一种基于非线性模式分解和自适应最优核的时频分析方法，采用非线性模式分解(Nonlinear Mode Decomposition，NMD)把待处理信号分解为一系列非线性模式分量，抑制噪声的干扰；然后将得到的各信号分量经过自适应最优核(AOK)分析方法处理抑制交叉项，由于自适应最优核的分析方法(AOK)能自动跟踪分析信号的变化，核函数能根据信号的变化而自适应变化，对不同的信号产生的核函数总是最优的，所以其时频分析性能比较好。

技术方案：为实现上述技术效果，本发明提供的技术方案为：

一种基于非线性模式分解和自适应最优核的时频分析方法包括步骤：

(1)定义待处理信号为s(t)，s(t)的采样频率为f_s、数据长度为N；采用非线性模式分解法将待处理信号s(t)分解为一组非线性模式分量，即：

其中，c_i(t)为s(t)的第i个非线性模式分量，n(t)表示噪声；

(2)对步骤(1)分解出的各非线性模式分量采用自适应最优核分析方法进行分析，包括步骤：

(2-1)构建径向高斯核函数

式中，表示用于控制径向高斯函数在径向角方向的扩展函数；

(2-2)以得到随着信号自适应变化的最优核为目标问题，构建优化问题模型为：

设置优化问题模型的约束条件为：

其中，为第i个非线性模式分量c_i(t)在极坐标下的模糊度函数，β为最优核的体积；在直角坐标系中的表达式为：

其中，A_i(t，θ，τ)为在直角坐标系中的映射量，s^*(t)，w^*(u)是s(t)，c_i(t)，w(u)的共轭复数形式；窗函数w(u)为以t的中心，宽度为2T的对称菱形窗函数，且|u|＞T时，w(u)＝0，变量τ和θ是一般模糊域{τ，u}的参数，|τ|＜2T；

(2-3)根据A_i(t，θ，τ)求解优化问题模型，得到非线性模式分量c_i(t)的最优核函数Φ_(i)opt(t，θ，τ)；

(2-4)计算非线性模式分量c_i(t)的自适应最优核的时频表示为：

式中，为非线性模式分量c_i(t)在t时刻的能量值；

(3)根据步骤(2)得到的所有非线性模式分量的自适应最优核的时频表示，计算时频分析的结果为：

进一步的，所述步骤(1)中采用非线性模式分解法将待处理信号s(t)分解为一组非线性模式分量的方法包括步骤：

(1-1)计算信号s(t)的小波变换表达式W_s(ω，t)：

其中，是s(t)的傅立叶变换，即s⁺(t)是s(t)信号的正频率部分，是小波函数，与互为傅立叶变换对，且满足条件ψ^*(t)，分别是ψ(t)，的共轭复数；ω_ψ表示小波峰值频率，ω_ψ＝1，f₀为分辨率参数，用来权衡在变换过程中时间和频率的分辨率；

(1-2)判断步骤(1-1)计算得到的小波变换是否为s(t)的最佳时频表示；若判断结果为否，则计算信号s(t)的加窗傅里叶变换表达式G_s(ω，t)：

其中，g(t)是加窗傅里叶变换的窗函数，为g(t)的傅里叶变换，满足以下条件：

且

(1-3)找出信号s(t)最佳时频表示的所有的脊曲线，并用脊方法重构谐波分量，其中，第h次谐波分量为：

x^(h)(t)＝A^(h)(t)cosφ^(h)(t)，h∈[1，2，…，N]

式中，A^(h)(t)、φ^(h)(t)分别是第h次谐波分量的幅度和相位；v^(h)(t)是第h次谐波分量的频率，y^(h)(t)≡φ′^(h)(t)；φ′^(h)(t)为φ^(h)(t)对时间t的导数；N表示谐波分量的最高次数，即数据长度；

(1-4)利用抗噪性替代检验方法对提取出的谐波分量的真伪进行鉴别，筛选出所有的真实的谐波分量，

(1-5)将所有的真实的谐波分量相加，得到一个非线性模式分量c₁(t)；

(1-6)从原信号s(t)中减掉非线性模式分解得到的c₁(t)，并且对残余分量重复步骤1-1到步骤1-6，得到所有的各个非线性模式分量c_i(t)；

最终，原目标信号s(t)可以表示为：

式中，n(t)表示噪声。

进一步的，所述步骤(1-3)中找出信号s(t)最佳时频表示的所有的脊曲线的方法为：

在时刻t_i，找出h个极大值点，并将找出的h个极大值点相连，形成时刻t_i的脊曲线为：

上式中i＝1，2，…，N，N为数据长度；H_s(ω，t)是信号s(t)的最佳时频表示，即为W_s(ω，t)或G_s(ω，t)；H_s(ω，t)中找出的所有时刻的脊点连线，即构成N条脊曲线。

进一步的，所述步骤(1-3)中用脊方法重构h次谐波分量x^(h)(t)＝A^(h)(t)cosφ^(h)(t)包括以下步骤：

若待处理的信号s(t)最佳时频表示是小波变换，则

若待处理的信号s(t)最佳时频表示是加窗傅立叶变换，则

式中，和是通过抛物线插值而得到的改进的离散化影响因子。

进一步的，所述步骤(1-4)中筛选真实的谐波分量的方法包括步骤：

(1-4-1)构建辨识统计量D(α_A，α_ν)，用于衡量各个谐波分量的幅值和频率的有序度；D(α_A，α_v)的表达式为：

式中，和分别表示A^(h)(t)和v^(h)(t)的谱熵；α_A，α_v为权值系数；

(1-4-2)为原信号s(t)创建N_s个傅立叶变换替代数据；其中，第j个傅立叶变换替代数据的表达式为：

(1-4-3)计算每一个傅立叶变换替代数据的最佳时频表示公式，并分别从其中提取出各次谐波分量，计算出各个傅立叶变换替代数据对应的辨识统计量为：

定义显著性水平指标式中为满足D_s＞D₀的替代数据的个数；D₀为原信号的h次谐波的有序度；

设置显著水平指标为p，当存在x个替代数据满足D_s＞D₀且x≥p×N_s时，判定对应的谐波分量不是噪声；否则，判定对应的谐波分量是噪声；

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

其中，

式中，w_A，w_φ，w_v代表的权值，ρ^(h)为幅度和相位一致性分配相等的权值ρ^(h)≡ρ^(h)(1，1，0)；

(1-4-5)定义综合度量值的阈值

(1-4-6)当一个谐波分量满足ρ^(h)≥ρ_min且显著性水平指标p≥95％时，判定此谐波分量通过检验，为真实的谐波分量。

进一步的，所述步骤(1-4-3)中判断一个替代数据是否满足D_s＞D₀的方法为：

选取任意三组不同值的参数(α_A，α_v)，分别计算三组参数对应的辨识统计量值，若存在任意一个辨识统计量值大于D₀，则判定该替代数据满足D_s＞D₀。

有益效果：与现有技术相比，本发明具有以下优势：

NMD方法不仅对于高斯白噪声具有良好的抑制效果，同时也能够有效抑制其他类型的噪声信号，具有很强的噪声鲁棒性和广泛的适应性；然而，通过NMD分解得到的非线性模式分量并不是像IMF一样的单分量，需用AOK方法的核函数能自适应随信号的变化而变化的特点来进行优化，AOK方法能够有效地、自适应地跟踪非平稳信号的变化。本发明基于NMD和AOK的新的时频表示方法，不仅利用了NMD对多分量信号有效的分解性能，同时还继承了AOK的优秀时频聚焦性能和有效抑制交叉项的能力，解决了EMD+AOK和EEMD+AOK的抗噪性能差以及仍含有一定量的交叉项的问题。

附图说明

图1为本发明的流程图。

具体实施方式

以下以一个多分量非平稳仿真信号为例，详细说明本发明的实施方式以及优越性。

假设s(t)为一个含有高斯白噪声的多分量信号：

s(t)＝s_d(t)+n(t)

s_d(t)＝cos(20πt)+sin(200πt)+sin(400πt)+sin(100π(t-0.5)²)

式中，n(t)表示高斯白噪声；s_d(t)为理想的多分量信号。设置采样频率为1kHz，采样时间1s，数据长度为1000。设置窗长度2T＝128，核函数体积限制为β＝5。

本发明的流程如图1所示，包括以下步骤：

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

步骤B：对信号s(t)进行NMD分析；

步骤B-1：计算信号s(t)的小波变换(Wavelet Transform，WT)W_s(ω，t)，

其中，是s(t)的傅立叶变换，即s⁺(t)是s(t)信号的正频率部分，ψ(t)是小波函数，与互为傅立叶变换对，且满足条件ψ^*(t)，分别是ψ(t)，的共轭复数。表示小波峰值频率。

ω_ψ＝1

其中，f₀是一个分辨率参数，用来权衡在变换过程中时间和频率的分辨率(通常情况下默认为1)。时间分辨率高则频率分辨率降低，反之亦然。

步骤B-2：检查步骤B-1得到的小波变换W_s(ω，t)是否为信号s(t)的最佳时频表示，若不是，则采用加窗傅里叶变换(Windowed Fourier Transform，WFT)G_s(ω，t)，WFT定义为：

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

步骤B-3找出信号s(t)最佳时频表示的所有的脊曲线这里是第h次谐波的脊曲线。所谓的脊曲线就是时频图上一些局部极大值点所连接而成的曲线。各个极大值点就叫脊点。

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

上式中i＝1，2，…，N，N为数据长度。H_s(ω，t)是信号s(t)的最佳时频表示，即为W_s(ω，t)或G_s(ω，t)。

将H_s(ω，t)中找出的所有时刻的脊点连线，即构成H条脊曲线。

步骤B-4：用脊方法重构第h次谐波分量为x^(h)(t)＝A^(h)(t)cosφ^(h)(t)，其中A^(h)(t)、φ^(h)(t)分别是第h次谐波分量的幅度和相位，ν^(h)(t)是第h次谐波分量的频率，ν^(h)(t)≡φ′^(h)(t)，φ′^(h)(t)为φ^(h)(t)对时间t的导数。

若待处理的信号s(t)的最佳时频表示为W_s(ω，t)，那么可以通过以下公式计算得到信号s(t)的第h次谐波分量x^(h)(t)。

若待处理信号s(t)的最佳时频表示为G_s(ω，t)，那么计算s(t)的h次的谐波相关分量x^(h)(t)为：

式中，和是通过抛物线插值而得到的改进的离散化影响因子。

步骤B-5：利用抗噪性替代检验方法对提取出的谐波分量的真伪进行鉴别，筛选出所有的真实的谐波分量，并且当连续三个谐波分量被判断为假时停止分解过程。包括以下步骤：

(1)从TFR中提取出谐波分量并计算出对应的辨识统计量D(α_A，α_v)；

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

其中，α_A，α_v是计算D(α_A，α_v)的权值系数。

(2)为原信号s(t)创建N_s个傅立叶变换替代数据；

令的幅值保持不变，相位变为均匀分布在[0，2π)上的N_s个随机相位这种随机分布的每个相位对应的傅里叶反变换即为s(t)的一个傅立叶变换替代数据s′_j(t)。

这里j＝1，2，…，N_s。

(3)计算与每一个替代数据相关的TFR，并分别从其中提取出各次谐波分量，从而计算出各个替代数据对应的辨识统计量

定义显著性水平指标式中是D_s＞D₀的替代数据的个数；D₀为原信号的h次谐波的有序度。

假设创建N_s个替代数据并且显著水平指标设置为p，即至少有p×N_s个替代数据满足D_s＞D₀才能认为该分量不是噪声，从而继续分解过程。

判断一个替代数据是否满足D_s＞D₀的方法为：选取任意三组不同值的参数(α_A，α_v)，分别计算三组参数对应的辨识统计量值，若存在任意一个辨识统计量值大于D₀，则判定该替代数据满足D_s＞D₀。

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

其中，

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

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

(6)当一个谐波分量满足ρ^(h)≥ρ_min且显著性水平指标p≥95％时，则认为此谐波分量通过检验，为真实的谐波分量。

步骤B-6：将所有的真实的谐波分量相加从而得到一个非线性模式分量c₁(t)。

步骤B-7：从原信号s(t)中减掉非线性模式分解得到的c₁(t)，并且对残余分量重复步骤B-1到步骤B-6得到所有的各个非线性模式分量c_i(t)。

最终，原目标信号s(t)可以表示为：

式中，n(t)表示噪声。

步骤C：对NMD方法分解出的各个非线性模式分量c_i(t)进行AOK分析，其过程如下：

步骤C-1：首先选择径向高斯核函数，其定义如下：

式中，表示用于控制径向高斯函数在径向角方向的扩展函数。

步骤C-2：为了得到随着信号自适应变化的最优核，应求解以下有约束的优化问题。

约束条件为：

其中，是极坐标下的模糊度函数，为每个非线性模式分量对应的β是最优核的体积。

在直角坐标系中的定义A_i(t；θ，τ)为：

其中，s^*(t)，w^*(u)是s(t)，c_i(t)，w(u)的共轭复数形式；窗函数w(u)为以t的中心，宽度为2T的对称菱形窗函数，且|u|＞T时，w(u)＝0；变量τ和θ是一般模糊域{τ，u}的参数，且|τ|＜2T。

步骤C-3：把代入步骤C-2中的(4)，(5)式，可以通过求解这个有约束的优化问题得到一个最优核函数它与短时模糊函数一样会随着时间变化。

步骤C-4：计算当前时间点(t时刻)某个非线性模式分量c_i(t)的自适应最优核的时频表示AOK TFR_i：

式中，为信号某一分量c_i(t)在某一个t时刻的能量值。

步骤C-5：重复步骤C-1至步骤C-4，得到所有分解出的非线性模式分量的自适应最优核时频表示。

步骤D：将得到的所有信号分量的AOK TFR_i(即)求和，得到最终的时频分析的结果P_NMD-AOK(t，ω)：

由上式即可以最终得到NMD+AOK时频分析的结果。

将本发明所提供的的技术方案与现有技术相比，可以得到以下分析结果：

只用AOK方法会产生一定量的交叉项，并且抗噪性能很差；EMD+AOK和EEMD+AOK方法能在一定程度上抑制交叉项，但是仍然对噪声过于敏感；而本发明所采用的NMD+AOK方法不论是在去除噪声还是抑制交叉项方面，明显优于其他的几种方法。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。