简介遥感影像地物分类是土地利用监测、环境评估、农业估产等领域的基础工作。传统方法依赖光谱特征与人工设计CNN卷积神经网络则能自动提取纹理、形状和上下文等深层特征分类精度与鲁棒性更优适用范围覆盖土地利用变化监测、作物分类与灾害评估等场景。这套完整工程面向遥感研究者、GIS开发者和深度学习入门者提供基于Landsat影像的地物分类Python解决方案。压缩包共10个文件含3个Python脚本分别完成影像切片、模型训练和新数据预测另有2个tif示例影像、1个h5训练好的模型以及xml、tfw辅助文件整体14.49MB结构紧凑、可直接运行。已有182人学习下载借助TensorFlow/PyTorch与OpenCV使用者无需从零训练加载h5模型即可快速验证分类效果结合README还能理清数据预处理、模型训练与预测的完整流程便于教学演示、科研实验和实际遥感项目二次开发。1. 基于CNN的Landsat地物分类到底比传统方法强在哪如果你手头有一幅30米分辨率的Landsat影像想分出水体、耕地、林地、建筑和裸土传统的随机森林或SVM也能做但你会发现两个问题一是靠人工挑的特征如NDVI、纹理决定了精度上限二是地物边界处总是一团糟。换成Python CNN深度学习之后同样的训练样本整体分类精度通常能提升8~15个百分点边界也更干净。这个标题里的“源代码训练好的模型”组合其实是告诉你一件事不需要从零训练一个大网络直接在预训练模型上微调或做推理是遥感从业者最常用的落地路径。这篇文章从数据预处理、网络选型、训练到推理制图把CNN做Landsat地物分类的完整链路讲清楚适合想真正跑通一版分类结果、又不想被复杂理论劝退的从业者。2. 把地物分类当成图像任务CNN的输入、输出与三种常见网络结构2.1 多光谱波段与RGB的差别通道数决定网络输入CNN处理普通照片时输入是三通道RGB而Landsat影像的陆地卫星系列通常有多个波段。以Landsat 8 OLI为例常见的处理输入是可见光波段蓝、绿、红、近红外、短波红外中的几个组合。你要做的第一个决定是到底用哪几个波段作为网络输入。我一般会优先选6个波段蓝B2、绿B3、红B4、近红外B5、短波红外1B6、短波红外2B7。为什么不把全部波段塞进去因为像气溶胶波段B1和热红外波段B10/B11要么受大气影响大要么分辨率不同直接输入反而会让网络把噪声当特征。如果你在USGS下载的是Landsat Collection 2 Level-2产品波段已经是地表反射率省去了大气校正但要注意每个波段的缩放因子。通道数的变化直接改网络第一层的卷积核深度。比如PyTorch里nn.Conv2d(in_channels6, out_channels32, kernel_size3)如果你的输入是RGB三通道这里就要改成in_channels3。很多新手拿着别人的RGB训练代码直接跑多光谱数据报错往往就出在这一行。严格来说Landsat的影像不是“照片”但CNN不关心物理含义它只把每个波段当成一个特征通道。你把NDVI、NDWI这些指数也拼进去当作额外通道也可以但要注意标准化方式必须与训练时一致。2.2 逐像素分类与语义分割两种思路对模型的要求用CNN对遥感影像做地物分类实际有两条路线一种是“逐像素分类”也叫图像块分类patch-based把每个像元周围的N×N窗口作为一个小图输入CNN输出该像元的类别。另一种是“语义分割”直接输入整幅或大块影像输出每个像元的类别典型模型是U-Net、SegNet。逐像素分类的好处是思路简单训练和推理都很直接很多老项目就是拿一个5层CNN在滑窗上跑。坏处是计算量大、有重复计算而且每个像元独立决策空间连续性不好——这正是你在结果里看到“椒盐噪声”的原因。语义分割则在网络结构里引入了编码-解码和多尺度特征合并输出直接是整张分类图边界保持好得多。现在的主流方案基本都转向语义分割尤其是U-Net。从落地角度我建议直接用U-Net。原因有三点第一U-Net在医学图像、遥感图像上被验证得最多网络结构小一张1080显卡就能训练第二它的跳跃连接能把浅层的纹理信息传到深层这对地物边界很重要第三代码工业级成熟随便搜一下就有几十个可复写的实现。后面我给的训练源代码也是精简版U-Net。2.3 怎么选网络U-Net、SegNet还是小CNN如果你的训练样本很少每类只有几百个标注像元逐像素分类的小CNN反而更稳因为参数少不容易过拟合。如果样本有几千块以上U-Net优势明显。SegNet和DeepLabV3等更大模型在遥感上也能用但显存和训练时间成倍增加对Landsat这种类别不多、地物相对简单的任务收益很小。还有一个很多教程不提的细节Landsat影像的类别分布极度不均比如大片农田和零星建筑。如果直接用交叉熵损失模型会把所有像元预测成农田因为这么做损失也不大。这时要改损失函数常见做法是加权交叉熵或Dice损失。我在自己的项目里就吃过这个亏——训练准确率99%全图分类结果却是一片绿。后面第5章会专门讲这个坑。3. 从Landsat原始影像到训练集数据预处理与标签制作3.1 下载与辐射定标L1还是L2缩放因子怎么用USGS EarthExplorer和GloVis是下载Landsat影像的常见去处。这里建议直接选Collection 2 Level-2产品因为L1级别是原始辐射亮度需要自己用ENVI或Python做辐射定标和大气校正而L2已经是地表反射率省掉一个最容易错的大气校正环节。下载时注意影像云量云量大于5%的景最好剔除尤其是训练数据云和云影会被模型学成一种“地物”推理时给你制造假异常。L2产品的每个波段有两个文件一个实际数据如*SR_B2.TIF一个质量评估文件如*QA_PIXEL.TIF。数据文件里存的是整型数值要得到真实地表反射率需要乘以0.0001。很多人直接拿原始DN值训练结果数值范围在0~65535之间和预训练模型期望的0~1输入差着数量级损失直接爆炸。所以我一般会先做一次预处理脚本import rasterio import numpy as np def load_landsat_band(band_path, scale0.0001): with rasterio.open(band_path) as src: band src.read(1).astype(np.float32) valid band ! src.nodata band[~valid] np.nan return band * scale # 缩放因子这个函数把Landsat L2波段读取为浮点型地表反射率无效像元置为nan。注意band ! src.nodata要先判断因为L2的填充值可能是0或65535不同影像不一致。读取之后你还需要做一步——把所有波段的nan像元统一掩膜不然输入到CNN里会出现“黑洞”。3.2 制作样本块滑动窗口裁剪与标签栅格对齐训练数据是“影像块 对应标签块”。如果你的标签是已有的土地分类产品如GlobeLand30或人工矢量化的地类边界需要先栅格化并重采样到和Landsat一样的30米分辨率、一样的地理范围。最常见的问题就是影像块和标签块错位半个像元导致模型学到错误边界。我一般用rasterio把标签和影像投影到同一坐标系用reproject做最近邻重采样避免插值产生新类别。做训练样本时我不会直接把整景影像输入网络。Landsat单景大约7000×7000像元显存放不下。常见做法是滑动窗口裁剪成256×256或128×128的小块。窗口大小取决于地物对象的尺度农田块整齐256合适城市建筑破碎128更合适。裁剪步长可以小于窗口大小做重叠采样让边缘地物也能被充分学习。import numpy as np from patchify import patchify # image: (H, W, C) 的 float32 数组, label: (H, W) 的 int 数组 patches_img patchify(image, (256, 256, image.shape[2]), step128) patches_lbl patchify(label, (256, 256), step128)这段代码用patchify把整景影像切成256×256小块步长128意味着相邻块有50%重叠。重叠采样会让训练样本数量翻倍但对边界地物来说利大于弊。注意patchify切出来的数组是个多维数组需要把前两维展平才能得到(样本数, 256, 256, 通道数)的标准数据集。如果你不想引入额外依赖也可以手动写双层循环切片效果一样。切好的小块按类别筛选如果某块里某一类占比低于1%这种块其实很难学混在训练集里反而干扰收敛我一般直接丢弃。另外建议按场景分景切块不要把同一景影像的块既放训练集又放验证集否则验证精度虚高换成新影像立刻翻车。3.3 数据增强旋转、翻转、随机亮度调整的合理范围遥感地物不像自然图片那样有多姿态变化但你依然可以做增强。常见的做法是随机水平/垂直翻转、90度/180度旋转、随机亮度抖动。尺度变化对Landsat这类固定分辨率数据意义不大不建议缩放增强因为缩放会改变地物纹理的真实尺度。import torchvision.transforms as T from PIL import Image transform T.Compose([ T.RandomHorizontalFlip(p0.5), T.RandomVerticalFlip(p0.5), T.RandomApply([T.ColorJitter(brightness0.2, contrast0.1)], p0.3), T.ToTensor(), ])这里是PyTorch常见的transform组合。brightness0.2表示亮度在0.8~1.2倍之间随机变化。对Landsat地表反射率来说这个范围算安全超过0.3就会把水体提亮到和裸土混淆。contrast我一般只给0.1遥感影像的对比度本身跨度不大增强过猛反而学不进去。一个容易忽略的要点增强必须只作用在影像上不能作用在标签上。所以上面的transform只传给影像分支标签块单独处理。如果你在写数据加载器时图省事把影像和标签拼成一个元组同时transform标签就会被插值变成非整数类别损失函数直接报错。4. 用Python PyTorch训练CNN地物分类模型源代码与参数调优4.1 网络结构定义一个可跑的U-Net精简版下面这段是U-Net的精简实现输入6波段图块输出5个地物类别水体、农田、林地、建筑、裸土。我刻意把通道数控制在16起步普通CPU训练也能勉强跑有GPU更好。这里用PyTorch因为它的动态图方便我们在训练中调试形状不匹配问题。import torch import torch.nn as nn class UNetSmall(nn.Module): def __init__(self, in_channels6, n_classes5): super().__init__() self.enc1 self._block(in_channels, 16) self.enc2 self._block(16, 32) self.enc3 self._block(32, 64) self.pool nn.MaxPool2d(2) self.up2 nn.ConvTranspose2d(64, 32, kernel_size2, stride2) self.dec2 self._block(64, 32) self.up3 nn.ConvTranspose2d(32, 16, kernel_size2, stride2) self.dec3 self._block(32, 16) self.out nn.Conv2d(16, n_classes, kernel_size1) def _block(self, in_c, out_c): return nn.Sequential( nn.Conv2d(in_c, out_c, 3, padding1), nn.BatchNorm2d(out_c), nn.ReLU(inplaceTrue), nn.Conv2d(out_c, out_c, 3, padding1), nn.BatchNorm2d(out_c), nn.ReLU(inplaceTrue), ) def forward(self, x): e1 self.enc1(x) e2 self.enc2(self.pool(e1)) e3 self.enc3(self.pool(e2)) d2 torch.cat([self.up2(e3), e2], dim1) d2 self.dec2(d2) d3 torch.cat([self.up3(d2), e1], dim1) d3 self.dec3(d3) return self.out(d3)这里说明几个关键点。in_channels6对应前面选的6个Landsat波段如果你换成了RGB三波段这里改成3就行。nn.ConvTranspose2d是转置卷积用来把特征图尺寸恢复到输入尺寸与普通上采样不同它能学习到更合适的插值方式。torch.cat把编码器对应层的特征拼到解码器这就是U-Net的“跳跃连接”高频边界信息能通过这条路直接传到输出层。如果分类类别大于5把n_classes改成实际类别数。但注意类别数变化后最后一层卷积输出维度变了训练好的模型不能直接跨类别复用。4.2 训练脚本数据加载、损失函数与checkpoint训练部分重点在损失函数和优化器。我推荐CrossEntropyLoss但要开启ignore_index参数把标签里的无效像元如-1跳过。类别不平衡的另一个做法是给每类设权重权重用1 - 类别频率计算具体在踩坑章讲。import torch.optim as optim from torch.utils.data import DataLoader, TensorDataset # train_inputs: (N, 6, 256, 256), train_labels: (N, 256, 256) dataset TensorDataset(torch.from_numpy(train_inputs), torch.from_numpy(train_labels)) loader DataLoader(dataset, batch_size8, shuffleTrue) model UNetSmall(in_channels6, n_classes5) weights torch.tensor([0.8, 0.5, 0.7, 1.5, 1.2]) # 按类频率反比设置 criterion nn.CrossEntropyLoss(ignore_index-1, weightweights) optimizer optim.Adam(model.parameters(), lr1e-3) for epoch in range(50): model.train() for imgs, lbls in loader: preds model(imgs) # (B, 5, 256, 256) loss criterion(preds, lbls.long()) # 标签必须是 int64 optimizer.zero_grad() loss.backward() optimizer.step() if epoch % 10 0: torch.save(model.state_dict(), flandsat_cnn_epoch{epoch}.pth)batch_size8是针对256×256输入在8GB显存下的经验值如果你显存小降到4或2。lr1e-3是Adam的默认起点跑20轮后可以降到1e-4做细调。每次保存的state_dict是模型权重文件也就是标题里说的“训练好的模型”。注意这里用的是model.state_dict()而非整个model因为只存权重更小并且加载时只要重新实例化同样结构的网络就能恢复。一个容易犯的错model(imgs)的输出形状是(B, 5, H, W)而CrossEntropyLoss期望输入是(B, C, H, W)标签是(B, H, W)的长整型。如果标签是one-hot编码或者形状是(B, 1, H, W)会报维度不匹配需要先squeeze。4.3 训练好的模型保存格式与精度评估训练结束后建议保存两种格式一种.pth用于继续训练或微调另一种.onnx用于部署和快速推理。转ONNX的代码很简单dummy torch.randn(1, 6, 256, 256) torch.onnx.export(model, dummy, landsat_unet.onnx, input_names[image], output_names[class_map], dynamic_axes{image: {0: batch}, class_map: {0: batch}})dynamic_axes声明batch维度是动态的这样推理时一次可以喂入任意张影像块。精度评估不能只看类别平均准确率我建议同时计算每类的F1-score和总体Kappa系数。因为水体面积大、农田面积大总体准确率会被多数类别拉高建筑这一小类做错了根本看不出来。验证时放在不同时相或不同区域的Landsat影像上测才能反映模型真实泛化能力。5. 踩坑与排查Landsat分类任务的5个常见问题5.1 训练精度高但验证精度低是过拟合还是标签错现象训练集准确率99%验证集只有70%且验证损失在下降后突然反弹。原因通常有两个一是模型过拟合到训练影像的纹理细节二是训练块和验证块来自同一景影像导致的空间相关性被误认为高精度。解决先检查数据划分是否按“景”隔离。如果同一个湖面被切成训练块和验证块模型相当于见过的题目的变体验证分自然虚高。我一般会在切块前按影像文件ID分组保证同一景的块只进一个集合。如果按景隔离后精度仍低再加数据增强或缩小网络通道数。5.2 全图被分成同一类别急着改网络先查损失和标签现象训练过程loss一直在降但推理结果整幅图都是农田。原因多半是类别极度不平衡交叉熵损失被多数类别主导也可能是标签图里只有农田这一类的像元被正确标注其他类全是背景。解决观察训练集每类样本比例如果农田占了80%至少要把损失权重设成权重 1 - 比例让模型对少数类更敏感。更有效的方式是改用加权Dice损失它不按像元数量而按类别交集计算天然处理不平衡class DiceLoss(nn.Module): def forward(self, preds, labels, eps1e-6): preds torch.softmax(preds, dim1) # (B, C, H, W) onehot torch.eye(preds.shape[1])[labels] # one-hot, 避开无效标签 onehot onehot.permute(0, 3, 1, 2).to(preds.dtype) intersection (preds * onehot).sum(dim(2, 3)) union preds.sum(dim(2, 3)) onehot.sum(dim(2, 3)) dice (2 * intersection eps) / (union eps) return 1 - dice.mean()上面这个DiceLoss实现里有个陷阱labels如果有-1值直接做one-hot会越界所以需要先把无效标签置0或在索引前排除。计算dice时取所有类的平均少数类的损失会被单独拉高效果立竿见影。5.3 投影与坐标不一致导致的“整幅偏移”现象训练和推理用的是不同源的数据推理结果在ArcGIS里叠到底图上整幅偏移了1~2个像元或完全错位。原因多半是训练标签用的是WGS84影像用的是UTM投影或者两个坐标系虽然名字一样但基准面不同。解决第一步先用rasterio检查两个数据的crs和transform必须完全一致。第二步重采样标签时使用nearest而不是bilinear因为双线性插值会把类别值变成小数。第三步做严格对齐测试——取影像中的道路或水岸线确认和标签边界没有半像元左右偏移再进入训练。5.4 显存溢出与内存爆炸batch size和图像块的博弈现象跑训练时程序报CUDA out of memory或者切数据时本机内存直接打满。原因很直接256×256×6的输入8个batch不算大但U-Net的中间特征图数量是输入的数倍加上梯度占用8GB显存跑不动很正常。解决优先把batch_size降到2再把输入块从256切成128。如果还想继续大块可以用梯度累积每2个小batch累加一次梯度等效于大batch但显存占用不变。内存爆炸则是因为一次性把所有训练块读入内存改成PyTorch的Dataset按需读取不要用TensorDataset塞全景。5.5 加载别人训练好的模型后效果变差预处理不一致现象拿到一个训练好的.pth模型加载后跑出来的分类结果全是黑色或噪声精度还不如随机分类。原因几乎都是预处理不一致——可能是影像没有乘以0.0001缩放因子也可能是没做均值方差归一化或用了全新的波段组合。解决办法是回看训练时的输入范围如果训练时用的是0~1反射率推理前必须保证影像块也是0~1。Landsat L2原始DN值在0~65535不缩放就输入模型所有激活函数都被打饱和。这里没有后悔药只能把预处理流程固定成一个函数训练和推理都调用它。6. 进阶用训练好的模型做整景推理输出带地理坐标的tif整景Landsat有7000×7000像元不能一次性塞进网络常见做法是滑窗推理再加回坐标信息。下面这段代码把前面保存的模型和预处理封装成完整推理流程并直接写出带投影信息的GeoTIFF省去你在ArcGIS里手工赋坐标的麻烦。import numpy as np import rasterio import torch from rasterio.transform import from_origin model UNetSmall(in_channels6, n_classes5) model.load_state_dict(torch.load(landsat_unet.pth, map_locationcpu)) model.eval() with rasterio.open(Landsat8_B2_B7.tif) as src: profile src.profile image src.read() # (6, H, W) transform src.transform # 归一化假设地表反射率范围0~1已经乘以scale h, w image.shape[1], image.shape[2] output np.zeros((h, w), dtypenp.uint8) overlap 32 step 256 - overlap for y in range(0, h, step): for x in range(0, w, step): y2 min(y 256, h); x2 min(x 256, w) block image[:, y:y2, x:x2] block_t torch.from_numpy(block).unsqueeze(0).float() with torch.no_grad(): pred model(block_t) # (1, C, y2-y, x2-x) cls torch.argmax(pred, dim1).squeeze(0).numpy() output[y:y2, x:x2] cls.astype(np.uint8)这里的关键参数是overlap32避免窗口边界处出现接缝。模型对每块的预测完全独立相邻块的重叠区域取后一次的结果即可要求不高时这样足够。如果希望更好可以对重叠区域按到中心的距离做加权平均但会慢不少。最后用原始投影信息写出结果with rasterio.open(result.tif, w, driverGTiff, heighth, widthw, count1, dtypeuint8, crsprofile[crs], transformtransform) as dst: dst.write(output, 1)写tif时务必带上crs和transform这两个值直接来自原始影像的profile。我吃过一次亏没有复制profile结果输出的tif在GIS软件里是“无坐标系”的还得手动二次配准非常麻烦。把推理脚本固化下来之后换一景新影像只需要改输入文件名半小时能出一整幅30米地物分类图。最后说一个习惯训练好的模型和预处理脚本永远放在同一个目录并给模型文件名加日期和输入波段数比如unet_6band_20240512.pth。不然三个月后你看着一个.pth文件根本不记得它是6波段还是3波段训练出来的查起来全是玄学。这个方向值得投入因为Landsat免费、重访周期短CNN分类模型一旦跑通你手上的历史影像都能批量变成土地利用数据。希望这些经验能帮到你。本文还有配套的精品资源点击获取