高光谱影像分类

有监督高光谱遥感数据分类
clc,clear,close all;
%% 数据说明
% Pavia University 数据是由德国的机载反射光学光谱成像仪（Reflective Optics Spectrographic Imaging System，ROSIS-03）
% 在 2003 年对意大利的帕维亚城所成的像的一部分高光谱数据。该光谱成像仪对 0.43-0.86μm 波长范围内的 115 个波段连续成像，
% 所成图像的空间分辨率为 1.3m。其中 12 个波段由于受噪声影响被剔除，因此一般使用的是剩下 103 个光谱波段所成的图像。
% 该数据的尺寸为 610×340，因此共包含2207400 个像素，但是其中包含大量的背景像素，包含地物的像素总共只有 42776 个，
% 这些像素中共包含 9 类地物，包括树、沥青道路（Asphalt）、砖块（Bricks）、牧场（Meadows）等。
% PaviaU是大学城地区数据，630340，一共有103个波段
% PaviaU_gr是大学城地区数据标签，一共有9个类别，加上空白区域是10个，空白区域直接填充
% PaviaU是大学城地区数据，630340，一共有103个波段
% PaviaCenter是中心城区数据，1094*492，一共102个波段
% PaviaCenter_gt是中心城区数据标签，一共有9个类别
%% 实验内容
% 利用原始高维的光谱特征作为影像特征，进行影像分类；
% 利用PCA降维方法对原始高维光谱特征进行降维，然后利用降维后的光谱特征进行影像分类；
% 在实验二的基础上，加入局部窗口的均值和方差特征，进行影像分类；
% 注：分类器选用支持向量机(SVM)
% %% 加载数据
% PaviaU=load(‘PaviaU.mat’);
% PaviaU_gr=importdata(‘PaviaU_gt.mat’);
% PaviaCenter=importdata(‘PaviaCenter.mat’);
% PaviaCenter_gt=importdata(‘PaviaCenter_gt.mat’);
%% 绘原始图
% PaviaU图
% 波段的像素范围不是0-255，很多都超出这个范围本身
PaviaU_RGB=zeros(610,340,3);
PaviaU_R=PaviaU.paviaU(:,:,57);
PaviaU_G=PaviaU.paviaU(:,:,34);
PaviaU_B=PaviaU.paviaU(:,:,3);
PaviaU_R_255 =Normalize(PaviaU_R); %归一化并取整
PaviaU_G_255 =Normalize(PaviaU_G);
PaviaU_B_255 =Normalize(PaviaU_B);
PaviaU_RGB_255(:,:,1)=PaviaU_R_255;
PaviaU_RGB_255(:,:,2)=PaviaU_G_255;
PaviaU_RGB_255(:,:,3)=PaviaU_B_255;
PaviaU_RGB_255=uint8(PaviaU_RGB_255);

% PaviaCenter图
PaviaCenter_RGB=zeros(1096,492,3);
PaviaCenter_R=PaviaCenter(:,:,57);
PaviaCenter_G=PaviaCenter(:,:,34);
PaviaCenter_B=PaviaCenter(:,:,3);
PaviaCenter_R_255 =Normalize(PaviaCenter_R); %归一化并取整
PaviaCenter_G_255 =Normalize(PaviaCenter_G);
PaviaCenter_B_255 =Normalize(PaviaCenter_B);
PaviaCenter_RGB_255(:,:,1)=PaviaCenter_R_255;
PaviaCenter_RGB_255(:,:,2)=PaviaCenter_G_255;
PaviaCenter_RGB_255(:,:,3)=PaviaCenter_B_255;
PaviaCenter_RGB_255=uint8(PaviaCenter_RGB_255);

%% 数据表及标签制作
% PriviaU原始数据
PriviaU_table_all=zeros(610340,103);
img_row=610;
img_col=340;
band_num=103;
m=1;
% 用reshape应该快一些
for i=1:img_row
for j=1:img_col
PriviaU_table_all(m,1:band_num)=PaviaU.paviaU(i,j,1:band_num);%一行一行索引
m=m+1;
end
end
% PriviaU原始数据标签
PriviaU_tag_all=zeros(img_rowimg_col,1);
m=1;
for i=1:img_row
for j=1:img_col
PriviaU_tag_all(m,1)=PaviaU_gr(i,j);
m=m+1;
end
end
% 删除背景，减少数据量
data_index=find(PriviaU_tag_all);%记录非0元素位置，便于后续恢复图像
PriviaU_table=PriviaU_table_all;
PriviaU_tag=PriviaU_tag_all;
PriviaU_table(PriviaU_tag_all0,:)=[];
PriviaU_tag(PriviaU_tag_all0,:)=[];
%% 分类样本随机选取
feature_num=size(PriviaU_table,2);%波段特征个数
% 十折交叉方式选取
times=20;
PriviaU_train=zeros(1,feature_num+2);
PriviaU_test =zeros(1,feature_num+2);
train_index=zeros(1,1);
test_index=zeros(1,1);
for k=1:9
% 分类选取
PriviaU_Index=logical(PriviaU_tag(:,1)k);% Logical矩阵可记录类别位置
PriviaU_class=PriviaU_table(PriviaU_Index1,:);
class_index=find(PriviaU_Index==1);
PriviaU_class(:,feature_num+1)=k;% 记录类别
PriviaU_class(:,feature_num+2)=class_index;% 记录位置
indices = crossvalind(‘Kfold’, size(PriviaU_class,1), times);
for i = 1:5 % 3/20，即15%作为训练集
train = logical(indices == i);% 记录选择数据位置
PriviaU_train_5 = PriviaU_class(train, 😃;
PriviaU_train=[PriviaU_train;PriviaU_train_5];% 串联数据矩阵
end
for j = 6:20 % 17/20，即85%作为测试集
test = logical(indices == j);
PriviaU_test_5 = PriviaU_class(test, 😃;
PriviaU_test=[PriviaU_test;PriviaU_test_5];
end
end
PriviaU_train(1,:)=[];
PriviaU_test(1,:)=[];
% 原始样本训练
PriviaU_train_lable=PriviaU_train(:,feature_num+1);
PriviaU_train_data=PriviaU_train(:,1:feature_num);
PriviaU_test_lable=PriviaU_test(:,feature_num+1);
PriviaU_test_data=PriviaU_test(:,1:feature_num);
cmd1 =[’-t 1’]; %线性核函数
cmd2 =[’-t 1 –c 100’]; %线性核函数
tic;
model = svmtrain(PriviaU_train_lable, PriviaU_train_data, cmd1);
predict_result=svmpredict(PriviaU_test_lable,PriviaU_test_data,model);
time=toc;
disp([‘time=’,num2str(time)]);
% 结果分析

% 进行各类样本的分类样本计数

% 计算Kappa系数

disp([‘kappa=’,num2str(kappa)]);
% 各类样本分类精度
accuracy_single=predict_class_num./single_num;
disp(accuracy_single);
disp(kappa_class);
disp(all_num);
%% 绘制分类图像
% 按索引恢复图像分类矩阵
normal_class_mat=zeros(img_row,img_col);
PriviaU_class_line=zeros(img_row*img_col,1);
PriviaU_class_tag=zeros(size(PriviaU_tag,1),1);
for m=1:size(PriviaU_train,1)
PriviaU_class_tag(PriviaU_train(m,feature_num+2),1)=PriviaU_train(m,feature_num+1);
end
for n=1:size(PriviaU_test,1)
PriviaU_class_tag(PriviaU_test(n,feature_num+2),1)=predict_result(n,1);
end
PriviaU_class_line(data_index,1)=PriviaU_class_tag;
% 按顺序恢复图像
m=1;
for i=1:img_row
for j=1:img_col
normal_class_mat(i,j)=PriviaU_class_line(m,1);
m=m+1;
end
end% 绘制原始标签图
ordinary_img=draw_picture(PaviaU_gr);
figure
imshow(ordinary_img);
title(‘PriviaU Data Ordinary Image’);
% 绘制分类全图
label_img=draw_picture(normal_class_mat);
figure
imshow(label_img);
title(‘PriviaU Data Classification Image’);
%% 用模型进行全图分类
all_predict_result=svmpredict(PriviaU_tag_all,PriviaU_table_all,model);
all_predict_mat=reshape(all_predict_result,img_col,img_row);
all_predict_img=draw_picture(all_predict_mat’);
figure
imshow(all_predict_img);
title(‘All PriviaU Data Classification Image’);