Python遥感图像分割与国土分类：从SLIC超像素到随机森林实践

简介面向计算机相关专业课程设计与期末大作业的实战项目针对卫星遥感图像国土分类中人工解译效率低的问题基于图像分割技术实现自动化分类涵盖PSPNet、DeepLabV3及DeepLabV3等主流分割模型适合需要完整项目源码进行二次开发、复现或拓展的学习者尤其适合深度学习初学者理解分割网络在遥感场景下的应用。压缩包共12个文件核心为8个Python源码文件覆盖数据加载、图像预处理含水体预处理、多种分割模型构建与训练入口同时附有2个训练日志便于对比验证收敛过程、1份说明文档和1张项目展示图整体仅2.49MB模块划分清晰。项目已经过严格调试下载即可运行省去环境配置与排错时间。已有158人学习下载特别适合正在完成课程设计或需要图像分割实战参考的学生快速上手也可作为遥感分类方向的基线代码扩展使用能显著节省从零搭建模型的时间。1. 一个课程设计题目为什么值得把分割讲清楚拿到“python实现基于图像分割对卫星遥感图像进行国土分类项目源码课程设计.zip”这个标题的同学多半正在做遥感方向的课设或者导师给了一个“用 Python 对遥感影像做地类识别”的题目搜了一圈发现网上要么是深度学习框架的 demo要么是 ArcGIS 里点几下就出结果的流程跟自己要交代码、交报告的需求对不上。这门课设真正难的不是“分类”两个字而是三个隐含问题第一卫星遥感图像和普通照片不一样一个像素可能对应地面上 30 米甚至更宽的区域直接拿像素做分类结果会碎得像撒了一地芝麻第二国土分类耕地、林地、水体、建设用地、裸地对边界的要求很高既要分得对还要分得干净第三作为课程设计你需要一套能讲清楚原理、能现场演示、能写进报告的技术路线而不是黑匣子一样的模型调用。我下面按自己平时做遥感图像处理的一套流程来讲从数据准备、图像分割、特征提取到分类器训练再到踩坑记录和成果整理。整个过程不依赖付费软件用 Python 开源库就能跑通适合课程设计也适合入门遥感图像分类的从业者。2. 国土分类第一步把卫星影像整理成能喂给算法的数据2.1 遥感影像的数据格式与读取方式卫星遥感图像最常见的格式是 GeoTIFF。它和普通 TIFF 的区别在于文件头里嵌入了地理坐标信息包括投影坐标系、像元大小、影像四角的经纬度范围。这意味着你在读取时不仅能拿到像素值还能拿到每个像素对应的真实地理位置后续做面积统计、成果出图都靠这套坐标信息。Python 里读取 GeoTIFF 有两个主流选择GDAL 和 Rasterio。import rasterio from rasterio.plot import show import matplotlib.pyplot as plt # 读取多光谱影像 dataset rasterio.open(landsat_clip.tif) print(波段数:, dataset.count) print(影像尺寸:, dataset.width, x, dataset.height) print(坐标系:, dataset.crs) print(地理范围:, dataset.bounds) # 读取全部波段为一个三维数组 (波段数, 高, 宽) img_array dataset.read() print(数组形状:, img_array.shape) # 用 RGB 三个波段先看一眼影像长什么样 rgb img_array[[3, 2, 1], :, :] # 按 Landsat 常见波段顺序取红绿蓝 show(rgb, transformdataset.transform, cmappink)这段代码做了三件事打开 GeoTIFF 文件、读取影像元信息、把指定波段组合成 RGB 预览图。dataset.read()返回的是(波段数, 高, 宽)的 numpy 数组这是后续所有处理的起点。注意波段索引从 0 开始[[3, 2, 1]]取的是第 4、3、2 波段对应红、绿、蓝波段这样预览图的颜色才接近真实地物。课程设计里经常遇到的一个问题是手头的影像可能包含十几个波段而分类任务只需要其中几个。需要先确认影像包含哪些波段再决定用哪几个组合。如果是为了区分植被和水体近红外波段通常在 Landsat 里是第 5 波段必须参与计算后面提到的 NDVI 指数就要用到它。2.2 影像裁剪、重投影与云遮挡预处理拿到手的卫星影像往往覆盖范围很大动辄上百公里跨度直接整幅跑分类算法内存和耗时都不现实。课程设计只需要一个示范区域所以第一步是裁剪出一个小范围比如一个县、一个镇或者一个水库周边区域。import rasterio from rasterio.windows import from_bounds from rasterio.crs import CRS from rasterio.warp import calculate_default_transform, reproject # 裁剪按经纬度范围切出一块实验区 bbox (116.0, 39.5, 117.0, 40.5) # (西经, 南纬, 东经, 北纬) with rasterio.open(landsat_full.tif) as src: window from_bounds(*bbox, transformsrc.transform) data src.read(windowwindow) profile src.profile.copy() profile.update(heightwindow.height, widthwindow.width, transformwindow.transform) with rasterio.open(landsat_clip.tif, w, **profile) as dst: dst.write(data) # 重投影统一到 UTM 投影方便后续按米计算面积 with rasterio.open(landsat_clip.tif) as src: dst_crs CRS.from_epsg(32650) # UTM 50N按实验区所在经度带选择 transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.profile.copy() kwargs.update(crsdst_crs, transformtransform, widthwidth, heightheight) with rasterio.open(landsat_clip_utm.tif, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_crssrc.crs, dst_crsdst_crs, resamplingrasterio.enums.Resampling.bilinear, )裁剪时用from_bounds把经纬度范围转成像素窗口这一步要求bbox的坐标系和影像本身的坐标系一致。如果影像已经是投影坐标系这里传的就是投影坐标而不是经纬度很多人在这步出错。重投影则要把坐标系从 WGS84经纬度转成 UTM 投影因为经纬度坐标系下每个像素的宽度随纬度变化没法直接算面积。云遮挡是卫星影像特有的问题。如果实验区有云分类结果里云和云影会被分到莫名其妙的类别里。课程设计的应对办法有两个一是选择季度内云量最低的影像下载页面一般有云量信息二是在分类之后把云掩模掉不参与精度统计。2.3 训练样本怎么准备别在分类阶段才想起标签国土分类属于监督分类必须有训练样本。课程设计常见的做法是拿 Google Earth 的高分辨率影像作为参考人工勾选样本区域。注意不要在待分类的影像上直接勾样本因为待分类影像分辨率通常不高30 米或 15 米等间隔取样本点后还要去高分辨率影像上确认地类。样本的组织形式推荐用 GeoJSON 或 Shapefileimport geopandas as gpd import pandas as pd # 人工勾选后得到的样本点/面文件 samples gpd.read_file(train_samples.shp) print(样本数量:, len(samples)) print(类别分布:\n, samples[class_name].value_counts()) # 查看字段确认类别都已经标注 print(samples[[class_name, class_id]].head()) # 检查是否有重叠或空几何避免后续提取像元时出错 assert samples.geometry.notna().all(), 存在空几何 assert samples[class_name].notna().all(), 存在未标注类别国土分类的类别体系建议控制在 5 类以内耕地、林地、水体、建设用地、裸地。类别太多比如把耕地细分成水田和旱地课程设计阶段样本量不够分类器学不好精度反而难看。类别太少只分 3 类报告又显得单薄。5 类是比较稳妥的选择既能展示多分类能力又不会因为样本不足翻车。样本数量上每个类别至少 50 个样本对象如果你做的是面向对象的分类每个样本是一个分割后的多边形每个对象里提取特征这样进入分类器的样本总量在 250 个以上训练一个随机森林分类器就够用了。3. 图像分割算法选型为什么课程设计首选 SLIC 超像素3.1 像素级分类和面向对象分类的本质区别很多第一次做遥感的同学会直接拿影像的每个像素去训练分类器这在课程设计中往往会得到一个“椒盐效果”严重的分类图一块完整的农田被分成几十个小碎片有些像素被分成耕地旁边紧挨着的像素被分成林地。造成这个问题的原因在于单个像素的光谱信息太单薄。一个像素只有几个波段的反射率数值相邻地物的光谱可能非常接近比如干燥的裸土和收割后的耕地在光谱上几乎没有区分度。这时候图像分割的价值就体现出来了先把影像切成一堆“对象”也叫像斑每个对象代表一块均质区域再用对象的均值、方差、纹理、几何特征去分类噪声被抑制了分类结果也更符合人对地物的认知。图像分割在国土分类中的角色用一句话概括不是分割本身要得到地类是分割把影像从“像素网格”转成“地物对象”为后续特征提取和分类提供载体。3.2 三种常见分割算法对比与课程设计选型Python 里能用的分割算法不少我列一下实际比较过的情况算法典型库特点适合场景SLIC 超像素scikit-image速度快、参数少、控制超像素数量直接课程设计首选分水岭OpenCV / scikit-image对边缘敏感需预处理梯度边缘清晰、地物边界分明的影像多尺度分割需要遥感专用库如 RSGISLib可分层分割效果好但配置复杂毕业设计级别时间充裕再说课程设计我推荐 SLIC原因是它只有一个核心参数n_segments控制整幅影像被切成的块数调参直觉非常直观块数越多对象越小块数越少对象越大。国土分类关注的是较大尺度的地物农田、水体、居民区不需要把单棵树都切出来所以n_segments取一个适中偏小的数量就行。分水岭算法更适合做边缘检测后处理它要先算影像梯度梯度图中“分水岭脊线”就是地物边界。这个方法对细小地物效果好但参数多梯度图的质量直接决定分割结果的好坏而梯度图怎么算又有一堆讲究。课程设计阶段不建议在这个上面花时间。3.3 用 SLIC 分割遥感影像的完整步骤分割的输入是前面准备的多光谱影像输出是一个和原始影像同尺寸的标签数组每个像素的值代表它属于哪个超像素块。from skimage.segmentation import slic, mark_boundaries from skimage.color import label2rgb import numpy as np # 读取预处理后的影像取需要的波段 with rasterio.open(landsat_clip_utm.tif) as src: img src.read() # (波段数, 高, 宽) # 转成 (高, 宽, 波段数) 并做标准化SLIC 对量纲敏感 img_hwc np.transpose(img, (1, 2, 0)).astype(np.float32) for b in range(img_hwc.shape[2]): band img_hwc[:, :, b] img_hwc[:, :, b] (band - band.min()) / (band.max() - band.min() 1e-6) # 执行 SLIC 分割 segments slic( img_hwc, n_segments2000, compactness10, sigma1, start_label1, ) print(分割出的对象数量:, segments.max())compactness是另一个关键参数。它控制分割结果的“方正程度”值越小分割越贴合地物边界对象形状越不规则值越大分割结果越接近规则方块地物边界会被裁直。对国土分类来说一般取 520 之间我习惯先取 10 看效果哪类地物边界碎了就调低哪类地物边界糊了比如水体边缘锯齿消失就调高。sigma是对影像做高斯预平滑的强度取 1 即可作用是抑制传感器噪声导致的细小分割块。n_segments2000这个数字是估的你可以根据影像尺寸调整一个常见经验是让平均每个对象的面积在 500 到 2000 个像素之间。分割完成后一定要目视检查。把分割边界叠加在原图上看看农田是否被切成一整块道路是否被单独分出来水体边界是否贴合河岸。这一步不能省分割质量决定了后面分类精度的上限。如果发现分割块太大精细地物小水塘、窄道路被并进了周围地物就增大n_segments如果分割块太碎同一种地物被切成几十块就减小。4. 特征提取与国土分类模型让分割后的对象开口说话4.1 光谱特征、纹理特征与几何特征每一类都要有物理意义分割完成之后每个对象超像素块就代替像素成为基本分析单元。接下来要做的是对每个对象提取特征向量然后把这个特征向量输入分类器。常用的特征分三类光谱特征每个波段内像元的均值、中位数、标准差。均值代表对象的光谱平均值标准差代表对象内部光谱均匀程度。农田和裸地可能均值接近但农田内部的种植行距会带来更高的纹理变化。指数特征NDVI归一化植被指数、NDWI归一化水体指数。NDVI 能把植被从土壤和建筑中分离出来一个对象如果在生长季 NDVI 均值超过 0.4大概率是农田而不是裸地。几何特征对象面积、周长、形状指数周长除以面积平方根。建设用地形状多为规则多边形自然林地形状不规则这个特征能有效区分。特征不是越多越好。课程设计的样本量通常只有几百个特征如果堆到几十维分类器会学到噪声。我建议控制在 10 个左右5 个波段均值 2 个指数 1 个标准差 2 个几何特征。from skimage.measure import regionprops import numpy as np # 对每个分割对象提取特征 def extract_features(img_hwc, segments, band_indices[3, 4, 5]): img_hwc: (高, 宽, 波段数) 的影像数组 segments: SLIC 分割结果与 img_hwc 同尺寸 band_indices: 要参与特征提取的波段索引 props regionprops(segments, intensity_imageimg_hwc[:, :, band_indices[0]]) features [] for prop in props: # 取该对象区域内的像素掩膜 mask segments prop.label # 光谱均值 spectral_means [] for b in band_indices: band_data img_hwc[:, :, b] spectral_means.append(band_data[mask].mean()) # 计算 NDVI假设第 4 波段近红外、第 3 波段红 nir img_hwc[:, :, 4][mask].mean() red img_hwc[:, :, 3][mask].mean() ndvi (nir - red) / (nir red 1e-6) # 几何特征面积、形状指数 area prop.area perimeter prop.perimeter shape_index perimeter / (2 * np.sqrt(np.pi * area) 1e-6) features.append(spectral_means [ndvi, shape_index, area]) return np.array(features) features extract_features(img_hwc, segments) print(特征矩阵形状:, features.shape) # (对象数量, 特征维度)regionprops是 scikit-image 里做对象属性提取的工具直接用分割标签数组就能拿到每个对象的面积、周长、惯性矩、强度统计等属性。这里的代码把光谱均值和几何特征手动拼装是为了后续加特征时思路清晰也方便在报告里列出每个特征的物理含义。这段代码里唯一需要调整的是band_indices和 NDVI 计算中用到的波段位置。不同传感器的波段定义不同写代码前先查清楚影像每个波段的波长再用真实配置改进。4.2 随机森林分类器训练与参数调节特征向量准备好后需要一个分类器。遥感领域做中低分辨率影像分类随机森林是稳定可靠的选择。它对特征量纲不敏感特征是均值、面积、指数混合在一起取值范围差异大随机森林不需要做标准化能处理多分类问题并且能输出特征重要性这正好用在报告里展示“哪个特征对分类贡献最大”。from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report, confusion_matrix # 样本特征和标签已经人工标注完成 X features # 每个训练对象的特征向量 y labels # 与特征对应的地类标签 # 划分训练集和验证集 X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.3, stratifyy, random_state42 ) # 训练随机森林分类器 rf RandomForestClassifier( n_estimators200, max_depth10, min_samples_leaf2, random_state42, ) rf.fit(X_train, y_train) # 在验证集上评估 y_pred rf.predict(X_test) print(classification_report(y_test, y_pred, target_names[耕地, 林地, 水体, 建设用地, 裸地]))随机森林的三个参数值得认真解释。n_estimators是树的数量。过小比如 50模型不稳定预测结果随随机种子变化大过大比如 500 以上训练时间变长边际精度提升很小。200 是一个性价比不错的默认值你可以在 100300 之间试。max_depth控制每棵树的深度。遥感影像特征之间往往存在相关性树太深会过拟合训练样本导致验证集精度虚高但实际应用翻车。课程设计数据量小深度设 10 左右能保留模型泛化能力。min_samples_leaf是指定叶子节点的最小样本数。设 2 的含义是如果某个叶子节点的样本少于 2 个就停止分裂。这个参数能防止分类器记住个别样本的噪声。数据量小的时候尤其重要。4.3 精度评估的标准操作混淆矩阵读法课程设计报告里不能只写“分类精度 85%”导师一定会问每个类别的精度分别是多少哪些类别容易混混淆矩阵是标准的答案。import matplotlib.pyplot as plt import seaborn as sns # 计算混淆矩阵 cm confusion_matrix(y_test, y_pred) # 转换为百分比便于报告展示 cm_percent cm.astype(float) / cm.sum(axis1, keepdimsTrue) * 100 fig, ax plt.subplots(figsize(8, 6)) sns.heatmap( cm_percent, annotTrue, fmt.1f, cmapBlues, xticklabels[耕地, 林地, 水体, 建设用地, 裸地], yticklabels[耕地, 林地, 水体, 建设用地, 裸地], axax, ) ax.set_xlabel(预测类别) ax.set_ylabel(真实类别) plt.tight_layout() plt.savefig(confusion_matrix.png, dpi200)混淆矩阵的行代表真实类别列代表预测类别。对角线是每个类别的分类正确率召回率。观察非对角线元素你会发现耕地和裸地经常互相混光谱相似建设用地和林地偶尔混阴影区域。这些信息都值得写进报告说明分类结果在哪类地物上不可靠以及改进方向是什么。补充一个常见问题整体精度是用验证集算的验证集必须和训练集独立。如果你用训练集回测精度会虚高在报告里会显得不专业。5. 国土分类避坑课程设计阶段最容易翻车的四个问题5.1 环境安装GDAL 装不上是初学者第一道坎现象pip install rasterio装到一半报错出现一堆红色输出常见的是Failed building wheel for GDAL或Microsoft Visual C 14.0 is required。原因rasterio、fiona、geopandas 这类库依赖底层的 GDAL 库GDAL 原生代码是用 C/C 写的pip 安装时需要在本地编译而 Windows 环境下缺少编译工具链就会直接失败。解决Windows 用户最稳妥的办法是使用 conda 安装地理空间库conda 会直接分发预编译的二进制包不需要本地编译。提示用 conda 创建独立环境不要往基础环境里装遥感库的依赖版本要求比较严格独立环境能避免库冲突导致的项目崩溃。conda create -n landcover python3.10 -y conda activate landcover conda install -c conda-forge rasterio geopandas scikit-image scikit-learn matplotlib如果只能用 pip先访问 PyPI 或第三方镜像站下载对应 Python 版本的预编译 wheel 文件再本地安装。注意 wheel 文件名里的cp310表示适配 Python 3.10版本不匹配会继续报错。5.2 中文路径读不了数据现象rasterio 打开以中文命名的 tif 文件时抛出file does not exist in the file system错误但文件明明就在那里。原因部分 GDAL 版本在 Windows 下处理中文路径时出现编码问题本质是底层 C 库的文件路径编码与 Python 字符串编码不一致。解决把整个项目路径全部改成英文字符包括文件名、文件夹名、工程目录路径。这是一劳永逸的办法。如果项目已经用了中文路径先复制一份所有数据到英文路径下再继续。5.3 分割结果过碎或过整分类精度都不高现象SLIC 分割后有些地类被切得稀碎分类器在这些碎片上的预测结果几乎完全错误另一些地类又被合并成一大块里面混了道路和房屋。原因n_segments或compactness设置有偏差。过碎是n_segments太大过整是compactness太大。解决不要一次只调一个参数而是交替观察效果。先固定compactness10从小到大试n_segments找到对象大小合适后再微调compactness让边界更贴合。每次调参后都把分割边界叠加到原图上目视检查不要把参数交给模型自动寻优课程设计阶段的影像数量不多手工调参是可控的。5.4 训练样本和验证样本不分家精度虚高现象验证集精度 90% 以上换一块区域应用分类模型结果完全不能看。原因训练样本和验证样本来自同一块区域甚至同一块农田上取了相邻的样本。由于空间自相关性验证样本和训练样本的光谱特征几乎一样模型相当于开卷考试。解决划分训练集和验证集时要按“区域”划分而不是按“单个样本”划分。比如用影像的东部区域做训练样本西部区域做验证样本。这样测试出的精度才接近模型的真实泛化能力。课程设计数据量有限时这个做法的代价是训练样本少但报告的说服力强得多。6. 把分割结果的矢量化和面积统计做成报告的加分项课程设计交了分类结果图通常还不够导师喜欢看到一个“能算数”的结果。比如分类完成后统计各土地类型的面积这在国土分类业务中是最常见的需求。做法是把分割分类后的栅格结果矢量化导出为 Shapefile再用 GeoPandas 做面积统计。from skimage.measure import find_contours import geopandas as gpd from shapely.geometry import Polygon import rasterio.features # 假设 final_label 是分类结果数组每个像素的值代表地类编码 # 将栅格分类结果转为矢量面 results [] with rasterio.open(landsat_clip_utm.tif) as src: transform src.transform for class_id in [1, 2, 3, 4, 5]: mask (final_label class_id).astype(uint8) for polygon, value in rasterio.features.shapes( mask, transformtransform ): results.append({ class_id: class_id, geometry: Polygon(polygon[coordinates][0]), }) # 转成 GeoDataFrame 并计算面积 gdf gpd.GeoDataFrame(results, crsEPSG:32650) gdf[area_m2] gdf.geometry.area area_stats gdf.groupby(class_id)[area_m2].sum() / 1e6 # 转为平方公里 print(area_stats)rasterio.features.shapes把每个连通区域转成一个矢量多边形同时把地类编码作为属性保留下来。矢量化之后面积计算就变得简单可靠——因为坐标系已经是 UTM 投影geometry.area直接得到平方米。最后一步是对面积统计做一个可视化画柱状图比直接贴数字直观得多。按面积从大到小排序后用不同颜色画柱状图每根柱子上面标注面积数值和占比。这一张图放进课程设计报告里比任何文字都有说服力。矢量化之后还能做一件有意思的事把矢量结果叠加到原始影像上生成对比图左半部分显示原始影像右半部分显示分类对象边界。这种图在答辩时特别受欢迎因为一眼就能看出分割对象和真实地物边界的关系也能直观展示分类的细碎程度。我做这个方向时养成的习惯是每一次参数调完先把结果图保存下来按日期和参数命名。这避免了反复对比时找不到之前的版本。课程设计的调试过程中参数要调很多轮不保存记录很容易调着调着就调糊了最后想回退都没法回退。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取