Landsat遥感影像地物分类实战：从预处理到CNN训练完整流程

简介基于PyTorch的CNN深度学习遥感影像地物分类Python源码聚焦Landsat数据处理与地物分类任务适合高校人工智能、遥感、自动化等专业学生及从业者用于毕业设计、课程设计或项目入门也可作为深度学习和遥感图像处理的基础练习。压缩包共10个文件大小约14.88MB包含3个Python脚本分别实现训练影像切片生成、模型训练、新数据预测、1个训练好的CNN模型文件.h5、2个GeoTIFF影像及配套xml/tfw辅助文件另有说明文档便于快速理解工程结构。已有86人学习下载。资源提供从数据预处理到模型训练与预测的完整流程代码经严格测试可正常运行内置示例影像和预训练权重可直接复现地物分类效果也便于在此基础上修改扩展适合实践学习与二次开发。1. 从Landsat到地物分类这个任务到底在解决什么用CNN做遥感影像地物分类很多人一上来就调网络结构但真正让分类精度掉到70%以下的往往是Landsat数据连预处理都没做干净。Landsat的DN值不是反射率云和阴影还混在波段里直接喂给PyTorch卷积核学到的是传感器噪声而不是地表特征。下面这条路径会把Landsat数据下载、辐射定标、切片建库、CNN训练到分类图输出的完整流程拆开所有代码都面向Python源码落地。适合刚接触深度学习和遥感、想在本地跑通第一个地物分类项目的工程师。你会看到预处理为什么决定了CNN的上限也能直接抄走最小实现。2. Landsat数据预处理把原始影像变成CNN能吃的张量Landsat预处理是CNN分类最重要的环节。几何校正和辐射定标、云掩膜等做不好后面全是白费。数据下载后拿到的是L1级产品先别急着读进模型要把波段结构、反射率换算和切片方式先理清楚。2.1 认识Landsat的存储结构波段顺序和元数据Landsat 8/9的一景L1TP压缩包内是多个单波段GeoTIFF命名如LC08_L1TP_xxx_B2.TIF。B2到B7分别负责蓝、绿、红、近红外和两个短波红外。地物分类常把这六个波段全部用上但为了匹配常见深度学习框架的输入规模很多人会先用B2到B5四个波段训练跑通后再加B6和B7。下表给出了这几个波段对地物的敏感度。波段波长µm对哪类地物更敏感B2 蓝0.452–0.512水体与建筑不透水面和水体的反射差异明显B3 绿0.533–0.590健康植被的反射峰常用于NDVI计算B4 红0.636–0.673裸土、道路和城区人工地物B5 近红外0.851–0.879植被高反射与水体和土壤对比最强B6 短波红外11.566–1.651土壤湿度、干燥植被B7 短波红外22.107–2.294矿物、火烧迹地与建筑阴影识别拿到数据后先用rasterio快速确认文件能正常读取python -c import rasterio; srcrasterio.open(LC08_L1TP_xxx_B2.TIF); print(src.width, src.height, src.crs)这段命令输出宽度、高度和坐标系。如果crs是EPSG:326XX说明是UTM投影L1TP已经做过几何校正不需要再处理几何偏移。rasterio是处理GeoTIFF的Python库用pip install rasterio安装。真正需要操心的是辐射预处理下一步就做这件事。2.2 辐射定标与简单大气校正从DN值到地表反射率Landsat原始文件存的是16位整数DN值要先读元数据MTL.txt里的RADIOMETRIC_RESCALING_MULT_BAND_x换算成表观反射率。如果不想对接6S大气模型常见做法是直接用L2级地表反射率产品但课题交付往往只给L1级源码里就得自己实现最小辐射定标。下面是一个用numpy处理单波段的简版import numpy as np def dn_to_reflectance(dn, mult, add, sun_elev): DN转反射率。mult/add来自MTL.txtsun_elev为太阳高度角。 radiance mult * dn add refl np.pi * radiance / (1361.0 * np.sin(np.deg2rad(sun_elev))) return np.clip(refl, 0.0, 1.0)这段代码的关键参数有三个mult和add是每个波段的缩放因子必须按波段从MTL.txt里单独取1361是太阳常数单位是W/m²·µm不同传感器版本可能需要查文档确认np.sin把太阳高度角换算成天顶角的余弦注意角度制到弧度制的转换。实际项目中气溶胶会让反射率整体抬高最简易的修正是DOS1暗像元减法def dos1(reflectance): dark np.percentile(reflectance, 0.1) return np.maximum(reflectance - dark, 0.0)暗像元取全波段反射率第0.1百分位用来近似大气程辐射贡献。在大多数陆地区域这比完全不做大气校正要稳。做完之后把各波段按(Band, H, W)叠起来下一条流程就会用到。2.3 制作训练切片用栅格窗口而不是全图读入一景Landsat影像普遍接近8000×8000像素直接读成numpy数组会吃穿内存。训练CNN时只需要局部patch。下面这段代码生成64×64的样本切片import rasterio from rasterio.windows import Window def iterate_patches(img_path, label_path, channels, patch64, stride32): with rasterio.open(img_path) as img_src, rasterio.open(label_path) as lbl_src: for r in range(0, img_src.height - patch 1, stride): for c in range(0, img_src.width - patch 1, stride): win Window(c, r, patch, patch) x img_src.read(channels, windowwin) y lbl_src.read(1, windowwin) yield x, y, r, cchannels是波段索引列表比如[2,3,4,5]对应B2到B5stride小于patch时生成重叠样本相当于数据增强label_path需要提前把地物矢量标签栅格化成与影像完全对齐的GeoTIFF分辨率也必须是30m。生成后要过滤掉标签区域全是背景值的块否则模型会学到大量空样本。这一步做好了训练集的平衡性才有保障。3. CNN模型构建与训练从零搭一个地物分类网络选择CNN其实是选择了感受野。传统逐像素分类只看单点的光谱值而CNN判断一个像素是不是建筑会看到周围几十米范围内的纹理、阴影和邻域结构。Landsat 30m分辨率下一栋房子可能只有一两像素但建筑区和裸地的邻域纹理完全不同这是CNN能比随机森林多拿几个点精度的原因。3.1 为什么用CNN而不是传统分类器传统方法先算NDVI、NDWI等指数再丢给SVM或随机森林。这些指数能区分高植被和水体但面对城市密集建筑与裸地的边界时很浑浊。CNN可以在卷积层自动学到类似“近红外减红”的特征也可以在更高层组合这些特征输出复杂空间关系。此外卷积的权重共享让模型对位置变化不敏感同一块地物出现在影像不同位置时输出不会因为像素坐标偏移而突变。PyTorch是目前做深度学习最常用的Python库装好CUDA环境就能开始。3.2 用PyTorch实现一个轻量CNN网络下面是一个输入(B,4,64,64)、输出(B,5)的小网络运行代码量很小CPU也能跑import torch.nn as nn class LandCoverCNN(nn.Module): Landsat地物分类网络输入多波段patch输出像元类别概率。 def __init__(self, in_channels4, num_classes5): super().__init__() self.features nn.Sequential( nn.Conv2d(in_channels, 32, 3, padding1), nn.BatchNorm2d(32), nn.ReLU(inplaceTrue), nn.Conv2d(32, 64, 3, stride2, padding1), nn.BatchNorm2d(64), nn.ReLU(inplaceTrue), nn.AdaptiveAvgPool2d(1), ) self.classifier nn.Linear(64, num_classes) def forward(self, x): return self.classifier(self.features(x).flatten(1))第一层卷积输入4个波段输出32个特征图padding1让空间尺寸保持64×64。第二层卷积stride2把特征图缩小到32×32通道升到64这样空间分辨率降低的同时通道数增加。AdaptiveAvgPool2d(1)把任意空间尺寸压成1×1向量避免后面全连接层依赖固定输入尺寸。全连接输出5类对应水体、植被、裸地、建筑、阴影。如果你想区分更多类别改num_classes并确保训练标签从0开始连续编号。这个网络的卷积层只有两层主要照顾Landsat这种30m分辨率的场景。地形纹理更复杂时可以在第二个卷积后面再接一个相同通道的残差块但训练时间会明显增加。3.3 训练参数设定把一套默认参数写成表格训练配置建议直接用下面这张表。参数建议值说明影像输入尺寸64×64覆盖约1.92km地面范围太大GPU内存不够波段数4B2到B5可按需改成6batch size648GB显存够用OOM就降到32学习率1e-3Adam默认值验证集震荡时改1e-4损失函数CrossEntropyLoss多分类标准训练轮数30观察验证集不再下降就早停优化器Adambeta10.9, beta20.999最小编码循环import torch import torch.nn as nn from tqdm import tqdm device torch.device(cuda if torch.cuda.is_available() else cpu) model LandCoverCNN(in_channels4, num_classes5).to(device) criterion nn.CrossEntropyLoss() optimizer torch.optim.Adam(model.parameters(), lr1e-3) for epoch in range(30): model.train() for img, lbl in tqdm(train_loader, descfepoch {epoch1}): img, lbl img.to(device), lbl.to(device) optimizer.zero_grad() out model(img) loss criterion(out, lbl) loss.backward() optimizer.step()train_loader是PyTorch的DataLoader每轮返回预处理后的张量(B,4,64,64)和标签(B,)。标签必须转成long类型因为CrossEntropyLoss不接受float当类别索引。每个epoch结束后要跑一遍验证集输出准确率根据情况调整学习率。4. 处理Landsat大影像的四个坑内存、拼接、样本与泛化训练完模型接着要跑到整景影像上。这一步比训练更容易翻车以下几个坑几乎每个Landsat分类项目都会遇到。4.1 内存溢出用patch-by-patch推理代替整图预测整景8000×8000像素如果不切块直接forward显存几乎立刻爆掉。常见做法是滑动窗口推理一次只处理一个64×64块推理时关闭梯度with torch.no_grad(): for r in range(0, height - patch 1, stride): for c in range(0, width - patch 1, stride): x extract_patch(r, c, patch) out torch.softmax(model(x.to(device)), dim1).squeeze().cpu().numpy() result[r:rpatch, c:cpatch] out.argmax(0)stride如果等于patch速度最快但块边界明显显存仍然不足时把batch size设为1不要在batch维度一次叠多个窗口。4.2 拼接边界效应重叠窗口加权平均直接把patch预测结果拼到result里会让patch边缘出现“马赛克”状跳变原因是相邻patch能看到的上下文不同。解决方法是让滑动窗口有重叠stride取patch的一半然后对重叠区域加权累加越靠近patch中心权重越大。下面用汉宁窗做加权融合import numpy as np def merge_patches(h, w, patch, stride, model, extract): acc np.zeros((h, w, num_classes), np.float32) weight np.zeros((h, w), np.float32) w_row np.hanning(patch).reshape(-1, 1) w_col np.hanning(patch).reshape(1, -1) win_w w_row w_col with torch.no_grad(): for r in range(0, h - patch 1, stride): for c in range(0, w - patch 1, stride): x extract(r, c, patch) prob torch.softmax(model(x), dim1).squeeze(0) # (C,H,W) prob prob.permute(1,2,0).cpu().numpy() # (H,W,C) acc[r:rpatch, c:cpatch] prob * win_w[..., None] weight[r:rpatch, c:cpatch] win_w return acc / weight[..., None]win_w是patch×patch的二维矩阵通过None索引广播成和概率图相同的形状。重叠区域预测值累加后除以权重边缘处自然平滑最终得到一张无明显接缝的概率图。4.3 类别不均衡加权损失与随机采样真实场景里农田和森林面积很大水体或建筑经常只占很小比例。直接用CrossEntropy模型会把小类别全叫成大类别。最直接的做法是统计训练集每个类别的像素占比计算权重from collections import Counter label_counts Counter(labels.flatten()) total sum(label_counts.values()) class_weights torch.tensor([total / max(label_counts[i], 1) for i in range(5)]).to(device) criterion torch.nn.CrossEntropyLoss(weightclass_weights)这是倒数加权某类只有1%像素时权重会接近100容易放大噪声所以常把权重归一化到0.5到3之间。除加权损失外也可以给DataLoader加WeightedRandomSampler对稀有类样本过采样两种方法能同时使用。4.4 过拟合与迁移数据增强和预训练权重训练集通常只有几千到几万张patch模型很容易记住训练区域的细节验证集正常但换一幅影像精度骤降。增强手段放在transform里from torchvision import transforms train_transform transforms.Compose([ transforms.RandomHorizontalFlip(), transforms.RandomVerticalFlip(), transforms.RandomRotation(15), transforms.Normalize(mean[0.5]*4, std[0.5]*4) ])注意Normalize的mean和std要从训练集统计不能固定写0.5。旋转会改变地理方向但地物分类中水体、农田对方向不敏感可接受。如果想用ResNet18预训练权重由于ImageNet输入是3通道Landsat是4或6通道常见做法是只取B4、B3、B2组成伪RGB并把网络第一层改成3通道。这会丢掉B5近红外信息因此很多项目干脆不加载预训练因为卫星波段分布和自然图像差异太大迁移收益有限。坑现象提前预防CUDA OOM训练或推理时显存耗尽patch改小、batch减半分类图呈块状预测掩模有方块边界stride减半重叠加权小类召回率低水体被划分为植被加权损失过采样验证好独立精度差训练和测试区域太近按区域划分train/test并增强5. 从源码到可复用脚本整理输出与分类后处理现在分类图已经预测出来但还只是numpy数组需要写回GeoTIFF做精度评估再封装成命令行工具。5.1 把概率图写成带地理坐标的GeoTIFF分类结果要能叠加到QGIS或ArcGIS里必须复制原始Landsat的transform和crs。下面把(H,W,C)的概率数组写入多波段GeoTIFFwith rasterio.open(LC08_L1TP_xxx_B2.TIF) as ref: profile ref.profile.copy() profile.update(countnum_classes, dtypefloat32) with rasterio.open(class_probs.tif, w, **profile) as dst: dst.write(probs.transpose(2,0,1).astype(float32))probs是(H,W,C)数组transpose成(C,H,W)才能被rasterio.write接收。硬分类图则新建一个count1、dtypeuint8的profile把argmax结果写进去。这样输出文件自带地理坐标能直接和原始影像叠加。5.2 用混淆矩阵和Kappa系数验证预测完独立的测试区域后需要用指标量化分类质量from sklearn.metrics import confusion_matrix, cohen_kappa_score y_true test_labels.flatten() y_pred test_pred.flatten() cm confusion_matrix(y_true, y_pred) kappa cohen_kappa_score(y_true, y_pred) print(cm) print(fkappa {kappa:.3f})Kappa大于0.75说明一致性较好混淆矩阵能看出哪两类常混分。Landsat分类中水体与建筑阴影很容易混因为阴影区在可见光波段的水体响应很接近。这时把B7短波红外加入通道能明显改善阴影误分。5.3 用argparse把流程封装成命令行工具可复用的源码包至少要暴露两个命令train.py负责训练predict.py负责推理。predict.py的参数可以这样设计import argparse parser argparse.ArgumentParser(descriptionLandsat land cover prediction) parser.add_argument(--img, requiredTrue, helpLandsat多波段堆叠后的tif) parser.add_argument(--checkpoint, requiredTrue, help训练好的模型权重) parser.add_argument(--out, defaultclass.tif, help输出分类图) parser.add_argument(--patch-size, typeint, default64) parser.add_argument(--stride, typeint, default32) args parser.parse_args()命令行调用变成python predict.py --img stacked.tif --checkpoint model.pt --out result.tif --stride 32stride默认取patch_size的一半正好对应重叠加权方案。把波段选择、归一化参数和模型类替换成你的版本这套结构就能批量处理多景Landsat影像。你的CNN结构、预处理函数、后处理权重都能在同一个仓库里找到对应文件切换传感器时只需要改波段索引和统计值。本文还有配套的精品资源点击获取