简介面向Landsat影像的地物分类需求这套基于CNN深度学习的Python源代码提供从训练到预测的完整流程适合遥感研究人员、GIS工程师及高校相关专业学生使用。压缩包共10个文件总大小约14.49MB主要包含3个Python脚本分别用于生成影像切片、训练模型、预测新数据、训练好的CNN模型h5格式、示例TIF影像及其配套的xml/tfw文件另有README说明文档。目前已有182人学习用户可直接加载h5模型对新影像进行分类无需从零训练网络节省大量时间。CNN能自动提取纹理、形状等深层特征相比传统方法在精度和鲁棒性上更有优势特别适用于土地利用变化监测、农作物分类、环境监测等场景。整套文件结构清晰既可作为实际项目的实施基础也可用于教学演示与二次开发为遥感地物分类实践提供便利工具。1. 从 Landsat 影像到地物分类为什么 CNN 方案能直接落地做遥感的人手里大概率都有那么一两套“历史遗留”的分类代码不是用的最大似然法就是随机森林跑出来的结果在城区边缘全是椒盐噪声。当你试过把一景 30 米分辨率的 Landsat 8 OLI 影像喂进 CNN 卷积神经网络后再回头去看那些传统分类器的结果就会明白什么叫代差。这套 Python 方案最核心的价值是把“源代码 训练好的模型”直接变成可复用的分类流水线——你不用自己去标注几万张样本、不用从零设计卷积网络结构只要按固定流程做辐射定标、切图、归一化把影像数据送到模型里就能得到建筑物、水体、耕地、裸地、林地这些地物类别。四个关键点值得先说清楚为什么咱非得用 Landsat 而不去用高分辨率影像为什么模型文件敢直接在本地跑而不用联网为什么精度能做到在混淆矩阵里 Kappa 系数稳定上 0.85为什么这份代码能迁移到别的区域。这篇内容我会从源码结构和推理逻辑入手带你把数据准备、模型训练、参数调节、还有那些常见的翻车场景全部捋一遍。适合已经跑通过几个深度学习的分类任务、但没正经处理过遥感影像的人也适合拿别人的模型做推理却经常在数据预处理上栽跟头的乙方单位。2. 看懂这份源码的模型结构与推理逻辑2.1 模型选型为什么是 CNN 而不是传统分类器LandSat 影像分类和做自然图像分类本质上不是一回事。一张 ImageNet 图片只包含一个语义主体而一景 Landsat 影像覆盖几百平方公里的地表里面同时存在水体、建筑、植被、裸土、云影等多种地物类别。传统分类器包括支持向量机和随机森林它们逐个像元地提取光谱特征完全忽略邻域信息所以结果在天际线阴影区域容易产生“同物异谱”错分。CNN 的核心竞争力在于它能用卷积核同时抽取光谱特征和空间上下文——比如一个 3x3 的卷积核对应地面上 90 米 x90 米的实际范围这就能把“房子旁边一定有道路、道路周边大概率有裸地”这类空间规律捕获到。在遥感影像分类这块业界做得最多、效果也最稳的是 U-Net 及其变体。它属于编码器-解码器结构编码器用逐级下采样把浅层特征从 512x512 缩小到 16x16 甚至更小的特征图解码器再用转置卷积逐级上采样恢复分辨率。最关键的是它中间加了跳跃连接把编码器每一层的特征图直接拼接给解码器对应层这样网络在最后做像素分类时既能保留高层的语义信息知道什么东西大概率出现在哪又不会丢失浅层精细的边缘细节比如建筑和裸地之间的严格边界。如果你拿到的源码不是 U-Net 而是简单的全卷积网络 FCN也没关系原理大同小异只是没有跳跃连接边缘精度会差一些。2.2 运行环境搭建与项目文件作用先把环境跑通再谈其他。我一般建议用 Python 3.8 或者 3.10配 CUDA 环境下的 PyTorch 1.10 以上版本。如果你只是拿训练好的模型做推理不需要显卡也能跑就是速度慢一点如果你打算自己重新训练一张 GTX 1660 级别以上、显存 6G 的显卡是入门线。这部分最容易踩坑的是 GDAL 和 Rasterio 的安装——Win 10 系统上直接用 pip 装 Rasterio 经常报缺少底层依赖的错误比你 python 安装教程里配环境变量还麻烦。我的建议是直接用 conda 创建环境conda create -n landsat python3.8 conda activate landsat conda install -c conda-forge rasterio gdal pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu118参数说明这里的--index-url指定的是 PyTorch 官方针对 CUDA 11.8 的编译版本GPU 版本的 PyTorch 不能用默认 pip 源安装否则会下载成纯 CPU 版本模型跑起来就是龟速。Rasterio 用 conda 而不是 pip 装是为了避免 GDAL 的动态链库冲突——用 pip 装的 Rasterio 在导入时经常报DLL load failed这属于根本没法用的状态。文件结构方面这份源码里最常见的是一份main.py或者train.py做训练入口、一份predict.py做推理入口、模型权重文件一般是unet_landsat.pth或者.pt结尾另外还有一个config.py里面记录了波段数、类别数、均值方差这些归一化参数。你拿到手之后第一件事不是去看代码逻辑而是打开config.py确认它读的是几个波段、模型输出是几类。2.3 加载训练好的模型并对单景影像做预测我拿到的这套源码里推理逻辑比较清晰没有一堆花里胡哨的类封装。核心步骤就是读取 TIF 文件、按训练时的套路做归一化、切成固定大小的 patch 逐个预测、最后再拼回完整的地理坐标系。关键代码如下import torch import torch.nn.functional as F import rasterio from rasterio.windows import Window import numpy as np from model import UNet # 假设模型定义在 model.py 中 # 加载模型 model UNet(in_channels7, out_channels5) # Landsat 8 选7个波段分5类 checkpoint torch.load(best_model.pth, map_locationcpu) # 没有GPU就强制用CPU model.load_state_dict(checkpoint[state_dict]) model.eval() # 读取影像元信息 with rasterio.open(input_landsat.tif) as src: width, height src.width, src.height profile src.profile all_bands src.read() # 形状 (7, height, width) # 归一化到 [0,1]注意必须和训练时一致 all_bands all_bands.astype(np.float32) for b in range(all_bands.shape[0]): all_bands[b] (all_bands[b] - np.min(all_bands[b])) / (np.max(all_bands[b]) - np.min(all_bands[b]) 1e-8) # 切块推理避免一次加载过大影像导致显存溢出 patch_size 256 stride 224 # 重叠一部分减少拼接缝隙 output np.zeros((5, height, width), dtypenp.float32) count np.zeros((height, width), dtypenp.float32) for y in range(0, height - patch_size 1, stride): for x in range(0, width - patch_size 1, stride): patch all_bands[:, y:ypatch_size, x:xpatch_size] patch_tensor torch.from_numpy(patch).unsqueeze(0) with torch.no_grad(): pred model(patch_tensor) # 形状 (1, 5, patch_size, patch_size) pred F.softmax(pred, dim1).squeeze(0).numpy() output[:, y:ypatch_size, x:xpatch_size] pred count[y:ypatch_size, x:xpatch_size] 1 output output / np.maximum(count, 1) class_map np.argmax(output, axis0).astype(np.uint8) # 写出结果保留地理参考信息 profile.update(dtypeuint8, count1, compresslzw) with rasterio.open(classification_result.tif, w, **profile) as dst: dst.write(class_map, 1)这段代码的逻辑说明输入影像有 7 个波段对应 Landsat 8 的可见光、近红外和两个短波红外通道。切块推理时用到了 224 的 stride让相邻 patch 之间重叠 32 个像素这样可以大幅缓解模型输出的拼接缝效应——投影在遥感这类纹理极其复杂的影像上时如果不做重叠拼接拼出来的图会有一条一条明显的人工边界。逐 patch 预测得到的五分类概率全部累加到一个输出矩阵中最后按 count 矩阵做归一化这相当于对重叠区域做平均投票结果比硬切硬拼稳定得多。输出前用profile把原始影像的地理坐标、投影信息原封不动地写进去这样结果可以直接丢进 ArcGIS 或 QGIS 里叠加分析。代码里最后没有用torch.argmax直接在模型输出上取结果而是先做了一次F.softmax拿到每个类别的概率值这样做的好处是中间概率图还可以用于后续的阈值筛选和人工复核。如果你拿到的源码里推理脚本写得比较乱建议你按照上面这个套路自己重构一版——把并行切块、概率融合、地理写入三个环节拆开每一块独立测试确认无误之后再拼起来。3. 制作自己的 Landsat 训练数据集从 TIF 到标签3.1 数据源选择与预处理辐射定标、大气校正很多人以为别人给了一份训练好的模型自己就可以无限地在别的区域预测。实际上一旦你换到植被覆盖或土壤底色差异极大的地区模型精度会断崖式下跌所以最终还是躲不开自己制备训练数据的环节。Landsat 数据通常从 USGS 官网或者国内的地理空间数据云下载拿到手的是 L1 级产品DN 值范围在 0 到 65535 之间。直接把这些原始数值喂给神经网络训练会导致 loss 在你还没调参的时候就崩掉了——因为不同影像的太阳高度角、大气条件、传感器增益不同DN 值的分布完全不可比。所以第一步一定要做辐射定标和大气校正。辐射定标把 DN 值转换成传感器入瞳处的辐射亮度值公式涉及 R 和 M 两个增益系数这两个系数直接写在每景影像的 MTL 元数据文件里。大气校正常见做法是用 ENVI 里的 FLAASH 模块或者 6S 模型但我在这份 Python 方案里更推荐用 envi 的ENVI_DOIT或者直接用Google Earth Engine导出地表反射率产品。如果你不想引入太多手动环节实际上很多开源数据集的制作套路是直接用 Landsat 的 Collection 2 Level-2 产品它自带地表反射率和云掩膜波段省掉了一大段物理反演的地狱难度。拿到之后需要做波段筛选Landsat 8 的 OLI 传感器共有 9 个波段但对地物分类来说真正有用的是 2、3、4、5、6、7 六个波段加海岸蓝波段 1大多数 CNN 项目的输入通道就是 6 到 7 个。再多的话训练时数据量需求会陡然上升模型参数也会膨胀性价比很低。3.2 用 Python 切分影像并生成训练样本含代码做训练样本时最重要的一件事是保证标签和影像在空间范围上严格对齐。标签图可能是你以前人工解译出来的 Shapefile也可能是一个类型栅格。这里会把标签栅格转成同尺寸的int数组类别从 0 开始编号比如 0 代表背景、1 代表水体、2 代表建筑、3 代表裸地、4 代表植被。针对 Landsat 影像的裁剪策略和普通图像任务不一样因为遥感影像中一些类别比如水体只占很小的面积比例如果完全随机在整景影像上滑窗采样训练集里会出现大量全背景 patch导致模型偏向把什么都预测成背景。常见的做法是先按类别做直方图统计计算每个类别的面积占比然后对占比很小的类别做定向采样——例如水体在影像边缘、河道附近这些位置多切一点 patch保证训练数据中每个类别至少占到 15% 以上。切出来的训练集不应该直接存成 .npz而是做成文件夹里面存着image目录和label目录各自存放按顺序排列的 .npy 文件import numpy as np import rasterio from sklearn.utils.class_weight import compute_class_weight label_arr np.array([0, 1, 2, 3, 4]) # 假设5类 # 计算各类别像素占比用于后续采样策略和损失函数加权重 num_classes 5 class_count np.zeros(num_classes) for c in range(num_classes): class_count[c] np.sum(label_arr c) total_pixel np.sum(class_count) class_weight total_pixel / (num_classes * class_count) class_weight np.clip(class_weight, 1.0, 10.0) print(类别权重:, class_weight)参数说明class_weight计算使用了中位数频率平衡法分子total_pixel除以类别数得到平均像素数再除以该类别的实际像素数得到比值。对于占比极小的类别这个权重会非常大比如水体只占 1%权重就可能到 8.9。后面在训练时用这个权重做加权交叉熵损失会在梯度更新时放大对少数类别的惩罚力度不然模型练出来等于一个高精度的“背景分类器”。切 patch 的时候我用rasterio.windows.Window按固定偏移滑窗每个 patch 大小定在 256x256同时检查 patch 里有效像素的比例如果 nodata 或云掩膜的面积超过 10%就跳过这块样本。可以落地的细节是Landsat 影像的几何分辨率是 30 米256 个像素在空间上覆盖 7.68 公里这对一个地物类别来说范围足够大大到 CNN 能学到“河道两侧有堤坝”这类上下文结构。4. 训练你自己的地物分类模型超参数与关键代码4.1 数据增强策略与训练集划分遥感影像的数据增强和自然图像不太一样直接照搬随机裁剪、翻转那套方法存在隐藏的坑。比如垂直翻转对城市影像来说相当于把北上广深的地物盖到了南半球变化剧烈模型难以收敛。更稳妥的增强策略是水平翻转、垂直翻转、旋转 90 度保持地物相对角度的一致性、和多光谱通道上加一个轻微的随机伽马校正模拟不同光照条件下地物光谱的偏移。不要把颜色增强做到 image caption 那种夸张的色相饱和度变换上去因为多光谱数据的通道有明确物理含义红色波段和近红外波段的比值就对应植被指数你一个随机的色相变换会把植被和裸地的光谱界限抹掉直接导致模型错分。训练集和验证集划分的时候不能随机划分到像素级而应该按“图幅”划分。怎么理解假设你有 5 景不同日期的 Landsat 影像把 4 景的所有 patch 拿来做训练留 1 景的 patch 做验证集。这么做是因为相邻像素间的强空间自相关会让随机划分的验证集虚高——打乱的像素里验证集里跟训练集紧挨着的样本几乎长得差不多模型“背题”都能拿高分。如果真的只有一景影像那就按空间位置把整景切成左右两半左半边训练、右半边验证千万别把同一个 patch 内的像素拆开。4.2 训练脚本中的核心参数说明训练脚本的核心配置我建议写在一个单独的 YAML 或者 dict 里而不要到处散落。以下是一个地表分类任务在实际操作中被验证过多次的参数组合train_cfg { epochs: 80, batch_size: 16, learning_rate: 1e-4, weight_decay: 1e-5, lr_scheduler: StepLR, lr_step_size: 30, lr_gamma: 0.1, optimizer: AdamW, loss: WeightedCrossEntropy, }batch_size这个参数看起来无害实际上一旦设大了显存直接溢出设小了则 BN 层统计的均值和方差波动很大模型不收敛。16 这个值在 256x256 输入下对应 6 G 显存的显卡刚好跑得动。learning_rate用1e-4而不是默认的1e-3因为遥感影像类别较多每个类别的样本量差异大过大的学习率会让梯度更新方向偏向大类别模型陷入局部最优解。AdamW相比传统的 Adam 多了一项解耦的参数衰减在训练后期抑制权重振荡实测能在验证集上多出 0.5 到 1 个百分点的 Kappa。配合StepLR在 30 个 epoch 和 60 个 epoch 时各把学习率降到原来的 1/10这样前期快速收敛、后期精细调整避免在最优解附近打转。4.3 训练过程监控与断点续训一句话训练监控只盯训练集 loss 和验证集 mIoU 两条曲线。每跑完一个 epoch把训练 loss 和验证集上的平均交并比打出来如果训练 loss 一直下降但验证集 mIoU 卡住不动那说明模型过拟合了应该把数据增强强度加上去或者降低模型容量。标准的训练循环长这样for epoch in range(train_cfg[epochs]): model.train() train_loss 0.0 for batch_x, batch_y in train_loader: batch_x, batch_y batch_x.to(device), batch_y.long().to(device) optimizer.zero_grad() outputs model(batch_x) loss criterion(outputs, batch_y, class_weight) loss.backward() optimizer.step() train_loss loss.item() # 验证阶段 model.eval() iou_score compute_iou(model, val_loader, device) print(fEpoch {epoch1}/{epochs} | Loss: {train_loss:.4f} | mIoU: {iou_score:.4f}) # 保存最高分模型 if iou_score best_iou: best_iou iou_score torch.save({state_dict: model.state_dict(), best_iou: best_iou}, best_model.pth)训练过程中有几个必须注意的细节。第一计算损失用WeightedCrossEntropy时class_weight不能一次性传给损失函数就完事在每一轮迭代中它会被重复计算梯度所以要把权重移到 GPU 上并且设置requires_gradFalse不然反向传播会把梯度积累到权重上造成数值溢出。第二验证集评估要放在torch.no_grad()上下文里验证时模型不更新梯度如果不加这一步BatchNorm 层里的 running_mean 会继续被更新导致模型在验证时的特征分布被污染。第三断点续训确实是个后悔药必须在每个 epoch 结束时同时保存模型权重和优化器状态否则中间停电或显卡崩溃你前面的训练就白熬了保存优化器状态可以保证你恢复训练时的学习率调度位置和 AdamW 的动量信息不丢。5. Landsat 分类模型落地避坑Kappa 系数、边界效应与类别不平衡5.1 现象模型把所有像素都预测成背景类这是新手遇到最多的症状跑完一景影像出图的结果基本是一片黑属于典型的“背景分类器翻车”。你去看验证集的 OA 指标可能高达 0.92但 Kappa 系数不到 0.3打开混淆矩阵预测类别集中在第 0 类而真实标签分布在 1、2、3、4 类上。原因是训练样本中 0 类背景占比太高可能占到 70% 以上而损失函数用的是普通的交叉熵。优化器发现了偷懒的捷径——把所有像素全判成背景总损失已经不高了。解决办法分两条路走第一条是在采样阶段就进行类别比例平衡强制每个批次里每个类别的 patch 数量平均第二条是在损失函数里加上类别权重给少数类别大梯度。但要注意权重不能设置到超过 10否则少数类别的梯度爆炸模型极端地全预测成少数类。5.2 现象推理速度极慢GPU 利用率不到 30%在跑预测脚本的时候GPU 利用率老是掉到 5% 到 30% 之间整景影像跑下来可能要好几个小时。这种问题的根源多半不在模型本身而在于你写推理脚本时把batch_size1。单张影像切块后每个 patch 单独过 GPU连续的矩阵运算之间插入了大量 CPU 的数据读取和预处理等待GPU 大部分时间在“等饭来”。血泪经验是先把数据集的 patch 都批量预读到内存再一次性将多个 patch 组成一个批次同时送入 GPU。具体做法就是构造一个DataLoader把全部 patch 加载进来batch_size调到 16在合成概率图之后再逐 patch 写回。你会发现推理速度直接快一个数量级。同时注意在推理循环外面把model.eval()和torch.no_grad()打开不要小看这一步它能省掉显存中保存中间激活值的开销尤其在你用 6G 小显存显卡的时候不开no_grad很容易在推理到一半时直接 OOM 掉线。5.3 现象预测图边界有严重的块状拼接痕迹拼接缝的成因是模型在 patch 边缘位置能看到的视野比较窄卷积核计算边界时缺少足够的上下文导致边缘预测不稳定。如果你用 3x3 卷积核网络层级越深有效感受野对应到原图上就越大但边界上始终有一圈特征被截断。解决思路是我在第 2 章推理代码里的重叠拼接方案切块的时候让 stride 小于 patch_size让相邻区域重叠 40 到 50 个像素推理完成后对重叠部分的概率软投票取平均。另一个参数指标是如果整景影像地物类别太多patch_size 可改大一些从 256 改成 512这样能缓解边界上下文中“建筑物只看到左侧没看到右侧”的问题。代价是 GPU 显存占用变成原来的 4 倍建议 12G 显存起步。5.4 现象Kappa 系数高但城区边缘地块轮廓毛刺严重这是因为模型已经学习了多数类别的空间特征但对边缘像素的敏感性不足。传统交叉熵只逐像素计算损失完全没有考虑相邻像素的标签一致性。在后处理步骤加一个多数滤波对分类结果做 3x3 窗口的众数平滑能肉眼修正一部分毛刺但也会牺牲掉狭窄的地物例如道路。比后处理更有效的是在训练时加一个边界损失项把标签图像用拉普拉斯算子卷积出边界权重图对于离边界越近的像素给更高的损失权重。实测这一步是让 Kappa 从 0.82 提升到 0.87 的核心操作属于不用换网络结构就见效最快的调参方案。6. 把模型用于大范围制图的进阶验证技巧当模型训练好、推理也通了你还要做最后一件事不要急于把结果直接交付给甲方先用验证数据集做一次严格“质检”尤其是当你把模型拿来做大范围制图时城市 A 训练的模型拿到城市 B 去推理精度可能自己就降几个点。我习惯在每次模型发布之前跑一份标准的验证报告里面包含三大块内容第一是全局指标——OA、Kappa、每类的 IoU 和 F1-score第二是按区块统计的精度混淆矩阵这样能发现模型在哪个空间区域表现失效第三是概率图叠加底图的人工抽检——把模型预测概率介于 0.4 到 0.6 之间的“犹豫区域”叠加到原始影像上人工看一眼这些地方到底是标签错了还是模型错了。你往往会在这种抽检中发现标签图本身存在大量噪声自己的标注和卫星影像对不上。写验证脚本的时候我会直接复用训练时的验证数据集抽取逻辑但注意不要用训练时的数据增强模块改而只用原始影像切块。用 sklearn 里的classification_report可以直接打印出每类的 F1-score然后单独计算 Kappa。建议在报告里把每个类别的混淆矩阵逐行打印因为遥感制图最怕的就是把建筑错分成裸地、把水体错分成阴影。多光谱影像里的阴影和水体在光谱特征上确实高度混淆如果能接受这种错分你需要在交付报告里坦诚地写清楚——这也是你专业性的体现。还有一个习惯训练完成后把模型在多个时相、不同季节的影像上各跑一遍统计精度方差。Landsat 影像随季节变化极大夏天植被是绿色冬天裸土可能和建筑屋顶混在一起。只在一个时相上验证通过的模型必须带着这个风险点去和用户沟通。希望这些踩坑和调参的经验能帮到你照着这套流程至少可以少走一大半弯路。评估无误后把模型权重、验证报告和概率图一起归档交付出去这才是整套 Python 遥感分类方案最完整的落地方式。本文还有配套的精品资源点击获取