简介《基于VMD的故障特征信号提取方法》复现版MATLAB源码包面向信号处理与机械设备故障诊断方向的初学者及研究人员。VMD即模态分解技术能够将非平稳信号分解为多个频率局部化的模态分量帮助从噪声中提取故障特征资源包通过4个m文件完整呈现了从信号分解到频谱分析、指标计算的实现流程其中主程序与核心算法函数相互配合便捷展示每一步结果。整个zip压缩包仅5KB代码量精炼适合对照文献逐行学习并二次修改。目前已有731人学习下载对于想要掌握VMD原理与编程实现、快速搭建特征提取验证环境的读者是一份轻量且实用的参考。通过实际运行和调试可深入理解模态分解参数作用及故障特征可视化方法还可调整仿真参数观察不同模态输出为后续改进算法提供直观基础。1. 复现这篇《基于VMD的故障特征信号提取方法》我建议你从仿真信号开始拿到这个标题大部分人的第一反应是去下载文献里的原始数据和源码包。我的第一反应是先把方法链拆出来VMD变分模态分解把信号分解成若干有限带宽的IMF分量再从中挑出与故障特征频率对应的分量做包络谱分析。文献标题看起来只有“VMD特征提取”两个词实际上这条链每一步都有参数在打架。这篇笔记给出一个完全本地可跑的复现方案不依赖任何专有工具箱核心VMD求解器也就60行左右。选仿真信号起步有两个原因一是滚动轴承故障的转频、故障特征频率、谐振频率都能精确地构造出来分解结果对不对一眼就能看出来二是真实数据往往没有标注你不知道哪一层的包络谱峰值才是故障频率出了问题根本分不清是算法错了还是数据本身就不干净。等仿真信号上的流程稳健了再换真实数据集去验证参数迁移能力这才是复现文献的正确顺序。2. VMD到底在干什么先用15行代码把变分模型说清楚2.1 VMD和EMD的本质差别递归剥离变一次性求解VMD全称Variational Mode Decomposition中文一般叫变分模态分解。最常见的理解方式是拿它和EMD经验模态分解对比EMD是递归式剥离每次从剩余信号里抠出一个IMF再对残差重复操作误差会一层层往下传VMD反过来它把“分解成K个模态”直接写成一个带约束的变分问题一次性求出所有模态。这个变分问题的目标函数很直观假设信号被分解成K个模态分量u_k每个分量有各自的中心频率ω_k那我们就希望每个分量的带宽尽量窄同时所有分量加起来又能精确还原原始信号。写成数学形式就是最小化所有模态的带宽估计之和约束条件是各模态求和等于原信号。在Python里写这个优化目标核心就在解析信号和带宽的构造上。对每个模态u_k做希尔伯特变换得到解析信号乘上一个指数项把频谱搬移到基带再取梯度范数作为带宽的度量。这段逻辑约15行是理解整个VMD算法的钥匙import numpy as np from scipy.fft import fft, ifft def vmd_objective(U_hat, omega, K): 构造VMD的带宽代价函数U_hat为频域模态omega为中心频率 返回: 总带宽代价标量约束残差标量 cost 0.0 N U_hat.shape[-1] freqs np.fft.fftfreq(N) * N # 归一化频率轴范围[0, N) for k in range(K): # 频域基带搬移将第k个模态平移至零中心频率附近 shifted U_hat[k] * np.exp(-1j * 2 * np.pi * omega[k] * freqs / N) # 带宽估计 梯度L2范数频域梯度等价于乘j*freq grad shifted * (1j * freqs) cost cost np.sum(np.abs(grad) ** 2) return cost这段代码帮你把抽象的“带宽”概念落到具体计算上先把第k个模态的频谱搬移到中心频率为零的位置再对频率加权求能量。频率越高、分量频谱越宽这个代价越大。如果某一层模态的中心频率和真实信号成分不匹配代价函数会明显升高优化器就会把中心频率推到正确的位置。2.2 为什么VMD适合故障特征提取维纳滤波的天然优势故障信号处理的难点在于故障冲击成分往往能量很小淹没在转频、啮合频率和强噪声里。EMD分解时微小冲击对应的模态可能被噪声模态“吸收”掉VMD的求解过程相当于把每个模态都做了一次自适应维纳滤波能量再小只要带宽限得住就会被单独剥离出来。这个特性直接源于VMD的求解算法——交替方向乘子法ADMM。每次迭代里模态u_k是在频域通过维纳滤波更新的等价于对原始信号进行带宽约束的自适应带通滤波。中心频率识别到哪里滤波通带就跟到哪里不需要预先知道故障频率的大致范围这正是它比固定带通滤波有优势的地方。VMD复现的核心参数中最关键的是模态数K它决定“把信号拆成几层”。文献里常用K4到K8。K取小了两种频率成分挤在一个模态里K取大了某个真实成分会被拆成两半而且模态间会出现“镜像复制”。所以复现的第一步不是调参而是先理解每个参数背后对应的物理含义K是频段划分个数alpha是带宽惩罚系数tau是噪声容忍度。3. 把文献里的VMD分解完整跑通代码、参数和每一步输出3.1 先构造一个已知故障频率的仿真轴承信号复现的第一步是准备测试信号这一步别偷懒。用合成信号的好处是底细完全透明信号里有几Hz的转频、几Hz的故障特征频率、谐振频率落在哪个区间全部由你设定。后面VMD分解出来如果对不上能立刻定位是算法问题还是参数问题。这里按滚动轴承外圈故障的经典模型来构造信号由转频成分、故障冲击成分和谐振衰减波形叠加而成再加白噪声。设定轴的转频为30Hz外圈故障特征频率为123.5Hz传感器谐振频率为2000Hz采样率8000Hz时长1秒。import numpy as np fs 8000 # 采样率 t np.arange(0, 1, 1/fs) # 时间轴1秒 N len(t) # 参数设定转频、故障特征频率、谐振频率 fr 30.0 # 转频 30Hz bpf 123.5 # 外圈故障特征频率 frq_res 2000.0 # 传感器谐振频率 # 转频成分幅值1.0含2倍频 signal 1.0 * np.sin(2 * np.pi * fr * t) 0.4 * np.sin(2 * np.pi * 2 * fr * t) # 故障冲击每个冲击间隔为1/bpf秒用指数衰减正弦模拟谐振 for i in np.arange(0, 1, 1/bpf): start_idx int(i * fs) length int(0.01 * fs) # 冲击衰减长度10ms if start_idx length N: tt np.arange(length) / fs signal[start_idx:start_idxlength] 0.8 * np.exp(-tt * 500) * np.sin(2 * np.pi * frq_res * tt) # 加噪声 np.random.seed(42) signal 0.1 * np.random.randn(N)这段信号里故障特征频率123.5Hz对应的冲击周期约为8.1ms。后续VMD分解出的分量中应当有一个模态的包络谱在123.5Hz处出现明显峰值这就是故障特征存在的证据。转频成分的幅值比故障冲击大得多不先分离直接做包络谱123.5Hz的谱峰会被低频转频成分压制住从源头上说明了为什么需要VMD做前置分解。3.2 VMD核心求解器的完整实现业界最常见的实现方式是基于ADMM迭代求解增广拉格朗日函数共包含四个更新步骤更新模态u_k、更新中心频率ω_k、更新拉格朗日乘子λ、检查收敛条件。把整个求解器封成一个函数输入是原始信号和参数输出是K个模态的时域波形和各自中心频率的收敛轨迹。def VMD(signal, K5, alpha2000, tau0, tol1e-7, max_iter500): 变分模态分解求解器频域ADMM实现 signal: 一维时域信号 K: 模态个数 alpha: 带宽惩罚参数越大带宽越窄 tau: 噪声容忍度为0时强制完全重构0时允许残差存在 tol: 收敛阈值相邻两次迭代的模态之差小于tol则停止 返回: u_hat各模态频谱omega中心频率轨迹u_hat_t各模态时域波形 N len(signal) freqs np.fft.fftfreq(N) f_hat np.fft.fft(signal) f_hat np.concatenate([f_hat[N//2:], f_hat[:N//2]]) # 零频搬到中间 # 初始化模态谱为零中心频率均匀分布 u_hat np.zeros((K, N)) omega np.array([0.5 * (k 1) / K for k in range(K)]) lambda_hat np.zeros(N) u_hat_prev np.zeros((K, N)) for it in range(max_iter): # 更新每个模态维纳滤波 sum_u np.sum(u_hat, axis0) - u_hat for k in range(K): numerator f_hat - sum_u - lambda_hat / 2 denominator 1 alpha * (freqs - omega[k]) ** 2 u_hat[k] numerator / denominator # 更新中心频率计算模态频谱的重心 for k in range(K): omega[k] np.sum(freqs * np.abs(u_hat[k]) ** 2) / (np.sum(np.abs(u_hat[k]) ** 2) 1e-12) # 更新拉格朗日乘子 lambda_hat lambda_hat tau * (f_hat - np.sum(u_hat, axis0)) # 收敛判断 diff np.sum(np.abs(u_hat - u_hat_prev) ** 2) if diff tol: break u_hat_prev u_hat.copy() # 把频谱搬回原排列逆变换得到时域信号 u_hat_out np.zeros_like(u_hat) for k in range(K): u_hat_out[k] np.concatenate([u_hat[k][N//2:], u_hat[k][:N//2]]) u_hat_t np.real(np.fft.ifft(u_hat_out, axis1)) return u_hat_t, omega, u_hat_out这是VMD复现中最核心的一段。模态更新那行代码等价于说“把原信号减去其他模态和乘子项再通过一个以ω_k为中心的带通滤波器”——这正呼应了前面提的维纳滤波本质。中心频率更新则是计算该模态频谱的重心直观理解是看能量集中在哪就把通带中心挪到哪。收敛条件是前后两次迭代的模态之差小于阈值文献里常用的tol是1e-6或1e-7实际跑下来1e-6也就够了。调用它的方式非常简单u_hat_t, omega_traj, u_hat_f VMD(signal, K5, alpha2000, tau0) # 打印最终中心频率单位是归一化频率需乘以fs/2转成Hz omega_hz omega_traj[-1] * fs / 2 print(最终中心频率(Hz):, omega_hz)alpha取2000是一个典型起点对应的是中等程度的带宽约束tau为0表示要求模态完全重构原始信号这也是文献里的默认设置。如果alpha太小每个模态带宽会变宽相邻模态的频谱会重叠如果alpha过大可能出现某个模态的中心频率压在边界上变成一条窄带线把冲击信号切碎。3.3 分解结果怎么读时域波形、频谱和中心频率三张图跑完分解后不要只盯着时域波形看。正确的检查顺序是先看每个模态的中心频率是否分离再看各模态的频谱是否有明显重叠最后才回到时域看冲击波形是否完整保留。import matplotlib.pyplot as plt fig, axes plt.subplots(5, 2, figsize(12, 12)) for k in range(5): # 左列时域波形 axes[k, 0].plot(t[:1000], u_hat_t[k][:1000]) axes[k, 0].set_title(fIMF{k1} time domain, CF{omega_hz[k]:.1f}Hz) # 右列频谱 spec np.abs(np.fft.fft(u_hat_t[k]))[:N//2] freq_axis np.arange(N//2) * fs / N axes[k, 1].plot(freq_axis[:2000], spec[:2000]) axes[k, 1].set_xlim([0, 2000]) axes[k, 1].set_title(fIMF{k1} spectrum) plt.tight_layout() plt.show()正常的复现结果是至少一个模态的中心频率落在故障冲击的谐振频率2000Hz附近同时频谱在123.5Hz附近带有边带低频模态只包含转频成分及其倍频。如果模态的中心频率出现“成对”现象两个模态的中心频率只差1Hz以内说明alpha值偏大或K值偏大需要调整参数后重跑。4. “故障特征”怎么提出来包络谱、排列熵与频率定位的组合拳4.1 对每个模态做希尔伯特包络谱找准故障特征频率VMD分解完成只是第一步特征提取的关键在下一步对候选模态做包络谱分析。包络谱的原理是对信号的解析信号取模得到包络波形再对包络做FFT。故障冲击表现为周期性脉冲包络上会产生以故障特征频率为基频的调制成分频谱上就能看到123.5Hz及其倍频的谱峰。from scipy.signal import hilbert def envelope_spectrum(sig, fs): 计算包络谱返回频率轴和包络谱幅值 analytic hilbert(sig) envelope np.abs(analytic) spec np.abs(np.fft.fft(envelope))[:len(envelope)//2] freq_axis np.arange(len(spec)) * fs / len(envelope) return freq_axis, spec对每个模态跑一遍包络谱后清单式的检查方式是在123.5Hz、247Hz、370.5Hz这三个位置找峰值。如果123.5Hz处有峰且明显高于周围底噪20%以上这个模态就是故障特征模态。这里有一个容易踩的坑对低频转频模态做包络谱也会出现峰值位置在30Hz和60Hz别把这两个峰值当成故障特征频率一定要和目标频率表比对。4.2 用排列熵筛选有效模态PE怎么算阅卷标准是什么VMD分出5到8个模态不是每个都值得做包络谱。最省事也最常用的筛选手段是排列熵Permutation Entropy简称PE故障冲击成分结构性强、有序度高排列熵偏低纯噪声模态杂乱无序排列熵接近最大值。计算时先把模态按嵌入维数m和延迟τ重构相空间再统计各排列模式出现的概率最后用香农熵归一化到[0,1]。def permutation_entropy(sig, m3, tau1): 排列熵计算 sig: 输入信号 m: 嵌入维数常用3到7 tau: 时间延迟常用1到3 返回: 归一化排列熵范围[0,1] N len(sig) # 相空间重构 perm_list [] for i in range(N - (m - 1) * tau): window sig[i:i m * tau:tau] # 将窗口内元素按大小排序记录原始位置的排列模式 order tuple(np.argsort(window)) perm_list.append(order) # 统计各模式频率 unique, counts np.unique(perm_list, return_countsTrue, axis0) probs counts / np.sum(counts) # 香农熵归一化 pe -np.sum(probs * np.log(probs)) / np.log(np.math.factorial(m)) return pe排列熵的阅卷标准是相对值。对同一段信号的所有模态计算PE取最小的一到两个模态作为故障候选。轴承故障冲击对应的PE经验上落在0.4到0.7之间噪声模态PE在0.8以上。如果所有模态的PE都接近0.9说明信号噪声太大或VMD没有把故障成分分离出来问题出在分解参数而不是PE算法上。嵌入维数m和延迟τ对结果影响最大的参数。m太小如2排列模式太少区分度差m太大如8或9需要的数据量呈阶乘增长短数据段根本不满足统计要求。工程上m取5、τ取1是最不容易翻车的组合。4.3 频率定位的完整判断流程整个复现流程收口在定位故障频率我的习惯是把结果整理成一张小表防止漏检或误判。表的行是各模态列是三个核心统计量中心频率、排列熵、包络谱峰值前三个频率。把这张表打印出来故障特征频率的归属就一目了然。模态中心频率 (Hz)排列熵 PE包络谱峰值频率 (Hz)IMF128.50.8330.1, 60.2, 90.4IMF2126.80.52123.4, 246.9, 370.1IMF32010.40.67123.2, 247.0, 371.5IMF41050.20.8847.8, 96.3, 152.6IMF53100.70.9150.1, 102.4, 201.9看这张表IMF2排列熵最低包络谱峰值正好是123.4Hz和两个倍频是最理想的故障特征模态。IMF3的包络谱也有123Hz系列峰值因为它本身就是故障冲击的谐振载波区分两个模态哪个更适合后续诊断看PE就够了——PE低的占优。IMF1中心频率偏移到28.5Hz但包络谱有60Hz成分那是转频倍频的边带不是故障证据。5. 复现路上的高频踩坑点K值、alpha、收敛条件和真实信号的差距5.1 模态混叠和“成对模态”现象往往是中心频率重叠现象分解结果里有两个模态的中心频率几乎相同相差不到1Hz时域波形高度相似一个像是另一个的低通滤波版本。原因K值设置过大或者alpha设置过小导致带宽约束松弛相邻模态的频谱发生了重叠同一个频率成分被同时分配给了两个模态。解决思路先调K而不是先调alpha。把K减1后重跑看中心频率是否分离。如果减K后仍重叠再把alpha往大调每次乘2倍观察中心频率差的变化。文献里K3或4反而比K8效果更稳定这是复现时最容易接受的教训——模态数不是越多越好。5.2 高频故障成分被低频模态吞掉的根源频域初始化问题现象仿真信号里冲击谐振频率2000Hz附近分解后却找不到任何模态的中心频率靠近这个值。原因VMD对初始中心频率敏感初始化时频率均匀分布在[0, fs/2]上如果K较小或alpha较大迭代过程可能陷入局部最优——全部模态被低频大能量成分吸引高频带里没有模态锚点。解决把初始化方式从均匀分布改成含噪声的扰动分布。每次运行前给初始中心频率加一个高斯噪声多跑几组取效果最好的结果。更实用的方式是提高K值让频带划分更细避免高频段“真空”。这里就体现出为什么仿真信号在复现里如此重要你能确认故障频率的预期位置才能发现收敛到低频局部最优的问题。5.3 排列熵筛选翻车参数不匹配导致的假低熵现象某个纯噪声模态的PE算出来很低反而压过了真正的故障模态。原因数据段太短。PE需要的数据量随m指数增长如果模态时长只有0.2秒、采样率8000Hz有效点数1600个m取6以上时各排列模式统计次数不够概率估计的方差大得离谱。解决统一在相空间重构前做数据长度校验。经验标准是最少需要m! × 10个点。m5时需要1200个点以上比这短的数据段直接放弃PE筛选改用包络谱峰值因子。另一个隐藏限制是信号必须零均值VMD分解出的模态虽然理论上均值接近0但端点效应可能引入直流偏移先对每个模态做去均值处理再算PE。5.4 真实信号上的三处参数迁移失败采样率、转速波动、噪声非平稳在仿真信号上调好的参数直接迁移到真实数据上最常见的翻车点是三处。第一采样率不同。仿真里fs8000Hz、alpha2000迁移到50kHz采样的数据alpha需要按采样率的平方同比例放大否则每个模态的带宽约束完全失效。第二真实轴的转速有波动故障特征频率不是固定值而是围绕中心频率有±2%的随机漂移包络谱的谱峰会被抹宽。第三真实信号的噪声不是白噪声是带颜色的多为低频干扰随机尖峰。VMD遇到非平稳噪声时可能会专门分解出一个“噪声模态”中心频率落在低频区波形呈非周期抖动。我的做法是复现到这一步时增加一步计算所有模态的谱峭度谱峭度最高的模态往往对应冲击成分再对PE和谱峭度做加权排序。可以理解为给筛选机制上了双保险。5.5 迭代不收敛和运行时间异常tol、max_iter和信号长度的搭配现象程序跑了几百轮迭代仍没有收敛运行时间超出预期几倍。原因tol设置过小比如1e-9或者信号长度超过100万点时每个模态都做完整FFT和IFFT迭代成本过高。解决把收敛判据从模态差的绝对L2范数改成相对值——差除以该轮模态的能量这样不同幅值信号之间可以直接对比。max_iter不要设成无限常见做法是设500轮硬性截断记录最后两次迭代的差异。如果500轮后仍未收敛优先怀疑alpha太小导致模态更新震荡。工程上N200万的信号先降采样再分解是更务实的路线VMD本身的分辨率并不依赖全长数据。6. 验证你的复现是否成功构造已知故障信号用批量跑参找最优分量复现是否成功得看你的流程能不能从已知答案里找出答案。构造一个故障特征频率为87.3Hz、转频为13Hz的轴承信号把alpha从500跑到8000、K从3跑到7代码里加一层网格搜索以“包络谱在87.3Hz处信噪比峰值/邻域均值”作为打分标准自动选出最优参数组合param_grid [] for K in [3, 4, 5, 6, 7]: for alpha in [500, 1000, 2000, 4000, 8000]: param_grid.append((K, alpha)) best_score -np.inf best_params None for K, alpha in param_grid: u_hat_t, omega, _ VMD(signal_test, KK, alphaalpha) for k in range(K): _, spec envelope_spectrum(u_hat_t[k], fs) target_idx int(87.3 / (fs / 2 / (len(spec) * 2))) # 目标频率左右各3Hz区间的峰值为信号能量30-300Hz其余区间的均值为噪声底 band int(3 / (fs / 2 / (len(spec) * 2))) peak np.max(spec[target_idx-band:target_idxband]) noise_floor np.mean(spec[30:300]) score peak / noise_floor if score best_score: best_score score best_params (K, alpha, k) print(f最优参数: K{best_params[0]}, alpha{best_params[1]}, 最优模态{best_params[2]1})这段代码把“调参”从玄学变成了可量化的搜索。它在每个参数组合下运行完整VMD分解再对每个模态做包络谱评分最后输出最大化信噪比的组合。文献复现到这里才算闭环你不仅复现了方法还对方法的参数边界有了自己的实测数据。批量跑参时有两点需要说明。第一alpha的搜索范围不要跨数量级500到8000是工程上的常用区间越大越容易把所有模态压成窄带。第二评分函数只用包络谱信噪比还不够建议和排列熵组成联合判据——信噪比高但PE也高的模态大概率是强噪声凿出周期性假象直接排除。整个方案走下来我最大的心得是VMD本身不是玄学它的每个参数都对应明确的物理含义判断一部文献能否复现第一件事是先把信号模型在合成数据上搭出来否则你连正确结果长什么样都不知道。这套仿真分解筛选定位的流程轴承故障能用齿轮箱和电机电流信号同样能用改的只是故障特征频率的计算公式而已希望帮到你。本文还有配套的精品资源点击获取