简介一份面向信号处理与数据分析人员的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 设大了同一个成分被劈成两半还会出现一堆没有物理含义的虚假模态。以前我调这两个参数靠手工试凑从 K3 试到 K12每组跑一遍对照效率极低结果还全凭经验。正余弦算法SCA是一种结构很简单的群体优化算法把 K 和 alpha 当成两个连续变量用包络熵做适应度函数让算法自己去搜最优组合。全套逻辑用 Python 实现从读数据到输出最优参数一共不过一百多行代码。这篇笔记适合正在做信号分解、故障诊断或者特征提取的读者照着跑通就能把人工试凑换成自动寻优。2. VMD 参数为什么要优化K 和 alpha 到底在控制什么2.1 变分模态分解的核心参数与取值边界先明确一个事实VMD 不是一劳永逸的分解工具它的效果完全取决于两个超参数。K 是模态个数决定把信号拆成几份alpha 是二次惩罚系数作用在带宽约束上。alpha 越小每个模态的带宽约束越松模态中心频率附近的频率分量更容易被吸收alpha 越大带宽越窄模态越“纯净”但也可能把同一个物理成分截断。工程上最常见的问题恰恰在这里——这两个参数互相牵扯K 变了alpha 的合适区间也跟着变。比如信号里有 50Hz 和 120Hz 两个正弦分量K2 同时把 alpha 从 500 拉到 3000分解出来的两个模态中心频率可能没变化但模态 1 的带宽明显收窄波形畸变更严重。我用过一段时间的固定参数经验表大致是这样的边界效应参数取值方向典型后果适用场景K偏小多个频率分量混叠模态频谱出现多峰频率成分少、分离度高的信号K偏大产生虚假模态相邻模态频率重叠信号成分多但每个都窄带alpha偏小(数百)模态带宽大噪声容易混进模态需要保留细节、容忍噪声alpha偏大(数千)模态整齐但可能截断真实的频率成分强周期成分、冲击特征明显这个表只作为起点参考真正的问题在于每个实际信号的“合适区间”都不同这就是为什么需要算法去搜。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, tau0, DC0, init1, tol1e-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)每个模态一行循环取出来算熵即可。tau0是 VMD 论文里的标准设置表示不考虑噪声容差。2.3 参数边界怎么定别拍脑袋先看信号再做范围适应度函数定了边界的设置是另一个影响结果的细节。我见过很多人直接把 K 定成 [1, 20]alpha 定成 [1, 100000]这种宽得离谱的区间让 SCA 白费大量迭代去探索完全没有意义的参数组合。实际经验是K 的下限低于 2 没有意义因为 VMD 至少应该把信号分成直流分量和主分量两路K 的上限按采样频率和信号长度来粗估对 1 秒 1000 点、成分不超过 4 个的信号K 上限设到 1012 已经足够。alpha 的边界更值得说小于 500 时模态带宽过大分解结果经常出现相邻模态互相“抢频段”大于 5000 时收敛变慢模态中心频率容易漂移出真实物理成分。8003000 是 VMD 论文里举例的常用区间落地时放宽到 5005000 保留寻优空间。还有一个工程经验在 SCA 初始化种群时塞几个经验解进去比如 (K8, alpha2000)剩下的个体随机生成。这比全随机种群收敛快因为正余弦算法本身没有很强的记忆能力一个好的起点能明显减少前期无效探索。种群规模和迭代次数的基本盘是 pop12、T60这个配置单条信号在普通笔记本上大约几十秒到两三分钟。如果想压时间优先降迭代次数而不是降种群因为种群太少容易让 SCA 直接陷进局部最优。3. 正余弦算法的核心逻辑四个随机参数和一个衰减因子3.1 SCA 位置更新公式拆解为什么它能兼顾探索和收敛正余弦算法是 Mirjalili 在 2016 年提出的原理解释起来比粒子群还简单每个候选解的位置更新只参照当前全局最优解用正弦或余弦函数来生成步进方向和步长。具体公式是当随机数 r4 0.5 时位置按正弦更新否则按余弦更新。每次更新里有四个随机数起作用r1 控制步长幅度它在迭代过程中先从 2 线性衰减到 0前期的探索步长大后期逐步收敛r2 是随机角度范围在 [0, 2π]决定移动方向r3 给最优解位置加一个随机权重范围 [0, 2]避免每次移动完全重复r4 则是决定这次用正弦还是余弦分支。这四个参数配合起来让 SCA 前期大范围跳跃、后期逐步逼近最优解算法结构没有任何微分和梯度依赖所以特别适合像 VMD 参数这样黑匣子式的目标函数。我把四参数的作用整理成一张表方便对照代码理解参数取值范围作用对搜索行为的影响r1a - t * a/T (a2)步长控制前期大范围探索后期精细收敛r2[0, 2π] 均匀随机方向控制决定位置朝最优解的外侧还是内侧移动r3[0, 2] 均匀随机最优解加权避免所有个体朝同一点靠拢r4[0, 1] 均匀随机分支选择平衡正弦分支与余弦分支的利用率理解 r1 的衰减最关键。它直接用当前迭代次数 t 和总迭代次数 T 计算a2 是标准配置。t 小的时候 r1 接近 2步长大t 接近 T 的时候 r1 接近 0每个个体在最优解周围微调。这个衰减逻辑很实用但要注意它没有自适应的能力如果 T 设得太小后期还没来得及收敛就停了如果 T 设得太大后半段有大量迭代在做无谓的微调。3.2 SCA 主循环代码从初始化到迭代更新SCA 的实现非常规整我就直接给出一个适用于 VMD 参数寻优的标准循环。为了保证可复现用numpy.random.default_rng管理随机源支持固定随机种子这个习惯能避免“每次运行结果不一样”的尴尬。def sca_vmd(signal, pop12, T60, seed42): 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里存的是连续值比如 K7.32、alpha2134.7要转成可执行的 VMD 参数在调用fitness时已经做了round和clip所以最终输出时要手动处理一遍。3.3 种群和迭代次数的工程经验先跑小规模验证再上全量SCA 的寻优效果对参数不算特别敏感但有一个原则优先保证种群多样性而不是盲目增大迭代次数。种群设 8 的时候我遇到过 SCA 在某个局部解附近反复打转的情况因为个体太少交叉覆盖能力太弱。种群到 1220 之后每增加一个个体的边际收益明显下降但耗时线性增加因为每个个体都要跑一次 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 tqdmvmdpy是一个轻量级的 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, tau0, DC0, init1, tol1e-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, pop12, T60, seed42): 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), descSCA 迭代): 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, endpointFalse) 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, pop12, T60) 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模态 {k1} 中心频率{omega[-1, k]:.2f} Hz)参数说明pop12表示每个种群有 12 组候选参数T60是总迭代次数。raw.std()归一化这一步很容易被忽略但很重要。实测如果不做归一化alpha 的寻优结果会受到信号幅值的直接影响同一段信号放大 10 倍后最优 alpha 可能翻倍归一化之后结果才具备传递性。VMD函数里的init1表示均匀初始化中心频率DC0表示不在第一个模态里保留直流分量。omega[-1, k]取的是最后一次迭代收敛后的中心频率这是判断模态是否分离成功的第一手证据。4.3 结果怎么判定收敛曲线和中心频率间隔是两条线程序跑完输出一组最优参数不是终点还要做两步验证。第一步看适应度收敛曲线把每一轮的best_fit记录下来画图如果收敛曲线在最后 10 代还在明显下降说明 T 设小了需要加迭代次数。第二步看中心频率间隔K5 时如果有一对相邻模态的中心频率差值只有 2Hz基本可以断定过分解哪怕包络熵算出来很低也要警惕。低包络熵可能是虚假冲击挣来的不代表分解有物理意义。从上述合成信号的运行结果看SCA 通常能稳定找到 K23、alpha10002000 之间的组合并且中心频率分别落在 50Hz 和 120Hz 附近。如果你用的是真实采集信号中心频率一定要对照设备的转频和故障特征频率这一步直接决定后续的故障诊断结论是否可信。5. SCA-VMD 落地避坑指南五个我踩过的常见问题这一章专门列几个基于真实使用经验的坑每一条我都按“现象 → 原因 → 解决”的路径写过。5.1 最优 K 永远顶在边界上先查范围设置再怀疑信号现象SCA 跑完输出 K12alpha5000正好落在边界最大值而且换几组随机种子结果都一样。原因分两种可能。第一种是参数范围本身定窄了真实最优可能确实在边界之外第二种是目标函数在边界区域单调递减算法一路顺着斜坡滑到边界上但斜坡尽头不是真正的谷底。处理办法是先放大范围再跑一次比如 K 上限到 15、alpha 上限到 10000如果结果仍然顶在边界上就手动扫一遍 K215 的包络熵曲线看边界处是不是一条明显下降的斜坡。如果是说明这个信号本身的合适 K 就很大比如包含十几个窄带分量这时候要回到信号本身去确认物理成分数量而不是继续加大边界。5.2 每次跑出来的最优 K 都不一样随机种子不是摆设现象同一个信号上下两次运行一次 K4一次 K8包络熵还都差不多。原因SCA 本身是随机优化算法np.random.default_rng没有固定种子时每轮初始种群和更新方向全不同。而 VMD 参数优化目标函数又带有不少局部极小不同的随机起点确实会收敛到不同的局部解。我的习惯是固定seed42作为默认配置跑三组不同种子交叉验证。如果三组结果都落在同一组参数附近说明这个解比较可靠如果差异很大优先怀疑信号里有明显混叠成分而不是算法有问题。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 较大时打印出的中心频率有两个值只差几赫兹明显是重复模态但包络熵确实比 K5 时更低。原因包络熵是单模态指标它完全不检查模态之间的独立性。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 一栏基本不再需要人工修改。希望帮到你。本文还有配套的精品资源点击获取