遥感影像监督分类和非监督分类与相关代码实现

...wd......wd......wd...遥感影像监视分类与非监视分类摘要：遥感影像的分类方法按照是否有先验类别可分为监视分类和非监视分类，这两种分类方法有着本质的区别但也有存在一定的联系。本文从分类原理和分类方法等不同角度分别介绍了监视分类和非监视分类方法，并对两种方法的分类结果进展了比照和分析。关键词：遥感；监视分类；非监视分类；ISODATA算法；贝叶斯分类算法。数据来源本文使用的数据是华盛顿广场上空由卫星拍摄的高光谱遥感影像。该幅影像使用的传感器系统覆盖0.4到2.4m的可见光到近红外的210个波段。由于0.9-1.4米光谱对应的区域上空大气透光性很差，因此将这个区域内的波段从影像中删除，最后得到191个波段。该数据集有1208个扫描行，每行307个像元，容量近似150MB。为了清晰地反映影像中地物的特征，本文从191个波段中选择了3个波段，分别是120、140和160。下面两幅图分别是全波段影像和三波段影像：图1全波段影像图2三波段影像遥感影像的监视分类2.1监视分类的原理监视分类〔supervisedclassification〕又称训练场地法，是以建立统计识别函数为理论基础，依据典型样本训练方法进展分类的技术。即根据训练区提供的样本，通过选择特征参数，求出特征参数作为决策规则，建立判别函数以对各待分类影像进展的图像分类，是模式识别的一种方法。要求训练区域具有典型性和代表性。判别准则假设满足分类精度要求，则此准则成立；反之，需重新建立分类的决策规则，直至满足分类精度要求为止。常用算法有：最大释然分类法、最小距离分类法、马氏距离分类法、平行六面体分类法、K-NN分类法。本文选用最大释然分类法对遥感影像进展监视分类。最大似然判别法.也称为贝叶斯(Bayes)分类,是基于图像统计的监视分类法,也是典型的和应用最广的监视分类方法.它建立在Bayes准则的基础上,偏重于集群分布的统计特性,分类原理是假定训练样本数据在光谱空间的分布是服从高斯正态分布规律的,做出样本的概率密度等值线,确定分类,然后通过计算标本(像元)属于各组(类)的概率,将标本归属于概率最大的一组.用最大似然法分类,具体分为三步:首先确定各类的训练样本,再根据训练样本计算各类的统计特征值,建立分类判别函数,最后逐点扫描影像各像元,将像元特征向量代入判别函数,求出其属于各类的概率,将待判断像元归属于最大判别函数值的一组。2.2最大释然分类算法的调试#include<iostream.h>#include<fstream.h>#include<stdlib.h>#include<math.h>#include<iomanip.h>constfloatPI=3.1415;floattrain1[60][3],train2[60][3],train3[60][3],train4[60][3];//存各类训练样本floatm1[3],m2[3],m3[3],m4[3],c1[3][3],c2[3][3],c3[3][3],c4[3][3];//各类均值向量及协方差矩阵floatp=0.25;//先验概率floattest[240][11];//检验样本各列分别存放：分类前类别，行号，列号，一波段灰度，二波段灰度，三波段灰度，分类后类别floata1,a2,a3,a4;//公式中常数项；floatqq,qq1,qq2,qq3,qq4;//分类精度floatk;//kappa值intq[4][4];//混淆矩阵floata[4][4];//精度指标矩阵各列分别为生产者精度漏分误差//用户精度错分误差intmain(){ voidgetdata1(char*filename); voidgetdata2(char*filename);voidcompute(floattrain[60][3],floatm[3],floatc[3][3],float*a); floatidentify(floatm[3],floatc[3][3],floata,floatband1,floatband2,floatband3); intclassify(floatp1,floatp2,floatp3,floatp4);voidtesting(); voidaccuracy(floatp[4][4]); voidoutput(char*filename); voidoutput1(char*filename); getdata1("05training.txt");//存放的是训练样本三个波段的灰度值 getdata2("05testing.txt");//检验样本行号列号类别代码三个波段的灰度值compute(train1,m1,c1,&a1);//计算各类均值向量及协方差矩阵 compute(train2,m2,c2,&a2); compute(train3,m3,c3,&a3); compute(train4,m4,c4,&a4);testing();//对检验样本分类 accuracy(a); output("05test-result.txt"); output1("05confusion-matrix.txt"); return0;}voidgetdata1(char*filename){ ifstreaminfile(filename,ios::in|ios::nocreate);//把每类样本的灰度值依次存入对应的数组if(infile==0) {cout<<"openinfileerror!"<<endl; exit(1); }inti,j; for(i=0;i<60;i++) for(j=0;j<3;j++) infile>>train1[i][j];for(i=0;i<60;i++) for(j=0;j<3;j++) infile>>train2[i][j];for(i=0;i<60;i++) for(j=0;j<3;j++) infile>>train3[i][j];for(i=0;i<60;i++) for(j=0;j<3;j++) infile>>train4[i][j]; infile.close();}voidgetdata2(char*filename){ ifstreaminfile(filename,ios::in|ios::nocreate);if(infile==0) { cout<<"openinfileerror!"<<endl; exit(1); }inti,j; for(i=0;i<240;i++)//检验样本的数量为240 for(j=0;j<6;j++)//行、列号、类别代码、三个波段灰度值 infile>>test[i][j];infile.close();}voidcompute(floattrain[60][3],floatm[3],floatc[3][3],float*a){ voidgauss(intf,floatn[3][3]);//求n矩阵的逆阵声明 inti,j;for(i=0;i<3;i++)m[i]=0; for(i=0;i<3;i++) { for(j=0;j<60;j++) m[i]+=train[j][i]; m[i]=m[i]/60; }floatTrain[3][60];//转置矩阵 for(i=0;i<3;i++) for(j=0;j<60;j++) Train[i][j]=train[j][i]; intk; for(i=0;i<3;i++)//求协方差 for(j=0;j<3;j++) {c[i][j]=0;for(k=0;k<60;k++) c[i][j]+=((Train[i][k]-m[i])*(train[k][j]-m[j])); c[i][j]=c[i][j]/60; } *a=c[0][0]*c[1][1]*c[2][2]+c[0][1]*c[1][2]*c[2][0]+c[0][2]*c[1][0]*c[2][1] -c[0][2]*c[1][1]*c[2][0]-c[0][0]*c[1][2]*c[2][1]-c[0][1]*c[1][0]*c[2][2]; *a=(-0.5)*log(*a)+log(p)-25*log(2*PI);gauss(3,c); }floatidentify(floatm[3],floatc[3][3],floata,floatband1,floatband2,floatband3){//计算待分样本的d(X) floatx[3];x[0]=band1;x[1]=band2;x[2]=band3;inti,j; for(i=0;i<3;i++) x[i]-=m[i];floatd[3]; for(i=0;i<3;i++) {d[i]=0; for(j=0;j<3;j++) d[i]+=(x[i]*c[j][i]); } doubleq=0; for(i=0;i<3;i++)q+=x[i]*d[i]; q=(-0.5)*q+a; q=exp(q); returnq;}intclassify(floatp1,floatp2,floatp3,floatp4)//根据后验概率将像元归为指定类别{floatp[4];p[0]=p1;p[1]=p2;p[2]=p3;p[3]=p4;inti,k=0;floatmax=p[0];for(i=1;i<4;i++)if(max<p[i]) {max=p[i]; k=i; }k++;returnk;}voidtesting(){ doublep1,p2,p3,p4,p5,p6,p7; floatkappa(intk[4][4]);inti,j; for(i=0;i<240;i++) { p1=identify(m1,c1,a1,test[i][3],test[i][4],test[i][5]);p2=identify(m2,c2,a2,test[i][3],test[i][4],test[i][5]); p3=identify(m3,c3,a3,test[i][3],test[i][4],test[i][5]);p4=identify(m4,c4,a4,test[i][3],test[i][4],test[i][5]);p5=identify(m5,c5,a5,test[i][3],test[i][4],test[i][5]);p6=identify(m6,c6,a6,test[i][3],test[i][4],test[i][5]);p7=identify(m7,c7,a7,test[i][3],test[i][4],test[i][5]); doubledx=p1+p2+p3+p4+p5+p6+p7; p1=p1/dx;p2=p2/dx; p3=p3/dx;p4=p4/dx;p5=p5/dx;p6=p6/dx;p7=p7/dx; test[i][7]=p1; test[i][8]=p2; test[i][9]=p3; test[i][10]=p4;test[i][11]=p5;test[i][12]=p6;test[i][113]=p7; test[i][6]=classify(p1,p2,p3,p4,p5,p6,p7); //cout<<setiosflags(ios::fixed)<<setprecision(3)<<p1<<""<<p2<<""<<p3<<""<<p4<<""<<p5<<""<<p6<<""<<p7<<endl; } for(i=0;i<7;i++) for(j=0;j<7;j++) q[i][j]=0;for(i=0;i<60;i++)//混淆矩阵中的各类别组合数量 { if(test[i][2]==1&&test[i][6]==1)q[0][0]++; if(test[i][2]==1&&test[i][6]==2)q[0][1]++; if(test[i][2]==1&&test[i][6]==3)q[0][2]++; if(test[i][2]==1&&test[i][6]==4)q[0][3]++;if(test[i][2]==1&&test[i][6]==4)q[0][4]++;if(test[i][2]==1&&test[i][6]==4)q[0][5]++;if(test[i][2]==1&&test[i][6]==4)q[0][6]++; }for(i=60;i<120;i++) { if(test[i][2]==2&&test[i][6]==1)q[1][0]++; if(test[i][2]==2&&test[i][6]==2)q[1][1]++; if(test[i][2]==2&&test[i][6]==3)q[1][2]++; if(test[i][2]==2&&test[i][6]==4)q[1][3]++;if(test[i][2]==2&&test[i][6]==4)q[1][4]++;if(test[i][2]==2&&test[i][6]==4)q[1][5]++;if(test[i][2]==2&&test[i][6]==4)q[1][6]++; } for(i=120;i<180;i++) { if(test[i][2]==3&&test[i][6]==1)q[2][0]++; if(test[i][2]==3&&test[i][6]==2)q[2][1]++; if(test[i][2]==3&&test[i][6]==3)q[2][2]++; if(test[i][2]==3&&test[i][6]==4)q[2][3]++;if(test[i][2]==3&&test[i][6]==4)q[2][4]++;if(test[i][2]==3&&test[i][6]==4)q[2][5]++;if(test[i][2]==3&&test[i][6]==4)q[2][6]++; } for(i=180;i<240;i++) { if(test[i][2]==4&&test[i][6]==1)q[3][0]++; if(test[i][2]==4&&test[i][6]==2)q[3][1]++; if(test[i][2]==4&&test[i][6]==3)q[3][2]++; if(test[i][2]==4&&test[i][6]==4)q[3][3]++;if(test[i][2]==3&&test[i][6]==4)q[3][4]++;if(test[i][2]==3&&test[i][6]==4)q[3][5]++;if(test[i][2]==3&&test[i][6]==4)q[3][6]++;}qq1=(float)q[0][0]/60.0;qq2=(float)q[1][1]/60.0;qq3=(float)q[2][2]/60.0;qq4=(float)q[3][3]/60.0;qq5=(float)q[4][3]/60.0;qq6=(float)q[5][3]/60.0;qq7=(float)q[6][3]/60.0;qq=(qq1+qq2+qq3+qq4+qq5+qq6+qq7)/4; qq=qq*100;k=kappa(q);// cout<<"1代表'房顶'，2代表'街道'，3代表'小路'，4代表'草地'，5代表'树木'，6代表'水体'，7代表'阴影'"<<endl;// cout<<"正确分类百分比分别为：";//cout<<qq1<<""<<qq2<<""<<qq3<<""<<qq4<<""<<qq5<<""<<qq6<<""<<qq7<<endl;//cout<<"总分类精度:"<<setprecision(2)<<qq<<"%"<<endl;// cout<<"Kappa值为:"<<setprecision(2)<<k<<"%"<<endl;}voidoutput(char*filename){ofstreamoutfile(filename,ios::out);if(outfile==0){cout<<"openoutfileerror!"<<endl;exit(1);}outfile<<"1代表'房顶'，2代表'街道'，3代表'小路'，4代表'草地'，5代表'树木'，6代表'水体'，7代表'阴影'"<<endl;outfile<<"正确分类百分比分别为：";outfile<<setprecision(3)<<qq1<<""<<qq2<<""<<qq3<<""<<qq4<<endl;outfile<<"总分类精度:"<<setprecision(4)<<qq<<"%"<<endl;outfile<<"row"<<"col"<<"before"<<"after"<<"proba1"<<"proba2"<<"proba3"<<"proba4"<<endl;inti;for(i=0;i<240;i++){ outfile<<setiosflags(ios::fixed)<<setprecision(0)<<test[i][0]<<""<<test[i][1]<<""<<test[i][2]<<""<<test[i][6]<<"";outfile<<setiosflags(ios::fixed)<<setprecision(3)<<test[i][7]<<""<<test[i][8]<<""<<test[i][9]<<""<<test[i][10]<<endl;} outfile.close();}voidgauss(intf,floatn[3][3]){inti,j,k;floatp,q,h[3];for(k=f-1;k>=0;k--) {p=n[0][0];for(i=1;i<f;i++) {q=n[i][0]; h[i]=((i>k)?q:-q)/p; for(j=1;j<=i;j++) n[i-1][j-1]=n[i][j]+q*h[j];}n[f-1][f-1]=1/p;for(i=1;i<f;i++) n[f-1][i-1]=h[i]; }for(i=1;i<f;i++)for(j=0;j<i;j++)n[j][i]=n[i][j];}floatkappa(intk[4][4]){ inti,j; floatK; intc[4],d[4];intN=0,M=0,D=0; for(i=0;i<4;i++) { c[i]=0; d[i]=0; } for(i=0;i<4;i++) for(j=0;j<4;j++) { N+=k[i][j]; if(i==j)D=D+k[i][j]; c[i]=k[i][j]+c[i]; d[j]=k[i][j]+d[j]; } for(i=0;i<4;i++) M=c[i]*d[i]+M; K=float((N*D-M))/(N*N-M); returnK;}voidaccuracy(floatp[4][4]){ inta[4],b[4]; floata1[4],a2[4],b1[4],b2[4]; inti,j; for(i=0;i<4;i++) { a[i]=0;b[i]=0; a1[i]=0;b1[i]=0; a2[i]=0;b2[i]=0; } for(i=0;i<4;i++) for(j=0;j<4;j++) { p[i][j]=0; a[i]+=q[i][j]; b[i]+=q[j][i]; } for(i=0;i<4;i++) { a1[i]=(float)q[i][i]/a[i]; a2[i]=1-a1[i]; b1[i]=(float)q[i][i]/b[i]; b2[i]=1-b1[i]; }for(i=0;i<4;i++) { p[i][0]=a1[i]*100;//生产者精度 p[i][1]=a2[i]*100;//漏分误差 p[i][2]=b1[i]*100;//用户精度 p[i][3]=b2[i]*100;//错分误差 }}voidoutput1(char*filename){ ofstreamoutfile(filename,ios::out);if(outfile==0){cout<<"openoutfileerror!"<<endl;exit(1);}outfile<<"Lakes"<<"Built-upLand"<<"ForestLand"<<"Cropland"<<endl;inti,j;for(i=0;i<4;i++){ switch(i) { case0: outfile<<"Lakes";break; case1: outfile<<"Built-upLand";break; case2: outfile<<"ForestLand";break; case3: outfile<<"Cropland";break; } for(j=0;j<4;j++) { outfile<<setiosflags(ios::left)<<setw(7)<<q[i][j]<<""; if(j!=0&&j%3==0)outfile<<endl; } }outfile<<endl;outfile<<endl;outfile<<endl;outfile<<""<<setw(11)<<"Producer's"<<setw(11)<<"Omission"<<"User's"<<"Commision"<<endl; for(i=0;i<4;i++) { switch(i) { case0: outfile<<"Lakes";break; case1: outfile<<"Built-upLand";break; case2: outfile<<"ForestLand";break; case3: outfile<<"Cropland";break; } for(j=0;j<4;j++) { outfile<<""<<setiosflags(ios::fixed)<<setprecision(2)<<setw(4)<<a[i][j]<<"%"<<""; if(j!=0&&j%3==0)outfile<<endl; } } outfile<<endl; outfile<<endl; outfile<<""<<"总分类精度:"<<setprecision(2)<<qq<<"%"<<endl; outfile<<""<<"Kappa系数:"<<setprecision(3)<<k<<endl; outfile.close();}2.3最大释然分类算法的执行在贝叶斯分类器中执行以上代码，得到的监视分类图如下所示：图三监视分类图利用最大释然分类算法对遥感影像进展监视分类，最后将图像上的地类分为7类，分别是：房顶，街道，小路，草地，树木，水体和阴影。从处理的图像来看，分类结果较为很理想，基本上可以识别出这7类地物。非监视分类非监视分类的原理非监视分类是以不同影像地物在特征空间中类别特征的差异为依据的一种无先验〔〕类别标准的图像分类，是以集群为理论基础，通过计算机对图像进展集聚统计分析的方法。根据待分类样本特征参数的统计特征，建立决策规则来进展分类。而不需事先知道类别特征。把各样本的空间分布按其相似性分割或合并成一群集，每一群集代表的地物类别，需经实地调查或与类型的地物加以比较才能确定。常用的非监视分类算法有：K-均值〔K-Means〕算法和迭代自组织数据分析技术〔ISODATA〕。本文使用ISODATA算法对遥感影像进展非监视分类。ISODATA是一种模糊迭代自组织数据分析算法,采用样本修正法调整样本所属类别来完成样本的聚类分析。它是一种非监视分类方法,可以自动的分类合并到达合理的分类数,兼顾了各空间实体之间相关信息,各个区域中特征的非均匀性可包含在初始图像分割处理中,且可以形成原始图像的特征空间,而不会像硬分割那样产生偏倚。ISODATA算法的调试#include"stdafx.h"#include"stdio.h"#include"math.h"#defineN10#defineeps0.00001structPointf{ intsequence; floatx; floaty;};structPointZ{ floatx; floaty;};floatCalDistancef(Pointfx1,Pointfx2){ returnsqrtf((x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y));}floatCalDistanceZ(PointZx1,PointZx2){ returnsqrtf((x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y));}floatCalDistancefZ(Pointfx1,PointZx2){ returnsqrtf((x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y));}intmain(intargc,char*argv[]){ Pointfpts[N]={ {0,0.0,0.0},{1,3.0,8.0},{2,2.0,2.0},{3,1.0,1.0},{4,5.0,3.0}, {5,4.0,8.0},{6,6.0,3.0},{7,5.0,4.0},{8,6.0,4.0},{9,7.0,5.0}, }; inti,j,m; printf("样本集为:\n"); for(i=0;i<N;i++) { printf("X%d(%.1f,%.1f)",i,pts[i].x,pts[i].y); if((i+1)%5==0) { printf("\n"); } } printf("\n"); printf("\n"); intNc=0; printf("pleaseinputNc(0-10):"); scanf("%d",&Nc); intZ[N]; for(i=0;i<Nc;i++) { printf("输入初始第%d聚类中心的序号(0-9):",i); scanf("%d",&Z[i]); }intNj[N];//记录每个类中元素的个数 PointZZArray[N]; PointfSAArray[N][N]; floatDjAv[N]; floatDeltaj[N][2]; floatDeltajmax[N]; intDeltajmaxCor[N]; floatDAv; intNreal=N; intcount=0; floatDij[N*N/2]; intDiji[N]; intDijj[N]; intq=0; intp=0;floatft; intit; intjt; intflag; intss=0; PointZZtp; PointZZArraytp[N]; intNctp; charch; intcur=0; for(i=0;i<N;i++) { Nj[i]=0; } //聚类中心的特征值 for(i=0;i<Nc;i++) { intihere=Z[i]; ZArray[i].x=pts[ihere].x; ZArray[i].y=pts[ihere].y; } intK,ThetaN; floatThetaS,ThetaC; intL,I;Step1: printf("输入预期聚类中心数目K:"); scanf("%d",&K); printf("输入每个聚类域中最少的样本数ThetaN:"); scanf("%d",&ThetaN); printf("输入同一聚类域中样本标准差的最大值:"); scanf("%f",&ThetaS); printf("输入不同聚类域距离最小值:"); scanf("%f",&ThetaC); printf("输入一次可以合并的聚类中心的最多对数:"); scanf("%d",&L); printf("输入最大迭代次数:"); scanf("%d",&I);Step2: for(i=0;i<Nc;i++) { Nj[i]=0; } printf("\n"); printf("这是第%d次归类\n",count+1); for(i=0;i<N;i++)//将模式样本归类 { if(pts[i].sequence==-1)continue;//假设该点的序号为－1则说明它是被剔除的 floatdis=1.0e+10; intxx=0; floatftemp; for(j=0;j<Nc;j++) { ftemp=CalDistancefZ(pts[i],ZArray[j]); if(ftemp<dis||fabs(dis-ftemp)<eps) { xx=j; dis=ftemp; } } SAArray[xx][Nj[xx]].x=pts[i].x; SAArray[xx][Nj[xx]].y=pts[i].y; SAArray[xx][Nj[xx]].sequence=pts[i].sequence; Nj[xx]=Nj[xx]+1; } for(i=0;i<Nc;i++) { printf("第%d个聚类中心是:(%.2f,%.2f)",i,ZArray[i].x,ZArray[i].y); printf("包含的元素有:",i); for(j=0;j<Nj[i];j++) { printf("X%d",SAArray[i][j].sequence); } printf("\n"); //printf("Nj(%d)is%d\n",i,Nj[i]); } count++;Step3: for(j=0;j<Nc;j++)//是否可以去掉一些数据 { if(Nj[j]<ThetaN) { for(i=0;i<Nj[j];i++) { pts[SAArray[j][i].sequence].sequence=-1; } i=j; inttr=j; Nreal-=Nj[j]; while(j<Nc-1) { for(m=0;m<Nj[j+1];m++) { SAArray[j][m].x=SAArray[j+1][m].x; SAArray[j][m].y=SAArray[j+1][m].y; SAArray[j][m].sequence=SAArray[j+1][m].sequence; } j++; } while(i<Nc-1) { Nj[i]=Nj[i+1]; i++; } Nc--; j=tr; } }Step4://修正各聚类中心 for(j=0;j<Nc;j++) { floattemx=0,temy=0; for(i=0;i<Nj[j];i++) { temx+=SAArray[j][i].x; temy+=SAArray[j][i].y; } ZArray[j].x=temx/Nj[j]; ZArray[j].y=temy/Nj[j]; }Step5://计算各聚类域中诸样本与聚类中心的平均距离 floattemp=0.0; for(j=0;j<Nc;j++) { for(i=0;i<Nj[j];i++) { temp+=CalDistancefZ(SAArray[j][i],ZArray[j]); } DjAv[j]=temp/Nj[j]; temp=0.0; }Step6://计算全部模式样本对应聚类中心的总平均距离 DAv=0; for(j=0;j<Nc;j++) { DAv+=Nj[j]*DjAv[j]; } DAv/=Nreal;Step7: if(count>=I)gotoStep14; if(Nc<=K/2)gotoStep8; if((count%2==0)||Nc>=2*K)gotoStep11;Step8://计算各聚类中样本距离标准差 for(j=0;j<Nc;j++) { floattemx=0.0,temy=0.0; for(i=0;i<Nj[j];i++) {temx+=(SAArray[j][i].x-ZArray[j].x)*(SAArray[j][i].x-ZArray[j].x); temy+=(SAArray[j][i].y-ZArray[j].y)*(SAArray[j][i].y-ZArray[j].y); }Deltaj[j][0]=sqrtf(temx/Nj[j]); Deltaj[j][1]=sqrtf(temy/Nj[j]); temx=0.0,temy=0.0; } Step9://求每个标准差向量中的最大分量 for(j=0;j<Nc;j++) { Deltajmax[j]=Deltaj[j][0]>Deltaj[j][1]?Deltaj[j][0]:Deltaj[j][1]; DeltajmaxCor[j]=Deltaj[j][0]>Deltaj[j][1]?0:1; }Step10://分裂判断和计算 for(j=0;j<Nc;j++) { if(Deltajmax[j]>ThetaS) { if((DjAv[j]>DAv&&Nj[j]>2*(ThetaN+1))||Nc<=K