灰狼算法优化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 为何需要 GWO：K 和 alpha 靠手调不靠谱

2.1 VMD 的核心参数与物理含义

VMD 把信号分解问题写成变分模型：寻找 K 个模态函数，使所有模态的带宽之和最小，同时这些模态之和能精确重构原信号。迭代求解时，alpha 是带宽惩罚因子：alpha 越大，各模态带宽越窄、频带越干净，但过大可能丢失瞬态冲击的能量；alpha 太小，相邻模态之间就会发生频谱重叠。除 K 和 alpha 外，tau 控制拉格朗日乘子的更新步长；DC=1 表示第一个模态包含直流分量；init=1 表示初始中心频率为均匀分布。

实际信号分析中，K 是最让人头疼的参数。K 设小了，多个频率成分叠在同一个模态里；K 设大了，又会出现相邻中心频率几乎重合的空模态或碎片模态。理论上可以通过观察中心频率是否聚集来判断 K 是否合适，但当信号包含噪声与多个边带时，肉眼判断很不靠谱。所以 GWO-VMD 的思路很直接：把 K、alpha、tau 一起放进一个优化目标里，用智能搜索代替手工试探。

2.2 为什么选灰狼优化算法而不是网格搜索或 PSO

网格搜索的第一问题是维度爆炸，K 取 2～10、alpha 取 200~2000 步长 100，组合数就是 9*19=171 次完整 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 最小 200，tau 最小 0 ub = [10, 2000, 0.3]; % K 最大 10 层，alpha 最大 2000，tau 最大 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 的搜索区间很小，如果为了减少搜索维度，可以固定 tau=0，把 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 行，此处省略 --- end

VMD 的五个参数里，注意 DC 参数对第一模态的影响。当 DC=1 时，第一个模态允许包含直流分量，适合处理含有明显偏置的传感器信号；DC=0 则强制所有模态围绕各自的中心频率带限分布。init=1 是均匀初始化，保证每次分解结果可复现；改成 init=0 从零开始迭代，会引入随机性，GWO 寻优时最好别用 init=0，否则同一组参数两次分解的包络熵会不一样，直接影响适应度对比。

4. 参数设置、收敛判据与运行提速的实用技巧

4.1 搜索边界和 GWO 参数怎么定：一张表说清推荐值

GWO-VMD 运行前要确定的参数包括 GWO 自身的参数和 VMD 参数的搜索范围。下表汇总了我在电力谐波、轴承故障、齿轮箱振动三种典型信号上使用的推荐值：

关于 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 寻优得到的最优参数；第二组是经验参数，比如 K=6、alpha=1000；第三组是随机参数，比如 K=4、alpha=500。分别用三组参数分解信号，取第一个模态做希尔伯特包络谱，观察特征频率处幅值。如果 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 的实际调试中非常管用，它能把纯数值寻优的结果与信号处理的物理直觉结合起来。

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