遥感影像监督分类与非监督分类及相关代码实现(31页).doc

《遥感影像监督分类与非监督分类及相关代码实现(31页).doc》由会员分享，可在线阅读，更多相关《遥感影像监督分类与非监督分类及相关代码实现(31页).doc（31页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、-遥感影像监督分类与非监督分类及相关代码实现-第 31 页遥感影像监督分类与非监督分类摘 要：遥感影像的分类方法按照是否有先验类别可分为监督分类和非监督分类，这两种分类方法有着本质的区别但也有存在一定的联系。本文从分类原理和分类方法等不同角度分别介绍了监督分类和非监督分类方法， 并对两种方法的分类结果进行了对比和分析。关键词：遥感；监督分类；非监督分类；ISODATA算法；贝叶斯分类算法。1. 数据来源本文使用的数据是华盛顿广场上空由卫星拍摄的高光谱遥感影像。该幅影像使用的传感器系统覆盖到的可见光到近红外的210个波段。由于米光谱对应的区域上空大气透光性很差，因此将这个区域内的波段从影像中删除

2、，最后得到191个波段。该数据集有1208个扫描行，每行307个像元，容量近似150MB。为了清晰地反映影像中地物的特征，本文从191个波段中选择了3个波段，分别是120、140和160。下面两幅图分别是全波段影像和三波段影像：图1全波段影像 图2 三波段影像2. 遥感影像的监督分类监督分类的原理监督分类 （supervised classification）又称训练场地法，是以建立统计识别函数为理论基础，依据典型样本训练方法进行分类的技术。即根据已知训练区提供的样本，通过选择特征参数，求出特征参数作为决策规则，建立判别函数以对各待分类影像进行的图像分类，是模式识别的一种方法。要求训练区域具有

3、典型性和代表性。判别准则若满足分类精度要求，则此准则成立；反之，需重新建立分类的决策规则，直至满足分类精度要求为止。常用算法有：最大释然分类法、最小距离分类法、马氏距离分类法、平行六面体分类法、K-NN分类法。本文选用最大释然分类法对遥感影像进行监督分类。最大似然判别法. 也称为贝叶斯(Bayes) 分类,是基于图像统计的监督分类法,也是典型的和应用最广的监督分类方法. 它建立在Bayes 准则的基础上,偏重于集群分布的统计特性,分类原理是假定训练样本数据在光谱空间的分布是服从高斯正态分布规律的,做出样本的概率密度等值线,确定分类,然后通过计算标本(像元) 属于各组(类) 的概率,将标本归属于

4、概率最大的一组. 用最大似然法分类,具体分为三步:首先确定各类的训练样本,再根据训练样本计算各类的统计特征值,建立分类判别函数,最后逐点扫描影像各像元,将像元特征向量代入判别函数,求出其属于各类的概率,将待判断像元归属于最大判别函数值的一组。2.2 最大释然分类算法的调试#include#include#include#include#includeconst float PI=3.1415;float train1603,train2603,train3603,train4603; /存各类训练样本 float m13,m23,m33,m43,c133,c233,c333,c433; /各类

5、均值向量及协方差矩阵float p=0.25; /先验概率float test24011; /检验样本 各列分别存放：分类前类别，行号，列号，一波段灰度，二波段灰度，三波段灰度，分类后类别float a1,a2,a3,a4; /公式中常数项；float qq,qq1,qq2,qq3,qq4; /分类精度 float k; /kappa值int q44; /混淆矩阵float a44; /精度指标矩阵 各列分别为生产者精度 漏分误差 /用户精度 错分误差int main()void getdata1(char *filename);void getdata2(char *filename); v

6、oid compute(float train603,float m3,float c33,float *a);float identify(float m3,float c33,float a,float band1,float band2,float band3);int classify(float p1,float p2,float p3,float p4); void testing();void accuracy(float p44);void output(char *filename);void output1(char *filename);getdata1(05traini

7、ng.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);return 0;void

8、getdata1(char *filename)ifstream infile(filename,ios:in|ios:nocreate); /把每类样本的灰度值依次存入对应的数组 if(infile=0) coutopen infile error!endl; exit(1); int i,j;for(i=0;i60;i+)for(j=0;jtrain1ij; for(i=0;i60;i+)for(j=0;jtrain2ij; for(i=0;i60;i+)for(j=0;jtrain3ij; for(i=0;i60;i+)for(j=0;jtrain4ij;void getdata2(ch

9、ar *filename)ifstream infile(filename,ios:in|ios:nocreate); if(infile=0) coutopen infile error!endl; exit(1); int i,j;for(i=0;i240;i+) /检验样本的数量为240for(j=0;jtestij;void compute(float train603,float m3,float c33,float *a)void gauss(int f,float n33); /求n矩阵的逆阵声明 int i,j; for(i=0;i3;i+) mi=0;for(i=0;i3;i

10、+)for(j=0;j60;j+) mi+=trainji; mi=mi/60; float Train360; /转置矩阵for(i=0;i3;i+)for(j=0;j60;j+)Trainij=trainji; int k;for(i=0;i3;i+) /求协方差 for(j=0;j3;j+) cij=0; for(k=0;k60;k+) cij+=(Trainik-mi)*(trainkj-mj); cij=cij/60;*a=c00*c11*c22+c01*c12*c20+c02*c10*c21 -c02*c11*c20-c00*c12*c21-c01*c10*c22; *a=(-0.

11、5)*log(*a)+log(p)-25*log(2*PI); gauss(3,c); float identify(float m3,float c33,float a,float band1,float band2,float band3) / 计算待分样本的d(X)float x3; x0=band1; x1=band2; x2=band3; int i,j;for(i=0;i3;i+)xi-=mi; float d3;for(i=0;i3;i+) di=0; for(j=0;j3;j+)di+=(xi*cji);double q=0;for(i=0;i3;i+) q+=xi*di; q

12、=(-0.5)*q+a; q=exp(q); return q;int classify(float p1,float p2,float p3,float p4) /根据后验概率将像元归为指定类别 float p4; p0=p1; p1=p2; p2=p3; p3=p4; int i,k=0; float max=p0; for(i=1;i4;i+) if(maxpi) max=pi; k=i; k+; return k;void testing()double p1,p2,p3,p4,p5,p6,p7;float kappa(int k44); int i,j;for(i=0;i240;i+

13、) p1=identify(m1,c1,a1,testi3,testi4,testi5); p2=identify(m2,c2,a2,testi3,testi4,testi5); p3=identify(m3,c3,a3,testi3,testi4,testi5); p4=identify(m4,c4,a4,testi3,testi4,testi5);p5=identify(m5,c5,a5,testi3,testi4,testi5);p6=identify(m6,c6,a6,testi3,testi4,testi5);p7=identify(m7,c7,a7,testi3,testi4,te

14、sti5); double dx=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; testi7=p1; testi8=p2; testi9=p3; testi10=p4; testi11=p5; testi12=p6; testi113=p7; testi6=classify(p1,p2,p3,p4,p5,p6,p7); / coutsetiosflags(ios:fixed)setprecision(3)p1 p2 p3 p4 p5 p6 p7endl;fo

15、r(i=0;i7;i+)for(j=0;j7;j+)qij=0; for(i=0;i60;i+) / 混淆矩阵中的各类别组合数量if(testi2=1&testi6=1) q00+;if(testi2=1&testi6=2) q01+;if(testi2=1&testi6=3) q02+;if(testi2=1&testi6=4) q03+; if(testi2=1&testi6=4) q04+; if(testi2=1&testi6=4) q05+; if(testi2=1&testi6=4) q06+; for(i=60;i120;i+)if(testi2=2&testi6=1) q10+

16、;if(testi2=2&testi6=2) q11+;if(testi2=2&testi6=3) q12+; if(testi2=2&testi6=4)q13+; if(testi2=2&testi6=4)q14+; if(testi2=2&testi6=4)q15+; if(testi2=2&testi6=4)q16+;for(i=120;i180;i+)if(testi2=3&testi6=1) q20+;if(testi2=3&testi6=2) q21+;if(testi2=3&testi6=3) q22+;if(testi2=3&testi6=4) q23+; if(testi2=

17、3&testi6=4) q24+; if(testi2=3&testi6=4) q25+; if(testi2=3&testi6=4) q26+;for(i=180;i240;i+)if(testi2=4&testi6=1) q30+;if(testi2=4&testi6=2) q31+;if(testi2=4&testi6=3) q32+;if(testi2=4&testi6=4) q33+; if(testi2=3&testi6=4) q34+; if(testi2=3&testi6=4) q35+; if(testi2=3&testi6=4) q36+; qq1=(float)q00/6

18、0.0; qq2=(float)q11/60.0; qq3=(float)q22/60.0;qq4=(float)q33/60.0;qq5=(float)q43/60.0;qq6=(float)q53/60.0;qq7=(float)q63/60.0; qq=(qq1+qq2+qq3+qq4+ qq5+qq6+qq7)/4;qq=qq*100; k=kappa(q);/cout1 代表 房顶 ， 2 代表 街道 ， 3 代表 小路 ，4 代表 草地 ，5 代表 树木，6 代表 水体，7 代表 阴影endl;/cout正确分类百分比分别为：; / coutqq1 qq2 qq3 qq4 qq5

19、qq6 qq7endl;/ cout总分类精度:setprecision(2)qq%endl;/coutKappa值为: setprecision(2)k%endl;void output(char *filename) ofstream outfile(filename,ios:out); if(outfile=0) coutopen outfile error!endl; exit(1); outfile1 代表 房顶 ， 2 代表 街道 ， 3 代表 小路 ，4 代表 草地 ，5 代表 树木，6 代表 水体，7 代表 阴影endl; outfile正确分类百分比分别为：; outfile

20、setprecision(3)qq1 qq2 qq3 qq4endl; outfile总分类精度:setprecision(4)qq%endl; outfilerow col before after proba1 proba2 proba3 proba4 endl; int i; for(i=0;i240;i+) outfilesetiosflags(ios:fixed)setprecision(0)testi0 testi1 testi2 testi6 ; outfilesetiosflags(ios:fixed)setprecision(3)testi7 testi8 testi9 te

21、sti10=0;k-) p=n00; for(i=1;ik)?q:-q)/p; for(j=1;j=i;j+) ni-1j-1=nij+q*hj; nf-1f-1=1/p; for(i=1;if;i+) nf-1i-1=hi; for(i=1;if;i+) for(j=0;ji;j+) nji=nij; float kappa(int k44)int i,j;float K;int c4,d4; int N=0, M=0, D=0; for(i=0;i4;i+)ci=0;di=0;for(i=0;i4;i+)for(j=0;j4;j+) N+=kij;if(i=j) D=D+kij;ci=ki

22、j+ci;dj=kij+dj;for(i=0;i4;i+)M=ci*di+M;K=float(N*D-M)/(N*N-M);return K;void accuracy(float p44)int a4,b4;float a14,a24,b14,b24;int i,j;for(i=0;i4;i+)ai=0; bi=0;a1i=0; b1i=0;a2i=0; b2i=0;for(i=0;i4;i+)for(j=0;j4;j+)pij=0;ai+=qij;bi+=qji;for(i=0;i4;i+)a1i=(float)qii/ai;a2i=1-a1i;b1i=(float)qii/bi;b2i=

23、1-b1i; for(i=0;i4;i+)pi0=a1i*100; /生产者精度pi1=a2i*100; /漏分误差pi2=b1i*100; /用户精度pi3=b2i*100; /错分误差 void output1(char *filename)ofstream outfile(filename,ios:out); if(outfile=0) coutopen outfile error!endl; exit(1); outfile Lakes Built-up Land Forest Land Cropland endl; int i,j; for(i=0;i4;i+) switch(i)c

24、ase 0:outfileLakes ;break;case 1: outfileBuilt-up Land ;break;case 2:outfileForest Land ;break;case 3:outfileCropland ;break; for(j=0;j4;j+) outfilesetiosflags(ios:left)setw(7)qij ;if(j!=0&j%3=0) outfileendl; outfileendl; outfileendl; outfileendl; outfile setw(11)Producers setw(11)Omission Users Com

25、mision endl;for(i=0;i4;i+) switch(i)case 0:outfileLakes ;break;case 1: outfileBuilt-up Land ;break;case 2:outfileForest Land ;break;case 3: outfileCropland ;break;for(j=0;j4;j+) outfile setiosflags(ios:fixed)setprecision(2)setw(4)aij% ; if(j!=0&j%3=0) outfileendl; outfileendl;outfileendl;outfile 总分类

26、精度:setprecision(2)qq%endl;outfile Kappa系数 :setprecision(3)kendl;2.3 最大释然分类算法的执行 在贝叶斯分类器中执行以上代码，得到的监督分类图如下所示： 图三 监督分类图 利用最大释然分类算法对遥感影像进行监督分类，最后将图像上的地类分为7类，分别是：房顶，街道，小路，草地，树木，水体和阴影。从处理的图像来看，分类结果较为很理想，基本上可以辨认出这7类地物。 3. 非监督分类3.1 非监督分类的原理 非监督分类是以不同影像地物在特征空间中类别特征的差别为依据的一种无先验（已知）类别标准的图像分类，是以集群为理论基础，通过计算机对图

27、像进行集聚统计分析的方法。根据待分类样本特征参数的统计特征，建立决策规则来进行分类。而不需事先知道类别特征。把各样本的空间分布按其相似性分割或合并成一群集，每一群集代表的地物类别，需经实地调查或与已知类型的地物加以比较才能确定。常用的非监督分类算法有：K-均值（K-Means）算法和迭代自组织数据分析技术（ISODATA）。本文使用ISODATA算法对遥感影像进行非监督分类。 ISODATA 是一种模糊迭代自组织数据分析算法, 采用样本修正法调整样本所属类别来完成样本的聚类分析。它是一种非监督分类方法,可以自动的分类合并达到合理的分类数, 兼顾了各空间实体之间相关信息, 各个区域中特征的非均匀

28、性可包含在初始图像分割处理中, 且可以形成原始图像的特征空间, 而不会像硬分割那样产生偏倚。3.2 ISODATA算法的调试#include #include #include #define N 10#define epsstruct Pointfint sequence;float x;float y;struct PointZfloat x;float y;float CalDistancef(Pointf x1,Pointf x2)return sqrtf(x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y);float CalDistanceZ

29、(PointZ x1,PointZ x2)return sqrtf(x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y);float CalDistancefZ(Pointf x1,PointZ x2)return sqrtf(x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y);int main(int argc, char* argv)Pointf ptsN=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

30、.0,8,6.0,4.0,9,7.0,5.0,int i,j,m;printf(样本集为:n);for(i=0;iN;i+)printf(X%d(%.1f,%.1f) ,i,ptsi.x,ptsi.y);if(i+1)%5=0)printf(n);printf(n);printf(n);int Nc=0;printf(please input Nc(0-10): );scanf(%d,&Nc);int ZN;for(i=0;iNc;i+)printf(输入初始第%d聚类中心的序号(0-9):,i);scanf(%d,&Zi); int NjN; /记录每个类中元素的个数PointZ ZArra

31、yN;Pointf SAArrayNN;float DjAvN;float DeltajN2;float DeltajmaxN;int DeltajmaxCorN;float DAv;int Nreal=N;int count=0;float DijN*N/2;int DijiN;int DijjN;int q=0;int p=0; float ft;int it;int jt;int flag;int ss=0;PointZ Ztp;PointZ ZArraytpN;int Nctp;char ch;int cur=0;for(i=0;iN;i+)Nji=0; /聚类中心的特征值for(i=

32、0;iNc;i+)int ihere=Zi;ZArrayi.x=ptsihere.x;ZArrayi.y=ptsihere.y;int K,ThetaN;float ThetaS,ThetaC;int L,I;Step1:printf(输入预期聚类中心数目 K :);scanf(%d,&K);printf(输入每个聚类域中最少的样本数ThetaN: );scanf(%d,&ThetaN);printf(输入同一聚类域中样本标准差的最大值: );scanf(%f,&ThetaS);printf(输入不同聚类域距离最小值: );scanf(%f,&ThetaC);printf(输入一次可以合并的聚

33、类中心的最多对数: );scanf(%d,&L);printf(输入最大迭代次数: );scanf(%d,&I);Step2:for(i=0;iNc;i+)Nji=0;printf(n);printf(这是第%d次归类n,count+1);for(i=0;iN;i+) /将模式样本归类if(ptsi.sequence=-1)continue; /若该点的序号为1则说明它是被剔除的float dis=1.0e+10;int xx=0;float ftemp;for(j=0;jNc;j+) ftemp=CalDistancefZ(ptsi,ZArrayj);if(ftempdis|fabs(dis-ftemp)eps)xx=j;dis=ftemp;SAArrayxxNjxx.x=ptsi.x;SAArrayxxNjxx.y=ptsi.y;SAArrayxxNjxx.sequence=ptsi.sequence;Njxx=Njxx+1; for(i=0;iNc;i+) printf(第%d个聚类中心是:(%.2f,%.2f) ,i,ZArrayi.x,ZArrayi.y); printf(包含的元素有:,i); for(j=0;jNji;j+) printf( X%d ,SAArrayij.sequence); printf(n);