企业官网建设流程全解析

简介：基于PyTorch的卷积神经网络遥感影像地物分类项目源码包，面向地理信息、遥感、人工智能等相关专业的高校学生、教师与科研人员，解决Landsat系列多光谱数据从样本制作、模型训练到地物分类的完整流程问题。包内三个Python脚本分别完成原始影像切片、七类地物模型训练与新增数据批量预测，配套的h5权重文件可直接加载，实现训练与推理的闭环；XML与TFW辅助文件记录影像元数据和地理坐标，便于分类结果在GIS中叠加上图。整个压缩包共10个文件，涵盖Python源码、GeoTIFF影像样例、XML配置、模型权重及说明文档等多种类型，大小仅14.88MB，结构紧凑、层次清晰，便于按模块对照学习。目前已有88人学习下载。代码经过严格测试、可稳定运行，不仅适合直接用于毕业设计、课程设计与项目初期演示，也可作为入门深度学习与遥感交叉应用的实践素材；理解算法流程后，可自行扩充地物类别并迁移到其他区域影像，拓展性较好。

1. 当Landsat影像遇上CNN：遥感地物分类为什么值得自己动手写一套Python源码

给县域做地表覆盖调查的同事跑过一个地物分类需求：面积不小，类别有耕地、林地、水体、建设用地，外加大片裸土和山地阴影。传统做法是先算NDVI、NDWI这类光谱指数再设阈值，一遇到农田边缘和阴影区，分类结果就碎成芝麻粒，补绘的时间比自动分类还长。换成“CNN深度学习遥感影像地物分类：Landsat数据处理Python源码”这套思路后，本质就一句话：让卷积神经网络从Landsat多光谱波段里自动学习地物的光谱-空间特征，替代手工阈值和人工勾绘。这条技术路线适合两类人：一类是要做区域尺度土地利用制图或生态参数反演的从业者，另一类是已经会Python和GDAL、想从随机森林一类传统分类器迁到深度学习的工程师。读懂这个方向，关键不是背网络结构，而是盯住一条完整链路：Landsat原始影像如何变成训练样本，CNN怎么训练和推理，预测结果如何输出成带地理坐标的专题图。

2. 读懂Landsat数据与CNN的适配逻辑：选型背后的三个关键理由

2.1 Landsat数据格式与波段物理意义：从DN值到反射率的换算

拿到Landsat数据的第一件事，不是急着把GeoTIFF塞进numpy，而是分清你手上是哪一级数据。Landsat 8/9 OLI传感器的每个波段以单独GeoTIFF文件存储，文件名类似LC08_L2SP_xxx_SR_B2.TIF，MTL元数据文件记录成像时间和缩放系数。Collection 2的L1级产品是未经大气校正的DN值，整型，范围通常在0到65535之间，数值大小同时受太阳高度角、大气散射和传感器增益影响，直接当反射率用会出大问题。L2级产品则是已经做过大气校正的地表反射率，但仍然以uint16存储，真实反射率被放大了10000倍左右，使用前必须按元数据里的scale factor缩放。

很多入门者把L2数据除以65535当成归一化，这是典型的翻车操作。L2表面反射率的有效记录范围不是0到65535，而是0到10000，除以65535会把反射率压到0.15以下，光谱间的区分度被严重压缩，CNN再强也学不出东西。正确的做法是乘以0.0001或元数据中给定的scale因子，把数值映射回0到1的物理反射率。

波段选择上，Landsat 8/9 OLI的反射波段有B2到B7六个30米波段：蓝光B2对水体敏感，绿光B3是植被反射峰所在，红光B4对应叶绿素吸收，近红外B5对植被整体反射强烈，短波红外B6和B7对土壤湿度、建筑物和矿物差异更敏感。做六波段输入是我最常用的配置，也就是B2到B7全部参与分类。六波段组合覆盖了可见光到短波红外的关键物理区间，CNN可以在卷积过程中自行组合出类似NDVI、NDWI甚至更复杂的特征表达式，不需要你手工预设。

2.2 为什么是CNN而不是随机森林或SVM：光谱-空间特征的关键差异

随机森林和SVM在遥感分类领域用了十几年，不是没有道理：样本量少的时候它们很稳，训练也快。但它们的常规输入是逐像元的光谱向量，也就是说模型看一个像元时只知道这个像元的六个波段数值，不知道它周围是什么。在30米分辨率的Landsat影像上，地物分类的难点往往不在光谱本身：耕地和草地的光谱曲线可能很接近，水体与山体阴影在可见光波段也容易混淆。真正能区分它们的线索，一部分在近红外和短波红外的数值组合里，另一部分就在空间上下文里——耕地有规则的田块纹理，草地相对均质，山体阴影通常连着山脊线。

CNN的卷积操作天然把邻域像素拉进了计算。一个卷积核在3x3窗口内滑动时，模型看到的是像元与其邻居的光谱联合模式；多个卷积层堆叠后，低层学到边缘和纹理，高层学到地块形状和组合关系。这种局部连接和参数共享的归纳偏置，恰好贴合遥感影像的地物空间结构，比把每个像元当孤立点看的SVM更合理，也比手工做GLCM纹理再喂给随机森林省去大量特征工程。

深度学习模型用在这个场景里也不需要堆得很深。Landsat的波段数只有六个，输入patch通常不超过15x15像素，信息量比ImageNet图像小几个数量级。我见过有人照搬ResNet50做Landsat分类，结果训练慢、过拟合明显，精度和轻量CNN几乎一样。这个方向真正有效的网络通常是几层卷积加池化的小模型，参数量从几万到几十万就够用了。遥感地物分类是一类“数据维度不高但空间关联重要”的任务，网络结构的重心应该放在如何组织空间上下文上，而不是追求层数。

2.3 从patch到像素：CNN输入数据组织方式

CNN不能像逐像元分类器那样直接输入一整景影像并输出每个像元的类别，常见做法是以待预测像元为中心裁剪一个正方形窗口，这个窗口叫做patch。模型输入patch，输出中心像元的类别概率。这样做的好处是模型既看到了中心像元的光谱值，也看到了它的空间上下文，代价是你得自己平衡patch的尺寸。

patch太小，比如3x3，上下文信息不够，模型退化成一个带邻域平滑的光谱分类器；patch太大，比如31x31，窗口内往往包含好几个地类，中心像元与远处像素的关系也不强，反而引入噪声。对Landsat这种30米分辨率影像，我一般把patch设在7到15之间。15x15的patch对应地面上450米x450米的范围，对耕地田块、居民区这类中等尺度地物已经足够；如果做山地或森林为主的区域，7到11更合适，避免混入过多相邻地类。

import numpy as np def extract_patches(array, patch_size=11, stride=5, valid_mask=None): """从多波段影像中按滑窗方式切patch，返回(n_patches, bands, h, w) array: 形状为 (bands, height, width) 的影像数组 patch_size: 窗口边长，通常取奇数，方便对齐中心像元 stride: 滑窗步长，小于patch_size时会得到重叠patch valid_mask: 有效像元掩膜，中心像元无效时跳过该patch """ bands, h, w = array.shape pad = patch_size // 2 # 边界填充，避免影像边缘像元无法裁剪 padded = np.pad(array, ((0, 0), (pad, pad), (pad, pad)), mode='reflect') if valid_mask is not None: padded_mask = np.pad(valid_mask, ((pad, pad), (pad, pad)), mode='constant', constant_values=0) patches = [] positions = [] for y in range(h): for x in range(w): cy, cx = y + pad, x + pad if valid_mask is not None and not padded_mask[cy, cx]: continue patch = padded[:, cy - pad: cy + pad + 1, cx - pad: cx + pad + 1] patches.append(patch) positions.append((y, x)) return np.stack(patches), positions

这段代码的逻辑是遍历每个像元位置，以它为中心裁取邻域。valid_mask用于跳过nodata或云掩膜区域，避免无效像元参与训练。stride参数在这个全遍历版本中体现为后续可以改成隔n个像元取一个样本，用于降低训练样本密度。如果做全图推断，通常stride设为1或取较小的重叠步长，推断结果再做投票平滑；如果做训练样本提取，stride可以放宽到5到10，减少相邻样本的冗余。patches和positions一一对应，positions记录每个patch左上角原影像坐标，最后拼接预测结果时要用到。

3. 环境配置与Landsat数据预处理：从GeoTIFF到numpy数组的落地流程

3.1 Python环境与核心库安装：最小依赖清单

做Landsat地物分类，Python环境不需要很复杂。核心库就三类：数值计算用numpy，地理数据读写用rasterio，深度学习框架用TensorFlow或PyTorch。rasterio是GDAL的Python封装，读写GeoTIFF比直接用GDAL的Python绑定顺手得多，而且能方便地保留坐标系、仿射变换等元数据。除了这三个，再加scipy做后处理时的滤波操作，scikit-learn做精度评估和样本划分，基本就能跑通全部流程。

conda create -n landcover python=3.10 -y conda activate landcover pip install numpy rasterio scipy scikit-learn pip install tensorflow

建议用conda新建独立环境，避免把系统Python装乱。python版本选3.10是一个稳妥选择，深度学习框架当前主流版本都支持。TensorFlow装CPU版也能跑这个任务，只是训练会慢一些；如果机器有N卡，先装CUDA和cuDNN再装GPU版。判断标准很简单：patch很小，模型很轻，训练集几十万样本规模，CPU训练一个epoch大概几分钟到十几分钟，不是不能等。

3.2 用rasterio读取Landsat影像并合成波段numpy数组

Landsat L2产品的每个波段是独立GeoTIFF，读取时要按固定顺序打开并堆叠成一个(bands, height, width)数组。波段顺序对CNN没有影响，但要保证训练和推断时顺序一致。我习惯按blue、green、red、nir、swir1、swir2排列，这个顺序对应Landsat 8的B2到B7。

import rasterio import numpy as np from pathlib import Path def read_landsat_sr(sr_dir): """读取Landsat L2地表反射率波段，堆叠为float32数组 sr_dir: 存放SR_B*.TIF文件的目录 返回: (bands, height, width)的顺序数组，以及来自第一个波段的地理元数据 """ band_names = { 'blue': 'SR_B2', 'green': 'SR_B3', 'red': 'SR_B4', 'nir': 'SR_B5', 'swir1': 'SR_B6', 'swir2': 'SR_B7', } bands = [] georef = None for key in ['blue', 'green', 'red', 'nir', 'swir1', 'swir2']: fpath = Path(sr_dir) / f"{band_names[key]}.TIF" with rasterio.open(fpath) as src: bands.append(src.read(1).astype('float32')) if georef is None: georef = (src.profile, src.transform) stack = np.stack(bands, axis=0) return stack, georef

读取时用float32而不是保留uint16，是为了后续归一化和CNN计算的精度。astype('float32')在波段数多、影像大时能省一半内存，6波段3000x3000的影像约占216MB，可接受。georef里存了投影坐标系和仿射变换参数，训练结束后把分类结果写回GeoTIFF时直接用这份元数据，保证输出与原始影像地理对齐。千万不要把每个波段分别读成python列表再循环拼接，那会有大量临时变量在内存里反复复制；一次性读进list再stack是最稳的做法。

拿到全波段数组后，按研究区范围裁剪能大幅减少后续计算量。裁剪很简单，直接用数组切片就行，但要知道切片后的地理坐标起点变了，输出结果时需要重新计算仿射变换参数。如果研究区边界是shapefile，更规范的做法是用rasterio.mask.mask按矢量边界裁剪；只是训练实验的话，直接切片够用。

# 假设研究区在原影像中的行列范围 x0, x1, y0, y1 = 1200, 4200, 800, 3800 img_crop = img[:, y0:y1, x0:x1] label_crop = label[y0:y1, x0:x1] # 用切片后的数组重新生成仿射变换 from rasterio.transform import from_origin new_transform = from_origin( georef[1][2] + x0 * georef[1][0], # 新的左上角x坐标 georef[1][5] + y0 * georef[1][4], # 新的左上角y坐标 georef[1][0], # 像素宽度 georef[1][4] # 像素高度 )

这里的关键是仿射变换的六个参数。原transform记录了左上角坐标、像素宽和高以及旋转项。裁剪后左上角坐标变了，必须通过原始行列号反算出新的左上角，否则后面写GeoTIFF时影像的坐标位置是错的。判断是否有地理错位有一个简单的自检方法：把原始影像和裁剪后影像分别转成png带坐标叠加，肉眼看到重叠一致才算通过。

3.3 训练样本制作：标签栅格与样本采样策略

CNN地物分类是有监督学习，需要每个训练patch对应一个真实类别标签。标签制作通常有两种来源：一种是已有土地利用矢量数据，栅格化成标签栅格；另一种是直接在影像上目视解译勾绘少量多边形，再转成栅格。对起步实验来说，后者的控制力更强，因为你能确保标签和影像在时间上基本一致。

from sklearn.model_selection import train_test_split import numpy as np def collect_training_samples(img_crop, label_crop, patch_size=11, samples_per_class=3000, seed=42): """从标签影像中随机抽取训练patch，按类别均衡采样 img_crop: (bands, h, w) 影像数组 label_crop: (h, w) 标签数组，类别从1开始编号，0为背景 """ rng = np.random.default_rng(seed) classes = np.unique(label_crop) classes = classes[classes > 0] X, y = [], [] pad = patch_size // 2 # 先检查边界，避免裁到影像外 for c in classes: ys, xs = np.where(label_crop == c) valid = [] for cy, cx in zip(ys, xs): if cy - pad >= 0 and cy + pad + 1 <= label_crop.shape[0] and cx - pad >= 0 and cx + pad + 1 <= label_crop.shape[1]: valid.append((cy, cx)) if len(valid) > samples_per_class: selected = rng.choice(len(valid), size=samples_per_class, replace=False) else: selected = np.arange(len(valid)) for idx in selected: cy, cx = valid[idx] X.append(img_crop[:, cy-pad: cy+pad+1, cx-pad: cx+pad+1]) y.append(c) return np.array(X), np.array(y)

这段代码按类别分别采样，每类最多取samples_per_class个patch，防止水体这类大而均质的类别在样本里占比过高。valid检查是容易漏掉的一步：Landsat边缘像元本身可能缺值，标签边界处的patch也可能跨不同类别，简单裁掉边界不像元最省事。如果把边界问题缓存给后续滑动窗口取patch，模型推理时会遇到同样问题，统一用reflect填充即可。

采样完成后，把数据集按8比2切分训练集和验证集。这里有一个常见坑：不能直接随机打乱切分，原因是相邻像元高度相关，随机切分会把同一块田地的像元同时分进训练和验证集，导致验证精度虚高。先按地理位置粗分为几大块，再把整块放进训练或验证，这才是可靠的做法。

X_train, X_val, y_train, y_val = train_test_split( X, y, test_size=0.2, stratify=y, random_state=42 ) # 标签从1到N转为0到N-1，便于softmax输出 y_train = y_train - 1 y_val = y_val - 1

stratify=y保证切分后各类别比例与原始数据一致，避免某类样本在验证集中缺失。类别数记得和模型输出层的n_classes参数对齐，后文中模型输出维度就是这里实际类别数量。

4. 搭建CNN地物分类模型：网络结构与训练参数详解

4.1 模型结构设计：轻量patch CNN还是U-Net

同一个Landsat分类任务，可以用patch CNN，也可以用U-Net这类全卷积语义分割网络。两者差别在于输出粒度：patch CNN输入一个patch，输出中心像元的类别；U-Net输入一整块影像，直接输出整块每个像元的分类结果。U-Net在边界精细度上有优势，但训练需要像素级标签，显存消耗大，实现也更复杂。

对绝大多数Landsat区域分类需求，我建议先从patch CNN开始。原因很直接：Landsat是30米中分辨率影像，一张patch里的空间信息并不需要像高分辨率遥感那样逐像素细化；patch CNN实现简单、训练快、调参直观，出图效果足够好。U-Net更像这些工作后期提升精度时考虑的方向。

from tensorflow.keras import layers, models import tensorflow as tf def make_patch_cnn(patch_size=11, n_bands=6, n_classes=5): """构建轻量patch分类CNN patch_size: 输入patch边长 n_bands: 输入波段数 n_classes: 地物类别数 """ inputs = tf.keras.Input(shape=(patch_size, patch_size, n_bands)) x = layers.Conv2D(32, 3, activation='relu', padding='same')(inputs) x = layers.MaxPooling2D(2)(x) x = layers.Conv2D(64, 3, activation='relu', padding='same')(x) x = layers.MaxPooling2D(2)(x) x = layers.Conv2D(128, 3, activation='relu', padding='same')(x) x = layers.GlobalAveragePooling2D()(x) x = layers.Dropout(0.5)(x) outputs = layers.Dense(n_classes, activation='softmax')(x) return models.Model(inputs, outputs) model = make_patch_cnn(patch_size=11, n_bands=6, n_classes=5) model.summary()

网络包含三个卷积层，卷积核都是3x3，通道数从32增加到128，中间穿插最大池化。GlobalAveragePooling2D把最后一个卷积层的特征图压缩成128维向量，替代Flatten加全连接的做法，能明显减少参数量并缓解过拟合。Dropout(0.5)在全连接前随机丢弃一半神经元，是这个小模型里最有效的正则化手段。

对于11x11的输入patch，经过两次池化后特征图缩到约3x3，第三个卷积层实际上主要承担特征抽象，信息保留已经足够。如果你用15x15的patch，网络可以不变，因为三层卷积的感受野大约是7x7到9x9，对15x15输入仍能覆盖中心区域的主要上下文。如果做多时相Landsat数据，比如输入包含两个时相的12个波段，只需把n_bands改成12，结构本身不需要调整。

4.2 训练脚本与参数设置：让loss真正下降

模型训练最核心的参数是学习率、batch size和epoch数。这个任务我用Adam优化器，初始学习率0.001，一般不需要额外调。batch size默认128，如果显存爆了降到32或64；经验法则是在模型收敛的前提下尽量用大batch，梯度更稳定。epoch设30，配合EarlyStopping看验证loss，通常十几个epoch就能收敛。

# 标签做one-hot编码 y_train_onehot = tf.keras.utils.to_categorical(y_train, num_classes=5) y_val_onehot = tf.keras.utils.to_categorical(y_val, num_classes=5) model.compile( optimizer='adam', loss='categorical_crossentropy', metrics=['accuracy'] ) history = model.fit( X_train, y_train_onehot, batch_size=128, epochs=30, validation_data=(X_val, y_val_onehot), callbacks=[ tf.keras.callbacks.EarlyStopping(patience=5, restore_best_weights=True) ] )

X_train的形状是(num_samples, 11, 11, 6)，TensorFlow要求的通道最后布局。如果你的数据从numpy stack直接来，形状是(num_samples, 6, 11, 11)，需要做np.transpose(0, 2, 3, 1)把波段轴移到最后。这是CNN实现里最常见的形状错误，报错通常是维度不匹配。

EarlyStopping的patience设为5，意味着验证loss连续5个epoch不下降就停止训练，并恢复到验证loss最低的权重。遥感样本往往有一定噪声，验证曲线不会单调下降，patience太小容易提前停，patience太大又浪费时间，5是一个比较稳的起点。训练结束后用model.evaluate(X_val, y_val_onehot)看整体精度，但这只是第一道检查，真正的精度判断在第六章用混淆矩阵做。

4.3 类别不平衡与损失函数：别硬训

地物类别天然不平衡。一个县域影像里水体可能只占5%，耕地占40%，如果按原始比例采样本，模型会把所有不确定像元都猜成耕地，因为这样整体准确率最高。解决采样层面的办法是像前面那样每类固定采样上限，另一个角度是在损失函数上给稀有类别更高的惩罚权重。

from sklearn.utils.class_weight import compute_class_weight import numpy as np classes = np.array([0, 1, 2, 3, 4]) # 与label编码后的类别对应 weights = compute_class_weight('balanced', classes=classes, y=y_train) class_weight = dict(zip(classes, weights)) print(class_weight) model.fit( X_train, y_train_onehot, batch_size=128, epochs=30, validation_data=(X_val, y_val_onehot), class_weight=class_weight, # 控制每类样本在loss中的权重 callbacks=[tf.keras.callbacks.EarlyStopping(patience=5, restore_best_weights=True)] )

compute_class_weight的'balanced'模式会自动按样本数量的反比计算权重，样本少的类别loss权重更大。这种方式比直接在采样时均衡更平滑，因为训练时仍然保留了各类别的光谱分布多样性。如果某类样本实在太少，少于100个patch，权重设得再高也难以学好，这时优先补充样本量，而不是依赖权重调整。

还有一类做法是用focal loss或dice loss处理不平衡。但从实践看，Landsat地物分类的类别不平衡严重程度远低于医学图像分割，修正采样本比例加class_weight已经能解决绝大多数问题，没必要引入复杂损失函数增加调参难度。

5. 遥感分类实战避坑：5个最常见的翻车现场与排查方案

5.1 叠加后地物错位严重：投影坐标系不一致

现象：分类结果和验证影像叠加后错位几十米甚至上百米，边界明显对不上。

原因：Landsat L2产品通常自带UTM投影坐标，但不同景可能落在不同UTM分带；标签数据如果来自其他来源，比如地方国土部门的shapefile用的是CGCS2000或WGS84经纬度，你和影像直接叠加必然错位。

解决：用rasterio打开两个文件，先看一眼坐标系再动手。

import rasterio with rasterio.open('landsat_sr.tif') as src: print(src.crs) # 输出影像坐标系 with rasterio.open('label.tif') as src_label: print(src_label.crs) # 输出标签坐标系

如果两者crs不一致，把标签重投影到影像坐标系。rasterio支持在读取时直接做重投影，但更省事的路径是用rasterio.warp.reproject把标签重采样到影像的transform和分辨率。做完之后随机选几个地物特征点，比如道路交叉口，在原始影像和重采样标签上看坐标是否对齐。这一步能省掉后面所有的错位烦恼。

5.2 loss曲线第一轮就飞到几十：像素没有归一化

现象：第一个epoch的loss值是三位数，accuracy却接近某一类别占比，比如始终0.6不动。

原因：影像数组仍停留在0到10000或0到65535的尺度，梯度计算时数值范围过大，尤其是还有nodata值为0混在里面时，模型对数值尺度极其敏感。

解决：训练前检查数组统计量，然后按波段分别归一化。

print(X_train.min(), X_train.max(), X_train.dtype)

如果最大值远大于1，先做反射率缩放：乘以0.0001把L2数据映射到0到1附近。如果不同波段数值范围差异大，再做逐波段标准化，每个波段减去均值除以标准差。我习惯的默认做法是先乘scale到0到1，再逐波段标准化。要特别注意影像中的0值，如果0是nodata伪值，归一化前要替换成有效的填充值或做掩膜，否则0会影响均值方差计算。

5.3 一训练内存就爆：整景影像一次性展开

现象：16G内存的机器，训练启动没多久就报MemoryError，或者进入swap后卡死。

原因：很多人在采样阶段把整景影像按滑窗切成了几十万甚至上百万个patch，每个patch在内存中独立存储，加上numpy的复制和过程中产生的临时变量，内存轻松冲破20G。

解决：训练时用生成器或分批从磁盘读取patch，不要把所有样本一次性堆进内存。如果样本集是几十万patch、每个patch只有6x11x11，单块数据约10G，勉强能放下，但训练过程中还会有批次数据转换和梯度计算的额外开销。

def batch_generator(X, y, batch_size=128): """轻量batch生成器，按epoch随机打乱后逐个返回batch""" n = len(X) while True: idx = np.random.permutation(n) for start in range(0, n, batch_size): batch_idx = idx[start:start + batch_size] yield X[batch_idx], y[batch_idx]

把batch_generator传给model.fit，配合steps_per_epoch设置为样本总数除以batch_size。这样每个epoch只在循环内处理一个batch的临时数据，内存占用稳定。如果连原始patch数组都放不内存，说明样本量太大，更彻底的方案是把影像按块存成npy文件，训练时按块读入。

5.4 验证精度高达96%，实际分类图却不忍直视

现象：混淆矩阵很漂亮，Kappa超过0.9，但生成的全图碎斑严重，行政边界处大面积错误，农田和草地互相混。

原因：样本切分时按像元随机打散，训练和验证样本来自同一地块的相邻像元，空间上高度相关。模型记住的是这块田的光谱纹理特征，而不是泛化出“农田”这个抽象类别，验证时恰好撞上记忆过的区域，精度自然虚高。

解决：按空间块划分训练验证集。把研究区分成1km x 1km的小块，随机拿八成块训练、两成块验证，或者更简单粗暴一点，用影像上一半区域训练、下另一半区域验证。虽然数值上验证精度会下降，比如从96%降到82%，但这个82%才是你真正部署时的预期精度。看到验证精度比全图肉眼效果高太多时，先怀疑采样方式，不要怀疑模型出了问题。

5.5 云和阴影没掩膜：把云当成了不透水面

现象：分类结果图上所有云区域都成了建设用地或裸地，云影区域全是水体。

原因：Landsat L2产品虽然做过了大气校正，但云和云影仍然存在，这些像元的反射率特征完全不体现地表信息。没有质量掩膜的模型只能按光谱强行归类，云高亮像元非常接近建筑物屋顶和裸土的光谱范围。

解决：用L2产品自带的SCL质量波段生成掩膜，只对真实地表像元训练和预测。SCL波段中地表覆盖有固定的值范围，云和云影另有编码，具体对照参考产品说明文档。读取SCL后生成bool掩膜，在训练样本采集中把掩膜为无效的像元直接跳过，推断完成后把无效像元重新写成nodata。

6. 精度验证与结果后处理：让CNN输出变成一张能用的专题图

6.1 混淆矩阵与Kappa系数：分类精度怎么算才可信

模型训练完，先别急着全图推理。拿一个独立验证集，注意是空间独立的验证集，计算混淆矩阵和Kappa系数。训练时打印的accuracy看着高是因为类别不平衡占便宜，混淆矩阵能逐类看出哪些地物互相混淆。

from sklearn.metrics import confusion_matrix, classification_report, cohen_kappa_score y_pred = model.predict(X_val) y_pred_class = np.argmax(y_pred, axis=1) # 真实标签如果是one-hot，先转回类别编号 # y_val_class = np.argmax(y_val_onehot, axis=1) cm = confusion_matrix(y_val_class, y_pred_class) print(cm) # 逐类精度报告 print(classification_report(y_val_class, y_pred_class, target_names=['耕地', '林地', '水体', '建设用地', '裸地'])) kappa = cohen_kappa_score(y_val_class, y_pred_class) print(f'Kappa: {kappa:.3f}')

混淆矩阵的行是真实类别，列是预测类别。对角线上是正确分类数量，非对角线上是误分情况。比如行列位置(0,1)数值大，说明耕地被大量误分成林地。classification_report里每个类别有precision、recall和f1-score，precision表示预测成该类别的像元里有多少是真该类，recall表示真该类的像元里有多少被找回来。对地物分类来说，recall更影响制图效果，因为它反映某一类地物是否被完整保存下来。

Kappa大于0.8算高度一致，0.6到0.8算中等。如果Kappa偏低，回到混淆矩阵看是哪两类在混。耕地和草地光谱接近是常态，解决方向是增加时相信息或者加入地形辅助特征，而不是盲目加深网络。

6.2 分类结果后处理与专题图输出：从逐patch到完整分类图

全图推理时不能整景直接塞进模型，还是要按patch来。推理的步长可以设置得比训练时更密，比如patch_size是11，stride取5，让相邻patch有重叠，预测结果做多数投票或只取中心像元值。重叠投票能明显减少斑块噪声，代价是计算量成倍增加。对区域制图来说，取中心像元的预测值加上一个小的中值滤波已经够用。

from scipy.ndimage import median_filter # 假设pred_full是推理得到的(height, width)分类编号数组 # 中值滤波去除孤立小斑块，3x3是比较温和的选择 smoothed = median_filter(pred_full, size=3)

中值滤波会把小于窗口尺寸的零散像元替换成周围多数类别，直接用分类编号做滤波是安全的，不需要one-hot。size取3只影响单像元噪点，5会开始抹掉细小线状地物。对建设用地这类边界复杂的类别，宁可保留一些碎点也别用太大的核，否则边界畸形。输出GeoTIFF才是关键，分类图没了地理坐标就是一张废图。

import rasterio from rasterio.transform import from_origin profile = { 'driver': 'GTiff', 'height': smoothed.shape[0], 'width': smoothed.shape[1], 'count': 1, 'dtype': 'uint8', 'crs': img_crs, 'transform': new_transform, 'nodata': 0, } with rasterio.open('landcover_result.tif', 'w', **profile) as dst: dst.write(smoothed.astype('uint8'), 1)

这里的crs和transform必须与采样时保持一致，丢失其中一个坐标信息都会错。写完先用rasterio打开，在GIS软件里叠加原始影像做一次目视抽查。我自己有个习惯：推理结果先从视觉上挑几个已知区域看类别对不对，再谈指标；数值再高、大图肉眼看着不对，一定是前面某些环节出了问题，这时候回头查预处理和样本采样，别急着调网络。还有一件经常做且值得做的事，是把验证集的预测概率保存下来，后续想调整类别阈值或者合并相似地类时不用重新训练模型，你等于给自己留了一颗后悔药。这套CNN地物分类流程里，真正决定上限的不是卷积层数有多深，而是你对Landsat数据、样本分布和地理坐标这三件事有没有控制住，希望帮到你。

本文还有配套的精品资源，点击获取