遥感数字图像处理编程实习报告
09031 2009302590123 吴煜晖
摘要
在这次遥感数字图像处理编程实习中，我选取了主分量变换、K均值聚类、ISODATA聚类以及小波变换四个专题来完成。程序我采用Visual C++ 6.0编译器完成，配合DIB和GDAL这两个类库完成对各种格式图像的读取和写入。整个程序我采用DIB句柄对图像进行操作，当用GDAL读入其他格式图像时我会将其数据转化为DIB句柄，这样无论打开哪种格式图像，都能进行我所完成的专题。 关键词
DIB GDAL 主分量变换 聚类 小波变换
一、原理介绍
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分析，蛋白质分析，气象学，通用信号处理，语言识别，计算机图形学，和多分形分析。
小波的一个用途是数据压缩。和其他变换一样，小波变换可以用于原始数据(例如图像)，然后将变换后的数据编码，得到有效的压缩。JPEG 2000 是采用小波的图像标准。
在这次实习中，我完成的是离散小波变换。
二、算法设计
1、主分量变化（k-L变换）
K-L变换是方差标准下对数据进行压缩的一种方法。该方法是先求出数据协方差矩阵E的本征值，并将其组成本征向量矩阵V1；然后对本征向量矩阵加以规范正交化，得到新矩阵V2,并用V2 作用于原始数据矩阵X,得到K-L变换后的矩阵Y；最后从矩阵Y中提取出若干个较大的本征值对应的向量，舍弃较小本征值对应的向量，即可达到压缩数据的目的。具体做法如下：
设X为N维随机向量，即 X = (X1, X2,……. XN,)T
Xi 的数学期望为 X?E(Xi) i=1,2,…….N。 X的协方差矩阵
ZX?E{(X X)(X X)T}
为一实对称矩阵，可求出N个本征值?i(i?1,2,......N)及其对应的本征向量V1i,(i=1,2,….N)。则有
ZXV1i??iV1i i=1,2,……N
对本征向量矩阵V1做正交变换（规范正交化），即：
?1当i?j时? V1T? （i, j = 1,2,……N） iV1j???0当i?j时?从而可得到规范正交矩阵V2 。因此K-L变换式为
Y?V2TX或Y?{Yi} i=1,2,……..N
矩阵Y即为K-L变换后的矩阵，它是原始矩阵X中各种信息的重新分配和组合。Y中的每个向量包含有X中所有向量的信息。其中协方差矩阵中的若干个较大的本征值对应的向量包含原始矩阵X中各向量的信息较多，且是主要的信息。因此，从中可把含有较多信息的向量提取出来，而舍弃那些含原来信息较少，亦即次要的向量，从而达到压缩数据的目的。压缩后的矩阵写为
Y?{Yi} i=1,2,…..k ；k 由于提取的若干个向量仍包含原来每个向量的信息，而且是主要的信息是，所以这样做相当于对原始数据进行了滤波，使信息变得更加分明，这对于聚类分析将是很有用的。 对于这次实习，我们把6波段图像的6个波段的灰度值作为6维向量，每个像素点为一个向量。 2、K均值分类 K均值算法能使聚类域中所有样品到聚类中心的距离平方和最小。其主要步骤主要如下： 11第一步：任选k个初始聚类中心：Z11,Z2,…,Zk(上角标记载为寻找聚类中 的迭代运算次数)。用数组classp[6*clsnumber]来存储类中心的值，一般可选定样品集的前k个样品作为初始聚类中心。但是考虑到这样做不太有利于后面的算法收敛。因此采用了最大最小距离选心法。该法的原则是使各初始类别之间，尽可能地保持远离。 任意选取M个初始中心，将其值存入iGrayValue[Bands*M]中；将第一个点X1作为第一个初始类别的中心Z1。计算X1与其它各抽样点的距离D。取与之距离最远的那个抽样点(例如X7)为第二个初始类别中心Z2则第二个初始类中心Z2=X7。 对剩余的每个抽2样点，计算它到已有各初始类别中心的距离 Dij(u=1,2,…，已知有初始类别数m),并取其中的最小距离作为该点的代表距离Dj,Dj=min(D1j,D2j…,Dmj)； 在此基础上，再对所有各剩余点的最小距离Dj进行相互比较，取其中最大者，并选择与该最大的最小距离相应的抽样点(如X11)作为新的初始类中心点，即 Z3=X11。 此时m=m+1; 如此迭代直到m>=clsnumber;即m = 0,1,2…,clsnumber; k第二步：设已进行到第k步迭代。若对某一样品X有|X－Zkj|＜|X－Zi|, kk则X∈Sk以此种方法将全部样品分配到j,其中Sj 是以Zj为聚类中心的样品集。 k个类中。即确定每个象素的类属k7中，如 k7= 3,即表示该象素属于第3类；并相应的对其赋值到array数组中，以便以后可以显示其分类结果。 第三步：计算各聚类中心的新向量值classo[i]; ?1Zkj＝ 1njx?skj?X （ j = 1,2,…,k） classo[i] = classo1[i]/NL[i]; 式中nj为Sj中所包含的样品数。classo1[i]表示所有属于第i类的像素的值的累加。 NL[i]表示属于第i类的像素总数。classo[i]重新分类后的聚类中心值。 因为在这一步要计算k个聚类中心的样品均值，故称为k均值算法。 ?1k第四步：若Zk≠Zjj, j = 1,2,…,k,则回到第二步，将全部样品n重新分?1类，重复迭代计算。若Zk＝Zkjj, j =1,2,…,k,则结束。在实现这步的时候， 根据需要设置了阈值thresholdc,如果改变前后的类中心的差别在阈值范围内则就可以结束。即 ((|classo[i]-classp[i]|)/ckassp[i]) 本次实习中，我们将六波段图像的六个波段值作为一个像素点的位置，取六维空间中的欧氏距离作为每两个点之间的距离。 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) （b）Nc≤K/2， 则将Zj分裂为两个新的聚类中心可加上kσjmax，其中0<≤k≤1； ,且加1。中相当于 中相当于 的分量， ，即Sj中样本总数超过规定值一倍以上， 的分量，可减去kσjmax。 如果本步完成了分裂运算，则跳回第二步；否则，继续。 第十一步：计算全部聚类中心的距离： ,i=1,2, …, -1 j=i+1, …, 第十二步：比较即 与θc值，将<θc的值按最小距离次序递增排列，