遙感原理實習報告

1、遙感原理課程設計專業(yè)：遙感科學與技術班級：學號：姓名：指導老師：一、原理介紹1.1 K-均值算法算法主要思想是先在需要分類的數(shù)據中尋找K組數(shù)據作為初始聚類中心，然后計算所有數(shù)據到各聚類中心的距離，將數(shù)據歸入與其距離最近的聚 類中心，之后再對這 K個聚類的數(shù)據計算均值，作為新的聚類中心，重復 以上步驟，直到新的聚類中心與上一次的聚類中心值相等時結束算法。算法 流程圖如下：否1.2 ISODATA （迭代自組織分析）通過設定初始參數(shù)而引入人機對話環(huán)節(jié)，并使用歸并與分裂的機制，當 某兩類聚類中心距離小于某一閾值時，將它們合并為一類，當某類標準差大 于某一閾值或其樣本數(shù)目超過某一閾值時，將其分為兩類。

2、在某類樣本數(shù)目 少于某閾值時，需將其取消。如此，根據初始聚類中心和設定的類別數(shù)目等 參數(shù)迭代，最終得到一個比較理想的分類結果。算法流程圖如下：1.3 IHS變換影像融合1.3.1 IHS變換影像融合的一般步驟為：1）將RGB圖像變換到IHS空間，得到I、H、S三個分量；2）對變換后的I分量與全色影像進行直方圖匹配，獲得修改后的I分量;3）最后用修改后的I分量結合H、S分量進行HIS逆變換，得到融合影 像。1.3.2 正變換：RGB to IHS1I （R G B）33RGBmin( R,G,B)=arccos(R_G) +(R_B)/2J(R_G)2 +(R_B)(G _B)1/21.3.3

3、逆變換：IHS to RGB當H在0,120之間B = I (1 -S)- ScosH 1R = I 1 +cos(60 - H )G =31 -(B R)當H在120,240之間R = I (1 -S)G=i 1Scos(H-120)IL cos(180 - H )B =31 -(R G)當H在240,360之間G =1(1 -S)一 Sco sH 2 40B =1 1-ILc o $3(0 0 H )R =31 -(G B)二、算法設計1.4 K均值算法1）讀取圖像，獲取圖像寬度 nWidth、高度nHeight、波段數(shù)nBand、數(shù)據首地 址pData、顯示數(shù)據指針pByte；2）將圖像

4、像素數(shù)據復制一份到 DataCopy中，之后用DataCopy處理，處理結果 重新賦值到pData中，刷新即可得到分類結果圖，從而避免了新建圖像這一 復雜過程。3）調用對話框，獲取分類類別數(shù)nClass,并定義中間變量參數(shù)：int *pixelClass=new intnWidth*nHeight;/ 判斷每個像素類別 double *mea nOld=new double nBand*n Class;/ 日的聚類中心 double *mea nNew=n ew double nBa nd* nClass;/新的聚類中心 double *gray=new double nBa nd; 某個像素

5、的各波段灰度值 double *dis=new double nClass;/某個像素到各個聚類中心的距離 in t *cou nt=n ew in t nClass;/ 每個類別像素數(shù)量4)5)6)7)8)1.57)8)9)10)11)12)13)14)確定初始聚類中心，選擇遙感影像前nClass個像素作為初始聚類中心Z Z2JJZ：(k=1,2,3,.,nClass；計算所有像素點到各個聚類中心的馬氏距離, 類，將所有像素劃分到對應類別；按照馬氏距離最小原則進行分計算各聚類中心的新向量值1k 1 丄 X k 1-Z j _ n j x Sj(j = 12 ,k)Z k 1 z k若Zj工Z

6、j (品n重新分類,現(xiàn)這步的時候,j = 1,2,nClass,k表示迭代次數(shù)），則回到第五步，將全部樣 Z k+ Z k重復迭代計算。若Zj _ Zj , j =1,2,nClass則結束。在實根據需要設置了閾值T,如果改變前后的類中心的差別在閾值k 1 k范圍內則就可以結束，即Zj -Zj nClassTempnClass/鮮且 迭代次數(shù)times為偶數(shù)，則執(zhí)行分裂操作16) 若當前類別數(shù) nClassTemp2*nClass或者 2*nClassnClassTempnClass/并且 迭代次數(shù)times為奇數(shù)，則執(zhí)行合并操作17) 計算各類協(xié)方差向量的最大值和最大值對應的維數(shù)18) 分裂

7、操作：如果某個類別的標準差大于閾值MSE，同時滿足以下條件之一：( 1 )類內平均距離大于所有類別平均距離，并且該類別像素數(shù)大于2*(NMin+1)； (2)當前類別數(shù)nClassTemp=nClass/2根據該類別協(xié)方差向 量最大值對應的維數(shù)S,將該類拆分為兩個新的類別，聚類中心在原來聚類 中心的基礎上 S 維度上加減協(xié)方差向量最大值的一半；19) 計算各類別之間的距離；20) 合并操作：若兩個類別之間的距離小于閾值 DisMin, 合并之，同時保證合并 對數(shù)小于 LMax ；21) 迭代次數(shù)達到閾值，停止迭代，刷新視圖，算法成功。1.6 IHS 變換影像融合1) 傳統(tǒng) HIS 變換算法IH

8、S 正反變換見算法原理介紹， 此處就 I 分量和高分辨率全色影像匹配再做詳細介紹。修改后的 I 分量與原始 I 分量以及全色影像的關系滿足下式：Inew=I*Pan+k2其中Inew修改后的高分辨率影像灰度值；Pan原始全色影像灰度值；k1 增益，由下式計算得到；k1=MSE(I)/MSE(Pan)其中MSE(I)、MSE(Pan)分別是I分量和全色影像的標準差。k2 偏移，由下式計算得到。k2=MEAN(I)-k1* MEAN(Pan)其中MEAN(I)、MEAN(Pan)分別是I分量和全色影像的灰度均值。2) 改進 HIS 變換算法 改進方法，在 HIS 正變換中 I 分量構造是采用 RG

9、B 三個波段的平均值 獲得，這一做法忽略了傳感器的光譜響應，可以采用三波段加權平均來重構I 分量，權重的計算可以采用多光譜影像各波段和全色影像的相關系數(shù)歸一 化獲取，或是采用多元線性回歸獲取。三、實現(xiàn)方法與過程1.7 K 均值算法代碼說明：本程序調用 GDAL 庫編寫，其中 CImage 類的定義此處省略， 具體定義詳見程序源碼。 本程序可以打開和處理各種格式圖像， 經改進后對 圖像波段數(shù)目沒有顯示， CNClassDlg 用來選擇輸入分類書目。為了便于顯 示，對 R、G、 B 三個波段進行顯示，故只默認這可以分 8 類，若要實現(xiàn)更 多類別，只需在以下代碼中分類后附色代碼稍作修改。1) K 均

10、值算法主要代碼 voidCGDALImageProcessDoc:OnKmean()inti,j,k;/ 循環(huán)變量nWidth=img.GetWidht();nHeight=img.GetHeight();intlinebytes=img.lLineBYTES;nBand=img.nBands;BYTE *pData=img.pData;BYTE *pByte=img.m_pDIBs;/ 獲取顯示數(shù)據指針BYTE *DataCopy=new BYTEnWidth*nHeight*nBand;for (k=0;knBand;k+)for (i=0;inHeight;i+)for (j=0;j8)

11、AfxMessageBox（ 輸入類別數(shù)過多 ）;intKMax;/ 循環(huán)次數(shù)上限KMax=1000;int delta;int *pixelClass=new intnWidth*nHeight;/ 判斷每個像素類別 double*meanOld=newdoublenBand*nClass;double *meanNew=new doublenBand*nClass;double *gray=new doublenBand;double *dis=new doublenClass;int *count=new intnClass;for (m=0;mnClass;m+)for (n=0;nn

12、Band;n+)meanOldm*nBand+n=DataCopyn*nWidth*nHeight+m;for (k=0;k0)for (m=0;mnClass;m+)for (n=0;nnBand;n+)meanOldm*nBand+n=meanNewm*nBand+n;/逐像素判斷memset(count,0.0,nClass*sizeof(int);for (i=0;inHeight;i+)for (j=0;jnWidth;j+)/求出各波段像素值，便于計算距離for (m=0;mnBand;m+)graym=DataCopyi*nWidth+j+m*nWidth*nHeight;for

13、 (n=0;nnClass;n+)disn=Distance(gray,&meanOldn*nBand);doubledismin=dis0;int t=0;for (int r=0;rnClass;r+)if (disrdismin)dismin=disr;t=r;/將對應像素復制到以上選擇的類別 ,同時給對應像素附上顏 pixelClassi*nWidth+j=t;countt+;/對各類別像素附上對應的顏色值，此處省略/end of j/end of imemset(meanNew,0.0,nBand*nClass*sizeof(double);for (n=0;nnBand;n+)fo

14、r (i=0;inHeight;i+)for (j=0;jnWidth;j+)int t=pixelClassi*nWidth+j;meanNewt*nBand+n+=(DataCopyi*nWidth+j+n*nWidth*nHeight);for (m=0;mnClass;m+)for (n=0;nnBand;n+)meanNewm*nBand+n/=(countm);delta=0;/delta 表示相鄰兩次迭代聚類中心的偏差for (n=0;nnClass;n+)for (m=0;mnBand;m+)delta+=abs(meanOldn*nBand+m-meanNewn*nBand+

15、m);if (delta0.01) break;/end of k/刷新試圖，此處省略/計算聚類中心的函數(shù)voidCGDALImageProcessDoc:MeanArray(double *mean,double *Ary)int *count=new int6; memset(mean,0.0,6*sizeof(double); memset(count,0.0,6*sizeof(int); intk,i,j;for (k=0;k6;k+)for (i=0;inHeight;i+)for (j=0;jnWidth;j+)if (ArynWidth*i+j+k*nWidth*nHeight!

16、=0)countk+=1;meank+=ArynWidth*i+j+k*nWidth*nHeight;meank/=countk; 計算歐式距離的函數(shù) double CGDALImageProcessDoc:Distance(double *a, double *b) doubledis=(a0-b0)*(a0-b0)+(a1-b1)*(a1-b1)+(a2-b2)*(a2- b2)+(a3-b3)*(a3-b3)+(a4-b4)*(a4-b4)+(a5-b5)*(a5-b5); return dis;1.8 ISODATA 算法程序主要代碼部分 /讀入圖像信息 inti,j,k,t;/ 循環(huán)

17、變量 nWidth=img.GetWidht();nHeight=img.GetHeight(); intlinebytes=img.lLineBYTES;nBand=img.nBands;BYTE *pData=img.pData;BYTE *pByte=img.m_pDIBs;/ 獲取顯示數(shù)據指針BYTE *DataCopy=new BYTEnWidth*nHeight*nBand;for (k=0;knBand;k+)for (i=0;inHeight;i+)for (j=0;jnWidth;j+)DataCopyi*nWidth+j+k*nWidth*nHeight=pDatai*nW

18、idth+j+k*nWidth*nHei ght; intn,m;/m n 分別為類別和波段數(shù)循環(huán)變量 int times=O;迭代次數(shù)intn ClassMax=20;最大類別數(shù)int *pixelClass=new intnWidth*nHeight;/ 判斷每個像素類別 double *meanOld=new doublenBand*nClassMax; double *meanNew=new doublenBand*nClassMax;double *gray=new doublenBand; /某一像素各波段像素值數(shù)組double *dis=new double nClassMax;

19、/某一像素到其對應類別聚類中心的 距離double *mse=new double nClassMax* nClassMax;/給類別協(xié)方差向量 double *mseMax=n ew double nClassMax;/記錄每一類最大協(xié)方差 int *S=new intnClassMax;/ 每一類中協(xié)方差向量的最大維數(shù)int *count=new intnClassMax;/ 某一類別樣本數(shù)目 double *disInter=new doublenClassMax*nClassMax; intnClass,NMin,LMax,IMax,nClassTemp;doubleMSE,DisMi

20、n;step1:for (n=0;nnBand;n+)meanOldn=DataCopyn*nWidth*nHeight;nClassTemp=1;step2:CIsodataParaDlgdlg;if (dlg.DoModal()=IDOK)nClass=dlg.m_K;NMin=dlg.m_Nmin;LMax=dlg.m_LMax;IMax=dlg.m_IMax;M SE=dlg.m_MSE;DisMin=dlg.m_Distance;step3:times+; memset(count,0,nClass*sizeof(int); memset(mse,0.0,nClass*nClass*

21、sizeof(double);/逐像素判斷for (i=0;inHeight;i+)for (j=0;jnWidth;j+)/求出各波段像素值，便于計算距離for (m=0;mnBand;m+)graym=DataCopyi*nWidth+j+m*nWidth*nHeight;if (times=1)for (n=0;nnClassTemp;n+)disn=Distance(gray,&meanOldn*nBand,nBand);elsefor (n=0;nnClassTemp;n+)disn=Distance(gray,&meanNewn*nBand,nBand); doubledismin

22、=dis0;int t=0; for (int r=0;rnClassTemp;r+)if (disrdismin)dismin=disr;t=r;/將對應像素復制到以上選擇的類別 ,同時給對應像素附上顏色 pixelClassi*nWidth+j=t;countt+;if (t=0)pDatai*nWidth+j=0;pDatai*nWidth+j+nWidth*nHeight=0;pDatai*nWidth+j +2*nWidth*nHeight=0;if (t=1)pDatai*nWidth+j=255;pDatai*nWidth+j+nWidth*nHeight=0;pDatai*nW

23、idth +j+2*nWidth*nHeight=0;/此處省略其他類別賦色語句/end of j/end of i /一次迭代完成step4:for (m=0;mnClassTemp;m+)if (countmNMin)/ 某一類別樣本數(shù)目過少，取消該類別for(n=0;nnBand;n+) meanNewm*nBand+n=KONGZHI; nClassTemp-=1;step5: /計算新的聚類中心 memset(meanNew,0.0,nBand*nClassTemp*sizeof(double); for (n=0;nnBand;n+)for (i=0;inHeight;i+)for

24、 (j=0;jnWidth;j+)int t=pixelClassi*nWidth+j;meanNewt*nBand+n+=(DataCopyi*nWidth+j+n*nWidth*nHeight);for (m=0;mnClassTemp;m+)for (n=0;nnBand;n+)meanNewm*nBand+n/=(countm);/ meanOldm*nBand+n=meanNewm*nBand+n;step6: memset(dis,0.0,nClassTemp*sizeof(double); /就算各類中樣本到其聚類中心平均距離 for (i=0;inHeight;i+)for (

25、j=0;jnWidth;j+)m=pixelClassi*nWidth+j;/ 判斷該像素屬于哪個類別for (k=0;knBand;k+)grayk=DataCopyi*nWidth+j+k*nWidth*nHeight;dism+=Distance(gray,&meanNewm*nBand,nBand);step7:/計算所有樣本到其相應聚類中心的距離的平均值double DIS=0,COUNT=0;for (m=0;mnClassTemp;m+)DIS+=dism;COUNT+=countm; dism/=countm;DIS=DIS/COUNT;step8:if (times=IMax

26、)DisMin=0;goto step12;if(nClassTemp=nClass/2&nClassTemp2*nClass)|(nClassTemp=(double)nClass/2&nClassTemp=2*nCla ss×%2=1)/ 合并條件滿足goto step12;step9:/計算各類協(xié)方差向量memset(mse,O.O,nClassMax*nClassMax*sizeof(double);/協(xié)方差向量初始 化for (i=O;inHeight;i+)for (j=O;jnWidth;j+)m=pixelClassi*nWidth+j;for (n=O;nnBand

27、;n+)t=DataCopyi*nWidth+j+n*nWidth*nHeight;msem*nBand+n+=(t-meanNewm*nBand+n)*(t-meanNewm*nBand+n);for (m=0;mnClassTemp;m+)for (n=0;nnBand;n+) msem*nBand+n/=countm; msem*nBand+n=sqrt(msem*nBand+n);step10:/計算各類協(xié)方差向量的最大值和最大值對應的維數(shù)for (m=0;mnClassTemp;m+)mseMaxm=msem*nBand;for (n=0;nmseMaxm) mseMaxm=msem

28、*nBand+n;Sm=n;step11:/分裂運算for (m=0;mMSE)if (dismDIS&countm2*(NMin+1)|nClassTemp=nClass/2)inttempS=Sm;for (n=0;nnBand;n+)/ 構造新的聚類中心meanNew(nClassTemp)*nBand+n=meanNewm*nBand+n; meanNewm*nBand+tempS+=0.5*mseMaxm; meanNew(nClassTemp)*nBand+tempS-=(0.5*mseMaxm); countnClassTemp=countm;n ClassTemp+=1;/新類

29、產生完畢/if(nClassTemp=nClassMax) goto step15;/end of if/end of if/end of forstep12:/計算各類之間的距離for (i=0;inClassMax;i+)for (j=0;jnClassMax;j+)if(i!=j)&(inClassTemp)&(jnClassTemp) disInteri*nClass+j=Distance(&meanNewi*nBand,&meanNewj*nBand,nBand);elsedisInteri*nClass+j=KONGZHI;step13:/對應書上 step13和 step14/合

30、并距離小于閾值的類i=0;for (t=0;tnClassTemp*nClassTemp;t+)if(disIntertDisMin)if(i=LMax) break;i+;int a=t/nClass;int b=t%nClass;a=min(a,b);b=max(a,b);/a 類聚類中心和數(shù)據更新for (n=0;nnBand;n+)meanNewa*nBand+n=(meanNewa*nBand+n*counta+meanNewb*nBand+ n*countb)/(counta+countb);counta=counta+countb;meanNewb*nBand+n=KONGZHI

31、;countb=0;if(b nClassTemp)/后面各類數(shù)據向前移動for (m=b;mnClassTemp;m+)if (m=nClassTemp-1)最后一類數(shù)據清空for (n=0;nnBand;n+)meanNewm*nBand+n=KONGZHI;/end of ifelse if(b nClassTemp-1)/其余各類數(shù)據前移countm=countm+1;for (n=0;nnBand;n+)meanNewm*nBand+n=meanNew(m+1)*nBand+n;/end/end of for/end of if(bnclasstemp)nClassTemp-=1;/

32、end of if distDIS/end of for tstep15:if (times=IMax)goto step16;elsegoto step3;step16:delete dis;delete meanNew; delete meanOld;deletemse; delete disInter; delete gray;deletemseMax;deletecount; deletepixelClass;for (k=0;k3;k+)for (i=0;inHeight;i+)for (j=0;jnWidth;j+)*(pByte+linebytes*(nHeight-1-i)+j

33、*3+2-k)=(pDatai*nWidth+j+k*nWidth*nHei ght);UpdateAllViews(NULL);1.9 IHS 影像融合(篇幅有些，此處只列出融合主函數(shù)和正變換函數(shù))template void CGDALImageProcessDoc:IHSfusion(T)CGDALImageProcessDoc*pDoc=(CGDALImageProcessDoc*)(CMainFrame*)AfxGetMainWnd()-GetActiveDocument();if(!pDoc-imgMulspec.IsValid()AfxMessageBox(open image e

34、rror);BYTE *HighData;HighData=(BYTE*)pDoc-imgFullColor.pData;/ 獲取全色影響數(shù)據T *lpData=(T *)pDoc-imgMulspec.pData;/ 獲取圖像數(shù)據指針BYTE *pByte=pDoc-imgMulspec.m_pDIBs;/ 獲取顯示數(shù)據指針DWORD nWidth=pDoc-imgMulspec.GetWidht();DWORD nHeight=pDoc-imgMulspec.GetHeight();DWORD linebytes=pDoc-imgMulspec.lLineBYTES;DWORD nBand

35、s=pDoc-imgMulspec.nBands;DWORD nWidth2=pDoc-imgFullColor.GetWidht();DWORD nHeight2=pDoc-imgFullColor.GetHeight();if(nWidth!=nWidth2)|(nHeight!=nHeight2)AfxMessageBox(影像大小不匹配);int i,j,k;float *MSData=new floatnWidth*nHeight*nBands;for(k=0;k3;k+)for(j=0;jnHeight;j+)for(i=0;inWidth;i+)ht;MSDataj*nWidth

36、+i+k*nWidth*nHeight=lpDataj*nWidth+i+k*nWidth*nHeigfloat *Intensity=new floatnWidth*nHeight;float *Hue=new floatnHeight*nWidth;float *Saturation=new floatnHeight*nWidth;TransformRGBtoIHS(MSData,Intensity,Hue,Saturation);Match(Intensity,HighData);TransformIHStoRGB(MSData,Intensity,Hue,Saturation);for

37、(k=0;k3;k+)for(i=0;inHeight;i+)for(j=0;jimg.datatype=GDT_Byte)*(lpData+i*nWidth+j+k*nWidth*nHeight)=max(0,min(ftemp,255);else *(lpData+i*nWidth+j+k*nWidth*nHeight)=(T)(ftemp);int nBandsShow=3;for(k=0;kimg.datatype=GDT_Byte)m_maxgray=255;m_mingray=0;elsefor(i=0;inHeight;i+)for(j=0;jnWidth;j+)if(*(lpD

38、ata+i*nWidth+j+k*nWidth*nHeight)m_maxgray)m_maxgray=*(lpData+i*nWidth+j+k*nWidth*nHeight);/else for(i=0;inHeight;i+)for(j=0;j GetActiveDocument();float r,g,b,I,H,S;int i,j;DWORD nWidth=pDoc-imgMulspec.GetWidht();DWORD nHeight=pDoc-imgMulspec.GetHeight(); for(i=0;inHeight;i+)for(j=0;j=b) H=float(acos

39、(f);else H=float(2*PI-acos(f);Intensityi*nWidth+j=I;Huei*nWidth+j=H;Saturationi*nWidth+j=S;四、實驗與分析K均值分類實驗如圖（a） - （c）,圖（b）、（c）分別是分5類和8類的 結果。圖（b）、（c）對比可知，分5類比分8類圖像中類別區(qū)域相對分布比 較集中，噪聲點較少，在應用時需要對圖（c）做簡單的高通濾波處理。ISODATA分類實驗如果（d）-（g）,圖（e）、圖（g）分別對用兩組ISODATA 算法參數(shù)的結果。這些參數(shù)對分類精度會產生一定的影響。HIS變換影像融合實驗如圖（h） - （k）,圖（j） （k、分別為經典IHS變 換融合結果和改進 HIS變換融合結果。由圖可知，改進后的融合算法在保 持影像空間信息的同時，能夠更好地保持圖像的光譜信息。通過計算圖像的 均方根誤差和光譜相關系數(shù)指標（篇幅有限，只考慮了這兩種評價指標），也證明了以上結論，評價結果見鞋標。加權系數(shù) W=0.5028，0.3203, 0.4774信息熵相關系數(shù)平均梯度