分水岭算法做果树图像分割与无人机路径生成完整实现

简介一份基于MATLAB的果树图像分割与无人机路径自动生成资源面向图像处理与无人机应用方向的开发者和学生。压缩包共八个文件包含六个m脚本和两张jpg果树图像整体仅272KB轻量易部署。m脚本搭建了完整处理链先通过高斯滤波去噪、直方图均衡化增强、Canny边缘检测等预处理突出果实边缘再设置种子标记、模拟水位上升完成分水岭变换将果树与背景分离随后识别目标位置借助A星或RRT等搜索算法规划避障路径并拆分航点、平滑轨迹以生成可执行航线配套图像数据可直接复现全流程。已有119人学习。代码组织清晰主流程与搜索函数分离便于拆解每个环节的参数结合样例图像能直观理解分水岭算法在农业场景中的运用以及无人机路径规划的设计思路是一份教学与实操兼顾的紧凑资源。1. 分水岭算法做果树图像分割生成无人机路径从树冠掩码到可飞航点的完整链路分水岭算法做果树图像分割再把分割出的树冠掩码接上无人机路径生成这套流程在植保喷药、变量施肥和产量估测场景里都能直接落地。相比U-Net这类深度图像分割模型分水岭不需要逐像素标注数据一张正射图加几个可调参数就能跑通代价是阈值和核参数要跟着树冠大小、光照场景微调。我去年在果园项目里用这个思路替代人工圈选树冠整片园子的路径生成从半天缩到十几分钟。这篇文章按「分割 → 掩码后处理 → 栅格地图 → 覆盖路径 → 航点导出」的顺序把完整链路讲透新手能照着跑熟手可以直接拿去对参数。2. 为什么是分水岭算法选型理由与整套流程的数据准备2.1 分水岭原理距离变换之后的树冠缝隙就是天然分水岭分水岭算法把图像看成一张地形图灰度值当作海拔高度。局部极小值的位置是「谷底」水从谷底向外漫漫到两个盆地的水要汇合时算法在交界处筑起一道「坝」这条坝就是把不同区域分开的边界线。果树分割里要处理的问题正好反过来——树冠之间是缝隙缝隙里通常是土壤或杂草颜色和纹理都和树冠不一样如果直接在原图上跑梯度分水岭噪声点会变成无数个极小值结果全是被过分割的碎块。所以常见的做法是不在灰度图上跑而是在距离变换图上跑。先做一次「树冠 vs 背景」的二值化再对前景做距离变换每个前景像素的灰度值变成它到最近背景点的距离。这样一来每棵树的中心区域最亮往外逐渐变暗树冠之间的缝隙变成暗色「山脊」分水岭的坝正好落在缝隙上。这套思路我在多张正射图上验证过只要树冠之间有肉眼可见的间隙效果比直接阈值分割好得多。而如果两张树冠真的挨在一起、像素级别都连上了再好的分水岭也没法拆开这种情况要靠后处理或更细分的算法兜底。2.2 与U-Net、阈值法、超像素对比分水岭的性价比边界做果树图像分割绕不开三个方案传统阈值分割、U-Net深度分割、分水岭。阈值分割最轻像Otsu或颜色范围都能做但它的输出永远是「树冠 vs 背景」一张二值图两棵树挨着就粘连不给个体编号后续想生成逐棵树的路径还得再做连通域分析而连通域分析在粘连情况下又必然出错。U-Net这类语义分割模型效果上限最高能处理阴影、光照不均、部分遮挡可它需要几百张带像素级标注的训练图果园场景里标注一棵树要几分钟数据准备成本直接把大多数项目劝退。超像素分割SLIC只能把图切成均匀的小块本身不产生「树」的语义还得再叠聚类规则。分水岭恰好站在两者中间不需要训练算得也快输出是带编号的连通域每棵树天然就是一个独立区域正好喂给后面的无人机路径规划。它最大的短板是阈值和核参数会跟着图像分辨率、树冠形态变化换一个果园可能要重新调 0.4 这个系数。我的选择标准是如果正射图质量稳定、树冠之间有间隙直接用分水岭如果图里严重遮挡或阴影大面积压过树冠才考虑收集数据去训练 U-Net。这个对比可以看下面这张表四项指标的差距一目了然。方案是否需要标注抗粘连能力输出单位落地成本阈值分割不需要弱整张二值图最低分水岭不需要中等带编号的树冠区域低超像素 SLIC不需要中等无语义的小块中U-Net需要上百张掩码强逐像素类别高2.3 整套流程的五步拆解与数据准备把完整流程拆开看实际上是五个串在一起的环节图像预处理、距离变换与种子标记、分水岭分割、掩码后处理、航点生成。前四步输出一棵棵树的二值掩码第五步把掩码转成航线。路径生成属于二维覆盖路径规划是无人机路径规划算法里最基础的一类先把二维跑稳再谈三维避障。这个流程用伪代码看会更清楚# 流程伪代码从正射图到无人机航点 # 1. 读取正射图 img # 2. 预处理ExG 超绿指数 Otsu - binary # 3. 距离变换 取 0.4 峰值做种子 - markers # 4. cv2.watershed(img_rgb, markers) - labels # 5. labels1 生成掩码形态学清理 - tree_mask # 6. 把 tree_mask 转栅格地图牛耕扫描 - waypoints # 7. 导出 CSV校验转弯半径数据准备上我一般会确认三类文件到位正射图无人机可见光拼接图TIFF 或 JPG 都行、相机的坐标参考信息用于像素坐标转 GPS、以及一个存放航点结果的空目录。如果有条件额外拍几张不同光照下的测试图用来验证阈值参数的稳定性——这一步看着不起眼实际是后面把所有坑都暴露出来的关键。代码完整、数据齐全的工程通常还会带一个标注文件夹但分水岭这条路不强制要求标注有的话用来验收分割精度没有也不影响跑通。3. 用OpenCV实现果树图像分割预处理、种子标记与分水岭调用3.1 超绿指数与Otsu阈值把树冠从土壤背景里拎出来分水岭之前必须先得到一张尽可能干净的前景/背景二值图。直接取绿色通道是最快的办法但土壤颜色偏黄、阴影偏蓝只靠单通道会把枯叶和阴影一起算进来。我一般用超绿指数 ExG 2G - R - B思路是同时压低红色和蓝色把绿色树冠的亮度抬到最高。import cv2 import numpy as np # 读入正射图OpenCV 默认是 BGR 顺序 img cv2.imread(data/orchard_rgb.tif) # 超绿指数绿色通道权重最高红蓝同时压低 b img[:, :, 0].astype(np.float32) g img[:, :, 1].astype(np.float32) r img[:, :, 2].astype(np.float32) exg 2 * g - r - b exg_norm cv2.normalize(exg, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) # Otsu 自动选阈值光照不均时比固定阈值稳 _, binary cv2.threshold(exg_norm, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU)Otsu 的阈值完全由直方图自动决定如果发现树冠里有大面积空洞可以改成固定阈值配合形态学闭运算如果地面杂草太茂盛把阈值拉低就在 Otsu 结果上做一个面积过滤把小于几百像素的碎块先丢掉。这一步骤的核心目标不是完美分割而是让距离变换的输入足够「胖」——宁可树冠边界外扩一点也不要缺一块因为缺了的地方在距离变换里会变成假种子直接引发过分割。3.2 距离变换与种子标记决定过分割还是粘连的关键参数拿到二值图后接下来不能直接上分水岭先做距离变换。这一步的目的很明确让每棵树的核心区域变成独立的亮斑这样分水岭才能从每棵树的中心开始漫水在缝隙处停下。如果跳过距离变换直接在二值图上灌水第一棵树的边缘会把相邻树冠一起淹没输出就是一片连通的大区域后面路径生成就全乱套了。# 欧氏距离变换每个前景像素保存到最近背景点的距离 dist cv2.distanceTransform(binary, cv2.DIST_L2, 5) # 取最大距离的 0.4 倍作为种子阈值 _, seeds cv2.threshold(dist, 0.4 * dist.max(), 255, cv2.THRESH_BINARY) seeds seeds.astype(np.uint8) # 3x3 膨胀一次让种子覆盖树冠核心而不是只有中心点 kernel_seed np.ones((3, 3), np.uint8) seeds cv2.dilate(seeds, kernel_seed, iterations1) # 连通域编号每个种子获得唯一标记背景从0开始 _, markers cv2.connectedComponents(seeds) markers markers 1 # 把0号留给背景 markers[binary 0] 0 # 背景区域强制设为0注意markers 必须是 int32 类型的单通道图用 uint8 会让分水岭返回错误标签。0.4 这个系数是我在多数果园正射图上调出来的经验值它表示只有距离峰值 40% 以上的核心像素才会变成种子。树冠小而密的果园要降到 0.3避免一颗小树没有种子而被漏掉树冠大而疏的果园可以调到 0.5减少粘连风险。膨胀这一步不是必须的但我发现膨胀一次能让分水岭起点更稳定避免单像素种子在漫水时被噪声带偏。如果发现两颗相邻树共享一个种子就说明 0.4 还不够大需要继续上调到 0.45 或 0.5。3.3 分水岭调用与掩码后处理去掉碎块、填掉树冠空洞分水岭本身是一次函数调用但它的输入要求是 8 位三通道图加一个标记图markers标记图里每个像素是一个整数编号0 是背景1-50 是已知前景区域分水岭会从这些已知区域往外扩展。这里最容易犯的错是直接把单通道灰度图喂进去函数会报错或产生垃圾结果。# 分水岭需要 8 位三通道输入这里转成 RGB 语义 img_rgb cv2.cvtColor(img, cv2.COLOR_BGR2RGB) labels cv2.watershed(img_rgb, markers) # labels1 是背景1 代表不同的树个体 tree_mask np.zeros_like(binary) tree_mask[labels 1] 255 # 闭运算填掉树冠内部的黑色空洞 kernel_close cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5, 5)) tree_mask cv2.morphologyEx(tree_mask, cv2.MORPH_CLOSE, kernel_close, iterations2) # 过滤过小的碎块正射图分辨率不同阈值要按树冠像素面积估 contours, _ cv2.findContours(tree_mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE) min_area 800 valid_contours [c for c in contours if cv2.contourArea(c) min_area] # 重新画掩码保证每个树冠区域是独立连通域 tree_mask np.zeros_like(tree_mask) cv2.drawContours(tree_mask, valid_contours, -1, 255, -1)labels 1 的含义是「已知前景区域扩展后形成的分水岭区域」背景因为前面把 0 号让出来所以是 1。如果你发现所有树冠最终连成一个区域大概率是预处理二值图里树冠之间的缝隙已经没有了如果你发现每棵树碎成好几块说明种子阈值 0.4 开得不够高或者膨胀操作 iterations 太多导致一个种子跨过了缝隙。后处理里的 min_area 也是个需要跟着分辨率走的参数0.05 米分辨率下 800 像素对应 2 平方米刚好能滤掉田埂上的小杂草块。掩码里树木通常是圆形轮廓椭圆核对它更友好矩形核容易留下方角。4. 从分割掩码到无人机路径栅格地图、牛耕覆盖与航点导出4.1 掩码转栅格地图像素坐标怎么换算成实际距离分割完成后掩码还是一张像素图无人机飞的是地理坐标中间隔着一条「像素怎么变成米」的换算。如果正射图是带地理参考的 GeoTIFF常见做法是用 GDAL 读取地理仿射变换参数直接得到像素中心对应的经纬度如果是普通 JPG就需要知道飞行高度和相机参数或者用图像里已知物体的实际尺寸做比例换算。我在项目里更常用第二种因为很多飞手导出的并不是严格配准的 DOM硬套六参数反而出错。# 已知正射图覆盖 200m x 200m图宽 4000 像素 cover_m 200.0 width_px 4000 res cover_m / width_px # 0.05 m/pixel def pixel_to_metric(row, col): x col * res # 水平方向单位米 y row * res # 垂直方向单位米 return x, y提示res 的值最好从地图工具里的实际地面距离量出来不要靠飞行参数估估出来的误差会在航点间距上成倍放大。这个缩放系数 res 是整个路径生成里最不能拍脑袋定的参数。你从飞控日志或地图工具里量出覆盖尺寸后写进代码如果 res 差个 10%航点间的间距和实际作业行距就全偏了。用 GeoTIFF 时我会额外做一个校验取任意两个航点用仿射变换算出距离再和地图工具里的实测距离对比误差超过 1% 就直接停下检查。4.2 牛耕覆盖路径生成逐行扫出每棵树的作业航点有了栅格地图和二值掩码路径生成就是一个经典问题怎么让无人机覆盖所有树冠同时尽量少飞无效距离。果园场景下最常用的是牛耕式覆盖路径也叫弓字形路径无人机沿固定方向一行一行飞到头转弯再折返。对树冠掩码来说每一行只需要在掩码的左右边界之间穿行没有树冠的区域可以直接跳过不必像扫全田那样无脑走直线。def boustrophedon(mask, spacing_m4.0, res0.05): h, w mask.shape step_pix int(spacing_m / res) # 行间距换算成像素 waypoints [] direction 1 # 1 表示从左往右-1 表示从右往左 for top in range(0, h, step_pix): row min(top, h - 1) cols np.where(mask[row, :] 0)[0] if len(cols) 0: continue # 这一行没有树冠跳过 left, right cols[0], cols[-1] if direction 1: waypoints.append([row, left]) waypoints.append([row, right]) else: waypoints.append([row, right]) waypoints.append([row, left]) direction * -1 return np.array(waypoints)这段代码的关键参数是 spacing_m它决定了作业行距实际应该由无人机的喷幅或传感器有效幅宽决定不是随便设的。行距设成 4 米而无人机有效喷幅只有 3 米中间就会漏喷设成 2 米作业时间几乎翻倍。另一个容易忽略的点是代码里的 left 和 right 是树冠边缘而不是树冠中心如果你只想对树冠中心喷药可以把边界改成整行树冠像素的均值坐标让航点落在树冠串的中心线上路径会平滑很多。4.3 航点导出CSV格式、转弯半径校验与GPS坐标换算航点算出来后最后一步是导出成无人机能认的格式。我一般导出 CSV保留 lat、lon 两列再根据需要加 alt作业高度和 speed。这里有两个必须做的校验一是相邻航点距离不能太小否则无人机还没到目标点就要拐弯二是转弯角度不能太尖锐植保机以大角度转弯会掉速度甚至触发飞控保护。import csv import math def export_waypoints(waypoints, out_csv, alt3.0): with open(out_csv, w, newline) as f: writer csv.writer(f) writer.writerow([lat, lon, alt]) for row, col in waypoints: lat, lon pixel_to_gps(row, col) # 由仿射或控制点换算 writer.writerow([lat, lon, alt]) def check_turn_angles(waypoints, max_deg90): bad_indexes [] for i in range(1, len(waypoints) - 1): p0 np.array(waypoints[i - 1], dtypefloat) p1 np.array(waypoints[i], dtypefloat) p2 np.array(waypoints[i 1], dtypefloat) v1 p1 - p0 v2 p2 - p1 cos_a np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2) 1e-6) angle math.degrees(math.acos(max(-1.0, min(1.0, cos_a)))) if angle max_deg: bad_indexes.append(i) return bad_indexespixel_to_gps 这一步如果是带地理信息的正射图直接用仿射变换的六个系数算如果没有用两到三个地面控制点的经纬度和像素坐标做线性映射。坐标级别的偏移可以靠校验来兜底但你不能把校验也省了——我见过最离谱的一次GPS 航点整体偏了 12 米原因是地图工具里量的是屏幕像素不是图像像素缩放时换算系数用了两倍所以每次我都会先导出一个航点拿到地面实测验证过再批量飞。如果要做更完整的无人机三维路径规划还要在航点上叠加数字高程模型给每个航点加上地形高度约束二维覆盖路径是基础三维那套得另起一套纵向规划。5. 避坑与排查分水岭分割和路径生成里最容易翻车的五个点分水岭这套流程看着代码不长真正跑到无人机起飞那一步处处是坑。我把几个反复踩过的问题按「现象 → 原因 → 解决」写出来都是可以直接对着排查的经验。做分割和路径生成的项目八成时间都花在这些边界情况上提前知道坑在哪比多写一百行代码都管用。5.1 过分割一棵树被切成了三块现象连通域统计结果里一棵成年果树被切成两个或三个区域最直接的表现是树冠数量比肉眼数出来的多 20% 到 50%生成的路径会绕同一棵树飞两遍。原因分水岭的种子点数量远超实际树冠数。树冠内部因为光照不均或叶片稀疏距离变换图上会出现多个局部高点每个高点都成为种子漫水还没到树冠边缘就先在内部汇合了或者种子阈值 0.4 太松膨胀次数不够导致种子没有连成完整核心。解决把 dist 阈值从 0.4 往 0.5 调膨胀 iterations 由 1 改成 2再看连通域数量变化。如果数量降下来了说明种子合并到位如果降不下来就到预处理里补一步椭圆核闭运算把树冠内部的暗斑填平。这个 0.4 调起来多少有点玄学但每次只改一个参数、对比连通域数量比凭感觉反复跑全流程高效得多。5.2 树冠粘连两棵相邻的树分不开现象两棵树的树冠在掩码里是一个连通域不管怎么调分水岭参数都拆不开。原因分水岭的坝只能筑在「已有缝隙」的位置如果预处理二值图里两棵树之间的像素已经全是前景分水岭根本找不到坝址。常见情况是杂草蔓延到树冠间隙或者阴影把间隙压暗后又被形态学操作填掉了。解决先可视化 binary确认间隙还在不在。间隙被杂草占掉就把 Otsu 阈值改成手动阈值比如 exg_norm 大于 120 才算树冠间隙被阴影盖掉就降低飞行高度重新采集提高分辨率让轮廓露出来。两颗树冠真正贴到一起的极端情况二维图像分割已经没有足够信息可用可以对二值图先做一次小核腐蚀把狭窄连接处断开再跑距离变换这个后门不稳定腐蚀核的大小很难一次调对只能作为备选方案。5.3 掩码毛刺让航点抖动现象树冠边缘锯齿严重同一行里左右边界点来回跳导出的航点序列像心电图。原因二值化和分水岭输出都带椒盐噪声而牛耕扫描是逐行取最左最右的白色像素边缘上任何一个孤立白点都会把边界点拉偏。图像分割最常见的翻车点就是后处理没做干净。解决掩码后处理里加 MORPH_OPEN 和 MORPH_CLOSE先开后闭开运算核用 MORPH_ELLIPSE尺寸不超过树冠半径的 1/20。5x5 椭圆核是一个起点锯齿还在就换 7x7但核太大会把细小的树冠直接抹掉。我一般在核大小和树冠面积之间取平衡树冠越小核就要越小宁肯保留一点毛刺也别把树抹掉。5.4 相邻航点转角过大无人机按航线飞不过去现象掉头处转角接近 180 度无人机实际飞的时候会画一个大圆弧航线在地面站里看着是尖角飞控日志里频繁报「转弯半径不足」。原因牛耕路径每一行的入口和出口正好在树冠边缘相邻两条航线方向相反转弯点天然就是锐角行距越小转弯越急。这是覆盖路径规划算法本身的结构问题不是参数偶然触发。解决在掉头点外扩一个转弯半径的距离插入一个圆弧过渡航点或者把掉头点从树冠边缘挪到行外 2 到 3 米的空白区域再拐。这个值跟无人机平台强相关多旋翼可以小半径原地转固定翼必须给足半径。我一般先查飞控的「最小转弯半径」参数再按它的两倍设置过渡航点偏移量转弯半径检查脚本跑完再导出。5.5 坐标整体偏移GPS航点落地偏了几米现象导出的航点导入地面站后位置整体偏移轻则几米重则十几米飞手拒绝起飞。原因pixel_to_gps 的换算关系错位。最常见的有三种把图像像素宽度当成地面覆盖宽度用了地图工具的屏幕坐标而不是图像坐标或者只用一个控制点做比例换算完全忽略正射图可能存在的旋转矫正。解决至少选三个不在同一直线上的地面控制点用仿射变换求解六个系数控制点要选道路交叉口、田埂拐角这类长期存在的地物不要选可能被收割的庄稼边界。第一个导出航点先拿去地图上目视比对误差在 0.5 米内再批量生成。坐标转换没有捷径这是全流程里最该较真的环节也是我拿一次整体返工换来的习惯。6. 进阶技巧用HSV预处理提升鲁棒性再用连通域统计验证分割结果6.1 HSV颜色阈值做分水岭前的增强ExG 超绿指数对纯果园场景足够但混种果园里有些果树叶片偏黄或偏红超绿指数会把它们直接压成背景。我一般在预处理里加一条 HSV 辅助通道把 H 在 35 到 85 之间的像素标记为绿色候选S 在 40 以上的像素标记为有颜色的植被两者取并集后再和 ExG 结果做交集。这样黄色叶片也能被覆盖到代价是多几个可调阈值。预处理不是一步定死的事先跑一遍分割把漏检的树在图上标出来再回头微调 HSV 范围这个循环跑两轮就够了。6.2 连通域统计脚本自动检查分割结果分割结果靠肉眼盯图容易漏我习惯在生成路径前加一个自动检查统计连通域数量、每个连通域面积的中位数和极值把面积异常小的连通域直接过滤掉。这一步能提前发现过分割和粘连比飞完再返工便宜得多。def validate_mask(mask, min_area800): num, labels, stats, _ cv2.connectedComponentsWithStats(mask, 8) areas [s[4] for s in stats[1:]] # 去掉背景 median_area np.median(areas) max_area np.max(areas) print(f树冠数量: {len(areas)}, 中位面积: {median_area:.1f}, 最大面积: {max_area:.1f}) # 面积差超过 8 倍说明大概率有粘连或漏检 if max_area 8 * median_area: print(警告: 存在大面积连通域疑似树冠粘连)这个 8 倍阈值是经验值成龄果园树冠大小通常比较均匀差出 8 倍基本就是两棵或更多树连在一起了。拿到警告后再单独看看那片区域比从头检查全图省时间。我自己的习惯是分割参数每次换果园都要重调一轮但坐标校验和转弯半径检查从来不动这两条是流程的底线。最开始做这套方案时我在没有地面控制点的情况下直接批量生成了全园航点结果整体偏移把飞手吓得不轻从那以后所有 GPS 换算都坚持「先验证一个点再放量」。这套流程走到今天已经足够稳定分水岭选择的分辨率边界、牛耕路径的行距取舍、坐标转换的核对顺序每一个环节都有据可查。希望帮到你。本文还有配套的精品资源点击获取