简介这是一份基于MATLAB的变分模态分解VMD算法程序包面向从事信号处理、故障诊断、数据分析的研究者与初学者。程序包含VMD核心函数、仿真信号脚本、包络谱与Hilbert分析工具以及多组轴承故障数据内圈、外圈等便于直接运行并观察非线性非平稳信号的分解与特征提取效果。压缩包共28个文件以m函数脚本、txt说明文档、mat数据文件和asv备份文件为主整体约9.62MB结构紧凑适合学习调试。目前已有1103人学习下载。通过调用VMD.m并调节模态数K与正则化因子α可灵活研究不同参数对分解结果的影响配套的仿真及实测数据脚本则展示了从信号加载、VMD分解到包络谱分析的一站式流程适合作为入门VMD算法与MATLAB实现的动手参考资料。1. 拆开 VMD 程序包之前先弄明白它和 EMD 差在哪处理滚动轴承振动信号、心电信号这类非线性非平稳数据时EMD经验模态分解是最先想到的工具但它有个老毛病模态混叠。一次冲击事件被拆到好几个 IMF 里包络谱上出现一堆假峰故障特征频率反而看不出来。VMD变分模态分解把「递归筛信号」改成「一次性求解约束优化问题」每个模态都被限制在各自的中心频率附近混叠问题被显著压制。这套 VMD.zip 里既有核心的 VMD.m也有 hua_fft1.m、hua_baoluo.m、hua_hilbert1.m 这些配套分析脚本还带着 1750内圈.mat、1750外圈.mat 等真实轴承故障数据适合做故障诊断、信号降噪、特征提取的工程师直接跑通一整套流程。接下来按「原理→文件→参数→实战→验证」的顺序逐一拆。2. VMD 的变分框架与程序文件逐项拆解2.1 变分问题如何构造带宽最小化与约束EMD 的思路是不断地找包络均值、剥出本征模态函数本质上是递归筛选误差会逐层累积。VMD 换了个思路先假设原始信号是由 K 个有限带宽的模态叠加而成每个模态都有一个中心频率然后构造一个约束优化问题——目标函数是全部模态的估计带宽之和最小约束条件是所有模态加起来等于原始信号。数学形式可以写成对每个模态 uk(t)先做 Hilbert 变换得到单边频谱再乘上指数项 e^{-jωk t} 把频谱搬到基带接着计算梯度范数的平方作为带宽估计。整个优化问题如下min Σₖ || ∂t [ (δ(t) j/πt) * uk(t) ] · e^{-jωk t} ||²约束 Σₖ uk(t) f(t)这里 ∂t 表示对时间求导* 是卷积ωk 是第 k 个模态的中心频率。直接求解带约束的变分问题不现实VMD 的做法是引入二次罚项 β 和拉格朗日乘子 λ(t)把约束罚进目标函数再通过交替方向乘子法迭代更新 uk、ωk 和 λ。每次迭代先固定其他变量更新第 k 个模态再固定模态更新对应中心频率最后更新乘子。迭代到满足收敛容差时停止输出的 uk 就是分解结果。理解了这一点你再看 VMD.m 里那些循环和矩阵运算就不会懵。程序里反复出现的u_hat、omega分别对应频域中的模态估计和中心频率迭代终止条件是前后两次更新的差值小于用户指定的 tol。这个框架天然支持对分解结果的数学解释而不是像 EMD 那样依赖经验筛选这也是它适合二次开发、能接进深度学习和故障诊断流程的原因。2.2 VMD.zip 内文件的角色划分拿到压缩包后不要急着运行先把文件清单过一遍搞清楚每个脚本在整条链路里的位置。下表是我拆包时整理的对应关系文件角色典型用途VMD.m核心分解函数输入信号和参数输出各模态及中心频率fangzhen.m / fangzhen.asv仿真入口脚本构造仿真信号、调用 VMD、画图仿真信号.m信号生成脚本定义 fs、时间轴和叠加分量hua_fft1.mFFT 频谱分析查看原始信号频域特征hua_hilbert1.mHilbert 变换求解析信号、瞬时幅值/相位hua_baoluo.m / hua_baoluo.asv包络谱分析解调出冲击成分特征频率Untitled.m测试/主控脚本串联数据读取、分解、绘图1750内圈.mat、1750外圈.mat故障数据轴承内圈/外圈故障振动信号72内.mat、1772内圈.mat附加数据其他工况或测点的样本6.txt、7.txt、8.txt、9.txt 等文本格式信号用于 load 直接读取的原始序列这些 .txt 文件通常是单列数值适合先用load(9.txt)读进来做快速验证.mat 文件则可能包含多个变量读入后先whos看一眼变量名和尺寸别默认变量名就叫 data。hua_baoluo.m和hua_hilbert1.m是配套的前者做包络谱时需要调用 Hilbert 变换得到包络信号后者把这一步封装好了单独拎出来也能用于瞬时特征提取。2.3 仿真信号入口脚本怎么读直接看fangzhen.m是最快的上手方式。常见做法是脚本先构造一个多分量信号比如低频正弦叠加高频调幅分量再加一点噪声然后调用 VMD.m。如果包里的仿真脚本没有按这个结构组织下面这个模板等价于它做的事fs 1000; N 1000; t (0:N-1)/fs; % 三个分量低频、中频调幅、高频正弦 x1 0.5 * cos(2*pi*20*t); x2 (1 0.3*cos(2*pi*5*t)) .* sin(2*pi*85*t); x3 0.3 * sin(2*pi*200*t); x x1 x2 x3 0.1*randn(size(t)); % 调用 VMDK3 对应三个分量 [u, ~, omega] VMD(x, 2000, 0, 3, 1e-7);这里 VMD 函数的形参按常见实现排列为信号、惩罚因子 alpha、噪声容忍度 tau、模态个数 K、收敛容差 tol。omega 返回的是归一化角频率要转成物理频率需要乘 fs/(2π)。如果某个文件的调用顺序不一样打开 VMD.m 看函数声明行就能确认不要凭记忆传参。跑通这个仿真后你就能验证 K 取 3 时三个模态是否对应 20Hz、85Hz、200Hz 三个分量这是检验程序包是否完整的第一步。3. 分解效果由参数决定alpha、K 的设置与副作用3.1 alpha 控制的是带宽不是分解个数alpha 是二次罚项的权重数值越大对模态带宽的限制越强每个模态被约束得越窄频域上各模态的重叠越少。这看起来是好事但 alpha 过大会让模态丢失细节把原本该包含的边频带切掉alpha 过小则模态带宽过宽两个模态会在频域上交叉出现类似 EMD 的混叠现象。对轴承故障信号故障冲击会产生一系列以共振频率为中心的边频带alpha 太大把这些边频削掉后包络谱里就只剩下孤零零的特征频率峰调制信息丢失。反过来alpha 太小则噪声会作为独立模态被分解出来。常见做法是先固定 alpha2000 跑一遍再看模态的频谱是否有明显重叠重叠明显就往大调边频缺失就往小调。实际项目中 alpha 落在 50010000 这个区间都是正常的不要迷信默认值。3.2 K 的选择方法中心频率观测与相关性校验K 是模态个数它的影响比 alpha 更直接。K 太小多个频率成分被塞进一个模态里模态频谱呈多峰形状K 太大会把单成分信号硬拆成两段出现频域相邻的虚假模态。判断 K 是否合适的核心手段是看中心频率的分离情况逐步增大 K把每次分解得到的中心频率打出来如果新增加的那个中心频率与已有某个中心频率非常接近说明 K 已经过大了。相关性校验也常用。分解出的每个模态与原始信号做 Pearson 相关相关系数高的模态代表真实成分相关系数突然跌到接近 0 的模态基本都是虚假分量。比如 K5 时第 5 个模态与原始信号的相关系数只有 0.02而前 4 个都大于 0.3那基本可以判定 K 应该取 4。注意相关性强弱和能量大小不是一回事一个高频小幅值成分相关系数可能不高但它确实是物理信号的一部分所以要结合频谱看不能只盯相关系数这一个指标。3.3 VMD.m 的调用方式与返回结果包里 VMD.m 的典型签名是[u, u_hat, omega] VMD(signal, alpha, tau, K, tol)其中 signal 是 1×N 的行向量输出 u 是 K×N 的矩阵每一行是一个模态u_hat 是对应模态的频域表示长度与信号相同omega 是 K 个中心频率的归一化角频率。tau 是噪声容忍度设为 0 表示完全强制重构精度设为非零值则允许一定残留噪声实测含噪信号时 tau 取 0.010.5 能避免把噪声强行掰成模态。% 对信号做预处理去均值 归一化 x x(:); x x - mean(x); x x / max(abs(x)); % 分解得到 4 个模态 [imf, ~, omega] VMD(x, 2000, 0.1, 4, 1e-7); % 换算成物理频率并打印 fs 1000; fprintf(中心频率(Hz): %.2f %.2f %.2f %.2f\n, ... omega*fs/(2*pi));预处理非常重要。去均值是防止直流分量占据一个模态归一化是避免幅值过大导致迭代不收敛。如果 VMD 分解后重构信号与原始信号误差很大优先检查输入是否做了这两步。另外 VMD 对数据长度敏感点太少解不出准确的带宽做短样本分析时至少保证 N 500否则结果不稳定。3.4 参数速查表参数常见范围经验初值表现出的问题与调整方向alpha500100002000模态频谱重叠→增大边频/细节丢失→减小K284模态多峰→增大出现相邻中心频率→减小tau00.50无噪声、0.1含噪分解结果忽大忽小→适当增大 tautol1e-91e-51e-7收敛太慢→放宽结果不收敛→调严这几个参数不是独立起作用的。alpha 和 K 共同决定模态的频域划分方式改 K 之后最好重新审视 alpha 是否还合适。调参的顺序我一般固定为先定 K靠中心频率观察再调 alpha靠频谱重叠和边频保留情况最后处理 tau 和 tol 这两个次要参数。4. 滚动轴承故障数据实战VMD 分解加包络谱识别4.1 读取 .mat 数据并做预处理包里的 1750内圈.mat、1750外圈.mat 是滚动轴承故障信号1750 通常对应 1750 r/min 的工况转速。先用 whos 查看变量再用plot看波形的冲击特征注意数据可能是列向量。% 载入内圈故障数据 load(1750内圈.mat); whos; % 确认变量名和尺寸 % 假设变量名为 data x data(:); fs 12000; % 西储大学轴承数据的常见采样率按实际填写 x x(1:fs*2); % 取前 2 秒数据避免数据量太大 % 去除直流和趋势项 x x - mean(x);如果 .mat 文件里装的是多个测点的信号选通道时优先选振动冲击最明显的那个。先做去均值和线性趋势去除原因很简单VMD 把趋势项当成低频模态输出浪费一个分解名额后续包络分析还要再剔除一次。数据量过大时 VMD 迭代速度明显变慢取 2 秒甚至 1 秒数据足够看到故障特征频率没必要把整段数据一次性丢进去。4.2 计算内圈与外圈的故障特征频率故障诊断的目标不是看波形有多乱而是找到冲击重复的频率。对内圈故障滚动体每经过一个缺陷点都会产生一次冲击冲击频率等于内圈故障特征频率 BPFI外圈故障对应 BPFO。计算公式如下BPFI n · fr · (1 d/D · cosφ) / 2BPFO n · fr · (1 − d/D · cosφ) / 2其中 n 是滚动体数量fr 是转频d 是滚动体直径D 是节圆直径φ 是接触角。这些几何参数属于具体轴承型号比如 6205 深沟球轴承的 n9接触角近似 0内圈特征频率系数约为 5.415外圈约为 3.585。nm 1750 r/min 时转频为 1750/60 ≈ 29.17 Hz代入系数fr 1750 / 60; % 转频 n 9; % 滚动体个数 BPFI n * fr * 5.415 / 9; % 换算到对应系数 BPFO n * fr * 3.585 / 9; fprintf(BPFI%.2f Hz, BPFO%.2f Hz\n, BPFI, BPFO);如果不知道轴承几何参数也可以用转速和故障特征系数近似估算分析时在包络谱上找 BPFI 或 BPFO 的整数倍频即可通常前 4 阶谐波都会出现。4.3 执行 VMD 分解并绘制模态对内圈故障信号K 取 4 比较合适一个低频旋转分量、一个共振带、一两个噪声分量。执行分解并观察各模态时域波形和频谱。K 4; alpha 1200; % 内圈故障边频带丰富alpha 偏小保留边频 [imf, ~, omega] VMD(x, alpha, 0.05, K, 1e-7); % 观察中心频率 disp(omega * fs / (2*pi)); % 画前 4 个模态的时域波形 for k 1:K subplot(K, 1, k); plot((0:length(imf(k,:))-1)/fs, imf(k,:)); ylabel([IMF, num2str(k)]); endalpha 取 1200 而不是默认 2000原因是内圈故障的调制边带比较宽太大的惩罚因子会把边频压掉。分解后重点看冲击最明显的模态——时域波形上呈现稀疏的大幅值脉冲这类模态往往对应故障冲击激起的共振成分。如果多个模态都有冲击用峭度来选峭度是四阶统计量冲击性强的信号峭度远大于 3。4.4 包络谱分析定位故障选定冲击明显的模态后做包络谱先求 Hilbert 变换得到解析信号取模得到包络对包络去均值后再做 FFT。hua_baoluo.m 封装的正是这个流程等价实现如下imf_sel imf(1, :); % 按峭度筛选出的模态 env abs(hilbert(imf_sel)); env env - mean(env); Nfft 2^nextpow2(length(env)); Y fft(env, Nfft); P abs(Y(1:Nfft/2)) * 2 / Nfft; f_axis (0:Nfft/2-1) * fs / Nfft; plot(f_axis, P); xlim([0 500]); xlabel(频率 (Hz)); ylabel(幅值);将峰值频率与 BPFI 对比如果峰值落在 BPFI 或其 2 倍、3 倍频附近并且谱峰明显高于周围噪声底就可以判定为内圈故障。外圈故障数据用同样的代码切换 BPFO 验证。注意包络谱前几个低频点经常有转频及其谐波这些不是故障特征区分方法是看它是否等于 fr 的整数倍而不是 BPFI/BPFO 的整数倍。VMD 分解在这个流程中的作用是把故障冲击模态从强噪声和旋转分量中分离出来让包络谱不被无关频率成分干扰。5. 用三个指标验证分解质量再谈参数自动寻优5.1 中心频率分离度判断 K 是否过大每次调完 K打印中心频率比较相邻中心频率的最小间隔。正常分解下各中心频率间距应明显大于 0.3 倍带宽如果最小间隔小到和带宽一个数量级说明两个模态混叠了。for K 2:8 [~, ~, omega] VMD(x, 2000, 0.05, K, 1e-7); cf sort(omega * fs / (2*pi)); min_gap min(diff(cf)); fprintf(K%d, 中心频率[%s], 最小间隔%.2f Hz\n, ... K, num2str(cf, %.1f ), min_gap); end观察哪个 K 对应的最小间隔突然缩小一个数量级K 的上限就在那个值之前。这是我在实际项目中用的第一道检查线比肉眼看频谱可靠。5.2 峭度与相关系数筛出有效模态分解完成后用峭度筛选冲击模态用相关系数剔除无效模态。峭度大于 3 表示信号比正态分布更尖锐故障冲击模态的峭度通常在 5 以上。相关系数则反映模态与原始信号的相似度太低说明这个模态主要来自噪声。K size(imf, 1); for k 1:K ku kurtosis(imf(k, :)); co abs(corrcoef(imf(k, :), x)); fprintf(模态%d: 峭度%.2f, 相关系数%.3f\n, k, ku, co(1,2)); end结合两个指标选模态峭度要高相关系数不能太低。只出现峭度高但相关系数很低的模态时先怀疑 K 偏大导致信号被拆碎而不是直接采信该模态。5.3 残余能量比衡量分解完整度把全部模态加起来与原始信号做差计算残余能量占比。这个指标直接反映分解是否损失了有效成分阈值的经验范围是 1%5%。resid x - sum(imf, 1); ratio sum(resid.^2) / sum(x.^2); fprintf(残余能量占比%.4f\n, ratio);tau 设为 0 时残余应该非常小除非迭代没收敛残余明显偏大时检查 tol 是否太宽松或者输入信号存在 NaN/Inf。tau 调大后残余能量占比允许略微上升但不应超过 5%。5.4 一条自动选择 K 的扫描思路把中心频率最小间隔、残余能量比、模态峭度串成一条自动扫描逻辑对 K2 到 8 分别执行 VMD计算每个 K 对应的最小中心频率间隔和最大模态峭度再绘制两者随 K 的变化曲线。曲线上的拐点就是合适的 K。这个思路类似 KMeans 聚类里用轮廓系数定类别数本质都是在分辨率与过拟合之间取平衡。扫描耗时可控时配合 alpha 在 8003000 范围内做 23 次粗扫基本能替代纯手工调参。把这套指标自动化的脚本保留下来后续换任何信号都能直接复用这是整个程序包里最值得沉淀的部分。本文还有配套的精品资源点击获取 VMD变分模态分解实战：从原理到轴承故障诊断
简介这是一份基于MATLAB的变分模态分解VMD算法程序包面向从事信号处理、故障诊断、数据分析的研究者与初学者。程序包含VMD核心函数、仿真信号脚本、包络谱与Hilbert分析工具以及多组轴承故障数据内圈、外圈等便于直接运行并观察非线性非平稳信号的分解与特征提取效果。压缩包共28个文件以m函数脚本、txt说明文档、mat数据文件和asv备份文件为主整体约9.62MB结构紧凑适合学习调试。目前已有1103人学习下载。通过调用VMD.m并调节模态数K与正则化因子α可灵活研究不同参数对分解结果的影响配套的仿真及实测数据脚本则展示了从信号加载、VMD分解到包络谱分析的一站式流程适合作为入门VMD算法与MATLAB实现的动手参考资料。1. 拆开 VMD 程序包之前先弄明白它和 EMD 差在哪处理滚动轴承振动信号、心电信号这类非线性非平稳数据时EMD经验模态分解是最先想到的工具但它有个老毛病模态混叠。一次冲击事件被拆到好几个 IMF 里包络谱上出现一堆假峰故障特征频率反而看不出来。VMD变分模态分解把「递归筛信号」改成「一次性求解约束优化问题」每个模态都被限制在各自的中心频率附近混叠问题被显著压制。这套 VMD.zip 里既有核心的 VMD.m也有 hua_fft1.m、hua_baoluo.m、hua_hilbert1.m 这些配套分析脚本还带着 1750内圈.mat、1750外圈.mat 等真实轴承故障数据适合做故障诊断、信号降噪、特征提取的工程师直接跑通一整套流程。接下来按「原理→文件→参数→实战→验证」的顺序逐一拆。2. VMD 的变分框架与程序文件逐项拆解2.1 变分问题如何构造带宽最小化与约束EMD 的思路是不断地找包络均值、剥出本征模态函数本质上是递归筛选误差会逐层累积。VMD 换了个思路先假设原始信号是由 K 个有限带宽的模态叠加而成每个模态都有一个中心频率然后构造一个约束优化问题——目标函数是全部模态的估计带宽之和最小约束条件是所有模态加起来等于原始信号。数学形式可以写成对每个模态 uk(t)先做 Hilbert 变换得到单边频谱再乘上指数项 e^{-jωk t} 把频谱搬到基带接着计算梯度范数的平方作为带宽估计。整个优化问题如下min Σₖ || ∂t [ (δ(t) j/πt) * uk(t) ] · e^{-jωk t} ||²约束 Σₖ uk(t) f(t)这里 ∂t 表示对时间求导* 是卷积ωk 是第 k 个模态的中心频率。直接求解带约束的变分问题不现实VMD 的做法是引入二次罚项 β 和拉格朗日乘子 λ(t)把约束罚进目标函数再通过交替方向乘子法迭代更新 uk、ωk 和 λ。每次迭代先固定其他变量更新第 k 个模态再固定模态更新对应中心频率最后更新乘子。迭代到满足收敛容差时停止输出的 uk 就是分解结果。理解了这一点你再看 VMD.m 里那些循环和矩阵运算就不会懵。程序里反复出现的u_hat、omega分别对应频域中的模态估计和中心频率迭代终止条件是前后两次更新的差值小于用户指定的 tol。这个框架天然支持对分解结果的数学解释而不是像 EMD 那样依赖经验筛选这也是它适合二次开发、能接进深度学习和故障诊断流程的原因。2.2 VMD.zip 内文件的角色划分拿到压缩包后不要急着运行先把文件清单过一遍搞清楚每个脚本在整条链路里的位置。下表是我拆包时整理的对应关系文件角色典型用途VMD.m核心分解函数输入信号和参数输出各模态及中心频率fangzhen.m / fangzhen.asv仿真入口脚本构造仿真信号、调用 VMD、画图仿真信号.m信号生成脚本定义 fs、时间轴和叠加分量hua_fft1.mFFT 频谱分析查看原始信号频域特征hua_hilbert1.mHilbert 变换求解析信号、瞬时幅值/相位hua_baoluo.m / hua_baoluo.asv包络谱分析解调出冲击成分特征频率Untitled.m测试/主控脚本串联数据读取、分解、绘图1750内圈.mat、1750外圈.mat故障数据轴承内圈/外圈故障振动信号72内.mat、1772内圈.mat附加数据其他工况或测点的样本6.txt、7.txt、8.txt、9.txt 等文本格式信号用于 load 直接读取的原始序列这些 .txt 文件通常是单列数值适合先用load(9.txt)读进来做快速验证.mat 文件则可能包含多个变量读入后先whos看一眼变量名和尺寸别默认变量名就叫 data。hua_baoluo.m和hua_hilbert1.m是配套的前者做包络谱时需要调用 Hilbert 变换得到包络信号后者把这一步封装好了单独拎出来也能用于瞬时特征提取。2.3 仿真信号入口脚本怎么读直接看fangzhen.m是最快的上手方式。常见做法是脚本先构造一个多分量信号比如低频正弦叠加高频调幅分量再加一点噪声然后调用 VMD.m。如果包里的仿真脚本没有按这个结构组织下面这个模板等价于它做的事fs 1000; N 1000; t (0:N-1)/fs; % 三个分量低频、中频调幅、高频正弦 x1 0.5 * cos(2*pi*20*t); x2 (1 0.3*cos(2*pi*5*t)) .* sin(2*pi*85*t); x3 0.3 * sin(2*pi*200*t); x x1 x2 x3 0.1*randn(size(t)); % 调用 VMDK3 对应三个分量 [u, ~, omega] VMD(x, 2000, 0, 3, 1e-7);这里 VMD 函数的形参按常见实现排列为信号、惩罚因子 alpha、噪声容忍度 tau、模态个数 K、收敛容差 tol。omega 返回的是归一化角频率要转成物理频率需要乘 fs/(2π)。如果某个文件的调用顺序不一样打开 VMD.m 看函数声明行就能确认不要凭记忆传参。跑通这个仿真后你就能验证 K 取 3 时三个模态是否对应 20Hz、85Hz、200Hz 三个分量这是检验程序包是否完整的第一步。3. 分解效果由参数决定alpha、K 的设置与副作用3.1 alpha 控制的是带宽不是分解个数alpha 是二次罚项的权重数值越大对模态带宽的限制越强每个模态被约束得越窄频域上各模态的重叠越少。这看起来是好事但 alpha 过大会让模态丢失细节把原本该包含的边频带切掉alpha 过小则模态带宽过宽两个模态会在频域上交叉出现类似 EMD 的混叠现象。对轴承故障信号故障冲击会产生一系列以共振频率为中心的边频带alpha 太大把这些边频削掉后包络谱里就只剩下孤零零的特征频率峰调制信息丢失。反过来alpha 太小则噪声会作为独立模态被分解出来。常见做法是先固定 alpha2000 跑一遍再看模态的频谱是否有明显重叠重叠明显就往大调边频缺失就往小调。实际项目中 alpha 落在 50010000 这个区间都是正常的不要迷信默认值。3.2 K 的选择方法中心频率观测与相关性校验K 是模态个数它的影响比 alpha 更直接。K 太小多个频率成分被塞进一个模态里模态频谱呈多峰形状K 太大会把单成分信号硬拆成两段出现频域相邻的虚假模态。判断 K 是否合适的核心手段是看中心频率的分离情况逐步增大 K把每次分解得到的中心频率打出来如果新增加的那个中心频率与已有某个中心频率非常接近说明 K 已经过大了。相关性校验也常用。分解出的每个模态与原始信号做 Pearson 相关相关系数高的模态代表真实成分相关系数突然跌到接近 0 的模态基本都是虚假分量。比如 K5 时第 5 个模态与原始信号的相关系数只有 0.02而前 4 个都大于 0.3那基本可以判定 K 应该取 4。注意相关性强弱和能量大小不是一回事一个高频小幅值成分相关系数可能不高但它确实是物理信号的一部分所以要结合频谱看不能只盯相关系数这一个指标。3.3 VMD.m 的调用方式与返回结果包里 VMD.m 的典型签名是[u, u_hat, omega] VMD(signal, alpha, tau, K, tol)其中 signal 是 1×N 的行向量输出 u 是 K×N 的矩阵每一行是一个模态u_hat 是对应模态的频域表示长度与信号相同omega 是 K 个中心频率的归一化角频率。tau 是噪声容忍度设为 0 表示完全强制重构精度设为非零值则允许一定残留噪声实测含噪信号时 tau 取 0.010.5 能避免把噪声强行掰成模态。% 对信号做预处理去均值 归一化 x x(:); x x - mean(x); x x / max(abs(x)); % 分解得到 4 个模态 [imf, ~, omega] VMD(x, 2000, 0.1, 4, 1e-7); % 换算成物理频率并打印 fs 1000; fprintf(中心频率(Hz): %.2f %.2f %.2f %.2f\n, ... omega*fs/(2*pi));预处理非常重要。去均值是防止直流分量占据一个模态归一化是避免幅值过大导致迭代不收敛。如果 VMD 分解后重构信号与原始信号误差很大优先检查输入是否做了这两步。另外 VMD 对数据长度敏感点太少解不出准确的带宽做短样本分析时至少保证 N 500否则结果不稳定。3.4 参数速查表参数常见范围经验初值表现出的问题与调整方向alpha500100002000模态频谱重叠→增大边频/细节丢失→减小K284模态多峰→增大出现相邻中心频率→减小tau00.50无噪声、0.1含噪分解结果忽大忽小→适当增大 tautol1e-91e-51e-7收敛太慢→放宽结果不收敛→调严这几个参数不是独立起作用的。alpha 和 K 共同决定模态的频域划分方式改 K 之后最好重新审视 alpha 是否还合适。调参的顺序我一般固定为先定 K靠中心频率观察再调 alpha靠频谱重叠和边频保留情况最后处理 tau 和 tol 这两个次要参数。4. 滚动轴承故障数据实战VMD 分解加包络谱识别4.1 读取 .mat 数据并做预处理包里的 1750内圈.mat、1750外圈.mat 是滚动轴承故障信号1750 通常对应 1750 r/min 的工况转速。先用 whos 查看变量再用plot看波形的冲击特征注意数据可能是列向量。% 载入内圈故障数据 load(1750内圈.mat); whos; % 确认变量名和尺寸 % 假设变量名为 data x data(:); fs 12000; % 西储大学轴承数据的常见采样率按实际填写 x x(1:fs*2); % 取前 2 秒数据避免数据量太大 % 去除直流和趋势项 x x - mean(x);如果 .mat 文件里装的是多个测点的信号选通道时优先选振动冲击最明显的那个。先做去均值和线性趋势去除原因很简单VMD 把趋势项当成低频模态输出浪费一个分解名额后续包络分析还要再剔除一次。数据量过大时 VMD 迭代速度明显变慢取 2 秒甚至 1 秒数据足够看到故障特征频率没必要把整段数据一次性丢进去。4.2 计算内圈与外圈的故障特征频率故障诊断的目标不是看波形有多乱而是找到冲击重复的频率。对内圈故障滚动体每经过一个缺陷点都会产生一次冲击冲击频率等于内圈故障特征频率 BPFI外圈故障对应 BPFO。计算公式如下BPFI n · fr · (1 d/D · cosφ) / 2BPFO n · fr · (1 − d/D · cosφ) / 2其中 n 是滚动体数量fr 是转频d 是滚动体直径D 是节圆直径φ 是接触角。这些几何参数属于具体轴承型号比如 6205 深沟球轴承的 n9接触角近似 0内圈特征频率系数约为 5.415外圈约为 3.585。nm 1750 r/min 时转频为 1750/60 ≈ 29.17 Hz代入系数fr 1750 / 60; % 转频 n 9; % 滚动体个数 BPFI n * fr * 5.415 / 9; % 换算到对应系数 BPFO n * fr * 3.585 / 9; fprintf(BPFI%.2f Hz, BPFO%.2f Hz\n, BPFI, BPFO);如果不知道轴承几何参数也可以用转速和故障特征系数近似估算分析时在包络谱上找 BPFI 或 BPFO 的整数倍频即可通常前 4 阶谐波都会出现。4.3 执行 VMD 分解并绘制模态对内圈故障信号K 取 4 比较合适一个低频旋转分量、一个共振带、一两个噪声分量。执行分解并观察各模态时域波形和频谱。K 4; alpha 1200; % 内圈故障边频带丰富alpha 偏小保留边频 [imf, ~, omega] VMD(x, alpha, 0.05, K, 1e-7); % 观察中心频率 disp(omega * fs / (2*pi)); % 画前 4 个模态的时域波形 for k 1:K subplot(K, 1, k); plot((0:length(imf(k,:))-1)/fs, imf(k,:)); ylabel([IMF, num2str(k)]); endalpha 取 1200 而不是默认 2000原因是内圈故障的调制边带比较宽太大的惩罚因子会把边频压掉。分解后重点看冲击最明显的模态——时域波形上呈现稀疏的大幅值脉冲这类模态往往对应故障冲击激起的共振成分。如果多个模态都有冲击用峭度来选峭度是四阶统计量冲击性强的信号峭度远大于 3。4.4 包络谱分析定位故障选定冲击明显的模态后做包络谱先求 Hilbert 变换得到解析信号取模得到包络对包络去均值后再做 FFT。hua_baoluo.m 封装的正是这个流程等价实现如下imf_sel imf(1, :); % 按峭度筛选出的模态 env abs(hilbert(imf_sel)); env env - mean(env); Nfft 2^nextpow2(length(env)); Y fft(env, Nfft); P abs(Y(1:Nfft/2)) * 2 / Nfft; f_axis (0:Nfft/2-1) * fs / Nfft; plot(f_axis, P); xlim([0 500]); xlabel(频率 (Hz)); ylabel(幅值);将峰值频率与 BPFI 对比如果峰值落在 BPFI 或其 2 倍、3 倍频附近并且谱峰明显高于周围噪声底就可以判定为内圈故障。外圈故障数据用同样的代码切换 BPFO 验证。注意包络谱前几个低频点经常有转频及其谐波这些不是故障特征区分方法是看它是否等于 fr 的整数倍而不是 BPFI/BPFO 的整数倍。VMD 分解在这个流程中的作用是把故障冲击模态从强噪声和旋转分量中分离出来让包络谱不被无关频率成分干扰。5. 用三个指标验证分解质量再谈参数自动寻优5.1 中心频率分离度判断 K 是否过大每次调完 K打印中心频率比较相邻中心频率的最小间隔。正常分解下各中心频率间距应明显大于 0.3 倍带宽如果最小间隔小到和带宽一个数量级说明两个模态混叠了。for K 2:8 [~, ~, omega] VMD(x, 2000, 0.05, K, 1e-7); cf sort(omega * fs / (2*pi)); min_gap min(diff(cf)); fprintf(K%d, 中心频率[%s], 最小间隔%.2f Hz\n, ... K, num2str(cf, %.1f ), min_gap); end观察哪个 K 对应的最小间隔突然缩小一个数量级K 的上限就在那个值之前。这是我在实际项目中用的第一道检查线比肉眼看频谱可靠。5.2 峭度与相关系数筛出有效模态分解完成后用峭度筛选冲击模态用相关系数剔除无效模态。峭度大于 3 表示信号比正态分布更尖锐故障冲击模态的峭度通常在 5 以上。相关系数则反映模态与原始信号的相似度太低说明这个模态主要来自噪声。K size(imf, 1); for k 1:K ku kurtosis(imf(k, :)); co abs(corrcoef(imf(k, :), x)); fprintf(模态%d: 峭度%.2f, 相关系数%.3f\n, k, ku, co(1,2)); end结合两个指标选模态峭度要高相关系数不能太低。只出现峭度高但相关系数很低的模态时先怀疑 K 偏大导致信号被拆碎而不是直接采信该模态。5.3 残余能量比衡量分解完整度把全部模态加起来与原始信号做差计算残余能量占比。这个指标直接反映分解是否损失了有效成分阈值的经验范围是 1%5%。resid x - sum(imf, 1); ratio sum(resid.^2) / sum(x.^2); fprintf(残余能量占比%.4f\n, ratio);tau 设为 0 时残余应该非常小除非迭代没收敛残余明显偏大时检查 tol 是否太宽松或者输入信号存在 NaN/Inf。tau 调大后残余能量占比允许略微上升但不应超过 5%。5.4 一条自动选择 K 的扫描思路把中心频率最小间隔、残余能量比、模态峭度串成一条自动扫描逻辑对 K2 到 8 分别执行 VMD计算每个 K 对应的最小中心频率间隔和最大模态峭度再绘制两者随 K 的变化曲线。曲线上的拐点就是合适的 K。这个思路类似 KMeans 聚类里用轮廓系数定类别数本质都是在分辨率与过拟合之间取平衡。扫描耗时可控时配合 alpha 在 8003000 范围内做 23 次粗扫基本能替代纯手工调参。把这套指标自动化的脚本保留下来后续换任何信号都能直接复用这是整个程序包里最值得沉淀的部分。本文还有配套的精品资源点击获取