高光譜圖像數(shù)學建模期末報告

模型建立與求解1.1高光譜圖像高光譜圖像在某種程序上比一般遙感圖像分類更精準，是因為很多礦物的光譜特征只能利用高光譜數(shù)據(jù)才能被探測到。高光譜數(shù)據(jù)的光譜分辨率比寬波段遙感高數(shù)十倍，在寬波段上無法反映出這些光譜特征，但在高光譜影像上很容易識別。1.2高光譜遙感不同于傳統(tǒng)遙感的新特點：⑴波段多：可以為每個像元提供幾十、數(shù)百甚至上千個波段；⑵光譜范圍窄：波段范圍一般小于10nm；⑶波段連續(xù)：有些傳感器可以在350~2500nm的太陽光內(nèi)提供幾乎連續(xù)的地物光譜；⑷數(shù)據(jù)量大：隨著波段數(shù)年的增加，數(shù)據(jù)量成指數(shù)增加；由目中提供的數(shù)據(jù)可以看出，光譜的波段較連續(xù)，所以此模型我采用高光譜圖像的方法進行求解。在光譜遙感圖像信息對地物進行精準的分類方法有很多，如叉相關(guān)光譜匹配，光譜吸收特征匹配，二值編碼匹配等。在這些分類中，光譜角度匹配分類是應(yīng)用最廣泛且最精確的分類方法之一。

接下來，我就使用光譜角度法進行計算。2.1光譜角度填圖光譜角度填圖方法的具體過程是將光譜數(shù)據(jù)視為多維空間的矢量 ，禾U用解析方法計算像元光譜與光譜數(shù)據(jù)庫光譜或像元訓練光譜之間矢量的夾角，根據(jù)夾角的大小來確定光譜間的相似程度 ，以達到識別地物的目的。由于與矢量的模無關(guān)，光譜角度對增益是不敏感的，從而能夠減弱因地形和照度引起的增益變化，因此訓練光譜可以直接和像元光譜相比較。光譜角度圖像的結(jié)果，是對像元賦予0°—180°的角度值,給出像元光譜與訓練光譜之間相似性的定量估計。2.2光譜角度填圖原理設(shè)有兩個n個波段的光譜向量，一個為光譜庫光譜或一個像元訓練光譜為T=(t1,t2,…tn),—個為實測地物光譜為R=(r1,r2,…rn),T、R不是零向量。它們的廣義夾角0定義為0=cos-1(TR)'TR/ntr即0二COS-1( 日t心)2iAL

ni=12iAL

ni=12irni=1以第1列數(shù)據(jù)為參考光譜，則將第2列到第28列中的光譜的幾何平均向量為類中心。設(shè)已知這 27列光譜中有27個點R2,R3,……,R28，則類中心為M=亦醫(yī)Ri2.3具體光譜角度分類步驟：⑴從光譜庫數(shù)據(jù)中選擇最感興趣的“最終成分光譜”，或選擇一個像元訓練區(qū)，并從中算出這個區(qū)域的幾何平均向量為類中心，形成“最終成分光譜”；⑵對“最終成分光譜”做重采樣，使“最終成分光譜”與實測地物的光譜分辨率和波長范圍一致。⑶計算最終成分光譜與圖像像元光譜兩個光譜向量 T=(t1,t2,…tn),R=(r1,r2,…rn)之間的角度0，以評價此兩個光譜向量相似性。9值域為(0,n/2),當0=0時，表示兩個光譜完全相似；當a=n/2時，則兩個光譜完全不同。3.最小二乘法設(shè)o(x),i(x),..., n(x)為n+1個線性無關(guān)(與向量的線性無關(guān)定義類似)的連續(xù)函數(shù)，門為o(x),i(x),..., n(x)所組成的n+1維線性n空間，即由其所有線性組合aJk(x),ak?R(k=0,1,...,n)構(gòu)成的k衛(wèi)集合，記作=spar<o(x),i(x),..., 「n(x)｝，任取p(x)?，則np(x) akk(x)，它是關(guān)于的線性函數(shù)。k=0對已知數(shù)據(jù)點(xi,yi)(i=0,1,..., m,在門中求p(x),使得m mnI [P(Xi)-yi]2 ['ak：k(Xi)-yi]2二mini=0 i二0k=0mn一2、('ak：k(Xi)-yi)](Xi)=0,j-0,1,...,n■aj i=0k=0m(jk)二'j(Xj)k(xi),j,k=0,1,...n,i=Q則得到如下的求解方程組

-1-1-他,％）(%,%)（%,粋）a0Wo,y)(%%)(%,%)（%,件）ai?,y)=??…(4)(%,%)（粋用1）（粋，粋）an%,y)?I.一把數(shù)據(jù)帶入方程組得到唯一解a=a0,a1,…，an則p（x）八akk（x）4.1帶入數(shù)據(jù)求解將第一列波得到得圖像作為已知曲線，利用matlab實現(xiàn)最小二乘,得到未知的27列地物發(fā)射信息的擬合曲線，再與已知曲線進行比較,得到他們之間的夾角0,由B的值判斷27列信息分別是哪種地物。得到的地物光譜曲線如下圖：第一類憫卿t-蛇紋石謁始光譜曲線御卩詭?蛇紋昏-包絡(luò)線去除后光譜陸線第二類

1chlorite-^扼石-眞皓光譜曲線chloritfr^混石“包結(jié)線上社忑無常曲域第三類epidote-^簾石"原箱北譜曲塢500 1000 1500 2000 2500epidote--^1chlorite-^扼石-眞皓光譜曲線chloritfr^混石“包結(jié)線上社忑無常曲域第三類epidote-^簾石"原箱北譜曲塢500 1000 1500 2000 2500epidote--^帝W--包略線1?除后北譜曲線3000第四類muscoM—&云母"原jtiSr曲寅1r021 i i i i i i0 500 10D0 1500 2000 25D0 3000musccvi-白云母-包終線去隱后允譜曲線第五類serpent-gt氨石-M姑北譜曲找0.4serpent-^紋石-包皓線去除后光詣曲線4.2最終結(jié)果第一類：5、6、7、11第二類：8 9、 13第三類：2、3、4第四類：10、12。14-26第五類：27、28、結(jié)果分析若輸入的是一組數(shù)據(jù)，將輸入的數(shù)據(jù)帶入這五條擬合曲線得出f,跟y做差方分析，若越小，則表示此地物跟相應(yīng)的類對應(yīng)，既可以判斷該地物是哪種。若輸入的是多組數(shù)據(jù)，則跟各類的主因子求最短距離，每組數(shù)據(jù)的最短距離相差不大的話，這幾組數(shù)據(jù)是一類，即這幾種地物是同一種。參考文獻[1]蘇紅軍,高光譜遙感圖像光譜特征選擇與特征提取的研究,遙感應(yīng)用與技術(shù)，2006.[2]張良培、張立福，高光譜遙感,武漢大學出版社，2005年10月[3]王志剛、朱振海、王紅梅、劉慶生，光譜角度填圖方法及其在巖性識別中的應(yīng)用,遙感學報1999年2月[4]王潤生,郭小方,王天興,等,成像光譜方法技術(shù)開發(fā)應(yīng)用研究,中國國土資源航空物探遙感中心,1999clearallloaddata;loadn1;loadn2;loadn3;loadn4;loadn5;A=alteredmineral_original;X=A(:,2:28);x=X';%聚類分析%j=zscore(x);y=pdist(j);m=linkage(y);t=cluster(m,5);stem(t,'DisplayName','t','YDataSource','t');figure(gcf)a=A';Y=A(:,1);%主成分分析%rn1=corrcoef(n1);rn2=corrcoef(n2);rn3=corrcoef(n3);rn4=corrcoef(n4);rn5=corrcoef(n5);[vn1,dn1]=eig(rn1);[vn2,dn2]=eig(rn2);[vn3,dn3]=eig(rn3);[vn4,dn4]=eig(rn4);[vn1,dn5]=eig(rn5);%求貢獻率%gz1=sum(dn1(:));gz2=sum(dn2(:));gz3=sum(dn3(:));gz4=sum(dn4(:));gz5=sum(dn5(:));gz11=dn1(1,1)/gz1;gz12=dn1(2,2)/gz1;gz13=dn1(3,3)/gz1;gz14=dn1(4,4)/gz1;gz21=dn2(1,1)/gz2;gz22=dn2(2,2)/gz2;gz23=dn2(3,3)/gz2;gz31=dn3(1,1)/gz3;gz32=dn3(2,2)/gz3;gz33=dn3(3,3)/gz3;gz51=dn5(1,1)/gz5;gz52=dn5(2,2)/gz5;gz41=dn4(1,1)/gz4;gz42=dn4(2,2)/gz4;gz43=dn4(3,3)/gz4;gz44=dn4(4,4)/gz4;gz45=dn4(5,5)/gz4;gz46=dn4(6,6)/gz4;gz47=dn4(7,7)/gz4;gz48=dn4(8,8)/gz4;gz49=dn4(9,9)/gz4;gz410=dn4(10,10)/gz4;gz411=dn4(11,11)/gz4;gz412=dn4(12,12)/gz4;gz413=dn4(13,13)/gz4;gz414=dn4(14,14)/gz4;gz415=dn4(15,15)/gz4;g1=[gz11gz12gz13gz14];stem(g1,'DisplayName','g1','YDataSource','g1');figure(gcf);Xlabel('因子');Ylabel('貢獻率');title(' 第1類');g2=[gz21gz22gz23];stem(g1,'DisplayName','g1','YDataSource','g1');figure(gcf);Xlabel('因子');Ylabel('貢獻率');title(' 第2類');g3=[gz31gz32gz33];stem(g1,'DisplayName','g1','YDataSource','g1');figure(gcf);Xlabel('因子');Ylabel('貢獻率');title(' 第3類');g4=[gz41gz42gz43gz44gz45gz46gz47gz48gz49gz410gz411gz412gz413gz414gz415];stem(g1,'DisplayName','g1','YDataSource','g1');figure(gcf);Xlabel('因子');Ylabel('貢獻率');title(' 第4類');g5=[gz51gz52];stem(g1,'DisplayName','g1','YDataSource','g1');figure(gcf);Xlabel('因子');Ylabel('貢獻率');title(' 第5類');%主因子擬合%x1=n1(:,1);x2=n2(:,1);x3=n3(:,3);x4=n4(:,1);x5=n5(:,1);p1=polyfit(x1,Y,2);p2=polyfit(x2,Y,2);p3=polyfit(x3,Y,2);p4=polyfit(x4,Y,2);p5=polyfit(x5,Y,2);y1=polyval(p1,Y);y2=polyval(p2,Y);y3=polyval(p3,Y);y4=polyval(p4,Y);y5=polyval(p5,Y);holdon;plot(x1,y1,'-');title(' 第1類擬合圖');figure;holdon;plot(x2,y2,'-');title（'第2類擬合圖'）;figure;holdon;plot(x3,y3,'-');title(' 第3類擬合圖');figure;holdon;plot(x4,y4,'-');title(' 第4類擬合圖');figure;holdon;plot(x5,y5,'-');title(' 第5類擬合圖');f=inline('a(1)*x+a(2)*x.A2.*exp(-a(3)*x)+a(4)','a','x');[xx1,res1]=lsqcurvefit(f,ones(1,4),x1,y1);disp('第1類的擬合函數(shù)是：');disp(poly2str(xx1,'x'));[xx2,res2]=lsqcurvefit(f,ones(1,4),x2,y2);disp('第2類的擬合函數(shù)是：');disp(poly2str(xx2,'x