VMD变分模态分解原理及Python实现：轴承故障特征提取与参数调优

简介一份面向机械故障诊断与信号处理学习者的VMD复现资源对应文献《基于VMD的故障特征信号提取方法》。代码包通过MATLAB脚本完整演示了振动模态分解在故障特征提取中的实现路径涵盖信号预处理、模态分解、频谱分析及特征可视化等环节适合具备一定信号处理基础和MATLAB编程能力的读者参考学习。压缩包共4个文件均为.m脚本整体仅5KB其中核心函数负责VMD分解的迭代优化与正则化过程主程序串联信号输入、分解调用与结果绘图另有辅助脚本用于功率谱计算等细节分析代码结构清晰、便于逐段研读和二次修改。目前已有731人学习下载代码包虽小但胜在简洁完整能帮助快速理解VMD从原理到落地的关键步骤复现过程不涉及具体轴承实验数据侧重理论算法实现通过运行代码可直观观察非平稳信号被分解为频率局部化模态分量的效果学会从噪声中剥离故障特征信息为后续设备状态识别与预测维护提供方法支撑。1. VMD 从论文到故障特征提取先搞清复现的是哪一步转子或轴承出现早期故障时故障特征频率往往埋在一堆背景噪声和多个谐波分量里。直接对原始信号做快速傅里叶变换频谱上通常找不到一条干净的故障谱线因为冲击成分的能量被展布在一个较宽的频带中转频及其倍频又在一旁“抢镜头”。变分模态分解VMD在这里的价值是把一个宽带信号自适应拆成若干个窄带模态每个模态围绕各自的中心频率展开故障特征所在的频带会被单独分离出来后续包络谱分析才有得可看。复现这篇文献不只是把论文公式抄成代码而是要把“分解—选模态—包络谱—找特征频率”这整条链路跑通。这条链路对做旋转机械故障诊断的工程师和正在读论文想落地验证的研究生都适用下面按我实际会用的顺序把每个环节说清楚。2. 用 Python 复现 VMD 最小实现核心迭代、代码与参数表2.1 复现前先看明白变分模型在算什么VMD 把观测信号 f 分解为 K 个本征模态函数 u_k每个模态被建模成一个调幅调频信号。目标函数是对每个模态的解析信号做频率搬移后求梯度二范数也就是在衡量模态带宽。约束条件是所有模态之和等于原始信号。写成优化问题就是min ∑_k ‖ ∂_t [ (δ(t) j/(πt)) * u_k(t) ] e^(−jω_k t) ‖₂²s.t. ∑_k u_k f∂_t 是对时间求导(δ(t) j/(πt)) 与 u_k 做卷积相当于构造解析信号再乘 e^(−jω_k t) 相当于把频谱搬到以中心频率 ω_k 为原点的坐标系里看。目标函数越小说明每个模态越窄带也就越接近“单一特征成分”。这个带约束优化问题用拉格朗日乘子和交替方向乘子法ADMM求解。一轮迭代里依次做三件事更新每个模态的频域表达、重估中心频率、更新拉格朗日乘子。写代码时不需要把整个拉格朗日函数展开记住三个更新式子就够用。模态谱的更新是一个维纳滤波结构分母是 1 2α(ω − ω_k)²α 越大偏离中心频率越远的分量被压得越狠模态就越窄中心频率按模态功率谱的重心重估拉格朗日乘子按分解残差累加。这三个式子互相耦合所以迭代收敛后得到的结果既能保持窄带特性又能尽量还原原始信号。理解了这一点对应到代码里每一行的作用就不会看晕。2.2 一个可复现的 numpy 版本复现 VMD 的常见做法是照着 ADMM 迭代自己写一版也可以移植论文作者公布的 MATLAB 实现。下面这个 numpy 版本把主循环压缩到 30 行以内适合先跑通再逐步扩展。信号长度建议取偶数避免奈奎斯特频率处共轭对称时丢点。import numpy as np from numpy.fft import fft, ifft def vmd(signal, K4, alpha2000.0, tau0.0, DCFalse, init1, tol1e-7): N len(signal) half N // 2 1 # 单边频谱点数 fhat fft(signal)[:half] # 单边原信号频谱 omega np.linspace(0, np.pi, half) # 数字频率轴 u_hat np.zeros((K, half), dtypecomplex) omega_k np.zeros(K) if init 1: # 中心频率在频带内均匀铺开 omega_k np.linspace(0.1 * np.pi, 0.9 * np.pi, K) else: # 所有模态中心频率从 0 开始 omega_k np.zeros(K) lam np.zeros(half, dtypecomplex) prev u_hat.copy() for it in range(500): for k in range(K): # 去掉其他模态贡献后的残差加上拉格朗日乘子的一半 residual fhat - np.sum(u_hat, axis0) u_hat[k] lam / 2.0 # 维纳滤波结构alpha 控制带宽 u_hat[k] residual / (1.0 2.0 * alpha * (omega - omega_k[k]) ** 2) # 按功率谱重心重估中心频率 power np.abs(u_hat[k]) ** 2 omega_k[k] np.sum(omega * power) / max(np.sum(power), 1e-12) if DC: # 直流分量单独固定到真实值 u_hat[0, 0] fhat[0] lam tau * (fhat - np.sum(u_hat, axis0)) eps np.sum(np.abs(u_hat - prev) ** 2) / max( np.sum(np.abs(prev) ** 2), 1e-12) if eps tol: break prev u_hat.copy() # 单边谱镜像回双边再反变换得到时域模态 u [] for k in range(K): full np.zeros(N, dtypecomplex) full[:half] u_hat[k] full[half:] np.conj(u_hat[k][1:half - 1][::-1]) u.append(np.real(ifft(full))) return np.array(u), u_hat, omega_k, it这段代码里有几个细节值得注意。更新某个模态时从 fhat 中减掉的是“其他模态当前的值”也就是说在同一个迭代轮次里先更新的模态会立刻影响后面模态的残差计算这种 Gauss-Seidel 风格的顺序更新比全部算完再统一更新的收敛速度更快也是复现版本里最常见的写法。中心频率重估用 np.sum(omega * power) 做加权平均物理含义就是找功率谱的重心不是找峰值位置。最后返回的三个结果里u 是时域模态u_hat 是单边频域模态omega_k 是迭代结束后的中心频率后面判断分解好坏全靠 omega_k。2.3 VMD 参数表新手照着设熟手按需改参数作用常见设置调参倾向K模态数量4 ~ 6偏大出现模态混叠偏小会漏频带alpha带宽惩罚因子500 ~ 3000偏大模态窄、抗噪强、易丢冲击细节偏小模态宽、易相互合并tau噪声容限0 或 fs/100含噪大时给正数可平滑重建信号但不是越大越好DC是否单独提取直流False信号有趋势项或零漂时置 Trueinit中心频率初始位置1均匀先均匀铺开看迭代结果不要直接猜故障频率tol收敛精度1e-7调大则提前停止中心频率会偏移调小迭代久K 和 alpha 是相互作用的一对参数。alpha 增大后每个模态带宽收窄同样的频带覆盖范围需要更多模态才够所以调 K 之前先固定 alpha不要两个参数同时乱调。tau 在原文默认取 0属“不允许重建误差存在”的硬约束实际带噪声信号复现时先保持 tau0 把主链路跑通再做噪声平滑时再给正数。2.4 用最小命令验证分解链路没写错先拿一个不含故障冲击的合成信号测试两个正弦加白噪声看 VMD 能不能把三个模态分清楚。这样能隔离“代码 bug”和“参数不合适”两类问题。fs 2000 t np.arange(0, 1, 1 / fs) test_signal (np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 200 * t) 1.2 * np.random.randn(len(t))) u, u_hat, omega_k, it vmd(test_signal, K3, alpha1500, tau0, DCFalse, init1, tol1e-7) print(中心频率(Hz):, omega_k * fs / (2 * np.pi))omega_k 是数字角频率范围在 0 到 π 之间换算成物理频率要乘 fs / (2π)。理想输出里三个中心频率应接近 50、200 和一个落在白噪声较高频段的数值。如果看到两个模态挤在同一个频率附近说明 K 给大了或 alpha 给小了。如果五次运行结果每次都不一样检查生成噪声时有没有固定随机种子VMD 本身没有随机性结果波动只来自输入信号。3. 复现文献实验流程构造轴承故障信号并用包络谱验证特征3.1 为什么先构造仿真信号而不是直接上实测数据直接拿现场采集的轴承数据复现最大的问题是不知道“正确答案”。故障特征频率理论值可以由轴承参数算出但现场转速波动、传递路径衰减、多点故障叠加都会让包络谱峰值不明显出了问题很难分清是 VMD 参数没选好还是数据本身太复杂。先用参数完全已知的合成信号把链路验证通再换实测数据才能定位问题出在哪一步。合成信号由三部分组成30 Hz 转频加上一倍频模拟转子不平衡周期性指数衰减冲击模拟轴承外圈局部故障高斯白噪声模拟背景干扰。冲击重复频率就是外圈故障特征频率 BPFO我们提前知道它等于 123 Hz最后用包络谱结果反推验证这就是复现文献里最常用的一种“闭环验证”做法。3.2 构造含冲击特征的一秒钟仿真信号import numpy as np from numpy.fft import fft fs 8192 N fs * 1 t np.arange(N) / fs # 转频及其二倍频模拟不平衡成分 rotor np.sin(2 * np.pi * 30 * t) 0.4 * np.sin(2 * np.pi * 60 * t) # 外圈故障特征频率理论给定 123 Hz BPFO 123.0 # 脉冲位置按时间轴计算取整到最近采样点 positions np.round(np.arange(0, 1, 1 / BPFO) * fs).astype(int) positions positions[positions N] impulse_train np.zeros(N) impulse_train[positions] 1.0 # 冲击响应指数衰减 2500 Hz 固有振荡模拟故障冲击波形 imp_len int(0.08 * fs) imp np.exp(-np.arange(imp_len) / 40) imp * np.sin(2 * np.pi * 2500 * np.arange(imp_len) / fs) fault np.convolve(impulse_train, imp, modesame) signal rotor 0.8 * fault / np.max(fault) 0.25 * np.random.randn(N)脉冲位置用 np.arange(0, 1, 1/BPFO) 生成的时间点乘以采样率再取整脉冲间隔会在 66 和 67 个采样点之间交替平均重复频率正好是 123 Hz。这在频谱上会表现为 123 Hz 及其倍频处有谱线只是能量会向相邻频点轻微扩散。冲击响应长度取 0.08 秒衰减系数 40 决定冲击的持续时间2500 Hz 是结构固有频率实际轴承故障冲击的高频共振频率通常在 2000 到 5000 Hz 之间取这个量级比较有代表性。0.8 是冲击分量的相对幅值0.25 是噪声标准差信噪比大约 10 dB 上下和现场测得的轴承信号强度相近。3.3 对信号做 VMD 并把频带与中心频率对应起来u, u_hat, omega_k, it vmd(signal, K5, alpha1500, tau0, DCFalse, init1, tol1e-7) freqs np.linspace(0, fs / 2, N // 2) for k in range(5): spec np.abs(fft(u[k]))[:N // 2] peak_freq freqs[np.argmax(spec)] center_hz omega_k[k] * fs / (2 * np.pi) print(f模态{k}: 频谱峰值 {peak_freq:.1f} Hz, 中心频率 {center_hz:.1f} Hz)运行后打印结果重点观察两列数值的关系。中心频率是迭代算出来的模态能量重心频谱峰值则是该模态中最强的单一频率。故障冲击模态的能量落在一个较宽的频带内中心频率可能停在 1800 Hz 附近而频谱峰值则出现在冲击振荡频率 2500 Hz 附近并不在 123 Hz 上。所以不要用“哪个模态的中心频率最接近 123 Hz”来选模态正确做法是对每个候选模态做包络谱之后再判断。这里看到的中心频率如果能和频谱峰值对应上说明模态分离没有出现交叉。3.4 包络谱里找 BPFO 及倍频选定冲击模态后做希尔伯特包络谱这一步是“故障特征信号提取”的落点。from scipy.signal import hilbert k_best 2 # 根据上一步打印结果选冲击信号所在的模态 env np.abs(hilbert(u[k_best])) env_spec np.abs(fft(env))[:N // 2] env_f np.linspace(0, fs / 2, N // 2) top_idx np.argsort(env_spec)[-5:] for i in top_idx: print(f{env_f[i]:.1f} Hz : {env_spec[i]:.2f})包络谱的原理是冲击性故障在时域上表现为重复的瞬态希尔伯特变换取模得到包络信号包络的频谱在冲击重复频率处会出现谱线。所以理论 BPFO 是 123 Hz包络谱峰值里应当出现 123、246、369 Hz 三条线。由于脉冲位置存在采样点取整123 Hz 处会看到一个小峰群取峰群中心位置做判断。包络谱实测峰理论位置判定约 123 HzBPFO 基频主峰必须明显约 246 Hz2 × BPFO应当存在约 369 Hz3 × BPFO弱峰或可见即可VMD 相对直接对原始信号做包络谱的优势在这一步体现得最明显。原始信号里 2500 Hz 的冲击振荡与 30 Hz、60 Hz 转频分量在包络谱中会互相干扰123 Hz 的峰可能被转频调制边带淹没。VMD 先把冲击模态单独分离再做包络谱噪声底更低123 Hz 峰与周围幅值对比更强。如果打印结果里 246 Hz 看不到多半是 alpha 太小导致模态带宽过大把相邻频率成分粘进来了。4. VMD 调参与排错K、alpha 和初始化怎么定4.1 K 的检查方法看中心频率是否“互相打架”K 是最容易拍脑袋给的参数。常见做法是从 K4 开始跑完后打印 omega_k按升序排列观察相邻中心频率的间距。如果两个中心频率之间的距离小于整个频带的 1% 左右说明这两个模态在挤同一个频带K 偏大。反过来如果某个模态频谱明显覆盖了另一个模态的频率范围且包络谱里原故障频率附近没有独立峰说明 K 偏小漏掉了频带。更系统的做法是做一个小循环让 K 从 2 变化到 8记录每个 K 下 omega_k 的最小间距和重构误差。当最小间距开始小于阈值时选前一个 K 作为最终值。注意这个循环只用来调试不要每次运行都做一遍否则 VMD 轮数会拖慢整体流程。提示中心频率是否重叠不能只看最后一次迭代结果要把迭代过程中的 omega_k 打印出来观察。只对比首尾两次容易漏掉模态交换因为 VMD 迭代过程中两个模态偶尔会交换顺序最终结果序是乱了但间距看起来“正常”。4.2 alpha 对带宽的影响和三个典型取值alpha 的物理含义是拉格朗日乘子对带宽约束的惩罚权重没有绝对最优值。文献复现时大多从 1500 或 2000 起步然后看包络谱的底噪和倍频清晰度再调整。三个典型取值对应三种场景alpha模态频谱形态适用场景主要风险500带宽较宽边带完整调制边带密集的齿轮箱信号相邻模态易合并1500频带和幅值平衡多数滚动轴承故障需配合 K 一起调3000带宽窄噪声抑制强强噪声下找单根故障谱线冲击能量被拆成多个模态调 alpha 时看包络谱的特征BPFO 峰旁边出现成对旁瓣说明模态带宽太宽把转频的调制边带包了进来应该加大 alpha如果 BPFO 峰找不到但中心频率稳定可能是 alpha 太大冲击能量被摊到两个相邻模态里应该减小 alpha。反复试两三次就能找到当前信号下的合适值。4.3 init、DC 和 tau 三个参数的排错细节init1 时中心频率在频带内均匀铺开迭代有可能收敛到局部最优。两个模态初始中心频率靠太近时迭代可能一直保持交叉而不交换顺序最后结果里 omega_k 不是按频率递增排列。遇到这种情况把 init 改成 0 全部从零开始重跑一次通常会得到另一组分解结果选模态间互相关系数更低的那组即可。DCTrue 只在信号有明显直流偏置时才需要。振动加速度信号经过零均值化处理后一般用不到。如果强行开启第一模态会被限制在 0 频附近当真实故障特征频率恰好低于 10 Hz 时会被这个设置坑掉。多数滚动轴承故障频率在几十到几百赫兹所以默认 False 更安全。tau0 是原文默认设置表示强约束分解残差为零。带噪数据上复现时给 tau 一个正数可以让重建信号更平滑但同时也会让每个模态不再严格累加回原信号。调试顺序建议先固定 tau0 把参数调好全部排错完成后再试 tau0 看包络谱底噪是否下降不要一开始就开这项否则会同时引入两个变量。4.4 一组通用排错流程按顺序执行以下检查比盲目组合 K 和 alpha 更快定位问题打印 omega_k确认所有中心频率单调递增且彼此间距合理。顺序乱掉先改 init间距过近则减小 K。把各模态累加得到重构信号与原始信号计算归一化均方误差。误差大于 1e-3 说明迭代没收敛增大最大迭代次数或把 tol 从 1e-7 调小到 1e-8。检查每个模态的频谱是否呈单峰形态。出现多个相隔很远的峰说明这个模态里混入了其他频带成分K 小了。对每个模态逐个做包络谱不要只查最像的那一个。有时故障模态的频谱峰值不高但包络谱峰值反而最明显只凭频谱选模态会漏掉它。全部通过后再上实测数据用轴承理论特征频率对照验证。这套流程把“两个参数对着调”变成按证据链定位每次改动只引入一个变量。复现论文时把这几步跑完再对比文献里的频谱图和包络谱图才算真正复现出结果而不是只看中心频率数量对上就结束。5. 验证特征提取结果的三个工程技巧5.1 用合成信号做盲测回代构造三到四组已知分量的混合信号把 VMD 输出和真值对齐再做相关性判断。先按中心频率从小到大排序然后计算每个真值分量和每个模态的相关系数相关系数大于 0.95 才认定该模态复现成功。这个回代测试能直接暴露一个问题中心频率张冠李戴。 例如已知分量 A 是 50 HzB 是 200 Hz如果最终模态 1 中心频率是 200 Hz、模态 2 是 50 Hz说明顺序乱掉了但不一定影响包络谱。把排序逻辑在代码里固定下来后面所有指标计算都基于排序后的索引避免人工看串行。5.2 用包络谱能量占比做量化指标单独比较高斯噪声底容易受频谱泄漏影响。定义一个更稳的指标取 0.9×BPFO 到 1.1×BPFO 频段内包络谱幅值之和除以 0 到 2000 Hz 内总幅值之和把这个比值作为故障特征强度指标。比值大于 0.3 说明故障特征非常明显小于 0.05 时要么故障不存在要么 K 和 alpha 还没调对。这个指标比单纯看峰高更稳因为它同时考虑了旁瓣和噪声底的贡献换数据时不会因为幅值缩放而失真。5.3 用 1/3 倍频程对比 VMD 前后的信噪比增益直接对比包络谱峰值与噪底容易受噪声估计方法影响常见做法是按 1/3 倍频程分带计算。以 123 Hz 为中心的倍频程带作为目标带在远离故障频率处取一段无特征频带作为噪底参考。target (env_f 110) (env_f 140) # BPFO 所在倍频带 noise (env_f 60) (env_f 90) # 无特征频带作噪底参考 gain 10 * np.log10(env_spec[target].sum() / env_spec[noise].sum())对 VMD 前后的包络谱分别计算这个增益。VMD 分解合理时目标带的信噪比增益比原始信号包络谱提升 6 dB 以上。这个提升量受 alpha 影响最大alpha 偏小提升量低alpha 偏大虽然单峰突出但冲击被展宽导致倍频丢失。所以最终参数不取增益最大的那一组而是取“增益超过 6 dB 且 2 倍频、3 倍频完整可见”的那一组。用这个双条件跑完再下结论说复现成功才站得住。本文还有配套的精品资源点击获取