简介面向机械故障诊断与信号处理研究者的MATLAB工具包提供基于灰狼优化算法GWO对变分模态分解VMD参数进行智能寻优的完整实现可有效解决VMD分解中惩罚因子与模态个数依赖人工经验设定的难题。程序内置四种适应度函数通过criterion参数灵活切换分别以排列熵、最小包络熵、信息熵或样本熵最小化为优化目标适应不同信号特征下的分解需求。压缩包共16个文件包括8个M脚本、5个Excel数据表、2个MAT数据文件和1张算法流程图主程序、目标函数、绘图模块一应俱全数据文件与表格可直接用于复现和结果对比。资源大小仅6.36MB已有322人学习适合具备一定MATLAB基础、希望优化VMD参数并快速开展实验验证的算法工程师与科研人员使用。1. VMD 分解的 K 和 alpha 定不住灰狼算法是现成的解法VMD 变分模态分解在 MATLAB 里实现并不难但真正卡住人的不是分解本身而是参数分解层数 K 和惩罚因子 alpha 怎么定。这两项直接决定结果是“过度细分”还是“模态混叠”工程信号里几乎没有一组固定的经验值能通吃——轴承故障信号好用的参数换到齿轮箱振动信号往往就失效。灰狼优化算法 GWO 解决的就是这个痛点用包络熵或排列熵作目标函数自动迭代搜索 K、alpha、tau 的最优组合把“试参数”变成可复现的自动化流程。对常年做故障诊断、振动分析或信号处理的人来说这套东西的价值在于它不要求你事先知道信号里到底有几个分量。本文从 GWO 为什么适合 VMD 讲起给出可在 MATLAB R2023b 上直接套用的 GWO-VMD 程序结构、核心代码和调参思路最后聊几个实际部署时容易踩的坑。2. VMD 为何需要 GWOK 和 alpha 靠手调不靠谱2.1 VMD 的核心参数与物理含义VMD 把信号分解问题写成变分模型寻找 K 个模态函数使所有模态的带宽之和最小同时这些模态之和能精确重构原信号。迭代求解时alpha 是带宽惩罚因子alpha 越大各模态带宽越窄、频带越干净但过大可能丢失瞬态冲击的能量alpha 太小相邻模态之间就会发生频谱重叠。除 K 和 alpha 外tau 控制拉格朗日乘子的更新步长DC1 表示第一个模态包含直流分量init1 表示初始中心频率为均匀分布。实际信号分析中K 是最让人头疼的参数。K 设小了多个频率成分叠在同一个模态里K 设大了又会出现相邻中心频率几乎重合的空模态或碎片模态。理论上可以通过观察中心频率是否聚集来判断 K 是否合适但当信号包含噪声与多个边带时肉眼判断很不靠谱。所以 GWO-VMD 的思路很直接把 K、alpha、tau 一起放进一个优化目标里用智能搜索代替手工试探。2.2 为什么选灰狼优化算法而不是网格搜索或 PSO网格搜索的第一问题是维度爆炸K 取 210、alpha 取 200~2000 步长 100组合数就是 9*19171 次完整 VMD 分解每次分解还要内部迭代几十轮在 MATLAB 里跑完一轮要几十分钟。第二问题是它无法感知连续性alpha 取 990 和 1000 的分解结果差异可能微乎其微网格却把它们当成两个完全独立的对象。灰狼优化算法的优势在于参数少、无需梯度。GWO 只依赖 a、A、C 三个系数完成对猎物的包围、狩猎和攻击模拟位置更新公式简单MATLAB 几十行就能实现。与粒子群 PSO 相比GWO 没有个体速度项不容易因速度过大而飞出边界与遗传算法相比GWO 没有编码解码过程直接以连续实数操作。在 VMD 参数寻优这类低维连续问题上GWO 收敛快且稳定。更关键的是GWO 天然支持混合整数——K 是整数alpha 和 tau 是连续量只需要在目标函数里对 K 做一次 round 取整即可。2.3 目标函数选包络熵还是排列熵GWO 的搜索方向完全由目标函数决定。包络熵度量模态包络的不确定性滚动轴承或齿轮产生冲击时包络信号中会出现稀疏的脉冲峰这种信号的包络熵较小如果模态被过度分解或混入噪声包络会变得杂乱包络熵变大。因此最小化包络熵可以得到“包络最稀疏、冲击特征最明显”的模态。排列熵则强调模态的时间序复杂性适合含有强随机噪声、需要抑制随机分量的场景。实际项目中我一般默认选包络熵因为它的物理意义明确、计算量小只有当信号本身没有明显冲击特征时才改用排列熵。下表给出常用目标函数的对比目标函数计算成本适用场景方向包络熵低轴承、齿轮等冲击型故障信号最小值排列熵中强噪声背景下的随机性分析最小值信息熵低普遍适用但灵敏度较低最小值模态混叠度高需要显式抑制中心频率重叠最小值包络熵的核心是希尔伯特变换后取模得到包络信号再求包络的能量分布熵值。它天然关注冲击成分的稀疏性这和大多数旋转机械故障特征高度吻合。所以目标函数这块没有太多悬念先用包络熵把 GWO 跑通再按需替换是性价比最高的路径。3. GWO-VMD 的 MATLAB 程序结构与核心代码实现3.1 顶层主脚本传数据、设边界、启动寻优GWO-VMD 的程序结构并不复杂主脚本负责加载信号和设置搜索边界调用 gwoVmd 寻优函数最后用最优参数执行一次完整 VMD 分解。注意rng(42) 这一行不能省因为 GWO 使用随机初始化灰狼种群固定随机种子才能保证可复现。%% main_gwo_vmd.m clc; clear; close all; rng(42); % 加载采样率为 fs 的信号 x列向量 load(bearing_fault.mat, x, fs); % GWO 参数设置 SearchAgents_no 25; % 灰狼数量推荐区间 20~30 Max_iter 30; % 最大迭代次数推荐区间 20~50 dim 3; % 待优化参数数量K, alpha, tau % 搜索边界行向量顺序与 dim 对应 lb [2, 200, 0]; % K 最小 2 层alpha 最小 200tau 最小 0 ub [10, 2000, 0.3]; % K 最大 10 层alpha 最大 2000tau 最大 0.3 % 调用 GWO-VMD 寻优 [bestPos, bestFitness, convergence] gwoVmd(... x, fs, SearchAgents_no, Max_iter, lb, ub, dim); % 输出最优参数并做最终分解 K_opt round(bestPos(1)); % VMD 层数必须是整数 alpha_opt bestPos(2); tau_opt bestPos(3); [u, u_hat, omega] VMD(x, alpha_opt, tau_opt, K_opt, 0, 1, 1e-7); fprintf(最优参数 K%d, alpha%.2f, tau%.4f\n, ... K_opt, alpha_opt, tau_opt);主脚本逻辑很直白但有两个易错点。一是 bestPos(1) 必须经过 round 取整之后才能传给 VMD否则会报维度错误。二是搜索边界并不是越大越好。K 上界设 10 是因为实际工程信号通常不会超过 10 个有效模态再多容易出现相邻中心频率合并的现象alpha 上界设 2000 是避免模态带宽过窄导致一个冲击被拆成多个窄带分量反而增加包络熵。tau 的搜索区间很小如果为了减少搜索维度可以固定 tau0把 dim 改为 2寻优速度会提升约三分之一。3.2 gwoVmd 核心函数包围、狩猎与攻击三阶段的位置更新灰狼优化算法的核心代码要处理三件事维护三只头狼Alpha、Beta、Delta的最优位置、用包围公式更新所有个体位置、线性衰减控制参数 a。下面是可直接运行的 gwoVmd 函数内部调用 vmdObjective 目标函数。function [bestPos, bestFitness, convergence] gwoVmd(... signal, fs, nAgents, maxIter, lb, ub, dim) % 灰狼优化算法搜索 VMD 最优参数 % signal: 待分解信号; fs: 采样率 % nAgents: 灰狼个数; maxIter: 最大迭代次数 % lb, ub: 参数下界与上界; dim: 优化维度 % 随机初始化灰狼种群 Positions zeros(nAgents, dim); for i 1:nAgents Positions(i, :) lb rand(1, dim) .* (ub - lb); end Alpha_pos zeros(1, dim); Alpha_score inf; Beta_pos zeros(1, dim); Beta_score inf; Delta_pos zeros(1, dim); Delta_score inf; convergence zeros(1, maxIter); for t 1:maxIter % a 从 2 线性衰减到 0控制探索与开发 a 2 - t * (2 / maxIter); for i 1:nAgents % 越界修正越界个体拉回边界 Positions(i, :) max(Positions(i, :), lb); Positions(i, :) min(Positions(i, :), ub); % 计算当前灰狼位置的适应度 fitness vmdObjective(signal, fs, Positions(i, :)); % 更新 Alpha、Beta、Delta 三只头狼 if fitness Alpha_score Alpha_score fitness; Alpha_pos Positions(i, :); elseif fitness Beta_score Beta_score fitness; Beta_pos Positions(i, :); elseif fitness Delta_score Delta_score fitness; Delta_pos Positions(i, :); end end % 更新所有灰狼的位置 for i 1:nAgents for j 1:dim r1 rand; r2 rand; A1 2*a*r1 - a; C1 2*r2; D_alpha abs(C1 * Alpha_pos(j) - Positions(i, j)); X1 Alpha_pos(j) - A1 * D_alpha; r1 rand; r2 rand; A2 2*a*r1 - a; C2 2*r2; D_beta abs(C2 * Beta_pos(j) - Positions(i, j)); X2 Beta_pos(j) - A2 * D_beta; r1 rand; r2 rand; A3 2*a*r1 - a; C3 2*r2; D_delta abs(C3 * Delta_pos(j) - Positions(i, j)); X3 Delta_pos(j) - A3 * D_delta; Positions(i, j) (X1 X2 X3) / 3; end end convergence(t) Alpha_score; end bestPos Alpha_pos; bestFitness Alpha_score; end代码里最关键的是 A 和 C 两个系数。A 的绝对值大于 1 时灰狼远离猎物对应算法的探索阶段A 绝对值小于 1 时逼近猎物对应开发阶段。a 从 2 线性衰减到 0使前期大范围搜索、后期局部精细收敛。C 是随机扰动项它在 [0,2] 之间波动作用是让狼群在接近最优解时不会完全停滞避免陷入局部极小值。三个头狼同时引导位置更新比单头狼引导的 PSO 有更好的种群多样性。越界修正放在适应度计算之前是因为 VMD 对负数参数会直接报错。alpha 一旦越界变成负数VMD 内部的带宽惩罚项就失去物理意义程序直接崩溃。有些人会在目标函数里做边界判断我却建议在主函数里处理这样可以减少 VMD 函数调用出错的机会。3.3 目标函数 vmdObjective每一次适应度评估都是一次完整 VMD目标函数是 GWO 和 VMD 之间的桥梁。编码上有点容易绕VMD 分解得到所有模态对每个模态算包络熵然后取均值作为适应度。为什么不取最小或最大单个模态的包络熵可能因为噪声出现极值均值更能反映整体分解质量。function fitness vmdObjective(signal, fs, param) % 计算 VMD 分解后的平均包络熵 % param [K, alpha, tau]K 为实数需取整 K round(param(1)); alpha param(2); tau param(3); if K 2 || alpha 0 fitness 1e6; % 非法参数直接给大惩罚值 return; end [u, ~, ~] VMD(signal, alpha, tau, K, 0, 1, 1e-7); [n, K] size(u); entropySum 0; for k 1:K % 希尔伯特变换得到解析信号取模得到包络 analytic hilbert(u(:, k)); envelope abs(analytic); p envelope ./ sum(envelope); % 包络熵公式加 eps 防止 log(0) entropySum entropySum - sum(p .* log(p eps)); end fitness entropySum / K; end这段代码有几个值得注意的细节。hilbert 是 MATLAB 信号处理工具箱自带函数执行一次希尔伯特变换后取模就得到包络。p envelope / sum(envelope) 是对包络做归一化使其满足概率分布和为 1 的条件。log 里加 eps 是为了防止包络为零时出现无穷大。适应度取所有模态的平均包络熵可以让搜索同时兼顾所有模态的质量而不是只盯着某一个分量。这里考虑一个实际问题GWO 在迭代过程中会反复调用 vmdObjective每次都要执行一次完整的 VMD 分解。如果信号长度为 10 万点灰狼数量 25迭代 30 次总共要执行 750 次 VMD。因此信号较长时建议先降采样到 2048 或 4096 点做寻优锁定最优参数后再用全采样率做最终分解。这个技巧能省下大量时间且不会对参数寻优结果造成明显影响。3.4 VMD 函数ADMM 迭代与五个参数的作用VMD 内核函数是整个流程的底层引擎它的完整实现有约 100 行代码核心是交替方向乘子法 ADMM 迭代。主循环中先用维纳滤波更新模态谱 u_hat再更新中心频率 omega最后更新拉格朗日乘子 lambda_hat。考虑到篇幅这里给出 VMD 函数的接口和各参数含义完整代码可直接参考原始论文仿写。function [u, u_hat, omega] VMD(signal, alpha, tau, K, DC, init, tol) % VMD 变分模态分解 % signal: 输入信号, 列向量 % alpha: 带宽惩罚因子, 越大模态带宽越窄 % tau: 拉格朗日乘子更新步长, 0 表示噪声为零 % K: 模态分解层数 % DC: 第一模态是否包含直流分量, 1 表示包含 % init: 初始中心频率方式, 1均匀分布, 0全零 % tol: 收敛精度 % u: 分解后的模态信号, 每列一个模态 % u_hat: 模态的频域表示 % omega: 各模态的中心频率 N length(signal); f (0:N-1) / N * 2 * pi; % 频率轴 % 镜像延拓处理边界效应 signal_mirror [signal(end:-1:2); signal; signal(end-1:-1:1)]; % --- ADMM 迭代主体约 60 行此处省略 --- endVMD 的五个参数里注意 DC 参数对第一模态的影响。当 DC1 时第一个模态允许包含直流分量适合处理含有明显偏置的传感器信号DC0 则强制所有模态围绕各自的中心频率带限分布。init1 是均匀初始化保证每次分解结果可复现改成 init0 从零开始迭代会引入随机性GWO 寻优时最好别用 init0否则同一组参数两次分解的包络熵会不一样直接影响适应度对比。4. 参数设置、收敛判据与运行提速的实用技巧4.1 搜索边界和 GWO 参数怎么定一张表说清推荐值GWO-VMD 运行前要确定的参数包括 GWO 自身的参数和 VMD 参数的搜索范围。下表汇总了我在电力谐波、轴承故障、齿轮箱振动三种典型信号上使用的推荐值参数搜索范围/取值推荐值设定理由K[2, 10]4~8K10 后模态中心频率聚集出现空模态alpha[200, 2000]1000 附近过大过小都会让包络熵增大tau[0, 0.3]0 或 0.01噪声小时取 0噪声大时取 0.01SearchAgents_no20~3025太少容易早熟太多耗时长Max_iter20~5030从收敛曲线观察是否充分收敛VMD 内部 tol1e-6 ~ 1e-81e-7精度足够再小影响不大关于 alpha 边界的设定我多说一句。alpha 太小模态带宽过宽相邻模态中心频率差很小的时候会直接混叠alpha 太大模态被压缩成极窄的谱线瞬态冲击的能量被分散到多个模态里。从包络熵角度看这两种情况都会让包络波形变复杂熵值升高搜索算法会自动避开。因此边界只要设置在一个合理范围就行不必追求完全贴合信号GWO 的连续搜索能力本来就能找到范围内的最优。4.2 从收敛曲线判断寻优是否正常收敛曲线 convergence 是判断寻优质量的最直接依据。把收敛曲线画出来主要有三种情况figure; plot(1:length(convergence), convergence, b-, LineWidth, 1.5); xlabel(迭代次数); ylabel(平均包络熵); title(GWO 收敛曲线); grid on;第一种是曲线单调下降后趋于平缓这是最理想的状态说明狼群逐步逼近最优解。第二种是前几代快速下降之后长时间不走对应的是算法陷入局部最优此时可以把 SearchAgents_no 从 25 增至 40或者把 lb/ub 范围缩小到已有最优解附近重新搜索相当于在局部区域做精细化寻优。第三种是曲线震荡无法收敛比如适应度值忽高忽低这往往是信号本身具有强随机性、或者 VMD 内部分解不稳定的信号可以用多次分解取平均包络熵的方式来平滑适应度。判断收敛的另一个做法是记录每代的 Alpha_pos如果连续 5 代 K 保持不变且 alpha 变化小于 5%基本可以提前终止迭代。这个早停策略在实际使用中可以省下约三分之一的计算时间且不会损失精度。4.3 用 parfor 并行加速多核 CPU 上的直接收益与边界GWO-VMD 的耗时集中在目标函数反复执行 VMD 分解。灰狼种群中每个个体的适应度计算相互独立完全可以用 parfor 做并行加速。修改方式只有一处把 gwoVmd 中计算适应度的 for 循环换成 parfor并把目标函数写成独立 m 文件。% gwoVmd 函数中 fitnessList zeros(SearchAgents_no, 1); parfor i 1:SearchAgents_no fitnessList(i) vmdObjective(signal, fs, Positions(i, :)); end直接换 parfor 有一个前提vmdObjective 必须能被 workers 访问即它必须是一个独立的 m 文件不能定义为嵌套函数signal 会被复制到每个 worker 上内存占用约等于 worker 数乘以信号长度。当信号超过 10 万点时8 个 worker 同时复制信号可能撑爆内存这种情况建议保留普通 for 循环或者把信号先做降采样再用 parfor 寻优。另注意parfor 内不能使用断点或 disp 输出调试时要切回 for。5. 用 GWO-VMD 跑完后的验证消融对比与包络谱分析GWO-VMD 输出最优参数后必须做的一步是验证分解结果的物理有效性。方法很简单用三组不同参数分别做 VMD再比较包络谱中的故障特征频率幅值。第一组是 GWO 寻优得到的最优参数第二组是经验参数比如 K6、alpha1000第三组是随机参数比如 K4、alpha500。分别用三组参数分解信号取第一个模态做希尔伯特包络谱观察特征频率处幅值。如果 GWO 参数组的特征频率峰值最高、噪声底较低说明寻优有效如果三组结果差距不大说明这个信号本身对参数不敏感直接使用 GWO 参数即可不必每次寻优。% 用最优参数分解 [u_best, ~, omega_best] VMD(x, alpha_opt, tau_opt, K_opt, 0, 1, 1e-7); % 用经验参数分解 [u_exp, ~, omega_exp] VMD(x, 1000, 0, 6, 0, 1, 1e-7); % 计算最优参数下第一模态的包络谱 analytic hilbert(u_best(:, 1)); env abs(analytic); N length(env); f_axis (0:N-1) * fs / N; env_spectrum abs(fft(env)); % 找前三个峰值对应的频率 [pks, locs] findpeaks(env_spectrum, f_axis, ... SortStr, descend, NPeaks, 3); disp(特征频率候选:); disp(locs(1:3));执行这段代码后如果 locs 中的频率与理论特征频率如外圈故障频率 BPFO接近说明分解结果可用如果峰值不明显或频率对不上就要回到目标函数把包络熵换成排列熵试试。还有一个技巧可以确认 K 是否选多检查最优参数对应的中心频率 omega计算相邻 omega 的差值。如果存在两个 omega 非常接近的模态如差值小于频域分辨率的 3 倍说明 K 偏大可以手动降低 K 上界后重新寻优。这个做法在 GWO-VMD 的实际调试中非常管用它能把纯数值寻优的结果与信号处理的物理直觉结合起来。本文还有配套的精品资源点击获取
灰狼算法优化VMD参数：MATLAB实现与故障诊断应用
简介面向机械故障诊断与信号处理研究者的MATLAB工具包提供基于灰狼优化算法GWO对变分模态分解VMD参数进行智能寻优的完整实现可有效解决VMD分解中惩罚因子与模态个数依赖人工经验设定的难题。程序内置四种适应度函数通过criterion参数灵活切换分别以排列熵、最小包络熵、信息熵或样本熵最小化为优化目标适应不同信号特征下的分解需求。压缩包共16个文件包括8个M脚本、5个Excel数据表、2个MAT数据文件和1张算法流程图主程序、目标函数、绘图模块一应俱全数据文件与表格可直接用于复现和结果对比。资源大小仅6.36MB已有322人学习适合具备一定MATLAB基础、希望优化VMD参数并快速开展实验验证的算法工程师与科研人员使用。1. VMD 分解的 K 和 alpha 定不住灰狼算法是现成的解法VMD 变分模态分解在 MATLAB 里实现并不难但真正卡住人的不是分解本身而是参数分解层数 K 和惩罚因子 alpha 怎么定。这两项直接决定结果是“过度细分”还是“模态混叠”工程信号里几乎没有一组固定的经验值能通吃——轴承故障信号好用的参数换到齿轮箱振动信号往往就失效。灰狼优化算法 GWO 解决的就是这个痛点用包络熵或排列熵作目标函数自动迭代搜索 K、alpha、tau 的最优组合把“试参数”变成可复现的自动化流程。对常年做故障诊断、振动分析或信号处理的人来说这套东西的价值在于它不要求你事先知道信号里到底有几个分量。本文从 GWO 为什么适合 VMD 讲起给出可在 MATLAB R2023b 上直接套用的 GWO-VMD 程序结构、核心代码和调参思路最后聊几个实际部署时容易踩的坑。2. VMD 为何需要 GWOK 和 alpha 靠手调不靠谱2.1 VMD 的核心参数与物理含义VMD 把信号分解问题写成变分模型寻找 K 个模态函数使所有模态的带宽之和最小同时这些模态之和能精确重构原信号。迭代求解时alpha 是带宽惩罚因子alpha 越大各模态带宽越窄、频带越干净但过大可能丢失瞬态冲击的能量alpha 太小相邻模态之间就会发生频谱重叠。除 K 和 alpha 外tau 控制拉格朗日乘子的更新步长DC1 表示第一个模态包含直流分量init1 表示初始中心频率为均匀分布。实际信号分析中K 是最让人头疼的参数。K 设小了多个频率成分叠在同一个模态里K 设大了又会出现相邻中心频率几乎重合的空模态或碎片模态。理论上可以通过观察中心频率是否聚集来判断 K 是否合适但当信号包含噪声与多个边带时肉眼判断很不靠谱。所以 GWO-VMD 的思路很直接把 K、alpha、tau 一起放进一个优化目标里用智能搜索代替手工试探。2.2 为什么选灰狼优化算法而不是网格搜索或 PSO网格搜索的第一问题是维度爆炸K 取 210、alpha 取 200~2000 步长 100组合数就是 9*19171 次完整 VMD 分解每次分解还要内部迭代几十轮在 MATLAB 里跑完一轮要几十分钟。第二问题是它无法感知连续性alpha 取 990 和 1000 的分解结果差异可能微乎其微网格却把它们当成两个完全独立的对象。灰狼优化算法的优势在于参数少、无需梯度。GWO 只依赖 a、A、C 三个系数完成对猎物的包围、狩猎和攻击模拟位置更新公式简单MATLAB 几十行就能实现。与粒子群 PSO 相比GWO 没有个体速度项不容易因速度过大而飞出边界与遗传算法相比GWO 没有编码解码过程直接以连续实数操作。在 VMD 参数寻优这类低维连续问题上GWO 收敛快且稳定。更关键的是GWO 天然支持混合整数——K 是整数alpha 和 tau 是连续量只需要在目标函数里对 K 做一次 round 取整即可。2.3 目标函数选包络熵还是排列熵GWO 的搜索方向完全由目标函数决定。包络熵度量模态包络的不确定性滚动轴承或齿轮产生冲击时包络信号中会出现稀疏的脉冲峰这种信号的包络熵较小如果模态被过度分解或混入噪声包络会变得杂乱包络熵变大。因此最小化包络熵可以得到“包络最稀疏、冲击特征最明显”的模态。排列熵则强调模态的时间序复杂性适合含有强随机噪声、需要抑制随机分量的场景。实际项目中我一般默认选包络熵因为它的物理意义明确、计算量小只有当信号本身没有明显冲击特征时才改用排列熵。下表给出常用目标函数的对比目标函数计算成本适用场景方向包络熵低轴承、齿轮等冲击型故障信号最小值排列熵中强噪声背景下的随机性分析最小值信息熵低普遍适用但灵敏度较低最小值模态混叠度高需要显式抑制中心频率重叠最小值包络熵的核心是希尔伯特变换后取模得到包络信号再求包络的能量分布熵值。它天然关注冲击成分的稀疏性这和大多数旋转机械故障特征高度吻合。所以目标函数这块没有太多悬念先用包络熵把 GWO 跑通再按需替换是性价比最高的路径。3. GWO-VMD 的 MATLAB 程序结构与核心代码实现3.1 顶层主脚本传数据、设边界、启动寻优GWO-VMD 的程序结构并不复杂主脚本负责加载信号和设置搜索边界调用 gwoVmd 寻优函数最后用最优参数执行一次完整 VMD 分解。注意rng(42) 这一行不能省因为 GWO 使用随机初始化灰狼种群固定随机种子才能保证可复现。%% main_gwo_vmd.m clc; clear; close all; rng(42); % 加载采样率为 fs 的信号 x列向量 load(bearing_fault.mat, x, fs); % GWO 参数设置 SearchAgents_no 25; % 灰狼数量推荐区间 20~30 Max_iter 30; % 最大迭代次数推荐区间 20~50 dim 3; % 待优化参数数量K, alpha, tau % 搜索边界行向量顺序与 dim 对应 lb [2, 200, 0]; % K 最小 2 层alpha 最小 200tau 最小 0 ub [10, 2000, 0.3]; % K 最大 10 层alpha 最大 2000tau 最大 0.3 % 调用 GWO-VMD 寻优 [bestPos, bestFitness, convergence] gwoVmd(... x, fs, SearchAgents_no, Max_iter, lb, ub, dim); % 输出最优参数并做最终分解 K_opt round(bestPos(1)); % VMD 层数必须是整数 alpha_opt bestPos(2); tau_opt bestPos(3); [u, u_hat, omega] VMD(x, alpha_opt, tau_opt, K_opt, 0, 1, 1e-7); fprintf(最优参数 K%d, alpha%.2f, tau%.4f\n, ... K_opt, alpha_opt, tau_opt);主脚本逻辑很直白但有两个易错点。一是 bestPos(1) 必须经过 round 取整之后才能传给 VMD否则会报维度错误。二是搜索边界并不是越大越好。K 上界设 10 是因为实际工程信号通常不会超过 10 个有效模态再多容易出现相邻中心频率合并的现象alpha 上界设 2000 是避免模态带宽过窄导致一个冲击被拆成多个窄带分量反而增加包络熵。tau 的搜索区间很小如果为了减少搜索维度可以固定 tau0把 dim 改为 2寻优速度会提升约三分之一。3.2 gwoVmd 核心函数包围、狩猎与攻击三阶段的位置更新灰狼优化算法的核心代码要处理三件事维护三只头狼Alpha、Beta、Delta的最优位置、用包围公式更新所有个体位置、线性衰减控制参数 a。下面是可直接运行的 gwoVmd 函数内部调用 vmdObjective 目标函数。function [bestPos, bestFitness, convergence] gwoVmd(... signal, fs, nAgents, maxIter, lb, ub, dim) % 灰狼优化算法搜索 VMD 最优参数 % signal: 待分解信号; fs: 采样率 % nAgents: 灰狼个数; maxIter: 最大迭代次数 % lb, ub: 参数下界与上界; dim: 优化维度 % 随机初始化灰狼种群 Positions zeros(nAgents, dim); for i 1:nAgents Positions(i, :) lb rand(1, dim) .* (ub - lb); end Alpha_pos zeros(1, dim); Alpha_score inf; Beta_pos zeros(1, dim); Beta_score inf; Delta_pos zeros(1, dim); Delta_score inf; convergence zeros(1, maxIter); for t 1:maxIter % a 从 2 线性衰减到 0控制探索与开发 a 2 - t * (2 / maxIter); for i 1:nAgents % 越界修正越界个体拉回边界 Positions(i, :) max(Positions(i, :), lb); Positions(i, :) min(Positions(i, :), ub); % 计算当前灰狼位置的适应度 fitness vmdObjective(signal, fs, Positions(i, :)); % 更新 Alpha、Beta、Delta 三只头狼 if fitness Alpha_score Alpha_score fitness; Alpha_pos Positions(i, :); elseif fitness Beta_score Beta_score fitness; Beta_pos Positions(i, :); elseif fitness Delta_score Delta_score fitness; Delta_pos Positions(i, :); end end % 更新所有灰狼的位置 for i 1:nAgents for j 1:dim r1 rand; r2 rand; A1 2*a*r1 - a; C1 2*r2; D_alpha abs(C1 * Alpha_pos(j) - Positions(i, j)); X1 Alpha_pos(j) - A1 * D_alpha; r1 rand; r2 rand; A2 2*a*r1 - a; C2 2*r2; D_beta abs(C2 * Beta_pos(j) - Positions(i, j)); X2 Beta_pos(j) - A2 * D_beta; r1 rand; r2 rand; A3 2*a*r1 - a; C3 2*r2; D_delta abs(C3 * Delta_pos(j) - Positions(i, j)); X3 Delta_pos(j) - A3 * D_delta; Positions(i, j) (X1 X2 X3) / 3; end end convergence(t) Alpha_score; end bestPos Alpha_pos; bestFitness Alpha_score; end代码里最关键的是 A 和 C 两个系数。A 的绝对值大于 1 时灰狼远离猎物对应算法的探索阶段A 绝对值小于 1 时逼近猎物对应开发阶段。a 从 2 线性衰减到 0使前期大范围搜索、后期局部精细收敛。C 是随机扰动项它在 [0,2] 之间波动作用是让狼群在接近最优解时不会完全停滞避免陷入局部极小值。三个头狼同时引导位置更新比单头狼引导的 PSO 有更好的种群多样性。越界修正放在适应度计算之前是因为 VMD 对负数参数会直接报错。alpha 一旦越界变成负数VMD 内部的带宽惩罚项就失去物理意义程序直接崩溃。有些人会在目标函数里做边界判断我却建议在主函数里处理这样可以减少 VMD 函数调用出错的机会。3.3 目标函数 vmdObjective每一次适应度评估都是一次完整 VMD目标函数是 GWO 和 VMD 之间的桥梁。编码上有点容易绕VMD 分解得到所有模态对每个模态算包络熵然后取均值作为适应度。为什么不取最小或最大单个模态的包络熵可能因为噪声出现极值均值更能反映整体分解质量。function fitness vmdObjective(signal, fs, param) % 计算 VMD 分解后的平均包络熵 % param [K, alpha, tau]K 为实数需取整 K round(param(1)); alpha param(2); tau param(3); if K 2 || alpha 0 fitness 1e6; % 非法参数直接给大惩罚值 return; end [u, ~, ~] VMD(signal, alpha, tau, K, 0, 1, 1e-7); [n, K] size(u); entropySum 0; for k 1:K % 希尔伯特变换得到解析信号取模得到包络 analytic hilbert(u(:, k)); envelope abs(analytic); p envelope ./ sum(envelope); % 包络熵公式加 eps 防止 log(0) entropySum entropySum - sum(p .* log(p eps)); end fitness entropySum / K; end这段代码有几个值得注意的细节。hilbert 是 MATLAB 信号处理工具箱自带函数执行一次希尔伯特变换后取模就得到包络。p envelope / sum(envelope) 是对包络做归一化使其满足概率分布和为 1 的条件。log 里加 eps 是为了防止包络为零时出现无穷大。适应度取所有模态的平均包络熵可以让搜索同时兼顾所有模态的质量而不是只盯着某一个分量。这里考虑一个实际问题GWO 在迭代过程中会反复调用 vmdObjective每次都要执行一次完整的 VMD 分解。如果信号长度为 10 万点灰狼数量 25迭代 30 次总共要执行 750 次 VMD。因此信号较长时建议先降采样到 2048 或 4096 点做寻优锁定最优参数后再用全采样率做最终分解。这个技巧能省下大量时间且不会对参数寻优结果造成明显影响。3.4 VMD 函数ADMM 迭代与五个参数的作用VMD 内核函数是整个流程的底层引擎它的完整实现有约 100 行代码核心是交替方向乘子法 ADMM 迭代。主循环中先用维纳滤波更新模态谱 u_hat再更新中心频率 omega最后更新拉格朗日乘子 lambda_hat。考虑到篇幅这里给出 VMD 函数的接口和各参数含义完整代码可直接参考原始论文仿写。function [u, u_hat, omega] VMD(signal, alpha, tau, K, DC, init, tol) % VMD 变分模态分解 % signal: 输入信号, 列向量 % alpha: 带宽惩罚因子, 越大模态带宽越窄 % tau: 拉格朗日乘子更新步长, 0 表示噪声为零 % K: 模态分解层数 % DC: 第一模态是否包含直流分量, 1 表示包含 % init: 初始中心频率方式, 1均匀分布, 0全零 % tol: 收敛精度 % u: 分解后的模态信号, 每列一个模态 % u_hat: 模态的频域表示 % omega: 各模态的中心频率 N length(signal); f (0:N-1) / N * 2 * pi; % 频率轴 % 镜像延拓处理边界效应 signal_mirror [signal(end:-1:2); signal; signal(end-1:-1:1)]; % --- ADMM 迭代主体约 60 行此处省略 --- endVMD 的五个参数里注意 DC 参数对第一模态的影响。当 DC1 时第一个模态允许包含直流分量适合处理含有明显偏置的传感器信号DC0 则强制所有模态围绕各自的中心频率带限分布。init1 是均匀初始化保证每次分解结果可复现改成 init0 从零开始迭代会引入随机性GWO 寻优时最好别用 init0否则同一组参数两次分解的包络熵会不一样直接影响适应度对比。4. 参数设置、收敛判据与运行提速的实用技巧4.1 搜索边界和 GWO 参数怎么定一张表说清推荐值GWO-VMD 运行前要确定的参数包括 GWO 自身的参数和 VMD 参数的搜索范围。下表汇总了我在电力谐波、轴承故障、齿轮箱振动三种典型信号上使用的推荐值参数搜索范围/取值推荐值设定理由K[2, 10]4~8K10 后模态中心频率聚集出现空模态alpha[200, 2000]1000 附近过大过小都会让包络熵增大tau[0, 0.3]0 或 0.01噪声小时取 0噪声大时取 0.01SearchAgents_no20~3025太少容易早熟太多耗时长Max_iter20~5030从收敛曲线观察是否充分收敛VMD 内部 tol1e-6 ~ 1e-81e-7精度足够再小影响不大关于 alpha 边界的设定我多说一句。alpha 太小模态带宽过宽相邻模态中心频率差很小的时候会直接混叠alpha 太大模态被压缩成极窄的谱线瞬态冲击的能量被分散到多个模态里。从包络熵角度看这两种情况都会让包络波形变复杂熵值升高搜索算法会自动避开。因此边界只要设置在一个合理范围就行不必追求完全贴合信号GWO 的连续搜索能力本来就能找到范围内的最优。4.2 从收敛曲线判断寻优是否正常收敛曲线 convergence 是判断寻优质量的最直接依据。把收敛曲线画出来主要有三种情况figure; plot(1:length(convergence), convergence, b-, LineWidth, 1.5); xlabel(迭代次数); ylabel(平均包络熵); title(GWO 收敛曲线); grid on;第一种是曲线单调下降后趋于平缓这是最理想的状态说明狼群逐步逼近最优解。第二种是前几代快速下降之后长时间不走对应的是算法陷入局部最优此时可以把 SearchAgents_no 从 25 增至 40或者把 lb/ub 范围缩小到已有最优解附近重新搜索相当于在局部区域做精细化寻优。第三种是曲线震荡无法收敛比如适应度值忽高忽低这往往是信号本身具有强随机性、或者 VMD 内部分解不稳定的信号可以用多次分解取平均包络熵的方式来平滑适应度。判断收敛的另一个做法是记录每代的 Alpha_pos如果连续 5 代 K 保持不变且 alpha 变化小于 5%基本可以提前终止迭代。这个早停策略在实际使用中可以省下约三分之一的计算时间且不会损失精度。4.3 用 parfor 并行加速多核 CPU 上的直接收益与边界GWO-VMD 的耗时集中在目标函数反复执行 VMD 分解。灰狼种群中每个个体的适应度计算相互独立完全可以用 parfor 做并行加速。修改方式只有一处把 gwoVmd 中计算适应度的 for 循环换成 parfor并把目标函数写成独立 m 文件。% gwoVmd 函数中 fitnessList zeros(SearchAgents_no, 1); parfor i 1:SearchAgents_no fitnessList(i) vmdObjective(signal, fs, Positions(i, :)); end直接换 parfor 有一个前提vmdObjective 必须能被 workers 访问即它必须是一个独立的 m 文件不能定义为嵌套函数signal 会被复制到每个 worker 上内存占用约等于 worker 数乘以信号长度。当信号超过 10 万点时8 个 worker 同时复制信号可能撑爆内存这种情况建议保留普通 for 循环或者把信号先做降采样再用 parfor 寻优。另注意parfor 内不能使用断点或 disp 输出调试时要切回 for。5. 用 GWO-VMD 跑完后的验证消融对比与包络谱分析GWO-VMD 输出最优参数后必须做的一步是验证分解结果的物理有效性。方法很简单用三组不同参数分别做 VMD再比较包络谱中的故障特征频率幅值。第一组是 GWO 寻优得到的最优参数第二组是经验参数比如 K6、alpha1000第三组是随机参数比如 K4、alpha500。分别用三组参数分解信号取第一个模态做希尔伯特包络谱观察特征频率处幅值。如果 GWO 参数组的特征频率峰值最高、噪声底较低说明寻优有效如果三组结果差距不大说明这个信号本身对参数不敏感直接使用 GWO 参数即可不必每次寻优。% 用最优参数分解 [u_best, ~, omega_best] VMD(x, alpha_opt, tau_opt, K_opt, 0, 1, 1e-7); % 用经验参数分解 [u_exp, ~, omega_exp] VMD(x, 1000, 0, 6, 0, 1, 1e-7); % 计算最优参数下第一模态的包络谱 analytic hilbert(u_best(:, 1)); env abs(analytic); N length(env); f_axis (0:N-1) * fs / N; env_spectrum abs(fft(env)); % 找前三个峰值对应的频率 [pks, locs] findpeaks(env_spectrum, f_axis, ... SortStr, descend, NPeaks, 3); disp(特征频率候选:); disp(locs(1:3));执行这段代码后如果 locs 中的频率与理论特征频率如外圈故障频率 BPFO接近说明分解结果可用如果峰值不明显或频率对不上就要回到目标函数把包络熵换成排列熵试试。还有一个技巧可以确认 K 是否选多检查最优参数对应的中心频率 omega计算相邻 omega 的差值。如果存在两个 omega 非常接近的模态如差值小于频域分辨率的 3 倍说明 K 偏大可以手动降低 K 上界后重新寻优。这个做法在 GWO-VMD 的实际调试中非常管用它能把纯数值寻优的结果与信号处理的物理直觉结合起来。本文还有配套的精品资源点击获取
 函数详解：语法、格式符与性能优化)
