简介：面向计算机、电子信息工程、数学等专业学生及信号处理初学者，这份材料围绕Matlab实现鲸鱼算法WOA优化VMD变分模态分解展开，重点解决模态数k与惩罚因子alpha依赖人工选取、分解效果不稳定的问题。资源包共1个PDF文件，约600KB，集中呈现算法实现说明、参数配置思路与分解结果图示，便于课程设计、期末大作业和毕业设计快速理解流程。文档给出边际谱、频率图、收敛曲线及超参数优化迭代过程图，可观察每次迭代的参数变化与分解效果；同时说明参数化编程、Excel数据替换和主程序运行方式，注释较清晰，适合作为时序分解方向的小型创新点或方法对比参考。已有1337人学习。

1. 从“K 值调到凌晨三点”说起：WOA-VMD 到底替你把哪一步自动化了

上周一个做旋转机械状态监测的同学发来消息，说他手里四十多组 Excel 振动数据，每组都要手动试 K 从 3 到 9、alpha 从 500 到 3000，调一晚上只跑完十几组，第二天起床还得重来一遍。VMD 这套时序分解方法本身不弱，弱在 K 和 alpha 这两个超参数全凭手感，网格搜索的开销随数据量线性上涨，越跑越像玄学。这个项目要解决的正是这最后一步：用鲸鱼优化算法 WOA 在二维空间里搜 (alpha, K) 组合，以包络熵为适应度，迭代收敛后直接给出分解结果，并输出边际谱、各模态时域波形、中心频率散布和收敛曲线。手上有机 Matlab、要做课程设计、毕业设计或者想在论文里补一个参数自适应创新点的开发者，跑通一遍大概十分钟，比手调省下不止一个晚上。

2. 工程骨架：文件清单、主程序入口与 Excel 数据的替换点

拿到这份源码，第一眼看到的是一堆.m文件，而不是一个点开就能跑的可执行程序。这不是缺点，恰恰是这类算法项目的常态——所有逻辑都摊在源码里，随时可以改。问题在于很多人第一次打开会愣一下：哪个是入口？哪个是工具函数？数据又该换个什么格式？

2.1 目录里都有什么，谁是入口

按这类项目的常见组织方式，压缩包里大致是下面这几类文件。不同版本文件名可能略差，但职责基本一一对应：

真正的入口只有一个：main.m。你在命令行里敲main回车，或者直接在编辑器里按 F5，剩下的调用链是main → WOA → fun → VMD → huatu，一层套一层。看代码的时候顺着这条链读，不会迷路。

2.2 main.m 从头到尾走一遍

主脚本是整个项目的配置中心，所有可调的东西都集中在这里。下面这份是我按原项目思路整理出来的结构，字段名和原版接近，注释做了补充：

%% main.m —— WOA-VMD 主程序 clc; clear; close all; rng(1); % 固定随机种子，保证结果可复现 %% 1. 读数据 data = xlsread('data.xlsx'); % 新版 Matlab 可换 readmatrix f = data(:,1)'; % 取第一列并转置成行向量 f = f - mean(f); % 去直流，避免 0 频分量干扰 %% 2. VMD 侧固定参数 tau = 0; % 噪声松弛，取 0 表示不做松弛 DC = 0; % 不强制第一个模态落在 0 频 init = 1; % 中心频率均匀初始化 tol = 1e-7; % 收敛容差 %% 3. WOA 参数 pop = 20; % 种群数量 Max_iteration = 30; % 最大迭代次数 lb = [800, 4]; % [alpha 下界, K 下界] ub = [2000, 10]; % [alpha 上界, K 上界] dim = 2; % 两个待优化参数 %% 4. 目标函数句柄，把固定参数带进去 fobj = @(x) fun(x, f, tau, DC, init, tol); %% 5. 优化 [~, Target_pos, curve, process] = WOA(pop, Max_iteration, lb, ub, dim, fobj); %% 6. 取整后重新分解一次 Target_pos(2) = round(Target_pos(2)); % K 必须是整数 [u, u_hat, omega] = VMD(f, Target_pos(1), tau, Target_pos(2), DC, init, tol); %% 7. 出图 huatu(process, 'VMD');

逐段说清楚关键点。第一步的f = data(:,1)'里那个转置是最容易漏的符号——xlsread读出来是 N×1 的列向量，VMD 内部是按行处理的，不转置的话fft会对整列做运算，结果就是一整条平线，看起来像“分解失败”，其实只是维度错了。去均值那一步同理，直流分量留下来会让中心频率计算在 0 附近打架。

第三步的lb和ub是最值得动手改的地方。原项目给的是[800, 4]到[2000, 10]，K 的范围 4～10 对多数机械振动信号合理，但 alpha 从 800 起步偏高；如果你的信号高频噪声重，alpha 需要更大，如果信号本身平滑，alpha 太大反而会把有用成分压成一条。后面第 5 章会专门讲这个区间调错之后会怎样。

第四步的fobj是一个匿名函数。它把f、tau、DC、init、tol五个固定参数包进了函数句柄，这样一来WOA.m里就只需要关心待优化的二维向量x，不必关心 VMD 那一堆参数怎么传。我第一次看到这种写法的时候还以为是多余的模块化，后来改成换优化器的时候才发现——换掉优化器完全不用动fun.m，多亏了这个句柄。

第五步的返回值process是三维数组，维度是Max_iteration × pop × dim，记录每一代每个个体的位置。它就是后面画“迭代过程图”的数据源，也是这项目比较讨喜的一个细节：别人只有收敛曲线，你有 alpha 和 K 两条各自随代数漂移的轨迹，论文里多一张图。

2.3 换数据：Excel 只需要改两件事

换自己的数据，实际改动量非常小。第一件是把data.xlsx放进当前工作目录，或者把文件名改掉；第二件是确认列的选取。如果你的 Excel 有多列，比如第一列是时间戳、第二列才是幅值，把data(:,1)改成data(:,2)即可。

需要注意两点。一是数据长度，VMD 内部做了镜像延拓，太短的信号（比如只有几百个点）延拓后同相分量太多，边界效应会很严重，一般建议至少两千点。二是数据的平稳性，WOA-VMD 适合分析非平稳、多分量的信号，比如轴承振动、肌电、心电、地震波，如果你拿一段纯正弦进去，K 优化到最后很可能收敛到 1 或 2，这不是算法出问题，是信号本身就没有多模态结构。

% 如果 Excel 是时间 + 幅值两列，改成这样 data = readmatrix('data.xlsx'); t = data(:,1); f = data(:,2)'; f = f - mean(f);

readmatrix是 Matlab R2019a 之后推荐替代xlsread的接口，速度更快，而且不需要装 Excel 或者装 COM 组件。老版本 Matlab 用xlsread也能跑，只是在没有 Excel 的服务器上会报错，这一点写论文跑服务器批处理的时候容易翻车。

3. 核心模块拆解：WOA 如何把 (alpha, K) 从二维空间里搜出来

优化器部分是这个项目最有含金量的地方，也是标题里“鲸鱼算法”四个字的实体。整段代码不到一百行，但要读懂它得先明白一件事：VMD 的参数优化本质上是一个没有解析梯度的黑箱问题——你没法对 alpha 求导，只能“给一组参数、跑一次分解、看看熵值多少”，所以只能用启发式群体算法。

3.1 搜索空间的定义：lb、ub 与 dim 三件套

优化器只认识三个东西：下界、上界、维度。对应到本项目就是：

dim = 2是因为只优化两个参数。原始代码里写死的lb = [800,4]、ub = [2000,10]是一组比较保守的默认值，能跑通，但换个信号质量差的数据集大概率要重新调。我的习惯是先把区间放宽到[300, 4000]和[3, 12]，跑一次看看最优解落在哪个角落，收窄范围后再跑一次正式的。这一步多花二十分钟，换来的收敛质量提升相当明显。

3.2 三段式捕食行为写成代码是什么样

WOA 的核心思想是把寻优过程拆成三种行为：包围猎物、螺旋气泡网攻击、随机搜索。代码里的分支结构就是这三种行为的直译：

function [Leader_score, Leader_pos, curve, process] = WOA(N, Max_iter, lb, ub, dim, fobj) % N 种群数量 % Max_iter 最大迭代次数 % lb/ub 搜索空间的下界/上界（1×dim 行向量） % fobj 适应度函数句柄，越小越好 % ---- 1. 随机初始化种群，落在 [lb, ub] 内 ---- Positions = rand(N, dim) .* (ub - lb) + lb; Leader_pos = zeros(1, dim); Leader_score = inf; curve = zeros(1, Max_iter); process = zeros(Max_iter, N, dim); % 记录每一代全部个体位置 for l = 1:Max_iter a = 2 - l * (2 / Max_iter); % 从 2 线性衰减到 0，控制探索→开发的切换 for i = 1:N % 越界处理：截断到边界内 Positions(i,:) = max(Positions(i,:), lb); Positions(i,:) = min(Positions(i,:), ub); % 评估当前个体 Fit = fobj(Positions(i,:)); if Fit < Leader_score Leader_score = Fit; Leader_pos = Positions(i,:); end r1 = rand(); r2 = rand(); A = 2 * a * r1 - a; % A 的绝对值小于 1 时开发，大于 1 时探索 C = 2 * r2; p = rand(); b = 1; % 螺旋形状常数，一般取 1 if p < 0.5 if abs(A) < 1 % ---- 2a. 包围猎物：朝当前最优个体靠拢 ---- D = abs(C * Leader_pos - Positions(i,:)); Positions(i,:) = Leader_pos - A .* D; else % ---- 2b. 随机搜索：随机挑一个同伴当参照 ---- idx = randi(N); Xr = Positions(idx,:); D = abs(C * Xr - Positions(i,:)); Positions(i,:) = Xr - A .* D; end else % ---- 3. 螺旋气泡网：绕最优个体做对数螺旋移动 ---- D = abs(Leader_pos - Positions(i,:)); Positions(i,:) = D .* exp(b) .* cos(2*pi*rand()) + Leader_pos; end end process(l,:,:) = Positions; % 存档 curve(l) = Leader_score; end end

几处参数值得单独拎出来说。a从 2 衰减到 0 是整段代码的灵魂，它决定了 A 的绝对值大小——A 在 (−1,1) 区间意味着“围绕当前最优解精细搜索”，在区间外意味着“去远处探索”。迭代前期a大，A 大概率越界，种群到处跑；迭代后期a小，A 收敛到小区间，种群开始局部开发。这个衰减规律是 WOA 区别于其他群智能算法的地方，也是它被论文引用最多的原因。

A = 2*a*r1 - a那一行里，r1是每次调用重新生成的随机数，所以同一个个体在不同迭代里的行为是随机的，这是群体算法的常态。如果你在调试时发现结果不稳定，第一反应应该是固定rng，而不是去改更新公式。

process(l,:,:) = Positions这一行是为了后面画迭代过程图。它会把每一代所有个体的位置都存下来，内存占用跟Max_iteration × pop × dim成正比，20×30×2 的数据量完全不用担心，但如果你把 pop 提到 200、迭代提到 500，就得考虑只存最优个体了，不然内存会被吃光。

3.3 适应度函数 fun.m：包络熵为什么能当尺子

优化器只负责“往哪走”，判断“好不好”的工作交给了fun.m。原始项目里用的准则是包络熵最小化，逻辑是——VMD 分解得好不好，就看每个模态的包络有多“稀疏”，包络越集中（熵越小），说明这个模态对应一个清晰的振荡成分，而不是噪声的叠加。

function fit = fun(x, f, tau, DC, init, tol) % x(1) = alpha，x(2) = K（连续值，需要取整） % 返回值越小，这组参数越优 alpha = x(1); K = round(x(2)); % 关键：在适应度函数内部就取整，不要拖到最后 % 取整之后可能越界，做一次硬保护 if K < 2, K = 2; end if K > 15, K = 15; end % 用当前参数跑一次 VMD，出错就返回一个大惩罚值 try [u, ~, ~] = VMD(f, alpha, tau, K, DC, init, tol); catch fit = 1e6; return; end % 逐模态计算包络熵，求和 fit = 0; for i = 1:K env = abs(hilbert(u(i,:))); % 解析信号 → 取模得到包络 p = env / (sum(env) + eps); % 归一化成概率分布 fit = fit - sum(p .* log(p + eps)); % 香农熵，负号保证熵越大值越小 end end

三段关键逻辑。第一段round(x(2))是必须在fun.m内部做的——如果你留到最后才取整，优化器看到的 K 是 4.1、4.5、4.9 三个不同值，但实际上它们跑的都是同一个 K=4 的分解，适应度景观就变成了阶梯状的平台，WOA 在平台上来回震荡却找不到下降方向，收敛曲线会拉出锯齿。把这个动作提前到适应度函数内部，优化器认知里的“好”和“坏”就与真实分解结果一致了。

第二段的try-catch是防止参数组合把 VMD 跑崩。某些 alpha 值配某些 K 值，ADMM 迭代可能不收敛，或者返回 NaN，包一层异常处理可以避免整个优化过程中断，代价是返回一个不可能被选中的大值。

第三段的包络熵公式：abs(hilbert(u))求解析信号模得到包络，然后归一化成概率分布，再算香农熵。因为原始熵的定义里前面有个负号，取负之后熵越大值越小，正好符合“最小化适应度”的方向。如果你直接抄其他项目的代码，经常会看到SSA_cg_curve = abs(SSA_cg_curve)这一行，它是把某些优化器返回的负值全部取正，让曲线朝下的方向统一——这行不影响结果，但会让收敛曲线上多出一段“已经收敛却在往上走”的假象，读图时要心里有数。

3.4 整数参数与连续参数的边界处理

K 是整数，alpha 是连续的，这两个参数放在同一个优化器里跑，边界处理比想象中麻烦。上界ub(2)=10意味着 K 只能到 10，如果取整之后是 10，包络熵偏大，优化器会慢慢把 K 往 9、8 推。这时候如果ub(2)再设成 8，9 就被掐掉了，最优的 K 可能正好在 8 和 10 之间，你永远找不到。

我的做法是把 K 的搜索区间设成比预期最优点宽两格，比如预判最优 K=6，就把lb(2)=3, ub(2)=10放开，让优化器自己去收敛。收敛之后如果发现 Target_pos 贴到了上界，说明上界还是设窄了，再往上放一次。这个操作两三次就到位。

4. VMD 本体与出图：从 u、u_hat、omega 到边际谱

优化器找出来的只是一组参数，真正完成分解的是VMD.m。这个函数在学术界有标准实现，代码结构相对固定，但项目对它做了一些适配，比如输入输出注释的中文化、镜像延拓的处理、收敛条件的简化。

4.1 VMD 函数签名逐参数说明

函数的输入部分在原项目里保留了英文注释，正好可以当参数说明用：

function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol) % Input and Parameters: % --------------------- % signal - 待分解的时域信号（行向量） % alpha - 惩罚因子，控制各模态的带宽 % tau - time-step of the dual ascent（0 表示不松弛） % K - 模态数，也就是分解层数 % DC - true 时第一个模态固定在 0 频 % init - 0 = 中心频率全部从 0 开始 % 1 = 均匀分布初始化 % 2 = 随机初始化 % tol - 收敛容差，一般 1e-6 ~ 1e-7 % % Output: % ------- % u - K×N 矩阵，每一行是一个模态的时域信号 % u_hat - K×N 矩阵，每一行是一个模态的频谱 % omega - K×2 矩阵，第一列是当前中心频率，第二列是上一轮的值

输入参数里最容易被忽略的是init。理论上 0/1/2 三种初始化都能收敛到相近结果，但在实际数据上，init=0（全从 0 频起步）有时会让两个模态的中心频率“抢”同一个位置，最后合成了一个；init=1（均匀分布）最稳，几乎是这类项目的默认选项；init=2（随机）只在做多次重复实验取平均的时候用得上。如果你懒得想，就直接用init=1。

输出方面，u是最终结果，omega是诊断工具。omega的数值是归一化频率，取值在 0 到 0.5 之间，乘以采样频率fs才是真实的中心频率（单位 Hz）。这一点很多人在写论文的时候会搞错，把 0.23 直接当成 0.23 Hz 写进表格里，审稿人一看就发现问题。

4.2 分解主循环：ADMM 在做什么

VMD 的核心是 ADMM 交替方向乘子法，把原问题拆成三个可分别求解的子问题。代码里对应的是三层嵌套循环，逻辑可以简化为下面这样：

% 1. 镜像延拓，抑制端点效应（两端各延一倍，总长变 3L） L = length(signal); f_mirror = [fliplr(signal), signal, fliplr(signal)]; T = length(f_mirror); freqs = (0:T-1)/T - 0.5; % 归一化频率轴，范围 [-0.5, 0.5) f_hat = fftshift(fft(f_mirror)); % 频域 % 2. 初始化中心频率和拉格朗日乘子 omega = zeros(K, 2); switch init case 1, omega(:,1) = (0.5/K) * (0:K-1)'; case 2, omega(:,1) = sort(rand(K,1)); otherwise, omega(:,1) = 0; end omega(:,2) = omega(:,1); u_hat = zeros(K, T); lambda_hat = zeros(1, T); % 3. ADMM 主循环 n = 0; uDiff = tol + eps; while uDiff > tol && n < 500 sum_omega = 0; for k = 1:K % 3a. 更新第 k 个模态的频谱 sum_uk = sum(u_hat,1) - u_hat(k,:); u_hat(k,:) = (f_hat - sum_uk - lambda_hat/2) ./ ... (1 + alpha * (freqs - omega(k,1)).^2); % 3b. 更新第 k 个模态的中心频率（频谱重心） num = sum(freqs .* abs(u_hat(k,:)).^2); den = sum(abs(u_hat(k,:)).^2) + eps; omega(k,2) = omega(k,1); omega(k,1) = num / den; omega(k,1) = min(max(omega(k,1), -0.5), 0.5); sum_omega = sum_omega + abs(omega(k,1) - omega(k,2))^2; end % 3c. 更新拉格朗日乘子 lambda_hat = lambda_hat + tau * (sum(u_hat,1) - f_hat); uDiff = sum_omega; n = n + 1; end % 4. 反变换并裁掉镜像延拓部分 u = real(ifft(ifftshift(u_hat, 2), [], 2)); u = u(:, L+1 : 2*L);

第 3a 步的公式是 VMD 最核心的更新式：分子是“剩余信号”（原始信号减去其他模态和拉格朗日项的贡献），分母是一个以omega(k,1)为中心的二次滤波响应，alpha越大，滤波器越窄，各模态分得越开。这也解释了为什么 alpha 需要优化——alpha 太小，模态之间混叠；alpha 太大，一个完整的模态可能被劈成两半。

第 3b 步用的是频谱重心法，本质上是该模态功率谱的加权平均频率，每次迭代重新算一遍。omega(k,2)存的是上一轮的值，用来判断本轮迭代是否已经稳定。

第 3c 步更新拉格朗日乘子，tau取 0 时这一项不生效，等价于不做约束，这也就是输入注释里说的“noise-slack”。实际项目中tau=0是主流做法，因为大多数振动信号本身带噪声，硬约束反而会拖慢收敛。

4.3 出图部分：边际谱、频率图、收敛曲线怎么组织

图是这个项目的卖点之一，原项目明确说了“包含边际谱、频率图、收敛曲线”。这些东西之所以值得单独做函数，是因为放在main.m里每次都要手动调整，最后画出来的图彼此不一致，论文插图的观感会被扣分。下面这段是按原思路重写的huatu和出图部分：

%% 出图：分解结果 + 频率图 + 边际谱 figure('Name','VMD分解结果'); subplot(K+1,1,1), plot(f,'k'), title('原始信号'), axis tight; for i = 1:K subplot(K+1,1,i+1), plot(u(i,:),'b'); title(['IMF ',num2str(i)]), axis tight; end figure('Name','各模态中心频率'); stem(omega(:,1) * fs, 'filled'); % 归一化频率 × 采样率 = 真实 Hz xlabel('模态序号'), ylabel('中心频率 / Hz'), grid on; figure('Name','边际谱'); hold on; for i = 1:K U = abs(fft(u(i,:))); fx = (0:length(U)-1) * fs / length(U); idx = 1:floor(length(U)/2); plot(fx(idx), U(idx), 'DisplayName', ['IMF',num2str(i)]); end xlabel('频率 / Hz'), ylabel('幅值'), legend show, grid on;

三张图各有各的诊断作用。模态波形图看的是什么？看每个模态是不是一个“规整的振荡”——如果某个模态看起来还是一团噪声，说明 K 给多了或者 alpha 给错了。中心频率图看的是什么？看 K 个频率点在频轴上是不是分得开——如果两个模态的频率在地图上几乎重合，K 就过多了。边际谱则是把所有模态的频谱叠在一张图上，看整体能量分布在哪些频段，跟原始信号的频谱对比，可以直接看出分解有没有漏掉成分。

我自己的习惯是把这三张图和原始信号的频谱放在一起看。如果边际谱的最高峰和原始信号的最高峰位置差超过 5%，那基本可以确定 alpha 或者 K 没调好，不要继续往下做后续分析了，先回头把参数重新跑一遍。

4.4 收敛曲线和迭代过程图的读法

huatu画的不是收敛曲线，而是两条“参数轨迹”：alpha 和 K 各自随迭代代数的变化，横轴是迭代次数，纵轴是参数值，两条线分别是种群平均值和种群最小值。下面这个简短版本捕获了它的核心信息：

function huatu(process, tag) % process: Max_iter × pop × dim，每代每个个体的位置 [Max_iter, ~, ~] = size(process); figure('Name',[tag ' 参数迭代过程']); subplot(1,2,1); plot(mean(process(:,:,1),2),'r-'); hold on; plot(min(process(:,:,1),[],2),'b-'); xlabel('迭代次数'), ylabel('alpha'); legend('种群均值','种群最优'), title([tag ' — alpha']); subplot(1,2,2); plot(mean(round(process(:,:,2)),2),'r-'); hold on; plot(min(round(process(:,:,2)),[],2),'b-'); xlabel('迭代次数'), ylabel('K'); legend('种群均值','种群最优'), title([tag ' — K']); end

读图重点在于“种群均值的走势”。如果均值在第 5 代左右就基本贴合最优点，说明种群多样性不足——要么初始值凑巧就落在了最优点附近，要么 pop 太小、随机性不够，后期基本靠最优个体带着走。如果到第 20 代均值还在明显波动，说明迭代次数不够，再加十几代看看。还有一种情况是 K 那条线一直贴着上界走，这几乎等于直接告诉你ub(2)设小了，K 的真实最优值在区间之外，回去改参数。

真正的收敛曲线（适应度值随迭代的变化）由优化器返回的curve提供，通常是单调不增的，中间可能出现小台阶。台阶的高度对应一次“找到了更好的解”，两段台阶之间的平坦期越长，说明搜索越困难。如果你看到曲线掉下去之后又翘上来，那一定是代码里对适应度值做了绝对值处理或者符号翻转，回去检查一下。

5. 避坑与排查：五条血泪经验

这段是我自己在这个项目上踩出来的坑，按“现象 → 原因 → 解决”的顺序整理。里面的每条几乎都有同学在群里问过。

5.1 参数与调用链上的坑

坑 1：函数名和注释对不上，也不知道到底在优化什么

现象：main.m里调用的是SSA(pop, ...)，注释里也写着“麻雀算法参数设置”，但项目标题、论文标题都明明白白写着“鲸鱼算法 WOA”。

原因：这类模板代码在不同群智能算法之间共用得非常多，作者复制粘贴之后改了标题，函数名没改，SSA.m里说不定装的其实是 WOA 的更新公式。

解决：打开优化器文件，看位置更新那段。如果出现A = 2*a*r1 - a、C = 2*r2以及exp(b).*cos(2*pi*rand())这个螺旋项，那它就是 WOA 无疑，只是文件名叫SSA。如果里面出现发现者-加入者-警戒者的三段结构，那才是真正的麻雀算法。论文或者答辩里必须把这一点理顺，不然被问到怎么答都别扭。

坑 2：目标函数返回的曲线是反向的，最优解反而不如图中最低点

现象：curve画出来是上升的，或者虽然看着在下降，但最终选出的Target_pos对应的是曲线上的一个高点。

原因：某些优化器接口把适应度定义为“越大越好”，而本项目是按“越小越好”写的。原项目里那句SSA_cg_curve = abs(SSA_cg_curve)就是为了把某些返回负值的结果统一翻成正的，但在别的优化器上直接照抄，符号就不对了。

解决：先把所有优化器统一改成“最小化适应度”，fun.m返回什么，优化器就按什么比。然后在WOA.m里用if Fit < Leader_score做判断，这一行写死，不要再去做取绝对值之类的处理。如果曲线露出来的形状还是不对，先单步打印前三次迭代的Fit值，对照包络熵数量级（通常在 3~15 之间）看一眼就能判断方向对不对。

坑 3：K 的取整放在最后一刻，收敛曲线全是台阶

现象：优化过程收敛得很慢，曲线一会儿降一会儿升，Target_pos出来之后取整，发现跟之前某一次试的 K 值一模一样。

原因：整数参数在连续空间里优化，round放在最后，优化器不知道 K 的“真实”值是多少，把 4.3 和 4.7 当成两个不同的解来比较，实际跑出来的却是同一个分解结果。适应度景观是阶梯状的，群智能算法在这种景观上表现极差。

解决：把取整放进fun.m，位置就在函数第一行。优化器看到的每个个体都是“取整之后实际跑的分解”，适应度的可比性就有保证。

5.2 数据与运行环境上的坑

坑 4：Excel 读进来方向错了，分解结果是平的

现象：VMD调用不报错，返回的u全是接近 0 的数值，omega里所有中心频率都堆在 0 附近。

原因：xlsread('data.xlsx')读到的是 N×1 列向量，如果不转置，VMD里面length(signal)得到的是 1，fft对整列做运算，结果全乱。或者是数据带直流没去，第一个模态把 0 频能量全部吸过去。

解决：f = data(:,1)';转置，然后再f = f - mean(f)去直流。两者都加，不要只做一个。判断方法很简单：size(f)打出来应该是1×N，如果看到N×1就是没有转置。

坑 5：种群只有 5 个，迭代只有 20 次，跑三次出三个不同结果

现象：每次清空工作区重新跑，得到的最优 K 分别是 5、7、6，alpha 也上下跳了快一百。

原因：原项目模板里pop = 5、Max_iteration = 20是演示配置，不是实际配置。5 个个体的种群在二维空间里几乎等于随机抽样，20 代还不到 WOA 该有的收敛速度，结果完全靠运气。

解决：pop提到 20～30，Max_iteration提到 30～50，并且加上rng(1)固定随机种子。改完之后再跑三次，如果结果还在大幅波动，那就不是种群小的问题，而是搜索区间设得不对，回去看第 3 章里那段关于 lb / ub 的说明。

6. 二次开发的几个下手点：换适应度、换优化器、批量跑

项目跑通之后想往论文方向走一步，改动量最大的往往不是 VMD 本体，而是三处边界：适应度函数、优化器、数据入口。每个地方的接口都很清晰，适合动手。

6.1 换适应度：从单一包络熵到复合指标

fun.m的返回协议就是“给一个数，越小越好”。想换成排列熵、样本熵、峭度，或者换成“平均包络熵 + K 的轻罚项”都可以，只要保证返回一个数值：

% 在 fun.m 末尾追加软惩罚，防止优化器把 K 一路推向边界 fit = fit + 0.002 * K;

加这一项的意义在于——原版的纯包络熵准则偏爱大 K，因为模态越多包络越集中，熵的绝对值往往越低。加上一个随 K 线性增长的惩罚，相当于在“分解得细”和“分解得少”之间找平衡，K 不会一路涨到上界。

6.2 换优化器：签名对齐就能无缝替换

把WOA.m换成一个灰狼优化器或者粒子群优化器，只要输入输出签名保持一致——输入(N, Max_iter, lb, ub, dim, fobj)，返回(Leader_score, Leader_pos, curve, process)，main.m一行都不用改。这是这个项目在结构上做得比较聪明的地方：把优化器当作可插拔模块，算法本身只是一个黑箱。

6.3 批量跑多个 Excel

main.m适合跑单个文件，要跑几十组，包一层循环即可。大致的骨架是：

files = dir('data/*.xlsx'); Results = zeros(numel(files), 3); % 每行存 [最优 K, 最优 alpha, 最终熵] for k = 1:numel(files) data = readmatrix(fullfile('data', files(k).name)); f = data(:,2)'; f = f - mean(f); fobj = @(x) fun(x, f, 0, 0, 1, 1e-7); [~, pos] = WOA(20, 30, [300,3], [4000,12], 2, fobj); Results(k,:) = [round(pos(2)), pos(1), fobj(pos)]; end writematrix(Results, 'summary.xlsx');

这段跑起来会比单文件慢一个数量级，因为每个文件要跑几十次 VMD。如果机子是多核的，把最外面这层for换成parfor能快三到四倍，前提是你已经装了并行计算工具箱，并且每个文件数据量都差不多、不会导致负载严重不均。

最后说个习惯。我早期做这类参数优化，跑出个好结果就直接进论文，从来没回头验一次。后来有次图省事，把alpha的搜索区间从 2000 直接砍到 1200，结果去跟另一组手调的结果对比时才发现，包络熵看起来更小，但模态混叠得厉害——熵小只是分布更集中，不代表分得对。从那以后我每次跑完优化，都会单独用最优参数手跑一次 VMD，对着原信号频谱目测一遍，确认模态频率没有并在一起，才敢往下做。这个动作多花不到五分钟，省下的却是整段重做的后悔药。希望帮到你。

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