基于U-Net的遥感图像分割与国土分类课程设计实践

简介：面向计算机相关专业课程设计和期末大作业场景，这份源码包实现了基于图像分割的卫星遥感影像国土分类项目。核心代码以 Python 脚本组织，覆盖数据加载、图像预处理、水体预处理以及 DeepLabV3、DeepLabV3+、PSPNet 等主流分割模型，整体结构清晰，便于借鉴和二次开发。压缩包共 12 个文件，以 8 个 .py 源文件为主，同时包含训练日志、数据说明和成果展示图，大小约 2.49MB；其中模型脚本与预处理脚本可直接衔接，结合日志可快速理解训练流程与分类效果。项目已经严格调试并能直接运行，下载后即可作为课程设计或期末大作业的完整参考，尤其适合希望从零搭建遥感图像分类任务、需要完整项目实战练习的计算机专业学生。目前已有 159 人学习，说明该资源在同类作业场景中具备一定参考价值。

1. 卫星遥感国土分类课程设计：它到底解决什么问题

做 python 的基于图像分割对卫星遥感图像进行国土分类项目源码，光看目录会以为难在模型，真跑一遍才知道，模型反而是最好办的部分。你会拿到一张几万乘几万的遥感 tif，想把里面的耕地、林地、建筑、水体、裸地、道路自动分出来，拼成一张完整的分类图——图像分割负责逐像素打标签，卫星遥感和国土分类决定数据怎么切、结果怎么拼。这个课程设计最劝退的点是：数据预处理占一半工作量，另一半是训练和推理的工程细节，模型只占很小一块。我下面按平时做一个完整课程设计项目的顺序讲，切片、掩膜、类别不平衡、训练、推理、避坑，一步步能复现。

2. 遥感影像预处理：窗口切片与掩膜制作是第一步

2.1 先看清楚你的 tif：波段顺序与地理信息

遥感影像不是普通照片，它通常是一个多波段 tif，常见的是 R、G、B 三波段，稍好一点的数据还带近红外波段，一共四个波段。尺寸也离谱，常见一条边就有五千到一万像素。直接整张图读进内存再喂模型，第一步就 OOM。

我拿到数据第一件事是用 rasterio 读元信息，而不是直接 read()。

import rasterio with rasterio.open("scene.tif") as src: print(src.width, src.height, src.count) # 宽、高、波段数 print(src.dtypes) # 每个波段的数据类型 print(src.transform) # 仿射变换参数 print(src.crs) # 坐标系

逻辑说明：width、height 决定切片范围，count 是波段数，dtypes 决定归一化方式。遥感影像很多是 16 位整型，最大值是 65535，不看 dtype 直接除以 255，结果会一片白。之后做掩膜时，transform 和 crs 更关键，切片、栅格化、拼接全要靠这两个参数对齐。

实际项目里我通常先看一眼波段组合，把影像按 RGB 或 R、G、B、NIR 的顺序确认好。不同来源数据波段顺序不一样，比如有的数据是 BGR，有的带 alpha 通道，写代码前先确认，不然训练和推理波段错位，精度会崩得莫名其妙。

2.2 重叠滑窗切片：tile_size 与 stride 怎么定

大图不能直接进 U-Net，常规做法是滑窗切片。tile_size 常见选 256、384、512，课程设计里我一般用 512，信息量足够，模型精度上限高一些，显存也还扛得住。

stride 是滑动步长，等于 tile_size 减去重叠长度。训练切片最好带重叠，因为地物会被切在边缘，比如一栋建筑正好一半在左片一半在右片，模型在边界处判断就很容易出错。推理时也要重叠，后面再讲。

import numpy as np import rasterio from rasterio.windows import Window from pathlib import Path def make_train_tiles(tif_path, out_dir, tile_size=512, overlap=128): stride = tile_size - overlap out_dir = Path(out_dir) out_dir.mkdir(exist_ok=True, parents=True) with rasterio.open(tif_path) as src: n_w = max((src.width - tile_size) // stride + 1, 1) n_h = max((src.height - tile_size) // stride + 1, 1) for row in range(n_h): for col in range(n_w): x = col * stride y = row * stride # 处理右侧和底部补不满的情况，直接贴边取最后一刀 x = min(x, src.width - tile_size) y = min(y, src.height - tile_size) win = Window(x, y, tile_size, tile_size) tile = src.read(window=win).astype(np.float32) np.save(out_dir / f"img_{row}_{col}_{x}_{y}.npy", tile)

逻辑说明：Window(x, y, tile_size, tile_size)表示从影像的(x, y)左上角位置取一块宽高都是tile_size的区域。stride = tile_size - overlap，overlap 设 128 时，512 大小的切片步长就是 384，相邻切片有 128 像素重叠。

这里有个小细节要注意，代码里min(x, src.width - tile_size)是为了防止最后一行或最后一列不够尺寸导致Window越界。代价是最后一次循环会产生一个与上一片重叠很大的切片，这在训练里没有大问题，增强一下就当多样本了。文件名里带上行列坐标，后面训练和拼接时能反推每个切片在原图上的位置。

tile_size 和 stride 的选择可以看下面这个表格：

2.3 矢量标注转掩膜：rasterize 让标注对齐影像

分类标签一般是 shapefile 多边形，比如一块耕地的 polygon、一片水体的 polygon。模型需要的是和影像像素一一对应的掩膜，也就是逐像素的类别标签。这一步用 rasterio.features 的 rasterize。

import geopandas as gpd from rasterio.features import rasterize def vector_to_mask(shp_path, ref_tif_path, out_npy_path, class_field="class_id"): with rasterio.open(ref_tif_path) as ref: shapes = gpd.read_file(shp_path).to_crs(ref.crs) mask = rasterize( [(row.geometry, row[class_field]) for row in shapes.itertuples()], out_shape=(ref.height, ref.width), transform=ref.transform, fill=0, all_touched=True, ) np.save(out_npy_path, mask.astype(np.uint8))

逻辑说明：rasterize 的第二个参数out_shape必须和影像尺寸一致，transform必须用同一张影像的，这样生成的掩膜每个像素都对齐到影像坐标。fill=0表示没有矢量的区域默认是背景类，all_touched=True表示只要多边形碰到的像素都标上，防止细线、窄面被栅格化后断掉。

这个函数里最容易翻车的是坐标参考不一致。shp 如果是经纬度、影像如果是投影坐标，to_crs(ref.crs)会强行转换，但如果 shp 本身坐标系定义错了，转换结果就全偏。所以我会在转换后打印一下 shp 和影像的 bounds 对比，确认多边形范围落在影像范围内。

2.4 类别不平衡与类别权重

国土分类数据集基本逃不开类别不平衡：水体、裸地经常占掉一半以上像素，道路、建筑只占几个百分点。如果直接用普通交叉熵，模型学到最后只会把大类别预测对，小类别几乎全丢。

正确做法是先统计类别像素占比，再算类别权重。

from collections import Counter import numpy as np def compute_class_weights(mask_paths, num_classes): counter = Counter() for path in mask_paths: mask = np.load(path) counter.update(np.unique(mask, return_counts=False).tolist()) total = sum(counter.values()) weights = np.zeros(num_classes, dtype=np.float32) for c in range(num_classes): if counter[c] == 0: weights[c] = 0.0 else: weights[c] = total / (num_classes * counter[c]) return weights

逻辑说明：权重公式是总像素数 / (类别数 * 该类像素数)。某类占比越小，权重越大。这个权重后面传给交叉熵损失。还有一类做法是先取中位数频率再归一化，但对课程设计来说，上面这种 inverse frequency 已经够用。

统计这一步建议在切完片之后做，因为切片的裁剪可能导致类别比例与全图不同。比如全图建筑只占 3%，但在切片里交通道路切片占比较高。我实际项目里都是直接统计所有训练切片，保证权重和训练分布一致。

2.5 Dataset 加载与数据增强

掩膜做完、切片落盘，接下来就是把它们喂给 PyTorch。这里要用 Dataset 把图片切片和掩膜切片配对起来。

import numpy as np import torch from torch.utils.data import Dataset class RSDataset(Dataset): def __init__(self, img_paths, mask_paths, max_value=65535, use_aug=True): self.img_paths = img_paths self.mask_paths = mask_paths self.max_value = max_value self.use_aug = use_aug def __len__(self): return len(self.img_paths) def __getitem__(self, idx): img = np.load(self.img_paths[idx]).astype(np.float32) / self.max_value mask = np.load(self.mask_paths[idx]).astype(np.int64) if self.use_aug: if np.random.rand() > 0.5: img = img[:, :, ::-1] mask = mask[:, ::-1] k = np.random.randint(0, 4) img = np.rot90(img, k, axes=(1, 2)) mask = np.rot90(mask, k) return torch.from_numpy(img), torch.from_numpy(mask)

逻辑说明：max_value取决波段类型，16 位数据是 65535，8 位数据是 255。归一化直接在 Dataset 里做，比训练前预处理好，不然切片落盘时归一化后数据类型和范围对不上，推理时还要再造一遍轮子。增强只做翻转和旋转，因为遥感影像地物方向不定，旋转 90 度不改变语义。水平翻转用[:, :, ::-1]，第一维是波段维，所以只翻转宽那一维。

熟手看到这里可能会问，为什么不用 albumentations？能用，但在课程设计里手写这几行就够了，还能少装一个依赖。想加颜色抖动之类再上 albumentations 也不迟。

3. U-Net 训练：模型结构、损失函数与超参数

3.1 为什么选 U-Net 而不是 DeepLabv3+

课程设计的数据量通常不算大，几千张切片都算多的。DeepLabv3+ 空洞卷积带来的大感受野，在数据不足时反而容易过拟合，训练时间也更长。U-Net 参数少、结构直观、对小样本友好，尤其是跳跃连接可以把浅层的边界细节直接送到解码器，这对分割耕地边界、道路边缘非常关键。

很多人一上来就去找预训练 backbone，其实在这个场景里，从零写一个 U-Net 就够了。我后面做对比实验时，torchvision 里的 DeepLabv3+ 在验证集上反而不如这个轻量 U-Net，因为数据量撑不起那么大的模型。

3.2 U-Net 的国土分类版：输出通道数改成类别数

U-Net 结构本身不复杂，就是下采样编码、上采样解码、跳跃连接。关键改动只有两处：输入通道数改成影像波段数，输出通道数改成分类类别数。

import torch import torch.nn as nn class DoubleConv(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, num_classes, features=(64, 128, 256, 512)): super().__init__() self.downs = nn.ModuleList() self.pool = nn.MaxPool2d(2) for f in features: self.downs.append(DoubleConv(in_ch, f)) in_ch = f self.bottleneck = DoubleConv(features[-1], features[-1] * 2) self.ups = nn.ModuleList() for f in reversed(features): self.ups.append(nn.ConvTranspose2d(f * 2, f, kernel_size=2, stride=2)) self.ups.append(DoubleConv(f * 2, f)) self.head = nn.Conv2d(features[0], num_classes, kernel_size=1) def forward(self, x): skips = [] for down in self.downs: x = down(x) skips.append(x) x = self.pool(x) x = self.bottleneck(x) for i in range(len(self.ups) // 2): x = self.ups[i * 2](x) x = torch.cat([x, skips.pop()], dim=1) x = self.ups[i * 2 + 1](x) return self.head(x)

逻辑说明：features控制每层通道数，课程设计用(64, 128, 256, 512)已经足够。最后self.head用 1x1 卷积把通道数映射到类别数，输出形状是(B, num_classes, H, W)。注意 U-Net 里下采样和上采样的步长都是 2，所以输入尺寸最好是 2 的整数次幂倍数，512 最省心。如果 w 和 h 偏 1 像素，torch.cat会报尺寸不一致，这点排查起来很烦，建议切片时就固定边长。

3.3 损失函数不是只有 CrossEntropy

很多课程设计直接nn.CrossEntropyLoss()就开训练了，然后发现小类目全是错的。原因前面说过，类别不平衡。CrossEntropy 的默认行为是大类主导梯度。

我一般用交叉熵加 DiceLoss 的组合。

def dice_loss(pred, target, smooth=1e-6): pred = torch.softmax(pred, dim=1) target_onehot = nn.functional.one_hot(target, num_classes=pred.shape[1]) target_onehot = target_onehot.permute(0, 3, 1, 2).float() inter = (pred * target_onehot).sum(dim=(2, 3)) union = pred.sum(dim=(2, 3)) + target_onehot.sum(dim=(2, 3)) dice = (2 * inter + smooth) / (union + smooth) return 1 - dice.mean()

逻辑说明：DiceLoss 对类别不平衡不太敏感，它衡量的是预测与标签的重叠度，而不是像素投票平均。one_hot之后要permute成(B, C, H, W)，不然形状对不上。总损失写ce_loss(out, target, weight=class_weights) + dice_loss(out, target)就行。weight 就是 2.4 算出来的类别权重。

这里有个玄学细节：两个 loss 的数值量级最好接近，不然大的那个会主导。常见做法是把 DiceLoss 乘 0.5 或直接把 CE 的 weight 调大。课程设计里先按 1:1 跑，看 loss 曲线如果震荡厉害再调。

3.4 训练超参数：学习率、批大小、epoch 怎么设

U-Net 训练不算重，课程设计用一块普通 GPU 就能跑。我经常用的参数是：

训练循环里我习惯每轮算一次验证 mIoU，然后按 mIoU 保存最优权重。

import torch optimizer = torch.optim.Adam(model.parameters(), lr=3e-4) scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau( optimizer, mode="max", patience=5, factor=0.5 ) best_iou = 0.0 for epoch in range(60): model.train() for imgs, masks in train_loader: imgs, masks = imgs.cuda(), masks.cuda() out = model(imgs) loss = ce_loss(out, masks) + dice_loss(out, masks) optimizer.zero_grad() loss.backward() optimizer.step() val_iou = validate(model, val_loader) print(epoch, loss.item(), val_iou) scheduler.step(val_iou) if val_iou > best_iou: best_iou = val_iou torch.save(model.state_dict(), "best_unet.pt")

逻辑说明：scheduler.step(val_iou)不要传 loss，mIoU 才是最终优化目标。validate函数会在 6.1 讲，但这里注意一点：如果验证集是按整张影像切过来的，mIoU 可以信任；如果是同一切片的随机划分，mIoU 可能虚高，后面第 5 章会专门讲这个坑。

4. 滑窗推理与拼接：把掩膜拼回完整分类图

4.1 滑窗预测：重叠区域取平均置信度

训练完模型，整张大图没办法一次推理，还是要滑窗。推理必须重复合，因为每个切片边缘地物不完整，直接硬切会让最终结果出现一条条拼接痕。常见做法是重叠区域把多次预测的概率平均，再 argmax。

import numpy as np import torch import rasterio from rasterio.windows import Window def predict_large_image(model, tif_path, num_classes, tile_size=512, overlap=128): stride = tile_size - overlap model.eval() with rasterio.open(tif_path) as src: h, w = src.height, src.width prob = np.zeros((num_classes, h, w), dtype=np.float32) count = np.zeros((h, w), dtype=np.uint8) for y in range(0, h - tile_size + 1, stride): for x in range(0, w - tile_size + 1, stride): win = Window(x, y, tile_size, tile_size) img = src.read(window=win).astype(np.float32) / 65535.0 img_t = torch.from_numpy(img).unsqueeze(0).cuda() with torch.no_grad(): p = torch.softmax(model(img_t), dim=1).cpu().numpy()[0] prob[:, y:y + tile_size, x:x + tile_size] += p count[y:y + tile_size, x:x + tile_size] += 1 prob /= np.maximum(count, 1) return np.argmax(prob, axis=0).astype(np.uint8)

逻辑说明：prob是类别数 x 高 x 宽的三维矩阵，每次预测结果累加进去。重叠区域被预测两次或更多次，概率值会更高，最后除以累加次数就得到平均置信度。count矩阵防止除零。这个代码没有处理右侧和底部的最后一条边缘，因为循环条件是range(0, h - tile_size + 1, stride)，非整整除的部分会被漏掉。我一般会在循环跑完后，手动对最右边和最底下的残块再补一次Window推理，逻辑和训练切片时的min修正一样。

4.2 特征图拼回原图：坐标映射与分块写入

推理代码里直接把概率累加到全图上，影像不大没问题。如果影像一条边超过一万像素，prob矩阵会非常大，这时候我会改用np.memmap或者逐行块写概率，避免内存翻车。

拼回来的分类结果要能导回 GIS，这一步必须保持变换参数。

with rasterio.open("result.tif", "w", driver="GTiff", height=h, width=w, count=1, dtype="uint8", transform=src.transform, crs=src.crs) as dst: dst.write(result, 1)

逻辑说明：transform和crs从原图读出来原样写回去，这样分类结果和原影像在 GIS 里完全重合。很多人在这步省事，只写数组不写 transform，结果导进 ArcGIS/QGIS 时位置全偏，这也是课程设计验收时容易被问住的地方。

4.3 分类后处理：去除孤立小块与面积估计

模型输出的掩膜通常有椒盐噪声，尤其在建筑、道路这类小类别上。我一般先做个小窗口模式滤波，再删掉面积过小的连通域。

from scipy.ndimage import median_filter from skimage.morphology import remove_small_objects # 3x3 中值滤波，去掉单点的跳变 result_clean = median_filter(result, size=3) # 按类别分别删掉小于阈值的孤立块 for c in np.unique(result_clean): mask_c = result_clean == c mask_c = remove_small_objects(mask_c, min_size=50) result_clean[mask_c] = c

逻辑说明：remove_small_objects对每个类别二值图单独处理，小于min_size的碎片会被标成 False，再赋回原类别号。注意不要直接对整张多类标签用remove_small_objects，它默认会把背景类当成删除对象，结果反而乱。滤波和碎块删除会让边界稍微平滑，对课程设计来讲是好事，因为最终要统计各类面积比例。

面积统计是需求文档里常有的功能，我一般这样输出：

labels, counts = np.unique(result_clean, return_counts=True) pixel_area = abs(src.transform[0] * src.transform[4]) # 一个像素对应的平方米 for lab, cnt in zip(labels, counts): print(lab, cnt * pixel_area)

逻辑说明：transform[0]是 x 方向分辨率，transform[4]是 y 方向分辨率，乘起来就是单个像素面积。cnt * pixel_area换成平方米，再除一万就是公顷。如果影像没有地理坐标，这一步就省了，只能给像素占比。

5. 避坑排查：五个让遥感分割翻车的细节

5.1 损失函数降得很低，预测结果却一片黑或全是背景

现象：训练 loss 掉到 0.1 以下，验证 mIoU 看着也还行，但推理整张大图时，结果大片大片的背景，只有零星几个类别碎片。原因一般是两种：一类是在切片落盘时 mask 和 img 的文件名排序没对应上，np.load出来的 label 都是同一个背景图；另一类是类别权重里稀有类被设成了 0，模型直接把那个类的梯度掐断了。解决：训练前写一个小脚本，随机抽 8 个 img 和 mask 对，用matplotlib可视化验证。别只 print 文件名，真去看一眼，彩色影像和标签轮廓对不对得上，这一步能省掉后面一整天的排查。

5.2 GPU 显存报错 OOM，连 512 都不行

现象：batch_size 设 8、tile_size 设 512，一跑就 CUDA OOM。原因不一定是显存小，可能是模型里ConvTranspose2d的中间特征图太大，加上概率矩阵也占显存。解决：batch_size 降到 4，如果还不行就把切片改成 384；再不行，把上采样替换成nn.Upsample(scale_factor=2, mode="bilinear", align_corners=False) + Conv2d，这个组合显存占用比转置卷积小 20% 到 30%。我一般先看nvidia-smi确认显存卡在哪个进程，再决定降 batch 还是降 tile_size。注意不要只改模型通道数，那会影响分割精度，课程设计犯不上。

5.3 切片边缘把耕地和道路切成两半，边界区域大片错分

现象：最终分类图有规律的棋盘格状错分，尤其在切片接缝处。原因：训练切片没有重叠，模型没见过跨边缘的上下文，推理时也硬切，边界上下文缺失就判错。解决：训练切片和推理切片都用重叠滑窗。训练时 overlap=128，推理时也设 128，重叠区域概率做平均再 argmax，棋盘格就会淡很多。如果你已经训练完了不想重来，可以在推理时用重叠滑窗，效果也会改善。

5.4 验证集 mIoU 高，换一张影像效果崩

现象：课程设计里常用一张大图切片训练验证，随机分 8:2，验证 mIoU 到了 0.85，自认为稳了，结果拿到另一张不同时间拍摄的影像上一测，mIoU 直接掉到 0.3。原因：随机划分切片的本质是同一张影像上的相邻区域同时进了训练和验证，模型早就见过验证切片的内容，这是典型的数据泄漏。解决：按影像 ID 分 fold，同一张影像的所有切片只能全在训练集或全在验证集。如果手里只有一张影像，那就按空间网格切，比如每隔 4 行取一行作为验证，保证验证切片和训练切片至少间隔一个滑动步长。这个教训我踩过一次后，现在每次分数据都先按影像序列号分。

5.5 掩膜和影像错位半像素，建筑边界飘了

现象：预测的边界和影像里建筑轮廓明显错开，有的地方差三四个像素。原因：rasterize时用了另一张影像的 transform，或者是 shp 与 tif 的坐标系没对齐就硬转。解决：先检查 shp 的 bounds 和影像的 bounds 是否大体一致，再gpd.read_file(...).to_crs(src.crs)后打印一遍坐标范围。最后用透明度图层把掩膜叠在影像上肉眼看一眼，轮廓贴合再进入训练。遥感这个过程不像自然图像能靠直觉发现问题，不核对坐标系，后面所有精度指标都是自欺欺人。

6. 验证技巧：mIoU 计算与训练集划分的进阶做法

最终验收不能只靠肉眼，mIoU 是课程设计里最常被问的指标。我一般这样算单个批次或整张图的 mIoU：

import numpy as np def compute_miou(pred, target, num_classes): ious = [] for c in range(num_classes): p = pred == c t = target == c inter = np.logical_and(p, t).sum() union = np.logical_or(p, t).sum() if union == 0: continue ious.append(inter / union) return np.mean(ious)

逻辑说明：mIoU 是每个类别交并比的平均值。union == 0说明这个类别在样本里不存在，跳过它，不要算成 0，那样会把分数拖低。如果你想严格一点，也可以在报告里注明“未出现类别不计入”，这是学术上比较常见的做法。

进阶做法是测试时增强，推理时对每张切片做水平翻转，把两次 softmax 概率平均后再 argmax。这个改动能让 mIoU 稳定提升 1 到 2 个百分点。代价是推理时间翻倍，但课程设计场景无所谓。数据划分上我现在的习惯是先把所有切片按来源影像分组，再在组级别做 train、val、test 划分，宁可让训练集少一点，也要保证测试集是模型没见过的区域。

这个项目做完后，我最深的感受是遥感分割的坑几乎全在数据层。模型代码半天就能写完，切片、掩膜、坐标对齐、数据划分这些环节才是真正花时间的地方。现在我做这类任务，会先花一个晚上把数据链路可视化确认清楚，再启动训练，而不是直接写模型。希望帮到你。

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