正余弦算法优化VMD参数：信号分解自动寻优方案

简介：一份面向信号处理与数据分析人员的Python实现包，聚焦正余弦算法(SCA)对变分模态分解(VMD)关键参数的自动优化。压缩包内共9个文件，包含4个xml工程配置、1个py核心算法脚本、1个txt示例数据，以及若干IDE辅助配置文件，整包仅631KB，轻量易用。VMD可将复杂信号分解为多个固有模态函数，但分解效果强烈依赖中心频率、罚项因子等参数；SCA通过初始化种群、计算适应度并利用正弦余弦函数更新位置，可快速搜索出使各IMF分量更纯净、彼此干扰更小的参数组合。代码开箱即跑，适合希望理解VMD参数寻优机制或开展时间序列分析的读者；同时借助工程配置文件，可在PyCharm中直接调试，通过对比不同参数下的频率特性与能量分布，直观验证优化效果，并观察SCA收敛过程及参数调整对分解结果的影响规律。目前已有661人学习，兼具教学演示、论文复现与工程预研价值，是一份轻量而完整的参考代码包。

1. 正余弦算法优化变分模态分解参数：信号分解别再做“盲人摸象”

拿到一段轴承振动信号，第一件事通常是做变分模态分解（VMD），把混在一起的频率成分拆开。可 VMD 的 K（模态数）和 penalty 参数 alpha 一旦设错，分解结果就全乱了：K 设小了，多个频率成分挤进同一个模态里；K 设大了，同一个成分被劈成两半，还会出现一堆没有物理含义的虚假模态。以前我调这两个参数靠手工试凑，从 K=3 试到 K=12，每组跑一遍对照，效率极低，结果还全凭经验。正余弦算法（SCA）是一种结构很简单的群体优化算法，把 K 和 alpha 当成两个连续变量，用包络熵做适应度函数，让算法自己去搜最优组合。全套逻辑用 Python 实现，从读数据到输出最优参数一共不过一百多行代码。这篇笔记适合正在做信号分解、故障诊断或者特征提取的读者，照着跑通，就能把人工试凑换成自动寻优。

2. VMD 参数为什么要优化：K 和 alpha 到底在控制什么

2.1 变分模态分解的核心参数与取值边界

先明确一个事实：VMD 不是一劳永逸的分解工具，它的效果完全取决于两个超参数。K 是模态个数，决定把信号拆成几份；alpha 是二次惩罚系数，作用在带宽约束上。alpha 越小，每个模态的带宽约束越松，模态中心频率附近的频率分量更容易被吸收；alpha 越大，带宽越窄，模态越“纯净”，但也可能把同一个物理成分截断。工程上最常见的问题恰恰在这里——这两个参数互相牵扯，K 变了，alpha 的合适区间也跟着变。比如信号里有 50Hz 和 120Hz 两个正弦分量，K=2 同时把 alpha 从 500 拉到 3000，分解出来的两个模态中心频率可能没变化，但模态 1 的带宽明显收窄，波形畸变更严重。

我用过一段时间的固定参数经验表，大致是这样的边界效应：

这个表只作为起点参考，真正的问题在于：每个实际信号的“合适区间”都不同，这就是为什么需要算法去搜。SCA 做的事情就是把 K 和 alpha 当作两个决策变量，在给定边界内找一组使目标函数最优的组合。因为正余弦算法是连续优化算法，K 在代码里需要取整后传给 VMD 函数，这是个容易翻车、但很容易躲开的细节，后面第 4 章会落到代码里。

2.2 适应度函数怎么选：包络熵是最稳妥的默认项

有了 K 和 alpha 还不够，优化还得有一个“好”的衡量标准。VMD 论文原始评价指标偏主观，做工程落地必须把它替换成数值目标。我试过排列熵、峭度、包络熵，最后日常方案固定在最小包络熵上。原理不复杂：对每个 IMF 做 Hilbert 变换得到包络，归一化后计算信息熵。如果分解结果里某个模态包含明显的冲击特征，它的包络会比较“稀疏”，熵值就低；如果模态里塞满了噪声，包络接近均匀分布，熵就高。取全部模态包络熵的最小值作为 SCA 的适应度函数，让算法去最小化它，方向上完全说得通。

有的读者习惯用峭度，但在 VMD 参数优化里，峭度有一个毛病：它只对冲击型信号敏感，对平稳信号或者调幅信号，峭度几乎不随 alpha 变化，目标函数太平，SCA 搜起来没梯度感。包络熵对信号形式不挑，是更通用的默认项。代码实现也就十几行：

import numpy as np from scipy.signal import hilbert from vmdpy import VMD def envelope_entropy(imf): # 希尔伯特变换得到解析信号，绝对值即包络 env = np.abs(hilbert(imf)) # 归一化，让包络序列变成一个概率分布 p = env / np.sum(env) # 包络熵：信息熵公式，加极小值防止 log(0) return -np.sum(p * np.log(p + 1e-12)) def fitness(params, signal, tau=0, DC=0, init=1, tol=1e-7): # params 来自 SCA 的位置向量：第一位是 K，第二位是 alpha K = int(np.clip(np.round(params[0]), 2, 12)) alpha = float(np.clip(params[1], 500, 5000)) # 调用 vmdpy 完成分解，返回三个值，只用第一个 u u, _, _ = VMD(signal, alpha, tau, K, DC, init, tol) # 对每个 IMF 算包络熵，取最小作为适应度 return min([envelope_entropy(u[k, :]) for k in range(K)])

这里有一个关键细节：np.clip把 K 按到 [2, 12]、alpha 按到 [500, 5000]，是为了防止 SCA 在探索过程中飞出边界。SCA 的位置更新本身不包含边界处理，如果不对每一代做裁剪，后面传给 VMD 的 K 可能是 200，直接导致分解耗时爆炸甚至u的维度异常。vmdpy库返回的u形状是 (K, N)，每个模态一行，循环取出来算熵即可。tau=0是 VMD 论文里的标准设置，表示不考虑噪声容差。

2.3 参数边界怎么定：别拍脑袋，先看信号再做范围

适应度函数定了，边界的设置是另一个影响结果的细节。我见过很多人直接把 K 定成 [1, 20]，alpha 定成 [1, 100000]，这种宽得离谱的区间让 SCA 白费大量迭代去探索完全没有意义的参数组合。实际经验是：K 的下限低于 2 没有意义，因为 VMD 至少应该把信号分成直流分量和主分量两路；K 的上限按采样频率和信号长度来粗估，对 1 秒 1000 点、成分不超过 4 个的信号，K 上限设到 10～12 已经足够。alpha 的边界更值得说：小于 500 时模态带宽过大，分解结果经常出现相邻模态互相“抢频段”；大于 5000 时收敛变慢，模态中心频率容易漂移出真实物理成分。800～3000 是 VMD 论文里举例的常用区间，落地时放宽到 500～5000 保留寻优空间。

还有一个工程经验：在 SCA 初始化种群时塞几个经验解进去，比如 (K=8, alpha=2000)，剩下的个体随机生成。这比全随机种群收敛快，因为正余弦算法本身没有很强的记忆能力，一个好的起点能明显减少前期无效探索。种群规模和迭代次数的基本盘是 pop=12、T=60，这个配置单条信号在普通笔记本上大约几十秒到两三分钟。如果想压时间，优先降迭代次数而不是降种群，因为种群太少容易让 SCA 直接陷进局部最优。

3. 正余弦算法的核心逻辑：四个随机参数和一个衰减因子

3.1 SCA 位置更新公式拆解：为什么它能兼顾探索和收敛

正余弦算法是 Mirjalili 在 2016 年提出的，原理解释起来比粒子群还简单：每个候选解的位置更新只参照当前全局最优解，用正弦或余弦函数来生成步进方向和步长。具体公式是，当随机数 r4 < 0.5 时，位置按正弦更新，否则按余弦更新。每次更新里有四个随机数起作用：r1 控制步长幅度，它在迭代过程中先从 2 线性衰减到 0，前期的探索步长大，后期逐步收敛；r2 是随机角度，范围在 [0, 2π]，决定移动方向；r3 给最优解位置加一个随机权重，范围 [0, 2]，避免每次移动完全重复；r4 则是决定这次用正弦还是余弦分支。这四个参数配合起来，让 SCA 前期大范围跳跃、后期逐步逼近最优解，算法结构没有任何微分和梯度依赖，所以特别适合像 VMD 参数这样黑匣子式的目标函数。

我把四参数的作用整理成一张表，方便对照代码理解：

理解 r1 的衰减最关键。它直接用当前迭代次数 t 和总迭代次数 T 计算，a=2 是标准配置。t 小的时候 r1 接近 2，步长大；t 接近 T 的时候 r1 接近 0，每个个体在最优解周围微调。这个衰减逻辑很实用，但要注意：它没有自适应的能力，如果 T 设得太小，后期还没来得及收敛就停了；如果 T 设得太大，后半段有大量迭代在做无谓的微调。

3.2 SCA 主循环代码：从初始化到迭代更新

SCA 的实现非常规整，我就直接给出一个适用于 VMD 参数寻优的标准循环。为了保证可复现，用numpy.random.default_rng管理随机源，支持固定随机种子，这个习惯能避免“每次运行结果不一样”的尴尬。

def sca_vmd(signal, pop=12, T=60, seed=42): rng = np.random.default_rng(seed) lb = np.array([2.0, 500.0]) # K 和 alpha 的下界 ub = np.array([12.0, 5000.0]) # K 和 alpha 的上界 # 种群初始化：一半随机，一半放经验参数 X = lb + rng.random((pop, 2)) * (ub - lb) X[0] = np.array([8.0, 2000.0]) # 经验起点，帮助早期收敛 X[1] = np.array([5.0, 1000.0]) # 另一组保守起点 fits = np.array([fitness(x, signal) for x in X]) best_idx = np.argmin(fits) best_pos = X[best_idx].copy() best_fit = fits[best_idx] a = 2.0 # r1 衰减初值 for t in range(T): r1 = a - t * (a / T) for i in range(pop): r2 = rng.uniform(0, 2 * np.pi) r3 = rng.uniform(0, 2) r4 = rng.random() if r4 < 0.5: X[i] = X[i] + r1 * np.sin(r2) * np.abs(r3 * best_pos - X[i]) else: X[i] = X[i] + r1 * np.cos(r2) * np.abs(r3 * best_pos - X[i]) # 边界裁剪：必须转回整数时要交给 fitness 处理 X[i] = np.clip(X[i], lb, ub) for i in range(pop): fits[i] = fitness(X[i], signal) if fits.min() < best_fit: best_idx = np.argmin(fits) best_pos = X[best_idx].copy() best_fit = fits[best_idx] return best_pos, best_fit

逻辑说明：每次迭代先对每个个体做一次位置更新，接着做边界裁剪，最后统一重新计算整个种群的适应度。与原版 SCA 每个个体更新完立即更新全局最优的做法不同，我用的是同步更新，好处是多个个体的适应度可以并行计算，后面第 6 章会说怎么接 joblib。best_pos里存的是连续值，比如 K=7.32、alpha=2134.7，要转成可执行的 VMD 参数，在调用fitness时已经做了round和clip，所以最终输出时要手动处理一遍。

3.3 种群和迭代次数的工程经验：先跑小规模验证，再上全量

SCA 的寻优效果对参数不算特别敏感，但有一个原则：优先保证种群多样性，而不是盲目增大迭代次数。种群设 8 的时候，我遇到过 SCA 在某个局部解附近反复打转的情况，因为个体太少，交叉覆盖能力太弱。种群到 12～20 之后，每增加一个个体的边际收益明显下降，但耗时线性增加，因为每个个体都要跑一次 VMD。迭代次数同理，60 次是我个人比较平衡的选择，少于 30 次往往还没完全收敛，多于 100 次则后半段的收敛收益极低。

如果信号特别长、VMD 单次分解需要几十秒，建议先用信号的前 0.5 秒做一次小规模预跑，确定一个大致的参数丘陵地形，再把它作为 SCA 初始种群的一部分。这个小技巧能省下大量时间，而且因为 VMD 参数优化对信号的全局统计特性更敏感、对局部波形不太敏感，短信号预跑的结果通常不会跑偏太多。如果你有明确的计算时间预算，就把种群和迭代先减半跑一版，再把结果作为经验解塞进第二轮的全量跑，这比直接跑全参数更实用。

4. 从信号读取到最优参数输出：完整可复现的 SCA-VMD 工程

4.1 工程依赖与最小文件结构

先把环境交代清楚。这个方案在 Python 3.8 以上版本能直接跑通，依赖四个库：numpy做数组运算，scipy提供 Hilbert 变换，vmdpy负责 VMD 分解，tqdm用来显示迭代进度。首次接触 Python 的读者，先确认python --version是否正常，再用虚拟环境装依赖，避免把系统 Python 环境搞乱。

python -m venv sca_vmd_env source sca_vmd_env/bin/activate # Windows 下用 sca_vmd_env\Scripts\activate pip install numpy scipy vmdpy tqdm

vmdpy是一个轻量级的 VMD 封装库，内核直接翻译自原作者论文里的 MATLAB 代码，接口只有一行VMD(f, alpha, tau, K, DC, init, tol)，非常适合被优化算法反复调用。不需要自己去复刻 VMD 的矩阵迭代逻辑，那是另外一条深坑路线，做参数优化只需要把它当黑匣子即可。

4.2 完整主程序：合成信号演示与结果输出

为了让你拿到就能跑，我用一个包含两个正弦分量加白噪声的合成信号做演示。实际应用时把raw = np.loadtxt("vibration.csv")换成你自己的数据源即可。完整代码如下：

import numpy as np from scipy.signal import hilbert from vmdpy import VMD from tqdm import tqdm def envelope_entropy(imf): env = np.abs(hilbert(imf)) p = env / np.sum(env) return -np.sum(p * np.log(p + 1e-12)) def fitness(params, signal, tau=0, DC=0, init=1, tol=1e-7): K = int(np.clip(np.round(params[0]), 2, 12)) alpha = float(np.clip(params[1], 500, 5000)) u, _, _ = VMD(signal, alpha, tau, K, DC, init, tol) return min([envelope_entropy(u[k, :]) for k in range(K)]) def sca_vmd(signal, pop=12, T=60, seed=42): rng = np.random.default_rng(seed) lb, ub = np.array([2.0, 500.0]), np.array([12.0, 5000.0]) X = lb + rng.random((pop, 2)) * (ub - lb) X[0] = np.array([8.0, 2000.0]) fits = np.array([fitness(x, signal) for x in X]) best_pos = X[np.argmin(fits)].copy() best_fit = fits.min() for t in tqdm(range(T), desc="SCA 迭代"): r1 = 2.0 - t * (2.0 / T) for i in range(pop): r2 = rng.uniform(0, 2 * np.pi) r3 = rng.uniform(0, 2) r4 = rng.random() if r4 < 0.5: X[i] = X[i] + r1 * np.sin(r2) * np.abs(r3 * best_pos - X[i]) else: X[i] = X[i] + r1 * np.cos(r2) * np.abs(r3 * best_pos - X[i]) X[i] = np.clip(X[i], lb, ub) for i in range(pop): fits[i] = fitness(X[i], signal) if fits.min() < best_fit: best_idx = np.argmin(fits) best_pos = X[best_idx].copy() best_fit = fits[best_idx] return best_pos, best_fit if __name__ == "__main__": fs = 1000 t = np.linspace(0, 1, fs, endpoint=False) raw = (1.2 * np.sin(2 * np.pi * 50 * t) + 0.6 * np.sin(2 * np.pi * 120 * t) + 0.1 * np.random.randn(fs)) # 去均值并做 z-score 归一化，避免幅值尺度影响 alpha 寻优 raw = (raw - raw.mean()) / raw.std() best_pos, best_fit = sca_vmd(raw, pop=12, T=60) K = int(np.round(best_pos[0])) alpha = round(float(best_pos[1]), 2) print(f"最优 K：{K}，最优 alpha：{alpha}，包络熵：{best_fit:.6f}") # 用最优参数做最终分解并打印中心频率 u, _, omega = VMD(raw, alpha, 0, K, 0, 1, 1e-7) for k in range(K): print(f"模态 {k+1} 中心频率：{omega[-1, k]:.2f} Hz")

参数说明：pop=12表示每个种群有 12 组候选参数，T=60是总迭代次数。raw.std()归一化这一步很容易被忽略，但很重要。实测如果不做归一化，alpha 的寻优结果会受到信号幅值的直接影响，同一段信号放大 10 倍后最优 alpha 可能翻倍，归一化之后结果才具备传递性。VMD函数里的init=1表示均匀初始化中心频率，DC=0表示不在第一个模态里保留直流分量。omega[-1, k]取的是最后一次迭代收敛后的中心频率，这是判断模态是否分离成功的第一手证据。

4.3 结果怎么判定：收敛曲线和中心频率间隔是两条线

程序跑完输出一组最优参数不是终点，还要做两步验证。第一步看适应度收敛曲线，把每一轮的best_fit记录下来画图，如果收敛曲线在最后 10 代还在明显下降，说明 T 设小了，需要加迭代次数。第二步看中心频率间隔，K=5 时如果有一对相邻模态的中心频率差值只有 2Hz，基本可以断定过分解，哪怕包络熵算出来很低也要警惕。低包络熵可能是虚假冲击挣来的，不代表分解有物理意义。

从上述合成信号的运行结果看，SCA 通常能稳定找到 K=2～3、alpha=1000～2000 之间的组合，并且中心频率分别落在 50Hz 和 120Hz 附近。如果你用的是真实采集信号，中心频率一定要对照设备的转频和故障特征频率，这一步直接决定后续的故障诊断结论是否可信。

5. SCA-VMD 落地避坑指南：五个我踩过的常见问题

这一章专门列几个基于真实使用经验的坑，每一条我都按“现象 → 原因 → 解决”的路径写过。

5.1 最优 K 永远顶在边界上：先查范围设置，再怀疑信号

现象：SCA 跑完输出 K=12，alpha=5000，正好落在边界最大值，而且换几组随机种子结果都一样。

原因分两种可能。第一种是参数范围本身定窄了，真实最优可能确实在边界之外；第二种是目标函数在边界区域单调递减，算法一路顺着斜坡滑到边界上，但斜坡尽头不是真正的谷底。处理办法是先放大范围再跑一次，比如 K 上限到 15、alpha 上限到 10000，如果结果仍然顶在边界上，就手动扫一遍 K=2～15 的包络熵曲线，看边界处是不是一条明显下降的斜坡。如果是，说明这个信号本身的合适 K 就很大，比如包含十几个窄带分量，这时候要回到信号本身去确认物理成分数量，而不是继续加大边界。

5.2 每次跑出来的最优 K 都不一样：随机种子不是摆设

现象：同一个信号，上下两次运行，一次 K=4，一次 K=8，包络熵还都差不多。

原因：SCA 本身是随机优化算法，np.random.default_rng没有固定种子时，每轮初始种群和更新方向全不同。而 VMD 参数优化目标函数又带有不少局部极小，不同的随机起点确实会收敛到不同的局部解。我的习惯是固定seed=42作为默认配置，跑三组不同种子交叉验证。如果三组结果都落在同一组参数附近，说明这个解比较可靠；如果差异很大，优先怀疑信号里有明显混叠成分，而不是算法有问题。

5.3 包络熵很低但分解结果明显不合理：评价指标被高频噪声带偏

现象：SCA 找到的 alpha 特别大，分解出的某个模态是一段高频噪声，包络熵却很低。

原因：包络熵对稀疏冲击敏感，如果某段 IMF 正好有尖锐脉冲，即使这个脉冲是没有物理含义的数值毛刺，包络仍然会显得稀疏，熵值大幅下降。这不是 SCA 的问题，是适应度函数选择的问题。解决方法是在fitness里加一个峭度约束，只有峭度大于某个阈值（比如 3.5）的模态才参与包络熵比较，否则直接给一个差评。也可以换成复合目标函数，比如最小包络熵除以峭度，牺牲一点纯净度换取更可靠的适用范围。

5.4 VMD 调用报错却不说清楚原因：先把参数类型和范围过一遍

现象：运行到一半弹TypeError或ValueError，甚至进程直接卡死。

原因：最常见的问题是 K 被传入浮点数，vmdpy里需要整数索引，SCA 的位置更新天然产出浮点值，必须显式int(np.round())；其次是 alpha 被传成 0 或负数，这是因为边界裁剪没有生效；还有一种隐蔽情况是信号长度太短，VMD 的矩阵维度在分解阶段就对不上。解决方式是在fitness入口处做一次防御性检查，打印出实际传入的 K、alpha 和信号长度。这里用错误的参数去调试 VMD 内部逻辑是在浪费时间，先把传入参数打印出来，九成问题当场就能定位。

5.5 中心频率分明叠在一起，SCA 却报成功：目标函数不感知过分解

现象：最优参数的 K 较大时，打印出的中心频率有两个值只差几赫兹，明显是重复模态，但包络熵确实比 K=5 时更低。

原因：包络熵是单模态指标，它完全不检查模态之间的独立性。K 过大时，两个相邻模态可能分别锁定了同一频率成分的两个相位片段，每个片段都比较“干净”，熵值反而下降，这是典型的过分解奖励。解决办法是在适应度函数里附加一个模态分离惩罚项：计算两两中心频率的最短间距，如果小于信号频率分辨率的某个倍数（比如 5 倍），就对该组参数施加加分惩罚。这一步能让 SCA 在搜索过程中避开过分解区域。

6. 让 SCA-VMD 真正进入工作流的三个进阶操作

把 SCA 跑通只是第一步，真正能落地的方案还得做交叉验证、性能优化和下游任务衔接。

先用网格搜索给 SCA 的结论做交叉验证。以当前最优 K 为中心，把 K 从 2 到 15 逐个固定，对每个 K 用 SCA 单独搜 alpha，画出 K-包络熵曲线。如果曲线在 SCA 给出的 K 附近有明显谷底，说明结论可靠；如果曲线是一条几乎平直的线，说明包络熵对这个信号不敏感，这时候加别的指标比继续调 SCA 更有意义。网格搜索的代价是几十次 VMD 分解，放在跑完 SCA 之后做，总耗时完全可以接受。

进一步看运行效率。SCA 的适应度计算完全独立，非常适合并行。用joblib把每代种群适应度计算分发到多进程，实测在四核机器上能提速到接近 3 倍。注意用进程池而不是线程池，因为vmdpy内部是 numpy 运算，存在 GIL 限制，线程并行反而更慢。还要留意内存：每个 VMD 分解都会保留一组u数组，进程数设成 CPU 核数的一半即可，避免大信号下内存涨爆。

最后说下游衔接。优化出来的 K 和 alpha 不是终点，它们的价值体现在后面两步：用最优参数做分解之后，对每个模态算希尔伯特包络谱，把包络谱里的峰值频率与故障特征频率做对照；如果包络谱里优势频率对应轴承外圈故障频率，再算这个频段的能量占比作为健康指标，送入后续的分类模型。SCA-VMD 在这里扮演的是特征提取前置模块，它解决的是“分解参数不再依赖人工拍板”这一环，但下游是否买账，还得看中心频率和包络谱是否长得符合物理机理。

我自己的做法是：在每次做完 SCA-VMD 后，固定把这些输出追加到实验记录里——最优参数、中心频率间隔、包络熵、以及一次手工 K 对比的结果。久而久之就攒出了自己的经验库，再遇到类似信号能直接给出起点参数，二次寻优时间能砍一半。这套脚本我已经在多个轴承数据集上验证过，K 一栏基本不再需要人工修改。希望帮到你。

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