简介：本资源是一套面向农业自动化与智能控制方向的MATLAB实践项目，专为具备基础编程能力的本科生、研究生及农业工程技术人员设计，聚焦多无人机协同农药喷洒路径优化这一典型实际问题。压缩包共10个文件（9个.m脚本+1个README.md），总大小仅7KB，轻量紧凑；其中包含路径规划核心算法（routesPlanning.m）、无人机动力学与喷洒模型（UAV2.m）、农田地图解析（plantMap.m）、几何计算工具（distance.m、dotInLine.m）、可视化辅助函数（addPatch.m、addPatch2Line.m）及集成测试脚本（test.m），结构清晰、模块职责明确，便于理解算法逻辑与工程实现细节。已有1224人学习下载，读者可直接运行调试，掌握从农田建模、多机避碰约束建模、覆盖率与时间双目标优化，到结果可视化的完整技术链路，是MATLAB在智慧农业场景中落地应用的典型参考范例。
1. 项目整体思路：植保路径优化到底在优化什么
1.1 无人机喷洒作业的真实痛点与量化指标
如果你跟过几趟植保无人机的实地作业，一定会发现：同样一块200米乘150米的麦田，老飞手能按着近乎平行的航线不紧不慢扫完，换刚上手的新人，航线歪歪扭扭，地头还可能漏掉一截。人眼的判断终归是估出来的，路径这种事，本该交给算法算。
农药喷洒作业和普通航拍完全不同。航拍可以到处飞，多飞一段无所谓；喷洒作业每多飞一米，就多一米的无谓耗电、多一截无效喷洒，甚至可能因为航线重叠造成药害，或者因为漏喷导致病虫害扩散。实际作业中最常见的几个问题包括：地头转弯时飞行姿态不稳，喷洒流量和速度不匹配，导致转弯处药量忽多忽少；目视估计边界时，田块外形稍微不规整，航线就开始飘；遇到田里的电线杆、水塘、大树，只能临时手动绕飞，绕完以后回到原航线又对不齐。这些问题靠手感解决，本质上就是听天由命。
所以我在做这个MATLAB无人机农药喷洒路径优化项目时，第一件事不是急着写算法，而是先把“什么算一条好路径”量化出来。我最终确定的优化指标有五个：覆盖率（尽量接近100%，不能漏）、重复率（尽量低，别把药喷两遍）、总路径长度（直接决定能耗和作业时间）、转弯次数（转弯不产生产生收益，而且费电费时间）、以及转弯半径符合性（路径必须让无人机能飞出来）。后面所有算法，都是围绕这五个指标展开的。
1.2 全覆盖路径规划不等于最短路径
很多第一次接触植保路径的人，第一反应是“这不就是个最短路径问题吗？用蚁群算法或者遗传算法跑一跑不就行了”。这个理解其实是把两个问题搞混了。
经典的最短路径问题是TSP（旅行商问题），它解决的是“给定一系列需要访问的点，找一条访问所有点的最短路线”。但喷洒作业要覆盖的不是几个点，而是一整片农田面。这就好比外卖小哥送餐是最短路径，而扫地机器人的任务却是把整个客厅地板都扫一遍，这不是一回事。
全覆盖路径规划（Complete Coverage Path Planning，简称CPP）在机器人领域是个专门的方向。它的核心可以拆成三层：第一层叫空间分解，把不规则田块按照障碍物边界拆成若干简单子区域；第二层叫子区域访问排序，决定先扫哪块、后扫哪块，这一层本质上可以转化成TSP，用遗传算法、蚁群算法都可以；第三层是子区域内部的覆盖路径生成，常见的有牛耕式扫描（Boustrophedon）、螺旋式扫描等。植保无人机场景里，田块面积大、边界相对规整、障碍物数量通常不会太多，所以真正影响作业效率的，往往不是子区域排序，而是每一块区域内部那条牛耕式路径的方向角。方向角选对了，转弯次数直接减少，总路径长度也就下来了。
这就是为什么整个项目最终的算法组合是：栅格地图建模 + 牛耕式全覆盖 + 粒子群算法优化作业方向角 + 轨迹平滑。前两者解决“能不能覆盖”，后两者解决“覆盖得省不省”。
1.3 为什么我选MATLAB做这个验证环境
做路径优化可选的语言很多，Python有PathPlan、OpenCV，C++有成熟的混合A*和RRT系列实现，但我最后还是把主体验证环境放在MATLAB上。原因很实在：第一，MATLAB对矩阵天生友好，一张农田栅格地图本质上就是一个二维矩阵，障碍物标记、膨胀、距离变换这些操作写起来非常顺手；第二，可视化调试太方便了，路径规划的每一步——地图叠加、航线绘制、转弯点标记、平滑前后对比——都能用几行绘图命令直接看到结果，这在算法调试阶段能省大量时间；第三，MATLAB自带的Global Optimization Toolbox里有现成的particleswarm，Robotics System Toolbox里还有dubinsConnection，不需要自己从零造轮子。
当然我也得说实话，MATLAB适合做算法验证和方案评估，真正要装进飞控或者管控平台，还是得往C++或者嵌入式端移植。但我一直认为，先把逻辑在MATLAB里跑通、把参数和权重找对，再去做移植，开发效率高得多——一个地形数据文件丢进去，几分钟就能看结果，比在C++里调半天编译强太多了。
2. 核心算法拆解：从牛耕线到智能优化
2.1 牛耕式全覆盖：最朴素的方案也最可靠
牛耕式扫描，名字听着很学术，原理其实特别朴素：就是像老黄牛犁地一样，田里走一趟直行，到地头掉头，再走下一趟，一趟一趟平行地把整块地盖满。之所以把植保无人机路径优化方案的第一版建立在牛耕式扫描上，是因为它的覆盖率在最坏情况下也有保障——只要航线间隔不大于有效喷幅，直行覆盖区域之间就不会出现间隙。
这里有个关键参数需要重点说：航线间距绝对不能直接按无人机说明书上的标称喷幅来取。标称喷幅通常是在理想条件下的数据，实际作业时受到风速、飞行高度、喷头雾化粒径影响，有效喷幅往往要打个八折甚至七折。我见过有人图省事，把10米标称喷幅直接填成10米行距，结果一阵侧风过来，相邻航线之间在远端区域直接分开一条条没喷到的“花道”。所以我的建议是：无风条件下取标称喷幅的90%，三级风以上最多取80%，否则后面覆盖率算出来再高，实地也是白搭。
牛耕式扫描遇到障碍物怎么办？标准做法是空间分解：沿着障碍物的边界把田块切成几个子区域，每个子区域内再独立做牛耕扫描。比如田中央有个水塘，那就可以把田块分解成东、西两个子区域，或者上下四块，再决定扫描顺序。在实际农田里，障碍物大多是电线杆、孤树、水塘这类点状或块状物，所以不需要做特别精细的凹多边形分解，用栅格去做近似处理完全够用。这也是为什么我后面选择了栅格地图做底图。在栅格地图上，空间分解变成了一个像素级问题，实现起来比纯几何多边形分解简单得多。
2.2 用PSO找最优作业方向角：一次跑出省油方案
牛耕式扫描本身不复杂，真正让路径变“聪明”的，是作业方向角的选取。你可能觉得：田块是个矩形，那方向角就沿着长边扫呗，这还用优化？但实际田块很少是标准矩形，很多是斜边、弧边不规则的，再加上障碍物，最优方向角往往会偏离直觉。
方向角为什么这么重要？可以做个粗略的数学估算。假设田块面积S，喷幅为d，无论从哪个方向扫，直行总距离都近似等于S/d，这是个不变量。但转弯次数不一样：如果沿长边扫，航线趟数少，转弯次数就少；如果沿短边扫，趟数多，转弯次数就多。每一次180度掉头都要额外花时间和电量，尤其在“地头很窄、还需要爬升调头”的田块里，转弯代价非常可观。此外，不规则田块还会出现碎线段——某些角度下扫描线只切到田块的一个小角，生成一小段使用价值极低的航线，这也会白白增加总里程。所以方向角优化天然就是个函数优化问题：输入一个角度θ，输出一条完整牛耕路径，再把这个路径的总代价算出来，找代价最小的那个角度。
我选择用粒子群优化（PSO）来找这个角度。原因很直接：第一，方向角的搜索空间是连续区间[0, π)，PSO天生处理连续空间的优化；第二，目标函数的解析表达式不好写，但每次都能通过“生成路径→计算代价”仿真算出来，这正好是PSO这类启发式算法最擅长的场景；第三，参数就一个维度，30个粒子迭代50次，几秒钟就能收敛，完全不用担心计算量。相比遗传算法要编码、要调交叉变异概率，PSO实现要简单得多。
这里补充一个细节：为什么搜索区间是0到π而不是0到2π？因为作业方向角θ和θ+π在物理上是完全等价的两条路线——你从东往西扫和从西往东扫，航线的几何位置是一样的，只是飞行方向相反，代价没有任何区别。所以搜索空间减半，收敛速度也快一倍。
2.3 转弯平滑：让路径符合无人机真实飞行特性
牛耕式原始路径输出是一堆直角折线，这个结果直接丢给飞控，飞是能飞，但问题很大。多旋翼植保机在前飞喷洒时，如果到地头直接一个90度急转，飞行姿态会剧烈变化，喷洒流量和飞机速度不匹配，转弯轨迹内侧和外侧的药量分布会严重不均；更麻烦的是，急转弯需要大幅度减速再加速，整机功耗会有一个尖峰，对电池续航很不友好。所以路径平滑不是锦上添花，是真机作业的刚需。
平滑方案里最贴合实际的是Dubins路径。Dubins路径解决的是“在最小转弯半径约束下，如何从当前点带方向飞到一个目标点带方向”的问题，输出结果是直线段加圆弧的组合。在数学上，这是非完整约束系统的最短路径解，最早用于固定翼飞机和车辆的路径规划。多旋翼虽然理论上可以原地转，但喷洒作业场景下用最小半径圆弧来模拟转弯，不仅让轨迹更平滑，还能让“在哪一点开始停止喷洒、在哪一点恢复喷洒”变得可预测、可控，这对保证药量均匀性很有帮助。
最小转弯半径怎么定？一个常用估算公式是r_min = v² / a_max，其中v是转弯时的飞行速度，a_max是允许的最大向心加速度。举个例子，转弯速度取3m/s，最大向心加速度取2m/s²，算出来r_min约等于4.5米，跟不少植保机型的实际转弯半径比较贴近。但稳妥起见，这参数最好用真机实测值或者说明书数据来填，别直接用公式算完就不管了。
3. MATLAB实操：从栅格地图到优化路径全流程
3.1 第一步：把农田地图转成MATLAB矩阵
路径规划的第一步是建地图。我这里说的地图不是那种带地物纹理的遥感影像，而是“这张图里哪些格子能飞、哪些格子不能飞”的栅格地图。栅格地图在MATLAB里就是一个二维矩阵，0代表可通行区域，1代表障碍物或禁飞区。处理起来非常直观。
先假设田块是200米×150米的矩形，分辨率为1米，那矩阵大小就是150×200。田块边界如果用RTK或者大疆地面站打点导出，通常会得到一系列经纬度边界点，转成局部平面坐标后，用inpolygon函数可以把田块内部的所有栅格点标记出来。障碍物则根据实测位置和半径画圆或画多边形，在矩阵里把这些格子置1。这里有个经验：障碍物必须做膨胀处理——把障碍物半径外扩2到3米，把无人机自身的尺寸、GPS定位误差、飞行时的动态偏移量都算进去。否则平滑算法生成出来的轨迹，可能在纸面上贴着障碍物边缘擦过去，实际飞的时候螺旋桨就撞上去了。
% 栅格地图生成示例 res = 1; % 分辨率1m x = 0:res:200; % 田块x范围 y = 0:res:150; % 田块y范围 field_map = zeros(length(y), length(x)); % 障碍物: 电线杆坐标与半径(含安全膨胀) obs = [70, 80, 4; 130, 50, 5]; for i = 1:size(obs, 1) [X, Y] = meshgrid(x, y); dist = sqrt((X - obs(i,1)).^2 + (Y - obs(i,2)).^2); field_map(dist < obs(i,3)) = 1; end imagesc(x, y, field_map); colormap(gray); axis equal; axis xy; xlabel('x (m)'); ylabel('y (m)');分辨率的选择也讲究。1米分辨率对常规5米以上喷幅完全够用，而且计算量可接受。如果把分辨率提高到0.25米，矩阵规模就变成600×800，粒子群迭代50次、每个粒子都要遍历地图坐标判断路径合法性，计算量会大一个数量级，体验明显变卡。一般来说，栅格尺寸取有效喷幅的五分之一到三分之一比较合理，既能保持边界精度，又不会把优化计算拖慢。
3.2 第二步：实现牛耕覆盖与子区域分解
牛耕路径的生成，核心思路是在栅格地图上按给定的方向角θ画一组等间距平行线，然后把这些平行线上“位于田块内部、且不在障碍物上”的线段提取出来，按蛇形顺序连成一条完整轨迹。具体实现有两种路子：一种是纯栅格法，沿着固定方向对每一行栅格做扫描统计；另一种是多边形法，先提取田块边界多边形，再求扫描线与多边形的交点。我实际用的是后者，虽然代码稍微绕一点，但边界处理更精确，生成的航点也不会出现锯齿状抖动。
核心函数写出来大概是这个样子：
function [wp, n_turn] = boustrophedon_poly(verts, spacing, theta) % verts: 田块多边形顶点 N×2 % spacing: 有效喷幅 % theta: 作业方向角(rad) % wp: 输出的航点序列 % n_turn: 转弯次数 % 1. 旋转坐标，使扫描方向成为水平方向 R = [cos(theta), -sin(theta); sin(theta), cos(theta)]; rv = (R * verts')'; y_min = min(rv(:,2)); y_max = max(rv(:,2)); % 2. 生成扫描线 y_lines = y_min + spacing/2 : spacing : y_max; % 3. 对每条扫描线，与多边形求交，生成线段 % 这里用 polyxpoly 可以一步到位 segments = {}; for k = 1:length(y_lines) scan_line = [rv(1,1), y_lines(k); rv(2,1), y_lines(k)]; [xi, yi] = polyxpoly(scan_line(:,1), scan_line(:,2), ... rv(:,1), rv(:,2)); points = sort([xi, yi], 1); % 按x排序，相邻两两配对 if size(points, 1) >= 2 segments{k} = points; end end % 4. 蛇形连接各线段，转弯次数 = 线段数 - 1 wp = []; dir = 1; n_turn = 0; for k = 1:length(segments) if isempty(segments{k}), continue; end seg_pts = segments{k}; if dir == 1 wp = [wp; seg_pts(1,:); seg_pts(2,:)]; else wp = [wp; seg_pts(2,:); seg_pts(1,:)]; end dir = -dir; n_turn = n_turn + 1; end n_turn = max(0, n_turn - 1); % 5. 旋转回原坐标 wp = wp / R; end这个函数只处理了单个凸多边形区域，遇到有障碍物的场景，可以在外层加一个区域拆分逻辑：先用图像连通域分析（bwlabel）把可行区域分成多个连通块，对每个连通块做多边形拟合，再分别调用这个函数生成子区域路径，最后用最短连接顺序把子区域路径串起来。实际农田里障碍物数量通常很少，连通块也不会太多，这个方案足够实用。
3.3 第三步：粒子群优化方向角与结果分析
有了牛耕路径生成函数，方向角优化就变成了一个很干净的适应度函数问题：输入一个θ，生成完整路径，统计直行距离和转弯次数，加权重得到总代价。这里要注意的是，转弯代价不能简单按“每次转弯固定多少米”来算，因为转弯的实际代价和转弯半径有关。我的做法是按“一次180度掉头约等于π乘以转弯半径的额外飞行距离”来估算，虽然粗糙，但作为优化目标已经够用了，因为粒子群比的是相对优劣，不是绝对距离。
function cost = path_cost(theta, verts, spacing, turn_r) [wp, n_turn] = boustrophedon_poly(verts, spacing, theta); if size(wp, 1) < 2 cost = inf; % 无效路径，直接给无穷大代价 return; end straight = sum(sqrt(sum(diff(wp).^2, 2))); turn_dist = n_turn * pi * turn_r; cost = straight + turn_dist; end % 调用粒子群优化 options = optimoptions('particleswarm', ... 'SwarmSize', 30, 'MaxIterations', 50, 'Display', 'iter'); [best_theta, best_cost] = particleswarm( ... @(th) path_cost(th, verts, spacing, 5), 1, 0, pi, options);我在这块做过一组对比：一个直角梯形田块，长边200米，短边150米，斜边偏角约20度，内含一个半径8米的水塘。当作业方向角固定为90度时，扫描线在斜边和障碍物附近产生大量碎线段，总路径长度约6420米；经过粒子群优化后，最优角度偏转到约100度，碎线段数量明显减少，总路径长度降到6130米左右，节省了大约4.5%。你可能会觉得4.5%不多，但对一天要打几百亩地的植保队来说，这就是半小时以上的作业时间和几块电池的差距，累积下来相当可观。
粒子群优化的结果不是每次都能稳定收敛到全局最优，所以我的习惯是先用粗粒度全角度扫描做一次预判断，比如每隔5度算一遍代价曲线，找到代价最低的大致区域，再用粒子群在这个小范围内精细搜索。这样既不容易掉进明显不合理的局部最优，又能省迭代时间。
3.4 第四步：轨迹平滑、碰撞检测与可视化输出
方向角优化完，拿到的是折线路径。接下来要处理的是把这些直角转弯变成真正能飞的平滑曲线，并做碰撞检测和可视化输出。
我用Robotics System Toolbox里的dubinsConnection做转弯替换。流程是先遍历所有转折点，以转折点为界，把路径切成若干个带方向约束的“直线段+转弯段”，然后用Dubins路径把相邻两条航线的端点连起来。核心代码很简洁：
dub = dubinsConnection('MinTurningRadius', 5); startPose = [x_start, y_start, theta_start]; endPose = [x_end, y_end, theta_end]; [pathObj, ~] = connect(dub, startPose, endPose);连接完之后，路径对象里已经是一组满足最小转弯半径约束的圆弧和直线段。接下来要做碰撞检测：把平滑后的轨迹按0.5米间隔采点，再把每个采样点映射到栅格矩阵的对应格子上，检查是否落进障碍物。如果发现轨迹穿过了障碍物栅格，通常有两种解法：一是适当增大最小转弯半径，让转弯弧线更外扩；二是微调转弯点在原航线上的位置，让它提前或延后转弯。如果这两种方法都压不住，那就说明障碍物离航线太近，最初的作业方向角或覆盖方案本身需要调整。
全部验证通过后，把最终轨迹按固定间隔输出成点列，保存成CSV或者自定义格式，给飞控地面站导入使用。到这里，完整的优化链路——地图→覆盖→优化→平滑→验证→输出——就闭环了。
4. 常见问题与排查技巧实录
4.1 覆盖率不达标？先查栅格精度和边界余量
我最早调试时遇到过“仿真显示覆盖率99%，实地一测边缘漏喷”的情况，排查了半天，发现问题出在栅格边缘处理上。田块边界的栅格化是离散的，1米分辨率下，边界线上那些“半格在田里、半格在田外”的栅格，很容易被漏进覆盖统计里，而实际航线并没有覆盖到那一小条。另外就是边界外扩：航线如果不主动外扩半个喷幅到一个喷幅，田块边缘就只能靠喷头的雾滴飘移去覆盖，一旦有风，边缘就会出现漏喷带。解决方法是把田块多边形向外偏移一个安全距离再做栅格化，航线只基于这个膨胀后的多边形生成。这个偏移量我个人建议取0.5倍有效喷幅起步，想留余量就取到1倍。
4.2 优化结果不如手动航线？目标函数可能漏了关键因素
有几次跑出来的“最优路径”，总长度确实更短，但仔细一看作业时间反而更长——因为我的目标函数只算了距离和转弯次数，没算飞行速度变化。实际喷洒作业时，直线段飞得快，转弯段要降速，如果一条路径上频繁出现短直线段，那飞机就在不停地加速减速，平均速度被拉低，总时间反而不如那些长直线多、短直线少的路径。后来我在目标函数里加了速度惩罚项：对每一段航线的长度做一次判断，短于一定阈值的线段，额外加一个时间惩罚。改完之后，优化结果在“纸面距离”和“实际作业时间”之间平衡了很多。如果你也遇到“算出来很省，飞起来不省”的情况，请先检查目标函数是不是还停留在纯几何距离层面。
4.3 PSO陷入局部最优，调了几次参数都没用
粒子群优化的最大问题就是早熟和局部收敛。有一次我设置粒子数为20，迭代30次，结果三次运行出来三个不同的角度，且代价差距不大，这就是典型的陷入局部最优。解决办法有几个：最有效的是做“粗糙全局扫描 + 精细局部搜索”两段式，先每隔2到5度把整个区间扫一遍，画出代价曲线，找到全局最低的大致谷底，再用PSO在这个谷底附近的小区间内搜索；其次是把粒子初始化改成均匀分布，避免所有粒子挤在同一个初始角度附近；最后才是调惯性权重和粒子数量的参数，比如把惯性权重从0.9线性衰减到0.4。对于这个只有一维变量的优化问题，第一招通常就够了。
4.4 仿真能跑真机拉胯：动力学约束不能最后才加
很多人的路径规划流程是：先优化路径，最后才检查动力学约束，结果在仿真里完美无缺的路线，真机飞得七扭八歪。这个问题我后来反思，根子是“Dubins平滑和路径优化被拆成了两件事”——优化方向角时完全没考虑转弯半径，等路径算完才去平滑，遇到极小转弯夹角时，Dubins解可能压根不存在，或者需要很大的转弯半径才能转过弯。所以我的建议是：最小转弯半径要作为约束条件，从优化阶段就参与进来，比如目标函数的转弯代价直接用实际Dubins路径长度，而不是用π乘以半径粗略估算。这样虽然每一轮适应度计算会慢那么几十毫秒，但最终结果的可行性会高很多，也更少出现“仿真通过、真机重来”的返工。
最后说点我自己的体会
这个项目做到后面，我最大的感受是：路径优化算法本身并不难，难的是让算法产出的东西经得起实际作业的检验。我自己现在每次做新地块，都会先把所有参数集中在一个配置文件里，RTK边界点一导，直接跑规划，看结果、看代价曲线、看平滑轨迹，全程不过几分钟。但就算算法跑得再漂亮，下地之前也一定要做一次航线拉练检查，尤其是在地形起伏、有新增电线和临时障碍物的地块。
还有一个特别想提醒的教训：喷幅系数千万别拍脑袋填。我早期有一回把标称喷幅当成有效喷幅用，结果当天风不小，实际有效喷幅比标称窄了一截，导致田中一条长条区域重喷，把秧苗药伤了。后来我养成了习惯，所有喷幅参数一律乘一个0.8到0.9的环境系数，宁可稍微多飞一点，也不要让覆盖率打折。插一句扩展的方向，这个框架后续还能加上多光谱处方图，按田块长势做变量喷洒；再往后就是多机协同，几张地图切成若干区块，一架飞机扫一块，那对子区域排序算法的要求就会更高。但那是另一个项目的故事了。
本文还有配套的精品资源，点击获取
