遥感影像分类这件事,很多人把精力全砸在选模型、调参数、跑训练上,等分类图一出来就以为大功告成。真正做过完整项目的人都知道,分类器输出的那张栅格图只是半成品,后面还有一大段路要走:碎斑要清理、类别要合并、精度要量化、误差要归因。我带过几个遥感解译的活儿,最后交付时甲方最关心的往往不是"你用了什么网络",而是"这块地的分类精度到底有多少、错在哪、能不能改"。这篇就围绕分类后处理与误差分析这条链路,把每一步为什么做、怎么做、容易踩什么坑讲透,适合已经跑通分类、准备做精度评估和成果整理的朋友参考。
1. 分类后处理到底在解决什么问题
1.1 分类器输出为什么不能直接用
监督分类或者深度学习分割出来的结果,本质上是逐像元的类别标签。模型只看单个像元(或局部感受野)的光谱和纹理特征,它不知道"一片耕地中间不应该孤零零冒出一个水体像元"这种空间上下文常识。于是结果图里经常出现大量椒盐噪声——专业上叫"麻点"或"碎斑",表现为零星散落的、与周围类别不一致的小图斑。
这些碎斑的来源主要有三类。第一是光谱混淆,比如阴影处的水体和深色沥青在部分波段反射率接近,模型在边界处来回横跳。第二是训练样本不均衡,某些类别样本少,模型对它的判别边界不稳,容易在边缘区域误判。第三是影像本身的噪声,比如传感器条带、云影残留、配准误差导致的混合像元。无论哪种来源,直接拿这张图去做面积统计或者变化检测,误差都会被放大。
我一般会先做一次目视检查,把结果图和原始影像叠加,重点看几类典型区域:地块内部是否干净、地类边界是否贴合、线状地物(道路、沟渠)是否连续。这一步不用工具,纯靠眼睛,但能快速判断后处理的工作量有多大。
1.2 后处理的目标不是"好看"而是"可信"
很多人误以为后处理就是把图修得漂亮,其实核心目标是让分类结果在空间上和语义上都更接近真实地表覆盖。具体拆成三个可量化的目标:
- 消除孤立碎斑:小于某个面积阈值的图斑,大概率是噪声而非真实地物,需要合并到邻域主导类别。
- 保持边界合理性:地类边界应该和影像上的实际边界(田埂、道路、水岸线)大致吻合,不能因为平滑操作把边界糊掉。
- 维持类别面积稳定:后处理前后各类别总面积不应发生剧烈变化,如果某个类别面积骤降,说明参数设置有问题。
这三个目标之间是有张力的。过度平滑会吃掉真实的小图斑(比如一块小菜地),平滑不足又清不干净噪声。所以后处理的参数需要结合影像分辨率来定,不能一套参数走天下。
1.3 后处理在整体流程中的位置
一个完整的遥感分类项目,流程大致是:数据预处理(辐射定标、大气校正、几何配准)→ 样本采集 → 特征工程 → 分类器训练与推理 →分类后处理→精度评价→ 成果输出与制图。后处理紧跟在推理之后,精度评价之前,这个顺序不能颠倒。
原因很简单:精度评价要评价的是"最终交付成果"的精度,而不是中间产物。如果你先评价再后处理,那评价结果和交付成果对不上,甲方验收时会出问题。我见过有人先算完混淆矩阵,后处理又改了一堆像元,最后报告里的精度和实际图对不上,返工重算,白白浪费两天。
2. 碎斑清理与类别合并的实操手法
2.1 众数滤波:最常用的碎斑清理手段
众数滤波(Majority Filter)的思路很直白:对每个像元,看它周围一个窗口(比如3×3或5×5)内哪个类别出现次数最多,就把中心像元改成那个类别。这相当于用邻域投票的方式压制孤立噪声。
在主流遥感软件里都有现成工具。以常见的栅格处理库为例,用Python实现一个简单的众数滤波逻辑大致是这样:
import numpy as np from scipy import ndimage def majority_filter(label_map, window_size=3): """ label_map: 二维整型数组,每个值代表类别 window_size: 奇数窗口,如3、5、7 """ def mode_func(values): vals, counts = np.unique(values, return_counts=True) return vals[np.argmax(counts)] result = ndimage.generic_filter( label_map, mode_func, size=window_size, mode='nearest' ) return result.astype(label_map.dtype)窗口大小的选择是关键。3×3窗口比较温和,只清理真正的孤立点;5×5会明显平滑,适合噪声多的结果;7×7以上就要谨慎了,容易把真实的小地物抹掉。我的经验是:先看影像分辨率。如果分辨率是10米级(比如Sentinel-2),一个3×3窗口覆盖30米×30米,相当于一个地块的最小单元,用3×3比较安全。如果是亚米级高分影像,3×3只覆盖几米,可以适当放大到5×5。
注意:众数滤波会改变像元类别,但不会改变类别数量。如果你的分类结果里有大量碎斑,滤波后面积统计会变化,务必在精度评价前完成,并记录滤波参数。
2.2 图斑合并:基于最小面积阈值的处理
众数滤波解决的是像元级噪声,但还有一种情况是成片的小图斑——比如一片林地中间有个几十个像元的小裸地,滤波窗口盖不住。这时候要用图斑级的处理:先做连通域分析,把相邻同类像元聚成图斑,然后统计每个图斑的面积,小于阈值的合并到相邻最大图斑。
from skimage.measure import label, regionprops import numpy as np def merge_small_regions(label_map, min_pixels=20): result = label_map.copy() for cls in np.unique(label_map): if cls == 0: # 假设0是背景 continue mask = (label_map == cls) labeled = label(mask, connectivity=2) for region in regionprops(labeled): if region.area < min_pixels: # 找到该图斑的边界邻域,取出现最多的其他类别 coords = region.coords neighbor_vals = [] for y, x in coords: for dy in [-1, 0, 1]: for dx in [-1, 0, 1]: ny, nx = y+dy, x+dx if 0 <= ny < label_map.shape[0] and 0 <= nx < label_map.shape[1]: if label_map[ny, nx] != cls: neighbor_vals.append(label_map[ny, nx]) if neighbor_vals: vals, counts = np.unique(neighbor_vals, return_counts=True) result[coords[:, 0], coords[:, 1]] = vals[np.argmax(counts)] return result最小面积阈值怎么定?这取决于你的制图比例尺和最小成图单元。比如做1:10000的图,最小图斑对应实地面积,换算到像元数就是阈值。一般经验值:10米分辨率影像,阈值设在10到50个像元之间;亚米级影像可以设到100以上。阈值太小清不干净,太大又会误伤真实小地物,需要反复试。
2.3 类别合并:从细分到应用导向
有时候分类体系本身太细,比如把"水田""旱地""菜地"分成三类,但应用方只关心"耕地"这一类。这时候需要做类别合并(Class Merging),把多个细类映射到一个大类。
合并的原则是语义一致性:合并后的类别在应用场景下应该有统一的含义。比如做城市扩张监测,把"高密度建筑""低密度建筑""工业用地"合并成"建设用地"是合理的;但如果把"林地"和"草地"合并成"植被",在做生态评估时可能就丢失了关键信息。
合并操作本身很简单,就是查表替换:
merge_dict = {1: 1, 2: 1, 3: 1, 4: 2, 5: 2, 6: 3} # 1,2,3合并为1;4,5合并为2;6保持为3 merged = np.vectorize(merge_dict.get)(label_map)但要注意:合并会改变混淆矩阵的结构。如果你先算了细类的混淆矩阵,再合并类别,需要把混淆矩阵的行列按映射关系相加,不能直接丢弃。这个后面讲误差分析时会细说。
2.4 边界保持:别把好边界也平滑掉
后处理最容易犯的错是"一刀切"平滑,把真实地类边界也模糊了。解决办法是引入边界约束:只在非边界区域做滤波,边界像元保持原样。
判断边界的方法是对分类图做梯度运算,梯度大的地方就是类别交界处。可以用简单的形态学操作:
from scipy.ndimage import binary_dilation def boundary_aware_filter(label_map, window_size=3): # 检测边界:某像元邻域内存在不同类别 boundary = np.zeros_like(label_map, dtype=bool) for dy in [-1, 0, 1]: for dx in [-1, 0, 1]: shifted = np.roll(np.roll(label_map, dy, axis=0), dx, axis=1) boundary |= (shifted != label_map) # 边界膨胀一圈,保护边界邻域 boundary = binary_dilation(boundary, iterations=1) filtered = majority_filter(label_map, window_size) result = label_map.copy() result[~boundary] = filtered[~boundary] return result这样处理下来,地块内部干净了,边界还保持锐利。代价是边界附近的碎斑清不掉,但那些碎斑往往就是混合像元造成的,强行清理反而可能引入新误差,保留它们更诚实。
3. 混淆矩阵:精度评价的核心工具
3.1 混淆矩阵的构造逻辑
混淆矩阵(Confusion Matrix)是精度评价的基石,它把"真实类别"和"预测类别"交叉列表,每一行代表真实类别,每一列代表预测类别(或者反过来,不同软件约定不同,报告里要写清楚)。矩阵对角线上的数值是正确分类的样本数,非对角线是错分样本数。
构造混淆矩阵的前提是有验证样本,而且验证样本必须和训练样本独立。这是最容易被忽视的一点。我见过有人图省事,直接拿训练样本当验证样本,算出来精度95%以上,结果实际应用一塌糊涂。这叫"自证清白",没有意义。验证样本要单独采集,最好用不同的来源(比如高分辨率影像目视解译、外业实测点),并且要覆盖各个类别和空间分布。
验证样本的数量也有讲究。统计学上,每个类别至少需要50个样本才能得到较稳定的精度估计,类别多的话建议每类100以上。样本太少,精度置信区间会很宽,结论不可靠。
3.2 从混淆矩阵能读出哪些指标
光看总体精度(Overall Accuracy)是不够的,它会被大类主导。比如一个区域90%是耕地,模型把所有像元都判成耕地,总体精度也有90%,但其他类别全错。所以要结合多个指标看:
| 指标 | 含义 | 关注点 |
|---|---|---|
| 总体精度 OA | 对角线和 / 总样本数 | 整体正确率,受大类影响大 |
| 生产者精度 PA | 某类正确数 / 该类真实总数 | 漏分情况,PA低说明漏分多 |
| 用户精度 UA | 某类正确数 / 该类预测总数 | 错分情况,UA低说明误判多 |
| Kappa系数 | 排除随机一致后的 agreement | 整体一致性,消除偶然因素 |
生产者精度和用户精度是一对容易混淆的概念。用生活化的例子:假设有100个真实的水体像元,模型只识别出70个,那生产者精度是70%,说明有30%的水体被漏掉了(漏分)。反过来,模型预测了80个水体像元,其中70个是真的,那用户精度是87.5%,说明模型说是水体的地方有12.5%是错的(错分)。做应用的人更关心用户精度——因为用户拿到分类图,看到"水体"标签,他希望这个标签是对的。
3.3 Kappa系数的计算与解读
Kappa系数衡量的是"分类结果比随机猜测好多少",公式是:
Kappa = (Po - Pe) / (1 - Pe)其中Po是观测一致率(就是总体精度),Pe是随机一致率(由各类别边际概率算出的期望一致率)。Kappa的取值和解读一般参考这个标准:
- 0.81 ~ 1.00:几乎完全一致
- 0.61 ~ 0.80:高度一致
- 0.41 ~ 0.60:中等一致
- 0.21 ~ 0.40:一般一致
- 0.00 ~ 0.20:微弱一致
但要注意,Kappa在类别极度不均衡时会有偏差,近年来有学者建议用其他指标(如F1、宏平均精度)补充。我的做法是报告里同时给出OA、Kappa和各类的PA/UA,让读者自己判断,不单独依赖某一个数。
3.4 用Python算混淆矩阵的完整代码
网上搜"python多分类混淆矩阵代码"能出来一堆,但很多没讲清楚标签对齐的问题。下面给一个完整可用的版本,重点处理了标签映射和指标输出:
import numpy as np import pandas as pd from sklearn.metrics import confusion_matrix, cohen_kappa_score def evaluate_classification(y_true, y_pred, class_names=None): """ y_true, y_pred: 一维数组,元素为类别标签(整数) class_names: 类别名称列表,顺序与标签对应 """ labels = sorted(np.unique(np.concatenate([y_true, y_pred]))) cm = confusion_matrix(y_true, y_pred, labels=labels) if class_names is None: class_names = [str(l) for l in labels] # 总体精度 oa = np.trace(cm) / np.sum(cm) # Kappa kappa = cohen_kappa_score(y_true, y_pred) # 逐类PA、UA rows = [] for i, name in enumerate(class_names): tp = cm[i, i] pa = tp / np.sum(cm[i, :]) if np.sum(cm[i, :]) > 0 else 0 ua = tp / np.sum(cm[:, i]) if np.sum(cm[:, i]) > 0 else 0 f1 = 2 * pa * ua / (pa + ua) if (pa + ua) > 0 else 0 rows.append({'类别': name, 'PA': round(pa, 4), 'UA': round(ua, 4), 'F1': round(f1, 4)}) df = pd.DataFrame(rows) print(f"总体精度 OA: {oa:.4f}") print(f"Kappa系数: {kappa:.4f}") print(df.to_string(index=False)) # 输出混淆矩阵 cm_df = pd.DataFrame(cm, index=class_names, columns=class_names) print("\n混淆矩阵:") print(cm_df.to_string()) return cm, oa, kappa, df这段代码有几个细节值得说。第一,labels参数显式指定了类别顺序,避免sklearn自动排序导致标签错位。第二,PA和UA的分母都做了零保护,防止某类样本为空时报错。第三,同时输出F1,方便对比。实际用的时候,把验证样本的真实标签和预测标签传进去就行。
提示:如果你的验证样本是点矢量,需要先做"点提取栅格值"操作,把分类图上的类别值提取到点上,再和真实标签对比。这一步在GIS软件里叫"Extract Values to Points"或类似名称。
4. 误差来源分析与精度提升思路
4.1 从混淆矩阵反推误差类型
混淆矩阵不只是给个精度数字,它还能告诉你误差的结构。看非对角线元素,哪两个类别之间错分最多,就说明这两个类别在特征空间里重叠严重。
比如混淆矩阵显示"耕地"和"草地"之间错分很多,可能的原因有:两者光谱曲线接近(尤其在生长季)、训练样本边界不清、纹理特征区分度不够。再比如"建筑"和"裸地"错分,往往是建筑阴影和裸地光谱相似导致的。
我习惯把混淆矩阵画成热力图,一眼就能看出误差集中在哪。虽然这里不画图,但思路是:行归一化后,对角线越亮越好,非对角线的亮块就是问题区域。
4.2 样本层面的误差排查
误差分析要往下钻到样本层面。把错分的验证样本挑出来,回到原始影像上看它们的位置和影像特征。常见的情况有:
- 混合像元:验证点落在两个地类的交界处,本身归属就模糊,这种错分不算模型的锅。
- 样本标注错误:验证样本本身标错了,或者标注时参考的影像和分类用的影像时相不一致。
- 类别定义模糊:比如"灌木"和"幼林"的界定标准不统一,不同人标注结果不同。
我一般会做一张错分样本分布图,看看错分是随机分布还是集中在某些区域。如果集中在某个区域,可能是该区域影像质量差(云影、阴影)或者训练样本没覆盖到;如果随机分布,那更可能是类别本身的可分性问题。
4.3 特征与分类器的改进方向
如果误差分析发现某两类持续混淆,可以从特征入手。比如加入纹理特征(GLCM)、植被指数(NDVI)、水体指数(NDWI)等,增强类别间的可分性。对于建筑和裸地的混淆,加入近红外波段和建筑指数往往有效。
分类器层面,随机森林和SVM对样本不均衡比较敏感,可以通过调整类别权重或者过采样来改善。深度学习模型可以尝试加入空间上下文信息(比如用U-Net代替逐像元分类),但要注意过拟合风险。
不过要提醒一句:精度不是越高越好,而是要和成本匹配。如果当前精度已经满足应用需求,没必要为了提升1个百分点去堆复杂的模型。我见过有人为了把OA从88%提到90%,多花了一周时间调参,结果应用方根本用不到这个精度。
4.4 精度报告怎么写才专业
最后交付的精度报告,不能只丢一个数字。一份合格的报告应该包含:
- 验证样本说明:来源、数量、采集方式、是否独立于训练样本。
- 混淆矩阵:完整的矩阵表格,标注行列含义。
- 精度指标表:OA、Kappa、各类PA/UA/F1。
- 误差分析:主要错分类型、可能原因、对应用的影响。
- 后处理说明:做了哪些后处理、参数是什么、对精度的影响。
这样甲方拿到报告,不仅知道精度多少,还知道误差在哪、能不能接受、后续怎么改进。这才是专业交付该有的样子。
5. 几个容易踩的坑和我的处理习惯
5.1 混淆矩阵"总合不唯一"是怎么回事
网上有人搜"yolo混淆矩阵总合不唯一",这个问题在遥感分类里同样存在。现象是:不同次运行、或者用不同工具算出来的混淆矩阵,总数对不上。原因通常有三个。
第一是标签顺序不一致。sklearn的confusion_matrix默认按标签值排序,如果你的标签是字符串或者不连续整数,不同工具的处理方式可能不同,导致矩阵行列错位,看起来总数就乱了。解决办法是显式传入labels参数。
第二是验证样本有重复或缺失。如果验证点矢量里有重复点,或者某些点的栅格提取值为NoData被丢弃,参与计算的样本数就会变化。我一般会在计算前先做一次数据清洗,去重、剔除无效值,并记录最终样本数。
第三是多分类与二分类的转换差异。有些工具在算多分类混淆矩阵时,会做one-vs-rest展开,导致总数膨胀。要确认你用的是标准的多分类混淆矩阵,而不是展开后的版本。
5.2 后处理顺序不能乱
正确的顺序是:分类推理 → 众数滤波 → 小图斑合并 → 类别合并 → 精度评价。如果先做类别合并再做滤波,滤波窗口内的类别数变少,投票结果会偏向大类,可能过度平滑。如果先做精度评价再做后处理,那评价的就是中间产物,没有意义。
还有一个细节:滤波和合并的参数要记录。我习惯在项目文档里写清楚每一步的参数,比如"众数滤波3×3,最小图斑阈值20像元,类别合并映射表见附件"。这样别人复现或者自己回头查都有据可依。
5.3 验证样本的独立性怎么保证
这是精度评价里最容易被糊弄的环节。严格来说,验证样本和训练样本应该来自不同的采样单元,最好在空间上也不重叠。实际操作中,如果样本量有限,至少要做到:训练样本和验证样本不共用同一个点,且验证样本的标注者最好和训练样本标注者不同(避免个人偏差)。
如果实在做不到完全独立,那在报告里要如实说明,并指出精度可能偏高。诚实比好看的数字更重要,这是我做了这么多项目最深的体会。
5.4 面积统计与精度脱节的问题
最后成果往往要给面积统计表,比如"耕地多少公顷、水体多少公顷"。这里有个陷阱:分类图上的像元数乘以像元面积,得到的是"分类面积",不是"真实面积"。如果某类用户精度只有80%,那这个面积统计就有20%的虚高或虚低。
严谨的做法是做面积校正,用混淆矩阵的错分比例反推真实面积。虽然校正后面积也不一定完全准确,但比直接用分类面积靠谱。如果项目要求不高,至少在报告里注明"面积为分类面积,未做精度校正"。
6. 写在最后的一点个人体会
遥感分类的后处理和误差分析,技术含量不比模型训练低,但往往被轻视。我见过太多项目,模型跑得漂亮,最后交付时因为精度报告不规范、后处理没做,被要求返工。其实这些工作并不难,关键是流程要完整、参数要记录、结论要诚实。
另外,工具的选择不用太纠结。QGIS、ArcGIS、ENVI、Python都能做这些事,核心是理解每一步在干什么。我个人的习惯是能用脚本就用脚本,因为可复现、可追溯,参数改起来也方便。但如果是单次项目、数据量不大,用现成软件点几下也完全没问题,没必要为了"显得专业"去写一堆代码。
最后分享一个小技巧:做精度评价时,先把混淆矩阵算出来,别急着看总体精度,先看非对角线。非对角线上的大数值才是真正要解决的问题,总体精度只是结果。把误差结构搞清楚了,改进方向自然就出来了。