☰
Python遥感图像道路提取:从U-Net到IoU评估的完整实践
2026/9/27 23:26:20 网站建设 项目流程

简介:这是一套基于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 格式,部分数据还附带有理多项式校正参数,需要在预处理阶段处理。

数据源全色分辨率多光谱分辨率波段数适用场景
高分一号 PMS2m8m3(R/G/B)大范围路网拓扑提取
高分二号 PMS0.8m3.2m3(R/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(window=window)) # 示例:从第 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_percent=2): """将多波段数组转为 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), dtype=np.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 路线 A:Canny 边缘检测 + 形态学闭运算做快速原型

传统路线我会用 OpenCV 完成,核心思路是边缘检测后用形态学闭运算把断裂的道路边缘连接起来,再用面积筛选去掉碎斑。这个流程对清晰的城市主干道有效,对乡村土路和阴影密集区效果有限,但作为 baseline 足够。

import cv2 import numpy as np # 读取上一章生成的 8 位 RGB 瓦片 img = cv2.imread("tile_8bit.png", cv2.IMREAD_GRAYSCALE) # 1. 高斯模糊抑制噪声,核大小取 5x5,sigma 默认 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, iterations=2) # 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 路线 B:U-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, padding=1), nn.BatchNorm2d(out_ch), nn.ReLU(inplace=True), nn.Conv2d(out_ch, out_ch, 3, padding=1), nn.BatchNorm2d(out_ch), nn.ReLU(inplace=True)) def forward(self, x): return self.conv(x) class UNet(nn.Module): def __init__(self, in_ch=3, out_ch=1): 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, stride=2) self.dec2 = ConvBlock(256, 128) self.up1 = nn.ConvTranspose2d(128, 64, 2, stride=2) 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], dim=1)) d1 = self.dec1(torch.cat([self.up1(d2), e1], dim=1)) return self.out(d1)

这个模型的关键设计在于跳跃连接——上采样后的特征图与同层编码器的输出在通道维度上拼接(torch.cat),让解码器同时拥有高层语义信息和低层边缘纹理。注意拼接后通道数翻倍,所以dec2和dec1的输入通道数分别是 256 和 128。ConvBlock里的BatchNorm2d在遥感影像这种分布波动大的输入上很重要,它能稳定训练过程,没有它的话模型可能需要三倍以上的 epoch 才能收敛。

训练时的参数选择直接影响结果好坏,我一般按下面的配置起步:

criterion = nn.BCEWithLogitsLoss() optimizer = torch.optim.Adam(model.parameters(), lr=1e-4) scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=10, gamma=0.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, threshold=0.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/目录的裁剪脚本开始跑,把数据走通一遍再进模型训练。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询