VS2010+GDAL+openCV的遥感图像K均值算法的实现和植被指数的计算

基于VS2010+GDAL+openCV的遥感图像K均

值算法的实现和植被指数的计算

摘 要：遥感图像可以理解为记录各种地物不同特征的图像，计算机处理的遥感图像都是数字图像，其中包含各种有价值的，或与地物相关的空间信息，所以对遥感影像的处理尤为重要。本文主要介绍基于VS2010+GDAL+openCV的遥感影像的K均值聚类算法实现和植被指数的计算算法的实现过程与结果分析。K均值算法是一种非监督分类算法，采用迭代、逐次移动聚类中心至得到最后分类结果，其聚类中心数已知，分类结果与初始聚类中心有关，与迭代次数也有关。植被指数的计算则是基于多光谱影像，针对植被等在不同波段下的特点，进行波段间的四则运算来反映关于植被指标变化及分布规律。 关键词：遥感图像处理，K均值，植被指数

1 原理介绍

1.1GDAL与openCV工具

GDAL是一个在X/MIT许可协议下的开源栅格空间数据转换库,可以很好的处理各种遥感影像栅格数据，比如在本次实习中的Erdas Imagine Images(img)格式；openCV则具有强大的处理各种基本数据类型的库，具有各种矩阵数组操作接口，而且可以快速的显示图像。结合GDAL和openCV使用VS环境，可以方便快速的处理各种格式的遥感图像，优化算法过程。

1.2K均值算法

K均值聚类算法是最常用的非监督分类算法之一，主要特点是高效简洁，在遥感影像的分类中广泛使用。K均值算法要求各类样本到聚类中心的距离平方和最小，它是在误差平方和准则的基础上建立起来的。该算法的主要过程为： (1) 任选K个聚类中心：z1(1),z2(1),...,zK(1)，程序中可以以默认的前K个样本作

为初始聚类中心。聚类中心的选择是任意的。如果迭代收敛，则与初始聚类中心的选择无关。

1

(2) 逐个将样本集｛x｝的每一个样本按照最小距离原则分配给K个聚类中心，

形成K类。在第m次迭代中，若：

||x?zj(m)||?|x?zi(m)||,i?1,2,...,K,i?j

则x?fj(m)，这里fj(m)表示第m次迭代时，以第j个聚类中心为代表的聚类域。 (3) 由步骤（2）重新计算积累中心，即：

1zi(m?1)?Nix?fi(m)?x,i?1,2,...,K

Ni表示第i个聚类域中的样本个数。其均值向量作为新的聚类中心。

(4) 如果在第m+1次迭代和第m次迭代过程中，两次的K个迭代中心均相等，则

计算完毕，输出迭代结果，否则返回步骤（2）。

1.3 植被指数

植被指数是指利用卫星不同波段探测数据组合而成的，能反映植物生长状况的指数。已用植被指数评估植被情况的基本原理是植物叶面在可见光红光波段有很强的吸收特性，在近红外波段有很强的反射特性，可以利用在这些波段附近的多光谱影像的波段间的四则运算来消除一些辐射影响，同时增强有用信息。

常用的植被指数有归一化差分植被指数（NDNI），比值植被指数（RVI），土壤调节植被指数（SAVI）,差值环境植被指数（DVI），垂直植被指数（PVI）等。各种植被指数的计算方法为：

(1) NDVI=(NIR-R)/(NIR+R)。用于检测植被生长状态、植被覆盖度和消除部分辐

射误差，可使植被从水和土中分离出来。NIR是近红外波段，R为红外波段。NDVI的值在[-1,1]之间，负值表示地面覆盖为云、水、雪等，对可见光高反射；0表示有岩石或裸土等，NIR和R近似相等；正值，表示有植被覆盖，且随覆盖度增大而增大。在进行计算时，可以利用这一特性，将不同NDVI值的结果赋以不同灰度值，可以直观的得出植被覆盖情况分布图。 (2) RVI=NIR/R。该植被指数能够充分表现植被在红光和近红外波段反射率的差

异，能增强植被与土壤背景之间的差异，但其对大气状况很敏感，而且当植被覆盖率低于50%时，它的分辨能力显著下降。

(3) SAVI=((NIR-R)/(NIR+R+L))(1+L)。其中，L是随着植被密度变化的参数，取

值范围从0-1，当植被覆盖度很高时为0，很低时为1。很明显，如果L=0,SAVI=NDVI。有研究表明，L取0.5时SAVI消除土壤反射率的效果较好。故在本次实习中计算SAVI时取L=0.5。

(4) DVI=NIR-R。该植被指数对土壤背景的变化极为敏感，有利于植被生态环境

的监测，故又被称为环境植被指数（EVI）。

2

(5) PVI=((S R-VR)^2+(SNIR-VNIR)^2)^(1/2)。S是土壤反射率，V是植被反射

率。PVI较好地消除了土壤背景的影响，对大气的敏感度小于其他VI。

2 实现方法与过程

首先构建单文档框架，使用一个悬浮窗口（非模态对话框）以方便直接对图像进行算法前后的操作：

打开图像，包括*.img，*.bmp，*.jpg， *png等格式。 文件路径

图像大小、波段数、位深度和图像格式 投影信息

图像左上角坐标 X，Y方向分辨率

分波段显示和三波段显示

植被指数计算与显示

K均值分类类别数设置与分类显示

2.1各种格式图像的读取及显示

实习采用的遥感图像格式包括Erdas Imagine Images(*.img)格式和*.bmp格式等，使用GDAL工具声明数据集对象指针读取数据：

GDALDataset *poDataset;//数据集对象指针

3

GDALAllRegister();//驱动注册

CPLSetConfigOption(\支持中文路径 poDataset = (GDALDataset*)GDALOpen((const char*)FileName,GA_ReadOnly);

for (int i=0;i

BYTE *pSrcData = new BYTE[m_width*m_height*m_depth/8]; pRasterBand[i]=poDataset->GetRasterBand(i+1);

pRasterBand[i]->RasterIO(GF_Read,0,0,m_width,m_height, pSrcData,m_width,m_height,iDataType,0,0); … }

可以获取图像的长宽、波段数（通道数）、图像格式、数据类型、图像位数深度、每个波段（通道）的数据域，以及对有地理信息的图像可以获取图像投影信息（*.img格式）、图像分辨率、图像左上角的坐标等。

获取到图像所有的基本信息之后，将多波段或单波段图像的每一个波段分别存储在使用openCV申明的IplImage*类型的矩阵之中，选择一个波段或三个波段进行显示，并使用CvvImage类绘制在视图中：

ASSERT_VALID(pDoc); if (!pDoc)return;

CvvImage &img=pDoc->m_img; CRect rc;GetClientRect(&rc);

img.DrawToHDC(pDC->GetSafeHdc(),&rc);

即所有格式的图像均由GDAL读出，保存在IplImage*矩阵中并显示出来。后续对各个波段进行聚类分析以及植被指数的计算中，可以方便直接的对存储在矩阵中的各波段数据进行操作。

针对不同数据类型的图像，对每一个波段，采用最大最小值拉伸的方法将其像素值线型化到[0,255]上：

double mm[2]={0.0};//最大最小值

cvMinMaxLoc(image_in[i],&mm[0],&mm[1]);//求得最小最大值

cvConvertScale(image_in[i],image_out[i],255.0f/(mm[1]-mm[0]),(255.0f-mm[0])/(mm[1]-mm[0]));//拉伸到[0,255]

2.2 K均值算法的实现

K均值算法支持对多波段、不同类型数据进行聚类分析。首先由读取过程中将像素值进行最大最小值拉伸到[0,255]上，再创建样本矩阵，将所有波段的影

4