类EMD的“信号分解方法”及MATLAB实现（第四篇）——VMD

重头戏来了。

在以往的应用经验里，VMD方法在众多模态分解方法中可以说是非常好的。从催更力度上看，这个方法也是格外受关注。笔者决定加快进度快一些写完这个方法，十月份了有些同学要开始做毕设，希望这篇文能帮上忙。

1. VMD（变分模态分解）的概念

VMD（Variational Mode Decomposition）即变分模态分解，与2014年由Dragomiretskiy[1]等人提出，虽然它也叫模态分解，但是和之前介绍过的EMD、EEMD、CEEMD、CEEMDAN这些模态分解方法在原理上有本质区别。VMD的整体思路为：

该方法假设任何的信号都是由一系列具有特定中心频率、有限带宽的子信号组成（即IMF）。
以经典维纳滤波为基础，通过对变分问题进行求解，得到 中心频率与带宽限制，找到 各中心频率在频域中对应的有效成分，得到模态函数。其模型构建涉及到维纳滤波、希尔伯特变换和解析信号等知识点。 [2]
VMD 的分解过程就是变分问题的求解过程，其算法主要包括 变分问题的构造和 变分问题的求解。 [3]VMD的求解过程主要包含两点约束：（1）要求每个模态分量中心频率的带宽之和最小；（2）所有的模态分量之和等于原始信号。 [4]

1.1 关于内涵模态分量

不同于黄锷先生提出的内涵模态分量（IMF）概念，VMD算法重新定义了约束条件更为严格的有限带宽的本证模态函数，该内涵模态分量被定义为调幅调频的分量模态函数，数学表达式为：

 为信号的包络幅值，为瞬时相位。需要注意该函数表征的分量也是同样满足EMD的约束条件的。

1.2 变分问题构造

由于变分问题涉及到泛函的相关知识，推导过程比较复杂，这里只介绍思路，对于详细过程有兴趣的同学可以直接看方法提出者的论文原文[1]。

所谓变分问题，就是求泛函的极值。在VMD中，泛函指的是VMD约束变分模型，而要求的极值，就是“每个模态分量中心频率的带宽之和最小”。

过去常遇到的是求函数极值，但有时我们需要对自变量也是函数的特殊函数求极值。

这种特殊函数即“函数的函数”，称为泛函，求泛函的极值问题称为变分问题。

构造出来的VMD约束变分模型是这样的[1]：

*上式表示的，就是求取“每个模态分量中心频率的带宽之和最小”时的模态函数和中心频率。 

1.3 变分问题求解

上式的求解过程应用到了比较多的数学工具，包括二次惩罚项、拉格朗日算子以及增广拉格朗日函数等等，详细过程[1]比较复杂，同样不展开说了。

2.VMD的显著特点

2.1 VMD思路

VMD方法是一种非递归变分模式的信号分解方法，其整体框架为一个变分问题。

2.2 VMD的独特优势

（1）可以指定想要得到的模态数。

（2）通过VMD方法分解出来的IMF都具有独立的中心频率，并且在频域上表现出稀疏性的特征，具备稀疏研究的特质。

（3）在对IMF求解过程中，通过镜像延拓的方式避免了类似EMD分解中出现的端点效应。

（4）有效避免模态混叠（K值选取合适的情况下）。

2.3 VMD的重要参数——模态数K

在使用EMD方法时，有很多同学就经常问要怎样指定分解出的IMF数。VMD方法就可以很好地满足这个要求，因为VMD的一个明显特点，就是信号在做VMD分解前，需要先设定模态分量的个数K。但是成也萧何败也萧何，某些应用场合下可指定IMF数属于优点，但是对于某些不预知信号隐含模态数的场景，怎样设置这个K值反而会让人挠头。

若设定的K小于待分解信号中有用成分的个数（欠分解），会造成分解不充分，导致模态混叠；若设定的K值大于待分解信号中有用成分的个数（过分解），就导致产生一些没有用的虚假分量。因此，K值的确定对于VMD就非常重要。

至于如何确定K值，可以从以下几个方面考虑：

（1）某些场景下可以预知模态数目，比如对于某齿轮箱振动数据，通过分析可知信号数据中主要包含齿轮的啮合频率和主轴转动频率，那么K可以设置为2。

（2）可以通过K值从小到大逐个尝试，并结合分解结果分析确定K值：随着K值增大各主要频率段数据能分布到不同IMF分量中，并且没有产生虚假分量，此时K值就是比较合适的。这个方法虽然比较高效且靠谱，但是不便于写到论文中，这就需要用到方法（3）了。

（3）可以使用一些辅助算法，比如使用峭度最大原理[5]、能量差值原则[6]等，或者结合一些寻优算法对K值（也可以同时对α）寻优。

2.4 VMD的重要参数——惩罚系数α

VMD 分解的过程中，预设的 K 值决定着 IMF 分量的个数，惩罚系数 α 决定着 IMF 分量的带宽。惩罚系数越小， 各 IMF 分量的带宽越大，过大的带宽会使得某些分量包含其他分量信号；α值越大，各IMF分量的带宽越小，过小的带宽是使得被分解的信号中某些信号丢失。该系数常见取值范围为1000~3000。

2.5 VMD的另一个参数——收敛容差tol

是优化的停止准则之一，即在连续两次迭代中，当向 IMF 收敛的绝对平均平方改进小于 tol时，优化停止。通常可以取 1e-6~5e-6。

3. VMD的应用案例

vmd函数在MATLAB2020a版本中内置到官方库了，不过按照“类EMD”系列的代码的统一风格，笔者重新进行了封装，一来使该方法可以应用于MATLAB2020以前版本，二来增加了绘制IMF分量与频谱对照的绘图函数，同时做了一些其他的适应性调整。

下边使用一个例子来展示VMD方法的突出效果：

首先生成一组测试信号，由二次趋势项、啁啾信号和分段余弦共同组成，采样频率为1000Hz。

%% 1.生产仿真信号
fs = 1e3;
t = 0:1/fs:1-1/fs;

x = 6*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...
    [cos(60*pi*(t(t<=0.5))) cos(100*pi*(t(t>0.5)-10*pi))];

figure('color','w')
tiledlayout('flow')
nexttile
plot(t,[zeros(1,length(t)/2) cos(100*pi*(t(length(t)/2+1:end))-10*pi)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,[cos(60*pi*(t(1:length(t)/2))) zeros(1,length(t)/2)])
xlabel('Time (s)')
ylabel('Cosine')

nexttile
plot(t,cos(4*pi*t+10*pi*t.^2))
xlabel('Time (s)')
ylabel('Chirp')

nexttile
plot(t,6*t.^2)
xlabel('Time (s)')
ylabel('Quadratic trend')

nexttile(5,[1 2])
plot(t,x)
xlabel('Time (s)')
ylabel('Signal')

四个信号分量及最终合成的Signal信号

这种程度的复合信号可以说是比较复杂的了，让我们看下对于该信号的VMD分解效果：

%% VMD分解参数设置
alpha=2000;  % alpha   - 惩罚因子
tol=1e-6;    % tol     - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6
K=4;         % K       - 指定分解模态数
[imf,CenFs] = pVMD(x,fs, alpha, K, tol); %调用函数实现分解与画图

Signal信号的VMD分解结果

上图是对于Signal信号的VMD分解结果，可以说是近乎完美了。

作为对比，EMD分解结果是下图这样的，就十分差强人意了。

Signal信号的EMD分解结果

在EEMD、CEEMD等改进方法中得到的结果也同样不理想，也正因此VMD方法可以说是笔者在解决工程问题时最常用的分解方法。

4.关于封装函数

上边案例中用到的函数“pVMD”就是笔者封装好的vmd画图程序，只需要输入待分解信号、采样频率、α值、K值、tol值，就可以得到信号的分解结果IMF和中心频率CenFs。这个函数的介绍如下：

[imf,CenFs]=pVMD(x,fs, alpha, K, tol);
% function [imf,CenFs] = pVMD(y,FsOrT, alpha, K, tol)
% 画信号VMD分解图
% 输入：
% data：待分解的数据(一维)
% FsOrT：采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% alpha   - 惩罚因子
% K       - 指定分解模态数
% tol     - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6

% 输出：
% imf：内涵模态分量，统一为n*m格式，其中n为模态数，m为数据点数。例如 imf(1,:)即IMF1，imf(end,:)即为残差
% 注意，为了与其他“类EMD”方法分解出来的imf分量保持一致，本程序内将imf排序进行了翻转
% 即保证imf排列是从高频向低频排列
% CenFs：即CentralFrequencies，各imf分量的中心频率
% 注意：在使用该代码之前，请务必安装工具箱：http://www.khscience.cn/docs/index.php/2020/04/09/1/

按照惯例，还有一个绘制VMD分解图及其频谱图的函数：

[imf,CenFs] = pVMDandFFT(x,fs, alpha, K, tol);
% function [imf,CenFs] = pVMDandFFT(y,FsOrT, alpha, K, tol)
% 画信号VMD分解与各IMF分量频谱对照图
% 输入：
% y：待分解的数据(一维)
% FsOrT： 采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% alpha   - 惩罚因子
% K       - 指定分解模态数
% tol     - 收敛容差，是优化的停止准则之一，可以取 1e-6~5e-6

% 输出：
% imf：内涵模态分量，统一为n*m格式，其中n为模态数，m为数据点数。例如 imf(1,:)即IMF1，imf(end,:)即为残差
% 注意，为了与其他“类EMD”方法分解出来的imf分量保持一致，本程序内将imf排序进行了翻转
% 即保证imf排列是从高频向低频排列
% CenFs：即CentralFrequencies，各imf分量的中心频率

% 例1：（FsOrT为采样频率）
% fs = 100;
% t = 1/fs:1/fs:1;
% y = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pVMDandFFT(y,fs,2000,2,1e-7);
% 例2：（FsOrT为时间向量，需要注意此时FsOrT的长度要与y相同）
% t = 0:0.01:1;
% y = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pVMDandFFT(y,t,2000,2,1e-7);

% 注意：在使用该代码之前，请务必安装工具箱：http://www.khscience.cn/docs/index.php/2020/04/09/1/

对于上边的案例，使用这个函数就可以得到各个imf分量的频谱图了：

Signal信号的VMD分解结果及其频谱图

上边的测试代码和封装函数（pVMD、pVMDandFFT等），包括工具箱都可以在下述链接中获取。

VMD画图工具（公开版） | 工具箱文档

EMD、EEMD、CEEMD以及HHT相关的程序也有，编程不易，感谢支持~

关于EMD、EEMD、CEEMD、CEEMDAN和HHT的相关介绍可以看这里：

Mr.看海：这篇文章能让你明白经验模态分解（EMD）——EMD在MATLAB中的实现方法

Mr.看海：希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位——Hilbert分析衍生方法及MATLAB实现

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第一篇）——EEMD

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第二篇）——CEEMD

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第三篇）——CEEMDAN

5. 结语

截至这篇，“类EMD”信号分解方法系列已经包含了EMD、EEMD、CEEMD、CEEMDAN和VMD这5篇，常用的方法基本已经涵盖，有部分比较小众的暂时不纳入，如果有哪种比较想看的可以在文章后留言。后续该系列将针对小波分解、小波包分解、SWT、EWT等方法介绍。

参考

1. ^abcd[1]Dominique, Zosso, Konstantin, et al. Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing A Publication of the IEEE Signal Processing Society, 2014. Variational Mode Decomposition | IEEE Journals & Magazine | IEEE Xplore
2. ^赵紫薇. 基于互信息VMD和SVM的煤矿流体管网泄漏检测算法与软件设计[D].山东科技大学,2020.
3. ^陶帅. 基于变分模态分解的语音增强算法的研究[D].哈尔滨理工大学,2021.
4. ^李士松. 基于改进VMD的单通道盲源分离滚动轴承复合故障诊断[D].哈尔滨理工大学,2021.
5. ^吴文轩, 王志坚, 张纪平,等. 基于峭度的VMD分解中k值的确定方法研究[J]. 机械传动, 2018, v.42；No.260(08):159-163.
6. ^宋玉琴,邓思成,路彦刚.K值优化的VMD在轴承故障诊断中的应用[J].测控技术,2019,38(04):117-121.