遥感原理实习报告

1、遥感原理课程设计专业：遥感科学与技术班级：学号：姓名：指导老师：一、原理介绍1.1 K-均值算法算法主要思想是先在需要分类的数据中寻找K组数据作为初始聚类中心，然后计算所有数据到各聚类中心的距离，将数据归入与其距离最近的聚 类中心，之后再对这 K个聚类的数据计算均值，作为新的聚类中心，重复 以上步骤，直到新的聚类中心与上一次的聚类中心值相等时结束算法。算法 流程图如下：否1.2 ISODATA （迭代自组织分析）通过设定初始参数而引入人机对话环节，并使用归并与分裂的机制，当 某两类聚类中心距离小于某一阈值时，将它们合并为一类，当某类标准差大 于某一阈值或其样本数目超过某一阈值时，将其分为两类。

2、在某类样本数目 少于某阈值时，需将其取消。如此，根据初始聚类中心和设定的类别数目等 参数迭代，最终得到一个比较理想的分类结果。算法流程图如下：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、波段数nBand、数据首地 址pData、显示数据指针pByte；2）将图像

4、像素数据复制一份到 DataCopy中，之后用DataCopy处理，处理结果 重新赋值到pData中，刷新即可得到分类结果图，从而避免了新建图像这一 复杂过程。3）调用对话框，获取分类类别数nClass,并定义中间变量参数：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;/ 每个类别像素数量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重新分类,现这步的时候,j = 1,2,nClass,k表示迭代次数），则回到第五步，将全部样 Z k+ Z k重复迭代计算。若Zj _ Zj , j =1,2,nClass则结束。在实根据需要设置了阈值T,如果改变前后的类中心的差别在阈值k 1 k范围内则就可以结束，即Zj -Zj nClassTempnClass/鲜且 迭代次数times为偶数，则执行分裂操作16) 若当前类别数 nClassTemp2*nClass或者 2*nClassnClassTempnClass/并且 迭代次数times为奇数，则执行合并操作17) 计算各类协方差向量的最大值和最大值对应的维数18) 分裂

7、操作：如果某个类别的标准差大于阈值MSE，同时满足以下条件之一：( 1 )类内平均距离大于所有类别平均距离，并且该类别像素数大于2*(NMin+1)； (2)当前类别数nClassTemp=nClass/2根据该类别协方差向 量最大值对应的维数S,将该类拆分为两个新的类别，聚类中心在原来聚类 中心的基础上 S 维度上加减协方差向量最大值的一半；19) 计算各类别之间的距离；20) 合并操作：若两个类别之间的距离小于阈值 DisMin, 合并之，同时保证合并 对数小于 LMax ；21) 迭代次数达到阈值，停止迭代，刷新视图，算法成功。1.6 IHS 变换影像融合1) 传统 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 分量，权重的计算可以采用多光谱影像各波段和全色影像的相关系数归一 化获取，或是采用多元线性回归获取。三、实现方法与过程1.7 K 均值算法代码说明：本程序调用 GDAL 库编写，其中 CImage 类的定义此处省略， 具体定义详见程序源码。 本程序可以打开和处理各种格式图像， 经改进后对 图像波段数目没有显示， CNClassDlg 用来选择输入分类书目。为了便于显 示，对 R、G、 B 三个波段进行显示，故只默认这可以分 8 类，若要实现更 多类别，只需在以下代码中分类后附色代码稍作修改。1) K 均

10、值算法主要代码 voidCGDALImageProcessDoc:OnKmean()inti,j,k;/ 循环变量nWidth=img.GetWidht();nHeight=img.GetHeight();intlinebytes=img.lLineBYTES;nBand=img.nBands;BYTE *pData=img.pData;BYTE *pByte=img.m_pDIBs;/ 获取显示数据指针BYTE *DataCopy=new BYTEnWidth*nHeight*nBand;for (k=0;knBand;k+)for (i=0;inHeight;i+)for (j=0;j8)

11、AfxMessageBox（ 输入类别数过多 ）;intKMax;/ 循环次数上限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/刷新试图，此处省略/计算聚类中心的函数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; 计算欧式距离的函数 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;/ 循环

17、变量 nWidth=img.GetWidht();nHeight=img.GetHeight(); intlinebytes=img.lLineBYTES;nBand=img.nBands;BYTE *pData=img.pData;BYTE *pByte=img.m_pDIBs;/ 获取显示数据指针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 分别为类别和波段数循环变量 int times=O;迭代次数intn ClassMax=20;最大类别数int *pixelClass=new intnWidth*nHeight;/ 判断每个像素类别 double *meanOld=new doublenBand*nClassMax; double *meanNew=new doublenBand*nClassMax;double *gray=new doublenBand; /某一像素各波段像素值数组double *dis=new double nClassMax;

19、/某一像素到其对应类别聚类中心的 距离double *mse=new double nClassMax* nClassMax;/给类别协方差向量 double *mseMax=n ew double nClassMax;/记录每一类最大协方差 int *S=new intnClassMax;/ 每一类中协方差向量的最大维数int *count=new intnClassMax;/ 某一类别样本数目 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)/ 某一类别样本数目过少，取消该类别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:/计算各类协方差向量memset(mse,O.O,nClassMax*nClassMax*sizeof(double);/协方差向量初始 化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:/计算各类协方差向量的最大值和最大值对应的维数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 类聚类中心和数据更新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)/后面各类数据向前移动for (m=b;mnClassTemp;m+)if (m=nClassTemp-1)最后一类数据清空for (n=0;nnBand;n+)meanNewm*nBand+n=KONGZHI;/end of ifelse if(b nClassTemp-1)/其余各类数据前移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 影像融合(篇幅有些，此处只列出融合主函数和正变换函数)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;/ 获取全色影响数据T *lpData=(T *)pDoc-imgMulspec.pData;/ 获取图像数据指针BYTE *pByte=pDoc-imgMulspec.m_pDIBs;/ 获取显示数据指针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类图像中类别区域相对分布比 较集中，噪声点较少，在应用时需要对图（c）做简单的高通滤波处理。ISODATA分类实验如果（d）-（g）,图（e）、图（g）分别对用两组ISODATA 算法参数的结果。这些参数对分类精度会产生一定的影响。HIS变换影像融合实验如图（h） - （k）,图（j） （k、分别为经典IHS变 换融合结果和改进 HIS变换融合结果。由图可知，改进后的融合算法在保 持影像空间信息的同时，能够更好地保持图像的光谱信息。通过计算图像的 均方根误差和光谱相关系数指标（篇幅有限，只考虑了这两种评价指标），也证明了以上结论，评价结果见鞋标。加权系数 W=0.5028，0.3203, 0.4774信息熵相关系数平均梯度