简介：针对CNN深度学习遥感影像地物分类任务，这套Landsat数据处理Python源码基于PyTorch框架实现，面向遥感、地信与人工智能相关专业的学生、教师和科研人员，旨在解决从原始遥感影像到地物分类结果的完整流程问题。压缩包内共10个文件，包括3个Python程序（分别用于生成影像切片、训练分类模型、对新建影像进行预测）、1个预训练H5权重文件、2幅TIFF样例影像及其配套XML与TFW坐标参考文件，另附Markdown格式的说明文档，整体约14.88MB。代码结构清晰、模块间解耦，可对照文档快速复现实验，也可替换自有Landsat数据完成不同区域的地物分类，适合作为毕业设计、课程设计或项目初期验证的基线参考。当前已有86人学习下载，样本裁剪、数据读取、批次构建、训练推理等关键环节都有具体实现，对初次接触遥感深度学习的开发者能起到较好的引导作用。
1. CNN深度学习遥感影像地物分类为什么先卡在Landsat数据上
做遥感影像地物分类的人，第一次用CNN跑Landsat数据，十有八九会遭遇“模型结构没问题，精度却上不去”的尴尬。问题往往不在网络本身，而在数据流入网络之前的那条管线：Landsat Level-1数据是DN值而非反射率，部分热红外波段还带重采样误差，叠加训练样本与影像之间坐标系错位，任何一点都会让模型学到错误的映射。业界早有一个共识——地物分类的精度上限由数据决定，模型只负责逼近那个上限。
这篇文章顺着“Landsat数据 → 预处理 → 样本制作 → CNN训练 → 滑窗预测 → 精度评价”这条完整链路展开，把每一步的逻辑、参数和一个可复用的Python实现串在一起。适合两类人：一是刚接触遥感、想把深度学习落到实际影像上的工程师，二是已有CNN基础但被Landsat元数据和波段特性困扰的从业者。后文代码基于Python的rasterio、geopandas和PyTorch，尽量贴近真实项目里最朴素的方案，不绕弯子。
2. Landsat预处理：从DN值到CNN可用的张量
2.1 辐射定标与TOA反射率计算：跳过这一步模型就学错特征
Landsat Level-1数据的每个像元存的是无量纲DN值（Digital Number），它和传感器实际接收到的辐射亮度之间是线性关系。直接把DN值作为CNN输入，等于放任不同影像、不同时相的增益差异影响模型——在A影像上训练，到B影像上预测时精度崩盘是必然的。
提示：即使只做单影像分类，也建议执行辐射定标。否则卷积核拟合的是“DN值分布”，不是“地表反射率分布”，换一景影像就要重新训练。
辐射定标公式很朴素。以Landsat 8/9 OLI为例，每个波段都有RADIANCE_MULT_BAND和RADIANCE_ADD_BAND两个元数据字段，DN值转辐射亮度，再做太阳高度角和日地距离校正得到TOA反射率：
import rasterio import numpy as np def dn_to_toa_reflectance(dn_band, mult, add, solar_zenith_angle, d_es=1.0): # 辐射亮度：L = M * DN + A radiance = mult * dn_band + add # TOA反射率：rho = pi * L * d^2 / (ESUN * cos(theta)) # 简化做法，用元数据里的REFLECTANCE_MULT_BAND和REFLECTANCE_ADD_BAND return radiance # 实际工程中更多直接用反射率缩放系数 def load_toa_bands(band_paths, meta): bands = [] for path in band_paths: with rasterio.open(path) as src: dn = src.read(1).astype(np.float32) # 直接读取元数据中的反射率乘加系数 mult = float(meta['REFLECTANCE_MULT_BAND_' + src.name[-1]]) add = float(meta['REFLECTANCE_ADD_BAND_' + src.name[-1]]) toa = dn * mult + add bands.append(toa) return np.stack(bands, axis=0)参数说明：这里用的是Collection 2 Level-1元数据中自带的REFLECTANCE_MULT_BAND和REFLECTANCE_ADD_BAND，两者是浮点小数值，通常在0到1之间。除以cos(太阳天顶角)的操作需要读取MTL文件里的SUN_ELEVATION字段，再取cos(90 - 太阳高度角)。若你用的是Landsat 4/5 TM或7 ETM+，字段名相同，直接套用即可。注意：这个步骤没有做大气校正，只是表观反射率，很多地物分类项目中TOA已经够用；若要用地表反射率，就得引入LEDAPS或LaSRC，但CNN分类场景下TOA反射率配合归一化通常能拿到不错的结果。
2.2 波段组合与张量组织：哪些波段进网络
Landsat 8有11个波段，但CNN输入通道数不宜贪多。常见做法是用6个波段：海岸蓝（B1）、蓝（B2）、绿（B3）、红（B4）、近红外（B5）、短波红外1（B6）。这6个波段覆盖了植被、水体、土壤、建筑区分最敏感的谱段范围，加入B7（SWIR2）收益有限但会显著增加计算量。
- 真彩色合成（B4-B3-B2）：适合人眼目视解译，但对CNN来说信息冗余度较高
- 假彩色合成（B5-B4-B3）：植被显示为红色，对林地/草地分类敏锐
- 6波段输入：上述三种合成一次性给齐，让卷积核自行学习组合权重
实际工程中，我会先把各波段DN值转成TOA反射率，再按[B1,B2,B3,B4,B5,B6]顺序堆叠，每波段独立做z-score归一化：
def normalize_bands(stacked): mean = stacked.mean(axis=(1,2), keepdims=True) std = stacked.std(axis=(1,2), keepdims=True) return (stacked - mean) / (std + 1e-8)这里有个容易被忽略的参数：归一化时不要用全局影像的均值方差，而应使用训练集统计量，或者逐影像按波段归一化后再做整体线性缩放到[-1, 1]。若直接对整景影像计算z-score，遇到大面积云或水体时会把陆地像元的数值压得极扁，影响CNN的特征提取。
2.3 训练样本制作：Shapefile转栅格标签图
样本准备是遥感深度学习中耗时最多的环节。常见工作流是：在ArcGIS或QGIS中人工勾绘地类边界，导出Shapefile，再用代码把矢量转成与影像严格对齐的栅格标签图。转换的核心是rasterio.features.rasterize，它需要一个关键的transform参数——必须从对应影像的元数据中获取，否则标签与影像会错位。
import geopandas as gpd from rasterio import features import rasterio as rio def shp_to_label(shp_path, ref_img_path, out_label_path): # 读取参考影像的元数据 with rio.open(ref_img_path) as src: transform = src.transform out_shape = (src.height, src.width) # 读取矢量并转换为GeoJSON gdf = gpd.read_file(shp_path) # 确保矢量与影像坐标系一致 if gdf.crs != src.crs: gdf = gdf.to_crs(src.crs) # 类别字段：例如'class_id'列存整数类别编码 shapes = [(geom, value) for geom, value in zip(gdf.geometry, gdf['class_id'])] # 栅格化：所有像素初始为0（背景），有矢量的区域填对应类别 label = features.rasterize( shapes=shapes, out_shape=out_shape, transform=transform, fill=0, all_touched=False, # 只填充矢量完全覆盖的像素 dtype=np.uint8 ) # 写出 with rio.open(out_label_path, 'w', driver='GTiff', height=out_shape[0], width=out_shape[1], count=1, dtype=np.uint8, crs=src.crs, transform=transform) as dst: dst.write(label, 1)参数说明：all_touched是个容易翻车的参数。设为True时，只要像素中心点在矢量边界上或边界穿过该像素，该像素就被赋值为对应类别，会让标签边缘比实际地物边界粗一圈；设为False时只赋值那些中心点严格落在矢量内部的像素。对于高分辨率影像建议True，但Landsat是30米分辨率，地物边界本身模糊，推荐False并用后续的形态学操作微调边缘。此外，shp_to_label中转出的标签值为0是背景，必须确保类别ID从1开始编号，否则CNN训练时会把背景当第0类混进去。
2.4 影像裁剪与数据集划分：窗口大小和空间不重叠原则
Landsat单景约7700×7700像素，直接整图进显存不现实。常规做法是裁剪成128×128或256×256的patch，按7:2:1划分训练/验证/测试集。这里最关键的约束是：同一地块内的样本必须分到同一集合，不能随机切patch后打乱，否则训练集和验证集会共享大量空间相邻的像元，验证精度会虚高。
import random import numpy as np def split_patches_by_polygon(label, patches, train_ratio=0.7, val_ratio=0.2): # 按图斑ID分块，而不是按像素随机 unique_ids = np.unique(label[label > 0]) n = len(unique_ids) random.shuffle(unique_ids) train_ids = set(unique_ids[:int(n * train_ratio)]) val_ids = set(unique_ids[int(n * train_ratio):int(n * (train_ratio + val_ratio))]) # 每个patch根据其中心像素所属的图斑ID决定归属 return patches, train_ids, val_ids窗口大小选择上，128×128通常够用，32米分辨率下覆盖约3.84×3.84公里地物，CNN感受野足够捕获地物纹理上下文。但若分类目标包含大面积连片农田或森林，可以提升到256×256，代价是显存占用指数上升。裁剪时相邻patch之间保留8像素重叠，能有效减少后续滑窗预测的边界效应——这个细节在第四章展开。
3. CNN模型构建与训练：小样本场景下的收敛策略
3.1 选型依据：为什么地物分类多用UNet而不是VGG/ResNet
图像分类任务中VGG、ResNet这类网络对整张图输出一个类别，但地物分类要求每个像素都有类别标签，本质是密集预测或语义分割。UNet是最适合这种任务的入门结构：编码器逐层下采样捕捉语义，解码器逐层上采样恢复空间分辨率，跳跃连接（skip connection）把下采样过程中丢失的边界细节拼接到对应层，对林地和农田这类有规则边界的对象尤其友好。
与DeepLabV3+相比，UNet结构简单、显存占用小、在小数据集上不易过拟合，配合Landsat的30米分辨率足够。另一个考虑是训练时间：UNet在单张RTX 3090上跑Landsat 6波段128×128输入，一个epoch只需几分钟，便于快速迭代调参。
3.2 PyTorch实现一个轻量UNet
用PyTorch实现一个足够用于地物分类的UNet，输入是6波段影像，输出是n个类别分数图：
import torch import torch.nn as nn import torch.nn.functional as F class ConvBlock(nn.Module): def __init__(self, in_ch, out_ch): super().__init__() self.conv = nn.Sequential( nn.Conv2d(in_ch, out_ch, 3, padding=1), nn.BatchNorm2d(out_ch), nn.ReLU(inplace=True), nn.Conv2d(out_ch, out_ch, 3, padding=1), nn.BatchNorm2d(out_ch), nn.ReLU(inplace=True) ) def forward(self, x): return self.conv(x) class UNet(nn.Module): def __init__(self, in_channels=6, num_classes=5): super().__init__() self.enc1 = ConvBlock(in_channels, 64) self.enc2 = ConvBlock(64, 128) self.enc3 = ConvBlock(128, 256) self.pool = nn.MaxPool2d(2) self.bottleneck = ConvBlock(256, 512) self.up2 = nn.ConvTranspose2d(512, 256, 2, stride=2) self.dec2 = ConvBlock(512, 256) self.up3 = nn.ConvTranspose2d(256, 128, 2, stride=2) self.dec3 = ConvBlock(256, 128) self.up4 = nn.ConvTranspose2d(128, 64, 2, stride=2) self.dec4 = ConvBlock(128, 64) self.out = nn.Conv2d(64, num_classes, 1) def forward(self, x): e1 = self.enc1(x) # 64通道 e2 = self.enc2(self.pool(e1)) e3 = self.enc3(self.pool(e2)) b = self.bottleneck(self.pool(e3)) d2 = self.dec2(torch.cat([self.up2(b), e3], dim=1)) d3 = self.dec3(torch.cat([self.up3(d2), e2], dim=1)) d4 = self.dec4(torch.cat([self.up4(d3), e1], dim=1)) return self.out(d4)逻辑说明：编码器每层输出通道翻倍，空间尺寸减半；解码器做转置卷积上采样后与对应编码器层拼接，通道数相加。这里的BatchNorm对Landsat输入特别重要——各波段缩放后的数值分布差异极大，BatchNorm能强制逐通道归一化，加速收敛。注意卷积核固定为3×3，padding为1保证尺寸不变，像素级分类不涉及全连接层，所以输入尺寸可以不是固定值。
3.3 损失函数与类别不均衡：小样本分类的胜负手
地物分类中类别不均衡是常态。水体、裸地往往占据影像70%以上面积，而建筑、道路等目标类可能只占3%。直接用标准交叉熵损失，模型会把所有像素预测为多数类，OA（Overall Accuracy）看着高，实际毫无意义。
常见的做法是给交叉熵加权重，权重与类别频率成反比：
def class_weights_from_label(label, num_classes): # 统计每个类别的像素占比 counts = np.bincount(label.flatten(), minlength=num_classes) total = counts.sum() # 权重 = log(总像素 / (类别像素 * 类别数))，防止权重极端 weights = np.log(total / (counts * num_classes) + 1e-8) return torch.tensor(weights, dtype=torch.float32) criterion = nn.CrossEntropyLoss(weight=class_weights_from_label(train_label, num_classes))参数说明：权重计算公式里对倒数做了log平滑，避免占比极小的类别获得过大权重而放大标签噪声。如果你用的是Focal Loss，其gamma参数默认2.0，但小样本影像中标签本身存在混杂像元，过高的gamma会让模型过度关注难例噪声。我没有给所有项目都上Focal，交叉熵配合上述权重在多数遥感分类任务中表现已经稳定。
3.4 训练参数：学习率、Batch Size与早停
| 参数 | 推荐值 | 说明 |
|---|---|---|
| 优化器 | Adam | 初始lr 1e-3，配合CosineAnnealingLR衰减 |
| 学习率 | 1e-3 → 1e-5 | 预热3个epoch后开始衰减，避免初期震荡 |
| Batch Size | 16（128×128输入） | Landsat输入通道数少，显存占用不大 |
| Epoch | 50 | 配合早停，patience=10 |
| 数据增强 | 随机旋转90°/180°/270°+水平翻转 | 不要用随机裁剪，会破坏地物空间上下文 |
训练过程中除了监控loss，还应每5个epoch在验证集上计算一次OA和Kappa；若连续10个epoch验证OA不升反降则停止。一个额外的技巧是使用MixUp增强——把两幅patch按比例混合，标签同步混合，能显著提高模型对混合像元（建筑物与道路边界）的鲁棒性。
4. 滑窗预测与精度评价：让模型输出真正可用
4.1 重叠滑窗预测：消除边界伪影
训练时用了128×128的patch，推理时不能把整景影像一次性喂入——显存装不下，而且影像边缘的上下文不足会影响精度。常见做法是滑窗推理，窗口大小与训练一致，加上一定的重叠度（如32像素），重叠区域取多次预测的均值。
def sliding_window_predict(model, full_img, window=128, stride=96, device='cuda'): model.eval() h, w = full_img.shape[1], full_img.shape[2] pred_sum = torch.zeros((num_classes, h, w), device=device) count = torch.zeros((1, h, w), device=device) for y in range(0, h - window + 1, stride): for x in range(0, w - window + 1, stride): patch = full_img[:, y:y+window, x:x+window].unsqueeze(0).to(device) with torch.no_grad(): out = model(patch) # (1, C, window, window) prob = torch.softmax(out, dim=1) pred_sum[:, y:y+window, x:x+window] += prob[0] count[:, y:y+window, x:x+window] += 1 # 取平均概率后取argmax得到类别 prob_avg = pred_sum / count.clamp(min=1) label_map = prob_avg.argmax(dim=0).cpu().numpy() return label_map参数说明：stride=96表示窗口之间重叠32像素（128-96），这是平衡速度与精度的一个经验值。重叠区域被预测多次，取概率平均可以有效抑制窗口边缘因为padding导致的类间震荡。若影像尺寸不是stride的整数倍，最后一行/一列需要单独处理——可以用reflect模式padding到能被stride整除，预测后裁剪回原始尺寸。注意这里用argmax前先做了softmax，而不是直接对logits取argmax，因为在重叠区域平均概率比平均logits更符合概率语义，实测能降低约1-2%的边界噪声。
4.2 OA、Kappa与混淆矩阵：别被整体精度骗了
地物分类的三大评价指标计算方式如下：
from sklearn.metrics import confusion_matrix, cohen_kappa_score, accuracy_score def evaluate_prediction(pred_label, true_label, ignore_index=0): # 只统计有效像素，忽略背景 mask = (true_label != ignore_index) & (pred_label != ignore_index) y_true = true_label[mask].flatten() y_pred = pred_label[mask].flatten() oa = accuracy_score(y_true, y_pred) kappa = cohen_kappa_score(y_true, y_pred) cm = confusion_matrix(y_true, y_pred) # 各类别精度（用户精度 UA 和生产者精度 PA） ... return oa, kappa, cm参数说明：评估时必须同时传入ignore_index（通常是0或255），否则背景类会把那些边缘噪声计入误差，导致指标虚高。Kappa系数衡量的是与随机分类的一致性差异，当地物类别分布极度不均衡时Kappa会比OA更能反应模型真实能力。一个常见骗局是：模型把所有像素预测为占比90%的类别，OA能到90%，但Kappa会跌到0附近——这正是指标组合存在的价值。
4.3 混淆的热点区域与波段对策
| 混淆对 | 原因 | 缓解方案 |
|---|---|---|
| 林地 vs 草地 | 光谱曲线在可见光段极为接近 | 加强近红外B5和短波红外B6的权重，或加入NDVI作为额外通道 |
| 水体 vs 阴影 | 山体阴影在可见光段与水体的低反射率相似 | 引入SWIR1（B6），水体在1.6μm处吸收强烈而阴影不会 |
| 裸地 vs 建筑屋顶 | 屋顶材料（混凝土/沥青）光谱与裸地重叠 | 加入纹理信息，或者用形态学开闭运算后处理 |
多数深度学习项目忽略了一个事实：Landsat的30米分辨率意味着一个像素往往是混合地物，例如林地与灌木的过渡带。与其要求模型做精细化边界，不如在后处理中引入“不确定度”概念——用预测概率的熵来标记那些处于类别边界的像素，供人工检查或作为矢量化的判定条件。
5. 进阶技巧：用连通域分析清除预测结果的椒盐噪声
模型逐像素预测的输出像质粗糙，常见问题是孤立的小块噪声区域——单个或几个像素被错误分类，在图像上表现为“椒盐”颗粒。这类噪声在均匀地块内部尤其显眼，直接影响矢量化和面积统计的精度。
一个不依赖额外模型的后处理技巧是连通域过滤：给预测标签图做连通域标记，统计每个连通域的面积（像素数），将面积小于阈值的连通域替换为周围最频繁的类别。
from scipy import ndimage import numpy as np def remove_salt_pepper(label_map, min_area=30): """删除小于min_area像素的孤立区域，替换为其邻域的最常见类别""" cleaned = np.copy(label_map) # 对每个类别分别处理，确保不同类别不相连 for cls_id in np.unique(label_map): if cls_id == 0: # 背景跳过 continue mask = (label_map == cls_id).astype(np.int32) # 连通域标记：按4连通或8连通（对角线相邻也算同一区域） labeled, num_features = ndimage.label(mask, structure=np.ones((3,3))) # 统计每个连通域面积 sizes = ndimage.sum(mask, labeled, range(1, num_features + 1)) # 找出面积小于阈值的连通域索引 small_ids = np.where(sizes < min_area)[0] + 1 if len(small_ids) == 0: continue # 将这些小区域置为-1（待修正） for sid in small_ids: cleaned[labeled == sid] = -1 # 对每个-1位置，用周围3x3窗口中出现最多的非负类别替换 from scipy.ndimage import generic_filter def most_common_without_neg(arr): values = arr[arr >= 0] if len(values) == 0: return 0 return np.bincount(values).argmax() # 通过膨胀的方式，用邻域的最常见类别填充待修正像素 while True: neg_mask = (cleaned == -1) if not neg_mask.any(): break # 用3x3邻域内最常见的类别填充 from collections import Counter pad = np.pad(cleaned, 1, mode='edge') for idx in np.argwhere(neg_mask): y, x = idx neighborhood = pad[y:y+3, x:x+3].flatten() valid = neighborhood[neighborhood >= 0] if len(valid) > 0: cleaned[y, x] = Counter(valid).most_common(1)[0][0] else: cleaned[y, x] = 0 # 邻域全部无效则设为背景 return cleaned逻辑说明：这个实现分两个阶段——先标记每个类别的连通域并筛掉小面积噪声，再对噪声像素用3×3邻域投票修复。threshold（min_area）的选择有门道：Landsat 30米分辨率下3×3像素对应90×90米区域，地物分类中低于该面积的目标（如小型建筑、窄道路）通常会与噪声混淆，所以设30可以保留大部分有效区域的同时滤除大部分椒盐点。如果研究的对象是碎片化景观，比如山区小规模梯田，threshold应降到10-15，否则真实地物会被一并消除。
验证这套后处理是否有效，可以计算处理前后预测结果的OA和Kappa，你会发现Kappa提升比OA更明显——因为连通域过滤主要修正了那些“小面积边界错分”，这些错分在OA里占比很低，但在Kappa的矩阵中显著影响对角线一致性。更进一步，把清理后的栅格和原始影像做半透明叠加，肉眼检查边界是否贴合实际地物轮廓，这一步往往比任何数值指标都更快发现问题。
本文还有配套的精品资源，点击获取
