★ FEATURED ARTICLE 简介这是一份面向信号处理研究者和工程技术人员的变分模态分解VMDMATLAB实现专门用于实测离散信号的模态分解可自适应处理非线性、非平稳信号与传统傅里叶、小波方法相比更能捕捉瞬态局部特征。压缩包仅含1个文件为约2KB的MATLAB脚本.m代码体积小、结构清晰便于直接阅读和二次修改由于只有单一脚本也能降低初学者定位核心算法的门槛。已有229人浏览学习说明该方法在信号处理实践中较受关注。读者拿到脚本后可导入自己的时间序列信号调试中心频率、迭代次数等关键参数得到若干具有不同频率特性的模态分量适用于噪声抑制、故障诊断、特征提取和信号恢复等场景脚本配合基本MATLAB操作即可运行且作者标明实测有效适合想快速上手VMD的入门者与实际数据分析人员。1. 拿到一段实测信号为什么我第一个想到的是VMD做故障诊断或者设备状态监测的工程师大概都经历过这种场景传感器采回来一段振动信号频谱图上有一片连续的宽峰不知道是轴承坏了、齿轮打了还是只是噪声用EMD分解各阶IMF在频谱上叠在一起模态混叠到怀疑人生滤波吧又怕把真正有用的特征滤掉。变分模态分解VMDVariational Mode Decomposition就是为了解决这类“信号成分分不开”的问题出现的。它能把一个非平稳、非线性的实测信号自适应地拆成若干个离散模态分量每个分量围绕一个中心频率、带宽受限对强噪声和间歇性干扰的容忍度比EMD高得多。适合做振动信号处理、故障特征提取、生物医学信号分析和地球物理数据处理的工程师与研究者。这篇笔记我会把VMD从原理、实现到实测信号落地的完整路径和踩过的坑一次讲清楚。2. 变分模态分解到底在做什么先理解“变分”和“模态”再动手2.1 VMD 在解决什么问题从 EMD 的模态混叠到自适应的变分框架在接触VMD之前多数人用的是经验模态分解EMD。EMD的思路很直接——把信号极值点包络的平均值当成“瞬时均值”不断用原始信号减掉这个均值筛出本征模态函数IMF。这个算法对简单信号很漂亮但实测信号一上来就露馅如果两个频率成分靠得近或者信号里有间歇性干扰筛出来的IMF在频谱上是叠在一起的这叫模态混叠另外EMD对采样率、噪声和端点极值非常敏感同样的信号换个时间窗分解结果可能都不一样。VMD换了个完全不同的思路。它把分解问题直接定义成一个约束优化问题在“所有模态加起来等于原始信号”的约束下让每个模态的估计带宽之和最小。这不是一个贪心的筛选过程而是一个全局寻优的过程因此模态与模态之间天然倾向于“各自占一段频带”。这种数学框架带来了两个直接好处第一模态之间的频谱重叠大幅减少混叠问题被约束条件压住第二对噪声有天然的鲁棒性——因为白噪声没有窄带结构会被分到一个宽频带的残余分量里而不会像EMD那样撕裂到各个IMF里。2.2 VMD 的数学模型最小带宽之和的约束优化VMD把信号 $f(t)$ 分解成 $K$ 个离散模态 $u_k(t)$每个模态都有一个中心频率 $\omega_k$。它的目标函数写出来是这个形式对每个模态做希尔伯特变换得到解析信号乘上 $e^{-j\omega_k t}$ 把频谱搬移到基带然后计算搬移后信号的梯度范数平方。这个值度量的是“这个模态绕着自己的中心频率有多分散”分散度越大带宽越大。VMD要做的就是让这 $K$ 个模态的带宽之和最小同时满足重构约束所有模态加起来等于原信号。这里的 $K$ 和惩罚参数 $\alpha$ 是人工设定的核心参数后面会专门讲怎么调。这个公式的工程含义很直接如果信号里有一个频率成分在150Hz附近VMD会倾向于生成一个中心频率接近150Hz、带宽尽可能窄的模态把它和别的成分分开。2.3 求解流程交替方向乘子法的一次次迭代这个约束优化问题没有闭合解常见做法是用交替方向乘子法ADMM迭代求解。迭代过程的核心是交替更新模态、中心频率和拉格朗日乘子。更新模态时把其他模态和乘子视作常量对 $u_k$ 求目标函数的极值。这个子问题在频域有解析解相当于对信号的维纳滤波——以当前中心频率为中心在频域做一个自适应带宽的滤波把属于这个模态的频带“挖”出来同时把不属于它的频带压下去。更新中心频率时把更新后的模态的功率谱质心当作新的中心频率也就是说中心频率会自动跑到模态能量最集中的地方。拉格朗日乘子则保证重构约束一步步收紧——迭代到后段重构残差会收敛到噪声水平。这个求解过程不是黑匣子。你可以看到每次迭代的中心频率 $\omega_k$ 变化轨迹如果两条模态的中心频率在迭代过程中互相追尾说明初始化或 $K$ 值有问题如果迭代几十次中心频率还在漂移说明惩罚参数太小、模态带宽约束不够强。把迭代信息打印出来看是排查VMD问题最直接的入口。3. 把 VMD 跑起来从测试信号到实测信号的最小实现3.1 先构造一个已知答案的测试信号最小可跑通的 MATLAB 脚本学习任何分解算法第一步永远是构造一个“你要什么答案它就要给你什么答案”的测试信号否则没法判断代码有没有跑对。下面这个脚本构造了一个三段已知频率合成信号加白噪声然后调用VMD函数分解%% 构造已知成分的测试信号 fs 1000; % 采样率 1kHz N 1024; % 数据长度 1024 点 t (0:N-1)/fs; s1 1.0 * cos(2*pi*60*t); % 60Hz 正弦分量 s2 0.6 * sin(2*pi*135*t); % 135Hz 正弦分量 s3 0.4 * cos(2*pi*260*t); % 260Hz 正弦分量 sig s1 s2 s3 0.15*randn(size(t)); % 叠加白噪声 %% 调用 VMD 分解 K 3; % 已知三个成分所以 K3 alpha 2000; % 带宽惩罚参数先取常用默认值 tau 0; % 噪声容忍系数0 表示严格重构 DC 0; % 第一个模态不设为直流分量 init 1; % 中心频率初始化为均匀分布 tol 1e-7; % 迭代收敛容差 [u, u_hat, omega] VMD(sig, N, alpha, K, tau, DC, init, tol); %% 查看分解结果和收敛中心频率 disp(最终中心频率 (Hz):); disp(omega(end, :) * fs / 2 / pi);这块代码逻辑不用一句一句背但有三个参数值得解释清楚。tau是噪声容忍系数设成0表示要求模态严格加起来等于原信号实测信号噪声大时可以放成0.1到0.5给重构留一点误差空间反而能让模态更干净。init是中心频率的初始化方式1表示在频带上均匀分布0表示全部从0开始统一优化实测信号建议用1否则低频分量可能会在迭代初期抢走过多能量。omega输出的是迭代过程中每一步的中心频率取最后一行就是收敛值。跑通这一段之后你就可以确信代码链路是通的。3.2 判断分解是否成功三个正经的检查点分解完不画图就等于没做。拿到u模态时域波形、u_hat模态频谱、omega中心频率迭代轨迹之后我一般按顺序检查三件事。第一看中心频率是否收敛。打印出来的omega(end,:)应该分别接近60、135、260单位换算成Hz后并且与真实成分的偏差在几个Hz以内。第二画各模态的时域波形和频谱图确认每个模态在时域上是一条干净的窄带振荡曲线频谱上只有一个孤立峰没有被相邻模态“咬掉”一部分。第三检查重构误差把sum(u,1)和原始信号相减残差的均方根应该接近你加的噪声幅度。如果残差明显高于噪声水平说明分解过程不收敛需要调大迭代次数或者检查输入参数类型。如果检查发现中心频率对但幅值偏小通常是初始化方式导致的能量分配不均如果中心频率对不上多半是K值给错了或者信号里有较强的高次谐波干扰。别急着调参数先把这三个检查在测试信号上跑通再去碰实测数据。这一步答应自己一定要做后面省掉的时间会值回票价。3.3 从测试信号到实测信号数据预处理这一步不能省把上面脚本里的sig换成实测信号之前有一项预处理不做的话VMD大概率会给你一堆“玄学”结果。实测信号通常有这几类问题一是直流偏置传感器或者采集卡会引入直流分量VMD会专门分出一个中心频率接近0的模态来装它白白浪费一个模态名额二是工频干扰50Hz的电网串扰会成为一个明显的窄带分量三是传感器安装共振引入的高频衰减或共振峰。我一般的处理顺序是先去均值再用陷波器滤掉工频干扰如果明显的话最后做一次带通滤波把分析频带限制在感兴趣的范围内。这个预处理对VMD的帮助非常大——VMD是全局优化方法信号里任何一个强窄带分量都会被当成一个模态单独分出来你不提前滤掉它就得用更多的K去承担它。实测信号里有时候还会碰到幅值很大的低频摆动比如转速波动、温度漂移这种分量也要先高通滤波去掉否则它会压低其他分量在目标函数里的权重。%% 实测信号预处理最小流程 sig_raw load(your_signal.mat); % 读入实测信号 sig_raw sig_raw(:); % 保证是列向量 sig sig_raw - mean(sig_raw); % 去直流 % 50Hz 陷波fs5120 的示例 fn 50 / (fs/2); % 归一化频率 [b, a] butter(2, [0.98*fn, 1.02*fn], stop); sig filtfilt(b, a, sig); % 零相位陷波 sig detrend(sig, linear); % 去除线性趋势filtfilt是零相位滤波这一点对后续分解很重要——普通filter会引入相位偏移VMD对时域波形是敏感的相位歪了会导致模态波形畸变。实测信号长度最好在1000点以上少于500点的信号VMD的边界效应会展不开分解出的端部模态基本都是假的。如果原始信号很长可以先截一段长度合适、工况稳定的片段再做分解而不是全丢进去。4. 把 VMD 用到实测信号核心参数的选择策略和调参手段4.1 K值怎么定先扫频谱再定模态数别靠猜VMD最头疼的参数就是模态个数 $K$。给少了几个频率成分挤在一个模态里分解等于没分解给多了一个真实的窄带成分会被劈成两三个相邻模态能量被撕碎。K值的初选我一般按三步走。第一步对预处理后的信号做FFT数一数频谱上有几个明显凸起的“峰群”。这些峰群就是天然模态的候选。第二步把峰的数量作为K的初值跑一次VMD看中心频率收敛值。第三步也是最关键的检查中心频率是否有两个值靠得很近比如只差5%以内如果有说明K太大了把K减一再来一次。反过来如果分解后某个模态的频谱宽得像个土堆说明这个模态里其实塞了不止一个成分K要加一。omega(end,:)这个收敛后的中心频率序列是判断K是否合适的直接证据比看时域波形直观得多。%% 用频谱峰个数辅助定 K [pxx, f] pwelch(sig, hann(512), 256, 1024, fs); mag 10*log10(pxx); % 挑出明显高于噪底的局部极大点阈值设为噪底10dB noise_floor prctile(mag, 30); [pks, locs] findpeaks(mag, MinPeakHeight, noise_floor10, ... MinPeakDistance, max(2, fs/512/2)); K_init length(pks); fprintf(建议初始 K %d\n, K_init);这里没有用复杂算法只是把“数峰”这一步自动化。注意MinPeakDistance至少设到两个频率分辨间隔以上避免同一个宽峰被数成两个模态。实测信号频谱往往是几个宽峰加许多小毛刺宁可先用少数明显主峰定K再用中心频率检查去迭代修正也不要第一次就把K设大。4.2 alpha 的脾气带宽惩罚参数决定模态“胖瘦”alpha是VMD里直接影响分解形态的参数。它控制约束条件里“带宽之和最小”这个目标的权重alpha越大每个模态越“瘦”带宽越窄中心频率分得越开alpha越小模态越“胖”带宽越宽容忍频率漂移的能力越强但相邻模态的频谱容易重叠。这个参数的经验取值区间非常大常见做法是试alpha 100/1000/2000/5000逐个比较频谱分离度。对振动信号我用2000起步对生物医学信号这类瞬态成分多的信号1000左右表现比较好如果信号里存在频率随时间漂移的成分比如变频器驱动的电机alpha取小一点反而能用一个模态把一段漂移的窄带包住。判断标准还是看频谱如果相邻模态的频谱在-20dB处不交叉alpha就算合适如果交叉得一塌糊涂alpha加倍再来。tau要不要调取决于你接不接受重构误差。严格重构tau0时VMD会把噪声也尽力塞进模态里模态波形会有些毛糙tau取0.1到0.3时算法允许一部分信号内容不被任何模态解释这部分正好是噪声模态时域波形会更干净。实测信号我多数情况用tau0直接上因为后续还要用包络谱提取特征模态干净与否影响不大但如果你的目的是看模态时域波形本身tau取0.2能明显好看。4.3 一次完整的调参流程从 K3 到收敛给你一条我实际用的调参路径照着走能少走弯路。第一步按4.1的频谱峰方法得到K_init。第二步固定 alpha2000、tau0、init1跑一次 VMD打印收敛中心频率。第三步检查中心频率序列如果两两之间最小间隔小于其中较小值的8%说明K过大如果存在中心频率接近信号的已知主频但频谱明显不对称、像被切了一刀说明K不足。第四步根据判断增减K再跑同一K下再试alpha为500、1000、2000、5000四档观察中心频率是否稳定。中心频率对alpha不敏感但频谱分离度对alpha敏感所以先定K再定alpha顺序不要反。第五步若分解结果依然古怪回看预处理和信号长度别在参数池里死磕。这套流程在测试信号上用不了两分钟在实测信号上可能要来回几次。注意每轮只改一个参数把每次的中心频率和alpha记下来否则两个参数一起动你根本不知道是谁搞坏了结果。5. 实测信号分解的避坑记录这些坑我替你踩过了5.1 信号两端出现明显振荡伪影边界效应现象分解后的每个模态在信号开头和结尾几十个点处都有大幅振荡中间正常频谱上表现为整个频带抬底。原因VMD内部要做希尔伯特变换和有限长度求导信号两端的数据窗不足以支撑窄带估计边界处模态定义被拉宽。数据越短边界效应占比越高。解决预处理时给信号两端各延长一段做完VMD再裁掉延长部分。延长方式最简单的是镜像延拓把信号前100个点翻转接到头部后100个点翻转接到尾部。不要用零填充零填充本身就是一种跳变会制造新的边界伪影。我用MATLAB的wextend函数做镜像延拓非常方便。5.2 K 取太大中心频率两个重叠在一起现象跑出来的K个中心频率里面有两个几乎一样对应的模态时域波形高度相似频谱峰也几乎重合。原因K取大了。VMD会尽量满足“K个模态”的约束当真实成分只有K-1个时它会把多余的那个模态硬生生插进某个已有模态里。但为什么中心频率不完全一样因为算法会在两个相同的中心频率上做带宽分配你看到的频谱实际上是同一个峰被切成两半。解决把K减1重跑。如果减1后重构误差明显变大那就不是K的问题回到预处理检查是不是有没滤掉的强窄带分量。这个坑最容易在实测信号里遇到因为实测频谱的峰不像测试信号那么清晰。5.3 实测信号噪声太大模态波形“毛刺化”现象分解出的模态时域波形高频毛刺严重频谱上除主峰外还有一串衰减的高频尾巴。原因tau0时算法必须严格重构信号噪声无法被拒在模态外就被均匀摊到每个模态上。测试信号噪声小看不出来实测信噪比低时非常明显。解决把tau设为0.1或0.2重新分解。你会看到模态时域波形光滑很多代价是重构残差变大——这没关系你要的是模态特征不是重构精度。如果改了tau还不够检查alpha是不是太小了alpha2000起步噪声大时用5000。5.4 强干扰分量把关注频段压得变形现象信号里有一个幅值很大的低频分量比如转频10Hz处幅值是关注频段100倍分解后关注频段所在模态幅值严重失真甚至中心频率偏离实际值。原因VMD的优化目标是最小化总带宽之和大振幅分量会占据主导权重迭代过程中的能量分配优先保大分量小分量被“挤”到收敛结果外。解决预处理阶段就把强干扰分量滤掉再分解不要指望VMD把它分到一个模态后再减掉。VMD的分辨率有限强分量和弱分量幅值差超过两个数量级时弱分量基本保不住。实测信号里我就会在预处理阶段把转频带高通滤除再跑VMD。5.5 分解结果跑一次一个样现象同样的信号、同样的参数跑两次VMD得到的模态时域波形不一样尤其体现在IMF的相位上。原因init参数控制了中心频率初始化默认的随机初始化会导致每次起点不同迭代收敛到的局部最优也不同。这不是VMD本身的bug而是所有交替优化算法的固有特性。解决固定init1频带均匀初始化并设置rand_seed如果有的话。只要初始化策略和种子固定VMD的结果就是确定性的。另外把中心频率的迭代轨迹打出来确认每一步没有异常跳变。做科研写论文的朋友记得把init固定状态写进代码注释里否则复现自己的结果都困难。6. 从仿真走向实测一套我一直在用的VMD落地验证流程前面写了这么多原理和参数最终你还是要拿着它对付一台真实设备的信号。我最后送你一套自己一直在用的落地流程它不复杂但每一步都是在真实项目里付出过代价换来的。拿到一段实测信号不要急着分解。先把整段信号的时域图、频谱图、包络谱三张图画出来把明显的主频、边带、噪底水平都看清。这一步能回答“信号里到底有什么”这个最关键的问题也决定了后面K的初值。然后预处理去直流、陷波、带通、去趋势再镜像延拓接着跑VMD最后做模态的包络谱分析——VMD分离出的每个模态本身是窄带信号用希尔伯特包络谱可以看到调制边带这是轴承故障和齿轮啮合故障诊断最常用的手法。%% 模态验证包络谱提取故障特征频率 for k 1:size(u, 1) analytic hilbert(u(k, :)); % 希尔伯特变换求包络 env abs(analytic); env env - mean(env); % 去直流 [pxx_env, f_env] pwelch(env, hann(256), 128, 512, fs); [~, loc] findpeaks(10*log10(pxx_env), SortStr, descend); fprintf(模态%d 包络谱峰值频率: %.2f Hz\n, k, f_env(loc(1))); end我自己的习惯是每次分解前在代码开头写一行注释记录“信号来源、采样率、分析频段、K和alpha从哪来”因为在实测项目里你很可能两周后回来看这段代码完全想不起来当初为什么选K6而不是K5。这行注释不花三秒钟但能在你怀疑人生时成为可靠的后悔药。最后说一句个人教训VMD调参最怕的就是“参数试一圈哪个图好看用哪个”。参数选择必须有依据依据就是中心频率的收敛值和频谱分离度不是视觉效果。把依据写下来、把每次实验的中心频率记下来你才能真正摸清这段信号的脾气。希望这篇笔记能帮你少走点弯路让VMD从“玄学”变成你手里的常规工具。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
