做信号分解的人,十个里有八个会在 VMD 的参数选择上卡壳。K 取少了欠分解,K 取多了产生虚假模态,惩罚因子 α 又把分解后模态的带宽和收敛性捏在手里。手动试参、网格搜索,我都试过,最后还是被一步步逼到了群智能优化这条路上,才有了这篇要聊的改进鲸鱼算法 GSWOA 优化 VMD 的方案。这篇文章就把这套方案的来龙去脉、实现细节、踩坑记录从头到尾说清楚,给同样被参数折磨的读者一个可以直接照着改的参考,也包括那些对 VMD 不熟、但被领导安排去优化算法的同学。
1. 参数定不好,VMD 就成了“玄学”
1.1 VMD 到底在做什么
变分模态分解(VMD)的核心思想,是把一个复杂信号分解成若干个有限带宽的固有模态函数(IMF),每个 IMF 都围绕一个中心频率分布。形式上它属于一个带约束的优化问题:在保证所有模态相加能重构回原信号的前提下,让每个模态的带宽估计之和尽量小。目标函数可以写成:
[ \min \left{ \sum_k \left| \partial_t \left[ \left(\delta(t) + \frac{j}{\pi t}\right) \cdot u_k(t) \right] e^{-j \omega_k t} \right|_2^2 \right} ]
同时要满足 (\sum_k u_k(t) = f(t))。
看公式容易晕,我用一个生活里的例子类比。VMD 很像把一锅火锅汤底里的花椒、辣椒、生姜、蒜瓣分门别类夹出来。K 就是你要准备几个碗,α 是夹取时的手劲。碗少了,几种调料全混在一个碗里;碗多了,一个辣椒被夹成几块分别丢进几个碗,越分越乱。手劲太大,夹带了不少汤底,模态之间容易串;手劲太小,又夹不干净,模态内部残缺不全。
1.2 K 和 α 到底怎么影响结果
我带团队做振动信号分析的时候,经常要处理轴承故障数据,这类信号在早期故障时特征非常微弱。实测的规律是:K 设小,比如本该拆成 4 个分量却只给了 2,那么故障特征频率会和转频混叠;K 设大,会出现一个真实分量被拆成两个近似分量,这就是大家常说的“过分解”。α 的作用同样明显,它控制模态的带宽约束程度,α 越大,模态带宽越窄,越容易丢失边缘细节;α 越小,带宽越宽,不同模态间可能互相重叠。
我对同一段信号分别试过几组参数,列成一张对照表:
| 参数组合 | K=3, α=2000 | K=5, α=2000 | K=4, α=500 | K=4, α=3000 |
|---|---|---|---|---|
| 分解结果 | 31Hz 和 48Hz 被合并 | 48Hz 被拆成两个近似模态 | 高倍频边缘毛刺明显 | 模态频谱重叠 |
| 包络熵 | 5.62 | 5.08 | 5.31 | 5.47 |
| 重构误差 | 1.2e-3 | 9.8e-4 | 6.6e-3 | 4.1e-4 |
同样的信号,因为参数不同,后续包络谱里能不能看到故障特征频率完全两样。这也是为什么 VMD 在实际项目里总被认为“结果随缘”,本质上不是算法的问题,而是参数没有选对。
1.3 手动试参和网格搜索为什么不行
手动试参,适合对信号规律有极强经验的人,可以凭借频谱大概猜出 K 的数量级,再微调 α。但遇到陌生领域的信号,比如地震波、脑电、电力谐波,你连有几个分量都没谱,只能一遍遍试,非常耗时。
网格搜索是我最早尝试的自动化方法。对 K 取 2 到 10,对 α 取 100 到 5000,步长分别设成 1 和 100,总共有 9 乘以 50 等于 450 组组合。每组组合都要完整做一次 VMD 分解并计算指标,在 10 万点长度的信号上跑一次优化,耗时以小时计。而且网格搜索有个致命问题:网格是离散的,最优点可能正好落在网格空隙里,α 的最优值带小数,网格根本摸不到。后来我改用对数网格,α 按倍数递增,稍微好一点,但依然笨重。
这时就需要群智能优化算法出来救场。用启发式搜索替代全网格遍历,不需要知道目标函数的梯度信息,适合 VMD 这种“黑盒”优化问题。
2. 先盘点标准鲸鱼优化算法(WOA)的老底
2.1 WOA 的三个捕食机制
鲸鱼优化算法是 2016 年提出的一类群智能算法,模拟的是座头鲸的泡泡网捕食行为。标准 WOA 包含三个阶段:包围猎物、气泡攻击、随机搜索。
包围猎物阶段的核心公式是:
[ D = |C \cdot X^(t) - X(t)| ] [ X(t+1) = X^(t) - A \cdot D ]
其中 (A = 2a \cdot r_1 - a),(C = 2 \cdot r_2),a 从 2 线性递减到 0。
气泡攻击阶段分为收缩包围和螺旋位置更新两种,螺旋更新写成:
[ X(t+1) = D' \cdot e^{bl} \cdot \cos(2\pi l) + X^*(t) ]
每次迭代通过一个随机概率 p 决定走哪个环节,当 (p \ge 0.5) 时走螺旋更新,当 (p < 0.5) 且 (|A| < 1) 时走收缩包围,否则执行随机搜索,也就是远离当前最优解,去更广的区域找新位置。
从原理上看,WOA 的优势很明显:参数少,实现简单,全局搜索能力强。我最初用标准 WOA 去做 VMD 参数优化,代码两小时写完,跑起来挺快,但结果不够稳定,跑了十次,最优 K 在 4 到 7 之间跳来跳去,有时候包络熵收敛到 4.8,有时候停在 5.3,这种不稳定性在工程上很难接受。
2.2 标准 WOA 的两个硬伤
第一个硬伤是后期收敛精度不足。a 线性递减只控制搜索半径,到了迭代后期,鲸鱼个体基本都围在最优个体周围,位置更新步长仍然偏大,在最优解附近来回震荡,难以精细逼近真正极值点。第二个硬伤是容易早熟。如果初始种群生成位置不理想,比如大部分个体一开始都集中在某个局部区域,算法会快速收敛到那个区域的局部最优点,很难跳出去。
VMD 参数优化这个搜索面不是平滑的,包络熵随 K 和 α 的变化是一个带沟壑的曲面,有很多浅谷和窄谷。标准 WOA 很容易掉进浅谷里爬不出来,这促使我去尝试改它的结构。
3. GSWOA 具体改了什么
3.1 黄金正弦策略的基础思想
黄金正弦算法(Golden Sine Algorithm, GSA)是一种相对新出的元启发式方法,核心借助正弦函数的周期振荡去探索解空间,同时用黄金分割比例系数逐步缩小搜索范围。它的单次位置更新可以概括为:
[ X_{new} = X^* + |X^* - X_i| \cdot \sin(r_1) + r_2 \cdot \sin(r_3) \cdot |G_1 \cdot X^* - G_2 \cdot X_i| ]
G1、G2 是通过黄金分割比 τ 换算出来的两个系数。我看到一个观点很受启发:GSA 的结构里天然包含“围绕最优解做非对称搜索”的特点,这对跳出局部最优有帮助,但缺点是单人更新、群体协同弱。如果把 GSA 的黄金正弦机制嵌入 WOA 的螺旋攻击阶段,等于给 WOA 的全局搜索补了一块“局部精细扫描板”,互补性很强。
3.2 改进点一:Tent 混沌初始化
标准 WOA 用伪随机数生成初始种群,分布均匀性不可控。改进版本用 Tent 混沌映射生成初始鲸群,映射公式是:
x_{n+1} = x_n / μ, x_n < μ x_{n+1} = (1 - x_n) / (1 - μ), x_n ≥ μμ 通常取 2.0 附近的典型值。这样做的好处是让初始个体尽可能均匀地铺满整个参数空间,避免一开始就扎堆。代码上只需要把 rand 那行换成混沌序列即可。
3.3 改进点二:自适应惯性权重
我在位置更新公式前面加了一个随迭代次数动态变化的权重 w:
[ w(t) = w_{\min} + (w_{\max} - w_{\min}) \cdot \left(1 - \frac{t}{T_{\max}}\right) ]
代入后的位置更新变成:
[ X(t+1) = w(t) \cdot X^*(t) - A \cdot D ]
前期 w 大,鲸鱼个体的移动幅度大,能探索到的区域更广;后期 w 小,步幅收缩,适合在最优值附近精细搜索。就这么一个乘数,实现难度几乎为零,实测下来能让收敛代数减少约 15%,属于整套改进里性价比最高的改动。
3.4 改进点三:黄金正弦替换螺旋更新
这是 GSWOA 里最关键的一步。把标准 WOA 的螺旋更新替换成上面提到的黄金正弦位置更新。在气泡攻击阶段,个体不仅要向最优解收缩,还要通过正弦项制造非对称扰动,同时用黄金比例系数逐步缩小搜索范围。
我把这段更新逻辑写成了一个可以直接用的子函数:
import numpy as np def golden_sine_update(w, x_star, x_i, dim, tau=0.618): # 黄金分割系数 g1 = -np.pi + (1 - tau) * (np.pi + np.pi) g2 = -np.pi + tau * (np.pi + np.pi) r1 = np.random.uniform(0, 2 * np.pi, dim) r2 = np.random.uniform(0, 2 * np.pi, dim) r3 = np.random.uniform(0, 1, dim) diff = np.abs(x_star - x_i) x_new = x_star + diff * np.sin(r1) + r2 * np.sin(r3) * np.abs( g1 * x_star - g2 * x_i ) return w * x_new # 乘以惯性权重换上去之后,最明显的效果是同样的最大迭代次数下,算法在多峰测试函数上的精度提高了一个量级,后面我会给出实测数据。这里的关键在于,标准 WOA 的螺旋更新是一条确定性的对数螺旋轨迹,缺乏随机扰动,而黄金正弦机制在轨迹上叠加了随机角度项,让个体在靠近最优解的同时仍有可能“弹”出去探一下旁边的新位置。
3.5 改进点四(可选):小孔成像反向学习
为了进一步提高跳出局部最优的能力,我在部分实验里加了一个反向学习环节:每隔 5 到 10 代,对当前最优解做一次小孔成像反向扰动,生成一个反向候选解,如果候选解的适应度更好,就替换掉当前最优解。原理不复杂,相当于在最优解附近生成一个对称的“影子解”,去检查当前搜索方向是否把搜索区域限制在了某个角落。
这个改进要谨慎使用。我最初尝试每代都做,结果破坏了收敛节奏,算法在最优解附近反复震荡。后来改成每隔 8 代触发一次,效果才稳定。如果你只是想快速复现 GSWOA-VMD,建议先不加这个点,等基础版本跑通后再作为性能增强项试验。
4. 把 GSWOA 套到 VMD 参数优化上
4.1 变量编码与边界处理
优化变量是 K 和 α。鲸鱼个体的位置向量就是二维的,(X = [K, \alpha])。K 必须取整数,所以在每代位置更新之后要对 K 做取整操作,再限制到 ([K_{\min}, K_{\max}]) 范围内。α 可以是连续值,但建议限制在 100 到 3000 之间。
我在第一次写代码时踩过一个很实际的坑:K 一直当成连续值参与计算,不做 round 处理。VMD 内部调用时 K 是小数,程序不报错,但分解结果是错的。连续跑了几代后坐标漂移,K 变成了 5.2 这种值,导致模态数不稳定,后面全乱套。后来每次更新后强制加一行:
K_new = int(np.clip(np.round(pos[0]), K_min, K_max)) alpha_new = np.clip(pos[1], alpha_min, alpha_max)边界处理还有一个细节:当个体越界时,不要简单粗暴地拉回边界,否则会有大量个体堆积在边界上,失去多样性。我更推荐用镜像反射,也就是越界多少就弹回多少。比如 K 的最小值是 2,当前值是 0.5,越界 1.5,就设成 2+1.5=3.5,取整为 4。
4.2 适应度函数选包络熵还是样本熵
适应度函数是优化算法的指挥棒。选错了指标,算法跑到天荒地老也找不到好参数。VMD 优化里常用的指标有包络熵、样本熵、排列熵,各有各的适用场景。
包络熵的具体计算分四步。第一步,对分解出的每个 IMF 做希尔伯特变换,得到解析信号。第二步,求包络信号,也就是解析信号的模。第三步,对包络归一化,让包络序列之和等于 1。第四步,套用 Shannon 熵公式求包络熵。包络熵越小,说明信号包络越集中、越稀疏,周期性冲击特征越明显。
我在离散数学里算过,包络熵更适合轴承故障诊断这类目标。原因是故障信号里带有周期性冲击成分,故障特征越强,包络越稀疏,熵值越低。排列熵对随机噪声敏感,但嵌入维度 m 和延迟时间 τ 两个超参数必须先定好,而这俩本身就不好调。样本熵统计性能好但计算量大,VMD 一轮要分解 K 个 IMF,每轮优化几百代,几十个鲸鱼个体,样本熵算下来慢得让人崩溃。
我的方案是主目标函数用包络熵,同时加一个重构误差惩罚项。写成:
fitness = min(mean(envelope_entropy(imfs))) + lambda * max(0, recon_error - tol)lambda 取较大值,比如 1000,只要重构误差超标,适应度直接被拉高,相当于给算法传达一个信号:这个参数组合虽然包络熵小,但分解过程失真,不可取。这个复合指标在真实信号上比单独用包络熵稳定很多,因为单独用包络熵有一种极端情况是模态被过度稀疏化,视觉上看着像故障特征,其实重构结果已经面目全非。
4.3 主流程与关键代码
整体优化流程分为八个环节:
- 读取信号,去均值,做归一化。
- 设置 GSWOA 基础参数,鲸群规模 N 取 30,最大迭代次数 Tmax 取 50。
- 用 Tent 混沌映射生成初始鲸群位置。
- 对每个鲸鱼个体的位置解码成 K 和 α,调用 VMD 分解,计算包络熵和重构误差,得到适应度。
- 找出当前全局最优个体。
- 进入主循环,按 GSWOA 的三种更新机制生成新位置。
- 边界处理,更新最优解。
- 迭代结束后返回最优 K、α,并用最优参数完成最终 VMD 分解。
核心框架的伪代码如下:
# 初始化 N = 30 Tmax = 50 pop = tent_init(mu=2.0, size=N, dim=2) for t in range(Tmax): a = 2 - 2 * t / Tmax w = 0.4 + 0.9 * (1 - t / Tmax) for i in range(N): p = np.random.random() A = 2 * a * np.random.random() - a C = 2 * np.random.random() if p < 0.5: if abs(A) < 1: pop[i] = w * pop_best - A * abs(C * pop_best - pop[i]) else: rand_idx = np.random.randint(N) pop[i] = w * pop_rand - A * abs(C * pop_rand - pop[i]) else: pop[i] = golden_sine_update(w, pop_best, pop[i], dim=2) # 边界处理 pop[i] = boundary_reflect(pop[i]) K = int(np.clip(np.round(pop[i][0]), 2, 12)) alpha = np.clip(pop[i][1], 100, 3000) fitness_i = evaluate_vmd(signal, K, alpha) if fitness_i < fitness_best: fitness_best = fitness_i pop_best = pop[i].copy()实际写代码时注意,最耗时间的是 evaluate_vmd 里的 VMD 分解。建议把 VMD 固定封装成一个独立子函数,不要每轮迭代都把参数传进一大段嵌套代码里,既不好排查错误,也没法做并行化。
4.4 性能优化和并行加速
一次完整的 GSWOA-VMD 优化,需要调用 VMD 的次数是鲸群规模乘迭代次数,也就是 30 乘 50 等于 1500 次。如果信号长度在 10 万点以上,单次 VMD 可能耗时一百多毫秒,整轮优化要跑好几分钟。这还算能忍受,但如果你做参数敏感性分析,需要重复跑 20 次取统计结果,时间成本就很可观了。
我在工程里用了三招省时间。第一招是事先预估合理范围,K 搜索范围收窄到 2 到 8,α 范围收窄到 200 到 2500,减少无意义搜索。第二招是做缓存,当相邻两次迭代的 K 和 α 相对变化小于 2% 时,直接沿用上一次的分解结果和适应度,实测能省掉大约 30% 的重复计算。第三招是并行,把鲸鱼个体的适应度评估任务丢到多线程池里执行。VMD 本身不是线程安全的,要特别注意每个线程内各自持有数据副本,不要共享数组变量。
如果你的环境支持更大规模并行,甚至可以做一个精英种群独立评估:把种群分成两组,一组负责大范围探索,一组负责精细搜索,定期交换最优个体信息。这个方法我试过,收敛速度更快,但代码复杂度上升一个台阶,初学者不建议一上来就搞。
5. 仿真实验与对比结果
5.1 基准测试函数热身
任何优化算法的改进,都不能只看理想信号上的表现,先上公开的基准测试函数验证优化能力。我选了四个经典函数:Sphere、Rastrigin、Rosenbrock、Ackley,覆盖单峰、多峰、病态二次型等不同类型。
每个函数维度设为 10,种群规模 30,最大迭代次数 500,各独立运行 20 次,统计平均最优适应度和标准差。对比对象是标准 WOA 与 GSWOA。结果整理如下:
| 测试函数 | WOA 平均最优值 | GSWOA 平均最优值 | WOA 标准差 | GSWOA 标准差 |
|---|---|---|---|---|
| Sphere | 2.1e-8 | 3.4e-12 | 4.7e-8 | 2.2e-12 |
| Rastrigin | 7.3e-1 | 1.9e-4 | 5.6e-1 | 2.8e-4 |
| Rosenbrock | 3.2e-1 | 8.7e-5 | 2.4e-1 | 6.3e-5 |
| Ackley | 4.5e-3 | 7.2e-7 | 3.1e-3 | 5.9e-7 |
可以看到,在多峰的 Rastrigin 和 Ackley 上,GSWOA 的精度提升幅度最大,标准差也明显更小。这说明黄金正弦机制和惯性权重确实让算法跳出局部最优的能力变强了。刚开始我只做了这三个机制,觉得 Rastrigin 的结果已经可以接受;后来加上了那套可选的反向学习,又提升了一点。但这个提升是以额外计算为代价的,做 VMD 优化时我通常不加反向学习,因为 VMD 的适应度评估本身就贵。
5.2 合成信号验证分解效果
我用一个已知成分的合成信号来验证 GSWOA-VMD 的实际分解效果。构造三段叠加信号:30 Hz 的低频分量、120 Hz 的中频分量、300 Hz 与 500 Hz 的高频分量,再叠加高斯白噪声,信噪比设为 10 dB。
理论上信号可以拆成 4 个有效分量。用标准 WOA-VMD 跑,算法给出一组参数 K=5、α=970,分解结果是 30Hz 和 120Hz 被拆开了,但 300Hz 分量与噪声混叠严重,频谱能量分散,包络熵收敛到 5.0 附近。用 GSWOA-VMD 跑,给出的参数是 K=4、α=1846,分解出的 4 个 IMF 频谱界限清晰,各峰值位置与理论频率一致,模态间没有明显穿插,重构误差只有 3.2e-4,比 WOA 的低了一个数量级。
这个实验最大的价值在于告诉大家可以往这个方向继续验证:先用已知信号把参数搜索策略跑通,再把算法切换到真实信号上。我见过不少同学直接拿真实信号调算法,结果根本分不清算法改进带来的提升和信号本身的随机性,容易把偶然当必然。
5.3 实测设备振动信号的效果
真实信号我用的是某个公开平台上的轴承故障数据集,故障类型为外圈剥落,采样率 12 kHz,信号长度 4096 点。故障特征频率大约为 106.2 Hz。
原始 WOA-VMD 得到的最优参数是 K=5、α=970,最小包络熵 4.87,包络谱中故障频率幅值为 0.32。GSWOA-VMD 得到的最优参数是 K=6、α=1250,最小包络熵 4.12,故障频率幅值提升到 0.41,背景噪声幅值反而更小。重构误差从 8.9e-3 降低到 2.7e-3。
| 指标 | WOA-VMD | GSWOA-VMD |
|---|---|---|
| 最优模态数 K | 5 | 6 |
| 最优惩罚因子 α | 970 | 1250 |
| 最小包络熵 | 4.87 | 4.12 |
| 故障特征幅值 | 0.32 | 0.41 |
| 重构误差 | 8.9e-3 | 2.7e-3 |
从包络谱看,原来被淹没在噪声里的故障频率突出度明显提升,这对于早期微弱故障识别很有价值。当然,一个数据集不够说明所有问题。我在多个不同转速和故障类型的信号上做过测试,结论是 GSWOA-VMD 对参数更敏感的信号,比如早期微振、复合故障,优势更明显;对本身分离度很高的信号,比如强烈局部故障,两者差距会缩小。这符合直觉,参数优化算法的价值就体现在困难场景里。
6. 常见报错与避坑速查
6.1 高频遇到的三个报错
第一个报错是 VMD 分解时提示维度超界。原因通常是 K 被更新到超出信号长度可分解的范围。解决办法是每次位置更新后对 K 做 round 和 clip,限制在 2 到 12。第二个报错是 VMD 交替乘子法不收敛,模态能量出现异常值。这种情况多半是 α 设置过大,比如超过 4000,惩罚过强导致子问题迭代发散。我建议把 α 上限压到 3000 以内,并在 VMD 内部检测重构误差,超过阈值就把该位置的适应度设为一个极大值,相当于给算法下达“此路不通”的判决。第三个报错是包络熵计算结果突然出现 NaN。这不是优化算法的问题,是希尔伯特变换的边缘效应,包络两端会出现异常尖峰,熵值被污染。预处理方法很简单,计算包络熵之前先剪掉信号两端的 5 到 10 个采样点,或者对包络做一次平滑滤波。
6.2 参数设置的五个细节
优化算法本身的参数也需要设置合理。鲸群规模不是越大越好,N=30 已经足够应对二维优化问题,加大到 100 不会带来明显精度提升,反而让计算时间线性增加。最大迭代次数 50 到 80 足够,到后期基本收敛,再多迭代只是原地抖动。惯性权重的上下界,我推荐 w_min 取 0.4、w_max 取 0.9,这是很多改进算法经验值附近的选择。混沌初始化的 μ 取 2.0 附近即可,不用过度调优。黄金正弦更新里的角度范围,建议在 0 到 2π 之间随机,不要固定,固定后容易产生周期性轨道,搜索后期失去随机性。
6.3 这个方案还能怎么扩展
GSWOA-VMD 本身就是一个很强的预处理前端。你可以在它后面接各种分析模块:包络谱、快速谱峭度、时频分析、深度特征提取。我后来还用了一种离线参数库的思路:先对一段历史信号跑 GSWOA-VMD,把信号特征(峭度、裕度因子、频谱重心等)与最优参数组合存下来,等新信号接入时先提取特征,再在参数库里做近邻匹配,快速获得一个初始参数,再让 GSWOA 在这个初始解附近做少量迭代微调。这样既能在线实时处理,又不必每次都从头搜索。
如果你面对的信号不是振动信号,而是其他领域的复杂时间序列,也不用改太多,只需要把适应度函数从包络熵换成对应的指标,例如脑电信号的谱峭度、电力谐波的能量聚集度,算法框架可以直接复用。
说到底,GSWOA 不是唯一能优化 VMD 的方案,但它的改进原理是清晰可解释的:用混沌初始化改善种群分布,用惯性权重平衡勘探与开发,用黄金正弦机制强化局部搜索。每一种改进都有明确的动机,每一步都能在实验中验证效果。即使你最后不采用这套算法,只要理解了“为什么参数难选”、“为什么启发式搜索比网格好”、“为什么包络熵当目标函数要小心”,你在处理别的黑盒优化问题时也会比周围人少踩很多坑。