@遥感图像光谱通道分类（Matlab）
基于遥感图像光谱通道的地物分类（Matlab）
一.实践数据与目的
1.高光谱影像：
高光谱影像的光谱通道数通常多达数十甚至数百个以上，而且各光谱通道间往往是连续的，因此能产生一条完整而连续的光谱曲线，对于每一个像素点，将其各个光谱通道上的灰度值作为其一个特征值，就可以建立一幅高光谱图像的像素点的特征空间。这就为以像素点的各光谱通道值作为特征，将像素点进行分类提供了理论基础，也就为高光谱影像的地物分类提供了原理。
此次试验的数据来源于ROSIS 高光谱传感器，影像的拍摄地点位于意大利，分别为：
（1）Pavia University，地点为意大利Pavia University，影像大小为610×340 pixels，影像共有103连续波段，影像的空间分辨率为1.3 米，主要地物有沥青路，草地，沙砾，树木，裸土地，柏油房等；
（2）Pavia Center，地点为意大利Pavia Center，影像大小为1096×715 pixels，影像共有103连续波段，影像的空间分辨率为1.3 米，主要地物有河水，树木，草地，砖块，裸土地，柏油房，沥青路等。
2.实验目的：
（1）实验内容：
利用以下三种方法进行分类处理，并比较这三种方法进行分类的精度比较：
1）利用原始高维的光谱特征作为影像特征，进行影像分类；
2）利用PCA降维方法对原始高维光谱特征进行降维，然后利用降维后的光谱特征进行影像分类；
3）在方法2）的基础上，加入局部窗口的均值和方差特征，进行影像分类；
（注：分类器统一选用支持向量机(SVM)）
（2）实验步骤：
1）从已标注的样本标签中每类选择X%的样本用于分类器训练，剩余的用于分类器精度测试（X 取 5%，10%，15%，20%，25%）；
2）评价指标：每类地物的正确分类精度，总体分类精度，kappa系数；
3）最后利用训练好的分类器对全图进行分类。
二.分类原理
1.样本划分：
对于模型的训练，需要对已知样本进行分类，将其分为训练集与测试集两部分，前者用于训练出模型，后者用于检测训练出模型的精确度，这当中的划分方法就决定了模型的训练精确度。此处采用留出法对样本进行划分，而划分的原理为：直接将整个数据集D划分为两个互斥的集合，其中一个集合S作为训练集，另一个作为测试集T。实验中对划分为两个互斥的集合的做法是：利用K均值算法将整个数据集强行划分为两个互斥的类，考虑到直接用这两个集合来分别进行训练和测试任务过于极端化，故按照一定的比例分别从两个类当中取出一定比例的数据作为训练集，剩下一部分作为测试集样本。
2.K均值算法：
该算法的思想是:基于使聚类性能指标最小化，所用的聚类准则函数是聚类集中每一一个样本点到该类中心的距离平方之和，并使其最小化。
3.主成分分析法（PCA）：
主成分分析法在地物分类中的应用主要体现在对于每个像素的所有像素通道，使用PCA方法查找出整个通道中起着重要作用的光谱层，这些光谱层就是描述这个像素最主要部分的代表性特征，我们也就可以将其看作区分不同地物的像素的主要特征。而PCA会将所有的光谱层特征向量综合压缩成一个数量更少的特征向量，在保留了特征值的区分度的同时，还降低了分类的工作量，使得分类工作更加快速有效的进行。
贡献率：使用主成分分析法时，每一个特征都会求出一个代表其特征效果的特征值，对于图像就表现为每一个像素点都有一个对应的特征值，其值越大代表其在图像中起到的特征效果越强，越可以用于代表整幅图像的特征，而某个特征值对于所有特征值之和的比率即代表了该特征值的贡献率。
4.Kappa系数：
Kappa系数是一种计算分类精度的方法。它是通过把所有地表真实分类中的像元总数(N)乘以混淆矩阵对角线(Xkk)的和，再减去某一类地表真实像元总数与被误分成该类像元总数之积对所有类别求和的结果，再除以总像元数的平方减去某一类中地表真实像元总数与该类中被误分成该类像元总数之积对所有类别求和的结果所得到的。
Kappa系数的计算结果为-1~1，但通常kappa是落在 0~1 间，可分为五组来表示不同级别的一致性：0.00.20极低的一致性(slight)、0.210.40一般的一致性(fair)、0.41~0.60 中等的一致性(moderate)、0.61~0.80 高度的一致性(substantial)和0.81~1几乎完全一致(almost perfect)。
三.实验代码
函数1：Chazhao.m
函数功能：输入一个m行n列的数字矩阵，查找每一行中非0元素的个数，返回一个m行的数量矩阵answer
函数代码：
function answer = Chazhao(aim)
zhong=0;
[m,n]=size(aim);
answer=zeros(m,1);
for i=1:m
for j=1:n
if aim(i,j)~=0
zhong=zhong+1;
end
end
answer(i,1)=zhong;
zhong=0;
end
end
函数2：Nummatching.m
函数功能：查找数字a在列矩阵aim的第多少行，返回行数answer，若不存在则返回0，若存在多个则返回第一次出现的行数
函数代码：
function answer = Nummatching(a,aim)
[m,~]=size(aim);
answer=0;
for i=1:m
if a==aim(i,1)
answer=i;
break;
end
end
end
函数3：Reorder.m
函数功能：排序函数，输入一个列矩阵a，将a中的元素按照要求进行升序（mark为2）或降序（mark为1）排列好放到result矩阵中，并将各元素在原来矩阵中的行数返回到name矩阵中依次对应进行保存
函数代码：
function [answer,name] = Reorder(a,mark)
[m,~]=size(a);%获得输入的数据的数量
answer=zeros(m,1);%初始化结果矩阵
name=zeros(m,1);
[mum,t1]=max(a);%先获得所有数据之中的最大值位置
[mim,t2]=min(a);%先获得所有数据之中的最小值位置
%先按降序进行处理
answer(1,1)=mum;name(1,1)=t1;
answer(m,1)=mim;name(m,1)=t2;
yuzhi=mim-1;%设置一个阈值方便改变元素值
a(t1,1)=yuzhi;a(t2,1)=yuzhi;
if m>2
for i=2:(m-1)
[data,num]=max(a);
answer(i,1)=data;
name(i,1)=num;
a(num,1)=yuzhi;
end
end
if mark==2%若要求升序排列
mid1=answer;
mid2=name;
for i=1:m
answer(i,1)=mid1(m+1-i,1);
name(i,1)=mid2(m+1-i,1);
end
end
end
函数4：KMeanvalue.m
函数功能：K均值聚类函数,输入一个含有m个样本,每个样本n个特征的m*n矩阵aim，再输入一个期望聚类中心数k，将返回一个k行的矩阵answer，其中每一行代表一个聚类，每一行中的非零元素表示被分配到这一聚类中的样本序号，0是为了填充矩阵所加入的，无实际意义; aim是输入的样本矩阵（需要按序号顺序输入），k是聚类中心数量，answer是返回的分类结果
函数代码：
function answer = KMeanvalue(aim,k)
[m,n]=size(aim);%m是输入的样本数量，n是聚类的特征数量
Julei=aim(1:k,1:n);%定义一个矩阵，用于储存聚类中心的特征值,并将输入的前k个样本作为初始的聚类中心
%针对初始聚类中心有重复的情况进行处理
gaizheng=1;
panduan=0;
while panduan==0
panduan=1;
for s1=1:k-1
a1=Julei(s1,1:n);
for s2=s1+1:k
a2=Julei(s2,1:n);
if isequal(a1,a2)==1
Julei(s2,1:n)=aim(k+gaizheng,1:n);
gaizheng=gaizheng+1;
panduan=0;
end
end
end
end %处理完成
Julei1=zeros(k,n);%定义一个矩阵，用于储存一次处理后聚类中心的特征值
Dist=zeros(m,k);%定义一个距离矩阵，用于后面储存样本与聚类中心的距离
result=0;
shumu=zeros(k,1);%定义一个矩阵储存各个聚类中的样本数量
Number=zeros(k,m);%定义一个矩阵用于储存各个聚类中的样本序号
while isequal(Julei,Julei1)==0%判断矩阵Julei与Julei1是否相等
Julei1=Julei;
Number=0*Number;%重置样本分类
for i=1:m
for j=1:k
for num=1:n
result = result + (aim(i,num)-Julei(j,num))*(aim(i,num)-Julei(j,num));
end
Dist(i,j)=sqrt(result);
result=0;
end
end
mid=transpose(Dist);
[~,u]=min(mid);%利用转置和min函数找出距离矩阵中每一行距离值最小值的列数，储存于行向量u之中
u1=transpose(u);
for r=1:k
midd=find(u1==r);%返回一个列向量，储存列向量u1中元素值等于r的元素的列数
[j,~]=size(midd);
shumu(r,1)=j;
midd1=transpose(midd);
Number(r,1:shumu(r,1))= midd1;
end
guodu=zeros(1,n);%建立中间矩阵，方便后面运算
for num1=1:k
for num2=1:m
if(u(1,num2)==num1)
guodu1=aim(num2,1:n);
guodu=guodu+guodu1;
end
end
Julei(num1,1:n)=guodu/shumu(num1,1);
guodu=0;
end
end
mixmum=max(shumu);
answer1=Number(1:k,1:mixmum);
answer=answer1;
end
函数5：Classify.m
函数功能：输入一个列向量，分辨出列向量中不同数字的数量m,并且将相同的元素归类;返回一个m行1列矩阵answer1，每一行元素是输入列向量中不同的数字，还有一个m行n列矩阵answer2，每行元素是answer1中对应行元素在输入列向量中的行号，不足处用0代替
函数代码：
function [answer1,answer2] = Classify(aim)
[m,~]=size(aim);%读取列向量的行数
mid=zeros(m,1);%过渡矩阵,储存不同的数字值
sumnum=zeros(m,1);%定义各不同元素的数目
%xuhao=zeros(m,m);%定义各不同元素对应的样本序号
xuhao=zeros(ceil(0.01*m),ceil(0.5*m));%当输入的样本数量过于大时，使用此式进行定义
difnum=1;%定义总共不同的数字的数目
mid(1,1)=aim(1,1);
sumnum(1,1)=1;
xuhao(1,1)=1;
for i=2:m %通过对比得出总的不同元素的数目以及各自元素的数目
panduan=0;
for j=1:difnum
if aim(i,1) == mid(j,1)
sumnum(j,1)=sumnum(j,1)+1;
xuhao(j,sumnum(j,1))=i;
panduan=1;
end
end
if panduan==0
difnum=difnum+1;
mid(difnum,1)=aim(i,1);
sumnum(difnum,1)=1;
xuhao(difnum,1)=i;
end
end
answer1=mid(1:difnum,1);
maxnum=max(sumnum);
answer2=xuhao(1:difnum,1:maxnum);
end
函数6：Assort.m
函数功能：样本分类函数，将输入的样本分为训练集和测试集，examples为输入的样本矩阵，每一行是一个样本，每一行的元素是样本的特征值；name是对应样本的列向量标签集，每一个元素与examples的每一行的样本对应；ratio为比例系数，表示理想的训练集占总样本数量的比例;train和test分别为返回的训练集与测试集矩阵,kind是总的k行1列样本标签，train和test的最后一列元素为1—k的整数m，分别表示该行样本的标签是kind中的第m行的标签由于此函数使用留出法对样本进行分离
函数代码：
function [train,test,kind] = Assort(examples,name,ratio)
[~,n]=size(examples);%获得输入的样本特征数量n
%利用Classify函数将样本进行聚类，相同标签的样本分到一起
%middle1为一个列向量，每一个元素表示对应行元素的标签
%middle2中每一行储存的是middle1中对应行的标签所包含的样本序号
[middle1,middle2]=Classify(name);
kind=middle1;
[zhong,~]=size(middle1);%获得总共的样本种类zhong
num=Chazhao(middle2);%获得每一种类样本的数量
su=1;
Examples=examples;%将样本按归类依次重新排序
for s=1:zhong
for t=1:num(s,1)
if middle2(s,t)~=0
examples(su,:)=Examples(middle2(s,t),:);
su=su+1;
end
end
end
number1=num;
for i=1:zhong
number1(i,1)=ceil(number1(i,1)*ratio);%储存对应样本应当拥有的训练集样本数目
end
number2=num-number1;%储存对应样本应当拥有的测试集样本数目
xunlian=sum(number1);%获得训练集样本数量
ceshi=sum(number2);%获得测试集样本数量
train=zeros(xunlian+10,n+1);%初始化训练集样本,+10是处于取整函数考虑，保留一定的容错空间
test=zeros(ceshi+10,n+1);%初始化测试集样本，+10是处于取整函数考虑，保留一定的容错空间
%下面用留出法将总样本分为训练集train与测试集test两部分
xu1=1;xu2=0;xu3=0;xun=0;ce=0;%定义中间选择系数
for t=1:zhong
xu2=xu2+num(t,1);
purpose=examples(xu1:xu2,:);
xu1=xu1+num(t,1);
if number2(t,1)~=0
answer=KMeanvalue(purpose,2);%利用K均值函数将每一类样本分为正反两个部分
a1=Chazhao(answer);%获得两部分的样本数量
a2=a1;
for i=1:2
a2(i,1)=ceil(a2(i,1)*ratio);%储存训练集样本数目
end
xuhao1=randi(a1(1,1),a2(1,1),1);
xuhao2=randi(a1(2,1),a2(2,1),1);
for p=1:a2(1,1)%给训练集赋值
xun=xun+1;
shu=xuhao1(p,1);
train(xun,1:n)=examples(answer(1,shu)+xu3,:);
train(xun,1+n)=t;
end
for p=1:a2(2,1)
xun=xun+1;
shu=xuhao2(p,1);
train(xun,1:n)=examples(answer(2,shu)+xu3,:);
train(xun,1+n)=t;
end
for q=1:a1(1,1)%给测试集赋值
if Nummatching(q,xuhao1)==0
ce=ce+1;
test(ce,1:n)=examples(answer(1,q)+xu3,:);
test(ce,1+n)=t;
end
end
for q=1:a1(2,1)
if Nummatching(q,xuhao2)==0
ce=ce+1;
test(ce,1:n)=examples(answer(2,q)+xu3,:);
test(ce,1+n)=t;
end
end
else
for i=1:num(t,1)
xun=xun+1;
train(xun,1:n)=purpose(i,:);
train(xun,1+n)=t;
ce=ce+1;
test(ce,1:n)=purpose(i,:);
test(ce,1+n)=t;
end
end
xu3=xu3+num(t,1);
end
train=train(1:xun,:);
test=test(1:ce,:);
end
函数7：Maincomonent.m
函数功能：主成分分析方法函数,data为输入的m*n的训练样本集，m行每一行为一个样本的特征值，n列是每一个样本的n个特征值;standard为主成分的理想贡献率，一般取大于0.8的数;answer为返回的一个结果矩阵，为t（t<=n)行n列,用于对样本进行K-L变换，t表示变换后的维度
函数代码：
function [answer,t]=Maincomonent(data,standard)
[m,n]=size(data);%先获得读入的样本矩阵尺寸
mid1=zeros(1,n);
for i=1:m
mid1=mid1+data(i,:);
end
M=mid1/m;%获得均值向量
mid2=zeros(n,m);%定义中间样本均值差矩阵
mid3=data';
for i=1:m
mid2(:,i)=mid3(:,i)-M';%将均值样本按列赋值
end
C=mid2*(mid2')/m;%得到协方差矩阵
[T1,T2]=eig(C);%利用eig函数求出协方差矩阵C的特征值矩阵T2与特征向量矩阵T1
spe=zeros(n,1);%定义特征值储存矩阵
for i=1:n
spe(i,1)=T2(i,i);
end
[special,line]=Reorder(spe,1);%利用Reorder函数将特征值按降序排列好，并将原本序号储存于line中
he=sum(special);%计算所有特征值的和
cal=0;%定义有实际作用的特征值的数量
midd=0;
for i=1:n
midd=midd+special(i,1);
cal=cal+1;
if (midd/he)>=standard
break;
end
end
final2=zeros(cal,n);%定义一个矩阵用于储存有效的特征向量
for i=1:cal
final2(i,:)=T1(line(i,1),:);
end
answer=final2;
t=cal;
end
函数8：Svmdealwith.m
函数功能：SVM模型训练函数，input1与input2分别为输入的训练样本集与测试样本集，其中每一行代表一个样本，其中每行样本的最后一列元素值表示该样本的标签序号，其他元素表示样本特征值; name为标签矩阵，每个序号对应相应行的标签名称； model为返回的训练模型,percent1为训练集测试精度矩阵，percent2为测试集测试精度矩阵； kappa1和kappa2分别为训练集和测试集样本的Kappa系数
函数代码：
function [percent1,percent2,kappa1,kappa2,model] = Svmdealwith(input1,input2,name)
[~,n]=size(input2);%获得样本的特征数量(n-1)
data=input1(:,1:n-1);%获得每个训练样本的特征值数据
feature1=input1(:,n);%获得每个训练样本的标签序号
test=input2(:,1:n-1);%获得每个测试样本的特征值数据
feature2=input2(:,n);%获得每个测试样本的标签序号
kind=size(name,1);%获得总的标签总数
cmd=('-t 1 -c 100');%函数类型可选参数，此处选择为线性核函数
model=svmtrain(feature1,data,cmd);%获得训练的模型
%下面利用训练的模型进行精度检测
%训练集
[lei1,xu1]=Classify(feature1);%获得各类训练样本的标签序号数量特征
number=Chazhao(xu1);%获得各标签序号的数量
testmark=svmpredict(feature1,data,model);%检验模型正确率
accarcy=zeros(kind,2);%建立一个训练样本正确率统计矩阵
hunxiao1=zeros(kind,kind);%建立一个混淆矩阵用于计算Kappa系数
for i=1:kind %依次将第一列赋值为标签序号
accarcy(i,1)=i;
end %依次将第二列赋值为分类正确的标签的数量
for i=1:size(data,1)
if feature1(i,1)==testmark(i,1)
accarcy(feature1(i,1),2)=accarcy(feature1(i,1),2)+1;
end
hunxiao1(feature1(i,1),testmark(i,1))=hunxiao1(feature1(i,1),testmark(i,1))+1;
end
for i=1:kind
j=Nummatching(i,lei1);
accarcy(i,2)=accarcy(i,2)/number(j,1);
end
p11=0;p21=0;
for i=1:kind
p11=p11+hunxiao1(i,i);%计算对角线元素之和
p21=p21+sum(hunxiao1(i,:))*sum(hunxiao1(:,i));
end
p31=p11/sum(number);
p41=p21/(sum(number)*sum(number));
kappa1=(p31-p41)/(1-p41);
percent1=accarcy;
%测试集
[lei2,xu2]=Classify(feature2);%获得各类测试样本的标签序号数量特征
number1=Chazhao(xu2);%获得各标签序号的数量
answer0=svmpredict(feature2,test,model);
answer=zeros(kind,2);%建立一个测试样本正确率统计矩阵
hunxiao2=zeros(kind,kind);%建立一个混淆矩阵用于计算Kappa系数
for i=1:kind
answer(i,1)=i;
end
for i=1:size(test,1)
if feature2(i,1)==answer0(i,1)
answer(feature2(i,1),2)=answer(feature2(i,1),2)+1;
end
hunxiao2(feature2(i,1),answer0(i,1))=hunxiao2(feature2(i,1),answer0(i,1))+1;
end
for i=1:kind
j=Nummatching(i,lei2);
answer(i,2)=answer(i,2)/number1(j,1);
end
p12=0;p22=0;
for i=1:kind
p12=p12+hunxiao2(i,i);%计算对角线元素之和
p22=p22+sum(hunxiao2(i,:))*sum(hunxiao2(:,i));
end
p32=p12/sum(number1);
p42=p22/(sum(number1)*sum(number1));
kappa2=(p32-p42)/(1-p42);
percent2=answer;
函数9：Design.m
函数功能：这是三种方法都有的一个主函数，这当中调用了前面所有的函数来完成一个整体的实验功能，包括数据处理，模型训练，精度测试以及结果图的绘制都将由此函数来完成，由于三种不同的方法中设计的Design函数的调用格式也略微存在一些不同，故下面将分开进行介绍
函数代码：
三种方法的Design函数都可以依次分为数据预处理、数据运算和结果展示三个
部分，其中第一种方法和第二种方法的数据预处理与结果展示两部分的代码完全
相同，方法三由于需要求出局部窗口均值与方差，所以数据预处理与数据运算部
分与前两种方法有一些不同，而结果展示部分也是相同的，具体代码如下所示：
①方法一和方法二的数据预处理：
[u,v,w]=size(aim);
num=0;%定义已知的样本数量
num1=0;%定义待确定的样本数量
for i=1:u %循环求出总的样本数量
for j=1:v
if biaoqian(i,j)~=0
num=num+1;
elseif biaoqian(i,j)==0
num1=num1+1;
end
end
end
data=zeros(num,w);%定义已知数据矩阵
mark=zeros(num,1);%定义已知数据标签矩阵
data1=zeros(num1,w);%定义未知数据矩阵
mark1=zeros(num1,1);%定义未知数据标签矩阵
station=zeros(num,2);%定义已知数据坐标矩阵
station1=zeros(num1,2);%定义未知数据坐标矩阵
%定义颜色矩阵
color=[0,0,255;0,255,0;255,0,0;255,255,0;255,10,155;80,116,189;148,140,82;144,238,144;231,100,55]/255;
t=0;%定义中间序号
t1=0;
for i=1:u
for j=1:v
if biaoqian(i,j)~=0
t=t+1;
data(t,:)=aim(i,j,:);
mark(t,1)=biaoqian(i,j);
station(t,1)=i;
station(t,2)=j;
elseif biaoqian(i,j)==0
t1=t1+1;
data1(t1,:)=aim(i,j,:);
mark1(t1,1)=0;
station1(t1,1)=i;
station1(t1,2)=j;
end
end
end
②三种方法的结果展示：
accuacy1=percent1;%将精度矩阵按照样本标签重新排列
accuacy2=percent2;
zhong=size(kind,1);%获得总的标签种类
for i=1:zhong
j=Nummatching(i,kind);
accuacy1(i,2)=percent1(j,2);
accuacy2(i,2)=percent2(j,2);
end
fprintf('\n下面是训练模型训练集精度矩阵\n');
disp(accuacy1);
fprintf('\n训练集Kappa系数为：');
disp(kappa1);
fprintf('\n下面是训练模型测试集精度矩阵\n');
disp(accuacy2);
fprintf('\n测试集Kappa系数为：');
disp(kappa2);
testanswer=svmpredict(mark1,data1,model);%调用模型进行标签预测
mark1=testanswer;%获得未知样本的标签序号
for i=1:num1 %将未知样本的标签序号转换为标签
mark1(i,1)=kind(mark1(i,1),1);
end
R=zeros(u,v);%定义合成的图像图层
G=zeros(u,v);
B=zeros(u,v);
for i=1:num %先给已知样本像素点赋色
R(station(i,1),station(i,2))=color(mark(i,1),1);
G(station(i,1),station(i,2))=color(mark(i,1),2);
B(station(i,1),station(i,2))=color(mark(i,1),3);
end
for i=1:num1 %再给未知样本像素点赋为背景白色
R(station1(i,1),station1(i,2))=1;
G(station1(i,1),station1(i,2))=1;
B(station1(i,1),station1(i,2))=1;
end
figure(1);
picimg1=cat(3,R,G,B);
imshow(picimg1);
%imwrite(picimg1,'testpic2.jpg');%将绘制结果输出
title('测试集像素分类结果');
for i=1:num1 %再给未知样本像素点赋色
R(station1(i,1),station1(i,2))=color(mark1(i,1),1);
G(station1(i,1),station1(i,2))=color(mark1(i,1),2);
B(station1(i,1),station1(i,2))=color(mark1(i,1),3);
end
figure(2);
picimg=cat(3,R,G,B);
imshow(picimg);
%imwrite(picimg,'answerpic2.jpg');%将绘制结果输出
title('总体像素分类结果');
figure(3);%单独画出地物标注图像
diwu={'沥青路','草地','沙砾','树木','金属板','裸土地','柏油房','砖块','阴影'};
%diwu={'河水','树木','草地','砖块','裸土地','沥青路','柏油房','瓦片','阴影'};
yanse=zeros(1,zhong);
xloca=[0,1,1,0];
for i=1:zhong
yloca=[i-1,i-1,i,i];
yanse(1,i)=fill(xloca,yloca,color(i,:));
hold on;
end
hold off;
legend(yanse,diwu,'Location','northeastoutside');
set(gcf,'position',[100,360,200,400]);
title('颜色分类标注');
③中间的数据处理为：
方法一：（直接调用原始数据）
[train,test,kind] = Assort(data,mark,bili);%分离出训练集train与测试集test
[percent1,percent2,kappa1,kappa2,model]=Svmdealwith(train,test,kind);%利用训练集train与测试集test对模型进行检测
方法二：（PCA变化原始数据后再处理）
sumdata=[data;data1];%将所有的样本综合起来
[answer,~]=Maincomonent(sumdata,gongxian);%获取贡献率得到的主成分分析变化矩阵
datas=(answer*data')';%获得主成分分析法变化得到的已知样本数据
data1s=(answer*data1')';%获得主成分分析法变化得到的未知样本数据
[train,test,kind] = Assort(datas,mark,bili);%分离出训练集train与测试集test
[percent1,percent2,kappa1,kappa2,model]=Svmdealwith(train,test,kind);%利用训练集train与测试集test对模型进行检测
④方法三的数据预处理与数据运算：
[u,v,w]=size(aim);
yuanshidata=zeros(u*v,w);%用于储存总的数据
t=0;
for i=1:u
for j=1:v
t=t+1;
yuanshidata(t,:)=aim(i,j,:);
end
end
[answer,~]=Maincomonent(yuanshidata,gongxian);%获取贡献率得到的主成分分析变化矩阵
maindata=(answer*yuanshidata')';%获得主成分分析法变化得到的原始样本数据
datanum=size(maindata,2);%获得主成分分析法变化得到的数据维数
newgetdata=zeros(u,v,datanum);%建立新矩阵还原主成分分析法变化得到的数据
t=0;
for i=1:u
for j=1:v
t=t+1;
newgetdata(i,j,:)=maindata(t,:);
end
end
num=0;%定义已知的样本数量
num1=0;%定义待确定的样本数量
for i=1:u %循环求出总的样本数量
for j=1:v
if biaoqian(i,j)~=0
num=num+1;
elseif biaoqian(i,j)==0
num1=num1+1;
end
end
end
data=zeros(num,datanum*3);%定义已知数据矩阵
mark=zeros(num,1);%定义已知数据标签矩阵
data1=zeros(num1,datanum*3);%定义未知数据矩阵
mark1=zeros(num1,1);%定义未知数据标签矩阵
station=zeros(num,2);%定义已知数据坐标矩阵
station1=zeros(num1,2);%定义未知数据坐标矩阵
%定义颜色矩阵
color=[0,0,255;0,255,0;255,0,0;255,255,0;255,10,155;80,116,189;148,140,82;144,238,144;231,100,55]/255;
t=0;%定义中间序号
t1=0;
for i=1:u %进行数据分类赋值
for j=1:v
if biaoqian(i,j)~=0
t=t+1;
data(t,1:datanum)=newgetdata(i,j,:);
mark(t,1)=biaoqian(i,j);
station(t,1)=i;
station(t,2)=j;
elseif biaoqian(i,j)==0
t1=t1+1;
data1(t1,1:datanum)=newgetdata(i,j,:);
mark1(t1,1)=0;
station1(t1,1)=i;
station1(t1,2)=j;
end
end
end
bianjie=(dimen-1)/2;
for i=1:num %对已知样本进行均值和方差计算
for j=1:datanum
s1=station(i,1)-bianjie;s2=station(i,1)+bianjie;%定义横轴方向坐标范围
t1=station(i,2)-bianjie;t2=station(i,2)+bianjie;%定义纵轴方向坐标范围
if s2>u %判断是否超出边界
s2=u;
end
if s1<1
s1=1;
end
if t2>v %判断是否超出边界
t2=v;
end
if t1<1
t1=1;
end
aimdata=zeros((s2-s1+1)*(t2-t1+1),1);%计算原始数据的均值与标准差
deltanum=0;
for p=s1:s2
for q=t1:t2
deltanum=deltanum+1;
aimdata(deltanum,1)=newgetdata(p,q,j);
end
end
data(i,datanum+j)=mean(aimdata);%均值赋值
data(i,datanum*2+j)=std(aimdata);%标准差赋值
end
end
for i=1:num1 %对未知样本进行均值和方差计算
for j=1:datanum
s1=station1(i,1)-bianjie;s2=station1(i,1)+bianjie;%定义横轴方向坐标范围
t1=station1(i,2)-bianjie;t2=station1(i,2)+bianjie;%定义纵轴方向坐标范围
if s2>u %判断是否超出边界
s2=u;
end
if s1<1
s1=1;
end
if t2>v %判断是否超出边界
t2=v;
end
if t1<1
t1=1;
end
aimdata=zeros((s2-s1+1)*(t2-t1+1),1);%计算原始数据的均值与标准差
deltanum=0;
for p=s1:s2
for q=t1:t2
deltanum=deltanum+1;
aimdata(deltanum,1)=newgetdata(p,q,j);
end
end
data1(i,datanum+j)=mean(aimdata);%均值赋值
data1(i,datanum*2+j)=std(aimdata);%标准差赋值
end
end
[train,test,kind] = Assort(data,mark,bili);%分离出训练集train与测试集test
[percent1,percent2,kappa1,kappa2,model]=Svmdealwith(train,test,kind);%利用训练集train与测试集test对模型进行检测
四.调参对实验结果的影响：
除了以上针对每一个实验不同输入参数的变化比较，还有一个三个实验中都存在的影响因素，那就是SVM模型参数的影响，此处在方法三的实验函数8：Svmdealwith.m后面添加了一些函数语句，得到了一些调参对结果的影响结果，并输出了一个表格用于展示：
添加的函数语句为：
%以下开始进行调参比较输出结果，并将输出结果输出为一个excel表格进行储存
str1='-t ';str2=' -c ';str3=' -g ';%定义三个参数前缀
shumu1=2;shumu2=180;shumu3=180;%定义三个参数的取值界限
excelshuchu=cell(shumu1*shumu2/30*shumu3/30+1,3+kind);%定义表格的表头
for i=1:kind
excelshuchu(1,i+1)={['地物',num2str(i),'精度']};
end
excelshuchu(1,1)={'向量机参数值'};
excelshuchu(1,2+kind)={'总体精度'};
excelshuchu(1,3+kind)={'Kappa系数'};
hangshu=1;
for p=1:shumu1
for q=30:50:shumu2
for r=30:50:shumu3
hangshu=hangshu+1;
midaccarcy=zeros(kind,2);%建立一个训练样本正确率统计矩阵
midhunxiao1=zeros(kind,kind);%建立一个混淆矩阵用于计算Kappa系数
designexcel=cell(1,3+kind);%用于储存计算结果
midcmd=[str1,num2str(p),str2,num2str(q),str3,num2str(r)];
designexcel(1,1)={midcmd};
midmodel=svmtrain(feature1,data,midcmd);%获得训练的模型
[midtestmark,midpercent1,~]=svmpredict(feature1,data,midmodel);%检验模型正确率
for i=1:kind %依次将第一列赋值为标签序号
midaccarcy(i,1)=i;
end %依次将第二列赋值为分类正确的标签的数量
for i=1:size(data,1)
if feature1(i,1)==midtestmark(i,1)
midaccarcy(feature1(i,1),2)=midaccarcy(feature1(i,1),2)+1;
end
midhunxiao1(feature1(i,1),midtestmark(i,1))=midhunxiao1(feature1(i,1),midtestmark(i,1))+1;
end
for i=1:kind
j=Nummatching(i,lei1);
midaccarcy(i,2)=midaccarcy(i,2)/number(j,1);
end
p11=0;p21=0;
for i=1:kind
p11=p11+midhunxiao1(i,i);%计算对角线元素之和
p21=p21+sum(midhunxiao1(i,:))*sum(midhunxiao1(:,i));
end
p31=p11/sum(number);
p41=p21/(sum(number)*sum(number));
midkappa1=(p31-p41)/(1-p41);%训练集Kappa系数
%下面进行测试集数据统计
[midanswer0,midpercent2,~]=svmpredict(feature2,test,midmodel);
midanswer=zeros(kind,2);%建立一个测试样本正确率统计矩阵
midhunxiao2=zeros(kind,kind);%建立一个混淆矩阵用于计算Kappa系数
for i=1:kind
midanswer(i,1)=i;
end
for i=1:size(test,1)
if feature2(i,1)==midanswer0(i,1)
midanswer(feature2(i,1),2)=midanswer(feature2(i,1),2)+1;
end
midhunxiao2(feature2(i,1),midanswer0(i,1))=midhunxiao2(feature2(i,1),midanswer0(i,1))+1;
end
for i=1:kind
j=Nummatching(i,lei2);
midanswer(i,2)=midanswer(i,2)/number1(j,1);
end
p12=0;p22=0;
for i=1:kind
p12=p12+midhunxiao2(i,i);%计算对角线元素之和
p22=p22+sum(midhunxiao2(i,:))*sum(midhunxiao2(:,i));
end
p32=p12/sum(number1);
p42=p22/(sum(number1)*sum(number1));
midkappa2=(p32-p42)/(1-p42);
accuacy1=midaccarcy;%将精度矩阵按照样本标签重新排列
accuacy2=midanswer;
for i=1:kind
j=Nummatching(i,name);
accuacy1(i,2)=midaccarcy(j,2);
accuacy2(i,2)=midanswer(j,2);
end
designexcel(1,kind+3)={[num2str(midkappa1),'//',num2str(midkappa2)]};
designexcel(1,kind+2)={[num2str(midpercent1(1,1)),'//',num2str(midpercent2(1,1))]};
for i=1:kind
designexcel(1,i+1)={[num2str(accuacy1(i,2)),'//',num2str(accuacy2(i,2))]};
end
excelshuchu(hangshu,:)=designexcel;
end
end
end
xlswrite('参数改变比对表',excelshuchu);
1074
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
