简介：一套基于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 个波段全丢给 CNN
Landsat 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, axis=0) print("合成后形状:", stacked.shape, "值范围:", stacked.min(), "-", stacked.max())代码逻辑说明：先把第一个文件的地理变换和投影信息保存下来，这是后面写回 GeoTIFF 时必需的。逐波段读入后堆叠成 (C, H, W) 形状，C 是波段数 6，H 和 W 是影像的高和宽。输出前先看一眼值范围——如果最大值超过 10000，说明这是反射率缩放了 10000 倍的 L2 产品，训练前要除以 10000 归一化到 0-1 区间。
归一化这步很关键，但很多新手会在这里翻车。Landsat L2 地表反射率的数值范围是 0 到 1，但实际存储时乘了 10000，直接扔给 CNN 会导致梯度爆炸。常规做法是除以 10000，或者对每个波段做 z-score 标准化。我倾向用全局最小最大值归一化而不是每景影像单独标准化，这样不同时间的数据在特征空间里是一致的，模型迁移时不需要重新训练。
参数说明：裁剪是为了让样本覆盖不同地理区域时保持尺寸一致。CNN 分类通常按固定大小的 patch 输入，比如 64×64 或 128×128。裁剪时要设置 stride，stride 小于 patch 尺寸可以产生重叠样本，增加训练数据量，但重叠太多会导致训练集和验证集空间上相关，精度虚高。我一般训练集 stride 取 patch 的 50%，验证集完全不相交。
3. CNN 模型结构与训练策略：从零搭一个能用的地物分类器
3.1 为什么遥感影像分类要用 CNN 而不是前馈神经网络
这是被问得最多的问题。全连接网络处理图像时要把每个像素拉成一维向量，一张 64×64×6 的 patch 展平后有 24576 个输入节点，第一层全连接的参数量轻松上百万，训练困难不说，还完全丢失了像素之间的空间结构。CNN 的卷积核在局部窗口内滑动卷积，天然假设相邻像素有相关性，这个归纳偏置让它用远少于全连接的参数量学会了"边缘""纹理""形状"这些中层特征。
遥感影像和自然图像还有一个重要区别：地物的尺度差异很大。一棵树的树冠在 30 米分辨率下可能只有几个像素，一片农田却占据上百个像素。CNN 通过堆叠卷积层和池化层形成层级特征——浅层关注边缘纹理，深层关注语义类别，这种多尺度表达能力是传统机器学习方法（比如随机森林 + 纹理特征）很难手工设计的。传统方法需要人为计算 GLCM 纹理、NDVI 指数、形状特征等，费时费力而且很难覆盖所有地物类型；CNN 把这些特征学习自动化了，这是它在这个领域快速替代传统方法的核心原因。
3.2 网络结构怎么选：小数据集用小网络，别一上来就上 ResNet
Landsat 影像和 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_bands=6, num_classes=6): super().__init__() self.features = nn.Sequential( # 第一层卷积，6 波段输入，32 个 3x3 卷积核 nn.Conv2d(num_bands, 32, kernel_size=3, padding=1), nn.BatchNorm2d(32), nn.ReLU(inplace=True), nn.MaxPool2d(2), # 第二层卷积，通道翻倍 nn.Conv2d(32, 64, kernel_size=3, padding=1), nn.BatchNorm2d(64), nn.ReLU(inplace=True), nn.MaxPool2d(2), # 第三层卷积，继续提取高层语义 nn.Conv2d(64, 128, kernel_size=3, padding=1), nn.BatchNorm2d(128), nn.ReLU(inplace=True), 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 卷积核加 padding=1，保证特征图尺寸只在池化层减半。BatchNorm2d 放在卷积和 ReLU 之间，对遥感数据尤其重要——Landsat 不同波段的数值范围差异很大，BatchNorm 能稳定每一层的输入分布。最后用 AdaptiveAvgPool2d(1) 做全局平均池化，把任意尺寸的特征图压缩成 128 维向量，这样训练时可以随意调整 patch 大小，模型结构不用改。
参数说明：num_bands=6 对应前面预处理合成的 6 个波段；num_classes=6 对应你要分的六个地类，比如耕地、林地、水体、建筑、裸地、草地。如果你有更多类别，只需改这个数字，最后一层全连接的输出维数会自动调整。池化做了三次（两个 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_samples=len(labels), replacement=True) # DataLoader 传入 sampler 即可 train_loader = torch.utils.data.DataLoader( train_dataset, batch_size=32, sampler=sampler )代码逻辑说明：WeightedRandomSampler 按照每个样本的权重进行有放回抽样，权重越大的类别被抽中的概率越高。class_weights 用类别样本占比的倒数，占比少的类别权重反而大，天然平衡了每个 epoch 内各类别出现的次数。replacement=True 允许同一个 patch 在一个 epoch 内被重复采样，这是必要的，否则样本少的类别凑不够一个 batch。
4. Python 源码搭建与踩坑排查：从环境配置到训练跑通的完整流程
4.1 环境准备与数据加载器：GDAL、PyTorch 和 Dataset 类的正确组合
技术栈可以简化为四个部分：GDAL 负责读写 GeoTIFF，NumPy 做数组运算，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_size=64, stride=64): # 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:y+patch_size, x:x+patch_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:y+self.patch_size, x:x+self.patch_size] lab = self.label[y+self.patch_size//2, x+self.patch_size//2] return torch.from_numpy(img), torch.tensor(lab, dtype=torch.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(p=0.5), T.RandomVerticalFlip(p=0.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(axis=1) # 每一类的用户精度：对角线值除以该类的预测样本总数 user_acc = cm.diagonal() / cm.sum(axis=0) 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 天就重访一次，收集不同时相的影像成本很低，这比调参涨的那一两个百分点要值钱得多。做遥感深度学习快六年，我最大的感受是：数据预处理决定精度的下限，模型结构决定精度的上限，而踩坑越多的工程师，出的分类图越干净。希望这篇笔记能帮你少走我走过的那些弯路。
本文还有配套的精品资源，点击获取
