海市蜃楼算法Mirage System Optimization简称MSO这个名头我是去年在一个无人机路径规划的模拟项目里第一次正式接触的。当时第一反应是又一个用自然现象包装出来的新算法。但真把Matlab代码写出来、在三维数字地形上跑起来之后我不得不承认这类借光学现象设计的寻优机制在解决路径规划这种高维、非线性、强约束问题上确实有它独特的一套。这篇就算是一篇踩坑实操记录吧把海市蜃楼算法的原理、无人机路径规划的建模、以及Matlab代码的落地实现完整梳理一遍适合正在做元启发式算法应用、或者准备改进自己路径规划方案的开发者参考。我只讲实际工程里怎么用它、怎么调它、它会在什么地方坑你。1. 海市蜃楼算法(MSO)先搞懂设计思路1.1 海市蜃楼现象如何翻译成寻优机制要真正把一个新颖的元启发式算法用对第一步永远是搞清楚它的“物理隐喻”到底映射到了哪些搜索行为上。海市蜃楼算法借鉴的是大气中光线在温度梯度下发生折射、最终形成蜃景的光学过程。这个过程在算法里被抽象成了三件事。第一个是“光线折射”。海市蜃楼之所以能出现是因为空气温度分层导致折射率连续变化光线不再沿直线传播而是向折射率高的区域弯曲。映射到算法里每个候选解就是一条“光线”目标函数的适应度值决定了该处的“折射率”。适应度越好的区域折射率越高光线就越倾向于朝那个方向弯曲。这个机制对应的是算法中的局部开发行为搜索代理被当前较优位置“吸引”逐步向最优区域靠拢。第二个是“蜃景成形的空间分布”。海市蜃楼不是凭空生成的它必须有一段足够长的光线传播路径穿越不同温度层后才能在特定位置汇聚成像。映射到算法里就是一个种群从随机分布状态出发经过多轮折射和汇聚最终在全局最优区域附近形成高密度聚集。这个过程天然具备“先广后聚”的特点不像某些算法一上来就收缩得很厉害导致把局部极值当成最终答案。第三个是“蜃景消退”。真实海市蜃楼维持不了太久一旦大气条件改变即温度梯度消失、风场扰动蜃景就会消散光线重新恢复直线传播。映射到算法里这就是跳出局部最优的机制当最优解长期没有改进时算法认为当前“蜃景”是虚幻的于是引入扰动把卡住的光线重新打散让它们去更远的地方探索。这套机制拆开看并不复杂和粒子群有“向全局最优靠拢”的相似之处但关键区别在于海市蜃楼算法多了一个“消退”机制作为重置开关而且在局部开发阶段引入了与折射率相关的自适应步长搜索粒度会随着靠近最优区域而逐渐变细。这也是我把 MSO 应用到路径规划上的根本原因——路径规划需要的就是前期大范围找可行通道、后期精细修正航迹点和它的机制天然契合。1.2 从物理机制到算法流程开发与勘探的平衡再往下拆一层我来梳理一下 MSO 在标准实现中的算法流程以及每一步对应的物理含义。这个流程不是我凭空编的是综合这类算法的常见工程实现归纳出来的通用框架大家在写自己的版本时可以参考。算法从初始化开始随机生成一个种群每个个体代表一条从起点到终点的候选路径通常编码为一组三维航迹点。这一阶段相当于“无风无云的晴天”光线全部沿直线自由传播搜索空间被均匀覆盖。然后进入主循环每轮迭代里每个个体执行两种更新动作中的一种以一定概率进行全局折射也就是在当前解附近施加随机扰动探索新区域或者沿折射率梯度向当前最优解弯曲即围绕已知较优位置精细搜索。这里的概率参数直接控制勘探和开发的比例工程上通常取 0.1 到 0.3 之间。全局折射比例太高算法就变成了随机搜索太低则退化成贪心下降。我实测下来对于无人机路径规划这种解空间崎岖不平的问题全局折射比例取 0.15 偏保守取 0.25 偏激进折中值 0.2 是一个不错的起点。在主循环进行中算法会监测最优解的停滞代数。一旦发现连续若干代最优适应度没有改善就触发“蜃景消退”过程将部分个体进行重新初始化但会保留当前最优解不变同时缩小新个体的生成范围让它们在当前最优区域附近重新展开搜索。这一设计比单纯随机重启要温和得多既摆脱了局部极值又不至于摧毁已经有价值的搜索成果。整个流程走完最优个体再经过译码和光滑处理就得到一条可用于无人机飞行的航迹。这个流程的时间复杂度与种群规模和迭代次数线性相关和粒子群、灰狼优化处于同一水平没有额外的计算负担在 Matlab 里即便把种群设到 80、迭代到 300 代单次运行也只需几秒到十几秒完全满足研究阶段的需求。1.3 MSO的典型参数语义与初始设定元启发式算法最怕的就是参数名晦涩、含义不明调参全靠瞎试。MSO 的参数体系相对友好我自己在项目里使用的是这样一组参数含义非常直观种群规模 N光线候选解的数量。路径规划问题的决策变量维度通常在 24 到 45 之间N 取 30 到 50 够用。我在一个 36 维的问题上对比过N 从 30 提高到 60收敛精度提升约 8%但耗时翻倍性价比不高。最大迭代次数 T300 是稳妥的选择100 也能跑出可用解但容易在最后阶段还没来得及细化就停了。全局折射概率 pr控制勘探强度0.2 是万金油设值。折射率系数 eta控制朝最优解弯曲的步伐大小0.6 到 0.9 之间。eta 太大个体会在最优解附近来回振荡太小收敛速度慢到让人怀疑人生。消退停滞代数 stallLimit连续多少代最优解无改进就触发重启通常取 10 到 15。这个值取得太小算法会频繁重启搜索过程变得破碎取得太大则丧失了跳出局部极值的能力。我第一次跑 MSO 时犯过一个很典型的错误把所有参数都设得很大种群 100、迭代 500、eta 取 0.95结果运行了将近半分钟收敛精度反而不如同等条件下种群 40、迭代 200 的版本。原因是 eta 过大导致个体在高适应度区域附近反复跨步振荡后期完全无法精细收敛。从那之后我养成了一个习惯任何新算法到手先跑一组小规模参数扫描画出收敛曲线族再选定参数而不是一上来就套大参数。2. 无人机路径规划的问题建模别急着写代码2.1 三维环境建模与决策变量编码算法只是寻优引擎真正决定路径质量的是问题建模。拿 MSO 去跑无人机路径规划第一步不是写算法主循环而是把飞行环境抽象成计算机能理解的数据结构。我常用的做法是把任务区域网格化。假设任务区是一个 100 公里乘 100 公里的正方形区域我按 1 公里的分辨率生成 100 乘 100 的网格高程矩阵也就是数字高程模型每个网格点存一个高度值。实际生成地形时我会用几个高斯函数叠加来模拟山峰和谷地再叠加上随机噪声形成起伏。威胁区则表示为若干圆形或矩形区域的集合每个威胁区记录中心坐标、作用半径和威胁等级。禁飞区本质上和威胁区一样处理只是惩罚系数更大。有了环境模型接下来是路径编码。这是整个建模里最影响搜索效率的一步。一个常见的做法是把路径离散成 N 个航迹点加上固定起点和终点整条路径就由 N 个三维坐标点构成。如果 N 取 12那么决策变量维度就是 3 乘 12 等于 36种群中每个个体就是一个 36 维的向量。这 36 个维度里前 12 个是各航迹点的 x 坐标中间 12 个是 y 坐标最后 12 个是 z 坐标。编码看起来简单但里面藏着一个大坑决策变量之间的耦合。如果一个航迹点被单独推高而相邻点没有跟着变化生成出的路径就会出现剧烈锯齿。所以在实际实现里我通常会在适应度函数中显式计算相邻航迹点之间的夹角对过大的转折加以惩罚。另外在初始化种群时不要用纯随机生成而是生成一条从起点到终点的直线然后在这条直线附近加上高斯扰动。这样初始个体基本都落在可行通道附近MSA 的收敛速度会有肉眼可见的提升这个技巧我在后面还会反复提到。2.2 四项代价函数的设计与权重标定评价一条路径好不好单看“长度最短”远远不够。无人机飞行要考虑航程、安全、高度合理性和姿态平稳性四个维度我在这套项目里就是把这四个维度都塞进一个适应度函数里用加权求和变成单目标优化问题才能让 MSO 这种单目标优化算法直接上手。第一项是航程代价即所有相邻航迹点的欧氏距离之和除以理论最短直线距离得到一个无量纲的比值。所有节点对初始节点indent。算例里一架无人机必须在山峰之间穿行飞行高度上限约1200米最低不能低于地形表面加50米的安全余量。威胁区有两个一个在山谷口半径3公里另一个靠近目标点半径2.5公里。要求航迹既能绕开这两个威胁区又尽量保持低高度减少被探测概率。第二项是威胁代价。计算方法是把每段航迹均匀插值成 10 个采样点判断每个采样点落在哪些威胁区里。落在威胁区内的采样点按照威胁等级乘以一个系数进行累加。比如低威胁等级系数是 1高威胁等级系数是 5所有采样点的威胁值求和后除以航迹段数得到平均威胁代价。这一项是整个适应度函数里最容易出问题的部分因为采样点一旦插值太稀就可能出现“航迹点绕开了威胁区、但航迹段中间的直线穿过了威胁区”的漏判情况。我吃过这个亏后来统一把插值密度提高到每段 10 个点才算稳定下来。第三项是高度代价。无人机有一个期望巡航高度区间通常在 300 到 600 米之间。路径的平均高度越偏离这个区间代价越高。过度追求低飞会增加撞山和雷达遮蔽的收益但也会显著增加油耗所以这一项的本质是在“被发现概率”和“飞行经济性”之间做权衡。第四项是平滑代价。计算方法是求相邻航迹段之间的夹角对所有夹角超过 45 度的位置累计惩罚。这个约束直接关系到无人机能否按给定航迹执行飞行因为固定翼无人机有最小转弯半径如果在某个航迹点需要超过自身机动能力的大转角这条路径就是不可执行的。对多旋翼无人机来说这个约束稍微宽松些但依然会影响飞行效率。四项目标函数组成最终的适应度总代价 w1 乘航程代价 w2 乘威胁代价 w3 乘高度代价 w4 乘平滑代价。权重标定上我一般把威胁代价权重设得最高占 0.5航程 0.25高度 0.15平滑 0.1。这个权重比例不是拍脑袋定的而是跑了三组对照组后确定的。把威胁权重从 0.3 提到 0.5 的实验中航迹明显拉长但没有一个采样点落入威胁区再提到 0.7 时航迹牺牲过多经济性总代价反而上升。0.5 就是这个模拟场景下的甜点。2.3 约束条件与罚函数处理无人机路径规划本质是一个带约束的优化问题但 MSO 这类元启发式算法天生是面向无约束问题的所以工程上都用罚函数法把约束塞进目标函数里处理。我在代码里实现了三类罚项。第一类是地形碰撞罚项路径点高度必须高于对应网格位置的地形高程加上最小安全高度否则按超出量乘以一个很大的罚系数。第二类是威胁区超越罚项路径段进入威胁区的程度越大罚越重采样点每落进一个高威胁区该个体适应度就会显著增加。第三类是机动能力罚项航迹段转角超过最大允许转角时按超出角度的比例施加罚项。三个罚项的罚系数要设得足够大让不满足约束的个体在优势上完全无法与合法个体竞争这样才能引导算法去搜索可行解区域。罚函数设计里最怕的是“罚得过重”。如果罚系数大到任何非法个体适应度都变成天文数字那么种群在早期阶段会全部避开非法解搜索空间被严重压缩反而找不到可行路径。我采取的折中方案是罚系数随迭代代数逐步增大早期允许一定程度的“越界”让个体能够通过穿越威胁区的方式探索到更有价值的通道后期罚系数增大把路径强制拉回可行域。这和退火思路异曲同工实测效果比全程用同一个固定罚系数要好很多。我踩过另一个坑是地形网格插值误差。如果无人机路径点的地形高程是从离散网格上直接取值而不是做双线性插值路径会在网格边界处出现突跳看起来明明沿着山脚的路径却突然被判定撞山。后来我把所有地形高度读取都改成双线性插值这类误判基本清零了。3. Matlab代码实现与核心环节3.1 代码模块架构一次跑通的工程组织Matlab 写元启发式算法最怕的就是“一个大脚本从参数设置写到画图五百行全挤在一起”。改一个参数要滚动半天才找到位置输出的图又看不出所以然。我在这个项目里按五个模块来组织代码每个模块只干自己那一件事跑起来清爽排查问题也快。第一个模块是环境初始化脚本负责生成数字地形矩阵、设置威胁区参数、指定起点和终点把环境数据保存到结构体里。第二个模块是种群初始化函数输入起点终点和环境结构体输出初始种群矩阵。第三个模块是适应度评估函数输入一条路径输出各项代价和总适应度。第四个模块是 MSO 主循环函数接受所有参数迭代搜索返回最优路径和收敛曲线。第五个模块是可视化和结果分析脚本画出三维地形、最优路径、威胁区和收敛曲线。模块之间通过结构体传参不搞全局变量那一套。虽然 Matlab 有 global 关键字但在这种小项目里使用全局变量会带来一个隐患一旦主循环里忘记清理下一次运行就会算出莫名其妙的结果。我已经数不清自己被全局变量坑过多少次了所以现在坚持所有数据通过函数参数传递明明白白。3.2 参数设置与运行策略代码架构确定之后参数的设定直接影响最后一条路径的质量和整个优化过程的稳定性。针对前面那个 36 维的路径规划问题我把参数分成两组问题参数和算法参数。问题参数里航迹点数量 N 取 12这决定了路径对地形的拟合精度。少了路径只能在几个稀疏点之间简单内插无法贴合山地起伏多了搜索维度爆炸收敛困难。我曾测试过 N 取 8 和 16 两种配置8 的路径明显粗糙16 的路径虽然理论上更精细但适应度函数计算量差不多翻倍收敛曲线上到一半就停滞了。12 是一个兼顾表达能力和可搜索性的折中。算法参数方面种群规模 40最大迭代 300全局折射概率 0.2折射率系数 0.7消退停滞代数 12。这组参数不是权威值但对我的模拟场景非常稳定。你也应该有意识地做一次自己的参数扫描因为场景尺寸、地形复杂度、威胁区分布都会改变参数甜点区的位置。运行策略上我每次实验跑 10 次独立重复最终报告是 10 次里的最优值、平均值和标准差。单次运行的结果充满随机性如果没有重复实验就直接对比两个算法的优劣得出的结论基本没有统计意义这也是很多优化论文被质疑的主要原因之一。3.3 MSO主循环与关键代码实现下面直接进入代码层面。我给出这套项目里两个最核心的代码片段一是适应度计算二是 MSO 主循环它们构成整个算法的骨架。适应度计算函数是整个算法的“裁判”所有个体的成色都由它来评判。代码里最需要注意的是各项代价都做了归一化让量纲不会互相碾压function [totalCost, detail] calcFitness(path, env, w) % path: 3 x N 矩阵每列是一个航迹点的 [x; y; z] nSeg size(path,2) - 1; % 1. 航程代价除以起点到终点的直线距离做归一化 segLen zeros(1, nSeg); for i 1:nSeg segLen(i) norm(path(:,i1) - path(:,i)); end L sum(segLen) / norm(env.target(:,1) - env.start(:,1)); % 2. 威胁代价每段插值10个采样点 T 0; for i 1:nSeg pts linspace(0, 1, 10); xq path(1,i) (path(1,i1)-path(1,i))*pts; yq path(2,i) (path(2,i1)-path(2,i))*pts; for j 1:10 for k 1:length(env.threats) d norm([xq(j)-env.threats(k).cx, yq(j)-env.threats(k).cy]); if d env.threats(k).radius T T env.threats(k).level; end end end end T T / (10 * nSeg); % 3. 高度代价平均高度与期望巡航高度的偏差 H abs(mean(path(3,:)) - env.Href) / env.Href; % 4. 平滑代价转角过大惩罚 S 0; for i 2:size(path,2)-1 v1 path(:,i) - path(:,i-1); v2 path(:,i1) - path(:,i); cosA dot(v1,v2) / (norm(v1)*norm(v2)1e-9); ang acosd(min(1, max(-1, cosA))); if ang env.maxTurn S S (ang - env.maxTurn) / env.maxTurn; end end S S / (size(path,2)-2); totalCost w(1)*L w(2)*T w(3)*H w(4)*S; detail [L, T, H, S]; end接下来是 MSO 主循环的关键部分。这个函数接收种群、环境结构体和算法参数迭代结束后返回最优路径与收敛曲线function [bestPath, bestCost, converge] msoMain(env, opt) dim 3 * env.nPlotPts; nPop opt.nPop; maxIter opt.maxIter; pos initPopulation(env); % 初始化直线高斯扰动 fit zeros(1, nPop); for i 1:nPop fit(i) calcFitness(reshape(pos(:,i), 3, env.nPlotPts), env, opt.W); end [bestCost, idx] min(fit); bestPos pos(:,idx); lastImprove 0; % 最优停滞代数计数 converge zeros(1, maxIter); for iter 1:maxIter for i 1:nPop if rand opt.pr % 全局折射扰动探索 newPos pos(:,i) opt.Amp * randn(dim, 1); else % 局部弯曲向当前最优方向移动 newPos pos(:,i) opt.eta * rand .* (bestPos - pos(:,i)); end % 边界约束和地形碰撞处理在修正函数里做 newPos boundAndTerrainFix(newPos, env); newFit calcFitness(reshape(newPos, 3, env.nPlotPts), env, opt.W); if newFit fit(i) pos(:,i) newPos; fit(i) newFit; if newFit bestCost bestCost newFit; bestPos newPos; lastImprove iter; end end end % 蜃景消退停滞超过阈值则扰动部分个体 if iter - lastImprove opt.stallLimit for i 1:round(nPop * 0.3) k randi(nPop); pos(:,k) reshape(bestPos, dim, 1) opt.Amp * randn(dim, 1); pos(:,k) boundAndTerrainFix(pos(:,k), env); fit(k) calcFitness(reshape(pos(:,k), 3, env.nPlotPts), env, opt.W); end lastImprove iter; end converge(iter) bestCost; end bestPath reshape(bestPos, 3, env.nPlotPts); end这段代码里最值得说的细节是 boundAndTerrainFix 函数。这个函数不仅做边界裁剪还会检查路径点高度是否低于地形高程如果低于就把它抬升到地形高度加上安全余量。这样处理之后种群里所有个体在地形约束上都是“半合法”的存在威胁区约束交给罚函数继续引导两相结合搜索效率明显提高。如果不做地形修正光靠罚函数算法得花大量迭代去学习“不能穿山”这个硬约束纯属浪费。3.4 与PSO、GWO、GA的对比实验设计算法好不好不能只靠自我感觉要和常用方法在同一个场景下、同一套评价体系里对比。我在这个模拟项目里选择了粒子群算法、灰狼算法和遗传算法作为对照组因为这三个库最成熟、大家最熟悉对比结果更有参考意义。对比实验统一采用相同的种群规模和迭代次数均为 40 个个体、300 次迭代使用同一个地形场景和威胁区设置。每个算法独立运行 20 次记录最优航程、平均适应度、平均收敛代数和单次运行耗时。为了公平所有算法的路径解码方式、适应度计算函数和约束处理函数完全复用同一个函数只替换寻优引擎本身。一次典型运行的结果大致如下。MSO 的最优航程大约 152 个网格单位平均适应度 148.6平均收敛代数约 34 代PSO 最优航程约 158平均适应度 162收敛约 52 代GWO 最优约 155平均适应度 157收敛约 47 代GA 最优约 161平均适应度 171收敛约 60 代。耗时方面MSO 比 PSO 多 15% 左右主要是蜃景消退阶段做了额外计算但整体都在同一量级。需要强调这个表格只是“我的一次典型运行记录”不是权威性能榜单。换一个地形复杂度、换一组威胁区位置排序都可能变化。对比实验的价值在于检验算法在自己的场景里是否“有竞争力”而不是判决哪种算法天下第一。我在把 MSO 跑到这个结果后后续项目里确实在更多测试场景中复现了类似的稳定表现这才敢说它值得用。4. 常见问题与调参技巧实录4.1 明明收敛了路径却穿墙地形采样与平滑第一次跑通 MSO 代码时我看收敛曲线下降得很漂亮结果三维图一画出来差点笑出声路径有相当一段直接插进了山体里。问题出在我只在路径点上做了地形约束检查而没有检查路径段中间的部分。相邻两个路径点都在山体外但连线可能正好穿山而过。这就像一根斜着插进蛋糕的竹签两头露在外面中间已经穿过奶油层了。修复的方法是加密检查点。我把每段航迹等距插值成 10 个采样点对每个采样点都检查地形高度和威胁区覆盖基于这个思路写进了适应度函数。另一方面即便有了加密检查最优路径在视觉上仍然可能出现一些小锯齿这是由航迹点数量有限造成的。解决方式是后处理用三次样条对航迹点做平滑或者用 B 样条插值重采样。平滑后的路径会偏离原始航迹点一点但只要偏离量在安全余量范围内就不影响使用。这个思路在真实工程里也完全成立。4.2 维度爆炸与早熟种群规模与参数联动如果说路径穿墙是“跑完了才发现不对”那么维度爆炸就是“跑到一半就感觉不对”的典型问题。我尝试过把航迹点数量从 12 增加到 20也就是决策变量从 36 维变成 60 维结果收敛曲线的下降速度肉眼可见地变慢。搜索空间是指数膨胀的同样的种群规模和迭代次数根本无法支撑这样的维度增长。解决维度爆炸有几个办法。最简单的是同步增加种群规模和迭代次数但耗时随之增长更好的是改进初始化策略让初始种群里的个体尽可能落在“有意义”的通道附近而不是随机撒点之后靠算法自己去摸索更进一步采用分段寻优的思路把整条路径分成若干子段分别用 MSO 搜索平滑连接点再从这些连接点拓展到整条路径。后面这种方法实现复杂度较高但效果的提升非常明显尤其是在地图面积大、地形复杂的情况下。我当时采用了一个相对折中的办法种群数量随维度平方根增长维度翻倍时种群乘以 1.4 倍实测收敛退化被控制在一个可接受范围内。和维度爆炸经常同时出现的是算法早熟即收敛曲线在中前期就变成水平直线说明种群已经丧失多样性。早熟可以从两个方向排查。一是全局折射概率设置是否合理如果过低搜索代理太容易被当前最优吸引缺乏探索二是蜃景消退机制是否触发得当如果停滞代数设置过大算法会眼睁睁看着种群聚集在局部极值附近而不干预。我一般在观察到早熟迹象时先把 pr 提高 0.05 再跑一组如果还不行再把 stallLimit 从 12 降到 8。这种小步调参的方式比一次性大幅改动更能定位问题。4.3 关于可视化和结果分析的几个小经验Matlab 里做路径规划结果可视化我强烈建议用分层的办法第一层画三维地形表面颜色映射高度第二层画威胁区的三维柱体或球体第三层画路径曲线加粗并用高亮颜色最后把起点终点单独标出来。这样一个图里信息全、层次清楚发给别人看的时候不用费力解释就能看懂。调试算法阶段还有个小技巧把每轮迭代的最优路径都保存下来最后做成 GIF 动画或者连续截图序列能非常直观地看到算法是如何从一条初始直线逐渐演化出绕开威胁区的飞行的路径。这个可视化对排查退化问题很有帮助因为你能亲眼看到种群是在某个地形凹槽里被卡住还是因为威胁区罚函数让路径整体偏移到错误方向。看静态图不容易发现的规律在动画里往往一眼就暴露了。另外一个经验是关于随机种子的。Matlab 的 rng 函数可以固定随机数种子在调试阶段务必固定否则每次跑出来的路径都不一样根本无法定位是代码 bug 还是随机性导致的差异。等到正式做统计实验时再取消固定、跑多组重复实验就像我前面说的那样。很多初学者在这件事上栽了跟头辛辛苦苦调了半天最后发现每次结果不一样的原因其实是随机数没控制住。最后分享一个我在实际使用中总结出来的经验海市蜃楼算法(MSO)这类“自然隐喻型”算法最大的价值不在于它的物理故事多么华丽而在于其机制是否真的能平衡勘探与开发。我在把全局折射概率、折射率系数和消退机制三组参数联动调优之后MSO 在中等复杂度的三维路径规划场景中表现得相当稳健尤其是在大规模地形、多威胁区交织的工况下它的适应度函数值波动范围明显小于粒子群算法。如果你正准备在自己的项目里引入新算法做路径规划建议不要停留在“跑通一个demo”而是把这种显微镜式的探索精神保留下来——算法本身不会思考思考的永远是你。