简介资源为基于粒子群优化算法PSO改进变分模态分解VMD的MATLAB实现方案面向信号处理、故障诊断及机器学习方向的研究者与工程师用于对非平稳、非线性信号进行自适应分解与特征提取。包内以4个m文件为主包含核心VMD算法函数、PSO优化调用及仿真脚本另含2个txt辅助说明文档和1个url链接便于理解参数设置与算法流程压缩包整体约2.29MB。已有221人学习下载。通过该代码可快速掌握PSO-VMD的耦合逻辑借助PSO全局搜索能力自动确定VMD惩罚因子与模态数有效避免人工调参的盲目性提升分解精度与运算效率仿真改动文件和辅助函数也提供了较好的二次开发基础适合需要对比不同分解方法、开展故障特征提取或学术研究的读者参考使用。1. 变分模态分解最难调的 K 和 alpha交给粒子群算法去定做旋转机械故障诊断的人大概率在变分模态分解VMD上卡过壳分解本身不复杂复杂的是 K模态数量和 alpha惩罚因子到底怎么定。K 设小了故障冲击和转频混在一个模态里K 设大了又会出现空模态或中心频率扎堆。alpha 更是直接影响频带宽度和 K 耦合在一起根本没法单独调。人工试凑一轮就要完整跑一次 VMD还要逐个看频谱耗时且没有一致性。把粒子群算法PSO套在 VMD 外面用包络熵最小作为寻优目标让粒子自己去搜索 K 和 alpha就是 pso-vmd 这类算法最主流的落地方式。这篇文章覆盖粒子群算法原理、适应度函数设计、完整可运行代码和工程参数边界适合信号处理方向的算法工程师和刚开始接触 VMD 参数优化的研究者。2. 变分模态分解为何对 K 和 alpha 敏感以及粒子群算法的适配逻辑2.1 变分模态分解的中心频率分配与惩罚因子耦合变分模态分解把信号看成 K 个调幅调频模态的叠加每个模态围绕自己的中心频率omega_k在频域内占据一段带宽。它的数学模型里有两个互相拉扯的目标一是所有模态的总带宽尽量窄二是 K 个模态加起来要能复现原信号。alpha 正是用来平衡这两者的权重alpha 偏大时模型把更多注意力放在重构精确上频带会变宽两个模态的中心频率一旦靠近就容易粘在一起这就是模态混叠alpha 偏小时频带被压得很窄重构误差变大部分真实分量会被削掉或拆散。K 和 alpha 从来不是独立变量。K 决定频谱要切几刀alpha 决定每一刀允许切多宽。同一个信号在 alpha2000 时 K4 可能很干净换成 alpha8000 之后 K4 可能就混叠了只能往上加到 K5 或 K6 才能重新分开。这就是为什么网格搜索很难奏效你要在一个二维平面上同时找两个互相牵制的参数网格步长稍微放大就会漏掉那块好参数区。一个更隐蔽的问题是中心频率的初始化。VMD 里init1表示中心频率均匀分布在频带上init0表示初始化为零。参数不合适时VMD 的迭代经常收敛到错误的局部频带也就是某个模态中心频率漂移到明显没有能量的位置。这种失败不是 VMD 本身的问题而是 (K, alpha) 的组合不合理。2.2 粒子群算法原理连续速度更新与离散 K 的兼容处理粒子群算法原理很适合这种场景它不依赖目标函数梯度只需要能算出这组参数比那组参数好就可以把解空间搜索起来。每个粒子是一个二维位置向量第 1 维是log10(alpha)第 2 维是 K。为什么 alpha 要对数化因为 alpha 在 500 到 10000 之间跨了两个数量级线性步长在小区间浪费精度在大区间跳得太快对数域搜索更接近人类的调参直觉。标准的速度-位置更新可以简写成下面这段# 第 i 个粒子的速度与位置更新 v[i] w * v[i] c1 * r1 * (pbest[i] - x[i]) c2 * r2 * (gbest - x[i]) x[i] x[i] v[i]其中 pbest 是单个粒子历史最优位置gbest 是整个种群的历史最优。惯性权重 w 控制粒子继承上一轮速度的程度前期希望它大一点去探索全局后期希望它小一点来精细搜索所以常见做法是让 w 从 0.9 线性衰减到 0.4。c1 和 c2 分别是自我认知和社会认知系数一般取 1.5 到 2.0。K 是整数而粒子的位移更新是连续的。常用的兼容处理是把更新后的 K 直接round()成整数再用边界语句拉回 [2, 10] 区间。PSO 不需要梯度信息所以离散化不会带来数值问题只会让粒子在整数点之间跳跃这在实际调参里完全够用。2.3 网格搜索、遗传算法与粒子群算法的代价对比VMD 本身是迭代算法单次分解耗时从几十毫秒到几秒不等参数寻优的核心指标是在尽量少的 VMD 调用次数内找到可用的 (K, alpha)。下面这张表是我在实际对比中的感受寻优方法评估次数量级整数变量处理实现成本主要问题网格搜索9×Malpha 分 M 档天然支持低档位稍密就爆炸漏掉非均匀分布的优区遗传算法种群×代数需要专门设计编码和解码中交叉变异算子对连续维的搜索效率偏低粒子群算法20×20 到 40×30速度更新后直接取整低早熟收敛风险随机种子敏感网格搜索最大的坑在于 alpha 档位怎么取。如果线性取 20 个点1000 到 10000 之间步长是 450精细结构会整个漏掉如果取 50 个点又要跑 450 次 VMD。遗传算法的交叉操作在离散 K 维上容易破坏掉已经找到的好组合。粒子群算法在二维问题上收敛很快20 个粒子跑 15 到 25 代通常就能找到可用的解代码量也最小。评估次数虽然也在几百次量级但因为每个粒子都保留自己的 pbest即使某一代全体方向带偏后期也能拉回来。3. 用 Python 把 pso-vmd 写成能直接跑的代码3.1 包络熵适应度函数与 vmdpy 的返回值pso-vmd 的代码组织通常是三层最外层是 PSO 主循环中间层是包络熵适应度函数最底层是 VMD 分解调用。VMD 不需要自己实现常见的 Python 封装是 vmdpy 库安装一行命令pip install vmdpy numpy scipyvmdpy 的完整调用参数是VMD(f, alpha, tau, K, DC, init, tol)返回值是u, u_hat, omega。其中u是分解出的模态数组shape 为(K, 信号长度)omega是每次迭代的中心频率取最后一行就是收敛后的中心频率。tau 是噪声容忍度通常设 0DC 设为 0 表示不单独提取直流分量init 设为 1 表示中心频率均匀初始化tol 用默认的 1e-7。建议拿到u先打印u.shape确认维度顺序不同封装的返回顺序偶尔有差异。适应度函数我用的是最小包络熵对每个 IMF 做 Hilbert 变换取包络归一化后计算信息熵。包络越小、冲击越集中熵值越小说明这个模态提取到的是干净的故障冲击而不是噪声。具体代码import numpy as np from scipy.signal import hilbert from vmdpy import VMD def envelope_entropy(imf): # 计算单个模态的包络熵值越小代表冲击特征越突出 analytic hilbert(imf) env np.abs(analytic) env env / (env.sum() 1e-12) return -np.sum(env * np.log(env 1e-12)) def cost_function(signal, alpha, K): # 调用 vmdpy 完成分解u 的形状是 (K, N) u, _, _ VMD(signal, alpha, 0, int(K), 0, 1, 1e-7) entropies [envelope_entropy(u[i, :]) for i in range(u.shape[0])] # 取最小包络熵作为适应度适合故障诊断场景 return min(entropies)取min而不是取平均是因为周期冲击通常集中在某一个模态上最小包络熵对应的就是那个最值得保留的分量。如果做的是信号去噪解法不同后面第 5 章会单独说。3.2 pso-vmd 主循环速度更新、整数化与边界处理下面是一个完整的 PSO-VMD 优化类可以直接放进脚本跑class PSO_VMD: def __init__(self, signal, pop_size20, max_iter20): self.signal signal self.pop_size pop_size self.max_iter max_iter # alpha 在对数域 [1000, 5000] 内搜索K 在 [2, 10] 内搜索 self.lb np.array([np.log10(1000), 2]) self.ub np.array([np.log10(5000), 10]) def cost_func(self, alpha, K): return cost_function(self.signal, alpha, K) def optimize(self): dim 2 w_max, w_min, c1, c2 0.9, 0.4, 1.5, 1.5 x np.zeros((self.pop_size, dim)) x[:, 0] np.random.uniform(self.lb[0], self.ub[0], self.pop_size) x[:, 1] np.random.randint(2, 11, self.pop_size) v np.zeros_like(x) pbest x.copy() pbest_cost np.full(self.pop_size, np.inf) gbest x[0].copy() gbest_cost np.inf trace [] for it in range(self.max_iter): w w_max - (w_max - w_min) * it / (self.max_iter - 1) r1 np.random.random((self.pop_size, dim)) r2 np.random.random((self.pop_size, dim)) v w * v c1 * r1 * (pbest - x) c2 * r2 * (gbest - x) x x v for i in range(self.pop_size): # 边界吸收越界位置拉回边界速度清零避免堆积 for d in range(dim): if x[i, d] self.lb[d]: x[i, d], v[i, d] self.lb[d], 0 if x[i, d] self.ub[d]: x[i, d], v[i, d] self.ub[d], 0 # K 必须是整数取整后保持在合法区间 x[i, 1] round(x[i, 1]) if x[i, 1] 2: x[i, 1] 2 if x[i, 1] 10: x[i, 1] 10 for i in range(self.pop_size): alpha 10 ** x[i, 0] cost self.cost_func(alpha, int(x[i, 1])) if cost pbest_cost[i]: pbest_cost[i] cost pbest[i] x[i].copy() if cost gbest_cost: gbest_cost cost gbest x[i].copy() trace.append(gbest_cost) print(fiter {it1}, cost{gbest_cost:.4f}, falpha{10**gbest[0]:.1f}, K{int(gbest[1])}) return 10 ** gbest[0], int(gbest[1]), gbest_cost, trace这段代码的关键点有三个。速度裁剪通过v 0来防振荡比限制速度幅值更直接避免粒子在边界频繁来回穿越。K 的取整放在速度更新之后保证 pbest 和 gbest 里存的都是合法整数。trace数组记录每一代全局最优适应度后面画收敛曲线用。3.3 pso-vmd 中三个影响寻优效果的关键参数粒子数量、迭代次数和搜索边界是 pso-vmd 里最值得调的三个参数其他参数保持经验值即可。参数常见取值对结果的影响pop_size20 到 40小于 15 容易早熟大于 50 评估次数翻倍max_iter15 到 30VMD 耗时高时优先保小迭代次数alpha 搜索范围log10 域 [1000, 5000]低于 1000 保真度差高于 5000 频带过宽K 搜索范围[2, 10]上限太高会出现空模态上限太低会漏频带w 衰减区间0.9 到 0.4前期全局探索后期局部细搜c1, c21.5 到 2.0两者相等时收敛最平稳alpha 搜索范围需要根据信号调整。采样率低、频带窄的信号适合把上限压到 3000高频成分多的信号可以把下限提到 2000。K 的上界强烈建议不超过 10工程上绝大多数信号分解到 8 个模态已经很难解释了K 搜到 10 还嫌不够时先检查信号里是不是混着直流偏移和趋势项。3.4 代码跑起来之后怎么判断收敛正常运行会看到每代打印的cost逐渐下降alpha 和 K 在前几代大幅跳动后面逐渐收敛到一个小范围。典型输出类似iter 1, cost2.3154, alpha4230.2, K7 iter 6, cost2.0418, alpha3020.7, K4 iter 12, cost1.9933, alpha2715.5, K5 iter 20, cost1.9881, alpha2560.1, K5判断收敛的标准连续 5 代cost变化小于 1%且 K 不再频繁切换。如果 alpha 一直在边界上跳多半是搜索范围给窄了如果 K 顶到上界不动说明信号本身不适合用 VMD 拆或者适应度函数选错了方向。注意vmdpy 的默认返回值顺序在一些版本里存在差异第一次调用时先打印 u.shape确认是 (K, N) 而不是 (N, K)再继续写后面的计算。4. pso-vmd 实战从仿真轴承信号中提取故障特征频率4.1 构造带故障冲击的仿真信号滚动轴承外圈故障的典型特征是周期性冲击每个冲击会激发系统固有共振并快速衰减。仿真信号可以用转频正弦分量叠加衰减冲击串再加噪声来构造import numpy as np fs 12000 # 采样率 12 kHz t np.arange(0.0, 1.0, 1.0 / fs) fr 30.0 # 转频 30 Hz bpfo 107.0 # 外圈故障特征频率约 107 Hz # 转频分量 signal 0.5 * np.sin(2 * np.pi * fr * t) # 周期性衰减冲击共振频率取 800 Hz interval int(fs / bpfo) for start in range(0, len(t), interval): idx np.arange(0, min(200, len(t) - start)) decay np.exp(-0.05 * idx) impact np.sin(2 * np.pi * 800 * idx / fs) * decay signal[start:start len(idx)] 1.5 * impact # 叠加小幅白噪声 signal 0.25 * np.random.randn(len(t))interval由采样率除以特征频率得到约 112 个采样点一个冲击。衰减系数 0.05 控制在 200 个采样点内冲击基本衰减完这与真实轴承外圈故障的衰减特性量级一致。噪声标准差取 0.25是为了让冲击仍然肉眼可见又不至于让分解太轻松。故障诊断里这种仿真信号常用来验证算法参数是否有区分度。4.2 运行 pso-vmd 并解读最优 alpha、K 与模态直接用上一章的 PSO_VMD 类跑opt PSO_VMD(signal, pop_size20, max_iter20) alpha_best, K_best, cost_best, trace opt.optimize() # 用最优参数重新分解保留中心频率结果 u, _, omega VMD(signal, alpha_best, 0, K_best, 0, 1, 1e-7) # 打印中心频率omega 最后一行是收敛值 freqs omega[-1, :] * fs / 2 for idx in range(K_best): print(fIMF{idx1}: center_freq{freqs[idx]:.2f} Hz) # 对每个模态做频谱找到最大谱峰位置 spectrum np.abs(np.fft.fft(u[1]))[:len(t) // 2] f_axis np.linspace(0, fs / 2, len(spectrum)) print(fIMF2 peak_freq{f_axis[np.argmax(spectrum)]:.2f} Hz)固定随机种子为 42 时一次典型运行的输出通常是 K4 或 K5alpha 落在 2000 到 3000 区间。K4 时一般会看到IMF1 是转频 30 HzIMF2 或 IMF3 是 800 Hz 附近的共振带107 Hz 的故障调制频率会出现在某个模态的包络谱上。如果输出的 K 总是顶到 10先检查信号是否混入直流或者趋势项。alpha 的通知意义在于它决定了共振带的宽度。alpha 偏低时800 Hz 共振峰会分裂成多个小峰分布在相邻模态alpha 偏高时转频和冲击串会被合到同一个模态。pso-vmd 给出的 alpha 在 2000 到 3000 区间对应的是共振带恰好完整放入一个模态的分界点这与人工调参时看到的临界位置基本一致。4.3 从中心频率分布识别模态混叠和空模态VMD 跑完不能只盯着包络熵看还要看中心频率的分布是不是合理。中心频率两两如果靠得太近说明两个模态在争抢同一段频带如果某个 IMF 能量占比不足 1%说明出现了空模态。这两类问题用一行代码就能检查diff_freqs np.diff(freqs) print(相邻中心频率差值, diff_freqs) # 计算每个模态的能量占比 energy_ratio np.array([np.sum(u[i] ** 2) for i in range(K_best)]) energy_ratio / energy_ratio.sum() print(模态能量占比, energy_ratio)检查项判定方法处理思路模态混叠相邻中心频率差小于 2 倍频率分辨率增大 alpha或者减少 K 后再试空模态能量占比小于 1%说明 K 偏大建议降 K中心频率漂移某个中心频率落在无明显频谱峰的位置调整 alpha 后重新运行频率分辨率在这个例子里是 1 Hz采样率除以信号长度所以相邻中心频率差小于 3 Hz 时就要警惕混叠。这些检查应当在 pso-vmd 每次寻优结束后作为固定校验步骤不能因为适应度低就直接采信。5. 让 pso-vmd 在工程里更稳的四个技巧5.1 把适应度从包络熵换成样本熵包络熵适合周期性冲击明显的信号也就是故障诊断的主场。如果待分解信号是非周期、非平稳的比如地震波或风电功率序列包络熵容易把噪声脉冲误判成有效成分。这时可以把 cost_function 里返回改为样本熵计算每个模态的样本熵后仍然取最小。改动只在 cost_function 内部PSO 主循环不用动因为 PSO 只认输入输出的数值大小。5.2 多随机种子跑三次取最优粒子群算法是随机优化一次运行的 gbest 不一定代表真实最优。工程上常见操作是固定三个随机种子各跑一次取 cost 最小的一组参数。代价是寻优时间变成原来的三倍但换来的是参数可复现。如果三次跑出来的 alpha 和 K 方差很大说明搜索空间设计可能有问题需要回头检查边界设置。5.3 搜到 K 顶在上界时先查信号质量K 顶到上界不是小事。先检查信号是否做了去均值有没有线性趋势有没有明显工频干扰。VMD 对直流和趋势项非常敏感这两类成分会独自占掉一个模态把 K 的有效容量吃掉。把趋势项用多项式拟合去掉后重新跑K 通常会从 10 回落到 5 附近。5.4 长信号先降采样再寻优VMD 每次迭代都要做 FFT信号长度翻倍单次分解时间也近似翻倍。长度超过 10 万点的信号直接做参数寻优很吃亏。通常先降采样到 2 kHz 左右完成寻优再用最优 K 和 alpha 对原始采样率数据做一次分解。注意降采样前加抗混叠滤波器alpha 与采样率强相关直接跨采样率套用 alpha 不一定成立按带宽比例折算后通常更接近最优值。按照带宽比例折算 alpha是长序列工程信号最省时间的处理路径。本文还有配套的精品资源点击获取