
简介这份资源围绕《基于VMD的故障特征信号提取方法》文献展开复现面向具备一定信号处理基础与MATLAB编程能力的读者尤其是从事机械设备故障诊断、非平稳信号降噪与特征提取方向的学习者和研究人员。资源包共4个文件均为m脚本压缩包约5KB其中核心算法、主程序流程与辅助分析函数分工明确便于按模块阅读与调试。已有731人学习下载说明其在VMD入门与文献复现场景中具有一定参考价值。通过运行主程序读者可完整走通信号分解、模态分量获取、频谱与时频可视化以及特征提取等环节直观观察VMD从噪声中分离故障特征的效果辅助函数则补充了功率谱密度计算、峰值检测等细节分析能力。整体代码结构紧凑适合作为理解VMD迭代优化与正则化流程的实践起点也可在此基础上修改参数、替换信号进一步验证和改进自己的信号处理方案。1. 从一段轴承振动信号说起VMD 故障特征提取到底在做什么手里有一段轴承外圈故障的振动信号采样率 12 kHz转频 30 Hz理论上外圈故障特征频率大概在 107 Hz 附近。你把它丢进 FFT频谱上确实能看到一堆峰但低频段被转频和它的倍频糊住高频段又混着共振调制故障特征频率那条线要么被淹没要么藏在边频带里根本认不出来。这时候很多论文会告诉你上 VMD把信号分解成若干模态故障特征就出来了。但真到自己复现《基于VMD的故障特征信号提取方法》这类文献时问题立刻变成三连K 取几alpha 取多少分解完怎么判断哪个模态才是故障模态VMD变分模态分解本质上是把信号分解问题写成一个带约束的变分问题每个模态被约束在中心频率附近通过交替方向乘子法迭代求解。它和 EMD 最大的区别是EMD 靠递归筛分模态混叠和端点效应是玄学VMD 靠频域迭代K 和 alpha 一旦定下来结果是确定的、可复现的。这也是为什么故障诊断领域这几年大量文献都往 VMD 上靠——它给了你一个可以调参、可以对比、可以写进论文的确定性框架。这篇笔记面向的是要真正把这类文献跑通的人手里有振动数据CWRU、XJTU-SY、自采都行想复现出「分解—选模态—包络解调—看到故障特征频率」这条完整链路并且想知道参数怎么设、坑在哪、结果不对时该看什么。下面按我实际复现的顺序讲代码用 Python核心库是 vmdpy 和 scipy。2. VMD 的数学骨架与复现前的环境准备2.1 变分问题的构造为什么它比 EMD 更「可控」VMD 的核心思想是把一个实信号分解成 K 个本征模态函数IMF每个 IMF 被定义为一个调幅调频信号围绕各自的中心频率。构造的变分问题是最小化所有模态的解析信号带宽之和约束条件是所有模态加起来等于原信号。为了求解引入二次惩罚因子 alpha 和拉格朗日乘子把约束问题转成无约束问题再用 ADMM 交替更新模态、中心频率和乘子。这里有两个参数直接决定结果K 是模态数alpha 是带宽惩罚。alpha 越大每个模态的带宽越窄模态之间越不容易混alpha 越小带宽越宽容易把多个频率成分塞进一个模态。K 决定了你能分出几个成分K 太小会把故障特征和转频混在一起K 太大则会产生虚假模态把噪声也拆成「模态」。理解这一点后面调参就有方向了K 和 alpha 不是拍脑袋而是由信号里实际有几个窄带成分决定的。2.2 环境搭建与最小可运行示例先装依赖。vmdpy 是 Python 里比较常用的 VMD 实现接口简单适合复现文献。pip install numpy scipy matplotlib vmdpy下面是一个最小可运行示例用一段合成的多分量信号验证 VMD 能不能把成分分开。import numpy as np import matplotlib.pyplot as plt from vmdpy import VMD # 构造合成信号三个分量 噪声 fs 12000 # 采样率 T 1.0 # 时长 t np.arange(0, T, 1/fs) f1, f2, f3 50, 300, 1200 # 三个中心频率 sig (np.cos(2*np.pi*f1*t) 0.6*np.cos(2*np.pi*f2*t) 0.3*np.cos(2*np.pi*f3*t) 0.1*np.random.randn(len(t))) # VMD 参数 alpha 2000 # 带宽惩罚 tau 0 # 噪声容忍0 表示严格保数据 K 3 # 模态数 DC 0 # 不含直流 init 1 # 中心频率初始化方式 tol 1e-7 # 收敛容差 u, u_hat, omega VMD(sig, alpha, tau, K, DC, init, tol) # u 形状为 (K, N)每行是一个模态 print(分解得到模态数:, u.shape[0]) print(各模态中心频率(Hz):, omega[-1] * fs)逻辑说明VMD 返回的 u 是 K 个模态的时域波形omega 是每次迭代的中心频率取最后一行就是收敛后的中心频率。参数上alpha2000 是文献里最常用的起点K3 对应我构造的三个分量。tau0 表示不允许噪声泄漏到模态外适合干净信号如果信号噪声大可以设 tau 为 0.1 到 0.3 之间让算法对噪声更宽容。跑完你会看到三个模态的中心频率大致落在 50、300、1200 Hz 附近说明分解有效。这一步跑通后面换成真实轴承信号只是换数据源。2.3 用 CWRU 数据替换合成信号复现文献时数据源通常是凯斯西储大学CWRU轴承数据集。假设你已经下载了 12k Drive End 的 mat 文件读取方式如下。from scipy.io import loadmat # CWRU 数据文件变量名通常是 X105_DE_time 这类 mat loadmat(105.mat) key [k for k in mat.keys() if DE_time in k][0] x mat[key].flatten() # 取一段做分析避免整段太长 x x[:12000] # 1 秒数据 x x - np.mean(x) # 去直流逻辑说明CWRU 的 mat 文件里变量名带 DE_time 的是驱动端加速度信号取 1 秒长度足够做一次分解。去直流是必须的否则 DC 分量会占据一个模态干扰判断。参数上采样率 12 kHz 对应 CWRU 的 12k Drive End 数据如果你用的是 48k 数据采样率要改成 48000后面算频率时同步改。3. 参数怎么定K 与 alpha 的选法及分解结果判读3.1 K 的确定从中心频率重复到峭度准则K 是 VMD 里最敏感的参数。K 太小故障特征频率和转频会被塞进同一个模态K 太大会出现中心频率相近的虚假模态。文献里常见的做法有两种观察中心频率法和峭度准则法。观察中心频率法从小到大试 K看最后一次迭代的中心频率。如果两个模态的中心频率非常接近比如相差不到 5%说明 K 取大了应该减小。这个方法直观但需要人工判断。峭度准则法峭度对冲击成分敏感故障冲击会让模态的峭度升高。计算每个模态的峭度选峭度最大的模态作为故障模态同时用峭度随 K 的变化判断 K 是否合适。下面是一个批量试 K 的脚本。from scipy.stats import kurtosis def try_k(x, fs, k_list, alpha2000): results {} for K in k_list: u, _, omega VMD(x, alpha, 0, K, 0, 1, 1e-7) freqs omega[-1] * fs kur [kurtosis(u[i]) for i in range(K)] results[K] {freqs: freqs, kurt: kur} print(fK{K}, 中心频率{np.round(freqs,1)}, 峭度{np.round(kur,2)}) return results res try_k(x, 12000, [3,4,5,6,7])逻辑说明这段脚本对每个 K 跑一次 VMD输出中心频率和每个模态的峭度。判断标准是如果某个 K 下出现两个中心频率几乎相同的模态说明 K 偏大如果峭度最大的模态在 K 增大时不再明显变化说明 K 已经够了。参数上alpha 先固定 2000等 K 定了再调 alpha。3.2 alpha 的调整带宽与模态混叠的权衡alpha 控制带宽。alpha 太小模态带宽宽容易混叠alpha 太大带宽窄可能把同一个故障特征拆到两个模态里。经验范围是 1000 到 5000文献里 2000 出现频率最高。调整方法是固定 K让 alpha 从 500 到 5000 变化看故障模态的包络谱里特征频率是否清晰。def try_alpha(x, fs, K, alpha_list): for alpha in alpha_list: u, _, omega VMD(x, alpha, 0, K, 0, 1, 1e-7) # 选峭度最大的模态 kur [kurtosis(u[i]) for i in range(K)] idx int(np.argmax(kur)) # 包络谱 env np.abs(scipy.signal.hilbert(u[idx])) env_spec np.abs(np.fft.rfft(env - np.mean(env))) freqs np.fft.rfftfreq(len(env), 1/fs) # 找包络谱峰值 peak_idx np.argsort(env_spec)[-5:] print(falpha{alpha}, 故障模态{idx}, 包络谱峰值频率{np.round(freqs[peak_idx],1)}) import scipy.signal try_alpha(x, 12000, 5, [500, 1000, 2000, 3000, 5000])逻辑说明对每个 alpha选峭度最大的模态做希尔伯特包络再对包络做 FFT 得到包络谱。故障特征频率会在包络谱上出现峰值。参数上如果包络谱峰值频率接近理论故障特征频率比如外圈 107 Hz说明 alpha 合适如果峰值杂乱或偏移说明 alpha 需要调整。3.3 故障模态的判读包络谱与理论频率对照分解完、选完模态最后一步是验证。以 CWRU 外圈故障为例理论故障特征频率计算公式是BPFO (n/2) * fr * (1 - (d/D)*cos(theta))其中 n 是滚珠数fr 是转频d 是滚珠直径D 是节径theta 是接触角。CWRU 的 6205 轴承参数n9d7.94mmD39.04mmtheta0。转频 30 Hz 时BPFO 约 107 Hz。在包络谱上你应该看到 107 Hz 及其倍频214、321 Hz的峰值。如果看到的是转频 30 Hz 的峰值说明选的模态是转频模态不是故障模态需要重新选模态或调整 K。提示包络谱的横轴范围建议限制在 0 到 500 Hz故障特征频率通常在这个范围内高频段是噪声和共振看了反而干扰判断。4. 复现时最容易翻车的五个地方4.1 现象分解出的模态全是噪声包络谱没有明显峰值原因信号没有去直流或者 alpha 太小导致模态带宽过宽把噪声也当成模态分解出来。另外如果原始信号里故障冲击本身很弱VMD 可能把冲击分散到多个模态里。解决先做去直流和简单带通滤波比如 500 到 3000 Hz把低频转频和高频噪声先压一压。alpha 从 2000 起步不要一上来就设 500。如果冲击弱可以先用包络谱确认原始信号里有没有故障特征再决定要不要 VMD。4.2 现象K 增大时中心频率出现重复但峭度最大的模态峭度反而下降原因K 过大VMD 把噪声拆成了虚假模态这些模态的中心频率和真实模态接近导致真实模态的能量被分散峭度下降。解决以中心频率不重复为第一准则峭度作为辅助。如果 K5 时中心频率开始重复就取 K4。不要盲目追求大 K文献里 K 取 3 到 6 最常见。4.3 现象包络谱峰值频率和理论故障特征频率对不上差了几赫兹原因转频估计不准。CWRU 数据里转频不一定是 30 Hz不同负载下转频会变。另外包络谱的频率分辨率是 fs/NN 是信号长度如果 N 太小分辨率不够峰值会偏移。解决先算准转频。CWRU 数据里可以用转速计信号或者从频谱里找转频峰值。包络谱分析时信号长度至少取 1 秒分辨率 1 Hz 以内。如果还是对不上检查轴承参数是否用错不同轴承型号参数不同。4.4 现象VMD 运行很慢或者报内存错误原因信号太长。VMD 的计算量和信号长度、K 都相关直接对 10 秒数据跑 K10 会非常慢。解决分段分析每段 0.5 到 1 秒分别做 VMD 和包络谱最后看特征频率是否稳定。另外vmdpy 的 tol 不要设太小1e-7 足够设 1e-10 会多迭代很多次。4.5 现象换一组数据后之前调好的 K 和 alpha 完全失效原因VMD 参数对信号敏感不同故障类型、不同负载、不同采样率下最优参数会变。文献里给的参数只针对特定数据。解决把 K 和 alpha 的搜索做成自动化流程每次换数据先跑一遍 K 搜索和 alpha 搜索用中心频率重复和包络谱峰值作为判据。不要指望一套参数打天下。5. 进阶把 VMD 和包络谱串成一条自动流水线复现文献的终点不是跑出一次结果而是把「分解—选模态—包络解调—频率对照」做成可重复的流程。我现在的习惯是写一个函数输入信号和轴承参数输出故障特征频率和包络谱图中间所有参数自动搜索。def vmd_fault_extract(x, fs, bpfo_theory, k_rangerange(3,8), alpha_range[1000,2000,3000]): best None for K in k_range: for alpha in alpha_range: u, _, omega VMD(x, alpha, 0, K, 0, 1, 1e-7) freqs omega[-1] * fs # 中心频率重复检查 if len(freqs) 1 and np.min(np.diff(np.sort(freqs))) 0.05 * fs / K: continue kur [kurtosis(u[i]) for i in range(K)] idx int(np.argmax(kur)) env np.abs(scipy.signal.hilbert(u[idx])) env_spec np.abs(np.fft.rfft(env - np.mean(env))) spec_freqs np.fft.rfftfreq(len(env), 1/fs) # 找包络谱里最接近理论故障频率的峰值 peak_idx np.argmax(env_spec) peak_freq spec_freqs[peak_idx] err abs(peak_freq - bpfo_theory) if best is None or err best[err]: best {K: K, alpha: alpha, peak_freq: peak_freq, err: err, u: u, idx: idx} return best best vmd_fault_extract(x, 12000, 107) print(f最优 K{best[K]}, alpha{best[alpha]}, 包络谱峰值{best[peak_freq]:.1f} Hz)逻辑说明这个函数遍历 K 和 alpha 的组合用中心频率重复做初筛用包络谱峰值和理论故障频率的误差做最终判据。参数上k_range 和 alpha_range 根据你的数据调整bpfo_theory 用轴承参数算出来。返回的 best 里包含最优参数和对应模态可以直接画图。验证方法拿 CWRU 的 105.mat外圈故障跑一遍看最优参数下包络谱峰值是否在 107 Hz 附近。再拿 130.mat内圈故障跑理论内圈故障频率约 162 Hz看是否对得上。如果两组数据都能对上说明流水线是可靠的。一个具体技巧包络谱画图时把理论故障频率及其倍频用竖线标出来一眼就能看出峰值是否对齐。这个习惯帮我省了很多反复对照的时间。注意自动搜索会跑很多次 VMD如果数据长、K 范围大耗时会明显增加。建议先用短数据0.5 秒粗搜再用长数据1 秒以上精搜。我踩过最深的一个坑是早期迷信文献里给的 K5、alpha2000换了一组数据后包络谱死活对不上折腾了一整天才发现是转频估错了导致理论故障频率算错后面所有判断都跟着错。从那以后我每次先花十分钟把转频和轴承参数确认清楚再动 VMD。希望帮到你。本文还有配套的精品资源点击获取