Python遥感图像道路提取：从U-Net到IoU评估的完整实践

简介这是一套基于Python的遥感图像道路提取算法课程设计源码面向高校计算机、遥感及相关专业学生也适合需要完成课程设计或期末大作业的编程学习者。压缩包共64个文件核心为35个py源码文件同时配有pyd/pyc编译模块、png测试图像、txt说明文档、打包配置文件及备份文件整体大小仅3.97MB目录按算法、界面、检测策略、测试等模块划分便于按需查阅。目前已有30人学习使用。项目完整覆盖从影像数据读取、灰度共生矩阵特征提取、K均值聚类、多种道路检测策略到结果可视化输出的全流程并提供了主窗口交互界面与辅助工具关键代码附有详细注释所有模块经过系统化测试可在典型环境中稳定运行。无论是用于理解遥感图像处理与机器学习算法还是作为毕业设计、课程实践的二次开发基础都能提供清晰完整的参考。1. Python 遥感图像道路提取一份能直接交课程设计的源码包遥感图像道路提取是高校计算机、测绘和地信专业课程设计里出现频率最高的题目之一。很多同学从 GitHub 拉一个 U-Net 工程跑通后把预测结果图往报告里一贴答辩时被老师问一句「训练样本怎么做的、精度指标怎么算的」就愣住了。也有同学自己从零写被 GDAL 安装、影像裁剪、标签制作几个硬骨头卡了一周还在原地。这份基于 Python 的遥感图像道路提取源码包把完整的算法链路串起来了高分影像预处理 → 训练样本制作 → OpenCV 传统算法快速原型 → U-Net 深度学习分割 → IoU 指标评估 → 栅格矢量化。它面向两类人一类是正在做课程设计、期末大作业、需要尽快跑通一条管线的本科生另一类是刚入门遥感图像分割、想找一个完整参考工程的算法爱好者。下面我会按「预处理 → 算法选型 → 避坑 → 验证加分」的顺序把每个环节的参数设置和踩坑点讲透。2. 影像预处理才是重头戏数据选型与 Rasterio 裁剪拼接很多课程设计失败不是死在模型上而是死在数据准备上。遥感影像不是普通的 JPG 图片它自带地理坐标系、辐射定标参数和多个波段如果直接拿 OpenCV 去读读出来的是一个形状height, width, bands的数组但像素值和真实地物反射率的关系你可能完全不知道。这一步不讲清楚后面所有算法都是建在沙子上。2.1 数据源怎么选高分影像参数与公开数据集课程设计里最常用的国产高分数据是高分一号的 PMS 传感器和高分二号的 PMS 传感器。高分一号全色分辨率 2 米、多光谱 8 米幅宽 60 公里适合做大范围路网提取高分二号全色 0.8 米、多光谱 3.2 米对单条道路的边缘细节保留得更好但幅宽只有 45 公里且单景数据量更大如果你的机器只有 8G 内存处理起来会比较吃力。两者的原始数据都是包含 RPC 文件的 TIF 格式部分数据还附带有理多项式校正参数需要在预处理阶段处理。数据源全色分辨率多光谱分辨率波段数适用场景高分一号 PMS2m8m3R/G/B大范围路网拓扑提取高分二号 PMS0.8m3.2m3R/G/B城市精细道路边缘分割Massachusetts Roads1.2m—3公开测试方便跑基线DeepGlobe Road0.5m—3竞赛数据样本质量高我给你的第一条建议是课程设计优先用公开数据集跑通全流程再用本地高分影像做验证。Massachusetts Roads 数据集是美国马萨诸塞州的航空影像标注是 1 像素宽的道路中心线DeepGlobe 道路提取挑战赛的数据是 0.5 米分辨率的卫星图标注更接近道路面。这两个数据集的共同好处是已经切割成 512×512 或 1024×1024 的瓦片省去了大影像裁剪的麻烦。如果导师要求必须用国产高分数据你需要向学校的数据中心申请或者用公开渠道下载已经做过正射校正的产品。拿到手后第一步是看元数据用 GDAL 打开看一眼投影坐标系、像元大小和波段数再决定是直接裁剪还是先做全色与多光谱融合。课程设计阶段我一般不推荐做融合因为 Gram-Schmidt 融合需要专门的软件或者用 GDAL 的 pansharpen 子集处理不好会出现光谱畸变反而影响后续分割。2.2 裁剪与波段合成Rasterio 实战拿到原始影像后第一步是按研究区或按网格裁剪出训练瓦片。这里我一般用 Rasterio 而不是 GDAL 命令行原因很简单Rasterio 的窗口读写机制可以按像素坐标直接切出一块不需要先把整景影像读进内存处理几个 G 的 TIF 也不会把内存占满。import rasterio from rasterio.windows import Window import numpy as np def crop_by_pixel(src_path, dst_path, col_off, row_off, width, height): 按像素坐标裁剪影像局部区域 col_off: 起始列row_off: 起始行width/height: 裁剪宽高 with rasterio.open(src_path) as src: window Window(col_off, row_off, width, height) transform src.window_transform(window) profile src.profile.copy() profile.update({ height: height, width: width, transform: transform }) with rasterio.open(dst_path, w, **profile) as dst: dst.write(src.read(windowwindow)) # 示例从第 1000 行 1000 列开始裁剪 512x512 的瓦片 crop_by_pixel(gf2_pms.tif, tile_1000_1000.tif, 1000, 1000, 512, 512)这个函数的核心是Window对象。注意Window的构造参数顺序是(col_off, row_off, width, height)先列后行很多人第一次写会写成(row_off, col_off)结果裁出来的区域在空间上错位了几百个像素而且因为都是整数不容易发现。裁剪时建议把瓦片尺寸设成 256 或 512 的整数倍后面喂给深度学习模型时不需要再做 resize避免分辨率变化带来的地物形变。裁剪出来的多波段 TIF 不能直接用于模型训练因为模型输入的通道数需要固定且像素值需要归一化到 0~255 或 0~1。下面的函数把多波段数据转成 8 位 RGB 图同时做了百分比线性拉伸——这是遥感图像显示里最常用的增强方式。def to_8bit_rgb(array_2d_percent2): 将多波段数组转为 8 位 RGB按百分比拉伸增强对比度 array: 形状为 (bands, height, width) 的 ndarray lower/upper: 拉伸百分位默认取 2%~98% bands array.shape[0] out np.zeros((array.shape[1], array.shape[2], 3), dtypenp.uint8) for i in range(min(bands, 3)): band array[i] p_low, p_high np.percentile(band, [lower, upper]) band_stretched np.clip((band - p_low) / (p_high - p_low 1e-6), 0, 1) out[:, :, i] (band_stretched * 255).astype(np.uint8) return out这里np.percentile的计算逻辑是先统计整景影像像素分布的 2% 和 98% 分位数然后把中间部分线性拉伸到 0~255。为什么要用百分位而不是直接用最大最小值因为遥感影像里常有云、高亮屋顶、水体等极端像元如果按整个动态范围拉伸整个图像会显得灰蒙蒙的道路和背景的对比度反而下降。如果你发现转出来的图片整体发黑或者亮成一片第一反应就应该是检查这个拉伸参数而不是去调模型。另一个常被忽视的预处理是影像归一化统计。不同景的影像因为拍摄时间、太阳高度角、大气条件不同亮度和对比度差异很大。常见的做法是对所有训练瓦片统计均值和标准差训练时用(img - mean) / std做标准化预测时也用训练集的统计值而不是每张测试图单独算。这是很多课程设计里测试集精度暴跌的根源后面避坑章节我会再展开。3. 两条算法路线从 OpenCV 形态学到 U-Net 深度学习道路提取的算法选型直接决定你报告的含金量和答辩时被提问的难度。如果只做传统算法工作量少但精度上限低如果做深度学习训练时间和机器要求高但指标好看。我的建议是两条线都做传统算法当 baseline深度学习当主结果这样的对比实验结构在课程设计答辩里几乎不会被问倒。3.1 传统算法为什么在遥感道路上失灵先说结论单独的 Canny 边缘检测、OTSU 阈值分割或 Hough 直线检测在遥感图像上都做不出完整的道路网。原因有三个。第一道路在影像上的光谱特征并不唯一。新建的沥青路是深灰色老旧水泥路偏白土路是黄褐色阴影下的道路和阴影下的农田光谱非常接近单纯按灰度阈值分割会把农田和建筑屋顶一起提出来。OTSU 这类全局阈值方法在道路占比小的影像上尤其不稳定因为算法假设图像只有前景和背景两类遥感图里的植被、水体、阴影、建筑会让直方图变成多峰分布。第二城市道路是弯曲的有交叉口、环岛和立交桥Hough 直线检测只能提取直线段遇到曲线路段需要分段拟合参数调起来非常繁琐。而且行道树遮挡会让边缘断裂断裂后的短线段数量巨大长度筛选阈值设多少都会顾此失彼。第三阴影区是传统算法的重灾区。高分辨率影像里建筑物阴影的灰度值和沥青路面接近边缘检测会在阴影边界产生大量虚假响应。如果你把 Canny 的阈值调高来压制阴影噪声道路边缘也会被一起滤掉阈值调低结果图里全是密密麻麻的短线段。所以传统方法的正确定位是作为快速原型验证数据质量以及作为深度学习结果的对比实验。跑通一条完整的传统算法只需要半小时能让你在训练深度学习模型之前确认影像预处理和标签是否对齐。3.2 路线 ACanny 边缘检测 形态学闭运算做快速原型传统路线我会用 OpenCV 完成核心思路是边缘检测后用形态学闭运算把断裂的道路边缘连接起来再用面积筛选去掉碎斑。这个流程对清晰的城市主干道有效对乡村土路和阴影密集区效果有限但作为 baseline 足够。import cv2 import numpy as np # 读取上一章生成的 8 位 RGB 瓦片 img cv2.imread(tile_8bit.png, cv2.IMREAD_GRAYSCALE) # 1. 高斯模糊抑制噪声核大小取 5x5sigma 默认 0 blur cv2.GaussianBlur(img, (5, 5), 0) # 2. Canny 边缘检测双阈值设为 50 和 150 edges cv2.Canny(blur, 50, 150) # 3. 形态学闭运算用 5x5 椭圆核连接断裂边缘 kernel cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5, 5)) closed cv2.morphologyEx(edges, cv2.MORPH_CLOSE, kernel, iterations2) # 4. 连通域分析删除面积小于 500 像素的碎块 num_labels, labels, stats, centroids cv2.connectedComponentsWithStats(closed) road_mask np.zeros_like(closed) for i in range(1, num_labels): if stats[i, cv2.CC_STAT_AREA] 500: road_mask[labels i] 255 cv2.imwrite(baseline_mask.png, road_mask)Canny 的两个阈值是这组参数里最关键的部分。阈值对50, 150的意思是梯度幅值超过 150 的像素确定为边缘低于 50 的丢弃介于两者之间的像素只有当它与高阈值边缘相连时才保留。遥感影像噪声明显时我会把低阈值抬高到 60高阈值抬到 180减少阴影边缘的虚检。形态学核的尺寸决定了断口连接的能力5×5 的椭圆核对 2 米分辨率的道路来说偏保守如果你想连更长的断口把核改成 7×7、迭代次数保持 2 次即可但要注意迭代太多次会把建筑物轮廓也连成片。面积阈值 500 像素意味着过滤掉所有小于约 2000 平方米的对象按 2 米分辨率算这个值要按你影像的分辨率折算。如果你用的是 0.5 米分辨率的 DeepGlobe 数据500 像素只相当于 125 平方米可能保留大量短线段需要调高到 2000 甚至 5000。这组参数没有万能值我在交付工程里用了配置文件的方式让每个数据集单独维护一组参数避免改数据就要改代码。3.3 路线 BU-Net 语义分割的实现要点深度学习路线我推荐用 U-Net 而不是 DeepLabV3 或 PSPNet原因是 U-Net 的编码器-解码器结构对中低分辨率影像的小目标分割更友好跳跃连接能把浅层边缘信息直接传给解码器训练样本需求也相对低。课程设计的样本量通常只有几百到一千张瓦片这个规模下 U-Net 比大模型更容易收敛。下面是一个精简版的 U-Net 实现输入是三通道 RGB 图像输出是单通道的道路概率图未经过 sigmoid 的 logits模型结构上只包含两个下采样块和一个上采样块实际工程里可以按这个模式继续加深到四层。import torch import torch.nn as nn class ConvBlock(nn.Module): def __init__(self, in_ch, out_ch): super().__init__() self.conv nn.Sequential( nn.Conv2d(in_ch, out_ch, 3, padding1), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue), nn.Conv2d(out_ch, out_ch, 3, padding1), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue)) def forward(self, x): return self.conv(x) class UNet(nn.Module): def __init__(self, in_ch3, out_ch1): super().__init__() self.enc1 ConvBlock(in_ch, 64) self.enc2 ConvBlock(64, 128) self.pool nn.MaxPool2d(2) self.mid ConvBlock(128, 256) self.up2 nn.ConvTranspose2d(256, 128, 2, stride2) self.dec2 ConvBlock(256, 128) self.up1 nn.ConvTranspose2d(128, 64, 2, stride2) self.dec1 ConvBlock(128, 64) self.out nn.Conv2d(64, 1, 1) def forward(self, x): e1 self.enc1(x) e2 self.enc2(self.pool(e1)) m self.mid(self.pool(e2)) d2 self.dec2(torch.cat([self.up2(m), e2], dim1)) d1 self.dec1(torch.cat([self.up1(d2), e1], dim1)) return self.out(d1)这个模型的关键设计在于跳跃连接——上采样后的特征图与同层编码器的输出在通道维度上拼接torch.cat让解码器同时拥有高层语义信息和低层边缘纹理。注意拼接后通道数翻倍所以dec2和dec1的输入通道数分别是 256 和 128。ConvBlock里的BatchNorm2d在遥感影像这种分布波动大的输入上很重要它能稳定训练过程没有它的话模型可能需要三倍以上的 epoch 才能收敛。训练时的参数选择直接影响结果好坏我一般按下面的配置起步criterion nn.BCEWithLogitsLoss() optimizer torch.optim.Adam(model.parameters(), lr1e-4) scheduler torch.optim.lr_scheduler.StepLR(optimizer, step_size10, gamma0.5) # 训练循环核心逻辑 for epoch in range(30): model.train() for img, mask in train_loader: pred model(img) # 输出 shape: (B, 1, H, W) loss criterion(pred, mask) # mask 为 0/1 的二值标签 optimizer.zero_grad() loss.backward() optimizer.step() scheduler.step()BCEWithLogitsLoss把 sigmoid 和交叉熵合并在一起计算数值上比「先 sigmoid 再 BCE」更稳定所以pred不需要手动过 sigmoid。学习率 1e-4 是 Adam 在图像分割任务里的常见起点不建议一上来就开 1e-3遥感影像的梯度噪声大学习率过高会让 loss 在前几个 epoch 剧烈震荡。StepLR 每 10 个 epoch 把学习率减半训练后期做精细调整。如果你的显存不够把输入 patch 从 256 降到 128 比强行压缩 batch 更有效因为 128 的输入显存占用只有 256 的八分之一。训练样本和标签的配对是另一门学问。前面提到 Massachusetts Roads 的标签是 1 像素宽的中心线这类标签直接训练会出问题道路在 ground truth 里占的面积比例极小模型会学成「什么都预测为背景」因为这样 loss 已经很低了。解决方式是训练前对标签做膨胀把中心线扩展成宽度约 5~7 像素的道路面或者在 loss 里给前景类别加权重。后面的避坑章节我会再详细说这个过程。4. 避坑排查道路提取中五个让人反复翻车的细节这部分是血泪经验。我拆过不下十份道路提取相关的课程设计源码也帮人调试过各种跑不通的工程最后发现翻车点高度集中在数据准备和结果后处理两个阶段。下面按「现象 → 原因 → 解决」的格式整理出五个最常见的坑。4.1 数据与训练阶段的三个坑第一个坑训练时显存溢出程序直接崩溃。现象是 U-Net 训练到第二个 batch 时报CUDA out of memory很多同学的第一反应是把 batch size 从 8 减到 2结果发现还是溢出。原因通常不是 batch 太大而是输入瓦片尺寸设成了 512×512在特征图最深的层显存占用会异常高。解决方法是优先把 patch 降到 256×256显存占用直接降到四分之一如果还不行再把 batch 降到 4并用torch.cuda.amp混合精度训练。我在源码包里给了默认的 256 patch 配置16G 显存的卡能稳定跑 batch 8。第二个坑模型在验证集上的 IoU 很高但换上另一景影像预测时结果一塌糊涂。现象是训练集和测试集来自不同传感器或不同季节的影像模型提取的道路边缘模糊、漏检严重。原因是遥感影像的辐射差异——不同成像时间的太阳高度角、大气条件、传感器增益都会让同一地物的像素值不同模型学到的颜色特征在新的辐射环境下失效。解决方法是训练时做数据增强对亮度、对比度做随机扰动模拟不同成像条件下的辐射差异预测前用训练集的均值和标准差对测试影像做标准化而不是让模型面对完全陌生的数值分布。第三个坑标签是 1~2 像素宽的中心线训练后模型预测结果几乎全黑。现象是 loss 下降很慢最终准确率高但以预测「全背景」实现。原因是正负样本极度不平衡道路像素占总像素比例不到 1%模型只要输出全零就能把 loss 压得很低。解决方法是训练前对标签做形态学膨胀用 3×3 或 5×5 的核把中心线扩展成道路面。遥感图像标注本身是费眼睛的活用 QGIS 手工标注一张 512×512 的图大概要 15 分钟如果标注的线细膨胀这一步省不掉。注意膨胀核的直径要和道路真实宽度匹配2 米分辨率的影像上 5×5 核膨胀后的道路宽度约为 10 米基本合理。4.2 结果后处理阶段的三个坑第四个坑预测结果图里道路断成一截截的中间有大量细碎孔洞。现象是模型输出的二值图里道路区块之间不连续看起来像虚线。原因是模型是逐像素分类的本身不感知道路连通性再加上训练样本里树木遮挡和车辆停泊造成的局部遮挡模型在这些区域会犹豫。解决方法是预测后先做形态学闭运算用 5×5 椭圆核把断口连接起来再做连通域分析只保留面积最大的前 N 个连通域。闭运算的迭代次数 2 次即可太多会把相邻的平行道路粘连成一块。第五个坑栅格转矢量后提取结果是一堆碎斑而不是光滑的路网线。现象是 ArcGIS 或 QGIS 打开生成的 shapefile发现多边形支离破碎边界锯齿严重。原因是预测概率图直接二值化时阈值选太低了比如 0.3把大量模棱两可的像素归为道路同时没有做平滑处理。解决方法是二值化阈值至少试 0.4、0.5、0.6 三个档观察哪个档的连通性最好转矢量之前先对二值图做一次中值滤波或形态学开运算去掉孤立噪点再做矢量化。如果后续要做拓扑分析还需要对矢量结果做「删除小于阈值的面」和「平滑边界」两步操作这在源码包里对应一个postprocess.py脚本。提示调阈值是玄学但有个笨办法最可靠——把预测概率图用灰度图保存下来在 QGIS 里用滑块动态调整显示阈值观察道路连通性和碎斑数量随阈值变化的规律选转折点处的值作为最终阈值。5. 交作业前的加分项IoU 评估与道路中心线提取第五个坑如果避开了你的结果图就已经是「道路区域」而不是「道路碎片」了。但课程设计答辩时老师最关心的往往是两个问题「精度多少怎么算的」和「结果除了图片还能不能导出矢量」。这一章讲怎么答好这两个问题。5.1 IoU 与精度指标计算道路提取这种二分类任务最常被问到的指标是 IoU交并比其次是精确率和召回率。IoU 的计算逻辑是预测道路区域与真实道路区域的交集除以并集取值范围 0 到 1越接近 1 说明预测和标注越吻合。import numpy as np def compute_iou(pred_mask, gt_mask, threshold0.5): 计算预测分割结果与标注之间的 IoU pred_mask: 模型输出的概率图 (H, W)值域 0~1 gt_mask: 标注图 (H, W)值为 0 或 255 threshold: 二值化概率的阈值默认 0.5 pred (pred_mask threshold).astype(np.uint8) gt (gt_mask threshold).astype(np.uint8) intersection np.logical_and(pred, gt).sum() union np.logical_or(pred, gt).sum() return intersection / (union 1e-6)注意gt_mask threshold这一步——如果标注图是 0 和 255 的格式直接和 pred 比较会全都判为背景所以先把标注也二值化。分母加1e-6是为了避免整张图都没有道路时 union 为 0 导致的除零错误。在源码包里我写了一个批量评估脚本它会遍历验证集所有瓦片输出每张图的 IoU、精确率、召回率以及全验证集的平均 IoU。课程设计报告里放一个这样的表格比放三张对比图有说服力得多。5.2 把预测结果变成能交差的东西栅格转矢量与道路中心线答辩时第二个高频问题是「你的结果在 GIS 里能不能用」。模型输出的 PNG 只是栅格图片但道路提取的最终交付物通常是矢量路网。转换的思路分两步先用 OpenCV 把二值栅格转为矢量多边形再用骨架提取得到道路中心线最后把中心线转成折线。骨架提取用 scikit-image 的skeletonize一句就能完成from skimage.morphology import skeletonize import cv2 import numpy as np # 假设 road_mask 已经是 0/255 的二值图 road_bool road_mask 0 skeleton skeletonize(road_bool) # 骨架是布尔数组转成 0/255 便于保存 skeleton_img (skeleton * 255).astype(np.uint8) cv2.imwrite(road_skeleton.png, skeleton_img)skeletonize的输入要求是布尔数组且前景用 1 表示。它会不断剥离边缘像素直到剩下单像素宽的中轴这个过程保留了道路的拓扑连通性。拿到骨架图后可以再用cv2.HoughLinesP把骨架线拟合出折线或者导入 QGIS 做矢量化。有些答辩老师还会追问中心线和道路实际宽度不一致的问题这时你解释清楚「骨架是拓扑抽象宽度信息已经从栅格结果中丢失做缓冲区分析时按平均宽度重建即可」这就体现了你真的理解了整个链路。经过这几个工程版本的迭代我养成了一个习惯不管任务多急预测完结果之后先强制自己走一遍「直方图检查 → 标签膨胀 → 连通域清洗 → IoU 三档阈值复核」这个流程。这条流程帮我挡住了不少翻车现场也让我在答辩时对每个结果图的来龙去脉都心里有数。这份源码包里的预处理脚本、模型训练脚本和后处理脚本都是按这套流程组织的下载后建议从preprocess/目录的裁剪脚本开始跑把数据走通一遍再进模型训练。希望帮到你。本文还有配套的精品资源点击获取