Landsat多光谱影像深度学习分类实战：从预处理到GeoTIFF输出

简介本资源是一套面向遥感图像处理研究者与深度学习初学者的实用工具包聚焦Landsat影像地物分类任务提供端到端的CNN建模与推理能力。资源包含10个文件3个Python脚本、2个TIFF遥感影像、2个XML元数据、1个H5训练模型、1个README说明及1个TFW地理配准文件总大小14.49MB其中trainModel.py和predictNewData.py构成训练-预测完整流程CNN_7class_3by3.h5为已训练好的7类地物分类模型example.tif与new_class.tif为示例影像数据便于快速验证与迁移应用。已有181人学习下载适合高校遥感/地信专业师生开展课程实验、科研原型开发或竞赛项目复现。用户可直接加载H5模型对新Landsat影像进行推理无需GPU环境即可完成地物分类水体、植被、建筑等配套代码结构清晰、注释完整并涵盖影像切片、数据增强、模型训练与结果可视化等关键环节显著降低深度学习在遥感领域的入门门槛。1. 用CNN对Landsat影像做地物分类不是调个库就能出图的事你手头有一套Landsat 8/9 的 Level-2 地表反射率数据.tif格式想自动区分水体、林地、农田、裸地、建成区——但直接扔进scikit-learn的RandomForest精度卡在72%上不去换ResNet50微调又发现输入尺寸不匹配、波段数对不上、训练时显存爆掉。这不是模型不行而是遥感影像的物理特性没被深度学习管道真正“吃透”Landsat有7个有效波段B2–B7空间分辨率30m辐射定标和大气校正已内置但像素值范围是0–65535整型不能像ImageNet图片那样直接归一化到[0,1]更关键的是单景影像动辄5000×5000像素全图送入CNN既不可行也不必要。真实落地的做法是把“遥感先验知识”编进数据流水线——比如用NDVI增强植被响应、按地理格网切块时保留邻域上下文、对多光谱通道做带权重的标准化。本文聚焦Python生态下可复现的端到端流程从原始Landsat .tar.gz解压开始到生成带地理坐标的分类GeoTIFF全程使用PyTorchRasterioGDAL不依赖ArcGIS或ENVI商业模块。适合已有遥感数据但缺乏深度学习工程经验的测绘、生态、地信从业者也适合作为高校遥感AI课程的实操基线。2. Landsat多光谱数据预处理为什么必须重写标准化与切块逻辑遥感影像的预处理不是图像处理的简单平移。Landsat Level-2产品虽已做大气校正但其DN值Digital Number仍需转换为物理量级的地表反射率0–1浮点且各波段动态范围差异极大近红外B5常达0.4以上而海岸带波段B1普遍低于0.1。若直接min-max归一化低反射率波段信息将被压缩至数值噪声层若统一除以65535又会丢失波段间辐射响应的相对关系。正确做法是按USGS官方文档推荐的缩放因子scale factor 0.0000275, offset -0.2进行线性变换并对每个波段单独计算统计量。2.1 解压与读取Landsat场景文件Landsat数据以.tar.gz压缩包分发内含多个波段TIFF及元数据MTL文件。需解析MTL提取辐射定标参数而非硬编码波段顺序import rasterio from pathlib import Path import numpy as np def parse_landsat_mtl(mtl_path): 从MTL文件提取辐射定标系数 with open(mtl_path, r) as f: lines f.readlines() params {} for line in lines: if RADIANCE_MULT_BAND_ in line or REFLECTANCE_MULT_BAND_ in line: key, val line.strip().split( ) params[key.strip()] float(val.strip().strip()) elif RADIANCE_ADD_BAND_ in line or REFLECTANCE_ADD_BAND_ in line: key, val line.strip().split( ) params[key.strip()] float(val.strip().strip()) return params # 示例读取B2-B7共6个波段Landsat 8去除了热红外B10/B11用于地物分类 bands_of_interest [B2, B3, B4, B5, B6, B7] # 对应蓝、绿、红、近红外、短波红外1、短波红外2 band_paths [Path(LC08_L2SP_123032_20220515_20220520_02_T1) / fLC08_L2SP_123032_20220515_20220520_02_T1_SR_{b}.TIF for b in bands_of_interest] mtl_path Path(LC08_L2SP_123032_20220515_20220520_02_T1) / LC08_L2SP_123032_20220515_20220520_02_T1_MTL.txt提示Landsat 9与8的波段编号一致但Landsat 7需跳过失效的B6热红外并注意B8全色不参与地物分类。务必通过MTL文件确认PRODUCT_CONTENTS字段避免波段错位。2.2 多光谱反射率计算与波段堆叠Landsat Level-2提供表面反射率SR产品但值域为0–10000整型需除以10000转为浮点反射率并应用USGS推荐的太阳天顶角校正cosθs提升跨季节一致性def read_and_calibrate_band(tif_path, mtl_params, sun_zenith_deg25.3): # 示例天顶角实际从MTL读取 with rasterio.open(tif_path) as src: band_data src.read(1).astype(np.float32) # Level-2 SR数据已做大气校正直接缩放官方说明中scale0.0001, offset0 reflectance band_data * 0.0001 # 太阳天顶角校正补偿不同采集时间的光照入射角差异 cos_theta_s np.cos(np.radians(sun_zenith_deg)) calibrated reflectance / cos_theta_s if cos_theta_s 0 else reflectance return calibrated, src.transform, src.crs # 堆叠6个波段为 (C, H, W) 张量 stacked_bands [] transforms, crs None, None for i, bpath in enumerate(band_paths): band_arr, transform, crs read_and_calibrate_band(bpath, mtl_params) stacked_bands.append(band_arr) if i 0: transforms, crs transform, crs X np.stack(stacked_bands, axis0) # shape: (6, 5000, 5000)2.3 按地理格网切块并保留空间上下文全图训练不现实但随机裁剪会破坏地物的空间连续性如道路、河流。采用滑动窗口重叠策略确保每个patch包含足够邻域信息def tile_with_overlap(array_3d, tile_size256, overlap32): 按重叠方式切块返回 (N, C, H, W) 和地理坐标列表 c, h, w array_3d.shape tiles, geoms [], [] # 计算起始行列避免边界不足 for i in range(0, h - tile_size 1, tile_size - overlap): for j in range(0, w - tile_size 1, tile_size - overlap): tile array_3d[:, i:itile_size, j:jtile_size] tiles.append(tile) # 计算该tile左上角地理坐标 left, top rasterio.transform.xy(transforms, i, j) geoms.append((left, top, left tile_size * transforms.a, top tile_size * transforms.e)) return np.array(tiles), geoms tiles, tile_geoms tile_with_overlap(X, tile_size256, overlap32) print(f原始影像尺寸: {X.shape} → 切得 {len(tiles)} 个256×256切片) # 输出: 原始影像尺寸: (6, 5000, 5000) → 切得 324 个256×256切片注意overlap32意味着相邻切片有32像素重叠推理时对重叠区域取平均可显著减少块效应blocky artifacts这是遥感分类区别于自然图像的关键技巧。3. 构建适配Landsat的CNN主干为什么标准ResNet要改输入层与初始化Landsat的6波段输入与ImageNet的3通道RGB存在本质差异蓝绿红近红外构成的“假彩色”组合无法直接套用RGB预训练权重短波红外波段B6/B7对土壤湿度、矿物成分敏感但其数值分布与可见光波段相差一个数量级。强行加载ImageNet权重会导致前几层梯度爆炸而从零训练又需要海量标注数据。折中方案是修改ResNet第一卷积层并采用遥感领域验证有效的初始化策略。3.1 自定义ResNet输入适配器PyTorch官方ResNet默认输入为3通道需替换conv1层并调整权重初始化import torch import torch.nn as nn from torchvision.models import resnet50 def build_landsat_resnet(num_bands6, num_classes5, pretrainedFalse): model resnet50(pretrainedpretrained) # 替换第一层卷积3→6通道保持kernel_size7,stride2,padding3不变 model.conv1 nn.Conv2d( in_channelsnum_bands, out_channels64, kernel_size7, stride2, padding3, biasFalse ) # 初始化新卷积层采用MSRA初始化适用于ReLU激活 nn.init.kaiming_normal_(model.conv1.weight, modefan_out, nonlinearityrelu) # 修改全连接层输出维度 model.fc nn.Linear(model.fc.in_features, num_classes) # 若启用pretrained仅初始化新层否则全部随机初始化 if not pretrained: for m in model.modules(): if isinstance(m, nn.Conv2d): nn.init.kaiming_normal_(m.weight, modefan_out, nonlinearityrelu) elif isinstance(m, nn.BatchNorm2d): nn.init.constant_(m.weight, 1) nn.init.constant_(m.bias, 0) return model # 实例化模型 model build_landsat_resnet(num_bands6, num_classes5, pretrainedTrue) print(f模型输入通道: {model.conv1.in_channels}, 输出类别: {model.fc.out_features}) # 输出: 模型输入通道: 6, 输出类别: 53.2 针对多光谱的BatchNorm优化标准BatchNorm在小batch size下不稳定而遥感切片训练常受限于显存256×256×6≈3MB/张。采用GroupNorm替代对通道分组归一化不受batch size影响def replace_bn_with_gn(model, num_groups4): 递归替换所有BatchNorm2d为GroupNorm for name, module in model.named_children(): if isinstance(module, nn.BatchNorm2d): model._modules[name] nn.GroupNorm( num_groupsnum_groups, num_channelsmodule.num_features ) else: replace_bn_with_gn(module, num_groups) return model model replace_bn_with_gn(model, num_groups4) # 验证替换结果 for m in model.modules(): if isinstance(m, nn.GroupNorm): print(f已替换为GroupNorm分组数: {m.num_groups}) break3.3 损失函数选择Focal Loss缓解地物样本不均衡遥感影像中水体、建成区常占比较小5%交叉熵损失易偏向多数类。Focal Loss通过调节难易样本权重提升少数类召回率class FocalLoss(nn.Module): def __init__(self, alpha1, gamma2, reductionmean): super().__init__() self.alpha alpha self.gamma gamma self.reduction reduction def forward(self, inputs, targets): ce_loss F.cross_entropy(inputs, targets, reductionnone) pt torch.exp(-ce_loss) focal_weight (1 - pt) ** self.gamma loss self.alpha * focal_weight * ce_loss if self.reduction mean: return loss.mean() elif self.reduction sum: return loss.sum() else: return loss criterion FocalLoss(alpha1, gamma2)参数说明gamma2是常用值增大gamma会进一步抑制易分类样本贡献alpha可设为各类别的逆频率如水体占比0.03则alpha_water1/0.03≈33但实践中常统一设为1并调gamma。4. 训练与验证如何用有限标注数据获得高精度分类结果遥感地物分类的最大瓶颈不是模型而是高质量标注数据稀缺。一套覆盖5类地物的Landsat切片人工勾绘1000张256×256样本需200小时以上。因此必须结合半监督与数据增强策略在有限标注下榨取最大性能。4.1 半监督训练Mean Teacher框架稳定教师-学生模型Mean Teacher利用无标签数据提升泛化能力学生模型参数实时更新教师模型参数由学生指数移动平均EMA得到二者预测一致性作为额外监督信号class MeanTeacherLoss(nn.Module): def __init__(self, consistency_weight1.0, consistency_rampup100): super().__init__() self.consistency_weight consistency_weight self.consistency_rampup consistency_rampup def forward(self, student_logits, teacher_logits, epoch): # 计算KL散度一致性损失 student_prob F.softmax(student_logits, dim1) teacher_prob F.softmax(teacher_logits, dim1) consistency_loss F.kl_div( torch.log(student_prob 1e-8), teacher_prob, reductionbatchmean ) # ramp-up权重前100轮线性增长至设定值 weight self.consistency_weight * sigmoid_rampup(epoch, self.consistency_rampup) return weight * consistency_loss def sigmoid_rampup(current, rampup_length): if rampup_length 0: return 1.0 current np.clip(current, 0, rampup_length) phase 1.0 - current / rampup_length return float(np.exp(-5.0 * phase * phase)) # 在训练循环中调用 consistency_criterion MeanTeacherLoss(consistency_weight1.0) consistency_loss consistency_criterion(student_output, teacher_output, epoch) total_loss supervised_loss consistency_loss4.2 遥感专用数据增强组合自然图像增强如水平翻转、色彩抖动对遥感无效——NDVI值会因波段顺序错乱而崩溃。必须设计光谱安全spectral-safe增强import albumentations as A from albumentations.pytorch import ToTensorV2 # 仅对空间维度操作保持波段间关系 train_transform A.Compose([ A.RandomRotate90(p0.5), A.HorizontalFlip(p0.5), A.VerticalFlip(p0.5), A.RandomScale(scale_limit0.1, p0.5), # 缩放不影响光谱值 A.GaussNoise(var_limit(0.001, 0.005), p0.3), # 添加微弱噪声模拟传感器误差 ToTensorV2(), ]) # 验证时仅归一化 val_transform A.Compose([ ToTensorV2(), ])关键约束禁用所有涉及颜色空间转换的增强如CLAHE,HueSaturationValue,RandomBrightnessContrast这些会破坏波段物理意义。NDVI (B5-B4)/(B5B4) 必须在增强前后保持数学一致性。4.3 验证指标混淆矩阵与Kappa系数必须地理加权Accuracy在遥感中具有欺骗性——若水体仅占1%模型全判为非水体也能得99%准确率。必须计算地理加权的IoU与Kappadef compute_geographic_iou(pred_mask, gt_mask, weightsNone): 按像素地理面积加权IoU考虑不同纬度像素面积差异 if weights is None: weights np.ones_like(pred_mask) iou_per_class [] for cls in np.unique(gt_mask): pred_cls (pred_mask cls) gt_cls (gt_mask cls) intersection np.sum((pred_cls gt_cls) * weights) union np.sum((pred_cls | gt_cls) * weights) iou intersection / (union 1e-6) iou_per_class.append(iou) return np.array(iou_per_class) # 计算Kappa系数需混淆矩阵 from sklearn.metrics import cohen_kappa_score kappa cohen_kappa_score(gt_flat, pred_flat) print(fKappa系数: {kappa:.3f} (0.8表示极强一致性))5. 推理与后处理从模型输出到可交付的GeoTIFF分类图训练完成的模型输出是(N, C, H, W)概率图但业务系统需要带地理坐标、符合GIS标准的GeoTIFF。这要求精确重建空间参考并对概率图做地理空间后处理如形态学滤波消除椒盐噪声。5.1 批量推理与地理坐标映射使用滑动窗口推理时需将每个切片的预测结果按原始地理坐标拼接回全图def predict_full_scene(model, tile_loader, tile_geoms, original_shape, device): 对全图切片推理并重建地理空间数组 model.eval() c, h, w original_shape # 初始化全图预测数组存储类别索引 full_pred np.zeros((h, w), dtypenp.uint8) # 初始化计数数组用于重叠区域平均 count_map np.zeros((h, w), dtypenp.uint8) with torch.no_grad(): for i, (tiles, _) in enumerate(tile_loader): tiles tiles.to(device) outputs model(tiles) # (B, C, H, W) preds torch.argmax(outputs, dim1).cpu().numpy() # (B, H, W) # 将第i个切片写入full_pred对应位置 tile_h, tile_w preds.shape[1:] left, top, right, bottom tile_geoms[i] # 根据transform反算像素行列号简化版实际需用rasterio.transform.rowcol # 此处假设tile_geoms已转换为像素坐标 start_row, start_col int((top - transform.f) / transform.e), int((left - transform.c) / transform.a) full_pred[start_row:start_rowtile_h, start_col:start_coltile_w] preds[0] count_map[start_row:start_rowtile_h, start_col:start_coltile_w] 1 # 取平均重叠区域 full_pred (full_pred / count_map).astype(np.uint8) return full_pred # 调用推理 full_prediction predict_full_scene(model, tile_loader, tile_geoms, X.shape, devicecuda)5.2 形态学后处理与矢量化导出概率图直接argmax会产生锯齿状边界需用地理感知的形态学操作平滑from scipy import ndimage import geopandas as gpd from shapely.geometry import Polygon def morphological_postprocess(pred_array, kernel_size3): 对分类图做开运算消除椒盐噪声 # 对每一类分别做形态学操作 processed np.zeros_like(pred_array) for cls_id in np.unique(pred_array): if cls_id 0: # 背景类跳过 continue binary_mask (pred_array cls_id) # 开运算先腐蚀再膨胀去除小噪点 kernel np.ones((kernel_size, kernel_size), dtypenp.uint8) opened ndimage.binary_opening(binary_mask, structurekernel) processed[opened] cls_id return processed smoothed_pred morphological_postprocess(full_prediction, kernel_size5) # 导出为GeoTIFF复用原始影像的transform和crs with rasterio.open( landsat_classification.tif, w, driverGTiff, heightsmoothed_pred.shape[0], widthsmoothed_pred.shape[1], count1, dtypesmoothed_pred.dtype, crscrs, transformtransforms ) as dst: dst.write(smoothed_pred, 1)5.3 矢量化地类图斑并导出为ShapefileGIS业务常需矢量图层使用rasterio.features.shapes直接栅格转矢量def raster_to_vector(raster_path, output_shp, class_namesNone): 将分类GeoTIFF转为带属性的Shapefile with rasterio.open(raster_path) as src: image src.read(1) mask image ! 0 # 掩膜非背景区域 # 生成矢量多边形 shapes list(rasterio.features.shapes(image, maskmask, transformsrc.transform)) # 构建GeoDataFrame polygons, values zip(*shapes) gdf gpd.GeoDataFrame({ geometry: [Polygon(poly) for poly in polygons], class_id: [int(v) for v in values] }, crssrc.crs) # 添加类别名称字段 if class_names: gdf[class_name] gdf[class_id].map(lambda x: class_names.get(x, unknown)) gdf.to_file(output_shp, driverESRI Shapefile) # 调用 class_mapping {0: Background, 1: Water, 2: Forest, 3: Cropland, 4: Bare, 5: Built-up} raster_to_vector(landsat_classification.tif, landcover.shp, class_mapping)验证技巧用QGIS打开生成的landcover.shp叠加原始Landsat真彩色影像B4-B3-B2合成目视检查道路、田埂、河岸线是否连续。若出现断裂需调大morphological_postprocess的kernel_size或增加训练epoch。本文还有配套的精品资源点击获取