简介这是一套基于CNN深度学习的遥感Landsat影像地物分类算法Python完整源码面向遥感、深度学习和计算机视觉方向的在校生、算法入门者及企业研发人员既可用于课程设计、毕业设计也可作为工程项目的初期演示。代码经测试运行成功流程完整覆盖训练样本制作、模型训练、新数据预测三个关键环节并附有预训练H5模型、TIF测试影像、XML/TFW辅助信息及Markdown运行说明能够帮助读者快速搭建实验环境、理解CNN地物分类的基本链路并直接查看预测效果。压缩包共10个文件主要包含Python脚本、遥感影像、模型权重和说明文档整体约14.88MB体积轻量、目录清晰方便本地复现与二次开发。目前已有470人学习浏览适合需要一套可运行基准来完成遥感图像分类实战或毕业设计选题的同学参考。1. “基于CNN深度学习的遥感landsat影像地物分类算法python完整源码”到底缺什么经常有人往工作群丢一个压缩包基于CNN深度学习的遥感landsat影像地物分类算法python完整源码.zip。下载的人不少真正跑通并出图的屈指可数。原因不是 CNN 多玄学而是前置链路更吃功夫Landsat 波段怎么排、标签坐标怎么对齐、类别不平衡怎么处理哪个环节缺了训练脚本再完整也等于白给。地物分类本身很直观就是把 30 米像元识别成水体、耕地、林地和不透水面难点在把原始 TIF 切成模型能吃的 patch。这篇笔记按 Python 工程落地的流程拆解数据准备、模型训练、推理出图和避坑适合已经装好 Python、想把手头 Landsat 影像真正变成分类图的工程师。2. 先立技术骨架Landsat 波段、标签与 CNN 模型的选型理由先澄清一个容易跑偏的地方地物分类不是把目标框出来而是给每个像元一个类别CNN 在这类任务里之所以管用不是因为“深度学习”听起来强而是卷积核天然合适。Landsat 像元本身只有光谱和邻域结构两种信息CNN 的 3×3、5×5 卷积正好同时编码“像元光谱特征”和“同质像元的空间纹理”。如果换用全连接层模型必须把整景影像展开后才能输入参数量和训练样本量都会跟着膨胀落地时你会怀疑人生。相比之下CNN 把 64×64×6 的小 patch 作为输入参数量小普通办公电脑也能跑这也是这类完整源码项目默认选 CNN 的真实原因。2.1 Landsat 8/9 哪些波段值得进模型Landsat 8 的 OLI 传感器给了 9 个多光谱波段看起来选择很多实际不能全塞给模型。我一般取六个波段B2 蓝、B3 绿、B4 红、B5 近红外、B6 短波红外 1、B7 短波红外 2。这六段把植被、水体和建筑裸土的关键光谱差异都覆盖了水面吸收近红外和短波红外的特征混凝土屋顶在短波红外上的响应和植被正好相反六个波段足够让一个浅层 CNN 学到这些可分性特征。波段组合常见做法适用场景B2B3B4真彩色输入主要用于目视检查影像不建议直接训练信息量偏少B4B5B6假彩色红外植被/水体的可分离性最好训练效果明显好过真彩色B2-B7 六通道标准分类输入城市与自然地表混合地物推荐作为模型起点B1/B9 卷云波段可选附加高云和海岸影像可能有帮助常规陆地任务建议扔掉另一个前提是产品级别。我一般直接用 Landsat Collection 2 Level-2 的表面反射率产品文件名里带_SR后缀因为 Level-1 的辐射值受大气影响很大没有做大气校正前直接训练模型学到的全是那一景当天的天气条件换一景就翻车。表面反射率产品省掉大气校正这道工序训练前只需按产品说明做尺度换算再按通道做归一化数值稳定性会好很多。2.2 选 patch CNN 还是 U-Net先看手里标签长什么样论文里常见的是全卷积网络做密集预测U-Net 也好DeepLab 也好但真实遥感项目里多数情况是从 patch CNN 起步的。我的判断标准是如果标签栅格来自目视解译或者从矢量面转栅格边界本身就比较粗一个 patch 记一个主导类别训练起来最稳patch CNN 在数据构造阶段就绕开了标签边界误差大的问题。U-Net 的优势是逐像元边界更干净损失是标签要足够精细。如果手里只有一个大面标签硬上 U-Net 反而会把边界噪声学进去训出来的模型在真实边缘上看着平滑实际上没有推理能力。如果手里的标签是高分影像手工逐像元勾出来的再换 U-Net模型骨架可以复用。无论选哪种特征输入都是同一套六通道真正影响进度的不是网络名字而是样本标签怎么构造。2.3 一份完整的 Python 源码包至少要凑齐五个文件拿到任何“完整源码”压缩包我做的第一件事不是看模型结构而是找五个文件数据预处理、数据集加载、模型定义、训练入口、推理输出。这五个缺一个要么跑不动要么跑通后无法换数据复现。常见尴尬是 model 定义写得漂亮训练脚本却硬编码了一张影像路径别人换数据就起不来。文件角色应该包含的内容验证标准数据预处理波段读取、重采样、标签对齐、切片输出的 npz/tfr 文件能单独检查数据集加载训练/验证划分、数据增强、归一化只改路径不改代码就能加载模型定义patch CNN 或 U-Net输入输出形状明确model.summary() 能打印且和数据集匹配训练入口超参、回调、日志、模型保存修改配置后能重新训练推理脚本滑窗推理、投票融合、GeoTIFF 写出输入整景影像输出带投影坐标的 tif如果压缩包按这个结构分清楚了哪怕没有文档我也愿意花时间去解。如果只有一个 ipynb 挤着全部逻辑建议直接换一套更清晰的骨架后面调数据、调超参、排查问题没有边界线改一步崩三处。3. 把 Landsat 原始影像切成 CNN 训练集重采样、标签对齐和滑窗切片Landsat 单景大概 7000×7000 到 8000×8000 像元不可能直接塞进 CNN。常见做法是把整景切成 64×64 或 128×128 的 patch。切之前有两件事必须先做统一六波段的像元网格把标签和影像的坐标对到同一个格子上。这两件事做错一步训练时看着每个 patch 数值正常推理时 patch 对应的地物位置其实已经错位了。3.1 六通道波段读入与统一重采样一个可复用的读图函数Landsat 产品里同为 30 m 分辨率的多光谱波段整体对齐良好但覆盖范围偶尔差几个像元最好统一用out_shape重采样。下面这个函数用 rasterio 把六个波段读成一个(h, w, 6)的浮点栈相当于做了一次内存级重采样。import glob import numpy as np import rasterio from rasterio.enums import Resampling BAND_TAGS { blue: B2, green: B3, red: B4, nir: B5, swir1: B6, swir2: B7, } def load_landsat_stack(src_dir, height256, width256): 把 Landsat 六个波段读成 (height, width, 6) 数组。 波段顺序固定为 B2/B3/B4/B5/B6/B7模型输入必须和这里保持一致。 bands [] for tag in BAND_TAGS.values(): matches glob.glob(f{src_dir}/*_{tag}_SR.TIF) if not matches: raise FileNotFoundError(f缺少波段 {tag} 的文件) with rasterio.open(matches[0]) as src: arr src.read( 1, out_shape(height, width), resamplingResampling.cubic, maskedFalse ).astype(np.float32) bands.append(arr) return np.stack(bands, axis-1) image, profile load_landsat_stack(./data/LC08_xxxx, height512, width512) print(image.shape, image.dtype, profile[crs])这段代码的取值逻辑out_shape写在read()里让 rasterio 在读盘过程中直接重采样不需要先把整景载进内存再缩适合单景几 GB 的场景。cubic适合表面反射率这种连续量但不能用于标签标签必须用nearest否则会插值出第 2.7 类这种不存在的类别。maskedFalse是刻意去掉掩膜nodata 先保留原始数值后续归一化时统一替换避免MaskedArray在训练时带来形状和类型上的隐患。3.2 标签栅格和影像对不齐先统一坐标系再裁 patch从整景影像落到带标签 patch最容易出的问题不是模型而是标签与影像差半像素。切片时正确的做法是用影像的 transform 推导窗口再用同一个窗口去标签文件里取数据而不是各自读整幅再凭直觉剪。如果标签和影像的 CRS 不一致标签常见是矢量转栅格后的投影坐标而影像本身是地理坐标先做一次重投影。import rasterio from rasterio.enums import Resampling from rasterio.warp import reproject import numpy as np def align_label_to_image(label_path, ref_src): 把标签栅格重投影/重采样到与参考影像完全相同的网格。 dst_array np.empty((ref_src.height, ref_src.width), dtypenp.uint8) with rasterio.open(label_path) as lab: reproject( sourcelab.read(1), destinationdst_array, src_transformlab.transform, src_crslab.crs, dst_transformref_src.transform, dst_crsref_src.crs, resamplingResampling.nearest ) return dst_array这里nearest是硬约束类别标签是离散值任何插值都会制造伪类别。重投影后立刻检查dst_array.shape是否等于ref_src.height, ref_src.width这样能发现因为边缘对齐产生的差几个像元问题。之后切片就可以用同一组行列范围同时取影像和标签坐标误差被控制在重采样允许的范围内。ref rasterio.open(./data/LC08_xxxx_B4_SR.TIF) label_aligned align_label_to_image(label.tif, ref) patch_size 64 for row in range(0, ref.height - patch_size 1, patch_size): for col in range(0, ref.width - patch_size 1, patch_size): img_patch image[row:row patch_size, col:col patch_size, :] lab_patch label_aligned[row:row patch_size, col:col patch_size] # 这里把 lab_patch 中 0 视为 nodata/背景不参与训练这段循环是切片主逻辑如果背景占比太高就跳过否则这个训练集会掺入大量无意义像元。切片后的lab_patch如果原有类别编码是 1 到 6进入模型前记得减一变成 softmax 的索引 0 到 50 这个类别值保留给 nodata不让它参与损失计算。3.3 滑窗步长、样本筛选与数据增强参数patch 的选取不能按顺序机械切。整景影像有相当一部分是背景、nodata 或云影模型收大量空白 patch分类能力会被稀释。我一般做两层筛选先统计 patch 内有效像元占比达到 90% 才保留再看类别分布对稀有小类别做重复滑窗让水体、不透水面这类低频类别至少积累一千个以上的 patch。参数/做法常用默认值说明patch size64×64覆盖约 1.9 km×1.9 km对地物纹理的表达足够训练切窗步长32前后 patch 重叠一半样本量更大推理切窗步长64 或 32推理重叠用于多数投票防止块状拼接背景剔除阈值90%有效像元占比不足 90% 的 patch 直接丢弃数据增强旋转 90/180/270、水平翻转对语义类别安全不会破坏光谱物理含义归一化各通道独立 mean/std 或分位数裁剪同景数据统一记录均值和方差推理时复用要特别提醒遥感分类的数据增强不要加入随机亮度、随机对比度、HSV 抖动这一类常用 CV 增强。Landsat 表面反射率本身有定量含义亮度抖动会让模型把同一块水体在不同光照下学成两个类别。光谱一致性比样本数量更重要。4. 训练 CNN 的 Python 骨架类别权重、模型定义与超参设置前面处理完数据训练环节有三件事要做算类别权重、定义 CNN、配训练参数。这套骨架只依赖 TensorFlow/Keras可以平移到自己的工程里如果模型输入形状和第 3 节输出的数组不一致先检查六通道顺序别急着怀疑网络结构。4.1 类别不均衡先统计标签再训练Landsat 数据里耕地、林地面积占比大水体和不透水面都是小类别。随机切片后水体可能只占 2%模型把一切预测成耕地整体准确率依然有 90%因为水体贡献的 loss 太小。常见做法是两路并进采样阶段让稀少类别多出 patch训练阶段给稀少类别更高的类别权重。下面这段代码统计标签栅格返回类别权重。def compute_class_weight(label_paths, n_classes6): freq np.zeros(n_classes, dtypenp.float64) for p in label_paths: with rasterio.open(p) as src: for _, win in src.block_windows(1): data src.read(1, windowwin) vals, counts np.unique(data, return_countsTrue) for v, c in zip(vals, counts): if 0 int(v) n_classes: freq[int(v) - 1] c freq freq / np.clip(freq.sum(), 1e-6, None) weights 1.0 / np.clip(freq, 1e-6, 0.5) return weights / weights.sum() class_weight compute_class_weight([train_label.tif, val_label.tif]) print(class_weight)这里的约定是标签文件取值 1 到 60 是 nodata/背景所以统计时写int(v) - 1。block_windows(1)让大影像按块逐个读避免整景一次性加载进内存。极端权重会被截断比如某个类占比 0.1% 时权重会很大直接放大到几十倍容易让模型把噪声也学进该类所以我用np.clip(..., 0.5)把最大权重限制在 2 倍附近。4.2 一个 64×64×6 的轻量 CNN 骨架模型定义不追求深关键是输入通道要和前面一致。这里给一个三卷积层加全局平均池化的轻量 CNN适合六通道小 patch 的分类。import tensorflow as tf from tensorflow.keras import layers def build_cnn(input_shape(64, 64, 6), num_classes6): inputs layers.Input(shapeinput_shape) x layers.Conv2D(32, 3, paddingsame, activationrelu)(inputs) x layers.BatchNormalization()(x) x layers.MaxPooling2D(2)(x) x layers.Conv2D(64, 3, paddingsame, activationrelu)(x) x layers.BatchNormalization()(x) x layers.MaxPooling2D(2)(x) x layers.Conv2D(128, 3, paddingsame, activationrelu)(x) x layers.BatchNormalization()(x) x layers.GlobalAveragePooling2D()(x) x layers.Dropout(0.3)(x) outputs layers.Dense(num_classes, activationsoftmax)(x) model tf.keras.Model(inputs, outputs) return model model build_cnn() model.summary()这段结构里GlobalAveragePooling2D把最后一张特征图直接平均比Flatten Dense的参数少一个数量级对这个规模的 patch 更稳。BatchNormalization在遥感输入上很有用能把不同时相、不同大气条件下的通道分布拉回统一尺度对跨景泛化有实际帮助。如果之后要换 U-Net前面的卷积块和归一化方式都可以继续沿用换的只是末端的输出结构。4.3 训练参数batch、learning rate、早停一步配齐训练脚本不复杂但参数选择要克制。遥感 patch 分类通常不需要上千 epoch更怕的是一路过拟合到单景影像。model.compile( optimizertf.keras.optimizers.Adam(learning_rate1e-3), losssparse_categorical_crossentropy, metrics[accuracy], ) callbacks [ tf.keras.callbacks.ReduceLROnPlateau( monitorval_loss, patience5, factor0.5, min_lr1e-6 ), tf.keras.callbacks.EarlyStopping( monitorval_loss, patience10, restore_best_weightsTrue ), ] history model.fit( train_dataset, validation_dataval_dataset, epochs60, class_weightclass_weight, callbackscallbacks, batch_size32, )这套参数的含义Adam 以 1e-3 起步验证 loss 连续 5 轮不降就把学习率减半最多降到 1e-6早停在 10 轮内没有提升就停并恢复历史最优权重。class_weight就是 4.1 节算出的权重数组。如果你已经把数据包装成了tf.data.Dataset且做了 batch就去掉batch_size传入原始 NumPy 数组时才需要保留。单张消费级 GPU 上64×64 patch、batch 32 基本没有压力如果把 patch 增到 128batch 要先降到 16 或 8。5. CNN 做地物分类最常踩的五个坑现象、原因与排查办法这一章集中写最常翻车的五个现场。每一条我都反复遇到过不少问题从精度指标上完全看不出来只有打开分类图那一刻才暴露。5.1 现象一训练集精度 95%换一景影像就崩现象训练验证精度很高模型换到相邻景或者不同月份的影像分类效果断崖式下跌。原因样本只来自一景、一个月或者一种天气条件。Landsat 影像里同一块耕地7 月偏绿10 月偏枯模型学到的可能是当天的大气、太阳高度角和土壤湿度而不是地物本身的属性。这是遥感分类最常见的“数据集漂移”。解决训练样本至少要覆盖不同季节、不同轨道号的三景以上暂时凑不齐多时相就在每景内部做通道归一化至少把绝对辐射水平压平。更重要的是一开始就留出一景完全独立的影像做测试不要用同景切片做验证。5.2 现象二水体消失了分类图里只剩耕地和裸土现象混淆矩阵整体精度 0.93打开分类图找不到那条河。原因水体占整景比例太低随机滑窗时水体 patch 数量不足或者类别权重算好了但没有在model.fit里真正传进去。水体作为稀少类别在损失函数里被耕地和森林稀释掉。解决切片阶段针对水体矢量做定向采样把每条主要水体的中心点附近的窗口额外滑一遍再把 4.1 节的class_weight传给训练接口。如果水体还是消失把水体设成必保留类训练集里强制保持至少 10% 的水体 patch。5.3 现象三一训练 GPU 显存就爆batch 降小了也不稳现象训练刚开始就报 CUDA out of memory把 batch 从 32 降到 8 仍然中断。原因patch 本身不大但数据 loader 可能用np.stack把全部 patch 一次性读进内存或者验证集是完整张量没有走流式。显存不是被模型吃掉的是被输入预处理吃掉的。解决把输入封装成tf.data.Dataset.from_tensor_slices(...).batch(...)的形式让数据按 batch 流式加载不要在训练前把所有 patch 都堆进一个 numpy 数组。显存仍然不够时把 patch 降到 48×48 或 batch 降到 8一般都能解决。5.4 现象四推理出图呈豆腐块边界锯齿非常重现象分类图能看出明显的 64×64 格子痕迹类别边界呈锯齿状。原因推理时切窗步长等于 patch 尺寸patch 之间没有重叠空间相邻的两个 patch 由不同的随机初始化权重给出不同判别自然拼出不连续。解决推理时把步长设为 patch 的一半例如 patch 64、stride 32重叠部分按 softmax 概率累加再取 argmax块状感会明显减弱。这个方案就是下一章给出的完整推理脚本如果应用场景对边界要求极高最终还是要换 U-Net。5.5 现象五换一台电脑rasterio 连 tif 都读不了现象新环境里import rasterio后调用rasterio.open报错或者提示找不到libgdal。原因用 pip 直接装的 rasterio 依赖系统里的 GDAL 动态库版本对不上就会失败换 Python 版本后有时连 wheel 都不匹配。环境问题是最不值钱的消耗但确实最磨人。解决用 conda 固定一个环境例如conda create -n remote python3.10再用conda install -c conda-forge gdal rasterio让 GDAL 和 rasterio 同源安装这样环境迁移完全可复现。装完先跑rio info x.tif自测不要等训练到一半才暴露。6. 整景影像推理出图用重叠滑窗把模型落成 GeoTIFF模型训练完最后要把整景影像落成一张带坐标的 GeoTIFF。下面的推理函数用重叠滑窗累加概率stride 取 patch 的一半正好解决第 5.4 节的豆腐块问题。def predict_full_scene(model, image_stack, patch_size64, stride32, num_classes6): h, w, _ image_stack.shape prob np.zeros((h, w, num_classes), dtypenp.float32) count np.zeros((h, w, 1), dtypenp.float32) for row in range(0, h - patch_size 1, stride): for col in range(0, w - patch_size 1, stride): patch image_stack[row:row patch_size, col:col patch_size] patch patch[None, ...] p model.predict(patch, verbose0)[0] prob[row:row patch_size, col:col patch_size] p count[row:row patch_size, col:col patch_size] 1 prob / np.maximum(count, 1.0) return prob.argmax(axis-1).astype(np.uint8)prob存每个像元每个类别的累积概率count记录每个像元被多少 patch 覆盖最后相除得到平均概率。这样重叠区域内不同位置预测的不一致性会被投票平滑掉。把输出写回 GeoTIFF 时沿用原始影像的 profile只改数据类型和波段数。profile dict(image_profile) profile.update(dtypeuint8, count1, compresslzw) with rasterio.open(cls_result.tif, w, **profile) as dst: dst.write(result, 1)写完后我会先不调任何指标直接把结果叠加在原始影像上看两块地方水体有没有断线城镇和裸土有没有混在一起。如果空间上明显不合理再回切片阶段补样本而不是急着调超参。这个方向值得坚持Landsat 免费且长时序把源码链路跑通以后换 U-Net 或做多时相模型都只改网络那一层。希望这篇笔记能给你省下那些我当初白耗的调试时间。本文还有配套的精品资源点击获取