简介针对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_es1.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, axis0)参数说明这里用的是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、短波红外1B6。这6个波段覆盖了植被、水体、土壤、建筑区分最敏感的谱段范围加入B7SWIR2收益有限但会显著增加计算量。真彩色合成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), keepdimsTrue) std stacked.std(axis(1,2), keepdimsTrue) 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( shapesshapes, out_shapeout_shape, transformtransform, fill0, all_touchedFalse, # 只填充矢量完全覆盖的像素 dtypenp.uint8 ) # 写出 with rio.open(out_label_path, w, driverGTiff, heightout_shape[0], widthout_shape[1], count1, dtypenp.uint8, crssrc.crs, transformtransform) 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_ratio0.7, val_ratio0.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, padding1), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue), nn.Conv2d(out_ch, out_ch, 3, padding1), nn.BatchNorm2d(out_ch), nn.ReLU(inplaceTrue) ) def forward(self, x): return self.conv(x) class UNet(nn.Module): def __init__(self, in_channels6, num_classes5): 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, stride2) self.dec2 ConvBlock(512, 256) self.up3 nn.ConvTranspose2d(256, 128, 2, stride2) self.dec3 ConvBlock(256, 128) self.up4 nn.ConvTranspose2d(128, 64, 2, stride2) 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], dim1)) d3 self.dec3(torch.cat([self.up3(d2), e2], dim1)) d4 self.dec4(torch.cat([self.up4(d3), e1], dim1)) return self.out(d4)逻辑说明编码器每层输出通道翻倍空间尺寸减半解码器做转置卷积上采样后与对应编码器层拼接通道数相加。这里的BatchNorm对Landsat输入特别重要——各波段缩放后的数值分布差异极大BatchNorm能强制逐通道归一化加速收敛。注意卷积核固定为3×3padding为1保证尺寸不变像素级分类不涉及全连接层所以输入尺寸可以不是固定值。3.3 损失函数与类别不均衡小样本分类的胜负手地物分类中类别不均衡是常态。水体、裸地往往占据影像70%以上面积而建筑、道路等目标类可能只占3%。直接用标准交叉熵损失模型会把所有像素预测为多数类OAOverall Accuracy看着高实际毫无意义。常见的做法是给交叉熵加权重权重与类别频率成反比def class_weights_from_label(label, num_classes): # 统计每个类别的像素占比 counts np.bincount(label.flatten(), minlengthnum_classes) total counts.sum() # 权重 log(总像素 / (类别像素 * 类别数))防止权重极端 weights np.log(total / (counts * num_classes) 1e-8) return torch.tensor(weights, dtypetorch.float32) criterion nn.CrossEntropyLoss(weightclass_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 Size16128×128输入Landsat输入通道数少显存占用不大Epoch50配合早停patience10数据增强随机旋转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, window128, stride96, devicecuda): model.eval() h, w full_img.shape[1], full_img.shape[2] pred_sum torch.zeros((num_classes, h, w), devicedevice) count torch.zeros((1, h, w), devicedevice) for y in range(0, h - window 1, stride): for x in range(0, w - window 1, stride): patch full_img[:, y:ywindow, x:xwindow].unsqueeze(0).to(device) with torch.no_grad(): out model(patch) # (1, C, window, window) prob torch.softmax(out, dim1) pred_sum[:, y:ywindow, x:xwindow] prob[0] count[:, y:ywindow, x:xwindow] 1 # 取平均概率后取argmax得到类别 prob_avg pred_sum / count.clamp(min1) label_map prob_avg.argmax(dim0).cpu().numpy() return label_map参数说明stride96表示窗口之间重叠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_index0): # 只统计有效像素忽略背景 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 阴影山体阴影在可见光段与水体的低反射率相似引入SWIR1B6水体在1.6μm处吸收强烈而阴影不会裸地 vs 建筑屋顶屋顶材料混凝土/沥青光谱与裸地重叠加入纹理信息或者用形态学开闭运算后处理多数深度学习项目忽略了一个事实Landsat的30米分辨率意味着一个像素往往是混合地物例如林地与灌木的过渡带。与其要求模型做精细化边界不如在后处理中引入“不确定度”概念——用预测概率的熵来标记那些处于类别边界的像素供人工检查或作为矢量化的判定条件。5. 进阶技巧用连通域分析清除预测结果的椒盐噪声模型逐像素预测的输出像质粗糙常见问题是孤立的小块噪声区域——单个或几个像素被错误分类在图像上表现为“椒盐”颗粒。这类噪声在均匀地块内部尤其显眼直接影响矢量化和面积统计的精度。一个不依赖额外模型的后处理技巧是连通域过滤给预测标签图做连通域标记统计每个连通域的面积像素数将面积小于阈值的连通域替换为周围最频繁的类别。from scipy import ndimage import numpy as np def remove_salt_pepper(label_map, min_area30): 删除小于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, structurenp.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, modeedge) for idx in np.argwhere(neg_mask): y, x idx neighborhood pad[y:y3, x:x3].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邻域投票修复。thresholdmin_area的选择有门道Landsat 30米分辨率下3×3像素对应90×90米区域地物分类中低于该面积的目标如小型建筑、窄道路通常会与噪声混淆所以设30可以保留大部分有效区域的同时滤除大部分椒盐点。如果研究的对象是碎片化景观比如山区小规模梯田threshold应降到10-15否则真实地物会被一并消除。验证这套后处理是否有效可以计算处理前后预测结果的OA和Kappa你会发现Kappa提升比OA更明显——因为连通域过滤主要修正了那些“小面积边界错分”这些错分在OA里占比很低但在Kappa的矩阵中显著影响对角线一致性。更进一步把清理后的栅格和原始影像做半透明叠加肉眼检查边界是否贴合实际地物轮廓这一步往往比任何数值指标都更快发现问题。本文还有配套的精品资源点击获取
