遙感影像監(jiān)督分類與非監(jiān)督分類及相關(guān)代碼實(shí)現(xiàn)

1、遙感影像監(jiān)督分類與非監(jiān)督分類摘 要：遙感影像的分類方法按照是否有先驗(yàn)類別可分為監(jiān)督分類和非監(jiān)督分類，這兩種分類方法有著本質(zhì)的區(qū)別但也有存在一定的聯(lián)系。本文從分類原理和分類方法等不同角度分別介紹了監(jiān)督分類和非監(jiān)督分類方法， 并對兩種方法的分類結(jié)果進(jìn)行了對比和分析。關(guān)鍵詞：遙感；監(jiān)督分類；非監(jiān)督分類；ISODATA算法；貝葉斯分類算法。1. 數(shù)據(jù)來源本文使用的數(shù)據(jù)是華盛頓廣場上空由衛(wèi)星拍攝的高光譜遙感影像。該幅影像使用的傳感器系統(tǒng)覆蓋0.4到2.4m的可見光到近紅外的210個(gè)波段。由于米光譜對應(yīng)的區(qū)域上空大氣透光性很差，因此將這個(gè)區(qū)域內(nèi)的波段從影像中刪除，最后得到191個(gè)波段。該數(shù)據(jù)集有1208個(gè)

2、掃描行，每行307個(gè)像元，容量近似150MB。為了清晰地反映影像中地物的特征，本文從191個(gè)波段中選擇了3個(gè)波段，分別是120、140和160。下面兩幅圖分別是全波段影像和三波段影像： 圖1全波段影像 圖2 三波段影像2. 遙感影像的監(jiān)督分類2.1監(jiān)督分類的原理監(jiān)督分類 （supervised classification）又稱訓(xùn)練場地法，是以建立統(tǒng)計(jì)識別函數(shù)為理論基礎(chǔ)，依據(jù)典型樣本訓(xùn)練方法進(jìn)行分類的技術(shù)。即根據(jù)已知訓(xùn)練區(qū)提供的樣本，通過選擇特征參數(shù)，求出特征參數(shù)作為決策規(guī)則，建立判別函數(shù)以對各待分類影像進(jìn)行的圖像分類，是模式識別的一種方法。要求訓(xùn)練區(qū)域具有典型性和代表性。判別準(zhǔn)則若滿足分類精

3、度要求，則此準(zhǔn)則成立；反之，需重新建立分類的決策規(guī)則，直至滿足分類精度要求為止。常用算法有：最大釋然分類法、最小距離分類法、馬氏距離分類法、平行六面體分類法、K-NN分類法。本文選用最大釋然分類法對遙感影像進(jìn)行監(jiān)督分類。最大似然判別法. 也稱為貝葉斯(Bayes) 分類,是基于圖像統(tǒng)計(jì)的監(jiān)督分類法,也是典型的和應(yīng)用最廣的監(jiān)督分類方法. 它建立在Bayes 準(zhǔn)則的基礎(chǔ)上,偏重于集群分布的統(tǒng)計(jì)特性,分類原理是假定訓(xùn)練樣本數(shù)據(jù)在光譜空間的分布是服從高斯正態(tài)分布規(guī)律的,做出樣本的概率密度等值線,確定分類,然后通過計(jì)算標(biāo)本(像元) 屬于各組(類) 的概率,將標(biāo)本歸屬于概率最大的一組. 用最大似然法分類,

4、具體分為三步:首先確定各類的訓(xùn)練樣本,再根據(jù)訓(xùn)練樣本計(jì)算各類的統(tǒng)計(jì)特征值,建立分類判別函數(shù),最后逐點(diǎn)掃描影像各像元,將像元特征向量代入判別函數(shù),求出其屬于各類的概率,將待判斷像元?dú)w屬于最大判別函數(shù)值的一組。2.2 最大釋然分類算法的調(diào)試#include#include#include#include#includeconst float PI=3.1415;float train1603,train2603,train3603,train4603; /存各類訓(xùn)練樣本 float m13,m23,m33,m43,c133,c233,c333,c433; /各類均值向量及協(xié)方差矩陣float p=

5、0.25; /先驗(yàn)概率float test24011; /檢驗(yàn)樣本 各列分別存放：分類前類別，行號，列號，一波段灰度，二波段灰度，三波段灰度，分類后類別float a1,a2,a3,a4; /公式中常數(shù)項(xiàng)；float qq,qq1,qq2,qq3,qq4; /分類精度 float k; /kappa值int q44; /混淆矩陣float a44; /精度指標(biāo)矩陣 各列分別為生產(chǎn)者精度 漏分誤差 /用戶精度 錯(cuò)分誤差int main() void getdata1(char *filename);void getdata2(char *filename); void compute(float

6、 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(05training.txt); /存放的是訓(xùn)練樣

7、本三個(gè)波段的灰度值getdata2(05testing.txt); /檢驗(yàn)樣本行號 列號 類別代碼 三個(gè)波段的灰度值 compute(train1,m1,c1,&a1); /計(jì)算各類均值向量及協(xié)方差矩陣compute(train2,m2,c2,&a2);compute(train3,m3,c3,&a3);compute(train4,m4,c4,&a4); testing(); /對檢驗(yàn)樣本分類accuracy(a);output(05test-result.txt);output1(05confusion-matrix.txt);return 0;void getdata1(char *fi

8、lename) ifstream infile(filename,ios:in|ios:nocreate); /把每類樣本的灰度值依次存入對應(yīng)的數(shù)組 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;infile.close();void getdata2(cha

9、r *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+) /檢驗(yàn)樣本的數(shù)量為240for(j=0;jtestij; infile.close();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+)

10、mi=0;for(i=0;i3;i+) for(j=0;j60;j+) mi+=trainji; mi=mi/60; float Train360; /轉(zhuǎn)置矩陣for(i=0;i3;i+)for(j=0;j60;j+)Trainij=trainji; int k;for(i=0;i3;i+) /求協(xié)方差 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-c

11、01*c10*c22; *a=(-0.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) / 計(jì)算待分樣本的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

12、;i3;i+) q+=xi*di; q=(-0.5)*q+a; q=exp(q); return q;int classify(float p1,float p2,float p3,float p4) /根據(jù)后驗(yàn)概率將像元?dú)w為指定類別 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

13、 i,j;for(i=0;i240;i+) 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

14、,a7,testi3,testi4,testi5); 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 p

15、3 p4 p5 p6 p7endl; for(i=0;i7;i+)for(j=0;j7;j+)qij=0; for(i=0;i60;i+) / 混淆矩陣中的各類別組合數(shù)量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(te

16、sti2=2&testi6=1) q10+;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&testi

17、6=4) q23+; if(testi2=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) q

18、36+; qq1=(float)q00/60.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正確分類百分比分別為：; / cou

19、tqq1 qq2 qq3 qq4 qq5 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; outfil

20、e正確分類百分比分別為：; outfilesetprecision(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)te

21、sti7 testi8 testi9 testi10=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;

22、if(i=j) D=D+kij;ci=kij+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-a

23、1i;b1i=(float)qii/bi;b2i=1-b1i; for(i=0;i4;i+)pi0=a1i*100; /生產(chǎn)者精度pi1=a2i*100; /漏分誤差pi2=b1i*100; /用戶精度pi3=b2i*100; /錯(cuò)分誤差 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;

24、 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+) outfilesetiosflags(ios:left)setw(7)qij ;if(j!=0&j%3=0) outfileendl; outfileendl; outfileendl; outfileendl; outfile setw(11)Producers

25、 setw(11)Omission Users Commision 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; outfil

26、eendl;outfileendl;outfile 總分類精度:setprecision(2)qq%endl;outfile Kappa系數(shù) :setprecision(3)kendl; outfile.close();2.3 最大釋然分類算法的執(zhí)行 在貝葉斯分類器中執(zhí)行以上代碼，得到的監(jiān)督分類圖如下所示： 圖三 監(jiān)督分類圖 利用最大釋然分類算法對遙感影像進(jìn)行監(jiān)督分類，最后將圖像上的地類分為7類，分別是：房頂，街道，小路，草地，樹木，水體和陰影。從處理的圖像來看，分類結(jié)果較為很理想，基本上可以辨認(rèn)出這7類地物。 3. 非監(jiān)督分類3.1 非監(jiān)督分類的原理 非監(jiān)督分類是以不同影像地物在特征空間中類

27、別特征的差別為依據(jù)的一種無先驗(yàn)（已知）類別標(biāo)準(zhǔn)的圖像分類，是以集群為理論基礎(chǔ)，通過計(jì)算機(jī)對圖像進(jìn)行集聚統(tǒng)計(jì)分析的方法。根據(jù)待分類樣本特征參數(shù)的統(tǒng)計(jì)特征，建立決策規(guī)則來進(jìn)行分類。而不需事先知道類別特征。把各樣本的空間分布按其相似性分割或合并成一群集，每一群集代表的地物類別，需經(jīng)實(shí)地調(diào)查或與已知類型的地物加以比較才能確定。常用的非監(jiān)督分類算法有：K-均值（K-Means）算法和迭代自組織數(shù)據(jù)分析技術(shù)（ISODATA）。本文使用ISODATA算法對遙感影像進(jìn)行非監(jiān)督分類。 ISODATA 是一種模糊迭代自組織數(shù)據(jù)分析算法, 采用樣本修正法調(diào)整樣本所屬類別來完成樣本的聚類分析。它是一種非監(jiān)督分類方法,

28、可以自動的分類合并達(dá)到合理的分類數(shù), 兼顧了各空間實(shí)體之間相關(guān)信息, 各個(gè)區(qū)域中特征的非均勻性可包含在初始圖像分割處理中, 且可以形成原始圖像的特征空間, 而不會像硬分割那樣產(chǎn)生偏倚。3.2 ISODATA算法的調(diào)試#include stdafx.h#include stdio.h#include math.h#define N 10#define eps 0.00001struct Pointfint sequence;float x;float y;struct PointZfloat x;float y;float CalDistancef(Pointf x1,Pointf x2)ret

29、urn sqrtf(x1.x-x2.x)*(x1.x-x2.x)+(x1.y-x2.y)*(x1.y-x2.y);float CalDistanceZ(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,

30、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,;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+)pri

31、ntf(輸入初始第%d聚類中心的序號(0-9):,i);scanf(%d,&Zi); int NjN; /記錄每個(gè)類中元素的個(gè)數(shù)PointZ ZArrayN;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;Poin

32、tZ ZArraytpN;int Nctp;char ch;int cur=0;for(i=0;iN;i+)Nji=0; /聚類中心的特征值for(i=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(輸入預(yù)期聚類中心數(shù)目 K :);scanf(%d,&K);printf(輸入每個(gè)聚類域中最少的樣本數(shù)ThetaN: );scanf(%d,&ThetaN);printf(輸入同一聚類域中樣本標(biāo)準(zhǔn)差的最大值: );

33、scanf(%f,&ThetaS);printf(輸入不同聚類域距離最小值: );scanf(%f,&ThetaC);printf(輸入一次可以合并的聚類中心的最多對數(shù): );scanf(%d,&L);printf(輸入最大迭代次數(shù): );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; /若該點(diǎn)的序號為1則說明它是被剔除的float dis=1.0e+10;int xx=0;float

34、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個(gè)聚類中心是:(%.2f,%.2f) ,i,ZArrayi.x,ZArrayi.y); printf(包含的元素有:,i); for(j=0;jNji;

35、j+) printf( X%d ,SAArrayij.sequence); printf(n); / printf(Nj(%d) is %dn,i,Nji); count+;Step3:for(j=0;jNc;j+) /是否可以去掉一些數(shù)據(jù)if(NjjThetaN)for(i=0;iNjj;i+)ptsSAArrayji.sequence.sequence=-1;i=j;int tr=j;Nreal-=Njj;while(jNc-1)for(m=0;mNjj+1;m+)SAArrayjm.x=SAArrayj+1m.x;SAArrayjm.y=SAArrayj+1m.y;SAArrayjm.s

36、equence=SAArrayj+1m.sequence; j+;while(iNc-1)Nji=Nji+1;i+;Nc-;j=tr;Step4: /修正各聚類中心for(j=0;jNc;j+)float temx=0,temy=0;for(i=0;iNjj;i+)temx+=SAArrayji.x;temy+=SAArrayji.y;ZArrayj.x=temx/Njj;ZArrayj.y=temy/Njj;Step5:/計(jì)算各聚類域中諸樣本與聚類中心的平均距離float temp=0.0;for(j=0;jNc;j+)for(i=0;iNjj;i+)temp+=CalDistancefZ(

37、SAArrayji,ZArrayj);DjAvj=temp/Njj;temp=0.0;Step6:/計(jì)算全部模式樣本對應(yīng)聚類中心的總平均距離DAv=0;for(j=0;j=I) goto Step14;if(Nc=2*K) goto Step11;Step8:/計(jì)算各聚類中樣本距離標(biāo)準(zhǔn)差for(j=0;jNc;j+)float temx=0.0,temy=0.0;for(i=0;iNjj;i+) temx+=(SAArrayji.x-ZArrayj.x)*(SAArrayji.x-ZArrayj.x); temy+=(SAArrayji.y-ZArrayj.y)*(SAArrayji.y-ZA

38、rrayj.y); Deltajj0=sqrtf(temx/Njj);Deltajj1=sqrtf(temy/Njj);temx=0.0,temy=0.0;Step9:/求每個(gè)標(biāo)準(zhǔn)差向量中的最大分量 for(j=0;jDeltajj1?Deltajj0:Deltajj1;DeltajmaxCorj=Deltajj0Deltajj1?0:1;Step10:/分裂判斷和計(jì)算for(j=0;jThetaS)if(DjAvjDAv&Njj2*(ThetaN+1)|Nc=K/2)float Garma=0.5;PointZ Zj1,Zj2;if(DeltajmaxCorj=0)Zj1.x=ZArrayj.x+Deltajmaxj*Garma;Zj1.y=ZArrayj.y;Zj2