基于CNN的Landsat遥感影像地物分类Python源码实战与避坑

简介基于CNN深度学习的遥感Landsat影像地物分类Python源码与项目说明面向计算机、人工智能、地理信息等相关专业学生和从业者适合课程设计、毕业设计以及深度学习入门实战练习。整个压缩包共10个文件约14.89MB主要包含3个Python脚本分别用于影像样本切片生成、CNN七分类模型训练和新数据地物分类预测附有已训练好的h5模型权重文件、示例遥感tif影像以及配套的坐标参考xml/tfw辅助文件另有Markdown格式的项目说明文档目录结构清晰方便按步骤复现实验。代码均经过测试运行成功可直接下载使用能帮助读者完整理解从遥感影像预处理、样本制作、模型训练到最终地物分类的深度学习流程也可作为课程作业和算法研究的基础参照。目前该项目已有970人学习对于希望快速上手遥感图像分类的开发者具有较高参考价值。1. 基于CNN深度学习的遥感Landsat影像地物分类这份源码到底能解决什么一景Landsat影像动辄几百万个像元想分清水体、农田、林地、建设用地靠单波段阈值和最大似然很快翻车换成CNN深度学习做遥感影像地物分类又卡在样本生成、模型训练和坐标对齐上。这份基于CNN深度学习的遥感Landsat影像地物分类Python源码包把影像切片、模型训练、新数据预测整条链路一起打包了还带了一个训练好的7类模型适合做课程设计、毕业设计或者想快速跑通深度学习遥感分类流程的入门者直接从example.tif复现。它解决的核心问题就一句话让一个没有遥感深度学习基础的人也能在半天内得到一张像样的地物分类图。2. 为什么是7类、3×3切片地物分类任务的输入设计与波段选择2.1 逐像素分类如何变成CNN的3×3输入Landsat影像中每个像元都有一个多光谱向量传统分类器只看这个向量本身属于“单点光谱分类”。问题是地物在30米分辨率下经常存在混合像元比如田地边缘、建筑与树木交错的地方单点光谱往往落在两个类别之间分类结果会出现大量椒盐噪声。解法是把单个像元扩展成一个小邻域让模型同时看到中心像元和周围像元的空间关系这就是切片chip的由来。项目里的3×3切片就是取目标像元为中心的3×3窗口一共9个像元全部波段叠加在一起作为模型输入。对Landsat 8 OLI的常见处理是使用6个波段那么单个样本的shape就是(3, 3, 6)。3×3窗口在30米分辨率下对应90米×90米的实际地面范围足够覆盖农田、裸土、水体这种面积较大的地物又不会把太多不同地物混进一个窗口。相比直接用单个像元这种带邻域的输入让CNN第一层就能学到纹理和边界特征。我一般会把这种输入理解成“微型图像分类”每个样本是一张3×3的小图通道数等于波段数。训练时模型在一个小图上做7分类预测时对整幅影像的每个像元都取一遍窗口再逐个分类。这个过程在1_createImageChips.py和3_predictNewData.py里分别对应样本生成和推理预测后面会逐步拆。2.2 波段组合与7类地物体系的选用逻辑Landsat 8/9的OLI传感器从可见光到短波红外有多个波段做地物分类时不会把全部波段都塞进去而是要选对地物光谱差异敏感的波段。常见组合是蓝(B2)、绿(B3)、红(B4)、近红外(B5)、短波红外1(B6)、短波红外2(B7)这6个波段。水体在近红外和短波红外强烈吸收植被在近红外高反射裸土和建设用地在短波红外有明显差别这套组合能把大多数地表覆盖类型区分开。项目里模型名是CNN_7class_3by3.h5明确就是7类、3×3窗口。7类的分类体系很常规我按常见做法理解为水体、农田、林地、草地、裸土、建设用地、阴影。阴影单独成一类是因为山体和高层建筑阴影在可见光波段看起来像暗色水体不单独分出来会导致水体面积虚高。如果你有自己的分类体系只要在2_trainModel.py里把标签映射表改一下即可网络结构不用大动。信息量上6个波段加上3×3邻域每个样本是9×654个数值。这个维度对CNN来说其实很小所以网络不需要很深很宽重点是样本质量和类别均衡。项目里没有额外计算NDVI等植被指数作为输入波段如果你需要区分更复杂的地物可以在此基础上把NDVI、NDWI拼进去代价是重训模型这个在第6章会说。2.3 源码包文件清单与整体流程拿到压缩包后先看文件清单。下面这几个文件和你的操作顺序直接相关文件作用使用顺序example.tif示例Landsat影像做切片和训练的底图第一步1_createImageChips.py从影像生成3×3训练切片第二步2_trainModel.py训练CNN模型并保存第三步CNN_7class_3by3.h5训练好的7类模型权重可直接使用new_class.tif待分类的新影像预测输入new_class.tif.aux.xmlGDAL辅助元数据文件自动生成无需管new_class.tfwTIF坐标配准文件辅助地理参考new_class.tif.xml影像XML元数据描述辅助信息3_predictNewData.py载入模型对new_class.tif分类最后一步README.md项目说明与运行步骤先读整个流程是线性的先对example.tif做切片得到训练样本集再用2_trainModel.py训练CNN得到CNN_7class_3by3.h5最后用3_predictNewData.py加载这个h5对new_class.tif逐像元预测输出一张分类图。如果你不想重训跳过前两步直接用现成的h5也可以但前提是待预测影像的波段顺序和数量与训练时一致这个坑在第5章专门讲。3. 从影像到样本集1_createImageChips.py切片生成实战3.1 滑动窗口生成3×3样本的核心逻辑1_createImageChips.py是整个流程的地基它做的事可以概括为读影像、读标签、滑动取窗口、存样本。项目压缩包里没有看到单独的标签文件这是遥感深度学习最常见的现实情况——标签通常需要自己准备。常见做法是用GlobeLand30、FROM-GLC这类公开土地覆盖产品作为标签或者自己在GIS里人工勾几个分类多边形栅格化。我根据这个项目的工作流给出一个可复现的切片脚本核心逻辑你拿到源码后对照着看更容易理解# 1_createImageChips.py 核心逻辑用rasterio读取Landsat影像和标签栅格 # 以每个有效标签像元为中心取3×3邻域作为CNN样本 import rasterio import numpy as np INPUT_TIF example.tif # 训练影像Landsat多波段 LABEL_TIF label.tif # 标签栅格类别编码1-7 CHIP_SIZE 3 # 窗口边长 PAD CHIP_SIZE // 2 # 窗口半径 with rasterio.open(INPUT_TIF) as src: img src.read() # shape: (bands, rows, cols) n_bands, rows, cols img.shape with rasterio.open(LABEL_TIF) as ds: label ds.read(1) # 读取第一波段标签 chips, labels [], [] for r in range(PAD, rows - PAD): for c in range(PAD, cols - PAD): cls label[r, c] if cls in [0, 255]: # 0和255通常是nodata跳过 continue chip img[:, r-PAD:rPAD1, c-PAD:cPAD1] # 取3×3多波段窗口 chip chip.astype(float32) / 10000.0 # 地表反射率缩放 chips.append(chip) labels.append(cls) chips np.array(chips) # (N, bands, 3, 3) labels np.array(labels) np.savez_compressed(chips.npz, chipschips, labelslabels) print(f生成样本 {chips.shape[0]} 个每样本 shape {chips.shape[1:]})这段代码里有两个参数值得说明。第一个是10000.0缩放Landsat Collection 2的地表反射率产品为了存成整数把0到1的反射率乘以了10000所以读出来是uint16整数训练前必须除以10000恢复到0到1区间否则数值太大CNN训练时会梯度爆炸。第二个是窗口边界设计从第1行循环到rows-1实际上是舍弃了影像最外一圈像元因为边缘没有完整邻域直接索引会得到形状不匹配的切片。保存下来的chips.npz里有样本和标签两个数组。chips的shape默认是(N, bands, 3, 3)而TensorFlow的卷积层默认使用通道在最后的布局所以在2_trainModel.py里加载后要转成(N, 3, 3, bands)或者直接用np.moveaxis调整轴顺序。这是新手最容易忽略的一步轴顺序错了训练时第一步就会报shape不匹配。3.2 样本均衡、随机抽样与标准化参数怎么设直接把全图所有像元都切下来是灾难。一景Landsat影像大约8000×8000像元即使只有一半是有效标签也会生成上亿个3×3×6样本内存先爆。我一般会在切片脚本里加两个机制随机抽样和类别配额。随机抽样的意思是如果标签有效像元太多就按比例随机选取一部分比如每类最多保留2万个样本。类别配额更重要因为地物类别天然不平衡建设用地可能只有几万个像元农田却有几百万个如果不限制CNN会被农田类主导其他类别学不出来。给一个实用的均衡采样片段直接放在3.1代码的循环结束之后# 类别配额每类最多保留MAX_PER_CLASS个样本超出部分随机丢弃 from collections import Counter import numpy as np MAX_PER_CLASS 20000 labels np.array(labels) # labels和chips一一对应这里直接统计并采样 counter Counter(labels) print(原始类别分布:, dict(counter)) keep_indices [] for cls, cnt in counter.items(): idx np.where(labels cls)[0] if cnt MAX_PER_CLASS: idx np.random.choice(idx, MAX_PER_CLASS, replaceFalse) keep_indices.append(idx) keep np.concatenate(keep_indices) chips chips[keep] labels labels[keep] print(均衡后类别分布:, dict(Counter(labels)))这段代码的逻辑并不复杂逐类别找到所有样本下标超过配额就随机抽取一部分最后把保留的下标拼起来做索引。MAX_PER_CLASS设多少要看你机器内存我常用的经验值是1万到3万样本太少模型欠拟合样本太多训练时间成倍增加。标准化参数也要保持一致。训练脚本里除以10000预测脚本里也必须除以10000这是最容易踩的不一致问题。另一个容易忽略的是波段顺序Landsat的B2到B7顺序必须固定你要是某天换了一景数据源波段顺序变了模型预测的准确率会直接崩掉这个在第5章里细说。3.3 切片阶段最容易出现的三个异常切片阶段有三个问题几乎每个人都会遇到。第一个是内存耗尽全图逐像元切片生成几十G的npz解决方法是像3.2那样限制每类样本量或者分块处理大影像一次只读一块子区域。第二个是array reshape错误报错信息类似“cannot reshape array of size 54 into shape (3,3,6)”原因通常是边界像元窗口不完整或者影像波段数与代码里写的不一致。第三个是标签和影像行列数对不上用ArcGIS导出的标签图如果投影不同要先重采样到同一分辨率再读。这三个问题解决起来都不难难在提前意识到。切片脚本跑之前先打印一下img.shape和label.shape确认两者行列数一致再打印类别分布确认没有奇怪的类别值比如标签图里混进了255当成有效类。这样能省下大量调试时间。4. 训练与预测2_trainModel.py 和 3_predictNewData.py 的复现路径4.1 小patch CNN的网络结构与超参数2_trainModel.py负责加载chips.npz并训练CNN。因为输入是3×3×6的小patch网络结构不需要很深传统的VGG式堆叠在这里会显得笨重。常见的做法是两层卷积加Dropout再直接接全连接分类我给出这个项目场景下说得通的网络定义# 2_trainModel.py 核心逻辑加载切片定义小patch CNN并训练 import numpy as np from tensorflow.keras.models import Sequential from tensorflow.keras.layers import Conv2D, Flatten, Dense, Dropout from tensorflow.keras.utils import to_categorical data np.load(chips.npz) X data[chips] # (N, bands, 3, 3) # 通道维移到最后一维变成 (N, 3, 3, bands) X X.transpose(0, 2, 3, 1).astype(float32) y data[labels].astype(int32) # 类别重映射到0-6 class_map {1:0, 2:1, 3:2, 4:3, 5:4, 6:5, 7:6} y np.array([class_map[v] for v in y]) y to_categorical(y, num_classes7) model Sequential([ Conv2D(32, (3, 3), activationrelu, paddingvalid, input_shape(3, 3, X.shape[-1])), Conv2D(64, (3, 3), activationrelu, paddingvalid), Flatten(), Dense(128, activationrelu), Dropout(0.5), Dense(7, activationsoftmax) ]) model.compile(optimizeradam, losscategorical_crossentropy, metrics[accuracy]) history model.fit(X, y, validation_split0.2, epochs30, batch_size64) model.save(CNN_7class_3by3.h5) print(训练完成模型已保存)这段代码里我特意没有加MaxPooling2D。原因是3×3输入经过2×2池化后直接变成1×1空间信息瞬间清零后续卷积等于在1×1的像素上做无意义计算这是很多人在小patch分类上翻车的点。3×3窗口能提供的信息量本来就有限两层3×3卷积已经能把邻域内的光谱组合特征提取得比较充分了。参数设置上batch_size 64、epochs 30、adam默认学习率对小数据集是稳妥起点。如果训练精度上不去先别调网络结构去看样本质量比如标签有没有错位、波段缩放是否一致。如果训练精度很高但验证精度低把Dropout放到0.5再增加样本量。小patch网络通常几分钟到十几分钟就能训完CPU也能跑这点比大分辨率影像分割友好得多。4.2 预测新影像3_predictNewData.py是怎么工作的预测脚本的核心是逐像元取3×3窗口喂给模型取softmax概率最大的类别写入新的TIFF。下面是核心逻辑的还原# 3_predictNewData.py 核心逻辑读取new_class.tif逐像元预测并输出分类图 import rasterio import numpy as np from tensorflow.keras.models import load_model model load_model(CNN_7class_3by3.h5) with rasterio.open(new_class.tif) as src: img src.read().astype(float32) # (bands, rows, cols) profile src.profile n_bands, rows, cols img.shape # 地表反射率缩放必须在预测时重复训练时的操作 img img / 10000.0 out np.zeros((rows, cols), dtypeuint8) PAD 1 for r in range(PAD, rows - PAD): for c in range(PAD, cols - PAD): chip img[:, r-PAD:rPAD1, c-PAD:cPAD1] # (bands, 3, 3) # 检测nodata或无效像元跳过不预测 if np.isnan(chip).any() or (chip 0).any(): continue chip chip.transpose(1, 2, 0)[np.newaxis, ...] # (1, 3, 3, bands) pred model.predict(chip, verbose0) out[r, c] np.argmax(pred[0]) 1 # 映射回类别1-7 profile.update(dtypeuint8, count1, nodata0) with rasterio.open(new_class_classified.tif, w, **profile) as dst: dst.write(out, 1) print(分类结果已写出: new_class_classified.tif)这段代码每行都有讲究。img / 10000.0必须在预测脚本里原样重复不然训练时模型见到的数值范围是0到1预测时喂进去的是0到10000的整数输出会完全错乱。chip[..., np.newaxis]加批次维是因为Keras的predict要求输入是四维即使只有一个样本也要补上batch那一维。out[r, c] np.argmax(pred[0]) 1里的1是把CNN输出的0到6索引还原成你定义的类别编码1到7。逐像元循环加上model.predict调用确实很慢一景影像可能要预测几十分钟甚至几个小时。优化方式很明确把整幅影像的patch一次性提取出来批量预测常见做法是用tf.image.extract_patches把影像重组成(N, 3, 3, bands)的数组然后model.predict(patchs, batch_size1024)一次跑完。缺点是会有大量内存开销但比单像元循环快两个数量级。4.3 复用CNN_7class_3by3.h5时先检查两件事包里已经带了训练好的CNN_7class_3by3.h5不少人拿到后直接load_model然后预测出问题才回头看。我建议复用时先检查两件事。第一是波段数如果训练时用的是6个波段待预测影像也是6波段模型才能跑通如果new_class.tif是5波段或者7波段模型加载不会报错但predict时数组shape直接不匹配。第二是坐标参考new_class.tif带有了tfw文件和aux.xml说明它是有地理参考的预测脚本用profile.update保留坐标信息输出的分类图才能和原影像准确叠加。另一点容易被忽略直接用现成h5意味着你接受它训练的7类定义和你的目标类别一致。如果项目里的7类和水体、农田、林地、草地、裸土、建设用地、阴影这套体系吻合那直接预测没问题如果分类体系不一致就不要偷懒必须用你自己的标签重新走一遍切片和训练流程。5. 避坑排查Landsat影像做地物分类的四个典型翻车场景5.1 分类结果整片黑色或大片为0现象预测输出的new_class_classified.tif打开后几乎全黑只有零星几个像元有类别。原因我用上面代码逐像元检查了chip 0就跳过如果new_class.tif的边缘背景是0值那么整幅影像的边界一圈全被跳过。更常见的原因是待预测影像中存在大量nodata像元或者影像本身的数值范围不是0到1而是原始的0到65535但预测脚本里忘了除以10000导致输入值普遍大于1模型softmax输出变得很平大量像元被分到概率最高的那一类其余被置为0。解决先单独读取new_class.tif统计像元值的最小值、最大值和nodata占比。如果nodata集中在影像外框可以在预测前用开运算或者掩膜把这些区域剔除只对有效区域循环。如果是缩放问题务必保证训练和预测用的是同一套标准化公式常见的两种是除以10000和除以255混用基本等于模型白训。5.2 模型总是把所有像元预测成同一个类别现象分类图上某一种地物占据了绝对主导比如整个影像都被预测成农田建筑物和道路完全没分出来。原因这是类别不平衡的典型症状切片阶段没有做类别配额某类地物样本量是其他类的几十倍CNN学会了“无脑输出多数类”来降低loss。另一个原因是标签图和影像没有严格对齐比如标签图整体偏移了十几个像元模型学到的窗口特征是错位的所有类别都会被扭曲到某个大类上。解决回到3.2的均衡采样逻辑把每类样本控制在接近数量级后再重训。同时检查标签图和影像的空间对齐在GIS软件里把两者叠加看道路和水体边界是否重合。标签错位是遥感深度学习里最隐蔽的问题它不会让训练报错只会在结果上表现为莫名其妙的偏置。5.3 加载h5模型时报shape不匹配或Unknown layer现象执行load_model(CNN_7class_3by3.h5)时报错有的提示Shape mismatch有的提示Unknown layer或Unknown loss function。原因Keras的h5文件把模型结构和权重一起打包但不同版本之间兼容性并不透明。如果训练脚本是在TensorFlow 2.4下保存的而你现在环境是TensorFlow 2.10以上部分旧格式可能触发警告或加载失败。另一个常见原因是模型里用到了自定义层或自定义lossload_model需要额外传custom_objects参数。解决我一般会固定环境压缩包里的README.md如果写了依赖版本直接用那个版本建虚拟环境最省事。没写的话先在命令行跑python -c import tensorflow as tf; print(tf.version)确认版本再尝试加载h5。如果实在加载不了退回一步不直接加载而是用2_trainModel.py重新训练一个模型反正数据和脚本都在。5.4 逐像元预测速度慢到无法接受现象3_predictNewData.py跑起来后每按一次回车都在刷进度预测一景几百MB的影像要好几个小时。原因问题出在循环里每次只预测一个样本。model.predict本身有批次开销单个样本也走一遍完整的预测流程相当于每次都在启动一次推理。在循环里调用predict而不是predict_on_batch开销会被放大到难以接受。解决改成批量预测一次性取几千个patch。具体做法是用np.lib.stride_tricks.sliding_window_view提取所有窗口或者直接用tf.image.extract_patches把整幅影像变成(N, 3, 3, bands)的数组再X X.transpose(0, 2, 3, 1)调整轴顺序后model.predict(X, batch_size1024)。我实测过批量预测能把几小时的推理压缩到几分钟这也是我后来每次跑预测脚本前必改的一个点。6. 进阶验证把example.tif和new_class.tif跑成一张分类图6.1 两条命令完成端到端验证拿到源码包之后我建议的验证路径是直接走预测分支最快看到结果# 先看依赖是否完整 python -c import rasterio, numpy, tensorflow # 直接预测验证预训练模型可用 python 3_predictNewData.py这一步只需要3_predictNewData.py和CNN_7class_3by3.h5、new_class.tif三个文件。跑通后打开输出的new_class_classified.tif和原影像做假彩色合成对比重点看水体边界、农田纹理、建设用地轮廓是否合理。如果输出明显异常回到第5章排查问题。想从零重训就把三步串起来# 1. 准备好标签图后切片 python 1_createImageChips.py # 2. 训练并覆盖h5 python 2_trainModel.py # 3. 对new_class.tif预测 python 3_predictNewData.py6.2 精度验证与多波段扩展方向没有独立测试集时精度验证只能靠目视解译这也是遥感分类的通用笨办法。我通常会做三件事第一把分类结果叠加到原始影像上透明度设0.4看边界是否锯齿严重第二随机选10个水体像元、10个植被像元回看它们原始的6波段光谱曲线确认类别归属合理第三对混淆最严重的阴影和水体做一个单独掩膜量一下两类面积比例是否离谱。如果7类分类已经满足项目目标可以尝试把输入从3×3扩大到5×5只改CHIP_SIZE和第一层input_shape其他不用动。更大的窗口能捕获农田纹理和道路走向但边缘细节会变模糊。也可以把NDVI和NDWI作为额外两个波段拼进去这样模型对植被和水体的辨识度会更好代价是切片脚本和预测脚本的波段数都要同步修改。我自己的习惯是每换一个数据集都会把“波段顺序、缩放系数、nodata定义”这三件套写进脚本头部写清楚才动手。这个习惯救了我很多次希望帮到你。本文还有配套的精品资源点击获取