遥感解译精度评定，混淆矩阵与Kappa系数
遥感解译结果的精度验证是遥感影像分类的一个重要过程。只有经过精度验证，我们才能知道解译结果是否可靠。
大部分遥感应用软件中带有遥感分类结果精度评定功能，如易康，ENVI中都有精度评定功能，能够自动计算混淆矩阵、Kappa系数等，但是ArcGIS目视解译结果往往无法使用自动的精度评定功能进行计算，那么这种情况下如何进行精度评定呢？
精度评定大致需要以下几个步骤：
选取验证点 提取分类属性到验证点 数据透视表制作混淆矩阵 精度评定参数计算
1 验证点的选取
验证点应随机均匀分布在解译区域内，建议精度评定的没种类别应至少有50个点，当区域很大或者分类种类较多时，每类最少数量应增加到75~100个，这个数量应根据具体情况有所调整。
验证点选取方法
在这里推荐三种验证点的选取方法：
野外采样
野外调查成本高，但是精度好，能准确的判断地物类型，采集所需的光谱、土地覆被信息。
谷歌地球
Google Earth提供的影像分辨率高，精度好，一般可以用于Landsat、Sentinel、MODIS等中低分辨率的遥感解译精度验证。但是由于谷歌地球也是遥感影像，因此需要一些先验知识，不能完全代替野外调查。
点击上方标题跳转详细介绍和网站
这是一个神奇的网站，整数经纬度的地物、土地覆被展示。这个计划的目标是访问这个世界每一个整数经纬度地点，在每个地方拍照，这些照片、到访这里的故事，会被展示到这个网站上。可以部分代替野外调查工作。
验证点制作
由于验证点可以野外采集获得，谷歌地球获得，整数经纬度网站获得，获取方式很多，具体方法不再一一介绍。
给大家看一个例子。
两个图层，yzd为验证点，点要素；class为解译结果，面要素。
从属性表中可以看到，验证点中存储真实地类的字段为DLMC地类名称字段。我们要想制作混淆矩阵，还需要给验证点中加上分类结果。
由于分类后的为矢量，面要素，所以我们需要使用空间挂接功能（Spatial Join）将分类结果属性赋到验证点中。
这样，验证点的准备工作就完成了，接下来就是令人头大的计算过程啦
2 混淆矩阵制作
基本概念
误差矩阵（Error Matrix）又称混淆矩阵（Confusion Matrix），是一个用于表示分为某一类别的像元个数与地面检验为该类别数的比较矩阵。通常，矩阵中的列代表验证数据，行代表由遥感数据分类得到的类别数据。有像元数和百分比表示两种。
EXCEL制作混淆矩阵
前面空间挂接的结果，可以直接使用EXCEL打开SHP中的dbf文件，制作数据透视表，生成混淆矩阵。
根据基本概念中的说明，“列”为验证数据，在这里是DLMC，“行”为解译结果，在这里对应为classname
3 精度评定参数计算
以下内容为网络内容整理，具体参考资料在参考文献中，如有问题欢迎留言指出
以下面3×3矩阵为例，介绍计算方法：
| Class | Water | Forest | Farm |
|---|---|---|---|
| Water | a | b | c |
| Forest | d | e | f |
| Farm | g | h | i |
对于上表来说
列为验证数据 行为解译结果 解译正确的部分为对角线，a,e,i
总体分类精度（Overall Accuracy, OA）
等于被正确分类的像元总和除以总像元数； 被正确分类的像元数目沿着混淆矩阵的对角线分布； 总像元数等于混淆矩阵的和。
Kappa系数（Kappa Coefficient）
概念说明
Kappa系数是一个用于一致性检验的指标 可以用于衡量分类的效果 对于分类问题，所谓一致性就是模型预测结果和实际分类结果是否一致 Kappa系数的计算是基于混淆矩阵的，取值[-1,1]，通常大于0 OA可以直接反映分类正确的比例，计算简单，但是由于分类过程中各个样本不平衡，OA高的时候部分类别可能很低，Kappa系数则能反映这一问题，有一类分类结果差，Kappa值就会明显降低
Kappa系数计算
它是通过把所有真实参考的像元总数（N）乘以混淆矩阵对角线（XKK）的和，再减去各类中真实参考像元数与该类中被分类像元总数之积之后，再除以像元总数的平方减去各类中真实参考像元总数与该类中被分类像元总数之积对所有类别求和的结果。
错分误差（Commission Error, CE）
指被分为用户感兴趣的类，而实际上属于另一类的像元，错分误差显示在混淆矩阵的行里面 以Water类为例：
漏分误差（Omission Error, OE）
指本属于地表真实分类，但没有被分类器分到相应类别中的像元数。漏分误差显示在混淆矩阵的列里 以Water类为例：
制图精度（Producer’s Accuracy, PA）
制图精度也称“生产者精度” 指分类器将整个影像的像元正确分为A类的像元数（对角线值）与A类真实参考总数（混淆矩阵中A类列的总和）的比率 对于Water类来说：
用户精度（User's Accuracy, UA）
正确分到A类的像元总数（对角线值）与分类器将整个影像的像元分为A类的像元总数（混淆矩阵中A类行的总和）比率 同样以Water类为例：
精度评定参数计算程序
上面可以用EXCEL来算，但是计算矩阵还是比较麻烦的，我还是喜欢R，一次劳动，重复使用。
废话少说，直接上代码
library(readxl)
AccData <- read_excel(path = "Test.xlsx", sheet = 1, col_names = T)
#EXCEL中存在空值，需要补零
AccData[is.na(AccData)] <- 0
AccMatrix <- data.matrix(AccData[,2:ncol(AccData)])
#提取矩阵主对角线
MatDiag <- diag(AccMatrix)
#提取总数
TotalNum <- sum(AccMatrix)
#总体分类精度
OA <- sum(MatDiag)/TotalNum
#Kappa系数
#函数包求Kappa
library(vcd)
K <- Kappa(AccMatrix)
#手动计算Kappa
colFreqs <- colSums(AccMatrix)/TotalNum
rowFreqs <- rowSums(AccMatrix)/TotalNum
p0 <- sum(MatDiag)/TotalNum
pe <- crossprod(colFreqs, rowFreqs)[1]
k2 <- (p0-pe)/(1-pe)
#制图精度
for (i in 1:nrow(AccMatrix)) {
PA <- AccMatrix[i,i]/sum(AccMatrix[,i])
print(paste(AccData[i,1], "制图精度为", PA*100, "%"))
}
#用户精度
for (i in 1:nrow(AccMatrix)) {
UA <- AccMatrix[i,i]/sum(AccMatrix[i,])
print(paste(AccData[i,1], "用户精度为", UA*100, "%"))
}
#计算结果输出
print(paste("总体分类精度为", OA*100, "%"))
print(paste("Kappa系数为", K[["Unweighted"]][["value"]]*100, "%"))
参考文献
赵英时. 遥感应用分析原理与方法[M]. 科学出版社, 2003. https://blog.csdn.net/momozhongxia/article/details/72599454 https://www.rdocumentation.org/packages/vcd/versions/1.4-7/topics/Kappa https://rdrr.io/cran/vcd/src/R/Kappa.R https://zhuanlan.zhihu.com/p/67844308
点击阅读原文查看培训班视频课程与讲解