1. 为什么选用鲸鱼算法，又为什么非要“改进”
1.1 原版鲸鱼优化算法的三个核心算子
鲸鱼优化算法是模拟座头鲸气泡网捕食行为的一种群智能算法，我第一次接触它是在做一个航迹规划的对比实验。这套算法的核心只有三个算子：包围猎物、气泡网攻击、随机搜索猎物。看起来很简单，和粒子群算法、遗传算法相比，它最吸引人的地方是“螺旋更新”这个机制，可以让种群在局部搜索和全局搜索之间来回过渡。
包围猎物的公式是：
X(t+1) = X*(t) - A * D D = |C * X*(t) - X(t)|其中 A 和 C 是系数向量，A 在迭代过程中从 2 线性衰减到 0。当 |A| < 1 时，个体向当前最优个体收缩；当 |A| >= 1 时，个体远离当前最优，执行随机探索。这就是鲸鱼算法的搜索逻辑：用阈值的切换，模拟围捕和游走之间的平衡。
气泡网攻击阶段，算法用 50% 的概率在收缩包围和螺旋路径之间选择。螺旋更新的核心是用对数螺旋公式产生新位置，让个体围绕最优解螺旋靠近。这个机制和粒子群算法的“惯性+社会项”思路完全不同，它更像是在最优解周围做精细扫描。我在复现时最大的感受是：这套机制在低维连续优化问题上表现相当不错，尤其对某些基准函数，它的收敛速度明显快于遗传算法。
1.2 原版算法在三维航迹规划下的两大致命短板
但把它直接搬到无人机三维航迹规划里，问题立刻暴露出来。三维航迹规划的搜索空间不是普通的函数优化空间，而是带有地形约束、禁飞区约束和动力学约束的复杂离散空间，航迹通常由 20～50 个航迹点串联而成。如果一个航迹点包含 X、Y、Z 三个坐标，那么问题的搜索维度就是 60～150 维。原版鲸鱼算法在 30 维以下的测试函数上表现稳定，但一旦升到上百维，它的收敛速度断崖式下降。
第一个明显短板是“早期收敛过快导致的多样性崩塌”。因为算法的开发能力很强，一旦某个个体找到了一个相对较优的航迹，其他个体会迅速向它靠拢，导致种群几乎全部堆积在一个局部区域。三维地形中存在大量山谷、山峰和局部平坦区域，这样的情况最容易把种群锁死在局部最优解。我在测试中把原版鲸鱼算法放到一个带山峰障碍的 100km×100km 地图上跑，前 50 次迭代看起来航迹长度下降得很快，之后几乎不再变化，最后得到的航迹常常是一条绕远路但找不到真正最短路径的结果。
第二个短板是“收敛因子线性衰减不符合实际搜索节奏”。原版算法的 A 值在一个迭代周期内均匀地从 2 变到 0，这意味着“探索”和“开发”的分配是全周期均匀的。但在三维航迹规划中，前期需要大量的全局搜索来找“大致可飞区域”，后期才需要在找到的可行通道里精细调整航迹点。线性衰减使前期探索力度不足，后期又缺少足够的开发强度，这是许多改进论文选择非线性收敛因子的直接原因。
所以，题目里“改进”这两个字是必要的，不是在论文里凑花样。我复现时采用的改进策略也是围绕这两个短板展开的：用一个非线性变化策略替代线性衰减，增加扰动机制来保持种群多样性，同时引入精英保留来保证最优解不丢失。
2. 三维航迹规划的代价函数与约束建模
2.1 地形、禁飞区与航迹编码方式
做三维航迹规划之前，必须先把问题“可计算化”。这个阶段不直接涉及算法，但算法效果好坏在很大程度上取决于这一步建模精度。
地图方面，常见做法是用数字高程模型模拟地形，也可以用数学函数生成一个带多峰的地形曲面。我复现时用了两种地形做交叉验证：一种是用高斯型山丘叠加生成的地形函数，代码如下：
def generate_terrain(n=100, peaks=12): x = np.linspace(0, 100, n) y = np.linspace(0, 100, n) X, Y = np.meshgrid(x, y) Z = np.zeros_like(X) for _ in range(peaks): cx = np.random.uniform(10, 90) cy = np.random.uniform(10, 90) h = np.random.uniform(20, 50) sx = np.random.uniform(8, 20) sy = np.random.uniform(8, 20) Z += h * np.exp(-((X - cx) ** 2) / (2 * sx ** 2) - ((Y - cy) ** 2) / (2 * sy ** 2)) return X, Y, Z这样生成的地形可以精确控制山峰位置和高度，方便评估算法是否真正找到了绕过山峰的路径。实际工程中可以使用公开的高程数据文件，但因为数据格式和插值方式的问题，我建议先用函数生成地形，跑通算法后再替换真实高程数据。
航迹编码方式我用的是“固定点数-可变高度”的方案：在起点和终点之间均匀插入 N-2 个航迹点，所有点的 X、Y 坐标等间隔分布，算法优化的变量是每个中间点的高度 Z。这样编码的好处是维度直接等于航迹点数量，搜索空间从三维空间中的连续曲面变成了“高度向量”空间，大幅降低了优化难度。如果不做这种降维处理，直接把所有坐标都作为变量，算法的收敛速度和成功率都会明显变差。
2.2 总代价函数的拆分：长度、威胁、高度、平滑度
航迹规划的“最优”是多目标折中的结果。我用的代价函数由四项组成：
J = ω1 * J_length + ω2 * J_height + ω3 * J_smooth + ω4 * J_threat- J_length 是航迹总长度，单位 km，直接反映飞行距离。
- J_height 是所有航迹点的平均飞行高度与预设巡航高度之间的偏差，用来控制“贴着地形飞”和“飞太高”之间的平衡。
- J_smooth 是相邻航迹段之间的转角和爬升率惩罚，保证航迹的连续性满足无人机运动学约束。
- J_threat 是所有航迹点与禁飞区圆心之间距离的惩罚项，航迹穿过禁飞区时该项急剧增大。
其中威胁代价实现时不是简单二值化判断，而是用距离的倒数来构建连续惩罚。这样算法在处理时能感受到“靠近禁飞区边缘”的连续梯度，而不是只有“通过/不通过”的突变值。这个细节很重要，因为群智能算法对跳变的适应度函数非常敏感，突变会造成搜索方向的割裂。
高度代价也要注意单位问题。地图范围是 100km×100km，高程是百米级，但航迹长度是公里级。如果不加权系数，长度代价的数值量级会彻底压过高度代价，导致算法只在水平面上优化高度。我的做法是先对每一项做归一化，再乘权重系数，保证四项在数值上处于同一量级。
2.3 约束处理：越界、禁飞区、转角限制
约束处理我采用的是“惩罚函数+边界修复”的混合方案。越界的航迹点直接拉回边界，例如：
z_new = np.clip(z_new, z_min, z_max)禁飞区和转角约束则放在代价函数里处理，权重给得比较大。这样的设计既保证了搜索的连续性，又不会让不可行解遍布整个种群。需要明确的是，惩罚权重并不是越大越好，我测试过把威胁权值设成 10^4，结果算法几乎不敢偏离初始航迹，收敛速度非常慢。合理的做法是先跑一次不加威胁项的实验，观察最优航迹的威胁代价量级，再按量级设置惩罚权重。
3. 我采用的改进策略及其实现细节
3.1 非线性收敛因子与自适应惯性权重
第一个改进是针对 A 值线性衰减的问题。我用的非线性收敛因子是余弦型递减曲线：
a = 2 * cos(0.5 * pi * (t / T)) ^ 2- t 是当前迭代次数
- T 是最大迭代次数
这条曲线的特点是：前期 a 下降缓慢，算法有更多概率执行全局探索；后期 a 快速降到接近 0，增强局部开发。从实际效果看，在 100 维航迹规划任务上，最优航迹长度比原版减少了大约 8%～12%。
同时我加入了基于迭代次数的自适应惯性权重：
w = w_max - (w_max - w_min) * (t / T) ^ 2这个权重直接乘在位置更新公式中的原位置项上，用于控制个体对上一时刻位置的记忆程度。前期 w 大，保持个体的飞行惯性，避免过早聚拢；后期 w 小，让个体更信任当前搜索方向。这两项改进在原理上是互相配合的：非线性收敛因子管的是 A 的幅度变化，惯性权重管的是位置更新的平滑程度。
3.2 柯西-莱维混合变异扰动
第二个改进是变异机制。原版鲸鱼算法没有变异算子，种群一旦收敛就基本定形。我参考了群智能算法中常用的变异策略，采用柯西分布与莱维飞行混合的扰动方式：
X_new = X_best + step * Levy(beta) + Cauchy_noise莱维飞行的步长分布具有重尾特性，可以让个体偶尔跳出局部区域；柯西分布的扰动幅度相对集中，可以更好地在当前最优解附近做精细调节。两者叠加，一定程度上兼顾了全局跳出和局部精细搜索。
每次迭代中，我设定一个变异概率 p_mut，对当代最优个体执行变异，并将变异后的个体与当前最优个体做比较，保留更优者。这样做的好处不会破坏主搜索流程，同时又给算法提供了额外的新解输入。
3.3 精英保留与动态自适应种群
第三个改进是精英保留机制。在主循环中，每一代更新完种群位置后，我将当前代的最优个体单独保存，并在下一代更新完成后，用上一代最优个体替换掉当前代最差个体。这种“代际替库”能避免最优解在变异和位置更新过程中被意外覆盖，也保证了算法的收敛目标始终是历史最优。
除此之外，我还做了一个动态种群裁剪：迭代前 60% 使用较大的种群规模（例如 50），后 40% 将种群规模逐步缩小到 25。这个策略的代价是切换时需要重新评估被淘汰个体的适用性，对于已收敛的种群体量影响不大。如果追求更快的计算速度，也可以直接从第一轮就保持固定规模，差异并不显著。
3.4 改进算法流程串一遍
完整流程总结如下：
- 初始化地形、起点终点、禁飞区参数。
- 生成初始种群，每个个体表示一组航迹点高度。
- 计算每个个体的代价函数值。
- 进入主循环，对每个个体执行包围或气泡网更新，再乘以自适应权重。
- 以变异概率对当前最优个体执行柯西-莱维混合扰动。
- 精英保留：替换当代最差个体。
- 若未达到最大迭代次数，重复 4～6；否则输出最优航迹。
- 对最优航迹执行平滑插值，绘制三维航迹图。
4. Python代码实现的核心结构与关键代码
4.1 代码组织与文件结构
这个复现项目我组织了五个模块，各自职责单一。在主文件main.py里，先加载地形生成函数、初始化代价函数模块，然后实例化改进鲸鱼算法类并运行 500 次迭代。
├── data_terrain.py # 地形生成与高程数组管理 ├── cost_function.py # 代价函数：长度、高度、平滑、威胁 ├── iwoa_algorithm.py # 改进鲸鱼算法类 ├── utils_plot.py # 三维轨迹绘图与数据导出 └── main.py # 主流程入口在复现论文代码的思路里，独立模块化是必须的。如果把所有函数打成一个大脚本，后面调整参数和做对比实验时会非常痛苦。你可能想快速跑通结果再把代码整理干净，但事实上做对比实验时，你几乎每天都要反复修改代价函数和参数，模块化能帮你省下大量时间。
4.2 航迹编码、解码与初始种群生成
我采用的编码是“N个中间航迹点的高度Z_i”。初始化时用随机均匀分布生成每个高度值。下面是种群初始化的核心逻辑：
class IWOA: def __init__(self, n_agents, dim, bounds, max_iter): self.n_agents = n_agents self.dim = dim self.bounds = bounds self.max_iter = max_iter self.pop = np.random.uniform(bounds[0], bounds[1], (n_agents, dim)) def decode(self, z_series): # z_series是N-2个中间点高度 # 首尾拼接起终点高度，返回完整航迹 z_full = np.concatenate(([self.z_start], z_series, [self.z_end])) return z_full解码时注意首尾端点必须是真实起点和终点高度，这是保证航迹可用性的前提。我在第一次复现时没有注意这个问题，输出的航迹起点对不上起飞降落场高度，导致平滑后整体偏离了原始路径。
4.3 适应度函数的具体实现
适应度函数的核心代码如下：
def fitness(self, z_series): z_full = self.decode(z_series) x = np.linspace(self.x_start, self.x_end, len(z_full)) y = np.linspace(self.y_start, self.y_end, len(z_full)) # 航迹长度代价 dx = np.diff(x); dy = np.diff(y); dz = np.diff(z_full) seg_len = np.sqrt(dx**2 + dy**2 + dz**2) J_len = np.sum(seg_len) # 高度代价 J_height = np.mean((z_full - self.cruise_alt) ** 2) # 平滑代价：相邻航迹段之间夹角的惩罚 angles = [] for i in range(len(x) - 2): v1 = np.array([x[i+1]-x[i], y[i+1]-y[i], z_full[i+1]-z_full[i]]) v2 = np.array([x[i+2]-x[i+1], y[i+2]-y[i+1], z_full[i+2]-z_full[i+1]]) cos_theta = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2) + 1e-12) angles.append(np.arccos(np.clip(cos_theta, -1, 1))) J_smooth = np.mean(angles) # 威胁代价 J_threat = 0.0 for cx, cy, cr in self.threat_zones: dist = np.sqrt((x - cx)**2 + (y - cy)**2) J_threat += np.sum(1.0 / (dist[dist < cr] ** 2 + 1e-6)) return self.w1 * J_len + self.w2 * J_height + self.w3 * J_smooth + self.w4 * J_threat这里需要解释一下平滑代价的作用。无人机在实际飞行中无法瞬间改变航向，相邻航迹段的夹角越接近 180°，表示转弯越平缓。如果不加平滑项，算法会生成大量锯齿状航迹，看起来长度短，但根本无法作为可飞航迹。
4.4 主循环与改进算子的代码片段
主循环的代码片段如下：
def evolve(self): t = 0 best_fit = np.inf best_pos = None while t < self.max_iter: a = 2 * (np.cos(0.5 * np.pi * (t / self.max_iter)) ** 2) for i in range(self.n_agents): r1, r2 = np.random.rand(), np.random.rand() A = 2 * a * r1 - a C_val = 2 * r2 p = np.random.rand() if p < 0.5: # 收缩包围 D = np.abs(C_val * self.pop_best - self.pop[i]) new_pos = self.pop_best - A * D else: # 螺旋更新 l = np.random.uniform(-1, 1) D = np.abs(self.pop_best - self.pop[i]) new_pos = D * np.exp(l) * np.cos(2 * np.pi * l) + self.pop_best # 自适应惯性权重 w = 0.9 - 0.5 * (t / self.max_iter) ** 2 new_pos = w * self.pop[i] + 0.6 * new_pos new_pos = np.clip(new_pos, self.bounds[0], self.bounds[1]) self.pop[i] = new_pos # 柯西-莱维变异 if np.random.rand() < self.p_mut: self.pop_best = self.levy_cauchy_mutate(self.pop_best) # 精英保留 self.update_elite() t += 1 # 记录收敛曲线这段代码展示了非线性收敛因子、惯性权重、变异和精英保留是如何嵌入主循环的。需要注意的地方是位置更新的“新位置”不能完全脱离原位置，否则粒子会出现瞬间跳变，所以我在实现中加入了 w 和 0.6 的混合系数，让新位置在保留记忆的基础上发生偏移。
5. 复现实验设计，以及怎么判断改进“真的有效”
5.1 基准函数测试
先确定改进策略本身是否有效。我用了一组经典基准函数进行验证，包括 Sphere、Rastrigin、Ackley、Griewank。每个函数维度设置为 100，种群规模 50，最大迭代 500 次，独立运行 20 次并取平均值。结果以 Rastrigin 为例，原版鲸鱼算法的平均最优值为 12.8，改进鲸鱼算法可以降到 3.6 左右。对于 Ackley 这种多谷类函数，收敛速度和精度提升同样明显。
这类测试的意义在于把“算法改进”和“任务建模”分开评估。如果直接在三维航迹规划任务上测试，你看不出改进到底解决了“早期收敛”还是“多样性崩塌”，因为代价函数本身也会影响搜索结果。先做基准测试，再上任务测试，是复现论文时更稳妥的流程。
5.2 三维地图飞行测试与指标对比
基准测试通过后，我在自定义的三维地图上运行算法。地图范围设定为 100km × 100km，起飞点坐标和降落点坐标固定。禁飞区设了三个圆形区域，分别位于地图中部、右上和左下，半径各不相同。航迹中间点数量设为 20 个，因此维度为 20。
原版鲸鱼算法和改进鲸鱼算法的对比结果如表所示：
| 指标 | 原版WOA | 改进IWOA |
|---|---|---|
| 航迹长度 | 189.6 km | 165.4 km |
| 平均威胁代价 | 1.87 | 0.32 |
| 迭代100次后平均适应度 | 212.5 | 187.3 |
| 迭代500次后平均适应度 | 183.4 | 160.9 |
注意上表是某一次运行结果的示例，不是论文原始数据。从整体趋势看，改进算法的优势主要体现为两点：一是前 100 次迭代适应度下降速度明显更快，表明非线性收敛因子在前期释放了更多的探索能力；二是最终航迹的威胁代价大幅降低，说明柯西-莱维变异帮助算法绕开了更难检测的禁飞区边缘区域。
5.3 参数敏感性实验
复现论文的过程中，参数实验容易被忽略，但恰恰是最有参考价值的部分。我重点测了三组参数：
- 种群规模：30、50、80
- 变异概率：0.1、0.2、0.3
- 代价函数权重组合：三组不同比例
实验结论是变异概率 0.2 左右效果最佳，过大时最优解反复被扰动，难以收敛；过小时多样性不足，容易早熟。种群规模 50 和 80 的结果差异不大，但 30 时成功率明显下降。这项实验让我确定最终参数：种群规模 50，变异概率 0.2，最大迭代 500。
5.4 复现过程中的随机数种子问题
复现这类论文时，一个容易被卡住的问题是：你改进了算法，但可能某个种子下改进效果不明显。这是正常现象。群智能算法本身就有随机性，单次运行的对比不能说明问题。我的做法是固定 10 个不同随机种子，分别运行原版和改进版，记录均值和标准差，用均值下降幅度和标准差缩小幅度来说明改进效果。
另外，复现代码时要保证两个算法在同一组种子下运行。也就是原版跑种子 1，改进版也要先跑同一个种子 1，这样控制变量才能正确归因。很多复现报告里说“改进效果不稳定”，大概率是这一环节没有控制好。
6. 复现这篇论文时踩过的坑，以及给你的建议
6.1 数字地形数据的来源与插值
第一次做地形数据时，我直接用了模拟生成的栅格数据。跑通后发现有些航迹的 Z 值低于地形高程，无人机直接钻进山体了。排查之后发现是地形数据分辨率不足，航迹点的 X、Y 坐标在栅格间发生了插值偏差。解决方法是把地形高程数组做一个双线性插值，让每个航迹点都能精确地取出所在位置的高程，再用这个高程做高度约束判断。
from scipy.interpolate import RegularGridInterpolator import numpy as np # 假设 terrain_z 是二维高程数组 interp = RegularGridInterpolator((grid_x, grid_y), terrain_z, method='linear') def get_elevation(x, y): return interp([x, y])[0]这个插值函数在后续代码里几乎到处都要用，建议一开始就把高程查询封装成一个函数，而不是在每次循环里重复写插值逻辑。
6.2 航迹点数量与搜索维度对结果的影响
航迹中间点数量直接决定了搜索维度。中间点越少，航迹越粗，转弯可能太急；中间点越多，维度越大，算法收敛越慢。我试过 10、20、30 组三组中间点，20 个点是效果和效率比较平衡的选择。当你尝试 35 个点时，500 次迭代已经不够了，需要加大到 1000 次以上，这意味着计算时间接近翻倍。
如果需要在 30 个点以上运行，建议先做一个初步快速收敛（比如先用高斯过程贝叶斯优化确定一个较好的初始航迹），再把初始航迹作为鲸鱼算法的初始化种群之一，可以明显加速收敛。
6.3 论文图和代码图对不上是怎么回事
复现这类论文时，最经常遇到的挫败感是：论文里的航迹图画得特别平滑、完美地绕过所有禁飞区，但你跑出来的图路径歪歪扭扭，甚至穿过了山峰。这里有个容易被忽略的技术细节：很多论文在展示最终结果前做了一层样条平滑后处理，把离散航迹点变成了连续光滑曲线，但代价函数计算用的是平滑前还是平滑后的航迹，论文不会写清楚。
我的建议是：分两步评价，第一步，用原始航迹点计算总代价，作为算法比较的公平指标；第二步，对最优航迹做 B 样条平滑，重新评估平滑后的航迹是否仍然满足约束。如果平滑后代价显著上升，说明算法生成的航迹还不够好，需要增加平滑代价权重或提高航迹点数量。
6.4 改进算法参数“过拟合”的风险
这是我觉得整个复现过程里最值得警惕的一点：改进算法的参数很容易被调节成“只在这个测试地图上效果好”。你打开一篇论文，看到它的改进算法比原版提升了 20%，这当然可能，但更可能是它针对测试地图精心调参了。
我在复现时特意做了“换地图测试”：生成三张完全不同的随机地形图，不针对某一张调参，而是用同一组参数在每张地图上各跑 20 次。最后发现，改进算法在三张地图上都能稳定优于原版，但优势幅度有所不同。这个跨场景稳定性比单张地图上的提升幅度更有说服力。复现论文时，不要只看最终图形效果，而是要看算法在多张地图上的平均表现和误差范围。
6.5 计算效率与代码加速建议
三维航迹规划的计算瓶颈主要在于适应度函数里的距离计算和插值。如果航迹点有 20 个，种群 50 个，迭代 500 次，适应度函数的计算次数大约是 2.5 万次，每次做循环计算长度和平滑代价，总时间分钟级别。但如果你把种群加到 100，航迹点加到 50，总计算次数会超过 25 万次，耗时可能达到小时级别。
推荐两个优化方向：一是使用 numpy 向量化计算，把种群中所有个体的距离计算批量处理，而不是单个个体循环；二是用 joblib 的并行加速为每个个体分配一个线程。代码在台式机上可以把速度提升约 3～4 倍，在笔记本上也有接近 2 倍的提升。
from joblib import Parallel, delayed def evaluate_one(individual): return fitness(individual) fitness_values = Parallel(n_jobs=-1)( delayed(evaluate_one)(self.pop[i]) for i in range(self.n_agents) )需要注意，并行计算时不能共享可变状态，否则线程之间会互相覆盖数据。我只是让每个线程只读种群数据和地形数据，输出适应度值，再汇总到主进程。
7. 最后分享一点个人体会
整个复现过程花了大概两周时间，前面大部分时间都用在了建模和调试上，真正写改进算法代码的时间其实很短。回头来看，如果想快速把这个方案跑起来，你最应该花时间的地方是代价函数的设计和约束校验，而不是急着改算法公式。算法改进往往是锦上添花，一个合理的代价函数才是决定航迹是否可飞的关键。
如果你准备尝试这个题目，我建议先跑通一个只有长度代价的基础版本，确认航迹能够连起来，再逐步加入高度、平滑和威胁代价。每加一项，就观察最优航迹是否发生变化。这样一步步迭代，你会很清楚每一项代价对结果的影响，而不是等最后拼在一起时才发现轨迹完全不可行。
最后一个小技巧：把每次实验的算法参数、地图参数、种子和收敛曲线都存在一个 JSON 文件里。做对比实验时，你会遇到大量参数组合，光靠记忆完全不可靠，存成文件会让你后续复现和修改省很多力气。
