简介面向遥感与深度学习交叉方向学习者的CNN地物分类实战资源对应Landsat影像的像素级分类任务适合计算机、地信、人工智能等专业学生用于课程设计或毕业设计。压缩包共10个文件14.89MB包含Python源码、训练好的模型权重h5、Landsat示例影像tif及配套地理配准文件tfw/xml以及项目说明文档。源码按流程拆分影像切片生成、模型训练、新数据预测三个脚本可直接对照README运行从数据准备到结果输出形成完整闭环。已有969人学习下载。通过该资源可掌握基于CNN的遥感影像地物分类基本流程了解训练数据制作、模型保存与加载、未知影像预测等关键环节也可作为进一步开展植被/水体/建筑等地物识别实验的起点。1. 遥感地物分类的现状为什么用 CNN 成了必然选择做 Landsat 影像地物分类很多人第一反应还是最大似然、随机森林这类传统分类器。遇到 30 米分辨率、7 个有效波段的 OLI 数据植被和农田、裸土和建筑、水体和阴影之间光谱重叠严重传统方法调参调到头也只能做到 85% 左右的总体精度。换成基于 cnn 卷积神经网络的深度学习方案之后同样的训练样本总体精度普遍能拉到 93% 以上水体和阴影这种常年翻车的类别也有了明显改善。这套「基于CNN深度学习的遥感landsat影像地物分类方法」源码包解决的就是从 Landsat 原始影像到最终分类成图这一整条链路的问题多波段数据怎么组织、标签怎么对齐、模型怎么训练、预测结果怎么拼回完整影像。适合手里有 Landsat 数据和一套标签、想把分类精度再往上提一提的从业者和研究生也适合刚接触遥感深度学习的新手把它作为第一个能完整跑通的项目。值得注意的是CNN 解决的远不只是「换个更好的分类器」这么简单它把特征提取和分类决策合并在了一个模型里这才让那些传统方法完全分不开的地物成为可能。2. 从影像到样本Landsat 数据预处理与训练切片2.1 波段选择不是所有波段都该进模型拿到 Landsat 8 OLI 的 Level-1 或者 Level-2 数据常见的是一个包含 11 个波段的 GeoTIFF或者分波段存储的一堆 TIF 文件。很多第一次做的人会想「波段越多信息越丰富」直接全波段往里塞。实际上 B1 海岸带波段气溶胶散射影响大B9 卷云波段主要是大气校正用的B10、B11 热红外是 100 米重采样到 30 米的空间细节本来就丢了。硬塞进去模型容量被无效信息占掉训练还变慢。我一般保留 B2蓝、B3绿、B4红、B5近红外、B6短波红外 1、B7短波红外 2这 6 个 30 米波段。B5 对植被叶绿素含量极其敏感B6/B7 对土壤湿度、建筑物材质区分度很高B2-B4 提供真彩色基准。下面是读取和筛波段的代码。import rasterio import numpy as np def load_landsat_bands(tif_path, bands[2, 3, 4, 5, 6, 7]): 读取 Landsat 8 OLI 多波段影像band 编号对应文件内的波段序号。 注意rasterio 读取时下标从 1 开始与 numpy 下标从 0 开始不同。 with rasterio.open(tif_path) as src: # 读取指定波段src.read 传入的是波段序号列表返回 shape (波段数, 高, 宽) data src.read(bands) profile src.profile transform src.transform # 转成 (高, 宽, 波段数) 的排列方便后续按像素窗口切片 data np.transpose(data, (1, 2, 0)) return data, profile, transform # 使用示例 image, profile, transform load_landsat_bands(LC08_L2SP_137039_20230501_20230509_02_T1_SR.TIF) print(image.shape) # 期望输出 (h, w, 6) print(image.dtype) # 期望输出 uint16Landsat SR 产品通常是 16 位整数这里有两处关键点第一src.read()接受的是波段序号列表开发者容易混淆的是 rasterio 的波段从 1 开始而 numpy 数组下标从 0 开始写代码时不要搞混。第二Landsat Collection 2 的 Surface ReflectanceSR产品虽然名义上是反射率但存储的是放大 10000 倍的整型dtype 是uint16直接丢进模型之前必须做归一化这点在后面 2.3 会具体处理。2.2 影像裁剪与归一化训练样本到底怎么切Landsat 单景影像大约是 7800×7800 像素GPU 显存根本放不下整景。常见做法是把影像和标签一起切成固定大小的 patch比如 128×128 或者 256×256。patch 太小模型看不到足够的空间上下文建筑物阴影和水体更难区分patch 太大类别占比失衡严重一个小村庄可能只占 20×20 像素其余全是农田。我做过对比128×128 在 30 米分辨率下对应 3.84 公里见方既能覆盖中等尺度的地物纹理又不会让单个类别占比极端化。归一化不能对整景影像用(x - min) / (max - min)因为单景影像里云、雪、亮色建筑会产生极端值把水体这种低值区域压到几乎为 0模型对水体基本学不动。更稳的做法是按百分位截断比如把每景影像的 2% 和 98% 分位当作最小值、最大值做线性拉伸。def normalize_percentile(image, lower2, upper98): 按分位数截断归一化逐波段独立计算。 输入 image shape (h, w, bands)dtype 任意数值型。 返回 float32范围 [0, 1]。 h, w, bands image.shape normalized np.zeros_like(image, dtypenp.float32) for b in range(bands): band_data image[:, :, b] p_low np.percentile(band_data, lower) p_high np.percentile(band_data, upper) # 防止影像中某波段是常数导致除零 if p_high - p_low 1e-6: normalized[:, :, b] 0.0 continue band_norm (band_data - p_low) / (p_high - p_low) normalized[:, :, b] np.clip(band_norm, 0, 1) return normalized.astype(np.float32)归一化时容易出现一个隐藏陷阱如果训练数据来自多景影像比如夏天一景、冬天一景那么分位数应该每景各自计算而不是全放一起算。不同季节的 Landsat 影像辐射差异本来就大放一起归一化等于把所有影像拉到同一个亮度基准反而破坏了季节差异带来的分类线索。我自己会在每景影像单独归一化之后再做一次全体均值和方差的统计确认各景之间分布没有特别离谱的偏移。2.3 标签生成与数据集划分空间自相关是最大隐患标签通常有两种来源一种是已有的矢量分类图比如土地利用调查的 shapefile需要栅格化成与影像同一个 grid另一种是人工目视解译画出来的多边形同样要栅格化。栅格化时最容易出问题的是像元对齐——矢量栅格化用的 transform 和影像的 transform 不一致导致标签整体偏移几个像素。最稳妥的方式是直接用rasterio.features.rasterize并且传入影像本身的 transform 和 shape。import rasterio.features import geopandas as gpd def vector_to_label(shp_path, profile, transform, class_fieldclass_id): 把矢量标签栅格化为分类标签图。 profile 和 transform 取自对应影像保证标签与影像严格对齐。 gdf gpd.read_file(shp_path) # 确保矢量与影像坐标系一致不一致时先重投影 if gdf.crs.to_string() ! profile[crs].to_string(): gdf gdf.to_crs(profile[crs]) shapes [(geom, int(value)) for geom, value in zip(gdf.geometry, gdf[class_field])] label_array rasterio.features.rasterize( shapes, out_shape(profile[height], profile[width]), transformtransform, fill0, # 0 作为背景类不参与训练 dtypenp.uint8 ) return label_array栅格化之后滑窗切片时必须保证影像 patch 和标签 patch 在空间位置一一对应最简单的方式是对影像和标签用同一个窗口起始坐标做切片不要单独生成各自的窗口。数据划分上踩过的大坑是真随机打乱像素。把整景影像的所有 patch 混在一起随机分训练集、验证集验证集精度能到 97%但换一景新影像直接掉到 70%。原因在于相邻 patch 之间有大量重叠像素或者同一块农田被切进两个集合模型「背」下了空间位置而非泛化能力。正确做法是按空间位置划分比如把影像从中间劈成两半一半做训练、一半做验证或者干脆用另一景完全不相邻的影像做验证。def split_patches_by_region(images, labels, val_ratio0.2): 按空间区域划分数据避免相邻 patch 被分到 train 和 val 两侧。 这里简单按行方向分块更好的做法是按地理坐标或行政区划划分。 n_patches len(images) split_index int(n_patches * (1 - val_ratio)) train_images, train_labels images[:split_index], labels[:split_index] val_images, val_labels images[split_index:], labels[split_index:] return (train_images, train_labels), (val_images, val_labels)这里的split_by_region看起来简单实际使用时必须保证传入的 patch 列表本身就是按空间顺序排列的切片生成顺序如果被打乱过这个函数就失效了。另一个更可靠的做法是在滑窗生成 patch 时记录每个 patch 的左上角像素坐标(row, col)然后按坐标范围而不是索引来做划分。3. 源码里的 CNN 模型核心结构与参数这样调3.1 模型骨架轻量 U-Net 比纯分类网络更适合地物制图Landsat 地物分类的输出是逐像素的地物类别所以这里不能用 ImageNet 那种输出单一类别的分类网络要做语义分割。最常见的骨干结构是 U-Net 及其变体它的编码器部分逐层下采样提取高层语义特征解码器部分通过跳连接把低层纹理细节融合回来最后由 1×1 卷积输出每个像素的类别概率。这个对称的编码器-解码器结构在遥感影像分割里几乎是默认选择。源码包里最常见的实现是简化的 U-Net编码器三层、解码器三层初始通道 32 或 64。通道数设太大会让参数量爆炸Landsat 才 6 个输入波段不需要 ResNet50 级别的容量设太小又学不动地物的纹理差异。30 米分辨率的地物分类32 起步、逐层翻倍到 128 就够用了。import torch import torch.nn as nn class SimpleUNet(nn.Module): 轻量 U-Net面向 Landsat 6 波段输入。 in_channels6out_channels类别数含背景就多一维。 def __init__(self, in_channels6, num_classes6, base_channels32): super().__init__() # 编码器两个卷积块 一次池化 self.enc1 self._block(in_channels, base_channels) self.enc2 self._block(base_channels, base_channels * 2) # 解码器转置卷积上采样 一个卷积块 self.up nn.ConvTranspose2d(base_channels * 2, base_channels, kernel_size2, stride2) self.dec1 self._block(base_channels * 2, base_channels) # 跳连接拼接后通道翻倍 self.out nn.Conv2d(base_channels, num_classes, kernel_size1) def _block(self, in_ch, out_ch): return nn.Sequential( nn.Conv2d(in_ch, out_ch, kernel_size3, padding1), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue), nn.Conv2d(out_ch, out_ch, kernel_size3, padding1), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue), ) def forward(self, x): e1 self.enc1(x) e2 self.enc2(nn.MaxPool2d(2)(e1)) d1 torch.cat([self.up(e2), e1], dim1) return self.out(self.dec1(d1))_block内每个卷积后面都接 BatchNorm这是遥感数据训练稳定的必要条件。Landsat 不同波段的量纲差异很大虽然做了归一化但模型内部特征值仍然会漂移BatchNorm 能显著缓解这种内部协变量偏移。另外注意跳连接拼接后通道数是原来的两倍dec1的输入通道要相应改成base_channels * 2很多人第一次搭 U-Net 都是在这里维度对不上报错。3.2 训练主循环损失函数与关键超参数的取舍地物分类的训练 loss 最常见的方案是CrossEntropyLoss。Landsat 地物类别通常是水体、植被、裸土、建筑、农田这样几类各类别像素占比极不均衡水体可能只占 3%裸土可能占 40%。这时候给CrossEntropyLoss传一个权重向量让占比小的类别得到更大的惩罚权重是见效最快的改进。import torch.optim as optim from torch.utils.data import DataLoader # 类别权重示例背景0不参与水体占比小权重给大 # class_weights 长度必须等于模型输出通道数 class_weights torch.tensor([0.0, 1.2, 0.8, 1.0, 0.7, 3.0], dtypetorch.float32).cuda() criterion nn.CrossEntropyLoss(weightclass_weights, ignore_index0) model SimpleUNet(in_channels6, num_classes6).cuda() optimizer optim.AdamW(model.parameters(), lr1e-3, weight_decay1e-4) scheduler optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max50) # 训练主循环伪代码batch 组织略 for epoch in range(50): model.train() running_loss 0.0 for images, labels in train_dataloader: images, labels images.cuda(), labels.cuda() optimizer.zero_grad() outputs model(images) # shape: (B, C, H, W) loss criterion(outputs, labels) # labels: (B, H, W) 长整型 loss.backward() optimizer.step() running_loss loss.item() scheduler.step()优化器选AdamW而不是老式的Adam。地物分类任务里 L2 正则对抑制边界过拟合有效Adam 的实现里 weight decay 和 L2 正则不是一回事AdamW 修正了这一点。学习率 1e-3 是很稳的起步值配合CosineAnnealingLR在 50 个 epoch 内从 1e-3 降到接近 0基本不需要手调学习率曲线。batch size 在 12GB 显存上设 16patch 128×1286 波段输入刚好能跑显存小就把 batch 降到 8 或者 patch 降到 96。3.3 类别不平衡不只是加权重那么简单类别不平衡是遥感分类里最容易让模型翻车的问题。水体在大部分区域影像里占比很低训练时模型发现「预测成植被」的 loss 降低最快就会倾向于把什么都预测成植被结果水体像碎片一样被吃掉这跟现实中专业人员的目视结果完全没法比。加权重是最简单的处理但权重怎么给有讲究。直接按「类别像素占比倒数」给权重往往会让模型对少量类别过度敏感把阴影、暗色屋顶全部分到水体里。我一般把权重上限压在 3.0下限拉到 0.7宁可少数类别漏分也不要误分一片。另一个补救办法是专门挖「难例」把训练中 loss 最大的那些 64×64 小块挑选出来追加训练这些小块通常集中在边界和水体边缘补一轮之后 IoU 提升非常明显。如果数据作者有精力做数据增强旋转、翻转、小角度的随机裁剪对 Landsat 这种平移不变性强的影像帮助很大。但注意不要做强光照增强Landsat 影像已经是地表反射率加亮度扰动等于伪造辐射信息真实应用时不可能有这种变化。4. 推理与成图从模型输出到完整分类结果4.1 滑窗预测与重叠拼接边缘错位问题的解法训练时切了 patch推理时同样要切 patch但推理的滑窗步长必须小于 patch 尺寸否则 patch 边界处的预测结果会因为感受野缺失出现明显的块状接缝。常用的做法是步长取 patch 的一半比如 patch 128、步长 64推理完每一块之后把中心 64×64 区域作为有效区域写入对应位置边界部分丢弃。下面是一个完整的滑窗推理示例输出为整景影像的类别图。def sliding_window_inference(model, full_image, patch_size128, stride64, batch_size8): 全图推理。stride patch_size 保证每个像素都至少位于一个 patch 的中心区域。 返回 uint8 类型的类别索引图。 model.eval() height, width, _ full_image.shape result np.zeros((height, width), dtypenp.uint8) count np.zeros((height, width), dtypenp.float32) # 按步长生成左上角坐标 rows list(range(0, height - patch_size 1, stride)) cols list(range(0, width - patch_size 1, stride)) if rows[-1] ! height - patch_size: rows.append(height - patch_size) if cols[-1] ! width - patch_size: cols.append(width - patch_size) with torch.no_grad(): for r in rows: for c in cols: patch full_image[r:rpatch_size, c:cpatch_size] # (128, 128, 6) patch_tensor torch.from_numpy(patch).permute(2, 0, 1).unsqueeze(0).float().cuda() logits model(patch_tensor) # (1, C, 128, 128) pred torch.argmax(logits, dim1).squeeze(0).cpu().numpy() # (128, 128) # 只保留中心区域边缘丢弃 start_r (patch_size - stride) // 2 start_c (patch_size - stride) // 2 end_r start_r stride end_c start_c stride result[rstart_r:rstart_cstride, cstart_c:cstart_cstride] pred[start_r:end_r, start_c:end_c] count[rstart_r:rend_r, cstart_c:cend_c] 1 # 理论上每个像素至少被覆盖一次count 可能为 0 的点用近邻填充兜底 unfilled count 0 if unfilled.any(): print(fwarning: {unfilled.sum()} pixels unfilled, applying nearest fill) # 简化处理对未填充坐标直接做最近邻插值 from scipy.ndimage import distance_transform_edt idx distance_transform_edt(unfilled, return_distancesFalse, return_indicesTrue) result[unfilled] result[tuple(idx[i][unfilled] for i in range(2))] return result推理的显存峰值是训练的一半左右batch_size 可以适当放大。预测量大的场景下逐 patch 推理的耗时主要花在 Python 循环和 GPU kernel 启动上一般的 RTX 3060 跑一景 7800×7800 影像大约需要 15 到 25 分钟这是正常水平。4.2 保存成 GeoTIFF坐标系和元数据别丢推理结果只是 numpy 数组要落回地理空间必须把原始影像的 transform、crs 一起写进输出文件。这一步如果漏了分类结果在 GIS 软件里就是一张没有地理定位的普通图片完全没法用。def save_result_tiff(result, profile, output_path): 把分类结果保存为 GeoTIFF。 profile 直接复用原始影像的改 dtype 和 count 即可。 out_profile profile.copy() out_profile.update(dtyperasterio.uint8, count1, compresslzw, nodata0) with rasterio.open(output_path, w, **out_profile) as dst: dst.write(result, 1) dst.update_tags(AREA_OR_POINTArea)用 LZW 压缩是因为分类结果大量连续区域是同一个值压缩比很高文件比原始 uint16 影像小很多。nodata 设成 0 是沿用背景类约定如果原始影像有真正的无值区域需要把那些区域在标签里就处理成 255 而不是 0避免和背景类混淆。5. 地物分类源码避坑训练到成图最常见的 5 个坑5.1 标签和影像错位几个像素模型精度原地不动现象训练 loss 正常下降验证精度也还行但把预测结果叠加到影像上肉眼检查发现地物边界整体偏向一侧或者边界处是宽窄不一的黑边。原因矢量栅格化时没有严格复用影像的 transform使用了矢量数据自己的分辨率或者坐标范围另一个常见来源是影像在预处理过程中被重投影过但标签还是老的坐标系。解决标签栅格化前先用gdf.crs.to_string()和profile[crs].to_string()做比对不同就直接to_crs(profile[crs])栅格化时显式传transformtransform和out_shape...。更简单的自查方法把标签数组和影像的近红外波段做成半透明叠加肉眼看一下河流边界对不对得上。5.2 验证集精度高换一景影像全部崩盘现象训练集和验证集是在同一景影像上随机采样切 patch 的训练时验证 mIoU 可达 90% 以上部署到相邻时相或者相邻轨道的新影像总体精度掉到 70% 以下水体区域尤其明显。原因空间自相关。同一景影像里相邻像素高度相关随机划分 patch 导致训练集和验证集互相「泄露」模型学到的更多是空间位置特征而不是地物光谱特征。出了这景影像位置特征失效。解决数据划分必须以空间区域为单位而不是以 patch 为单位。最严格的做法是用完全独立的影像做验证如果数据有限至少按遥感影像的行方向切分保证训练和验证之间没有空间相邻的 patch中间留出至少一个 patch 宽度的缓冲带更稳。5.3 阴影和暗色屋顶被全部划成水体现象模型在水体和阴影上反复横跳明明是建筑物阴影预测结果清一色是水体或者反过来水体边缘大量归属为阴影类。原因水体在可见光波段是低反射率阴影区域同样低反射率两者光谱曲线在 B2-B4 上高度相似CNN 虽然能用纹理区分但光谱相似的 patch 占比一大模型偏向用光谱捷径决策。解决第一尽量在 Landsat 影像的 B5 近红外波段上验证水体近红外吸收很强、反射极低阴影地表的近红外反射通常比纯水高这个波段是区分水体和阴影的关键第二给水体类别加更高的 loss 权重是双刃剑权重过高会引发误分建议同时做难例挖掘第三如果阴影太多考虑在预处理阶段加一个基于 B5/B2 比值的阴影掩膜把高置信阴影区域在 loss 计算时置为 ignore。5.4 直接对 16 位整数影像做训练loss 曲线乱跳现象数据加载后没有归一化就直接进模型训练 loss 前几个 epoch 降不下来甚至出现 NaN或者把 uint16 影像直接除以 255所有波段被压到 0 到 60 的窄区间模型特征提取失效。原因Landsat SR 产品是 uint16数值范围 0-65535除以 255 得到的不是 0-1 的反射率而是 0-257 的错误区间。模型对输入量纲极其敏感BatchNorm 虽然能在一定程度上兜底但这种量级错误会让第一层卷积的梯度极不稳定。解决统一走 2.2 的分位数归一化流程或者用官方公式reflectance DN / 10000缩放到真实反射率范围然后再做标准化。至少确保模型输入的数值范围在 0-1 左右RGB 波段均值在 0.05-0.3 之间是健康的状态。5.5 训练速度慢、显存溢出不是模型的锅现象一开训练12GB 显存的卡直接 OOM勉强跑起来一个 epoch 要 20 分钟50 个 epoch 得十几个小时。原因patch 开到 256×256、batch 开到 8 或 16再加上全 11 个波段输入显存需求指数上涨。很多人以为 patch 越大效果越好Landsat 30 米分辨率下 128×128 的感受野已经足够覆盖大型地物256 带来的提升微乎其微显存和耗时代价却很大。解决显存优化的优先级是先把输入波段砍到 6 个 → patch 降到 96 或 128 → batch size 降到 4 或 8 → 开启 PyTorch 的torch.cuda.amp.autocast混合精度训练显存直接减半。混合精度在 V100/Turing 架构之后的显卡上对精度影响很小遥感地物分类这种任务可以放心用。6. 精度验证与扩展看完这些指标再决定要不要深挖训练完成的模型不能只看 loss 曲线要产出一份可追踪的验证报告。针对保留的独立验证影像跑完推理后计算混淆矩阵、总体精度 OA、Kappa 系数和各类别的 IoU。其中 IoU 比 OA 更有参考价值因为 OA 会被植被这种占比大的类别拉高水体 IoU 即使只有 0.3OA 也可能显示 90%。from sklearn.metrics import confusion_matrix, cohen_kappa_score def evaluate_prediction(label_true, label_pred, class_names): label_true/label_pred: 展平后的 1D 数组仅包含参与评估的类别。 cm confusion_matrix(label_true, label_pred) oa cm.diagonal().sum() / cm.sum() iou_per_class [] for i in range(cm.shape[0]): intersect cm[i, i] union cm[i, :].sum() cm[:, i].sum() - intersect iou_per_class.append(intersect / union if union 0 else 0.0) kappa cohen_kappa_score(label_true, label_pred) # 逐类别 IoU 打印方便定位是哪一类在拖后腿 for name, iou in zip(class_names, iou_per_class): print(f{name}: IoU{iou:.4f}) print(fOA{oa:.4f}, Kappa{kappa:.4f}) return cm, oa, kappa, iou_per_class如果验证结果中某一类 IoU 明显低于其他类不要急着调模型结构先回去看这个类别的训练样本数和样本纯净度。我遇到过裸土 IoU 一直上不去的 case最后发现是标签里把收割后的农田错分成了裸土标签噪声导致模型学到的是农田纹理验证时却要求它输出裸土。先清洗标签再谈调参常常比改模型更有效。模型跑通之后源码包的扩展空间主要在三个方向。其一是输入从单时相改成多时相把不同月份的两景影像按波段维拼接输入通道变成 12可以让模型学会「同一块地在不同季节的光谱变化规律」对农作物分类精度提升显著。其二是把 U-Net 编码器替换成预训练的 EfficientNet 或 ResNet 骨干利用 ImageNet 预训练权重做迁移学习在训练样本少于 5000 个 patch 的场景下收益明显。其三是完全跳开逐像素分类尝试基于 Transformer 的 SegFormer 模型对地物边界细节的保持更好代价是显存占用更高、推理更慢。我自己的经验是先跑通当前这套 U-Net 流程把数据、标签、验证体系都理顺再考虑上述升级否则一步跳到新模型只会让问题排查变得复杂。希望这套基于 cnn 卷积神经网络做 landsat 影像地物分类的流程能帮你在自己的数据上少踩几个坑。本文还有配套的精品资源点击获取