企业官网建设流程全解析

简介：一套用于复现《基于VMD的故障特征信号提取方法》文献算法的MATLAB代码包，面向机械设备故障诊断与信号处理方向的研究者、工程师及研究生。代码包围绕VMD降噪与特征提取技术，将复杂非平稳信号分解为多个频率局部化的模态分量，帮助识别隐藏在噪声中的故障特征，进而判断设备运行状态。包内共有四个MATLAB脚本文件，压缩后体量仅约五千字节，结构完整清晰，包含VMD核心分解函数、主运行程序以及用于频谱分析和指标计算的辅助脚本，方便逐模块阅读、调试和二次开发。目前已有七百三十一人学习下载，适合具备信号处理基础与MATLAB编程能力、希望通过代码复现经典方法的读者。通过运行并研读这套代码，可以清晰理解VMD分解中的迭代优化与正则化策略，掌握从原始振动信号中提取特征分量的完整流程，为实际工程中的故障诊断与预测性维护提供可靠的技术参考。

1. 变分模态分解与故障特征信号提取：为什么照着文献复现 VMD 总在第一步就翻车

去年我拿到一篇学位论文，标题就是“基于 VMD 的故障特征信号提取方法”，参数表写得清清楚楚：模态数 K=6，惩罚因子 α=2000。我把这两个数字原样抄进代码，跑出来的结果却是一堆中心频率挤在一起的怪波形，包络谱里连一个像样的故障特征峰都找不到。问题不在论文，而在变分模态分解（VMD，Variational 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 必须等于或大于实际分量数，小于会让两个分量挤进同一个模态，大于会拆出虚假模态。alpha=2000 是大多数场景的起步值，它越大每个模态的带宽越窄，越小越宽。tau 一般保持 0，只有信号本身噪声极重时才调高。DC 在大多数振动信号场景设 0，除非你要专门分离趋势项。

2.3 怎么确认迭代真的收敛：画中心频率曲线比盯损失值更直观

很多人只取 omega 最后一列，但 VMD 默认最大迭代次数有限，如果没收敛，最后一列也是“算完”的数值，并不代表它是对的。我的习惯是把中心频率随迭代次数的曲线画出来，一眼判断收敛质量：

import matplotlib.pyplot as plt for k in range(K): plt.plot(omega[k, :], label=f'IMF{k+1}') 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 值别照抄文献：先跑一遍中心频率扫描，看“打架”的位置

文献里写 K=6，那是它自己的信号成分决定的。你换轴承型号、换转速、换测点，采集到的信号分量数完全不一样。把文献的 K 当真理，是我见过最多人犯的错。正确做法是从 K=2 开始往上扫，打印每个 K 下收敛后的中心频率，观察相邻值之间的距离：

for K in range(2, 8): u, u_hat, omega = VMD(sig, alpha=2000, tau=0, K=K, DC=0, init=1, tol=1e-7) freqs_hz = np.sort(omega[:, -1] * fs / (2*np.pi)) print(f'K={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 的顺序，我一般会先固定 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, m=3, delay=1): 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, alpha=2000, tau=0, K=K, DC=0, init=1, tol=1e-7) pes = [permutation_entropy(u[k]) for k in range(K)] print(f'K={K}, 平均PE={np.mean(pes):.4f}, 各分量={np.array(pes).round(4)}')

逻辑说明：PE 计算的是每个模态的排列熵，m=3 是最常用的嵌入维度，delay=1 表示用连续采样点。K 从 2 升到 3 时，平均 PE 通常显著下降；再继续增加 K，平均 PE 趋于平稳或回升，新增模态的单个 PE 值会明显高于老模态。取平均 PE 最低或回落点的 K 作为最终模态数。

参数说明：m 不要取太大，序列长度有限时 m 太大统计不充分，一般 m=3 或 4。排列熵对幅值不敏感，这既是优点也是缺点——它只关心序列的序关系，两个幅值差异很大但形状相似的信号会得到相近的 PE 值。

4. VMD 复现避坑：模态混叠、端点毛刺与参数玄学的 6 个现场

4.1 现象：K 设 6 出来 4 个模态中心频率挤在 96~105Hz，包络谱出现三条兄弟峰

照抄文献 K=6 时，分解结果里两个或三个模态的中心频率会挤在一个很窄的频带内，包络谱在这个频带附近出现多个峰值，很难判断哪个是真正的故障特征频率。

原因：K 超过信号真实分量数，过分解导致同一个频带被多个模态瓜分。ADMM 迭代时，能量会随机分配到相邻模态里，谁先收敛谁拿大头，结果不可预测。

解决：用上一章的中心频率扫描脚本，从 K=2 开始逐个看，找到中心频率开始“打架”的那个 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_len+1][::-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 的 init=1 表示中心频率均匀初始化，init=0 表示全部从 0 开始，这两个选项都试一次，取收敛更平滑的那个。很多博客会把这两个初始化写反，遇到结果异常时两个都跑一遍最稳妥。

4.5 现象：打印出来的中心频率是 0.157、0.377，和文献里的 50Hz、120Hz 对不上

文献里给的是 Hz，vmdpy 输出的是归一化角频率，两者差一个换算系数。第一次复现的人经常拿着这个数字怀疑自己代码写错。

原因：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 先筛掉“废模态”：峭度 + 相关系数双指标选 IMF

VMD 输出 K 个模态不是每个都有诊断价值。实测信号里常有两个模态是噪声或趋势项。我的筛选习惯是用一个简单循环给所有模态打分：

from scipy.stats import kurtosis for k in range(K): r = np.corrcoef(sig, u[k])[0, 1] kurt = kurtosis(u[k], fisher=False) # 经典峭度定义，正态分布为3 print(f'IMF{k+1}: 相关系数r={r:.3f}, 峭度={kurt:.2f}')

逻辑说明：相关系数衡量模态与原信号的线性相关程度，能量占比高的模态 r 自然大；峭度衡量波形冲击性，故障冲击信号峭度远大于 3，白噪声峭度接近 3。两个指标一起看：r 高但峭度低，说明这是大能量正弦成分；峭度高但 r 低，可能是噪声尖峰；两者都高的才优先进入包络谱分析。

参数说明：具体的筛选阈值要参考背景噪声水平。噪声低时 r>0.3 且峭度>3 是比较稳的组合；噪声大时把 r 阈值放宽到 0.1，但峭度阈值不要放宽，否则噪声模态混进来。

5.2 包络谱对齐：先算故障特征频率，再在谱图里找峰，别靠肉眼“看着像”

选好 IMF 之后做包络谱分析。包络谱的原理是：故障冲击会通过高频共振幅值调制表现出来，用希尔伯特变换取出包络，再对包络做 FFT，故障特征频率就会在低频段形成谱峰。

def bearing_freqs(fr, n_balls, d, D, contact_angle=0): 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, prominence=0.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、bpfo+fr。用 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:i+len(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 类方法，第一步永远是构造一个已知特征的合成信号跑通整条链路，再换真实数据。这个习惯帮我省掉了无数次无效调参。希望帮到你。

本文还有配套的精品资源，点击获取