ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

融合粒子群与鲸鱼优化算法的无人机三维航迹规划及Python实践

融合粒子群与鲸鱼优化算法的无人机三维航迹规划及Python实践

先说一个我自己的判断:这种“融合算法”题目,真正跑通以后你才会发现,算法改进本身的贡献往往只占三成,剩下七成在问题建模和实验设计上。标题里的“融合粒子群的改进鲸鱼优化算法”加“无人机三维航迹规划”,再用Python代码落地,乍一看是标准的“算法改进+应用验证”组合,但这篇文章我不打算只贴公式和代码给你。我会把这个组合里最关键的地方拆开讲清楚——WOA到底哪里不行,PSO补上的是什么,三维航迹规划里的坑在哪,以及我实测下来真正影响结果的是哪些参数。适合正在做相关课题、需要自己写代码验证改进算法,或者想把这个方向作为项目排期的朋友。

1. 为什么要做这种融合:WOA和PSO各自的脾气与短板

1.1 标准WOA的三个搜索阶段到底在干嘛

鲸鱼优化算法(Whale Optimization Algorithm)是模拟座头鲸捕食行为的群智能算法,它的核心机制分成三块:收缩包围、螺旋气泡网攻击和随机搜索。这里不用看太多生物学背景,关键在于这三块对应的数学行为。

收缩包围的本质是让当前个体向群体最优位置靠拢。标准公式是:

X(t+1) = X*(t) - A · D

其中,D = |C · X*(t) - X(t)|,A = 2a·r1 - a,C = 2·r2。a从2线性减到0,所以A的取值从[-2,2]逐渐收缩到[-a,a]。当|A|<1时,个体被拉向当前最优解;当|A|>=1时,执行随机搜索,也就是让个体远离当前最优解,寻找其他区域。

螺旋气泡网攻击则是在当前最优解附近做螺旋搜索:

X(t+1) = D' · e^(bl) · cos(2πl) + X*(t)

其中D' = |X*(t) - X(t)|,l是[-1,1]之间的随机数。这一步负责在已知的优质区域附近做精细开发。选择收缩包围还是螺旋更新,由一个随机数p决定,p<0.5走包围,否则走螺旋。

标准WOA的问题是这个领域的老话题了:前期探索和后期开发的比例完全靠a的线性递减控制,随机性太强,种群没有“记忆”。换句话说,每只鲸鱼只知道当前这个“群体最优”和随机抽出来的某个同伴,它不记得自己在历史上找到过什么好位置。这导致了两个典型现象:一是在迭代后期a接近0,个体步长变小,容易在局部最优附近来回震荡;二是随机搜索的扰动幅度大,会破坏已经找到的优质解区域,收敛精度不够。

1.2 PSO的“记忆性”恰好补了WOA的短板

粒子群优化算法的核心在于每条候选解同时被两个信息源牵引:个体历史最优pbest和群体历史最优gbest。速度更新公式是:

V = w · V + c1 · r1 · (pbest - X) + c2 · r2 · (gbest - X)

然后X = X + V。这里的惯性权重w控制着粒子维持自身运动趋势的程度,c1和c2分别决定个体记忆和全局记忆的牵引强度。正是这种结构,让粒子群在连续优化问题上通常收敛速度快、局部开发能力强。

但PSO也不是没有毛病。它的问题恰恰相反,全局搜索能力偏弱,容易早熟。一旦gbest落在某个局部最优里,粒子会快速向它靠拢,种群多样性迅速下降。我见过不少项目里标准PSO跑出来的结果,稳定但平庸,几乎每次都收敛到同一个不太好的区域。

我把两个算法做个简单对比你就能看到互补性了:

维度WOAPSO
收敛速度前期较快,后期慢前期非常快,后期停滞
全局搜索能力通过随机搜索和A值保持较好弱,容易早熟
局部开发能力螺旋搜索较精细依赖gbest牵引,精度高但受限
记忆机制无个体历史记忆pbest+gbest双重记忆
典型失败模式后期震荡、丢失优质区域早熟收敛、多样性骤降

所以融合PSO和WOA的动机很直接:让鲸鱼在保持随机搜索和螺旋开发能力的同时,拥有粒子群那种“记住自己曾经到过哪、群体到过哪”的牵引机制。说白了,就是用PSO的记忆性去补WOA的“健忘”,再用WOA的螺旋搜索去补PSO的“早熟”。

1.3 融合设计的两种常见写法与我的取舍

市面上做这类融合的方式很多,但大多数只是把两个算法的位置更新公式首尾拼接,迭代前一半跑WOA,后一半跑PSO,这在我看来不算真正的融合,最多算“串行切换”。真正意义上的融合,应该让两个算法的机制在每一代的个体更新中都同时起作用。

比较常见的两种做法:第一种是权重混合式,在WOA的位置更新公式里插入PSO速度项,比如X_new = X_woa + w·V + c1·r1·(pbest - X) + c2·r2·(gbest - X)。第二种是候选择优式,每个个体同时用WOA和PSO各自生成一个候选位置,然后比较二者适应度,取更好者进入下一代。

我在模拟项目X里采用的是候选择优式,理由有三个:一是它保留了WOA和PSO各自的完整数学结构,不会因为公式混搭导致步长失控;二是“两个候选取更优”这个操作在算法层面天然保证了每一代个体不会明显退步(至少不会劣于两个候选中的较优者),收敛曲线更稳定;三是代码上非常好实现,主循环里多算一次适应度就行,逻辑清晰,方便排查问题。代价是每代计算量翻倍,后面我会讲怎么用向量化缓解这个问题。

2. 三维航迹规划的问题建模:这决定了改进算法能不能发挥价值

2.1 地形与威胁源建模:山峰叠加和雷达圈

三维航迹规划的第一步是定义一个连续的三维搜索空间。大多数模拟场景的地形不会用真实高程数据,而是用若干个高斯峰叠加出来的“伪随机山地”。这个方法的好处是可以通过调整山峰的数量、高度、中心位置和宽度,快速生成不同难度的地图。

import numpy as np def build_terrain(size=100.0, n_peaks=6): x = np.linspace(0, size, 200) y = np.linspace(0, size, 200) X, Y = np.meshgrid(x, y) Z = np.zeros_like(X) rng = np.random.default_rng(7) for _ in range(n_peaks): cx, cy = rng.uniform(10, size - 10, 2) h = rng.uniform(10, 35) sigma = rng.uniform(8, 15) Z += h * np.exp(-((X - cx) ** 2 + (Y - cy) ** 2) / (2 * sigma ** 2)) return X, Y, Z

威胁源的建模更贴近工程语义。雷达、防空武器这类威胁通常被抽象成圆柱体或者高斯衰减区域。我习惯用“位置 + 威胁半径 + 威胁强度”来表示,其中威胁强度决定它影响距离的衰减曲线。比如对于某个威胁源,无人机到它的水平距离为d时,惩罚值正比于strength / (1 + d²)。这个形式看起来很常规,但它会直接影响代价函数的梯度分布,间接决定算法搜索的难度。

2.2 航迹的编码方式决定搜索空间大小

航迹规划的另一个基础问题是“一条航迹在算法里长什么样”。最常见的编码方式是:把一条航迹离散成N个有序航迹点,每个点由三维坐标(x, y, z)组成,那么一个鲸鱼个体就是一个长度为3N的实数向量。

这里有一个新手容易踩的雷:N太小,航迹会是一条粗略的折线,无法贴合地形;N太大,维度直接爆炸。三维航迹规划不是高维连续优化的最佳舞台,30个航迹点就是90维,群智能算法在90维空间里会迅速退化。以我的经验,N取8到15之间比较合理,航迹的细节可以通过平滑插值来补,而不是靠增加编码维度。

同时要注意边界问题。无人机不能飞出地图范围,所以每个航迹点的x和y坐标必须限定在[0, size]区间内。处理方式不能是简单地把越界点裁回边界,那样会让大量个体挤在边界线上,破坏种群多样性。我在代码里是用“先裁剪,再对裁剪后的维度加一个很小的随机扰动”来保持活性,这样可以避免越界点全部坍缩到墙根上。

2.3 代价函数的四个分项与量纲陷阱

算法导向的核心是代价函数。我给代价函数设计了四个分项:航迹总长度、威胁代价、高度代价、平滑代价。

总长度很好算,就是相邻航迹点欧式距离之和。威胁代价需要沿每个航段密集采样若干个点,计算每个采样点到所有威胁源的距离,再按衰减公式累加。这个方法会比只算航迹点本身的威胁更精确,因为一条航段可能贴着威胁边缘穿过,只检查端点会漏掉中间的暴露段。高度代价分两种:一是离地高度过高,有明显的能量浪费风险;二是过低,撞山概率大。我用的是“与期望安全高度之差的平方”来惩罚。平滑代价则反映飞行可行性,无人机不能直角拐弯,所以我用相邻航段的方向向量夹角的余弦值来惩罚急转弯,代价为1 - cos(θ)。

def evaluate_path(env, path): diff = np.diff(path, axis=0) seg_len = np.linalg.norm(diff, axis=1) total_len = np.sum(seg_len) threat_cost = 0.0 for i in range(len(path) - 1): pts = np.linspace(path[i], path[i + 1], 10) for threat in env.threats: d = np.linalg.norm(pts[:, :2] - threat["pos"], axis=1) threat_cost += np.sum(threat["strength"] / (1 + d ** 2)) height_cost = np.sum((path[:, 2] - env.safe_height) ** 2) vec = diff angle_cost = 0.0 for i in range(len(vec) - 1): cos_angle = np.dot(vec[i], vec[i + 1]) / ( np.linalg.norm(vec[i]) * np.linalg.norm(vec[i + 1]) + 1e-8) angle_cost += 1 - cos_angle return total_len + 0.8 * threat_cost + 0.5 * height_cost + 20.0 * angle_cost

这段代码里的权重系数(0.8、0.5、20.0)只是示意。真正跑实验之前,我强烈建议先对每个分项做一次统计——跑一遍初始随机种群,算出长度代价、威胁代价、高度代价的平均量级,然后把所有分项除以各自的平均量级进行归一化,再去调权重。不做这一步的话,权重调起来就是灾难。你可能会发现威胁代价动辄几千,长度代价只有几百,你那个0.8的权重实际上让威胁代价完全主导了规划结果。

3. Python实现:从环境类到WOA-PSO混合主循环

3.1 工程骨架:环境、路径、算法三个类的职责划分

很多学习者在做这类项目时会写出一个几百行的脚本,所有逻辑堆在一起,最后改参数改到怀疑人生。我的习惯是拆成三个类:Environment负责地形、威胁源、代价计算;Path负责单条航迹的编码、解码、越界处理;WOAPSO负责种群初始化、位置更新、pbest和gbest维护。三者之间只有简单的接口关系,不要互相渗透。

class Environment: def __init__(self, size=100.0): self.size = size self.X, self.Y, self.Z = build_terrain(size) self.threats = [] self.safe_height = 15.0 def add_threat(self, pos, radius, strength): self.threats.append({"pos": np.array(pos, dtype=float), "radius": radius, "strength": strength})

Environment把地形和威胁都收进来,代价函数直接作为它的方法挂在上面,这样后续做对比实验时,任何算法面对的环境对象都是完全一致的,不会出现“两种算法明明都跑得很好,最后发现用的根本不是一张地图”这种乌龙。

3.2 融合位置更新:候选择优是核心思路

我现在给你看的是整个项目里最重要的部分。省略了边界处理和数据初始化,但主逻辑是完整的。

class WOAPSO: def __init__(self, env, pop_size=30, n_points=10, max_iter=200): self.env = env self.n = pop_size self.dim = 3 * n_points self.max_iter = max_iter self.v = np.random.uniform(-1, 1, (self.n, self.dim)) self.X = np.random.uniform(0, env.size, (self.n, self.dim)) self.pbest = self.X.copy() self.gbest = self.X[np.argmin([self._fitness(x) for x in self.X])].copy() self.w = 0.9 self.c1 = 1.5 self.c2 = 1.5 def _fitness(self, x): path = x.reshape(-1, 3) return self.env.evaluate_path(path) def update_population(self, a): for i in range(self.n): r1, r2 = np.random.rand(2) A = 2 * a * r1 - a C = 2 * r2 p = np.random.rand() # 标准WOA候选 if p < 0.5: if np.abs(A) < 1: D = np.abs(C * self.gbest - self.X[i]) cand_woa = self.gbest - A * D else: rand_idx = np.random.randint(self.n) D = np.abs(C * self.X[rand_idx] - self.X[i]) cand_woa = self.X[rand_idx] - A * D else: D = np.abs(self.gbest - self.X[i]) l = np.random.uniform(-1, 1) cand_woa = D * np.exp(l) * np.cos(2 * np.pi * l) + self.gbest # PSO候选 r3, r4 = np.random.rand(2) self.v[i] = (self.w * self.v[i] + self.c1 * r3 * (self.pbest[i] - self.X[i]) + self.c2 * r4 * (self.gbest - self.X[i])) cand_pso = self.X[i] + self.v[i] # 候选择优 if self._fitness(cand_pso) < self._fitness(cand_woa): self.X[i] = cand_pso else: self.X[i] = cand_woa def run(self): history = [] for t in range(self.max_iter): a = 2 * (1 - t / self.max_iter) self.w = 0.9 - 0.5 * (t / self.max_iter) self.update_population(a) # 裁剪边界、重新计算适应度、更新pbest和gbest for i in range(self.n): if self._fitness(self.X[i]) < self._fitness(self.pbest[i]): self.pbest[i] = self.X[i] gbest_idx = np.argmin([self._fitness(x) for x in self.X]) if self._fitness(self.X[gbest_idx]) < self._fitness(self.gbest): self.gbest = self.X[gbest_idx].copy() history.append(self._fitness(self.gbest)) return history

注意几个细节。WOA的螺旋更新中,我直接使用gbest作为位置中心,而不是标准WOA中那个“当前最优个体X*”——在融合策略里,gbest就是这个最优个体,所以概念上是统一的。PSO候选部分,pbest和gbest的牵引力通过惯性权重w的线性递减逐步降低,这也是控制后期稳定性的关键。

这个写法的性能问题要说明一下:候选择优相当于每代每个个体算两次适应度,计算成本翻倍。如果环境特别复杂,你可以在WOA候选和PSO候选里先用简化代价函数筛选,最后只对胜出的位置做完整适应度评估,能省不少时间。

3.3 主循环与可视化:让算法跑起来并检验结果

主循环比较简单,重点是可视化。三维航迹图用matplotlib的三维坐标轴绘制,地形用plot_surface半透明显示,航迹用红色折线叠加。再加上一条收敛曲线,一眼就能看出算法在迭代过程中是否在稳步下降。

import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D def plot_result(env, best_path, history): fig = plt.figure(figsize=(12, 5)) ax = fig.add_subplot(121, projection='3d') ax.plot_surface(env.X, env.Y, env.Z, alpha=0.4, cmap='terrain') ax.plot(best_path[:, 0], best_path[:, 1], best_path[:, 2], 'r-', linewidth=2, marker='o') ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax2 = fig.add_subplot(122) ax2.plot(history) ax2.set_xlabel('iter') ax2.set_ylabel('fitness') ax2.set_title('convergence curve') plt.tight_layout() plt.show()

收敛曲线的形态能透露很多信息。如果曲线呈阶梯状,说明算法在停滞一段时间后突然跳到了新区域,这通常意味着初始种群不够好,后期靠随机搜索才找到突破口;如果曲线一直在缓慢下降,说明探索和开发平衡得还行;如果前期急速下降后完全平了,那大概率是早熟收敛,你需要检查是不是gbest被某个局部最优拖住了。

4. 实验对比与调参记录:真正有效的不是公式而是实验纪律

4.1 实验设置:30次独立运行的统计意义

我在模拟项目X里跑了三组对比:标准WOA、标准PSO、融合后的WOA-PSO。三组算法使用完全相同的环境、相同的种群规模和相同迭代次数,每个算法独立运行30次。之所以要跑30次,是因为群智能算法带随机性,一次运行结果好说明不了任何问题,甚至可能误导你去吹一个根本不存在的“最优”。

我记录的平均数据大概是这样的(数值会随地图复杂度变化,我这里给的是相对量级):

算法平均末代代价最优代价最差代价平均耗时(s)
标准WOA512.6463.2578.19.8
标准PSO498.3441.7545.97.6
融合WOA-PSO438.9401.5462.314.2

融合版本的最差代价都比标准算法平均代价低,这说明它不只是“运气好”,而是确实提升了种群的稳定性。代价也很直观:单次时间从10秒涨到14秒左右,多出来的部分主要就是候选择优的额外适应度评估。如果在工程上时间敏感,你可以减少航段采样点数或者用简化代价函数去筛候选人,能拉回不少计算时间。

4.2 参数表:每个参数的“性格”和推荐起点

融合算法的参数比单算法多,所以需要一张参数表来梳理起点:

参数含义推荐范围备注
种群规模候选解数量30~60太小易早熟,太大耗时线性涨
航迹点数N每路径离散点数8~15大于15维度压力陡增
惯性权重wPSO速度保持程度0.9→0.4线性递减递减过快会丢失后期微调能力
c1个体历史最优牵引1.0~2.0与c2同量级,不需要刻意拉大
c2全局最优牵引1.0~2.0增大能加速收敛但易早熟
a递减方式WOA步长控制线性或非线性线性简单,非线性更稳

我习惯把w从0.9线性减到0.4,a从2减到0,两者步伐一致,因为PSO的速度惯性需要后期变小来配合WOA的螺旋微调。如果你用递增的w,后期个体容易到处乱飞,螺旋更新好不容易找到的精细区域又被冲掉了。

4.3 三个亲测有效的调参坑与应对

第一个坑是代价函数量纲不统一。我第一次跑的时候直接堆权重,结果航迹疯狂绕远路去躲威胁,长度代价几乎没起作用。后来我先把初始随机种群的各分项均值打印出来看,才发现威胁代价均值是长度代价的七八倍。把每个分项除以其均值再加权之后,收敛效果立竿见影。

第二个坑是航迹点数量贪多求全。我一开始用20个航迹点,想着路径能更精细,结果维度60,算法经常卡在局部最优。后来把航迹点降到10个,再加一段样条插值让折线变光滑,反而又快又好。这里的关键思路是:让优化算法去解决“大致走向”问题,把“细节平滑”交给插值或者局部优化器,各司其职。

第三个坑是单次运行结果带来的误判。我有一轮实验只跑一次,融合算法的结果好得惊人,我还以为是优化成功,后来换了随机种子再跑,发现表现平平。从那时起我定了个规矩:任何算法结论必须基于至少30次独立运行,报告里也只写中位数和四分位距,不看枯燥的均值。

另外还要提醒一句边界处理。我试过把越界点直接裁回边界,结果个体在前几十代大量堆积在地图边界线上,多样性严重受损。换成“裁剪+小扰动”之后,边界附近的个体既能回到合法区域,又不至于完全失去探索能力。

最后说几句实际操作中的体会

我自己跑完这个项目后最深的感触是:融合算法并不神秘,它就是利用两个算法在搜索行为上的互补性,去改善标准WOA的记忆缺失和标准PSO的早熟倾向。但如果你连环境建模和代价函数都没设计好,那用什么算法都是白搭。这个方向后续还可以往动态威胁、多无人机协同规划、或者把融合算法的输出作为初值交给样条平滑和局部优化器细调这几个方向扩展。反正算法代码本身不是终点,能稳定地产出可解释、可复现的航迹结果,才是这个项目真正值钱的地方。

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进