处理过遥感影像的同学应该都有体会项目组或者导师递过来一张高分辨率遥感影像第一句话大概率是“把这几个地类给我分出来”。以前我的第一反应是打开 GIS 软件手动勾绘或者直接架起深度学习模型慢慢做语义分割。但有一次项目周期只有三天既没有现成的标注样本也没有 GPU我改用 k-Means 聚类算法对遥感影像做语义分割结果出图效果意外地能打。这篇文章就完整记录我怎么用 scikit-learn 里的 k-Means 算法把一张多光谱遥感影像快速分成植被、水体、建筑、裸地并最终在 ArcGIS Pro 里生成地类图斑矢量图层的全过程。适合遥感、GIS、城市规划方向的初学者和从业者特别是短期项目里需要快速出结果、又不想折腾语义分割数据集的场景。1. 项目概述为什么用 k-Means 做遥感影像语义分割1.1 项目背景与核心需求我接到这个任务时手头是一幅某区域的高分辨率多光谱遥感影像空间分辨率约 2 米波段包含蓝、绿、红、近红外四个波段范围覆盖大概几十平方公里。目标很明确把影像里的地表覆盖类型自动划分出来生成一份地类图斑矢量数据供后续规划分析使用。以前处理这种需求要么人工目视解译一个大范围区域就要熬好几个通宵要么用深度学习但需要先制作语义分割数据集逐像素标注大量样本还要准备 GPU 环境短时间根本完不成。k-Means 的优势恰好在这里它是无监督算法不需要任何标注样本实现极其简单sklearn 里几行就能调用计算开销小普通 CPU 就能处理几百万像素。对“快速出结果、精度要求并非极致”的业务场景k-Means 是最省事的突破口。当然它也有明显短板比如对噪声敏感、类别数需要人工指定、只考虑光谱距离而忽略空间上下文这些在后面我会详细展开。这套流程本质上回答了一个问题在没有深度学习和标注数据的条件下如何用最小成本完成一次可交付的遥感影像语义分割任务。1.2 无监督聚类与深度学习分割的取舍很多人一听到“语义分割”就默认要上 DeepLabV3、U-Net 这些深度学习模型。确实在标注数据充足、训练充分的前提下深度模型的精度上限远高于 k-Means。我在另一个项目里也用过 DeepLabV3 系列做建筑物提取效果确实好但前提是花了整整一周做标注还得有一台显存够用的 GPU 机器。而 k-Means 走的是另一条路用像素在光谱空间中的距离自然聚类同类别像素在光谱特征上天然接近。它把“每个像素属于哪个类别”这个问题转化成“像素点在特征空间里离哪个簇中心最近”的问题。对于地表覆盖类型差异明显的影像比如水体、植被、裸地之间的光谱响应差异非常大k-Means 能轻松分出来。但如果地类之间光谱非常接近比如干土和水泥地或者阴影和深色水体单靠 k-Means 就会混成一团。我的经验是k-Means 适合做预分类、快速制图和辅助标注深度学习适合在数据条件允许时追求更高精度两者不是对立关系而是一条工作流里的不同环节。尤其是制作语义分割数据集的时候先用 k-Means 聚类出初稿再在 GIS 软件里人工修正标注效率能提升一个量级。2. 数据准备与影像预处理2.1 遥感影像数据选择思路直接用原始影像跑聚类不是不行但效果通常很一般。首先我会确认影像的格式和波段情况如果是直接下载的多光谱 L2A 级产品一般已经做过大气校正可以直接用如果只是原始 DN 值影像最好先做一遍辐射定标消除传感器响应差异和大气散射的影响。其实对快速分割这种需求只要波段之间可比DN 值和反射率的差别对结果影响没有想象中那么大但有一个细节必须注意各个波段的量纲和直方图分布要相近否则值域大的波段会在欧氏距离中占据主导聚类结果基本被这一个波段带偏。另外一个容易被忽略的点是影像的坐标系和范围。不同来源的影像可能是不同投影、不同分辨率直接叠加比对会错位。我在项目里会先用 gdalinfo 查看影像元信息确认投影、分辨率、波段数然后把研究区统一裁剪到同一范围。如果只是单景影像这一步可能多余但如果是多时相或多源数据融合统一坐标系和分辨率就是必须做的功课。2.2 预处理流程辐射定标、大气校正与裁剪我这次使用的影像是经过正射校正的多光谱数据所以重点做两件事裁剪研究区和波段整理。裁剪用 GDAL 的 gdal_translate 命令最方便输入影像范围和输出范围对上就行。命令行大致是这样gdal_translate -projwin 118.20 31.10 118.50 30.80 input.tif output_clip.tif-projwin 后面的四个参数分别是左上角 X、左上角 Y、右下角 X、右下角 Y注意坐标系必须与影像投影一致否则截出来是乱的。如果你用的是 WGS84 经纬度坐标系就直接填经纬度范围内的四至坐标。裁剪完之后我会顺手用 gdalinfo 看一眼波段数和数据类型确认输出是 4 波段或更多方便后面构建特征。大气校正这一步如果用的是 Level-2 级产品一般可以省掉如果拿到的是原始 Level-1 数据又急着出结果我建议至少做一次波段归一化让每个波段都缩放到 0~1 区间。这其实比严格的大气校正对聚类结果的稳定贡献更大因为 k-Means 基于欧氏距离特征尺度统一是前提。如果再讲究一点可以用 6S 或 FLAASH 这类模型做辐射校正但耗时长、参数多对快速分割来说性价比不高。注意如果影像存在 NoData 值通常是 -9999 或 0一定要在聚类前剔除否则这些无效值会被当成一个特殊类别参与聚类结果会出现大片异常色块。2.3 特征工程光谱特征与纹理特征纯光谱波段可以作为基础特征但只用四个波段做聚类类别边缘往往非常破碎尤其在城市区域阴影和屋顶点混在一起。我一般会额外构建两个特征NDVI 和纹理特征。NDVI 的计算公式是 (NIR - Red) / (NIR Red)能有效突出植被信息把植被和水体、裸地的区分度拉大。纹理特征我常用灰度共生矩阵GLCM里的对比度窗口选 3×3 或 5×5太小噪声大太大边缘会糊实测 3×3 在 2 米分辨率影像上性价比最高。构建特征的做法是把每个像元对应波段的数值、NDVI、纹理值一起拼接成一个多维向量。假设原始影像 4 个波段加 NDVI 再加 1 个纹理值就变成 6 维特征。所有像素摊平成一个 N×6 的矩阵N 是所有有效像元数量然后丢给 k-Means 去聚类。这里的核心理念是数据表达的质量决定了算法效果的上限k-Means 再怎么调参也救不了糟糕的特征表达。好的特征能让类别在空间中自动分开差的特征只会让聚类结果变成一锅粥。3. k-Means 聚类原理与参数解析3.1 k-Means 算法核心逻辑k-Means 的基本思路特别直白随机挑选 K 个初始质心然后反复迭代——每个像素归入最近质心对应的簇更新质心为簇内所有像素的均值直到质心位置不再变化。这个过程可以类比成“班级里按考试成绩分组先随便选几个组长其他人加入离自己成绩最近的组选完再重新算各组平均水平直到稳定”。对于遥感影像来说像素的特征向量就是“成绩”特征空间中的距离就是组与组之间的相异度。sklearn 里调用 k-Means 非常简单from sklearn.cluster import KMeans kmeans KMeans(n_clusters5, random_state42, n_init10, max_iter300) labels kmeans.fit_predict(feature_matrix)fit_predict 之后labels 数组里每个值 0~4 就是该像素的类别编号。看似简单背后有三个参数值得解释n_init 表示用多少个不同的随机初始质心去跑最终取误差最小的那次结果这能避免随机初始化导致的局部最优max_iter 是单次迭代上限一般 300 次足够收敛random_state 固定随机种子确保每次运行结果可复现。我在项目里固定 random_state42方便结果对比和排查问题。3.2 K 值选定方法K 值是 k-Means 里最让人头疼的超参数选少了类别分不开选多了同类被切碎。我常用的方法是肘部法则对 K 2~10 分别计算簇内误差平方和画一条 K-SSE 曲线找出曲线像手肘一样的拐点这个拐点对应的 K 就是相对合理的类别数。但实测下来遥感影像的拐点往往不太明显尤其在复杂地类区域曲线平滑下降根本找不到明显的“肘”。这时候我会结合业务经验来定 K比如我知道研究区大致有植被、水体、建筑、裸地四类就先设 K4再根据聚类结果局部调整。还有一种思路是用轮廓系数评价聚类效果值越接近 1 说明簇内紧凑、簇间分离好。我一般在 K 已经大概确定后用轮廓系数验证一下比如 K4 和 K5 哪个轮廓系数更高就选哪个。注意轮廓系数计算量大像素超过几十万时建议先随机抽样一部分像素来算不然内存和时间都吃不消。我实际做的时候先对影像做了降采样把几千万像素抽样到 50 万左右算轮廓系数几秒钟就能出结果。还有一个技巧不要只看单个指标最好把分类后的影像在 GIS 里打开目视检查各类的空间分布是否符合地学常识。比如水体应该是连通的面状分布不应该星星点点散落在山体上如果出现这种情况说明 K 值或者特征设计有问题单纯靠数字指标是发现不了的。3.3 特征标准化与聚类前的最后准备如果直接输入 DN 值或者反射率特征一定要做标准化。我用 StandardScaler 对特征矩阵进行 z-score 标准化让每个特征的均值为 0、方差为 1。为什么要做因为近红外波段的值域通常比蓝波段大很多如果不标准化近红外在欧氏距离中的权重会被无限放大聚类结果几乎只取决于近红外的值其它波段等于白提了。可以类比成比较两个人的“综合实力”一个人年龄差折算成米一个人身高差折算成岁单位都不一样没法直接比。标准化的作用就是把所有指标拉到同一把尺子上。标准化之后特征矩阵就可以直接丢进 KMeans 了。但还有一步容易被忽略如果影像异常大直接用全量像素聚类会很慢甚至内存溢出。我的做法是先随机抽样一部分像素比如 10 万到 50 万个在样本上训练 k-Means 得到簇中心再用 kmeans.predict 对全量像素进行分类。这样既保留了聚类精度又把时间和内存开销控制住了。4. 核心实现Python 完整流程4.1 环境搭建与依赖库我的开发环境是基于 conda 管理的 Python 3.9核心库是 rasterio、numpy、scikit-learn、scikit-image。rasterio 负责读写 GeoTIFFnumpy 负责数组运算sklearn 提供 KMeansscikit-image 提供纹理特征计算。GDAL 功能虽然更全但安装相对麻烦rasterio 对大多数遥感影像读取、裁剪需求已经足够。安装命令直接一行conda install rasterio scikit-learn scikit-image -c conda-forge如果电脑没有 conda用 pip 装也一样。我建议用 conda因为 rasterio 在 Windows 上依赖比较多conda 会自动处理底层库可以少踩很多坑。另外建议在 Jupyter Notebook 里先做交互式探索把影像读取、特征构建逐步跑通再整理成独立脚本开发和调试效率会高很多。4.2 影像读取与特征矩阵构建核心的读取函数是这样import numpy as np import rasterio def read_image_tensor(path): with rasterio.open(path) as src: data src.read() meta src.meta.copy() return data, metadata 是 (波段数, 高度, 宽度) 的三维数组。接下来要把这个三维数组转成“像素 × 特征”的二维矩阵。同时要注意处理无效值比如遥感影像的 NoData 值通常为 -9999 或 0这些像元不能参与聚类否则会造成离谱的类别。做法是先构建一个有效像元遮罩提取有效像素的索引聚类完再把这些索引位置填成 NoData。构建完整特征矩阵时我封装了一个函数def build_feature_matrix(data, nodata-9999): bands, h, w data.shape valid_mask np.all(data ! nodata, axis0) pixels data[:, valid_mask].T features [pixels] # NDVI近红外与红光波段的归一化差值 nir pixels[:, 3] red pixels[:, 1] denom np.maximum(nir red, 1e-6) ndvi (nir - red) / denom features.append(ndvi[:, None]) # GLCM 纹理基于灰度图像的对比度 from skimage.feature import graycomatrix, graycoprops gray np.mean(data, axis0) gray ((gray - gray.min()) / (gray.max() - gray.min()) * 255).astype(np.uint8) glcm graycomatrix(gray, distances[1], angles[0], levels256, symmetricTrue, normedTrue) contrast graycoprops(glcm, propcontrast)[..., 0, 0] features.append(contrast[valid_mask, None]) X np.concatenate(features, axis1) return X, valid_mask这段代码里的 NDVI 计算已经加了防零操作因为像元值可能出现 NIR Red 0 的情况除以 0 会出现 inf。GLCM 计算我取了距离为 1、方向为 0 的共现统计更精细的做法是取四个方向的均值能削弱方向性纹理偏差但耗时也翻倍。对于快速分割单个方向已经够用。4.3 聚类执行与结果映射特征矩阵准备好后聚类本身非常快from sklearn.preprocessing import StandardScaler from sklearn.cluster import KMeans scaler StandardScaler() X_scaled scaler.fit_transform(X) kmeans KMeans(n_clusters4, random_state42, n_init10) labels kmeans.fit_predict(X_scaled)在约 50 万有效像素上这段代码在普通笔记本 CPU 上跑完大约需要几秒到十几秒。如果影像非常大比如上亿像素先用 random.sample 抽 10% 像素做聚类再把全部像素输入 kmeans.predict 来预测速度能提升不少精度损失很小。聚类结束后要把一维 labels 映射回二维影像平面。做法是先建一个全 NoData 的二维数组然后把 valid_mask 对应的位置填上聚类标签label_img np.full((h, w), nodata, dtypenp.uint8) label_img[valid_mask] labels到这一步你已经得到了一个语义分割的栅格结果。但注意这个结果还是像素级的直接看会有点“花”因为每个像素单独赋类缺少空间连贯性。这也是下一节要处理后处理的原因。4.4 后处理噪声去除与图斑矢量化初步聚类结果通常是“椒盐状”的单个像素类别跳变严重因为 k-Means 只考虑了光谱距离完全没有空间上下文。我惯用的处理是众数滤波也就是用一个小窗口统计中心像素邻域内出现最多的类别把中心像素替换成众数类别。3×3 窗口去噪一次图面会干净很多但如果去噪窗口开太大比如 7×7小地物会被抹掉边界也变圆滑反而损失准确性。我实测 3×3 处理一次的效果最好像孤立噪点能消掉一大半。矢量化这一步可以直接用 rasterio 导出标签栅格再到 ArcGIS Pro 里用“栅格转面”工具转成面要素。考虑到很多人最后要在 ArcGIS Pro 里做后续编辑我更喜欢先在 Python 里存成 GeoTIFF再在 ArcGIS Pro 里转换这样步骤更透明出了问题也好排查。导出代码with rasterio.open(seg_result.tif, w, **meta) as dst: dst.write(label_img[np.newaxis, :, :])meta 是第一步读取时存下来的元数据但记得修改 count1 和 dtypeuint8否则写多波段或者浮点型会报错。导出后可以在 GIS 软件里叠加原影像做质量检查。5. 在 ArcGIS Pro 中加载与出图5.1 加载影像与聚类结果把遥感影像和聚类结果 GeoTIFF 直接拖进 ArcGIS Pro 的目录窗口地图视图里马上就能看到。加载完成后建议先右键聚类结果图层打开图层属性在“源”里确认像素类型和 NoData 值是否正确不然影像会出现整片黑色或白色。接下来把聚类结果的渲染方式改成“唯一值”并按聚类类别数设置分类。这一步我会给不同类别手动指定颜色植被用绿色水体用蓝色建筑用灰色裸地用土黄色。这样一来成果图立刻变得可读。ArcGIS Pro 里还提供了其他非监督分类工具比如“ISO 聚类非监督分类”它本质上和 k-Means 一脉相承只是内置了更完善的流程。如果你不想写代码直接在工具箱里选波段、指定类别数就能跑出结果。但用 Python 控制流程的好处是透明、可复现、参数可调而且可以批量处理多景影像。我自己的习惯是正式项目用 Python 脚本临时看一眼用工具箱。5.2 创建地类图斑矢量图层与符号化聚类结果栅格在 GIS 软件里更适合做分析底图但很多项目要求提交矢量面图层。ArcGIS Pro 里的操作是搜索“栅格转面”工具输入聚类结果栅格字段选择 Value输出要素类就是每块同类区域合并成一个多边形。这里有一个关键细节ArcGIS 的“栅格转面”默认会合并相邻且值相同的区域所以输出的图斑面数量不会太多但是边界呈锯齿状符合栅格数据特征。如果觉得锯齿太严重可以后续使用“平滑面”工具做一次平滑但要小心平滑过头导致边界偏离实际地物。生成面要素后我习惯再新建一个“地类图斑”矢量图层用 ArcGIS Pro 的“创建要素”功能绘制补充或修正区域比如把 k-Means 错分的区域手动改过来。在目录中新建要素类时选择面类型和坐标系进入编辑状态后就能手工绘制或修改。但既然已经有聚类结果矢量化出的面我更推荐优先用矢量化结果只有局部修正时才需要手工画。毕竟手工画费时费力用聚类结果作为底图逐类检查修改效率高一个量级。6. 常见问题与排查技巧实录6.1 典型问题与解决方案速查表我在实际操作中整理了以下常见问题做成速查表方便大家对照排查问题现象可能原因解决办法聚类结果一片混沌类别无规律特征未标准化或 NoData 参与计算先做 StandardScaler并用掩膜剔除无效值水体被分成好几类K 值过大或阴影被单独成簇减少 K或将阴影区域与水体合并处理所有像元几乎聚成一类特征值域差距过大某一波段主导距离标准化后重新聚类结果噪点非常多k-Means 没考虑空间上下文3×3 众数滤波或叠加纹理特征ArcGIS Pro 里显示全黑NoData 设置不对或拉伸方式错误检查 NoData改用唯一值渲染内存溢出像素太多特征矩阵过大抽样聚类predict 全量6.2 独家避坑经验第一不要迷信 K 值越大越好。我见过很多新手把 K 设成 8、10结果建筑物屋顶被按光照角度切成好几个类后期合并反而更麻烦。K 值宁小勿大类别不够可以分开后在 GIS 里合并类别多了归并起来就非常痛苦。第二NoData 处理要放在聚类之前。很多初学者直接读数据丢给 KMeans结果 NoData 区域被单独聚成一类或者聚成几个奇怪的类然后要花很久去处理这些异常区域。在读取时就把 NoData 掩膜做掉整个流程会干净很多。第三标准化之后记得保留 scaler 对象后续如果要把新影像输入同样的模型做预测需要用同一个 scaler 做转换否则特征分布不一致预测结果会偏。这个坑我踩过一次当时换了一景影像直接 predict结果类别严重错位排查了半天才发现是忘了重新标准化。第四很多时候“语义分割”并不需要一步到位。先用 k-Means 结果作为底图在 ArcGIS Pro 里做局部修正其实是制作语义分割训练数据集的最高效路径。深度学习模型需要逐像素标注样本手工画太累拿 k-Means 初分类结果当预标注人工只需要修边界和错分区域标注效率能提升一个数量级。这是我从实际项目中得出的体会特别适合后续想做 DeepLabV3 等深度模型但苦于没有训练数据的情况。7. 个人实操体会与扩展思路跑完整个流程我的感受是 k-Means 在遥感语义分割里的定位不是替代深度学习而是“快速启动”和“数据预标注”。项目工期短、没有标注样本、只有 CPU 机器这三条只要中了两条k-Means 基本就是最优解。而且它得到的聚类结果并不粗糙配合众数滤波和人工小修完全可以达到业务交付的及格线。后续还可以扩展的方向不少。比如把 k-Means 的簇中心当作初始值喂给 GMM 高斯混合模型或者谱聚类继续精化能更好地处理光谱分布不规则的类别。也可以用简单线性迭代聚类SLIC先做超像素分割再用 k-Means 对超像素区域聚类这样结果的空间连续性好很多边缘也更平滑。我在另一个实验里试过超像素加 k-Means 的组合噪点明显减少边界质量接近简单 CNN 的效果。最后说个实在的建议如果你只是临时用一次直接在 ArcGIS Pro 自带的“ISO 聚类非监督分类”工具里也能完成类似效果但用 Python 控制流程的好处是透明、可复现、参数可调。把特征工程、聚类和后处理写成脚本存起来下次换一张影像把路径改一改就能直接用这才是工程上最值得做的事。我后来把这套流程封装成了一个脚本任何新影像进来十分钟就能出结果效率比手动操作高了不知道多少倍。希望这套流程也能帮你少踩一些坑快速把影像变成能用的分类成果。