☰
城市场景无人机三维路径规划:NMOPSO多目标粒子群算法实现解析
2026/10/21 8:50:35 网站建设 项目流程

拿到“基于城市场景下无人机三维路径规划的导航变量的多目标粒子群优化算法NMOPSO研究(Matlab代码实现)”这个题目,很多人可能第一眼觉得信息量很大,实际上拆分下来就三件事:城市环境怎么建、路径怎么编码、多目标怎么搜。我以前做二维路径规划时,只要绕开楼层就够了;但三维城市场景里,一架无人机要去执行任务,光给一串航路点根本不够,还要考虑高度层、障碍物形状、机动力约束、任务时间、整段能耗,这些约束和目标放在一起,传统单目标寻优完全招架不住。NMOPSO就是在这种背景下被拿出来的方案,它在粒子群优化框架里引入非支配排序、外部存档、拥挤距离修剪,最后输出的是整条Pareto前沿,而不是孤零零的一条“最优路线”。

这篇文章不会写成论文综述,我会以实际做实验的视角,把你需要知道的环境建模思路、NMOPSO核心机制、Matlab实现结构和调参经验全部串起来,尽量让新手能照着落盘,也让已经有基础的人能少踩几个坑。

1. 项目为什么难:城市场景下的三维路径规划到底在优化什么

1.1 无人机在城里飞,先过环境建模这一关

城市环境跟野外山区最大的区别是“人造障碍物密集且形状规整”。高层建筑、商超裙楼、信号塔、跨江桥梁,这些障碍物在三维空间里可以被简化成不同高度的垂直棱柱或长方体。实际工程里,我们不需要像测绘那样精细到每块墙体,而是按飞行安全缓冲半径把障碍物膨胀一圈。

Matlab里做这种建模,早期我习惯用栅格地图,把地图切成一格一格的立方体voxel。缺点是分辨率低时路径很“锯齿”,分辨率高了内存爆炸。后来改成了多边形体素加碰撞检测函数:每个建筑物用中心坐标、底边长度、顶高三个参数描述,碰撞检测时就判断路径点是否落在任意建筑物的膨胀包围盒内。

% 障碍物定义示例:格式 [x_center, y_center, half_width, half_depth, z_min, z_max] obs = [ 20, 30, 8, 6, 0, 45; 45, 60, 12, 8, 0, 80; 70, 20, 10, 10, 0, 35; ];

这种方式比体素地图好用得多,因为三维空间里碰撞判断的次数跟路径点数量线性相关,而不是跟体素数相关。路径规划迭代几百步、每步评估几百个粒子时,这个差异会被放大得特别明显。

1.2 导航变量的提法:不是只给几个路径点就行

题目里反复强调“导航变量”,这个词比“路径点”准确得多。一条无人机三维路径,如果只看位置,就是一串(x, y, z);但如果要交给飞控执行,还得有偏航角ψ、俯仰角θ、飞行速度v、各航段的时间间隔Δt。把这些变量合在一起,才是完整的导航变量集。

NMOPSO里每个粒子的位置向量就是一条完整候选路径的编码。假设路径由N个航路点组成,每个航路点包含x、y、z、ψ四个核心量,那么粒子的维度就是4N。再加上速度v或时间间隔,维度会上升到5N甚至6N。这里就引出“高维”问题——粒子群算法在处理这种几十甚至上百维的决策变量时,收敛速度和多样性之间的平衡会变得很微妙。

我在项目中采用的做法是:路径点数量固定为10个,两端分别固定到起点和终点,中间8个点为待优化节点。这样粒子维度是8×3=24维,如果要体现具体航向,再额外加上8个偏航角自由度,总维度32维。这个规模对MOPSO而言算是中等偏上,既能体现高维多目标难度,又不至于让计算时间失控。

1.3 高维多目标的“高维”体现在哪里有

所谓“多目标优化”,至少两个目标以上,比如路径长度和威胁度。但无人机在城市里飞,实用场景下目标远不止这两个。我实际建模时至少用了6个目标,这也是“高维”的一个体现:

  • 总路径长度:直接影响飞行时间和能耗。
  • 障碍物安全裕度:路径点离建筑物边界的最近距离,越远越安全。
  • 高度变化代价:频繁爬升和下降非常耗电,而且机体姿态会剧烈变化。
  • 转弯角度代价:相邻三段路径的夹角之和,过大的转弯角说明航路“拧巴”。
  • 航向一致性代价:无人机与参考巡航方向之间的偏差,便于维持稳定飞行。
  • 飞行总时间估算:每段距离除以该段速度后累加,发电机限制下的时间约束。

多目标粒子群算法要做的事情就是把所有目标拉进同一个优化框架,寻找一组不能被彼此支配的解。所谓“支配”,就是解A在所有目标上都不差于解B,并且至少有一个目标严格优于B,那A就支配B。最终保留的是互不支配的Pareto解集。

2. NMOPSO核心机制拆解:从标准PSO到非支配多目标粒子群

2.1 从PSO到MOPSO:粒子怎么同时记住多个最优

标准粒子群里,每个粒子有两个“记忆”:个体历史最优pbest和全局最优gbest。速度更新公式就一行:

v(i,:) = w * v(i,:) + c1*r1*(pbest(i,:) - x(i,:)) + c2*r2*(gbest(i,:) - x(i,:)); x(i,:) = x(i,:) + v(i,:);

单目标下,gbest就是历史上适应度最好的那一个,pbest也是自己遇到的最好适应度。但多目标下,“最好”不再是一个标量。解A可能在路径最短上优于解B,解B在安全裕度上又优于解A。这时如果只把某个单一目标当作比较标准,就会把问题退化成加权单目标,失去多目标意义。

所以MOPSO要做两件事:第一,每个粒子的pbest不再是一路追着单个目标走,而是记录自己在非支配关系下无法被淘汰的偏好解;第二,全局最优gbest不再只有一个,而是从Pareto最优解集里动态选出引导粒子。

2.2 NMOPSO的改进点:非支配排序、外部存档、拥挤距离

把NMOPSO拆开看,核心就是三个机制:

  • 非支配排序:种群中所有粒子两两比较支配关系,将不被任何粒子支配的个体划到第一层,去掉第一层后再划出第二层,直到全部分层完毕。第一层就构成当前Pareto前沿。
  • 外部存档:单独维护一个数组archive,用来存放迄今为止发现的最优非支配解。每代更新时,把当前种群的Pareto前沿合并进存档,再对存档做一次非支配排序,剔除掉被支配的旧解。
  • 拥挤距离修剪:如果存档中非支配解数量超过设定上限,需要删除一部分,否则迭代几百代后存档会膨胀。标准做法是计算每个解在目标空间里的拥挤距离,优先删除距离小、周围解扎堆的个体,这样能在前沿上保留均匀分布的端点解。

这里我要多说一句:为什么常用拥挤距离而不用网格法。网格法把目标空间划分成固定大小的格子,一个格子内只保留一个解,这在高维多目标场景下会遇到“格子爆炸”问题,目标维度增加到6个以上时,哪怕每个维度只切10等分,网格数也是10的6次方,绝大多数格子是空的。拥挤距离只依赖当前目标空间样本的倒角,计算复杂度跟种群规模相关,不会随维度指数膨胀。

2.3 为什么在高维目标下,拥挤距离比档案网格更稳

实际测试中,我用6目标模型跑过对比,网格法在后半程会产生很多空白区域,导致Pareto前沿出现断层;拥挤距离修剪因为是对已有个体计算最近邻欧氏距离,所以不会出现假空白。

但拥挤距离也有缺陷:在目标维度特别高时,欧氏距离会被高维稀疏性稀释,所有点的距离都比较接近,修剪起来区分度下降。解决办法不是扔掉拥挤距离,而是配合支配强度函数修饰距离——如果某个局部区域已经极密集,给该区域的支配强度加惩罚权重。NMOPSO里这种策略属于“加权拥挤距离”,比原始版本稳定不少。

还有一个细节容易被跳过:外部存档的引导选择。gbest不能只从存档里随机挑,否则会丢失“拉开前沿”的能力。更稳妥的做法是先计算每个存档个体被选中的概率,概率与其稀疏程度成正比。稀疏位置的解能牵引粒子往尚未覆盖的方向飞。简化的选择实现如下:

% 按拥挤距离从大到小排序,排在前面的有更高概率成为gbest [~, idx] = sort(CrowdDist, 'descend'); candidate = arc_idx(idx(1:max(2, round(length(idx)*0.1)))); gbest = archive_individual(candidate(randi(length(candidate))));

3. Matlab实现主体框架与关键代码

3.1 总体文件结构与数据流

我习惯将项目拆成以下几块:

  • main_runner.m:程序入口,读入地图和任务参数,循环调用PSO更新。
  • initSwarm.m:初始化种群和外部存档。
  • evaluate.m:计算每个粒子的6个目标值,并做碰撞约束检查。
  • updateArchive.m:合并非支配解,修剪存档。
  • updateParticle.m:速度、位置更新,pbest和gbest选择。
  • plotResults.m:绘制路径轨迹与Pareto前沿图。

整个数据流是:initSwarm生成初始粒子 → evaluate计算目标 → updateArchive整理Pareto集合 → updateParticle更新位置 → 循环直到迭代结束 → plotResults输出结果。

注意,evaluate是整个程序里最耗时的部分。因为每个粒子都是一条完整路径,路径上的任何一段碰到障碍物,整个粒子都要被判罚。如果障碍物检测写得低效,600代×60个粒子的实验要跑几个小时。

3.2 个体编码与初始化

采用“两头固定、中间可变”的编码方式。设起点为S,终点为G,中间节点按比例插入8个变量点。粒子向量X的维度是24维,对应8个点的x、y、z坐标。

dim = num_inner_points * 3; lb = [x_min * ones(1,N), y_min * ones(1,N), z_min * ones(1,N)]; ub = [x_max * ones(1,N), y_max * ones(1,N), z_max * ones(1,N)]; pop = lb + rand(pop_size, dim) .* (ub - lb);

初始化时不能完全随机,至少要保证起点和终点之间的直线走廊附近有足够初始样本,否则第一步迭代很多粒子都在障碍物内部,惩罚值完全一样,导致算法前期没有任何梯度信息。我特别加了短距离偏移初始化:

for i = 1:pop_size for j = 1:num_inner_points ratio = j / (num_inner_points + 1); base = S + ratio * (G - S); pop(i, j) = base(1) + randn * init_std; pop(i, j+N) = base(2) + randn * init_std; pop(i, j+2N) = base(3) + randn * init_std; end end

init_std设为地图短边尺寸的10%。这样初始种群既均匀覆盖了走廊附近,又不至于全部叠在一条直线上。

3.3 目标函数与约束处理

目标函数这块我建议写成一个独立的函数,返回一个行向量f,并单独返回violate标记。

function [f, feasible] = evaluateOnce(x, S, G, obs, params) pts = decodePath(x, S, G); f1 = totalPathLength(pts); f2 = safetyMargin(pts, obs); f3 = climbingCost(pts, params.max_climb); f4 = smoothnessCost(pts); f5 = cruiseHeightChange(pts, params.ref_height); f6 = timeEstimate(pts, params.speed_profile); f = [f1, f2, f3, f4, f5, f6]; feasible = checkFeasible(pts, obs, params.min_turn_radius); end

约束的处理不要直接丢给惩罚函数了事。我做实验时发现,惩罚函数太容易把Pareto前沿的多样性毁掉,因为不约束的粒子会把大量目标值拉得极差,导致Pareto前沿收敛到某个不切实际的角落。推荐的方案是“硬约束+修复机制”:如果粒子路径穿过禁飞区,就给该粒子全部目标值都加上一个远超正常范围的惩罚向量,并在下一步强制将其拉回可行域附近。约束越严格,修复越重要。

3.4 主迭代:速度更新、档案维护、选择引导

主迭代主体不算复杂,但坑很多。精简版本如下:

for iter = 1:max_iter for i = 1:pop_size % 选择gbest:从存档中按拥挤度权重随机选 gbest = chooseGuide(archive, crowd_dist); % 更新速度与位置 r1 = rand(1, dim); r2 = rand(1, dim); if rand < 0.8 v(i,:) = w*v(i,:) + c1*r1.*(pbest_x(i,:)-pop(i,:)) ... + c2*r2.*(gbest_x - pop(i,:)); pop(i,:) = pop(i,:) + v(i,:); else pop(i,:) = mutateByPolynomial(pop(i,:), lb, ub, iter/max_iter); end end % 边界处理:越界分量拉回 pop = clipToBound(pop, lb, ub); % 评估 for i = 1:pop_size [f(i,:), feas(i,:)] = evaluate(pop(i,:)); end % 合并归档、更新pbest archive = updateArchive(archive, pop, f, feas); [pbest_x, pbest_f] = updatePbest(pbest_x, pbest_f, pop, f, feas); end

这里我用了20%概率的多项式变异来保证后期多样性。很多人把MOPSO做成纯“粒子群+Shuffled”,实际跑下来会发现粒子在连续解空间里极为容易聚集到单个Pareto点。

关于惯性权重w,我用的是线性递减策略,初期0.9,后期0.4。这种策略在多目标下比常数w稳定,因为前期需要大范围搜索,后期需要小范围局部精调。

3.5 路径还原与可视化

Pareto前沿上每个解都是一组航路点,还原路径时要做一次平滑插值。我用三次样条插值把稀疏航路点加密到每0.5米一个点,再次检查是否穿障碍物。这一步是为了确保最终结果能被飞控用于轨迹跟踪。

可视化用Matlab的plot3画路径线,用patch画障碍物立体包围盒,再用scatter3输出Pareto前沿的分布。

figure; for k = 1:length(selected_path) p = pts_cell{k}; plot3(p(:,1), p(:,2), p(:,3), 'LineWidth', 1.2); hold on; end patch(obs_vertices, 'FaceColor', [0.6 0.6 0.6], 'EdgeColor', 'k'); scatter3(pf(:, 1), pf(:, 2), pf(:, 3), 40, 'filled');

这个图我一般要出两张,一张是Pareto前沿三维投影,一张是城市环境中所有非劣路径叠加显示。后者对判断“解是否具有工程可用性”特别重要,很多解虽然数学上非劣,实际飞起来却会贴着楼顶旋翼半径范围内擦过,工程上不可用。

4. 参数标定与实验设计:一版能跑的配置和评估方法

4.1 种群规模、迭代次数、三维空间分辨率怎么定

参数这块没有万能答案,但有一个合理起点。我用过的典型配置如下:

参数取值说明
种群规模60太小则Pareto前沿覆盖不足,太大则计算开销过大
最大迭代2006目标场景200代能基本稳定,继续增加边际收益低
外部存档上限100与目标维度相关,6目标时100个解已足够描述前沿
惯性权重0.9到0.4线性递减前期探索,后期开发
c1/c21.8/2.2稍微偏重社会引导,让种群更快向存档区聚拢
变异概率0.2高维度下保持多样性
网格分辨率5米障碍物膨胀后5米大小比较均衡

空间分辨率不是越细越好,5米栅格意味着航路点的最小有效移动粒度是5米。如果设成1米,粒子维度不变,但障碍物检测次数增加5~10倍,实验时间成倍增长,而路径质量的提升往往并不明显,除非任务窗口极其狭窄。

4.2 城市场景障碍物建模参数

城市区域我用的是模拟的“商务区网格”:地图为1000米×800米×120米的空域,其中均匀分布若干10米到80米高的长方体建筑。障碍物设置需要遵循两个原则:一是起点和终点之间至少存在一条经过人工预检查的可行通道;二是建筑密度不能过低,否则三维路径退化成二维。

具体来说,我在场景里放了18个障碍物,分成低层裙楼、高层塔楼、高架桥三类。高层塔楼的高度达到80米,导致无人机必须限制在特定的高度走廊内飞行。这个场景对NMOPSO的挑战在于:不同高度层的可行路径长度和安全裕度是强冲突的,能同时凸显算法在高维多目标下的对抗能力。

4.3 实验结果怎么评:超体积、IGD、运行时间

光看Pareto前沿图画得漂亮不够,还要量化比较。我这边常用三个指标。

超体积(Hypervolume, HV):Pareto前沿与参考点围成的目标空间体积,HV越高越好。参考点怎么选?取每个目标的最大值或稍微超过最大值的点。注意各目标要先归一化到[0,1]或[0,10],否则量级大的目标会主导HV计算。

反世代距离(IGD):先构造一组参考前沿点,然后计算算法得到的Pareto前沿离参考前沿的平均最小距离,越低越好。没有标准参考前沿时,可以用多次独立运行合并后的Pareto前沿作为近似参考。

运行时间:不可忽略。NMOPSO每个迭代里要求更新存档并做非支配排序,复杂度约O(MN²)。我给两组对比实验做过计时:种群60、迭代200,算6目标24维,单次运行耗时在Matlab里大约90秒左右;如果把目标数降到3个,不改变其他参数,运行时间可以压缩到35秒。这说明目标数增加带来的不只是质量解释难度的增加,更是全局排序计算量的压力。

评估时一定要跑至少10次独立的随机种子实验,取平均值和标准差。只跑一次就下结论的做法,在粒子群这类随机优化算法里完全没意义。

注意:Matlab里如果用contains、ismember这类函数嵌套在目标函数内层循环,千万要小心预分配数组。能向量化的地方一定要向量化。比如f1路径长度完全可以对所有粒子一次性算,不要写成逐粒子循环里的while。

5. 实际工程中那些坑:我踩过的和帮你避免的

5.1 目标函数数值量级不一致导致的效果漂移

最常见的坑:路径长度是以“千米”为单位,安全裕度是“米”,转弯代价是“度”。这几个值如果直接放在同一个Pareto支配关系里比较,支配关系本身不受影响,因为支配只关心每个目标的相对好坏。真正受影响的是拥挤距离计算——欧氏距离在量级差异大的目标上会被某一维完全主导,其他维度的差异淹没在数值噪声里,导致前沿分布严重偏向某一目标。

我踩过一次非常深刻的坑:路径长度大约500米,安全裕度大约5到15米,高度变化约30到60米,三者数值差一个数量级。跑出来的Pareto前沿在长度维度几乎没有区分,看起来所有解长度都很短,而其他维度被拉得很散。这不是算法错了,而是距离度量失效。

解决方法是目标归一化。每一代拿到所有粒子的目标后,对每个目标做min-max归一化,再用归一化后的值计算拥挤距离和支配强度。归一化之后的Pareto前沿分布明显均匀了很多。

5.2 无人机动力学约束跟路径点矛盾时的处理

这里要特别提醒:三维路径规划里优化的导航变量指的是航路点,但它不等于无人机能直接飞出来的轨迹。无人机有最小转弯半径、最大爬升角和速度限制。即使路径点本身不碰障碍物,连接路径点之间的平滑线段也可能因为转弯半径不够而撞墙。

我在第一次做较复杂地形时就是这样,算法跑出来的Pareto前沿异常漂亮,每条路径都绕开了建筑中心,但把三段路径衔接处的轨迹画出来,转弯半径远小于无人机实际参数。解决办法是在目标函数里额外补了一个“动态可行性代价”:对每个转折点计算实际转弯半径,如果小于最小允许转弯半径,就返回一个大罚值。

实际中,转弯代价turningCost可以和高度变化代价组合成一个目标,否则目标数量越加越多,问题变得不必要地难处理。做项目要有“目标聚合”的自觉。

5.3 NMOPSO与MOPSO/NSGA-II对比的测试结果

同场景下我做过三组对比:经典MOPSO(有网格法维护外部存档)、NMOPSO(非支配排序+拥挤距离修剪)和NSGA-II。

NSGA-II在低目标数时表现很稳定,3目标下Pareto前沿均匀度和收敛性都很好。但目标升到6个之后,NSGA-II的交叉变异算子在高维决策空间里的搜索效率明显下降,需要很多代才能追平NMOPSO。

经典MOPSO用网格法时,前期收敛很快,但后期多样性不足,多次运行的前沿解之间重叠率很高,容易被局部分布困住。

NMOPSO的优势主要体现在中等高维目标上,6目标24维决策变量的配置下,单调了近20%,花费时间增加的幅度还能接受。如果目标数量真的上升到10个以上,建议考虑加入目标降维或偏好诱导机制,否则任何MOPSO都很难保证前沿质量。

5.4 一个容易被忽略的Matlab性能细节

更新外部存档时,很多人的代码用判断add、save来暴力分配数组,用顶点3次就慢得可怕。要改成提前分配max_archive容量的矩阵,只用索引记录当前有效个数,修剪时用逻辑掩码一次性替换。不要在有隐藏快速循环的地方使用size反复计算,尤其不要把历史的最优解全量复制后重排。

还有就是在每次迭代里做约束修复时,如果遇到不可行解,不要完全丢弃,直接重置到某类边界值会让种群丧失引导信息。我倾向于做局部回弹:只把越界的分量拉回边界,并给该粒子的速度乘一个0.5的衰减系数,这样它下一轮会朝可行方向飞回,而不是原地停滞。

6. 一些更贴近实际任务的扩展经验

当前实现的目标函数是通用的,但真要在城市场景落地,导航变量肯定不止这些。复杂项目里,我会把“时间窗口”加进去,比如无人机必须在指定时间段完成侦察或配送任务,这时目标函数要增加时间窗违反代价。还有通信遮挡问题,城市楼群会引起信号遮挡,可以把“链路维持概率”作为目标加进去。

这些扩展的难点反而在目标定义上,不在NMOPSO本身。如果能把工程目标精确量化,算法框架不需要大改。

还有一个少有人提的实用技巧:不要让无人机贴着障碍物顶部飞行。“安全裕度”目标用的是路径点与障碍物边界的最短距离,但如果路径点离楼顶太近,又一整条路径绕了很远,距离目标会直接判定安全,但实际雷达回波和尾流影响会很大。我给安全裕度目标加了楼顶上空域的惩罚锥:超过建筑物高度一定范围,仍视为风险高度。这样优化出来的路径不会紧贴楼顶割草,更符合实操。

如果你也想复现这套代码,建议先做一个小规模场景,把三目标版本打通,再增加到六目标。直接上高维版本容易把时间都耗在Debug收敛问题上。我当初第一次跑完整代码,三个小时没跑完,最后定位到是障碍物检测函数里的for循环嵌套过深。把检测逻辑向量化后,运行时间降了三分之二。

NMOPSO这套框架在无人机三维路径规划场景下给我的整体印象是:它不像教课书里那么简单,也不像某些论文写的那样玄乎。核心还是环境建模精度和参数调优,算法本身只是把Pareto寻优做得足够稳固。跑通一批结果之后,你会发觉很多所谓“最新的高维多目标优化”,拆开都是经典机制的合理组合——包括非支配排序、拥挤距离、修复策略。真正考验人的,是你能不能把这些机制放到Matlab的数值通道里顺畅跑起来。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询