简介一套用于复现《基于VMD的故障特征信号提取方法》文献算法的MATLAB代码包面向机械设备故障诊断与信号处理方向的研究者、工程师及研究生。代码包围绕VMD降噪与特征提取技术将复杂非平稳信号分解为多个频率局部化的模态分量帮助识别隐藏在噪声中的故障特征进而判断设备运行状态。包内共有四个MATLAB脚本文件压缩后体量仅约五千字节结构完整清晰包含VMD核心分解函数、主运行程序以及用于频谱分析和指标计算的辅助脚本方便逐模块阅读、调试和二次开发。目前已有七百三十一人学习下载适合具备信号处理基础与MATLAB编程能力、希望通过代码复现经典方法的读者。通过运行并研读这套代码可以清晰理解VMD分解中的迭代优化与正则化策略掌握从原始振动信号中提取特征分量的完整流程为实际工程中的故障诊断与预测性维护提供可靠的技术参考。1. 变分模态分解与故障特征信号提取为什么照着文献复现 VMD 总在第一步就翻车去年我拿到一篇学位论文标题就是“基于 VMD 的故障特征信号提取方法”参数表写得清清楚楚模态数 K6惩罚因子 α2000。我把这两个数字原样抄进代码跑出来的结果却是一堆中心频率挤在一起的怪波形包络谱里连一个像样的故障特征峰都找不到。问题不在论文而在变分模态分解VMDVariational Mode Decomposition本身是一个需要和你的数据“对齐”的变分框架K 和 α 只要没对上信号的成分结构后续提取特征全是白做。这篇笔记想给你的是一条能直接走通的路VMD 到底在算什么、最小复现代码怎么写、K 和 α 怎么定、复现时常在哪几个坑里翻车以及最后怎么用合成信号做闭环验证。适合在做机械设备故障诊断、电力暂态分析、生物医学信号处理的工程师照着文献做 VMD 时卡在参数或结果对不上的人。2. VMD 原理不啃论文把带约束变分问题拆成三个能看懂的步骤2.1 一分钟读懂 VMD 在干什么和 EMD 比它升级在哪VMD 的基本假设是一段一维振动信号 x(t) 由 K 个“有限带宽”的模态组成每个模态围绕自己的中心频率振荡。它要解的优化问题很简洁所有模态加起来能还原原信号同时每个模态的带宽之和最小。这是一个带约束的变分问题约束是“模态求和等于原信号”目标函数是“各模态带宽总和最小”。带宽在 VMD 里怎么算先对每个模态做希尔伯特变换得到解析信号再乘一个指数项把中心频率搬到基带附近最后用梯度 L2 范数的平方估计带宽。整个问题用交替方向乘子法ADMM迭代求解每轮迭代交替更新模态 u_k、中心频率 ω_k 和拉格朗日乘子 λ直到满足收敛容差。你不需要背住全部推导但必须记住一个结论VMD 是一次性把 K 个中心频率和 K 个模态同时解出来而不是像 EMD 那样逐层从高频筛到低频。EMD 的最大痛点是模态混叠基本不可控端点处三次样条拟合误差还会逐层污染后面的分量。VMD 把分解问题变成了一个可调参的优化问题K 控制模态数量α 控制带宽惩罚强度。这让它从一个“黑匣子”变成了一个你手上能拧的旋钮代价是你必须理解这两个旋钮否则结果比 EMD 还难看。2.2 最小复现用 vmdpy 把一个三频仿真信号拆开常见做法是用 vmdpy 这个开源库它是 VMD 作者 Matlab 思路的 Python 移植接口非常收敛。安装之后一个函数调用就能拿到分解结果。先构造一个三频信号的仿真数据跑通整个调用链路from vmdpy import VMD import numpy as np fs 1000 t np.arange(0, 1, 1/fs) # 仿真信号50Hz 正弦 120Hz 正弦 300Hz 正弦互不重叠 sig 0.6 * np.sin(2*np.pi*50*t) 0.4 * np.sin(2*np.pi*120*t) 0.2 * np.sin(2*np.pi*300*t) K 3 # 模态数和信号分量数一致 alpha 2000 # 惩罚因子控制每个模态的带宽 tau 0 # 噪声容忍度0 表示不做松弛 DC 0 # 不让第一个模态特殊化为直流分量 init 1 # 中心频率均匀初始化 tol 1e-7 # 收敛容差 # u 是分解出的 K 个模态omega 是每次迭代的中心频率 u, u_hat, omega VMD(sig, alpha, tau, K, DC, init, tol) print(归一化角频率:, omega[:, -1]) print(每个模态的能量占比:, [np.sum(u[k]**2) / np.sum(sig**2) for k in range(K)])逻辑说明这段代码把三个分量用 VMD 拆开u 的形状是 (K, N)每一行对应一个模态信号omega 的第二维是迭代次数最后一列是收敛后的中心频率。打印出来的中心频率是归一化角频率范围在 0 到 π 之间对应实际频率用 f omega * fs / (2*pi) 换算50Hz、120Hz、300Hz 分别对应 0.157、0.377、0.942 左右。参数说明K 必须等于或大于实际分量数小于会让两个分量挤进同一个模态大于会拆出虚假模态。alpha2000 是大多数场景的起步值它越大每个模态的带宽越窄越小越宽。tau 一般保持 0只有信号本身噪声极重时才调高。DC 在大多数振动信号场景设 0除非你要专门分离趋势项。2.3 怎么确认迭代真的收敛画中心频率曲线比盯损失值更直观很多人只取 omega 最后一列但 VMD 默认最大迭代次数有限如果没收敛最后一列也是“算完”的数值并不代表它是对的。我的习惯是把中心频率随迭代次数的曲线画出来一眼判断收敛质量import matplotlib.pyplot as plt for k in range(K): plt.plot(omega[k, :], labelfIMF{k1}) plt.xlabel(迭代次数) plt.ylabel(中心频率归一化角频率) plt.legend() plt.grid(True) plt.show()逻辑说明每条曲线代表一个模态的中心频率收敛过程。正常情况是前几十轮快速调整后段变得平直。如果某条曲线在后段还在明显振荡说明 alpha 或 K 设置不当或者信号本身不适合直接分解。这时哪怕 omega[:, -1] 打印出来数值“合理”也不能放心用。参数说明曲线在后 1/3 保持水平或波动小于 1% 就算收敛。若不收敛先加大 tol 到 1e-6 试试如果还振荡通常不是你计算精度不够而是 K 比真实分量多某两个模态在抢同一段频带。3. 复现文献的关键K 和 alpha 怎么定中心频率曲线说了算3.1 K 值别照抄文献先跑一遍中心频率扫描看“打架”的位置文献里写 K6那是它自己的信号成分决定的。你换轴承型号、换转速、换测点采集到的信号分量数完全不一样。把文献的 K 当真理是我见过最多人犯的错。正确做法是从 K2 开始往上扫打印每个 K 下收敛后的中心频率观察相邻值之间的距离for K in range(2, 8): u, u_hat, omega VMD(sig, alpha2000, tau0, KK, DC0, init1, tol1e-7) freqs_hz np.sort(omega[:, -1] * fs / (2*np.pi)) print(fK{K}: , .join([f{f:.1f}Hz for f in freqs_hz]))逻辑说明随着 K 增加频谱会被切得更细。K 偏小时两个相近分量被并成一个模态中心频率会落在两者之间K 偏大时多出来的模态会和相邻模态抢频带表现为两个中心频率非常接近比如 98Hz 和 103Hz 同时出现。参数说明我的经验判断标准是两个相邻中心频率的差小于较大频率的 10% 时就说明已经过分解了。比如 98Hz 和 103Hz 这种差 5% 不到明显是同一个频带被切开。此时把 K 减 1 再看一遍。注意中心频率接近是过分解的充分信号但不是唯一信号。就算中心频率拉得开还要看一眼对应模态波形是否畸变比如出现明显的高频毛刺或“假正弦”状这也要回退 K。3.2 alpha 是带宽旋钮从 2000 起步按 2 倍步进去试alpha 的名字叫惩罚因子直观理解是“每个模态允许有多宽的频带”。alpha 越大带宽越窄模态越接近纯单频alpha 越小带宽越宽越能保留冲击类信号的宽频成分但不同模态之间越容易重叠。alpha 范围带宽表现常见问题500~1000宽带宽保留冲击细节模态混叠中心频率分离但波形互相串扰2000默认带宽适中大多数旋转机械振动信号可接受4000~8000窄带宽谐波分得很干净冲击波形可能被削成近正弦包络谱模糊调 alpha 的顺序我一般会先固定 K然后按 500、1000、2000、4000 这样二倍步进去试比较每个 alpha 下目标 IMF 的包络谱。注意分解图好看不等于有用判断标准是后续包络谱里故障特征频率峰是否突出。alpha 太小包络谱会出现一排“兄弟峰”alpha 太大包络谱主峰变矮变宽。实际工程里我有 90% 的场景 alpha 落在 1000 到 3000 之间超过 5000 的场景非常少。3.3 用排列熵PE自动选 K把“看中心频率”变成一条可量化的曲线中心频率扫描靠肉眼总有点“玄学”。如果想让选 K 变得可复现可以用排列熵Permutation Entropy, PE来辅助也就是文献里常说的 PE-VMD 思路。排列熵衡量时间序列的随机性噪声分量熵高纯谐波熵低混叠模态熵介于两者之间。当 K 合适时各模态的排列熵相差不大且整体偏低K 过大时新增的虚假模态熵明显偏高。from itertools import permutations from math import factorial def permutation_entropy(x, m3, delay1): n len(x) perms list(permutations(range(m))) perm_to_idx {p: i for i, p in enumerate(perms)} cnt np.zeros(len(perms)) for i in range(n - (m - 1) * delay): pattern tuple(np.argsort(x[i:i m * delay:delay])) cnt[perm_to_idx[pattern]] 1 p cnt / cnt.sum() p p[p 0] return -np.sum(p * np.log(p)) for K in range(2, 7): u, _, _ VMD(sig, alpha2000, tau0, KK, DC0, init1, tol1e-7) pes [permutation_entropy(u[k]) for k in range(K)] print(fK{K}, 平均PE{np.mean(pes):.4f}, 各分量{np.array(pes).round(4)})逻辑说明PE 计算的是每个模态的排列熵m3 是最常用的嵌入维度delay1 表示用连续采样点。K 从 2 升到 3 时平均 PE 通常显著下降再继续增加 K平均 PE 趋于平稳或回升新增模态的单个 PE 值会明显高于老模态。取平均 PE 最低或回落点的 K 作为最终模态数。参数说明m 不要取太大序列长度有限时 m 太大统计不充分一般 m3 或 4。排列熵对幅值不敏感这既是优点也是缺点——它只关心序列的序关系两个幅值差异很大但形状相似的信号会得到相近的 PE 值。4. VMD 复现避坑模态混叠、端点毛刺与参数玄学的 6 个现场4.1 现象K 设 6 出来 4 个模态中心频率挤在 96~105Hz包络谱出现三条兄弟峰照抄文献 K6 时分解结果里两个或三个模态的中心频率会挤在一个很窄的频带内包络谱在这个频带附近出现多个峰值很难判断哪个是真正的故障特征频率。原因K 超过信号真实分量数过分解导致同一个频带被多个模态瓜分。ADMM 迭代时能量会随机分配到相邻模态里谁先收敛谁拿大头结果不可预测。解决用上一章的中心频率扫描脚本从 K2 开始逐个看找到中心频率开始“打架”的那个 K往后退一位。同时可以把 alpha 提升到 3000 左右窄带宽能让各模态更“守本分”。不要试图用后处理去合并模态那种做法又麻烦又不可靠。4.2 现象alpha 设 500低频冲击成分被高频模态吃掉BPFO 找不到了alpha 偏小时每个模态的带宽很宽低频和高频模态在频域里重叠。迭代过程中能量可能被某个高频模态“吸走”导致低频模态里只剩一点残渣包络谱里根本看不到外圈故障特征频率。原因alpha 太小模态带宽竞争失去约束。VMD 本质上按频带划分能量带宽重叠等于没有划分。解决把 alpha 提到 2000 以上重新分解。如果冲击信号本身频带很宽且含强噪声先对原始信号做带通滤波把分析频带限定在传感器共振频带附近再做 VMD。直接拿全频带信号硬分解任何参数都救不了。4.3 现象IMF 两端长出大尾巴包络谱低频出现一串毛刺分解出的第一个和最后一个模态在信号两端出现明显幅值放大包络谱低频段多出一堆没有物理意义的峰。原因VMD 在频域里迭代有限长度信号的边界不连续会影响希尔伯特变换结果。这不是 VMD 独有EMD 也有只是 VMD 的表现相对轻。解决分解前做镜像延拓把信号左右各拼一段翻转数据分解完成后只取中间的原始长度部分def mirror_extend(sig, ext_len): left sig[1:ext_len1][::-1] right sig[-ext_len-1:-1][::-1] return np.concatenate([left, sig, right]) ext_len int(0.1 * len(sig)) sig_ext mirror_extend(sig, ext_len) u_ext, _, _ VMD(sig_ext, alpha, 0, K, 0, 1, 1e-7) u u_ext[:, ext_len:ext_len len(sig)] # 截掉延拓段逻辑说明镜像延拓让边界处导数连续减小端点突变。ext_len 一般取信号长度的 5%~10%太长会引入虚假周期太短没效果。截取后边的模态前几十个点仍可能有轻微畸变做包络谱时可以用窗函数加权。4.4 现象中心频率曲线第一次迭代就跳到 3.0然后持续振荡打印 omega 曲线时某条曲线从初始值瞬间跳变到接近 π之后来回振荡tol 设多小都收敛不了。这种一般是数据本身的问题。原因信号没去均值幅值量纲太大比如振动加速度以 g 为单位ADMM 更新时梯度数值过大。也可能是 init 初始化方式和数据不匹配。解决分解前先对信号减去均值再除以标准差做标准化。vmdpy 的 init1 表示中心频率均匀初始化init0 表示全部从 0 开始这两个选项都试一次取收敛更平滑的那个。很多博客会把这两个初始化写反遇到结果异常时两个都跑一遍最稳妥。4.5 现象打印出来的中心频率是 0.157、0.377和文献里的 50Hz、120Hz 对不上文献里给的是 Hzvmdpy 输出的是归一化角频率两者差一个换算系数。第一次复现的人经常拿着这个数字怀疑自己代码写错。原因VMD 在频域计算时使用归一化角频率范围 0 到 π 对应 0 到 fs/2。不理解这一点就会在结果解读上卡住。解决统一用 f_hz omega[:, -1] * fs / (2 * np.pi) 换算成 Hz 再分析。画中心频率收敛曲线时横轴用迭代次数纵轴用换算后的 Hz这样和文献对比时不会错位。4.6 现象同样一段信号EMD 分得还行VMD 反而分出一条接近零的模态不是 VMD 不如 EMD而是你把趋势项、噪声和故障冲击一股脑喂了进去。VMD 会把微弱但占一个模态的噪声单独分出来也可能是 K 设大后多出的一条“空模态”。原因VMD 假设每个模态都有非零带宽纯随机噪声也满足这个假设。信噪比很低、K 又偏大时噪声会被“包装”成一个看似合理的模态。解决先做带通滤波或小波阈值降噪再确定 K。分解后如果出现接近全零的模态把 K 减 1 重跑不要硬留。滤波时注意保留目标故障频带否则等于白做。5. 从分解到特征值用峭度、相关系数和包络谱把故障特征频率揪出来5.1 先筛掉“废模态”峭度 相关系数双指标选 IMFVMD 输出 K 个模态不是每个都有诊断价值。实测信号里常有两个模态是噪声或趋势项。我的筛选习惯是用一个简单循环给所有模态打分from scipy.stats import kurtosis for k in range(K): r np.corrcoef(sig, u[k])[0, 1] kurt kurtosis(u[k], fisherFalse) # 经典峭度定义正态分布为3 print(fIMF{k1}: 相关系数r{r:.3f}, 峭度{kurt:.2f})逻辑说明相关系数衡量模态与原信号的线性相关程度能量占比高的模态 r 自然大峭度衡量波形冲击性故障冲击信号峭度远大于 3白噪声峭度接近 3。两个指标一起看r 高但峭度低说明这是大能量正弦成分峭度高但 r 低可能是噪声尖峰两者都高的才优先进入包络谱分析。参数说明具体的筛选阈值要参考背景噪声水平。噪声低时 r0.3 且峭度3 是比较稳的组合噪声大时把 r 阈值放宽到 0.1但峭度阈值不要放宽否则噪声模态混进来。5.2 包络谱对齐先算故障特征频率再在谱图里找峰别靠肉眼“看着像”选好 IMF 之后做包络谱分析。包络谱的原理是故障冲击会通过高频共振幅值调制表现出来用希尔伯特变换取出包络再对包络做 FFT故障特征频率就会在低频段形成谱峰。def bearing_freqs(fr, n_balls, d, D, contact_angle0): cos_a np.cos(np.radians(contact_angle)) bpfo 0.5 * n_balls * fr * (1 - d / D * cos_a) bpfi 0.5 * n_balls * fr * (1 d / D * cos_a) return bpfo, bpfi from scipy.signal import hilbert, find_peaks env np.abs(hilbert(imf)) spec np.abs(np.fft.rfft(env)) freqs np.fft.rfftfreq(len(env), 1/fs) peaks, _ find_peaks(spec, prominence0.05 * spec.max()) for p in peaks: if abs(freqs[p] - bpfo) 0.01 * bpfo: print(f在 {freqs[p]:.2f}Hz 找到外圈故障特征峰)逻辑说明bearing_freqs 里的 fr 是转频n_balls 是滚动体数量d 是滚动体直径D 是节圆直径接触角一般取 0 或查轴承手册。这些参数决定理论上的外圈 BPFO 和内圈 BPFI。之后用 find_peaks 找包络谱的显著峰和理论值做匹配。参数说明容差取理论值的 1% 到 2%。频率分辨率由 fs/N 决定N 是信号点数点数越多容差可以设得越小。如果峰值离理论值偏差超过 2%先检查 fr 是否算错再检查轴承参数是否查对最后才怀疑分解参数。5.3 成组确认基频、倍频、边频带三个位置一起看单峰匹配很容易被噪声尖峰骗过去。轴承故障谱的典型特征是故障特征频率及其 2 倍频、3 倍频同时转频会产生边频带也就是 f ± fr 的位置出现调制峰。判断故障时这三个位置成组出现才可靠。实际操作是把理论频率列表扩展开bpfo、2bpfo、3bpfo、bpfo-fr、bpfofr。用 find_peaks 找到所有显著峰后逐个检查这些位置附近是否有匹配峰。成组出现的判定比单峰出现的判定可靠得多尤其在转速波动、负载变化引起边带模糊时边频带的出现反而能从侧面印证故障诊断结论。6. 上真实数据前的最后一步合成信号闭环验证锁定参数真实数据没有标准答案调参调得再辛苦你也无法判断结果对错。所以我复现文献里 VMD 类方法时一定会先构造一个已知故障频率的冲击信号跑闭环测试参数能在仿真信号上找回设定频率才拿真实数据去碰运气。# 构造已知故障特征频率 97.3Hz 的周期性冲击信号 fs 8192 t np.arange(0, 2, 1/fs) imp np.zeros_like(t) for i in range(0, len(t), int(fs/97.3)): n np.arange(int(0.02*fs)) imp[i:ilen(n)] 0.8*np.exp(-60*n/fs)*np.sin(2*np.pi*800*n/fs) sig imp 0.01*np.random.randn(len(t)) u, _, _ VMD(sig, 2000, 0, 4, 0, 1, 1e-7)这段代码构造的是每隔约 10.3ms 出现一次的衰减振荡脉冲代表一个故障特征频率为 97.3Hz 的冲击序列。分解后选峭度最大的 IMF用第 5 章的包络谱代码找峰如果检出频率和 97.3Hz 的相对误差小于 1%就说明这套 K、alpha 组合对冲击型信号是有效的。我第一版复现就是直接拿实测数据调参调了三天都被噪声和边带干扰搞得“看着像又不敢确定”。改成合成信号闭环之后一个下午就把 K 和 alpha 锁定了再回到真实数据时只需要微调。之后的习惯是凡是从文献里复现 VMD 类方法第一步永远是构造一个已知特征的合成信号跑通整条链路再换真实数据。这个习惯帮我省掉了无数次无效调参。希望帮到你。本文还有配套的精品资源点击获取
原理、实测数据处理与排障指南)
