简介本资源面向计算机、电子信息工程、数学等专业的学生与研究人员提供一套Matlab实现的WOA-VMD鲸鱼算法优化变分模态分解完整方案可用于课程设计、期末大作业与毕业设计也适合作为信号处理方向的创新点。资源利用鲸鱼优化算法自动寻优VMD中的模态数k与惩罚因子alpha输出边际谱、频率图、收敛曲线及超参数迭代过程图直观呈现每次迭代的变化趋势。代码采用参数化编程注释清晰数据为Excel格式便于替换运行主程序main即可直接执行。压缩包共1个文件为1份pdf文档大小约600KB内容涵盖算法原理、程序设计思路与参考资料便于快速理解整体框架。目前已有1334人学习下载读者可借此掌握WOA与VMD结合的核心流程、参数寻优策略及分解效果评估方法为后续信号分析与特征提取工作提供可复用的实现参考。1. 当VMD的K值靠猜时WOA-VMD鲸鱼算法优化变分模态分解到底在做什么如果你做过轴承故障诊断、心电去噪或者任何跟振动信号打交道的活大概率被VMD变分模态分解坑过分解层数K和惩罚因子α这两个参数调一次跑一次K大了模态混叠K小了欠分解α偏了直接过平滑。很多人靠经验拍脑袋或者拿网格搜索硬扫扫一晚上出来的结果还不一定稳。WOA-VMD鲸鱼算法优化变分模态分解干的就是把这两个参数交给鲸鱼优化算法WOA自动搜用包络熵或者样本熵当适应度函数让算法自己找到一组让分解结果最“干净”的K和α。这套组合在旋转机械故障特征提取里已经算常见做法适合手上有一维时序信号、想做模态分解但不想手动调参的人。下面从原理到代码把这条链路拆开讲清楚。2. WOA和VMD为什么能凑成一对参数搜索的底层逻辑2.1 VMD的两个参数到底卡在哪VMD的核心思想是把一个实值信号分解成K个有限带宽的模态分量IMF每个模态围绕自己的中心频率。它的变分问题构造里有两个关键参数模态数K和二次惩罚因子α。K决定你允许分解出多少个分量α决定每个分量的带宽约束强度。K设小了原本独立的频率成分被塞进同一个模态叫欠分解K设大了同一个频率成分被拆到两个模态里叫过分解还会引入虚假分量。α太小模态带宽过宽噪声混进来α太大带宽过窄模态失真。传统做法是固定α只扫K看哪个K下没有中心频率相近的模态就停。但α和K是耦合的固定一个搜另一个找到的是局部最优。这就是为什么需要把两个参数一起优化。2.2 鲸鱼优化算法凭什么适合搜这两个参数WOA是模仿座头鲸气泡网捕食的群智能算法核心行为有三步包围猎物、气泡网攻击、随机搜索。它的位置更新公式里有一个随迭代从2线性降到0的参数a控制探索和开发的切换。相比粒子群和遗传算法WOA的优势是参数少、收敛快、不容易早熟——当然这是相对而言后面避坑章节会讲它什么时候也会翻车。把WOA用在VMD参数搜索上映射关系很直接每条鲸鱼的位置是一个二维向量[K, α]K取整数范围比如2到10α取实数范围比如200到5000。适应度函数用包络熵对每个IMF算包络再算包络的熵值取所有IMF包络熵的最小值作为该位置的适应度。包络熵越小说明模态越稀疏、特征越集中分解质量越好。2.3 适应度函数选包络熵还是样本熵包络熵的计算步骤是对每个IMF做Hilbert变换取包络对包络做归一化得到概率分布再算熵。样本熵则是衡量时间序列复杂度值越小表示序列越规律。两者在VMD优化里都有人用区别在于适应度函数计算量对噪声敏感度适用场景包络熵中等较敏感故障冲击特征明显的信号样本熵较大一般复杂度变化明显的信号排列熵小较敏感短数据、实时性要求高我一般先用包络熵因为故障诊断里冲击成分的包络特征最直接。如果信号长度短于2000点包络熵的Hilbert变换边缘效应会明显这时候换排列熵更稳。3. 用Python把WOA-VMD跑通从适应度函数到主循环3.1 环境准备与VMD实现选择Python里没有scipy级别的官方VMD常见做法是用vmdpy这个包或者自己按原始论文的ADMM流程写。我倾向自己写一份精简版因为优化过程中要反复调用VMD第三方包的接口封装有时会拖慢速度。下面这份VMD实现只保留核心ADMM迭代够用且可控。import numpy as np from scipy.signal import hilbert from scipy.fft import fft, ifft def vmd(signal, K, alpha, tau0, DC0, init1, tol1e-7, max_iter500): 精简版VMD实现 signal: 输入一维信号 K: 模态数 alpha: 惩罚因子 tau: 噪声容限对含噪信号可设0 DC: 是否包含直流分量 init: 中心频率初始化方式1为均匀分布 tol: 收敛容差 max_iter: 最大迭代次数 N len(signal) t np.arange(N) / N # 频域初始化 freqs t - 0.5 - 1/N f_hat fft(signal) f_hat_plus f_hat.copy() f_hat_plus[:N//2] 0 # 只保留正频率部分 # 初始化模态频谱和中心频率 u_hat_plus np.zeros((max_iter, K, N), dtypecomplex) omega_plus np.zeros((max_iter, K)) if init 1: for i in range(K): omega_plus[0, i] (0.5 / K) * i else: omega_plus[0, :] np.sort(np.exp(np.log(0.5) np.log(N) * np.random.rand(K))) lambda_hat np.zeros((max_iter, N), dtypecomplex) u_diff tol 1 n 0 while u_diff tol and n max_iter - 1: for i in range(K): # 计算残差 sum_others np.sum(u_hat_plus[n, :, :], axis0) - u_hat_plus[n, i, :] # 更新模态频谱 numerator f_hat_plus lambda_hat[n, :] / 2 - sum_others denominator 1 alpha * (freqs - omega_plus[n, i]) ** 2 u_hat_plus[n1, i, :] numerator / denominator # 更新中心频率 power np.abs(u_hat_plus[n1, i, N//2:]) ** 2 omega_plus[n1, i] np.sum(freqs[N//2:] * power) / (np.sum(power) 1e-16) # 更新拉格朗日乘子 lambda_hat[n1, :] lambda_hat[n, :] tau * (np.sum(u_hat_plus[n1, :, :], axis0) - f_hat_plus) # 收敛判断 u_diff np.sum(np.abs(u_hat_plus[n1, :, :] - u_hat_plus[n, :, :]) ** 2) / (np.sum(np.abs(u_hat_plus[n, :, :]) ** 2) 1e-16) n 1 # 重构时域模态 u np.zeros((K, N)) for i in range(K): u[i, :] np.real(ifft(u_hat_plus[n, i, :])) return u, omega_plus[n, :]这段代码的关键点f_hat_plus只保留正频率因为实信号的负频率是共轭对称的省一半计算量。omega_plus的初始化用均匀分布比随机初始化更稳。收敛判断用模态频谱的相对变化量tol设1e-7是常用值再小收益不大。3.2 包络熵适应度函数的实现细节包络熵是WOA-VMD里最容易被写错的地方。常见错误是直接对IMF取绝对值当包络正确做法是Hilbert变换取解析信号的模。def envelope_entropy(imf): 计算单个IMF的包络熵 imf: 一维模态分量 # Hilbert变换取包络 analytic hilbert(imf) envelope np.abs(analytic) # 归一化 envelope envelope / (np.sum(envelope) 1e-16) # 避免log(0) envelope envelope[envelope 0] # 计算熵 entropy -np.sum(envelope * np.log(envelope)) return entropy def fitness_function(position, signal): WOA适应度函数 position: [K, alpha] signal: 原始信号 K int(np.round(position[0])) alpha position[1] # 边界约束 K np.clip(K, 2, 10) alpha np.clip(alpha, 200, 5000) try: u, _ vmd(signal, K, alpha) # 计算所有IMF包络熵的最小值 entropies [envelope_entropy(u[i, :]) for i in range(K)] return np.min(entropies) except Exception: # 分解失败返回大值 return 1e10这里取np.min(entropies)而不是均值是因为我们关心的是最优模态的稀疏性均值会被差模态拉高导致优化方向偏。K用np.round取整因为VMD要求整数。边界约束用np.clip硬截断比惩罚函数简单有效。3.3 WOA主循环与参数设置WOA的标准流程不复杂但有几个参数直接影响搜索效果。def woa_optimize(signal, dim2, n_whales20, max_iter30, lbNone, ubNone): WOA主循环 dim: 优化维度这里是2K和alpha n_whales: 鲸鱼种群规模 max_iter: 最大迭代次数 lb, ub: 下界和上界数组 if lb is None: lb np.array([2, 200]) if ub is None: ub np.array([10, 5000]) # 初始化种群 positions np.random.uniform(lb, ub, (n_whales, dim)) fitness np.array([fitness_function(p, signal) for p in positions]) # 记录最优 best_idx np.argmin(fitness) best_pos positions[best_idx].copy() best_fit fitness[best_idx] for t in range(max_iter): # a从2线性降到0 a 2 - 2 * t / max_iter for i in range(n_whales): r1 np.random.rand() r2 np.random.rand() A 2 * a * r1 - a C 2 * r2 p np.random.rand() if p 0.5: if abs(A) 1: # 包围猎物 D abs(C * best_pos - positions[i]) positions[i] best_pos - A * D else: # 随机搜索 rand_idx np.random.randint(n_whales) D abs(C * positions[rand_idx] - positions[i]) positions[i] positions[rand_idx] - A * D else: # 气泡网攻击 D abs(best_pos - positions[i]) b 1 l np.random.uniform(-1, 1) positions[i] D * np.exp(b * l) * np.cos(2 * np.pi * l) best_pos # 边界处理 positions[i] np.clip(positions[i], lb, ub) # 重新计算适应度 new_fit fitness_function(positions[i], signal) if new_fit fitness[i]: fitness[i] new_fit if new_fit best_fit: best_fit new_fit best_pos positions[i].copy() print(fIter {t1}/{max_iter}, Best fitness: {best_fit:.6f}, fK{int(np.round(best_pos[0]))}, alpha{best_pos[1]:.1f}) return best_pos, best_fit参数设置上n_whales20和max_iter30是性价比比较高的组合。种群再大适应度函数调用次数线性增长但VMD本身计算量不小20条鲸鱼跑30代就是600次VMD分解普通笔记本大概几分钟。如果信号很长超过10万点建议先降采样或者截取一段。a的线性递减策略是WOA的标准做法但有些改进版用非线性递减在VMD优化里区别不大因为搜索空间只有二维不算复杂。3.4 完整调用示例与结果解读拿一段模拟的轴承故障信号来跑通整条链路。if __name__ __main__: # 生成模拟信号三个频率成分噪声 fs 1000 t np.arange(0, 1, 1/fs) # 低频成分 s1 0.8 * np.sin(2 * np.pi * 30 * t) # 中频成分 s2 0.5 * np.sin(2 * np.pi * 120 * t) # 高频冲击成分模拟故障 s3 0.3 * np.sin(2 * np.pi * 300 * t) * (1 0.5 * np.sin(2 * np.pi * 5 * t)) # 噪声 noise 0.1 * np.random.randn(len(t)) signal s1 s2 s3 noise # 运行WOA优化 best_pos, best_fit woa_optimize(signal, n_whales15, max_iter20) K_opt int(np.round(best_pos[0])) alpha_opt best_pos[1] print(f\n最优参数: K{K_opt}, alpha{alpha_opt:.1f}) print(f最优包络熵: {best_fit:.6f}) # 用最优参数做最终分解 u, omega vmd(signal, K_opt, alpha_opt) print(f中心频率: {omega})跑完之后看两个东西一是中心频率有没有两个模态靠得很近如果靠得近说明K还是偏大二是最优包络熵对应的IMF是不是集中在冲击成分所在的频段。如果分解出来的模态里有一个明显对应300Hz附近的冲击说明参数搜对了。4. 避坑与排查WOA-VMD翻车的五个典型场景4.1 适应度函数值一直不降输出全是初始值现象跑完WOA最优适应度和第一代差不多K和α几乎没动。原因最常见的是VMD内部迭代没收敛就返回了导致每次适应度计算都是随机值。检查vmd函数里的max_iter默认500在信号长的时候可能不够。另一个原因是tol设得太小1e-7在数值精度边缘改成1e-6更稳。解决在fitness_function里加一个判断如果VMD返回的模态能量和原始信号能量差超过10%直接返回大值避免无效分解污染搜索方向。4.2 K被优化到边界值10但实际信号只有3个成分现象WOA搜出来的K总是贴着上界。原因包络熵取最小值时K越大总有一个模态的包络熵会很小这是包络熵作为适应度函数的固有偏向。K10时至少有一个模态非常稀疏熵值自然低。解决适应度函数改成包络熵的均值加上一个对K的惩罚项比如fitness mean_entropy 0.05 * K。或者把K的上界根据信号频谱的峰值个数来定先做FFT看有几个明显峰K上界设峰值数加2。4.3 每次运行结果都不一样参数波动大现象同一段信号跑两次WOAK和α差很多。原因WOA本身是随机算法初始种群随机p和r1、r2都是随机的。如果适应度函数在最优附近比较平坦随机性就会被放大。解决固定随机种子np.random.seed(42)至少保证可复现。另外把n_whales加到30种群大了之后最优解的稳定性会明显提升。如果还不行说明适应度函数区分度不够换样本熵试试。4.4 VMD分解报错或返回全零模态现象fitness_function里抛异常或者返回的u全是零。原因α设得太小比如接近0时分母1 alpha * (freqs - omega)**2接近1更新公式退化成简单平均模态会塌缩。另外如果信号本身是常数或者接近零VMD的频域初始化会出问题。解决α的下界不要低于100我一般设200。信号做VMD之前先减去均值去掉直流分量。如果信号幅值特别小比如1e-6量级先归一化到[-1, 1]再分解。4.5 优化出来的K和α在测试集上表现差现象在训练信号上包络熵很低换一段信号用同样的K和α分解结果一塌糊涂。原因WOA-VMD优化的是特定信号的参数不是通用参数。每段信号的频率成分不同最优K和α自然不同。解决如果要做批量处理要么每段信号单独优化计算量大但准确要么先对同类信号做统计分析取一个折中的K和α作为固定参数。我一般对同类工况的信号取10段分别优化后取K的众数和α的中位数。5. 让WOA-VMD真正可用的三个进阶技巧第一个技巧是适应度函数的混合策略。单纯用包络熵最小值容易偏向大K我现在的做法是先算所有IMF包络熵的均值再算最小值适应度取0.7 * mean 0.3 * min。这样既考虑整体稀疏性又不完全忽略最优模态。实测在轴承故障信号上这个混合策略搜出来的K比纯最小值策略平均小1到2更接近真实成分数。第二个技巧是WOA的种群初始化用拉丁超立方采样代替均匀随机。均匀随机在二维空间里容易聚集拉丁超立方保证初始点在整个搜索空间里分布更均匀。实现上就是把np.random.uniform换成scipy.stats.qmc.LatinHypercube代码改动不超过五行但收敛速度能快20%左右。第三个技巧是分阶段搜索。先固定α2000用WOA只搜K搜到一个粗略的K值然后固定K用WOA搜α最后两个参数一起微调。这样做的好处是避免二维搜索时K的整数跳变干扰α的连续优化。分阶段的总计算量比直接二维搜索少因为每阶段的搜索空间是一维的。验证优化结果是否靠谱我习惯用中心频率重叠度来判断把优化后的K个中心频率排序如果相邻两个的差值小于频率分辨率的1.5倍就认为存在模态重叠需要减小K或者增大α。这个判断比看包络熵直观得多。最后说一个我踩过的坑有次跑WOA-VMD适应度函数忘了做边界截断K被优化到负数VMD直接报错但异常被try-except吞了返回1e10结果WOA以为找到了最差位置反而把种群往负数方向推。排查了半天才发现是边界处理漏了。从那以后我在fitness_function第一行就写np.clip不管后面怎么改先把参数锁在合法范围里。希望帮到你。本文还有配套的精品资源点击获取