遥感数字图像处理实习报告

遥感数字图像处理编程实习报告090312009302590123吴煜晖摘要在这次遥感数字图像处理编程实习中，我选取了主分量变换、K均值聚类、ISODATA聚类以及小波变换四个专题来完成。程序我采用VisualC++6.0编译器完成，配合DIB和GDAL这两个类库完成对各种格式图像的读取和写入。整个程序我采用DIB句柄对图像进行操作，当用GDAL读入其他格式图像时我会将其数据转化为DIB句柄，这样无论打开哪种格式图像，都能进行我所完成的专题。关键词DIBGDAL主分量变换聚类小波变换一、原理介绍1、主分量变换（K-L变换）在特征维数较多的情况下，选取其中有利于分类的特征，以减少特征维数，减少内存负担，加快运算速度，提高效率，是自动分类中的一个重要问题。从n维特征选取m维特征，删去的n-m维特征不一定就是无用的信息，如何在信息损失最小的情况下选取特征，在理论上就显得更严密些。通常采用正交变换，得到新的经变换的模式，以保证信息损失最小情况下获得有利于分类的特征。离散K-L变换（Karhunen-Loeve变换）就是常用的方法。2、K均值分类为了将各个样本进行分类，首先确定若干初始聚类中心，然后依一定算法改变或调整这些中心，使它们逐步趋于合理。K均值算法要求各类样本到聚类中心的距离平方和最小，它是在误差平方和准则的基础上建立起来的3、ISODATA分类ISODTA意为迭代自组织数据分析技术。这个算法与K均值算法相似，也是以均值迭代确定聚类中心，但它加进了人机对话环节，可以调整参数，并且引入了归并与分裂的机制，即当某两类别中心间距小于某一阈值时，将它们合并为一类，而在某类样本标准差大于某一阈值时，或其样本数目超过某一阈值时，则将它分为两类，在类别数目少于某一阈值时，也实行分裂。另外，在某类样本数目少于某阈值时，又需要将其消除。因此，这种方法灵活性更强，更为合理。4、小波变换传统的信号理论，是建立在Fourier分析基础上的，而Fourier变换作为一种全局性的变化，其有一定的局限性。在实际应用中人们开始对Fourier变换进行各种改进，小波分析由此产生了。小波分析是一种新兴的数学分支，它是泛函数、Fourier分析、调和分析、数值分析的最完美的结晶；在应用领域，特别是在信号处理、图像处理、语音处理以及众多非线性科学领域，它被认为是继Fourier分析之后的又一有效的时频分析方法。小波变换与Fourier变换相比，是一个时间和频域的局域变换因而能有效地从信号中提取信息，通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析，解决了Fourier变换不能解决的许多困难问题。小波变换分成两个大类：离散小波变换(DWT)和连续小波转换(CWT)。两者的主要区别在于，连续变换在所有可能的缩放和平移上操作，而离散变换采用所有缩放和平移值的特定子集。通常来讲，DWT用于信号编码而CWT用于信号分析。所以，DWT通常用于工程和计算机科学而CWT经常用于科学研究。小波变换现在被大量不同的应用领域采纳，经常取代了傅里叶变换的位置。很多物理学的领域经历了这个範式的转变，包括分子动力学，天文物理学，密度矩阵局部化，地震地质物理学，光学，湍流，和量子力学。其他经历了这种变化的学科有图像处理，血压，心率和心电图分析，DNA分析，蛋白质分析，气象学，通用信号处理，语言识别，计算机图形学，和多分形分析。小波的一个用途是数据压缩。和其他变换一样，小波变换可以用于原始数据(例如图像)，然后将变换后的数据编码，得到有效的压缩。JPEG2000是采用小波的图像标准。在这次实习中，我完成的是离散小波变换。二、算法设计1、主分量变化（k-L变换）K-L变换是方差标准下对数据进行压缩的一种方法。该方法是先求出数据协方差矩阵E的本征值，并将其组成本征向量矩阵V1；然后对本征向量矩阵加以规范正交化，得到新矩阵V2,并用V2作用于原始数据矩阵X,得到K-L变换后的矩阵Y；最后从矩阵Y中提取出若干个较大的本征值对应的向量，舍弃较小本征值对应的向量，即可达到压缩数据的目的。具体做法如下：设X为N维随机向量，即X=(X1,X2,…….XN,)TXi的数学期望为)(iXEXi=1,2,…….N。X的协方差矩阵XEZX{(XX)(X})T为一实对称矩阵，可求出N个本征值),......2,1(Nii及其对应的本征向量V1i,(i=1,2,….N)。则有iiiXVVZ11i=1,2,……N对本征向量矩阵V1做正交变换（规范正交化），即：时当时当jijiVVjTi0111（i,j=1,2,……N）从而可得到规范正交矩阵V2。因此K-L变换式为}{2iTYYXVY或i=1,2,……..N矩阵Y即为K-L变换后的矩阵，它是原始矩阵X中各种信息的重新分配和组合。Y中的每个向量包含有X中所有向量的信息。其中协方差矩阵中的若干个较大的本征值对应的向量包含原始矩阵X中各向量的信息较多，且是主要的信息。因此，从中可把含有较多信息的向量提取出来，而舍弃那些含原来信息较少，亦即次要的向量，从而达到压缩数据的目的。压缩后的矩阵写为}{iYYi=1,2,…..k；kN由于提取的若干个向量仍包含原来每个向量的信息，而且是主要的信息是，所以这样做相当于对原始数据进行了滤波，使信息变得更加分明，这对于聚类分析将是很有用的。对于这次实习，我们把6波段图像的6个波段的灰度值作为6维向量，每个像素点为一个向量。2、K均值分类K均值算法能使聚类域中所有样品到聚类中心的距离平方和最小。其主要步骤主要如下：第一步：任选k个初始聚类中心：11Z,12Z,…,1kZ(上角标记载为寻找聚类中的迭代运算次数)。用数组classp[6*clsnumber]来存储类中心的值，一般可选定样品集的前k个样品作为初始聚类中心。但是考虑到这样做不太有利于后面的算法收敛。因此采用了最大最小距离选心法。该法的原则是使各初始类别之间，尽可能地保持远离。任意选取M个初始中心，将其值存入iGrayValue[Bands*M]中；将第一个点X1作为第一个初始类别的中心Z1。计算X1与其它各抽样点的距离D。取与之距离最远的那个抽样点(例如X7)为第二个初始类别中心Z2则第二个初始类中心Z2=X7。对剩余的每个抽2样点，计算它到已有各初始类别中心的距离ijD(u=1,2,…，已知有初始类别数m),并取其中的最小距离作为该点的代表距离Dj,jD=min(jD1,jD2…,mjD)；在此基础上，再对所有各剩余点的最小距离jD进行相互比较，取其中最大者，并选择与该最大的最小距离相应的抽样点(如11X)作为新的初始类中心点，即3Z=11X。此时m=m+1;如此迭代直到m=clsnumber;即m=0,1,2…,clsnumber;第二步：设已进行到第k步迭代。若对某一样品X有|X－kjZ|＜|X－kiZ|,则X∈kjS,其中kjS是以kjZ为聚类中心的样品集。以此种方法将全部样品分配到k个类中。即确定每个象素的类属k7中，如k7=3,即表示该象素属于第3类；并相应的对其赋值到array数组中，以便以后可以显示其分类结果。第三步：计算各聚类中心的新向量值classo[i];1kjZ＝jn1kjsxX（j=1,2,…,k）classo[i]=classo1[i]/NL[i];式中jn为jS中所包含的样品数。classo1[i]表示所有属于第i类的像素的值的累加。NL[i]表示属于第i类的像素总数。classo[i]重新分类后的聚类中心值。因为在这一步要计算k个聚类中心的样品均值，故称为k均值算法。第四步：若1kjZ≠kjZ,j=1,2,…,k,则回到第二步，将全部样品n重新分类，重复迭代计算。若1kjZ＝kjZ,j=1,2,…,k,则结束。在实现这步的时候，根据需要设置了阈值thresholdc,如果改变前后的类中心的差别在阈值范围内则就可以结束。即((|classo[i]-classp[i]|)/ckassp[i])threshold即可以结束算法。本次实习中，我们将六波段图像的六个波段值作为一个像素点的位置，取六维空间中的欧氏距离作为每两个点之间的距离。3、ISODATA分类与刚介绍的K均值分类类似，这里的ISODATA分类也是一种非监督分类方法，只是步骤有所不同：第一步：将个模式样本{,i=1,2,3,…,}读入。第二步：将N个模式样本分给最近的聚类，假如Dj=min(‖x-zj‖,i=1,2,…,),即‖x-zj‖的距离最小，则x∈Sj。第三步：如果Sj中的样本数NjΘn,取消样本子集。第四步：修正聚类中心值j=1,2,…,第五步：计算各聚类域Sj中诸聚类中心间的平均距离：第六步：计算全部模式样本对其相应聚类中心的总平均距离：第七步：判别分裂、合并及迭代运算等步骤：①如迭代运算次数已达I次，即最后一次迭代，置θc=0,跳到第十一步，运算结束。②如≤K/2，即聚类中心的数目等于或不到规定值的一半，则进入第八步，将已有的聚类分裂。③如迭代运算的次数是偶次，或≥2K，不进行分裂处理，跳到第十一步；如不符合以上两个条件（即既不是偶次迭代，也不是≥2K），则进入第八步，进行分裂处理。分裂处理：第八步：计算聚类样本距离的标准差向量：第九步：求每一标准差向量{，σj=1,2,…，}中的最大分量，以{σj=1,2,…，}代表。第十步：在任一最大分量集{σj=1,2,…，}中，如有θS（该值给定），同时又满足以下二条件中之一：(a)，即Sj中样本总数超过规定值一倍以上，（b）Nc≤K/2，则将Zj分裂为两个新的聚类中心,且加1。中相当于的分量，可加上kσjmax，其中0≤k≤1；中相当于的分量，可减去kσjmax。如果本步完成了分裂运算，则跳回第二步；否则，继续。第十一步：计算全部聚类中心的距离：,i=1,2,…,-1j=i+1,…,第十二步：比较与θc值，将θc的值按最小距离次序递增排列，即{,,…,}式中，…。第十三步：如将距离为Di1j1的两个聚类中心zi1和zj1合并，得新中心为l=1,2,…，L式中，被合并的两个聚类中心向量，分别以其聚类域内的样本数加权，使为真正的平均向量。第十四步：如果是最后一次迭代运算（即第I次），算法结束。否则GOTO第一步——如果需由操作者改变输入参数；或GOTO第二步——如果输入参数不变。4、小波变换首先我们定义一些需要用到的信号及滤波器。x[n]：离散的输入信号，长度为N。g[n]：lowpassfilter低通滤波器，可以将输入信号的高频部份滤掉而输出低频部份。h[n]：highpassfilter高通滤波器，与低通滤波器相反，滤掉低频部份而输出高频部份。Q：downsamplingfilter降采样滤波器，如果以x[n]作为输入，则输出y[n]=x[Qn]。此处举例Q=2。清楚规定以上符号之后，便可以利用阶层架构来介绍如何将一个离散信号作离散小波变换：架构中的第1层(1ststage)架构中的第2层(2ndstage)可继续延伸架构中的第α层(αstage)离散小波变换如下：此时的输入信号变成x[m,n]，而转换过程变得更复杂，说明如下：首先对n方向作高通、低通以及降频的处理接着对v1,L[m,n]与v1,H[m,n]延著m方向作高低通及降频动作经过以上两个步骤才算完成离散小波变换的一个stage。三、实现方法与过程1、主体框架设计及本程序核心程序由于我在先期设计时计划采用DIB类对位图进行操作，所以主体采用DIB，但在后期为了实现打开其他格式图像（诸如img、jpg）等格式，使用了GDAL库，由于已经进入后期，不想对前期的程序进行修改，于是在用GDAL时，每打开一幅图像将其自动写为bmp格式，退出时不保存则不存储在硬盘内。而由于我又在聚类及主分量变换中实现了打开多幅图像，并对其进行操作的功能，因为每打开一幅新的图像，就相当于新建了一个文档视图类，于是我将数据（仅限于聚类及主分量分析时）存在App类中，方法如下：（1）分别在工作区中建立lib和include文件夹，将DI