简介：面向遥感图像中的道路提取任务，这份MATLAB实现提供了从预处理、特征提取到分割与后处理的完整参考流程，适合遥感、GIS及计算机视觉方向的初学者和研究者学习和二次开发。资源包共5个文件，均为.m脚本，包含主程序、不同测试图像的处理脚本及辅助验证函数，压缩包仅3KB，代码轻量直观，便于逐行理解道路识别的核心思路。目前已有537人学习下载。借助Canny边缘检测、形态学膨胀与腐蚀等经典算子，使用者可在MATLAB中快速复现道路区域分割效果，并参照脚本注释调整参数，观察不同地物背景下提取结果的变化。这份资料不仅是课程实验与毕业设计的实用示例，也可作为后续接入深度学习模型、提升高分辨率影像道路检测精度的研究起点。
1. 遥感道路提取不是边缘检测：这套MATLAB思路值得先看
拿到无人机正射影像或高分辨率卫星图，很多人的第一反应是“把道路提取出来”，然后下意识地打开Canny算子找边缘。结果往往是林缘、阴影、田埂、建筑轮廓全都被提了出来，路网反而埋在噪声里。这套基于MATLAB的道路提取资源走的不是边缘检测路线，而是更符合遥感影像特性的处理链：灰度化、滤波、阈值分割、形态学后处理、连通域几何过滤。每一段都有对应的脚本和参数注释，能跑通一条完整的“影像输入到路网mask输出”流程。适合两类人：一类是课程设计或毕设需要一套基线代码的在校学生，另一类是刚接触图像分割、想搞明白阈值和形态学实操组合的从业者。整个过程不依赖神经网络工具箱，普通MATLAB安装加Image Processing Toolbox就能跑，复现门槛很低。
2. 预处理这一步别省：灰度化、滤波与局部对比度增强
预处理是整条链路里最容易被跳过的环节，但绝大多数提取失败其实都发生在这一步。拿到一张遥感影像，先别急着调阈值，按顺序做三件事：降维、滤波、必要时局部增强。顺序反了效果会差很多，下面逐个展开。
2.1 灰度化：为什么RGB三通道不如单通道稳妥
道路在可见光影像里的灰度通常比周边植被和农田高，灰度化的目的就是把这种亮度差异集中到一个通道上。MATLAB里用rgb2gray一行就完成，它的底层权重是0.299R + 0.587G + 0.114B，绿色权重最高，这和人眼对亮度的感知曲线基本一致。
我见过不少新手直接对R、G、B三个通道分别做阈值分割，然后想融合结果。这个思路听起来更“充分利用信息”，实际工程里很麻烦：计算量变成三倍，每个通道的分割阈值和形态学参数都要单独调，三个结果不一致时还得分情况讨论要不要取交集。遥感道路提取任务里，道路相对背景的亮度差异通常足够大，一个灰度通道就能承载主要信息。真正需要回头用颜色信息的是那种道路和植被亮度接近但色差明显的场景，可以先在HSV空间看S通道或V通道的分离效果，再决定是否额外加一路分割。
% 读取影像，自动判断灰度还是彩色 img = imread('roads_scene.png'); if size(img, 3) == 3 gray = rgb2gray(img); % 彩色转灰度，权重0.299/0.587/0.114 else gray = img; % 本来就是灰度图就直接用 end这里的size(img, 3) == 3判断通道数，避免把单波段影像硬当RGB处理报错。转灰度后建议立刻看灰度直方图，确认是不是双峰分布：一个峰是背景（农田、林地），另一个峰是道路和人工地表。如果直方图已经明显双峰，后面OTSU分割会非常省事；如果拖出一条长尾巴，说明光照不均或阴影干扰严重，需要做2.3节的局部增强。
2.2 滤波：中值滤波在道路任务里优于高斯滤波
遥感影像的噪声来源很杂：传感器热噪声、JPEG压缩块效应、传输误码造成的椒盐点。这些噪声在阈值分割阶段会被成片放大，所以滤波不是可选项。任务里我一般直接用中值滤波，不推荐高斯滤波。
原因在于中值滤波保边，高斯滤波磨边。道路边缘在后续形态学处理里要靠结构元素“撑住”，边缘一旦被高斯模糊成渐变带，阈值分割后的边缘位置会漂移，骨架化时更容易产生毛刺。中值滤波对椒盐噪声的抑制效果突出，同时不会把边缘的锐度吃掉，这是它在道路提取场景里的核心价值。
% 中值滤波，窗口5×5 filtered = medfilt2(gray, [5 5]);窗口尺寸的选择逻辑：[3 3]能滤掉孤立噪点，但影像若压缩严重或噪声密集，滤不干净；[7 7]开始侵蚀道路边缘细节，细小的乡村道路可能在滤波这步就被抹平了。[5 5]是绝大多数中分辨率遥感影像的安全选择。如果图像尺寸特别大，比如1万像素级别的正射影像拼图，可以考虑先降采样再做中值滤波，效率高很多。
2.3 光照不均：局部直方图增强处理阴影段
高分辨率影像里建筑阴影和树荫经常把道路盖住一段，全局阈值分割时这些暗部区域会被直接划到背景类，路网断口从此产生。这里有一个简单但有效的处理：对滤波后的图像做局部直方图均衡，把阴影里被压暗的灰度重新拉开。
% 局部直方图均衡，增强暗部细节 enhanced = adapthisteq(filtered, 'ClipLimit', 0.02, 'NumTiles', [8 8]);adapthisteq把图像分成若干小块分别做直方图均衡，ClipLimit控制对比度增强强度，默认0.01，激进一点用0.02；NumTiles设为[8 8]意味着分成64块，块越小局部适应性越强，但块太小会让同质区域出现假纹理。这个函数对光照缓慢变化的影像效果明显，处理完再看直方图，原来挤在低灰度区的道路像素会被摊开，和背景的峰更容易分开。
需要注意一个红线：局部增强必须在滤波之后做。如果先adapthisteq再中值滤波，噪声也被同步放大，中值滤波反而要加大窗口才能压住，边际成本不划算。这个顺序问题我在项目里翻过车，后面会专门写进避坑清单。
3. 分割与形态学三件套：把道路从影像里真正切出来
预处理做完，接下来的主流程就是分割加形态学处理。这套流程是资源的绝对核心，或者说整份代码的价值都在这一章里。思路是OTSU全局阈值、开闭运算去噪补断口、连通域几何过滤，三层下来道路掩膜基本像样。
3.1 OTSU阈值：全局阈值是基线，但不是终点
OTSU又称最大类间方差法，核心思路是找一个灰度阈值，让分割后前景和背景两类之间的方差最大，也就是说这个阈值两侧的像素尽可能“类内抱团、类间分离”。MATLAB里用graythresh计算归一化阈值，配合imbinarize完成二值化。
% OTSU阈值分割 level = graythresh(enhanced); bw = imbinarize(enhanced, level);graythresh返回一个0到1之间的归一化值，不是绝对灰度值。imbinarize接收它之后会把图像中大于阈值的像素置为前景。这里有个版本坑：imbinarize是R2016a引入的函数，老代码里用的是im2bw(x, level)，如果脚本报“未定义函数”一类错误，优先查这处替换关系。OTSU在处理光照均匀的影像时足够好，但遇到全局光照梯度大的影像，单个阈值必然顾此失彼。这种情况我会退回2.3节把局部增强参数调强，让阈值有更大可能分开两类；如果还不行，就得考虑分块阈值，把图像切成网格分别计算局部阈值再拼接。
3.2 形态学开闭运算：结构元素尺寸要贴合道路形状
阈值分割后的二值图可用性很差，零散噪点、破碎道路段、孤立亮块交织在一起。形态学操作是这一段的真正主力。标准的做法是“先开后闭”：开运算去掉小亮块和碎点，闭运算接回道路断口。
% 开运算去碎点，用矩形结构元素贴合道路形状 seOpen = strel('rectangle', [15 3]); bwOpen = imopen(bw, seOpen); % 闭运算接路网断口，用圆形结构元素 seClose = strel('disk', 5); bwClose = imclose(bwOpen, seClose); % 删除面积小于500像素的孤立碎块 bwClean = bwareaopen(bwClose, 500);结构元素的选择是这段代码里最值得讲的地方。开运算用的是[15 3]矩形，长边沿道路方向、短边沿垂直方向。道路的特点是长条状，矩形结构元素在腐蚀阶段不会把道路主干腐蚀断，但能把圆形的建筑亮块、小面积噪点消掉。如果用同尺寸的disk结构元素，腐蚀半径3就会把宽度只有几个像素的乡间小路整段消没，这是很多提取任务“路没了、房子还在”的直接原因。
闭运算用disk(5)而不是矩形，是因为断口方向不确定，圆形在任意方向都能桥接。这里有个参数联动关系：disk(5)能接回的断口宽度上限大约是10个像素，超过这个宽度桥不上。如果阈值分割后断口普遍很宽，优先回去调阈值或增强，不要硬把disk半径加大到10以上，那会把旁边的非道路区域也糊进来。
bwareaopen的500像素阈值是基线，实际使用时根据影像分辨率调整。无人机影像的建筑物动辄几千像素，500以下的基本是噪声；但如果是高分辨率卫星图，道路像素宽度成倍增加，阈值也可以相应调高。建议分段统计regionprops的面积分布后再定，不要拍脑袋。
3.3 连通域几何过滤：用手工规则把非道路踢出去
形态学做完，剩下的连通域里通常还混着河流、田埂、建筑群。这些地物在阈值分割阶段因为灰度接近进入了道路掩膜，此时需要几何特征过滤。
% 连通域标记，8邻域连接 cc = bwconncomp(bwClean, 8); stats = regionprops(cc, 'Area', 'Eccentricity', 'Solidity', ... 'MajorAxisLength', 'MinorAxisLength', 'Orientation'); % 面积大于800像素，且偏心率大于0.9，才认为是道路段 idxKeep = find([stats.Area] > 800 & [stats.Eccentricity] > 0.9); keepMask = ismember(labelmatrix(cc), idxKeep);关键参数是两个：面积阈值和偏心率（Eccentricity）。偏心率描述区域接近直线的程度，0是正圆，1是线段。道路段通常是细长条，偏心率轻松超过0.9；而建筑群、广场这类块状地物偏心率低，这一刀能切掉不少干扰。第二个可用的特征是Solidity，区域面积和凸包面积的比值，道路段因为边缘相对平滑，Solidity一般也较高。如果结果里混入太多短粗目标，可以再加MajorAxisLength / MinorAxisLength > 3的长宽比约束。
这段跑完，路的雏形就出来了。但离“可以用的路网”还有距离，骨架化和断点缝合是下一章的内容。
4. 从掩膜到路网：骨架化、毛刺清理与断口拼接
分割结果在视觉上是“一条一条的亮带”，但无法直接用于GIS矢量化和道路长度统计，因为没有人想数着一堆像素宽度做分析。这章的核心是把亮带降维成单像素宽度的网络骨架，同时处理毛刺和断口，让路网更接近真实拓扑结构。
4.1 骨架化：thin和skel不是一回事
MATLAB里bwmorph有两个常用操作：'thin'和'skel'。'skel'生成的是最大内切圆圆心轨迹，保留了更多形态学信息，但在道路场景里容易出多余的分叉；'thin'是迭代去除边缘像素但不破坏连通性，结果更干净，更适合路网提取。我习惯优先用'thin'，迭代次数给Inf让它跑到收敛。
% 提取路网骨架，Inf表示迭代到不能再变细为止 skel = bwmorph(bwClean, 'thin', Inf); % 去掉长度不超过5个像素的毛刺分支 skelClean = bwmorph(skel, 'spur', 5);'spur'的第二个参数是允许去除的分支长度阈值，单位是像素操作迭代次数。设成5表示“长度小于等于5个像素的小分叉直接去掉”。这个值不宜过大，因为道路交叉口的短连接长度可能也就十来个像素，设到20以上会把交叉口切碎，路网变成一堆断线。
骨架化顺序有个硬规则：必须等前序操作全部完成后再骨架化，因为bwmorph的thin是破坏性操作，之后再做开闭运算很可能会把已经细化的骨架再次撕裂。如果骨架化后发现毛刺仍然严重，说明输入掩膜边缘不够光滑，应当回头调整闭运算的结构元素尺寸，而不是一味加大spur参数。这是我反复强调的一点，很多人卡在毛刺上反复调spur，却忽略了根源是掩膜边缘毛糙。
4.2 断口缝合：方向和距离两个约束缺一不可
骨架化之后依然会有断口：树荫遮挡、阈值分割遗漏、后处理误删，都会让路网断成若干段。一个工程化的做法是检测端点、配对、补线。端点检测在MATLAB里很简单：
% 检测骨架端点 endPoints = bwmorph(skelClean, 'endpoints'); [rows, cols] = find(endPoints);拿到端点坐标之后，配对逻辑不能只看距离，还要看方向。两个端点的连线方向和各自局部的道路走向一致，才值得缝合；平行但错开很远的断口不应该强行连。一个合理的判断条件是：两点距离小于N个像素（N建议取30到50，按影像分辨率折算），且连线方向与两点局部切线方向的夹角都小于30度。
这里不展开完整缝合代码，因为在真实项目里我通常只缝合能量化的断口，剩下的交给人机交互手动补。骨架自动缝合的边界条件太多，田间小路断口毫无规律时，盲目缝合会画出虚构道路，这比留一段缺口更影响后续分析的可信度。
4.3 结果验证：用IoU量化提取效果
很多人在这一步就停了，看结果“感觉还行”就交了。工程上需要量化验证，最常用的指标是IoU，即预测掩膜和人工标注掩膜的交并比。
% 计算与人工标注的IoU inter = sum(maskGt(:) & maskPred(:)); union = sum(maskGt(:) | maskPred(:)); iou = inter / union;maskGt是手工标定的标准mask，maskPred是算法提取的mask。IoU的计算结果在0.5以上就可以看作基线可用的提取效果，0.7以上已经具备实际应用价值。要留意的是，计算IoU前确认两个mask尺寸一致、坐标系对齐，否则算出来的数字没有意义。另外建议把IoU拆成“完整性”和“正确性”两个方向看：查全率低说明断口多，精度低说明误检多，不同问题对应的修复手段完全不同。
5. 避坑指南：五个把结果做砸的典型现场
这条路我走过的坑不算少，下面五条是反复出现的，每一条都按“现象 → 原因 → 解决”来写，直接可以对照排查。
坑一：河流被当成道路提了出来
现象：分割结果里河流和道路同时在mask里，且河流段长、连续性高，几何特征和道路几乎无法区分。
原因：水体在可见光影像中和浅色混凝土路面灰度接近，OTSU没有区分能力；加上Eccentricity和长宽比都是细长形状，几何过滤也拦不住。
解决：如果源数据是多光谱，直接额外计算NDWI水体指数把河流抑制掉；如果只有RGB，用局部方差区分——水面灰度极其均匀、纹理几乎为零，而道路路面有纹理变化。把stdfilt的局部标准差作为辅助特征，凡是标准差极低的细长连通域优先判定为水体。
坑二：阴影把道路切成碎块
现象：道路掩膜呈现“虚线”状态，断口全部集中在建筑物阴影或树荫覆盖区。
原因：OTSU全局阈值把阴影归入背景类，被阴影覆盖的道路像素亮度整体塌陷，不会出现在前景里。
解决：分两步。先回到2.3节做局部增强，提高阴影区的动态范围；再做闭运算接断口。如果断口宽度超过了结构元素能力，考虑对阴影区域单独拉低阈值做第二次分割，把两次结果做并集，能明显改善完整性。
坑三：田埂和地块边界混进来了
现象：农田区域密集的田埂被提取成密集平行线，路网被淹没。
原因：田埂的宽度、亮度、方向特性和乡村道路相似，纯几何特征无法区分。
解决：加纹理特征。道路周边的局部灰度方差通常高于田埂，或道路两侧存在明显的人工地表过渡带。也可以利用宽度一致性：同一条真实道路的宽度在影像上变化不大，而田埂宽度过于均匀。用bwdist计算骨架到掩膜边界的距离分布，方差小的细长区域更可能是田埂，这一招在我的项目里救回来好几次。
坑四：脚本在另一台机器上直接报错
现象：代码在自己电脑上跑通，换到同事机器上imbinarize报错“未定义函数”。
原因：R2016a之前没有imbinarize，之前版本都是im2bw。另外部分机器没装Image Processing Toolbox，regionprops、bwmorph全部不可用。
解决：脚本开头加功能介绍和版权说明不解决实际问题，真正要做的是写一层兼容封装：检测MATLAB版本号，旧版本自动走im2bw，或者统一替换成img > level * 255这种原始判断。新脚本里我一般直接写bw = enhanced > level * 255，绕开版本差异，这在任何版本都稳定。
坑五：骨架化后毛刺多到没法看
现象：骨架线跟树枝一样，分叉密布，路网拓扑完全不可读。
原因：掩膜边缘不光滑，微小突起在thin迭代时被放大成分叉，spur阈值设小没有清干净。
解决：先回掩膜层面做一轮闭运算，把边缘平滑掉再骨架化。spur参数从5开始往上试探，但不能一刀切得过大。如果毛刺是有规律的出现在同一方向，比如影像中农田边界方向一致，考虑停止自动骨架化，改用Hough变换提取主方向线段重建路网，这个降级方案虽然费事，但结果可靠得多。
6. 让提取结果更稳：连通域特征向量与轻量分类决策
形态学和几何过滤都做完后，结果质量基本取决于阈值定得准不准。如果追求更高的鲁棒性，有一个不依赖深度学习也能明显提升效果的做法：把人工规则的“多条件AND判断”升级成轻量分类器。做法是先用regionprops把每个连通域的多维特征导出成一张表，再用手工标注几十个样本喂给KNN或线性判别分类器，代替写死的阈值规则。
% 对每个连通域提取形态学特征向量 cc = bwconncomp(bwClean, 8); stats = regionprops(cc, 'Area', 'Eccentricity', 'Solidity', ... 'MajorAxisLength', 'MinorAxisLength', 'Orientation'); feat = table(); feat.Area = [stats.Area]'; feat.Ecc = [stats.Eccentricity]'; feat.Solidity = [stats.Solidity]'; feat.LenWidRatio = [stats.MajorAxisLength]' ./ [stats.MinorAxisLength]'; % 训练KNN分类器，trainIdx是手工标注过的样本索引 mdl = fitcknn(feat(trainIdx, :), label(trainIdx), 'NumNeighbors', 5); % 在测试样本上验证效果 pred = predict(mdl, feat(testIdx, :)); acc = mean(pred == label(testIdx));fitcknn的NumNeighbors取5是常见折中，邻域太小容易过拟合，太大则会把边界样本分错。手工标注训练样本时，正负样本比控制在1:1附近，每类标20到30个就够启动。跑完看分类准确率，如果低于85%，回头查特征分布，找出重叠区域再做一次特征组合调整。实际项目里这一招能把IoU再拉高5到10个百分点。
这套资源的代码链跑通之后，我个人的一个习惯是：每次换新影像，先跑一遍预处理，再分别保存阈值分割前、形态学后、过滤后三张中间图，逐张检查是哪一个环节损失了道路信息。这种“逐层检查中间产物”的做法帮我避免了很多盲目调参的玄学时间。记得有一次只要检查到某一步就定位问题，而不用把阈值从头试到尾，从那以后我每次换数据集都会强制走一遍这个流程，希望帮到你。
本文还有配套的精品资源，点击获取
