CN107506330A

CN107506330A - 一种基于粒子群算法的变分模态分解算法参数优化方法

Info

Publication number: CN107506330A
Application number: CN201710693495.4A
Authority: CN
Inventors: 于雪莲; 曲学超; 申威; 李海翔; 周云
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2017-08-14
Filing date: 2017-08-14
Publication date: 2017-12-22

Abstract

发明公开了一种基于粒子群算法的变分模态分解算法参数优化方法,属于非平稳信号时频分析方法的范畴.该方法主要流程：首先定义能够反映非平稳信号特征模态函数(Intrinsic Mode Function，IMF)分解效果的特征量即IMF傅立叶谱的熵，以IMF的傅立叶谱熵函数作为粒子群优化算法的适应度函数；再对PSO算法进行初始化，然后开始搜索适应度函的数最优解；最后，把得到的最优解作为变分模态分解算法的参数。本发明提出的方法，即通过粒子群算法自动求解出变分模态分解算法K，α的最优参数组合，求出的参数K和α不仅最优，而且快速高效。

Description

一种基于粒子群算法的变分模态分解算法参数优化方法

技术领域

本发明属于非平稳信号时频分析方法的范畴，应用于对非平稳信号进行变分模态分解时, 确定其所对应的变分模态分解算法的参数K和α的最优组合。其中，K代表将非平稳信号分解为单分量IMF(特征模态函数)的个数，α为二次惩罚参数。

背景技术

一直以来，经验模态分解(Empirical Mode Decomposition，EMD)算法广泛用于分解非平稳信号来获取其特征模态函数(Intrinsic Mode Function,IMF)，分解得到的IMF具有明确的物理意义，常用于非平稳信号的时频分。由于EMD算法采用循环筛分剥离的处理方式获得IMF 分量，所以EMD算法存在如下缺陷:缺乏有力的数学理论支撑；算法采用递归筛选逐次分解，无法反向误差校正；通常存在模态混叠现象或模态丢失等，这些缺陷严重的制约了EMD的应用。针对EMD存在的上述问题，Konstantin Dragomiretskiy和MominiqueZosso在2014年提出变分模态分解(Variational Mode Decomposition,VMD)算法。该方法在获取非平稳信号的IMF 分量时摆脱了EMD算法处理信号的模式，提出了将信号分解过程转移到变分的框架内的方法，通过引入增广Lagrange函数将有约束转化为无约束，用ADMM算法求最优解来实现信号自适应分解，在这个迭代求解的过程中，每个IMF分量的频率中心及带宽不断更新，最终实现信号频带的自适应剖分及得到若干窄带IMF分量。该算法成功的克服了EMD算法的诸多不足。因此变分模态分解算法广泛应用于非平稳信号的时频分析、雷达目标识瞬时频率特征的提取、传动系统故障的信息的诊断与分析、数据挖掘等领域。但是应用变分模态分解算法时针对特定非平稳信号进行特征模态函数提取时，需要事先指定相应的变分模态分解参数，不同的非平稳信号其对应的变分模态分解参数往往也不一样。因此若想正确提取非平稳信号的特征模态函数，必须保证该信号对应的变分模态分解参数设置的合理。比如，非平稳信号含有的IMF个数K，二次惩罚参数α，对偶上升步长μ，判别精度ε等。特别是参数K，α对IMF的提取结果影响很大，指定不合理的话，会发生模态混叠或模态丢失等问题。而目前确定待处理的非平稳定信号对应的模态个数K，惩罚参数α的常用方法是随机确定一系列K，α组合，然后对比分析各个组合参数下的分解的结果，进而确定比较合理的组合参数。这种方法对模态数目特别少时还可以应用，但是当模态较多时将会非常耗时，而且很难找到最优K，α组合。

对偶上升步长μ一般选取0到1小数即可，对分解结果影响不大，判别精度ε通常取1x10^-6即可。参数K及α对分解结果影响比较大而且很难确定。鉴于此，本发明通过构建IMF的傅立叶谱熵函数，作为粒子群算法(Particle Swarm Optimization，PSO)算法的适应度函数，通过粒子群算法自动求解出K，α的最优参数组合的方法。不仅求解出了最优参数K，α组合，而且快速高效。

发明内容

本发明的目的在于提出一种通过粒子群优化算法快速、准确的寻找出不同非平稳信号通过变分模态分解算法分解出特征模态函数时其所对应的变分模态分解算法的参数K和α最优组合的方法。

为了解决上述算法问题，本发明公开了一种基于粒子群算法的变分模态分解算法参数优化方法，该方法包括：

步骤1:对多分量非平稳信号y进行傅立叶变换得到多分量非平稳信号y的频谱，统计出频谱波峰的个数并记为M，认为信号y含有的单分量IMF的个数为M；

步骤2:初始化粒子群优化算法的参数：

粒子群优化算法：

其中，x＝[K,α]，K代表将非平稳信号分解为单分量IMF(特征模态函数)的个数，α为二次惩罚参数，影响分解出来的单分量的带宽，表示第i个粒子在第n次迭代时的速度，表示第i个粒子在第n次迭代时的位置，n为迭代次数，下标i代表第i个粒子，v代表粒子速度， x代表粒子位置，b_i代表第i个粒子搜索到的局部最优解，g_best代表整个粒子群迄今为止搜索到的最优位置，w表示惯性权重，e₁、e₂为学习因子，ξ为约束因子，λ每次迭代时为随机整数；

位置x初始化：x_i＝round(10rand(1,1)(x_max-x_min))+x_min (c1)

速度v初始化：v_i＝-1+round(2rand(1,1)) (c2)

round实现对某个数四舍五入取整，rand是在0到1之间随机取数，x_max、x_min分别代表自变量x＝[K,α]的最大值和最小值；

设置粒子群算法的粒子个数为N_p，粒子群搜索次数为L，其中(M-M/2)≤K≤(M+M/2)， 10≤α≤6000；

步骤3:采用步骤2每次迭代完成后得到的x＝[K,α]，作为变分模态分解算法的参数。然后使用变分模态分解算法对多分量非平稳信号y进行分解得到多个单分量IMF；

步骤4:计算步骤3得到的所有IMF模态混叠度系数C₁:

C₁＝Imf_w_sum-K+1

其中，Imf_w_sum表示所有单分量IMF中频率值的个数；

步骤5:计算步骤3得到的所有IMF模态重复度系数C₂

C₂＝Imf_rw_sum+K

其中：Imf_rw_sum表示所有IMF中不同频率值重复次数的和；

步骤6:构建粒子群优化算法的适应度函数：

其中：E(p)表示粒子群优化算法的适应度函数，C＝C₁·C₂，β(j)是第j个IMF的傅立叶谱；

步骤7:将步骤2中粒子群优化算法迭代L次，得到L个粒子群优化算法的适应度函数值，从这些值中寻找出最大值，认为该最大值对应的[K,α]为期望的最优解；

步骤8：将步骤7得到的最优解作为变分模态分解算法的参数，再对多分量非平稳信号y 进行变分模态分解得到多个IMF；判断分析IMF时域和频域波形，当不同IMF之间没有相同频率，而且同一IMF没有出现两个及以上的频率时，则认为步骤7最优解判定正确；否则，返回步骤1重新计算最优解。

进一步的，所述步骤2中惯性权重w设置为1，e₁、e₂均设置为2，约束因子ξ设置为1，随机系数λ每次迭代时取0到10之间的随机整数，粒子总数为[8，12]，迭代总次数为K和α最大值的乘积。其中，惯性权重w设置为1，保障了粒子群下一次迭代的速度合理地继承上一次迭代的速度；e1、e2均设置为2，权衡了粒子群算法中某个粒子的最优解和全局最优解对所有粒子速度和位置的影响程度，随机系数λ确保了粒子群求出的最优解为全局最优解，而不是局部最优解。

进一步的，所述步骤4中Imf_w_sum的计算方法为：通过傅立叶变换对步骤S3得到的多个IMF分别进行傅立叶变换，得到每个IMF的频谱并分别记为Imf_n_fft，简记为β；求取每个 IMF的傅立叶频谱β的均值FFT_mean和最大值FFT_max，以频谱均值FFT_mean的十倍Mean_ten作为阈值，若频谱不存在大于Mean_ten的值，则以频谱的最大值FFT_max的十分之一作为阈值；单独搜索出各IMF傅立叶频谱的所有波峰，将大于阈值的谱峰所对应的频率记为相应IMF的频率值Imf_n_fft_wi，并统计该IMF频率的个数Imf_n_w_num；再求取所有 IMF的频率值个数Imf_n_w_num之和Imf_w_sum。其中，阈值的确定方法，有效地避免了把 IMF的傅立叶频谱中由噪声引起的波峰所对应的频率记为该IMF的频率值的现象，进而确保了正确地统计出该IMF所包含的频率成分及个数。

进一步的，所述步骤5中Imf_rw_sum的计算方法为：选取不同IMF中的频率值进行比较，若选取的两个频率值差的绝对值小于等于0.0001则判定这两个频率值相等，统计所有IMF中某一频率值的个数P，则该频率值的重复次数为P-1，记为Imf_wi_num，再求取所有不同频率值重复次数Imf_wi_num之和Imf_rw_sum。其中，提出的判断两个频率是否相等的准则，在噪声干扰环境下，可以合理地判别某两个频率是否相同。

通过实施例的图2到图6可以证明本发明能够正确、高效地求出非平稳信号分解为多个单分量IMF所对应的变分模态分解算法的参数K和α。

附图说明

此处所说明的附图用来加深对本发明的进一步理解，构成本发明的一部分。

图1为本发明基于粒子群算法的变分模态分解算法的参数优化方法一个实施例的流程框图；

图2为仿真例sig信号的时域波形；

图3为仿真例IMF1的时频波形；

图4为仿真例IMF2的时频波形；

图5为仿真例IMF3的时频波形；

图6为仿真例IMF4的时频波形；

图7为仿真例各个模态的频谱信息度P_j。

具体实施方式

以下将配合实施例来详细说明本发明的实施方式，以便对如何应用本发明技术手段来解决技术问题有更加深刻的理解，以期达到良好地解决实际问题目的，并据以实施。实施例用到的非平稳信号是包含四个IMF的仿真多分量非平稳信号sig，具体参数现具体说明：

构建非平稳信号sig。

sig＝sig_1+sig_2+sig_3+sig_4+η

其中，sig_1＝sin(2π×30t)；sig_2＝sin(2π×50t)；sig_3＝sin(2π×60t)；sig_4＝sin(2π×130t)；η是均值为0，方差0.09的高斯白噪声，采样率800hz，图2给出了仿真信号sig的时域波形。

本发明实施例步骤按照本发明基于粒子群算法的变分模态分解算法的参数优化方法的步骤实施，即如图1所示步骤流程框图。具体实施如下：

步骤1:通过傅立叶变换得到仿真信号sig的频谱，并统计出频谱的波峰的个数M为5，预估计信号sig含有的单分量IMF的个数为5；

步骤2:初始化粒子群优化算法的参数：

粒子群算法中，惯性权重w设置为1，e₁、e₂均设置为2，约束因子ξ设置为1，λ每次迭代时取0到10之间的随机数即round(10·rand(1，1))。粒子群的适应度函数自变量分别为K和α取值范围分别为，2≤K≤8，10≤α≤6000。粒子的个N_p数设置为8，速度v范围设置为[-1,1]。粒子群搜索次数L设置为10000；

步骤3:使用变分模态分解算法对仿真信号sig进行分解得到其多个单分量IMF：

粒子群算法的位置参数K和α随机取值，作为变分模态分解算法的参数，并设定变分模态分解算法的对偶上升步长μ为0.2，判别精度ε取1x10^-6。然后，仿真信号sig进行分解得到其多个单分量IMF(特征模态函数)，并分别记为Imf_n；步骤S4:计算步骤S3得到的所有IMF模态混叠度系数C₁:

通过傅立叶变换对步骤S3得到的多个IMF分别进行傅立叶变换，得到每个IMF的频谱并分别记为Imf_n_fft，简记为β。求取每个IMF的傅立叶频谱β的均值FFT_mean和最大值FFT_max，以频谱均值FFT_mean的十倍Mean_ten作为阈值，若频谱不存在大于Mean_ten的值，则以频谱的最大值FFT_max的十分之一作为阈值。搜索出IMF傅立叶频谱的所有波峰，将大于阈值的谱峰所对应的频率记为相应IMF的频率值Imf_n_fft_wi，并统计该IMF频率的个数Imf_n_w_num。再求取所有IMF的频率个数Imf_n_w_num之和Imf_w_sum。最后，可求得所有IMF模态混叠度系数C₁＝Imf_w_sum-K+1。其中，K值为IMF的个数。

步骤4:计算步骤S3得到的所有IMF模态重复度系数C₂

当不同IMF的频率值Imf_n_fft_wi差的绝对值小于等于0.0001则判定这两个频率值相同，即这两个IMF的重复，并且将不同IMF相同的频率重复次数记为Imf_wi_num。再求取所有不同频率重复次数Imf_wi_num的和Imf_w_sum。最后，可求得所有IMF模态重复度系数C₂＝Imf_w_sum+K。其中，K值为IMF的个数。

步骤5:应用粒子群算法搜所到K和α作为变分模态分解算法的参数，对仿真信号sig进行分解，并求出粒子群的适应度函数AdaptFunc的函数值E_p：

步骤6:获取变分模态分解算法参数K和α的最优组合：

通过10000次粒子群搜索，得到10000组K和α组合，E_p取到最大值时对应的K和α分别为4和590，把该组合作为仿真信号y进行变分模态分解时的最优参数。

步骤7:把步骤S7获得的最优K和α的组合作为变分模态分解算法的参数，再仿真信号sig 进行变分模态分解得到多个IMF，分析分解结果。

在以步骤7得到的[K,α]＝[4,590]组合作为变分模态分解算法的参数时得到4个IMF，图3 为IMF1的时频波形，可知IMF1为正弦信号，频率为30hz，图4为IMF2的时频波形，可知 IMF2为正弦信号，频率为50hz，图5为IMF3的时频波形，可知IMF3为正弦信号，频率为60hz，图6为IMF4的时频波形，可知IMF4为正弦信号，频率为130hz，图7画出了步骤S6 定义的sig信号每个模态的频谱信息度p_j，各个p_j比例几乎相同，这与sig中各个平稳信号比重相同一致。综上，可知原始非平稳信号sig，在以PSO算法提供最优解作为变分模态分解的参数情况下，确实正确地分解出了其包含的所有平稳信号。IMF1对应sig_1，IMF2对应sig_2，IMF3对应sig_3，IMF4对应sig_4。得知[4，590]确实是[K，α]的最优组合。

上述说明示出并描述了发明应用的实施例，但如前所述，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述发明构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围，则都应在发明所附权利要求的保护范围内。

Claims

1.一种基于粒子群算法的变分模态分解算法参数优化方法，该方法包括：

步骤2:初始化粒子群优化算法的参数：

粒子群优化算法：

其中，x＝[K,α]，K代表将非平稳信号分解为单分量IMF的个数，α为影响分解单分量的带宽参数，表示第i个粒子在第n次迭代时的速度，表示第i个粒子在第n次迭代时的位置，n为迭代次数，下标i代表第i个粒子，v代表粒子速度，x代表粒子位置，b_i代表第i个粒子搜索到的局部最优解，g_best代表整个粒子群迄今为止搜索到的最优位置，w表示惯性权重，e₁、e₂为学习因子，ξ为约束因子，λ每次迭代时为随机整数；

位置x初始化：x_i＝round(10rand(1,1)(x_max-x_min))+x_min (c1)

速度v初始化：v_i＝-1+round(2rand(1,1)) (c2)

设置粒子群算法的粒子个数为N_p，粒子群搜索次数为L，其中(M-M/2)≤K≤(M+M/2)，10≤α≤6000；

步骤4:计算步骤3得到的所有IMF模态混叠度系数C₁:

C₁＝Imf_w_sum-K+1

其中，Imf_w_sum表示所有单分量IMF中频率值的个数；

步骤5:计算步骤3得到的所有IMF模态重复度系数C₂

C₂＝Imf_rw_sum+K

其中：Imf_rw_sum表示所有IMF中不同频率值重复次数的和；

步骤6:构建粒子群优化算法的适应度函数值：

其中：E(p)表示粒子群优化算法的适应度函数值，C＝C₁·C₂，β(j)是第j个IMF的傅立叶谱；

步骤7:将步骤2中粒子群优化算法迭代L次，得到L个粒子群优化算法的适应度函数值E_p，从这些E_p寻找出最大值，认为该最大值对应的[K,α]为最优解；

步骤8：将步骤7得到的最优解作为变分模态分解算法的参数，再对多分量非平稳信号y进行变分模态分解得到多个IMF；判断IMF时域和频域波形，当不同IMF之间没有相同频率，而且同一IMF内没有出现两个及以上的频率时，则认为步骤7最优解判定正确；否则，返回步骤1重新计算最优解。

2.如权利要求1所述的一种基于粒子群算法的变分模态分解算法参数优化方法，其特征在于所述步骤4中Imf_w_sum的计算方法为：通过傅立叶变换对步骤S3得到的多个IMF分别进行傅立叶变换，得到每个IMF的频谱并分别记为Imf_n_fft，简记为β；求取每个IMF的傅立叶频谱β的均值FFT_mean和最大值FFT_max，以频谱均值FFT_mean的十倍Mean_ten作为阈值，若频谱不存在大于Mean_ten的值，则以频谱的最大值FFT_max的十分之一作为阈值；单独搜索出各IMF傅立叶频谱的所有波峰，将大于阈值的谱峰所对应的频率记为相应IMF的频率值Imf_n_fft_wi，并统计该IMF频率的个数Imf_n_w_num；再求取所有IMF的频率值个数Imf_n_w_num之和Imf_w_sum。

3.如权利要求1所述的一种基于粒子群算法的变分模态分解算法参数优化方法，其特征在于所述步骤5中Imf_rw_sum的计算方法为：选取不同IMF中的频率值进行比较，若选取的两个频率值差的绝对值小于等于0.0001则判定这两个频率值相等，统计所有IMF中某一频率值的个数P，则该频率值的重复次数为P-1，记为Imf_wi_num，再求取所有不同频率值重复次数Imf_wi_num的和Imf_rw_sum。

4.如权利要求1所述的一种基于粒子群算法的变分模态分解算法参数优化方法，其特征在于所述步骤2中惯性权重w设置为1，e₁、e₂均设置为2，约束因子ξ设置为1，系数λ每次迭代时取0到10之间的随机整数，粒子总数为[8，12]，迭代总次数为K和α最大值的乘积。