简介一套基于PyTorch的CNN深度学习遥感影像地物分类项目源码面向人工智能、遥感、自动化、电子信息等专业的高校师生与从业者适用于毕业设计、课程设计或项目初期演示。代码经严格测试可正常运行包含数据切块、模型训练、新影像预测三个Python脚本并配有说明文档便于快速上手与二次开发。压缩包共10个文件约14.88MB除脚本外还包括2个Landsat影像TIFF及配套xml/tfw坐标信息以及1个已训练好的h5模型权重整体结构清晰覆盖从数据处理到预测分类的完整链路。目前已有87人学习下载。通过它可系统了解Landsat遥感影像地物分类的实践流程包括训练样本制作、CNN模型搭建与参数保存、新数据预测等关键环节适合作为深度学习与遥感交叉方向的学习参考也可直接用于相关课程设计或项目演示。1. 拿到一景 Landsat 影像怎样让 CNN 替你完成地物分类做遥感的人大概都有过这种经历领导丢过来一景 Landsat 8 影像说下周把这块区域的地类图给我你打开 ENVI 开始一波一波地选样本、跑最大似然分完一看农田和裸地糊成一片阴影里的水体直接消失。这个场景正是CNN 深度学习遥感影像地物分类Landsat 数据处理 Python 源码这个标题想解决的问题——把传统人工目视解译和基于像元的统计分类换成卷积神经网络自动提取空间特征让分类结果更干净、边界更完整、重复性更好。这篇笔记面向两类人一类是刚接触深度学习的遥感从业者手里有 Landsat 数据但不知道从哪一步开始另一类是学过 CNN 分类但被数据预处理折腾过的人比如波段怎么选、样本怎么打、精度为什么上不去。你可以把下文当成一份能直接照着做的技术方案从 Landsat 数据的波段特性讲到 CNN 建模再到训练避坑和精度验证最后落到结果怎么导回 GIS 里用。整个流程用的都是 Python 生态GDAL 处理栅格、PyTorch 搭网络验证集上的 Kappa 系数能做到 0.85 以上这是一条成熟、可复现、不需要超算也能跑的路线。2. Landsat 数据从下载到预处理模型还没建坑已经埋了三层2.1 Landsat 8/9 的波段组合为什么不能把 11 个波段全丢给 CNNLandsat 8 和 9 的 OLI 传感器有 11 个波段但其中沿海波段Band 1、卷云波段Band 9和两个热红外波段Band 10/11用途相对特殊做常规地物分类时很少全部参与建模。常见做法是把 Band 2、3、4 作为蓝、绿、红真彩色组合再加上 Band 5 近红外、Band 6 短波红外 1、Band 7 短波红外 2组成 6 个波段参与分类。这六个波段对植被、水体、建筑物、裸地的区分度最好也是 USDA 等机构做土地覆盖产品时的标准输入。有人会问既然 CNN 能自己学特征是不是波段越多越好实际不是。热红外波段空间分辨率是 100 米重采样到 30 米后边缘模糊反而干扰浅层卷积核学习锐利特征。卷云波段主要用来检测薄云参与分类只会引入噪声。我一般先用 6 波段组合跑一个基线如果某些类别混淆严重再针对性加入短波红外或计算指数特征而不是一上来全波段堆进去——输入维度上去后显存翻倍、训练时间变长精度却不涨甚至下降。2.2 L1 级产品必须做大气校正跳过这一步等于白干从 USGS EarthExplorer 下载的 Landsat 数据分两个层级L1 级是辐射校正后的原始 DN 值L2 级是已经做过大气校正的地表反射率产品。如果拿到的是 L1 级直接拿 DN 值去训练模型会遇到一个非常隐蔽的问题——不同时相影像的 DN 值范围受太阳高度角和大气状态影响模型在一个季节学的特征在另一个季节完全不通用。这不是模型鲁棒性差而是输入数据本身不在同一个物理量纲上。正确做法是做辐射定标加大气校正。辐射定标是把 DN 值转成辐亮度公式是 L ML × DN AL其中 ML 和 AL 来自影像头文件MTL 文件里的 RADIANCE_MULT_BAND_x 和 RADIANCE_ADD_BAND_x 字段。大气校正把辐亮度转成地表反射率常用的工具是 ENVI 的 FLAASH 模块也可以用开源的 Py6S 库在 Python 里完成。L2 级产品可以直接跳过这一步它已经提供了 SR 波段这也是为什么我建议新手优先下载 L2 级数据把精力放在后面更关键的环节。2.3 波段合成、裁剪和归一化的具体操作预处理落到代码层面第一步是把单波段 GeoTIFF 合成多波段文件。用 GDAL 读取六个波段然后叠加成一个数组同时保留地理参考信息。这里有一个容易踩的坑GDAL 读栅格的顺序和波段编号不一定一致必须先遍历 GetRasterBand 确认或者直接用 gdal_translate 按文件名匹配。from osgeo import gdal, gdal_array import numpy as np # 打开六个单波段文件按顺序合成多波段数组 files [B2.TIF, B3.TIF, B4.TIF, B5.TIF, B6.TIF, B7.TIF] bands [] ds gdal.Open(files[0]) geo_transform ds.GetGeoTransform() projection ds.GetProjection() cols, rows ds.RasterXSize, ds.RasterYSize for f in files: ds gdal.Open(f) band ds.GetRasterBand(1).ReadAsArray().astype(np.float32) bands.append(band) ds None # 堆叠成 (C, H, W) 顺序和 PyTorch 输入格式对齐 stacked np.stack(bands, axis0) print(合成后形状:, stacked.shape, 值范围:, stacked.min(), -, stacked.max())代码逻辑说明先把第一个文件的地理变换和投影信息保存下来这是后面写回 GeoTIFF 时必需的。逐波段读入后堆叠成 (C, H, W) 形状C 是波段数 6H 和 W 是影像的高和宽。输出前先看一眼值范围——如果最大值超过 10000说明这是反射率缩放了 10000 倍的 L2 产品训练前要除以 10000 归一化到 0-1 区间。归一化这步很关键但很多新手会在这里翻车。Landsat L2 地表反射率的数值范围是 0 到 1但实际存储时乘了 10000直接扔给 CNN 会导致梯度爆炸。常规做法是除以 10000或者对每个波段做 z-score 标准化。我倾向用全局最小最大值归一化而不是每景影像单独标准化这样不同时间的数据在特征空间里是一致的模型迁移时不需要重新训练。参数说明裁剪是为了让样本覆盖不同地理区域时保持尺寸一致。CNN 分类通常按固定大小的 patch 输入比如 64×64 或 128×128。裁剪时要设置 stridestride 小于 patch 尺寸可以产生重叠样本增加训练数据量但重叠太多会导致训练集和验证集空间上相关精度虚高。我一般训练集 stride 取 patch 的 50%验证集完全不相交。3. CNN 模型结构与训练策略从零搭一个能用的地物分类器3.1 为什么遥感影像分类要用 CNN 而不是前馈神经网络这是被问得最多的问题。全连接网络处理图像时要把每个像素拉成一维向量一张 64×64×6 的 patch 展平后有 24576 个输入节点第一层全连接的参数量轻松上百万训练困难不说还完全丢失了像素之间的空间结构。CNN 的卷积核在局部窗口内滑动卷积天然假设相邻像素有相关性这个归纳偏置让它用远少于全连接的参数量学会了边缘纹理形状这些中层特征。遥感影像和自然图像还有一个重要区别地物的尺度差异很大。一棵树的树冠在 30 米分辨率下可能只有几个像素一片农田却占据上百个像素。CNN 通过堆叠卷积层和池化层形成层级特征——浅层关注边缘纹理深层关注语义类别这种多尺度表达能力是传统机器学习方法比如随机森林 纹理特征很难手工设计的。传统方法需要人为计算 GLCM 纹理、NDVI 指数、形状特征等费时费力而且很难覆盖所有地物类型CNN 把这些特征学习自动化了这是它在这个领域快速替代传统方法的核心原因。3.2 网络结构怎么选小数据集用小网络别一上来就上 ResNetLandsat 影像和 ImageNet 不一样它只有 6 个波段空间分辨率 30 米地物边界相对平滑不需要特别深的网络去拟合极端复杂的纹理。我踩过的坑是一开始直接套 ResNet50 做迁移学习结果在样本量只有几千个 patch 的情况下严重过拟合训练精度 99%验证精度卡在 78% 上不去。后来换成自己搭的一个 4 层卷积网络验证精度反而涨到 88%。这里的原因其实不复杂预训练模型的权重是为三通道自然图像设计的Landsat 的 6 个波段在数值分布和物理含义上差异巨大迁移时要么把预训练权重的前三层丢掉重训要么修改输入层适配 6 通道。两种做法都意味着大量参数要从头学而遥感分类的样本量远不如 ImageNet参数越多越容易翻车。下面是一个针对 Landsat 设计的轻量分类网络参数量不到 50 万普通显卡几分钟就能训练一轮。import torch import torch.nn as nn class LandsatCNN(nn.Module): def __init__(self, num_bands6, num_classes6): super().__init__() self.features nn.Sequential( # 第一层卷积6 波段输入32 个 3x3 卷积核 nn.Conv2d(num_bands, 32, kernel_size3, padding1), nn.BatchNorm2d(32), nn.ReLU(inplaceTrue), nn.MaxPool2d(2), # 第二层卷积通道翻倍 nn.Conv2d(32, 64, kernel_size3, padding1), nn.BatchNorm2d(64), nn.ReLU(inplaceTrue), nn.MaxPool2d(2), # 第三层卷积继续提取高层语义 nn.Conv2d(64, 128, kernel_size3, padding1), nn.BatchNorm2d(128), nn.ReLU(inplaceTrue), nn.AdaptiveAvgPool2d(1), # 全局池化自适应输入尺寸 ) self.classifier nn.Linear(128, num_classes) def forward(self, x): x self.features(x) x x.view(x.size(0), -1) return self.classifier(x)模型结构说明三个卷积层都用了 3×3 卷积核加 padding1保证特征图尺寸只在池化层减半。BatchNorm2d 放在卷积和 ReLU 之间对遥感数据尤其重要——Landsat 不同波段的数值范围差异很大BatchNorm 能稳定每一层的输入分布。最后用 AdaptiveAvgPool2d(1) 做全局平均池化把任意尺寸的特征图压缩成 128 维向量这样训练时可以随意调整 patch 大小模型结构不用改。参数说明num_bands6 对应前面预处理合成的 6 个波段num_classes6 对应你要分的六个地类比如耕地、林地、水体、建筑、裸地、草地。如果你有更多类别只需改这个数字最后一层全连接的输出维数会自动调整。池化做了三次两个 MaxPool2d 加一个 AdaptiveAvgPool输入 64×64 的 patch 经过两个池化层后变成 16×16再经过全局池化直接降成 1×1计算量非常小。3.3 训练参数怎么调学习率、批次大小和类别不均衡训练参数是 CNN 分类最容易出现玄学的地方。先说学习率遥感影像分类任务我用得最多的是初始学习率 0.001、配合余弦退火调度。0.01 太大Loss 在前几个 epoch 会震荡不收敛0.0001 太慢训练 50 个 epoch 精度还在缓慢爬升。一个经验做法是先用单个 batch 试跑看 Loss 是否在下降然后按 10 倍步长搜索合适区间。批次大小受显存限制我一般设 32 或 64。批次太小BatchNorm 统计量不稳定批次太大一个 epoch 内参数更新次数少收敛慢。这里有个容易忽略的点类别不均衡。如果把水体、林地、农田、建筑四类按原始比例切 patch建筑类样本往往只有农田的十分之一模型会对稀有小类别基本无视——准确率看着有 95%建筑的召回率只有 40%。解决不均衡的常规做法是加权采样或加权损失。加权损失更简单直接计算每个类别的样本占比把占比的倒数归一化作为权重传给 CrossEntropyLoss 的 weight 参数。我实际用下来这个方案对小类别召回率的提升非常明显而且不改变模型结构训练代价为零。from torch.utils.data import WeightedRandomSampler # labels 是所有训练 patch 的类别标签数组 class_counts np.bincount(labels) class_weights 1.0 / class_counts sample_weights class_weights[labels] sampler WeightedRandomSampler(sample_weights, num_sampleslen(labels), replacementTrue) # DataLoader 传入 sampler 即可 train_loader torch.utils.data.DataLoader( train_dataset, batch_size32, samplersampler )代码逻辑说明WeightedRandomSampler 按照每个样本的权重进行有放回抽样权重越大的类别被抽中的概率越高。class_weights 用类别样本占比的倒数占比少的类别权重反而大天然平衡了每个 epoch 内各类别出现的次数。replacementTrue 允许同一个 patch 在一个 epoch 内被重复采样这是必要的否则样本少的类别凑不够一个 batch。4. Python 源码搭建与踩坑排查从环境配置到训练跑通的完整流程4.1 环境准备与数据加载器GDAL、PyTorch 和 Dataset 类的正确组合技术栈可以简化为四个部分GDAL 负责读写 GeoTIFFNumPy 做数组运算PyTorch 负责搭建和训练 CNN最后用 scikit-learn 计算精度指标。装环境时最容易出问题的是 GDAL 和 PyTorch 的版本共存——GDAL 依赖的 numpy 版本太老会导致 import 报错这种情况常见于直接 pip install 的场景别去手动降级 numpy否则 PyTorch 也会跟着出问题。写数据加载器时核心是把读取 patch和读取标签封装进 torch.utils.data.Dataset 的子类。这里有个容易被忽略的细节把整景影像一次性读入内存再在getitem里按索引切片。如果每调用一次getitem就用 GDAL 重新读一次磁盘训练速度会被 IO 拖垮到原来的十分之一。影像大小在几千乘以几千像素时6 个波段的 float32 数组也就几百 MB完全放得下内存。class LandsatPatchDataset(torch.utils.data.Dataset): def __init__(self, image_array, label_array, patch_size64, stride64): # image_array: (C, H, W) 归一化后的影像数组 # label_array: (H, W) 标签数组0 表示无数据区 self.image image_array self.label label_array self.patch_size patch_size self.patches [] h, w label_array.shape for y in range(0, h - patch_size 1, stride): for x in range(0, w - patch_size 1, stride): patch_label label_array[y:ypatch_size, x:xpatch_size] # 仅保留中心像素有有效标签的 patch if patch_label[patch_size // 2, patch_size // 2] ! 0: self.patches.append((y, x)) def __len__(self): return len(self.patches) def __getitem__(self, idx): y, x self.patches[idx] img self.image[:, y:yself.patch_size, x:xself.patch_size] lab self.label[yself.patch_size//2, xself.patch_size//2] return torch.from_numpy(img), torch.tensor(lab, dtypetorch.long)代码里的一个技巧用中心像素的标签代表整个 patch 的类别。这样做避免了 patch 边缘落在两种地物交界处导致标签不确定的问题同时让预测结果更平滑。如果训练时用整个 patch 的众数做标签模型会学到多数类吞并少数类的倾向边界处的地物会被系统性地抹掉。数据加载器的参数需要注意patch_size 和 stride 相同表示不重叠切块适合样本量充足的情况想要增广数据可以把 stride 调成 patch_size 的一半但验证集千万别这样用否则精度会虚高。4.2 训练循环与样本增广三行代码的事别自己重复造轮子训练循环本身并不复杂PyTorch 的标准写法是每个 epoch 遍历 DataLoader、前向传播、计算 Loss、反向传播、更新权重。但遥感影像分类有一个独特问题过拟合。Landsat 影像同质区域的 patch 高度相似如果不做增广模型会在训练集上死记硬背验证精度停滞不前。常规做法是随机水平翻转、垂直翻转和旋转 90 度。这三个操作完全符合遥感影像的物理特性——地物朝向不固定翻转旋转不会改变类别语义。不推荐随机裁剪和颜色抖动前者会改变 patch 中心对应的标签位置后者会破坏地表反射率的物理一致性。import torchvision.transforms as T train_transform T.Compose([ T.RandomHorizontalFlip(p0.5), T.RandomVerticalFlip(p0.5), T.RandomRotation((90, 90)), # 只在 90 度的倍数上旋转不产生锯齿 ]) def train_one_epoch(model, loader, optimizer, criterion, device): model.train() total_loss, correct, total 0.0, 0, 0 for images, labels in loader: images, labels images.to(device), labels.to(device) # 增广需要 4D 张量 (B, C, H, W)直接调用即可 images train_transform(images) outputs model(images) loss criterion(outputs, labels) optimizer.zero_grad() loss.backward() optimizer.step() total_loss loss.item() * images.size(0) correct (outputs.argmax(1) labels).sum().item() total images.size(0) return total_loss / total, correct / total训练逻辑说明train_transform 作用于每个 batch即同一个 batch 内的样本每次 epoch 看到的翻转状态都不相同等效于把训练集扩大了数倍。注意这里没有做标准化——因为数据在预处理阶段已经归一化过了再套 ImageNet 的 mean/std 反而会把反射率分布搞歪。这是一个很容易被迁移学习教程误导的地方。训练时的设备选择和精度设置也有讲究。GPU 可用时把模型和数据搬到 CUDA用混合精度训练torch.cuda.amp能省一半显存尤其当 patch 尺寸调到 128 时优势明显。CPU 训练不是不行但一个 epoch 可能要跑十几分钟迭代调参会非常痛苦建议至少用一块入门级 GPU。4.3 避坑专题五次翻车记录每一条都在源码里真实发生过现象一训练 Loss 下降很快但验证集精度一直在 60% 上下徘徊上不去。原因是样本标签本身有误——打标签时把阴影里的林地标成了草地模型学到的是暗的就是草地这种错误规律。解决方法是检查训练集里每个类别的 patch 可视化随机抽 20 张图看一眼有问题重新标注。现象二预测结果出现明显的条带噪声地物边界呈锯齿状不是平滑的曲线。原因是 stride 过大导致 patch 之间重叠不足相邻 patch 的预测结果在边界处不连续。解决方法是预测时用滑窗加投票即 stride 设为 patch_size 的四分之一多个 patch 对同一像素的预测结果取众数。现象三验证集精度很高但整景影像的预测图上出现了大量椒盐噪声单个像素被分成不同的类别。原因是逐 patch 预测时没有考虑空间上下文CNN 对孤立像素的预测不稳定。解决方法是加一个简单的多数滤波3×3 窗口取众数噪声基本能消除想更精细可以做 CRF 后处理但多数滤波在大多数场景下已经够用。现象四模型在 A 区域的影像上精度很高但换到相邻 B 区域就明显下降。原因是两个区域的物候期不同植被的反射率差异大模型没见过 B 区域的特征分布。解决方法是在训练集里混入多个时相的影像或者干脆别跨时相用同一个模型Landsat 每 16 天重访一次模型按季节各自训练反而更省事。现象五程序跑着跑着内存爆掉进程被系统杀掉。原因是整景影像以 float64 类型读入内存6 个波段 8000×8000 像素就是 3 GB。解决方法是读入后立刻转成 float32 并除以 10000同时在 Dataset 初始化时用 del 释放中间变量必要时把影像按块读取、分块预测再拼接结果。这五条踩坑记录里最值得重视的是第三条和第四条。椒盐噪声可以用后处理解决但跨时相泛化问题如果不提前规划模型上线后返工成本极高。我的建议是做项目前先问清楚影像的时相范围若跨越两个季节务必每个季节都准备样本。5. 精度验证与结果导出混淆矩阵、Kappa 系数和 GeoTIFF 输出模型训练完不是终点遥感地物分类项目要交付的是带地理坐标的分类图不是一组预测数字。最后一个环节有两件事必须做一是用一套严格的空间不重叠验证集计算精度二是把预测结果写回 GeoTIFF 文件保留地理参考信息以便在 ArcGIS 或 QGIS 里叠加显示。精度验证的核心指标不是整体准确率而是每一类的生产者精度和用户精度以及整体 Kappa 系数。生产者精度衡量的是这个类别被正确识别出来的比例用户精度衡量的是被分到这个类别的像元里真正属于该类别的比例。两者差异大说明某一类容易被误分为其他类比如建筑和裸地的混淆在 30 米分辨率下非常常见。from sklearn.metrics import confusion_matrix, cohen_kappa_score # y_true 是验证集真实标签y_pred 是模型预测标签 cm confusion_matrix(y_true, y_pred) kappa cohen_kappa_score(y_true, y_pred) # 每一类的生产者精度对角线值除以该类的真实样本总数 producer_acc cm.diagonal() / cm.sum(axis1) # 每一类的用户精度对角线值除以该类的预测样本总数 user_acc cm.diagonal() / cm.sum(axis0) print(Kappa:, round(kappa, 4)) print(类别 生产者精度 用户精度) for i in range(len(class_names)): print(f{class_names[i]:6} {producer_acc[i]:.3f} {user_acc[i]:.3f})评估代码说明混淆矩阵的行是真实类别列是预测类别。对角线上的值越大说明分类越准但真正要关注的是非对角线上的混淆模式——比如建筑被分到裸地的数量是否过多这决定了你是否需要引入新的特征或增加这两类的训练样本。Kappa 系数小于 0.8 时分类图拿去出报告会比较勉强优先回头增加训练样本而不是调模型。结果导出到 GeoTIFF 是一个容易出问题的地方。预测时模型输出的是每个 patch 中心像素的类别你需要按 patch 的坐标位置把类别填回一个与原始影像尺寸一致的数组里然后把数组和前面保存的 geo_transform、projection 一起传给 GDAL 写盘。这一步写错会导致分类图和影像对不上像元整体偏移。一个值得沿用的习惯是预测时用一个独立于训练集的非重叠窗口并且保存一张置信度图——即模型输出的 Softmax 概率最大值。置信度低的区域往往是地物边界、阴影或混合像元所在的位置后续人工检查时直接按置信度阈值筛出可疑区域能省下大量目视核查时间。我自己每次出图后都会把置信度低于 0.6 的像元标出来叠加在真彩色影像上那种感觉就像给自己的模型做了一次彻底的健康检查。最后分享一个经验如果项目预算允许尽量在训练集里加入不同季节、不同云量的影像哪怕每个季节只有几百个 patch也能显著提升模型在实际工作流里的生存能力。毕竟 Landsat 时隔 16 天就重访一次收集不同时相的影像成本很低这比调参涨的那一两个百分点要值钱得多。做遥感深度学习快六年我最大的感受是数据预处理决定精度的下限模型结构决定精度的上限而踩坑越多的工程师出的分类图越干净。希望这篇笔记能帮你少走我走过的那些弯路。本文还有配套的精品资源点击获取