遥感分类精度评价方法--混淆矩阵和kappa系数

遥感影像分类之后需要进行分类精度评价，精度评价方法中最常见的就是混淆矩阵和kappa系数。现把指标列举如下:
混淆矩阵（confusion matrix）
误差矩阵（error matrix）又称混淆矩阵（confusion matrix），是一个用于表示分为某一类别的像元个数与地面检验为该类别数的比较阵列。通常，阵列中的列代表参考数据，行代表由遥感数据分类得到的类别数据。有像元数和百分比表示两种。
总体分类精度（Overall Accuracy）
等于被正确分类的像元总和除以总像元数。被正确分类的像元数目沿着混淆矩阵的对角线分布，总像元数等于所有真实参考源的像元总数
Kappa系数（Kappa Coefficient）
它是通过把所有真实参考的像元总数（N）乘以混淆矩阵对角线（XKK）的和，再减去各类中真实参考像元数与该类中被分类像元总数之积之后，再除以像元总数的平方减去各类中真实参考像元总数与该类中被分类像元总数之积对所有类别求和的结果。
错分误差(Commission Error)
指被分为用户感兴趣的类，而实际上属于另一类的像元，错分误差显示在混淆矩阵的行里面。
漏分误差（Omission Error）
指本属于地表真实分类，但没有被分类器分到相应类别中的像元数。漏分误差显示在混淆矩阵的列里。
制图精度（Producer’s Accuracy）
制图精度或生产者精度是指分类器将整个影像的像元正确分为A类的像元数（对角线值）与A类真实参考总数（混淆矩阵中A类列的总和）的比率。
用户精度（User’s Accuracy）
是指正确分到A类的像元总数（对角线值）与分类器将整个影像的像元分为A类的像元总数（混淆矩阵中A类行的总和）比率。
现举例说明：
研究区中有7种类别地物，基于面向对象和SVM结合的方法得到分类结果。
SVM采用LibSVM工具包实现，LibSVM可以得到.predict文件。在ArcGIS中把.predict结果加入到eCognition的分割结果.shp中，得到分类结果可视化图；然后把人工矢量化参考影像，根据共有字段 join到分割结果.shp图层中，至此，参考数据和预测数据准备完毕。
利用matlab实现上述指标，直接上代码。

clc;
clear all;
close all;
%% read reference data
t_fn = 'D:\...\area1intersect200.xls';
t_ft = xlsread(t_fn);
t_label = t_ft(:,1);

%% read predicted data
pred_fn = 'D:\...\predict200and0.1.xls';
pred_ft = xlsread(pred_fn);
pred_label = pred_ft(:,1);
%pred_label = dlmread(pred_fn);
%% 
%找到lael1
ind1 = find(t_label == 1);
label1 = t_label(ind1);
%label_1=label(label==1);
%ft1 = ft(ind1,:);

%找到lael2
ind2 = find(t_label == 2);
label2 = t_label(ind2);
%ft2 = ft(ind2,:);

%找到lael3
ind3 = find(t_label == 3);
label3 = t_label(ind3);
%ft3 = ft(ind3,:);

%找到lael4
ind4 = find(t_label == 4);
label4 = t_label(ind4);
%ft4 = ft(ind4,:);

%找到lael5
ind5 = find(t_label == 5);
label5 = t_label(ind5);

%找到lael6
ind6 = find(t_label == 6);
label6 = t_label(ind6);

%找到lael7
ind7 = find(t_label == 7);
label7 = t_label(ind7);



%% calculate confusion matrix
cls = unique(t_label);
num_ref = zeros(length(cls),1);
cf_mt = zeros(length(cls));
for i=1:length(cls)
    ind = find(t_label == cls(i));
    num_ref(i) = length(ind);
    for j = 1:length(cls)
        ind1 = find(t_label == cls(i) & pred_label == cls(j));
        cf_mt(i,j) = length(ind1);
    end
end
%% 计算PA
PA1 = cf_mt(1,1)/sum(cf_mt(1,:));
PA2 = cf_mt(2,2)/sum(cf_mt(2,:));
PA3 = cf_mt(3,3)/sum(cf_mt(3,:));
PA4 = cf_mt(4,4)/sum(cf_mt(4,:));
PA5 = cf_mt(5,5)/sum(cf_mt(5,:));
PA6 = cf_mt(6,6)/sum(cf_mt(6,:));
PA7 = cf_mt(7,7)/sum(cf_mt(7,:));
%% 计算UA
UA1 = cf_mt(1,1)/sum(cf_mt(:,1));
UA2 = cf_mt(2,2)/sum(cf_mt(:,2));
UA3 = cf_mt(3,3)/sum(cf_mt(:,3));
UA4 = cf_mt(4,4)/sum(cf_mt(:,4));
UA5 = cf_mt(5,5)/sum(cf_mt(:,5));
UA6 = cf_mt(6,6)/sum(cf_mt(:,6));
UA7 = cf_mt(7,7)/sum(cf_mt(:,7));

cf_mt(isnan(cf_mt))=0;
[r,c] = size(cf_mt);

RT_A = sum(cf_mt,2); % row total
CT_A = sum(cf_mt,1); % column total

N_A = sum(sum(cf_mt));

PO_A = trace(cf_mt)/N_A; % PO: percentage overall accuracy of A

PC_A = 1/N_A^2 * sum(RT_A.*CT_A'); % PC: proportion of pixels that were classified due to chance agreement

%% Calculate kappa statistics
Kappa = (PO_A-PC_A)/(1-PC_A);