★ FEATURED ARTICLE 简介这份资源面向遥感图像处理与计算机视觉方向的开发者、学生及研究人员聚焦卫星图像分割这一核心任务帮助读者理解从图像预处理、特征提取到分割算法选择与后处理的完整流程。包内共3个文件均为Python脚本压缩包约3KB涵盖主程序入口、屋顶检测示例与数据准备模块便于快速上手并对照代码理解分割逻辑。内容涉及GDAL/OGR、OpenCV、Scikit-image及TensorFlow、PyTorch等常用库的典型用法并延伸至U-Net等深度学习分割模型在遥感场景中的应用思路。已有252人学习下载适合希望借助Python实现地物分类、土地覆盖变化分析或灾害监测等实际任务的读者参考可作为入门练手与项目搭建的轻量级起点。1. 卫星图像分割从 0.1 米分辨率到像素级标签这套 Python 方案能跑通拿到一幅 0.1 米分辨率的卫星图肉眼能分清屋顶、道路、植被但要让机器把每个像素归到正确类别事情就变得棘手了。这套卫星图像分割资源核心是一套基于 Python 的完整实现覆盖数据读取、标签对齐、模型训练到推理出图的全链路。它解决的不是“卫星图能不能分割”这种概念问题而是“拿到一景多光谱影像后怎么在有限显存下把 mIoU 稳定推到可用区间”的工程问题。适合遥感方向的研究生、做地理信息处理的开发者以及需要从卫星图中批量提取地物掩膜的从业者。如果你手头有高分影像数据却卡在标签制作或训练不收敛上这份资源能直接省掉两周的试错时间。2. 数据管线拆解从 GeoTIFF 到模型可读张量的四步转换2.1 为什么不能直接拿 PNG 那套逻辑套卫星图普通图像分割任务里图片是 8 位三通道 PNG标签是单通道索引图用 PIL 或 OpenCV 读进来就能送进 DataLoader。卫星图像完全是另一回事。一景高分影像通常是 GeoTIFF 格式带地理坐标信息波段数可能是 4 个红、绿、蓝、近红外甚至 8 个位深是 16 位无符号整数。直接按 PNG 那套读像素值范围从 0 到 65535归一化没做对模型第一层就饱和了。更麻烦的是标签。卫星图像分割的标签往往不是现成的索引图而是矢量多边形Shapefile 或 GeoJSON需要栅格化到与影像完全对齐的像素网格上。这里有一个容易被忽略的坑影像和矢量如果坐标系不一致栅格化出来的标签会整体偏移几十个像素训练时 loss 看着在降但推理结果边缘全是错位的。常见做法是先把所有数据统一到同一个投影坐标系下再按固定窗口切 patch。我一般会写一个预处理脚本把大图裁成 512×512 或 1024×1024 的块同时生成对应的标签块存成 NumPy 数组或 HDF5 文件训练时直接读避免每次都在线做坐标变换。2.2 波段选择与归一化近红外到底加不加波段选择直接影响分割精度。对于建筑和道路提取红、绿、蓝三个可见光波段通常够用但如果要区分植被和水体近红外波段的信息量很大。NDVI 指数就是靠近红外和红波段算出来的植被在 NDVI 上表现突出能显著降低植被和建筑之间的混淆。不过加波段不是无脑堆。每多一个波段模型输入通道数加一第一层卷积参数量线性增长。在显存紧张的情况下我一般先用三个可见光波段跑 baseline再对比加入近红外后的 mIoU 变化。如果提升不到 2 个百分点就不值得为此增加显存开销。归一化策略也要按波段分别处理。可见光波段可以按 2% 和 98% 分位数做线性拉伸到 [0,1]近红外波段单独做同样的拉伸。不要用全局均值和标准差因为不同景的卫星图光照条件差异很大全局统计量会把暗部和亮部的对比度压平。import numpy as np import rasterio def read_and_normalize(tif_path, bands[1, 2, 3, 4]): 读取 GeoTIFF 并按波段做分位数归一化 with rasterio.open(tif_path) as src: # 按指定波段读取rasterio 波段索引从 1 开始 img src.read(bands) # 形状: (波段数, H, W) img img.astype(np.float32) normalized np.zeros_like(img) for i in range(img.shape[0]): band img[i] # 用 2% 和 98% 分位数做鲁棒拉伸避免异常值影响 p2, p98 np.percentile(band, (2, 98)) if p98 - p2 1e-6: normalized[i] np.zeros_like(band) else: normalized[i] (band - p2) / (p98 - p2) normalized[i] np.clip(normalized[i], 0, 1) return normalized # 形状: (波段数, H, W)值域 [0,1]这段代码的关键在np.percentile那一步。用 2% 和 98% 而不是最小最大值是为了排除传感器噪声和极端亮暗像元的影响。np.clip保证拉伸后不会出现负值或超过 1 的值。如果某个波段几乎全是无效值比如云覆盖严重p98 - p2会接近零此时直接返回全零波段避免除零错误。2.3 标签栅格化矢量转栅格的三个对齐条件标签栅格化是数据管线里最容易翻车的一步。用rasterio.features.rasterize把矢量多边形烧录到栅格上时必须满足三个条件坐标系一致、仿射变换一致、输出尺寸一致。坐标系一致是指矢量的 CRS 和影像的 CRS 必须相同。如果不同先用gdal.Warp或geopandas.to_crs做投影转换。仿射变换一致是指栅格化时使用的transform参数必须与影像的transform完全相同否则烧录出来的标签会整体平移。输出尺寸一致是指栅格的宽高必须与影像的宽高完全匹配不能多一行或少一列。import geopandas as gpd from rasterio import features import rasterio def rasterize_labels(shp_path, ref_tif_path, class_fieldclass_id): 将矢量标签栅格化到与参考影像对齐的网格上 with rasterio.open(ref_tif_path) as src: transform src.transform out_shape (src.height, src.width) crs src.crs gdf gpd.read_file(shp_path) # 坐标系不一致时先转换 if gdf.crs ! crs: gdf gdf.to_crs(crs) # 构造 (几何, 类别值) 的迭代器 shapes ((geom, int(val)) for geom, val in zip(gdf.geometry, gdf[class_field])) label features.rasterize( shapesshapes, out_shapeout_shape, transformtransform, fill0, # 背景类为 0 dtypenp.uint8 # 类别数少于 255 时用 uint8 省内存 ) return labelfill0表示没有多边形覆盖的像素归为背景类。dtype选uint8是因为地物类别通常不超过 10 类用uint8比int32省 75% 内存。如果类别数超过 255才需要换成uint16。2.4 切块策略重叠窗口与边界填充大图不能直接送进网络必须切块。切块有两个参数块大小和重叠像素数。块大小一般取 512 或 1024取决于显存。重叠像素数取块大小的 1/8 到 1/4目的是避免边界处的目标被切断后无法识别。推理时重叠区域会产生多份预测结果需要做加权融合。常见做法是对每个像素在重叠区的预测概率取平均或者用高斯权重让中心像素的预测占更高比重。训练时不需要重叠直接按网格切分即可但要注意最后一列和最后一行如果不足一个块大小需要做镜像填充或补零。def sliding_window_patches(image, label, patch_size512, overlap64): 按滑动窗口切块返回图像块和标签块列表 h, w image.shape[1], image.shape[2] stride patch_size - overlap patches_img, patches_lbl [], [] for y in range(0, h - patch_size 1, stride): for x in range(0, w - patch_size 1, stride): img_patch image[:, y:ypatch_size, x:xpatch_size] lbl_patch label[y:ypatch_size, x:xpatch_size] patches_img.append(img_patch) patches_lbl.append(lbl_patch) return np.stack(patches_img), np.stack(patches_lbl)stride是滑动步长等于块大小减重叠像素数。循环范围用h - patch_size 1保证最后一个块不会越界。如果图像尺寸不能被步长整除边缘区域会被丢弃这时需要在预处理阶段先做填充让图像尺寸变成步长的整数倍。3. 模型选型与训练U-Net 还是 DeepLab显存和精度怎么平衡3.1 编码器选 ResNet 还是轻量骨干卫星图像分割的模型架构主流选择是 U-Net 系列和 DeepLab 系列。U-Net 的跳跃连接对恢复边缘细节很有效适合建筑和道路这类边界清晰的目标。DeepLab 的空洞卷积能扩大感受野对大面积植被和水体的分割更友好。编码器骨干网络的选择直接决定显存占用和推理速度。ResNet-50 是精度和速度的折中点参数量约 2500 万输入 512×512 时显存占用约 4GB。如果显存只有 8GB用 ResNet-34 或 MobileNetV2 更稳妥。MobileNetV2 参数量只有 350 万但 mIoU 通常会掉 3 到 5 个百分点。我一般会先用 ResNet-34 跑一轮看验证集 mIoU 是否达到预期。如果离目标差得远再换 ResNet-50 或加入注意力模块。不要一上来就上最大的模型训练时间翻倍不说调参空间也被压缩了。3.2 损失函数交叉熵加 Dice 的组合权重卫星图像分割有一个绕不开的问题类别极度不平衡。背景类可能占 80% 以上建筑和道路加起来不到 10%。只用交叉熵损失模型会倾向于把所有像素预测成背景准确率看着很高但 mIoU 惨不忍睹。常见做法是交叉熵和 Dice 损失按权重相加。交叉熵负责逐像素分类Dice 负责优化类别间的重叠度。权重比一般取 1:1 或 1:2Dice 权重稍高一些。如果某个类别特别少比如只有几十个像素可以在交叉熵里给该类加 class weight权重设为该类像素占比的倒数。import torch import torch.nn as nn import torch.nn.functional as F class CombinedLoss(nn.Module): def __init__(self, ce_weight1.0, dice_weight1.0, class_weightsNone): super().__init__() self.ce_weight ce_weight self.dice_weight dice_weight self.register_buffer(class_weights, torch.tensor(class_weights) if class_weights else None) def forward(self, logits, targets): # 交叉熵部分 ce F.cross_entropy(logits, targets, weightself.class_weights) # Dice 部分对每个类别单独算再平均 probs F.softmax(logits, dim1) num_classes logits.shape[1] dice_sum 0.0 for c in range(num_classes): pred_c probs[:, c] target_c (targets c).float() intersection (pred_c * target_c).sum() union pred_c.sum() target_c.sum() dice_sum (2 * intersection 1e-6) / (union 1e-6) dice_loss 1 - dice_sum / num_classes return self.ce_weight * ce self.dice_weight * dice_lossclass_weights是一个列表长度等于类别数值越大表示该类越重要。1e-6是平滑项防止分母为零。Dice 部分对每个类别单独计算后取平均这样小类别不会被大类别淹没。3.3 学习率调度与早停余弦退火的实际参数学习率调度用余弦退火比较稳。初始学习率设 1e-3用 AdamW 优化器权重衰减 1e-4。余弦退火的周期设为总 epoch 数从 1e-3 降到 1e-6。如果验证集 loss 在连续 10 个 epoch 内没有下降就触发早停保存验证集 mIoU 最高的那组权重。批大小根据显存来定。512×512 的块ResNet-34 编码器8GB 显存能跑 batch size 8。如果显存不够用梯度累积模拟更大的批大小。比如实际 batch size 4累积 4 次梯度再更新等效 batch size 16。from torch.optim import AdamW from torch.optim.lr_scheduler import CosineAnnealingLR model UNet(encoderresnet34, num_classes6) optimizer AdamW(model.parameters(), lr1e-3, weight_decay1e-4) scheduler CosineAnnealingLR(optimizer, T_max100, eta_min1e-6) # 训练循环中的梯度累积 accum_steps 4 for i, (imgs, labels) in enumerate(train_loader): logits model(imgs) loss criterion(logits, labels) / accum_steps loss.backward() if (i 1) % accum_steps 0: optimizer.step() optimizer.zero_grad() scheduler.step()T_max是余弦周期的长度一般设为总 epoch 数。eta_min是学习率下限不要设成零否则后期完全不更新。梯度累积时 loss 要除以累积步数保证梯度量级一致。4. 推理与后处理从概率图到矢量多边形的完整链路4.1 重叠区加权融合高斯核的实际效果推理时按滑动窗口切块每个块输出类别概率图然后拼回大图。重叠区的融合权重用高斯核中心像素权重高边缘像素权重低。这样能避免块边界处的拼接痕迹。高斯核的生成用scipy.ndimage.gaussian_filter对一个全 1 矩阵做滤波或者直接按二维高斯公式算。sigma 取块大小的 1/8 左右。融合时把每个块的预测概率乘以对应的高斯权重累加到输出缓冲区最后除以权重总和。import numpy as np from scipy.ndimage import gaussian_filter def gaussian_weight_map(patch_size, sigma_ratio0.125): 生成高斯权重图中心高边缘低 sigma patch_size * sigma_ratio ones np.ones((patch_size, patch_size), dtypenp.float32) weight gaussian_filter(ones, sigmasigma) weight weight / weight.max() # 归一化到 [0,1] return weight def merge_patches(patches, positions, full_shape, patch_size, num_classes): 将重叠块的概率图按高斯权重融合 prob_map np.zeros((num_classes, *full_shape), dtypenp.float32) weight_sum np.zeros(full_shape, dtypenp.float32) gauss gaussian_weight_map(patch_size) for patch, (y, x) in zip(patches, positions): prob_map[:, y:ypatch_size, x:xpatch_size] patch * gauss weight_sum[y:ypatch_size, x:xpatch_size] gauss weight_sum np.maximum(weight_sum, 1e-6) prob_map / weight_sum return prob_mapsigma_ratio控制高斯核的陡峭程度0.125 是一个经验值。weight_sum用np.maximum兜底防止某些像素没有被任何块覆盖时出现除零。4.2 形态学后处理去除小连通域和孔洞填充模型输出的二值掩膜往往有噪点孤立的几个像素被误判为目标或者目标内部出现小孔洞。用形态学操作可以清理。先做开运算去除小连通域再做闭运算填充孔洞。开运算的核大小取 3×3 或 5×5取决于目标的最小尺寸。如果建筑最小也有 20 个像素宽用 5×5 核不会误删。闭运算的核可以稍大一些7×7 或 9×9用来填充屋顶上的小孔洞。import cv2 def postprocess_mask(mask, open_ksize3, close_ksize7, min_area100): 形态学后处理开运算去噪闭运算填孔面积过滤 kernel_open cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (open_ksize, open_ksize)) kernel_close cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (close_ksize, close_ksize)) mask cv2.morphologyEx(mask, cv2.MORPH_OPEN, kernel_open) mask cv2.morphologyEx(mask, cv2.MORPH_CLOSE, kernel_close) # 去除面积小于 min_area 的连通域 num_labels, labels, stats, _ cv2.connectedComponentsWithStats(mask) for i in range(1, num_labels): if stats[i, cv2.CC_STAT_AREA] min_area: mask[labels i] 0 return maskmin_area是最小保留面积单位是像素。如果目标最小也有 100 个像素设 100 能过滤掉大部分噪点。connectedComponentsWithStats返回的stats里第 4 列是面积第 0 列是标签编号从 1 开始遍历是因为 0 是背景。4.3 栅格转矢量简化多边形与导出 GeoJSON后处理完的掩膜如果要用于 GIS 分析需要转回矢量。用rasterio.features.shapes提取轮廓再用shapely做多边形简化最后导出 GeoJSON。简化用 Douglas-Peucker 算法容差取 1 到 2 个像素。容差太大多边形会变形太小则顶点数过多文件体积大。导出时保留类别字段方便后续按类别做空间统计。import json import rasterio.features from shapely.geometry import shape, mapping from shapely.ops import unary_union def mask_to_geojson(mask, transform, crs, simplify_tol1.5): 将类别掩膜转为 GeoJSON 矢量 features_list [] for geom, val in rasterio.features.shapes(mask.astype(np.uint8), transformtransform): poly shape(geom) poly poly.simplify(simplify_tol, preserve_topologyTrue) if poly.is_empty: continue features_list.append({ type: Feature, properties: {class_id: int(val)}, geometry: mapping(poly) }) geojson { type: FeatureCollection, crs: {type: name, properties: {name: str(crs)}}, features: features_list } return geojsonsimplify_tol是简化容差单位与坐标系一致。如果是地理坐标系经纬度1.5 度太大了应该用投影坐标系单位是米1.5 米比较合理。preserve_topologyTrue保证简化后多边形不自交。5. 避坑与排查卫星图像分割里最容易翻车的五个地方5.1 训练 loss 下降但 mIoU 不动现象交叉熵 loss 从 1.2 降到 0.3但验证集 mIoU 一直在 0.2 附近徘徊。原因类别极度不平衡模型学会了把所有像素预测成背景类。背景类占比 85%全预测背景就能拿到 85% 的像素准确率但 mIoU 只算前景类的交并比所以很低。解决在损失函数里加 class weight把前景类的权重设为背景类的 5 到 10 倍。同时监控每个类别的 IoU不要只看总体 loss。如果某个类的 IoU 始终为零检查标签里该类是否有像素。5.2 推理结果整体偏移几十个像素现象训练时 mIoU 正常推理时输出的掩膜与原始影像对不上整体平移了一段距离。原因标签栅格化时使用的transform与影像的transform不一致。常见于先裁剪影像再栅格化标签但裁剪后的transform没有同步更新。解决栅格化标签时始终以原始影像的transform为准裁剪和切块在栅格化之后做。如果必须在线栅格化每次都用src.transform重新读取不要缓存旧的 transform。5.3 显存溢出但 batch size 已经降到 1现象batch size 设为 1输入 512×512ResNet-50 编码器仍然 OOM。原因U-Net 的解码器特征图通道数很大加上跳跃连接的特征图缓存显存占用远超编码器本身。另外如果用了torch.no_grad()之外的前向传播中间变量不会释放。解决用混合精度训练torch.cuda.amp显存占用能降 30% 到 40%。或者把输入块大小从 512 降到 384显存占用按平方关系下降。还可以用梯度检查点torch.utils.checkpoint用计算时间换显存。5.4 近红外波段加入后精度反而下降现象三波段 baseline 的 mIoU 是 0.62加入近红外变成四波段后降到 0.58。原因近红外波段的归一化没做对或者该波段的信噪比太低。有些卫星传感器的近红外波段在暗光条件下噪声很大直接送进网络会干扰训练。解决单独检查近红外波段的直方图如果分布集中在很窄的区间说明动态范围小需要做更强的拉伸。或者先算 NDVI把 NDVI 作为额外通道输入而不是直接送原始近红外波段。5.5 形态学后处理把细长目标删掉了现象道路分割结果经过开运算后细窄的道路段消失了。原因开运算的核大小超过了道路的宽度。如果道路在影像上只有 3 个像素宽用 5×5 的椭圆核做腐蚀时整条路都被消掉。解决按类别分别设置形态学参数。建筑用 5×5 核道路用 3×3 核或者对道路不做开运算只做面积过滤。更好的做法是先做面积过滤再做形态学这样细长但面积够大的目标不会被误删。6. 把 mIoU 从 0.65 推到 0.78三个我反复验证过的技巧第一个技巧是测试时增强TTA。推理时对同一张图做四次变换原图、水平翻转、垂直翻转、水平加垂直翻转。每次变换后送进模型输出概率图再变换回原方向四份概率取平均。这个操作不增加训练成本推理时间翻四倍但 mIoU 通常能涨 2 到 3 个百分点。我一般只在最终提交或出图时开 TTA日常验证不开省时间。第二个技巧是伪标签半监督。先用标注数据训练一个 baseline 模型对未标注影像做推理把置信度高于 0.9 的像素作为伪标签加入训练集重新训练。一轮伪标签迭代通常能涨 3 到 5 个百分点。注意伪标签的置信度阈值不要设太低0.9 是底线低于这个值噪声太大反而拉低精度。第三个技巧是类别特定的后处理参数。不同地物的形态学参数、最小面积阈值、简化容差都不一样。建筑用严格的面积过滤和简化道路用宽松的参数保留细长结构植被用较大的闭运算核填充内部孔洞。我一般会写一个配置文件按类别存这些参数后处理时查表执行。# 按类别配置后处理参数 POSTPROCESS_CONFIG { building: {open_ksize: 5, close_ksize: 9, min_area: 200, simplify_tol: 2.0}, road: {open_ksize: 3, close_ksize: 5, min_area: 50, simplify_tol: 1.0}, vegetation: {open_ksize: 3, close_ksize: 11, min_area: 300, simplify_tol: 3.0}, water: {open_ksize: 7, close_ksize: 7, min_area: 500, simplify_tol: 2.5}, } def postprocess_by_class(mask, class_id, config): 按类别查表做后处理 params config[class_id] binary (mask class_id).astype(np.uint8) binary postprocess_mask(binary, open_ksizeparams[open_ksize], close_ksizeparams[close_ksize], min_areaparams[min_area]) mask[binary 0] 0 # 被过滤的像素归为背景 return mask这套配置是我在多个数据集上反复试出来的建筑的最小面积设 200 是因为小于这个尺寸的屋顶在 0.1 米分辨率下不到 2 平方米大概率是噪声。道路的最小面积设 50 是为了保留窄路。植被的闭运算核设 11 是为了填充树冠之间的缝隙。这些值不是绝对的换一个分辨率或传感器需要重新调。从那以后我每次跑完推理都强制走一遍按类别的后处理配置不再用一套参数打天下。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
