SRTP结题报告

西南交通大学第五期大学生科研训练计划（SRTP）遥感与气象观测数据耦合的地表干旱状态监测结题报告2010年4月至2011年4月目录目录11绪论21. 1项目背景21. 2项目介绍31. 3 开源库GDAL与HDF文件格式介绍42数据源及GDAL库的配置72. 1 数据源72. 2 GDAL库的配置93数据处理原理及实现103. 1 遥感数据处理103.2 气象数据处理与降水空间插值法比较203. 3 指数耦合263. 4 程序设计274实验结果分析275不足与期望296项目感想31参考文献32附录1 程序主要源代码341 绪论1. 1项目背景我国是一个旱灾非常严重的国家，旱灾给农业、农村和农民造成了巨大的损失。据统计，19501990年间，我国共有11年发生了重特大干旱，发生频次为26.83，因干旱造成粮食损失占粮食总产量的4.02。而19912009年间，我国共有8年发生重特大干旱，因干旱造成粮食损失占粮食总产量的6.09。近年来，我国年年有干旱，平均不到3年发生一次重特大旱灾，尤其经常发生区域性特大旱灾。同时我国旱灾分布面积广，过去，我国旱灾高发的区域主要在干旱缺水的北方地区，特别是西北地区。近几年，在传统的北方旱区旱情加重的同时，南方和东部多雨区旱情也在扩展和加重，目前旱灾范围已遍及全国。与此同时，旱灾影响范围已由传统的农业扩展到工业、城市、生态等领域，工农业争水、城乡争水、超采地下水和挤占生态用水现象越来越严重。因此实现对旱灾的监测，有利于实施抗旱措施，同时合理分配水资源，节约用水，对促进农业生产、保障粮食安全和区域可持续发展具有重要的现实意义胡其峰.科学理性对待西南旱灾减灾需立足长远. 光明日报.2010.4。传统的干旱监测采用气象数据，气象数据由气象站观测得到，通过建立区域气象站网同时辅以水文、社会和经济等数据实施干旱监测。这种方法应用十分广泛，在单点精度很高，但其性质是以点带面来监测，这样很难给出不同干旱区域的分界线，同时整体精度也相对较低，难以适应大范围的干旱监测。随着卫星遥感技术的迅速发展，国内外已在借助遥感手段进行大范围的干旱监测方面开展了大量的研究和实际应用。基于遥感手段的干旱监测能频繁和持久的提供地表特征的面状信息，具有宏观，高时效，经济等特点，且点点俱到，适应大范围的干旱监测。当然，由于气象卫星运行特点和遥感传感器本身性能的限制，该方法也存在一定的缺点和问题，比如目前相应的基于遥感数据的一些反演模型不够成熟，并且其受大气影响较大，在很多有云的情况下，数据不能使用，其监测的结果在小范围内也不如传统的气象监测精确。针对单一遥感数据在进行地表干旱状态监测时的不足，考虑降水对土壤水分的阶段性影响，我们应该把两种方法结合在一起考虑，形成地表温度、植被状态、空间降水等参数耦合的干旱监测模型。这样能提高干旱监测的准确性和及时性。1. 2项目介绍项目基于对地观测MODIS数据与气象站点实测资料，建立起多源数据耦合下的地表干旱监测方法与相关指标体系，在开源GDAL库的基础上，开发可用于大范围区域地表干旱状态监测的应用系统。项目中涉及四个方面的内容：MODIS预处理算法程序设计、气象观测资料处理、热红外地表温度反演模型程序设计、多元数据耦合模型参数化方案及设计。其中，MODIS预处理算法程序设计部分包括HDF格式数据读取、MODIS几何校正与辐射标定处理。气象观测资料处理包括气象观测站点降水资料读取、空间插值与栅格化。热红外地表温度反演程序设计主要实现利用分裂窗法的地表温度反演模型。耦合模型部分则包括了地表温度参数与降水参数的耦合模型参数计算、模型程序实现。本项目主要分为三个部分：遥感数据处理，气象数据处理以及干旱监测指数的耦合。遥感数据处理：利用开源库GDAL提取MODIS 1B的数据，在此基础上利用MODIS数据中的相应波段反演地表温度（LST）和归一化植被指数（NDVI），得到遥感干旱监测因子作物供水指数（VSWI）。气象数据处理：利用ArcGIS的二次开发平台ArcGIS Engine实现气象数据降雨量的多种空间插值算法。在比较不同空间插值算法误差的基础上，综合多年数据选取最佳的空间插值算法进行降雨量的插值运算，并计算出气象干旱监测因子综合降水距平指数（MSRI）。干旱监测指数的耦合：建立作物供水指数和综合降水距平指数的耦合模型，得到农业旱情监测指数，并划分干旱等级。项目流程图如图1所示：遥感数据-MODIS1B数据地面观测数据-降雨量云检测NDVILSTVSWINDVI/LST单日数据内插旬合成当旬降水量距平综合降水距平指数旬合成综合旱情监测指数几何校正图1-1项目流程图1. 3 开源库GDAL与HDF文件格式介绍1. 3. 1 GDAL库介绍在MODIS数据的处理中，我们没有选用HDF Group开发的HDF函数库，而是选用了十分流行的开源库GDAL。因为开源库GDAL在处理栅格图像格式有很大优势，其不仅读取数据较快，而且能够把结果转换成任一其支持的栅格格式。1998年末，加拿大的Frank Warmerdam开始了GDAL(Geospatial Data Abstraction Library)项目的编写工作。该项目得到了许多个人和团队的支持。GDAL是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。OGR是GDAL项目的一个分支，功能和GDAL类似，只不过他提供对矢量数据的支持。因此GDAL能同时提供栅格和矢量数据的操作。有很多著名的GIS类产品包括ESRI的ARCGIS，Google Earth和跨平台的GRASS GIS系统都使用了GDAL/OGR开源库的若干模块。利用GDAL/OGR库，可以使基于Linux的地理空间数据管理系统提供对矢量和栅格文件数据的支持。GDAL提供对多种栅格数据的支持，包括Arc/Info ASCII Grid(asc)，GeoTiff (tiff)，Erdas Imagine Images(img)，ASCII DEM(dem) 等格式。提供读取、写入、转换、处理各种栅格数据格式（有些特定的格式对一些操作如写入等不支持）操作。GDAL使用抽象数据模型(Abstract Datamodel)来解析它所支持的数据格式，抽象数据模型包括数据集(Dataset)，坐标系统，仿射地理坐标转换(Affine GeoTransform)，大地控制点(GCPs)，元数据(Metadata)，栅格波段(Raster Band)，颜色表(ColorTable)，子数据集域(Subdatasets Domain)，图像结构域(Image_StructureDomain)，XML域(XML:Domains)。GDAL库采用ANSI C和C+语言编写，能通过所有c/c+编译器的编译。其核心类框架如图2所示GDALMajorObjectGDALDatasetGDALDriveGDALDriverManagerGDALRasterBand图12 GDAL核心类框架其中：GDALMajorObject类：带有元数据的对象。GDALDataset类：通常是从一个文件中提取的相关栅格波段数据的集合，以及其相关元数据。GDALDriver类：文件格式驱动类。GDALDriverManager类：文件格式驱动管理类。GDALRasterBand类：栅格波段数据类。1. 3. 2 HDF文件格式HDF（Hierarchical data format）格式是美国国家高级计算机应用中心为了满足各种领域研究需求而研制的一种高效存储和分发科学数据的层次数据格式，主要用来存储由不同计算机平台来产生的各种类型科学数据，便于在不同的计算机平台上扩展。HDF被设计为：自述性：对于一个HDF文件里的每一个数据对象，有关于该数据的综合信息。通用性：许多数据类型都可以被嵌入在一个HDF文件里。灵活性：HDF允许用户把相关的数据对象组合在一起，放到一个分层结构中，向数据对象添加描述和标签。它还允许用户把科学数据放到多个HDF文件里。扩展性：HDF极易容纳将来新增加的数据模式，容易与其他标准格式兼容。跨平台性：HDF是一个与平台无关的文件格式。HDF文件无需任何转换就可以在不同平台上使用。目前HDF有两个版本，即HDF4和HDF5，两个版本的数据结构模型变化不大，但其具体实现方法变化很大，因此这两个版本是不兼容的。由于GDAL库采用一个抽象的数据模型来支持遥感栅格图像，因此使用GDAL读取这两个版本的HDF数据的方法基本相同。一个HDF文件包括一个文件头，至少一个数据描述块以及若干个数据元素（图1-3）。HDF文件像一本带目录的书，数据描述块可以看成书的目录，而书的每一章内容则是HDF文件的数据元素。数据元素HDF文件头空数据描述符数据描述符空数据描述符空数据描述符HDF文件数据描述块图13 HDF文件格式本项目采用的MODIS 1B数据采用HDF文件格式储存。图1-4 是通过HDFView软件浏览器查看某MODIS 1km分辨率的文件列表示意图。波段：5*2040*1354卫星成像参数：408*271波段：15*2040*1354波段：2*2040*1354图14 HDF（MODIS）文件列表示意图2 数据源及GDAL库的配置2. 1 数据源2. 1. 1遥感数据源中等分辨率成像光谱仪MODIS数据采用HDF(Hierarchical data format)数据格式存储元数据，提供36个波段的地球综合信息，分布在0.414m的电磁波谱范围内，对开展自然灾害与生态环境监测、全球环境和气候变化研究以及进行全球变化的综合性研究等有重要意义。MODIS仪器的地面分辨率为250m、500m和1000m，扫描宽度为2330km。在对地观测过程中，每秒可同时获得6.1兆比特的来自大气、海洋和陆地表面信息，每日或每两日可获取一次全球观测数据。本项目研究的数据来源于NASA公开下载的MODIS 1B数据。本研究中用1000m分辨率的数据，主要在于1000m分辨率的数据包含了36个波段全部数据，而且其分辨率大小适中。所获取的数据是2009年8月1日到2009年8月30日白天数据，范围覆盖四川省。2. 1. 2气象数据源气象数据主要有四川省12个降雨量观测站2000年到2009年每月3旬的气象观测数据，站点分布如下图所示：图2-1气象数据站点分布图2. 1. 3 其它数据源四川省省界来源于国家基础地理信息中心（National Geomatics Center of China），另外还有中国区域的数字高程模型（DEM，分辨率为1km，来源于NASA）、投影全部采用埃尔伯斯等积圆锥投影（Albers Conical Equal Area），这些基础数据用于生成旱情监测系统的基础数据库。2. 2 GDAL库的配置GDAL库默认不支持项目所用的HDF文件格式。要获得支持HDF文件格式的GDAL，需要向GDAL源代码中，添加HDF文件驱动重新编译。配置过程如下：a) 从官网下载开源库GDAL (1.6版本)和HDF4函数库(42r4版本)。HDF4函数库放在”C:GDAL1.6HDF4”目录下。b) 打开gdal文件夹下的nmake.opt，修改GDAL_HOME = C:GDAL1.6，把路径改到需要把gdal安装的地方。c) 添加HDF4驱动。编辑gdal根目录下的nmake.opt，找到“# Uncomment the following and update to enable NCSA HDF Release 4 support.”这一行。把下面两行前面的#去掉，然后改成：HDF4_DIR = C:GDAL1.642r4-winHDF4_LIB = $(HDF4_DIR)dllhd424m.lib $(HDF4_DIR)dllhm424m.lib.$(HDF4_DIR)libhd424.lib $(HDF4_DIR)libhm424.lib完成后保存并关闭。d) 打开控制台，输入: ”D:Program FilesMicrosoft Visual Studio .NET 2003Vc7binvcvars32.bat 注册vc的编译环境。依次运行：nmake /f makefile.vcnmake /f makefile.vc installnmake /f makefile.vc devinstall如果编译不成功，这时有可能是GDAL代码自带的BUG。需要修改gdalfrmtslevellerlevellerdataset.cpp文件171行： “?, kPI / 180.0, UNITLABEL_DEGREE 将“ “? ”修改为“ ”? “ ”，保存后重新编译。e)编译好后，把C:GDAL1.6bin添加到系统环境变量中。并在项目预定义内容中添加以下内容：#include C:GDAL1.6includegdal_priv.h#include C:GDAL1.6includegdal.h#include C:GDAL1.6includecpl_string.h#include C:GDAL1.6includecpl_conv.h#pragma comment (lib,C:/GDAL1.6/lib/gdal_i.lib) /编译时链接GDAL的lib库这样在项目中就能使用支持HDF格式的GDAL库了。3 数据处理原理及实现3. 1 遥感数据处理干旱是降水时空分布不均，长时间缺乏降水而形成的气象现象。这样就会造成局部地区作物生长过程中因供水不足，阻碍作物的正常生长而发生的水量供应不足的现象。一般来讲，植被的生产状况主要与水分有关，水分供应程度便成了作物生长的关键因素，水分供应充足，植被生长良好，反之生长变差。植被指数的时空变化与土壤水分状况有一定的相关性，因此植被指数可以作为干旱监测的因素。植被指数是卫星传感器不同通道探测数据的线性或非线性组合，能够反映绿色植物生长和分布的特征指数。本项目采用十分常用的归一化植被指数NDVI(Normalized Difference Vegetation Index)，其定义为：NDVI=(CH1-CH2)/(CH1+CH2) (1)其中，CH1和CH2分别表示近红外波段和红光波段的反射率。地表温度(Land Surface Temperature, LST)也可用于干旱监测。地表温度是控制地球表面大多数物理、化学和生物过程的参数之一。影响地表温度变化的因素也比较多，比如地表湿度、气温、光照强度、地表材质（比如是草坪还是裸露土地，还是水泥地面，或者是沥青地面）等。综合考虑上述的两个干旱监测因子，本项目采用植被供水指数VSWI作为遥感干旱监测因子，其定义为：VSWI=NDVI/Ts (2)VSWI的物理意义：当植被供水正常时，卫星遥感的植被指数在一定的生长期内保持在一定的范围，而卫星遥感的植被冠层温度也保持在一定的范围，如果遇到干旱，植被供水不足，一方面植被的生长受到影响，卫星遥感的植被指数将降低，另一方面当植被受旱胁迫时，为减少水分损失，叶面的气孔会部分关闭，导致叶面温度的升高，从而植被冠层温度将升高。因此可以用植被供水指数来评价旱情状况农业科学院农业资源与农业区划研究所. 全国旱情监测改进方法.。遥感干旱监测模块流程图如下所示：MODIS 第1和第2波段植被指数NDVI地表比辐射率31和32大气水汽含量MODIS 第2和19波段大气透过率31()和32()MODIS第31和32波段分裂窗算法地表温度LST亮度温度T31和T32热辐射强度I31和I32植被供水指数VSWI几何校正遥感干旱监测指数图图3-1遥感干旱监测计算流程图2. 1. 1遥感数据的预处理为了减小数据量的大小，MODIS把接收的数据转转换为整型数据存储。要获取真实的反射亮度值和辐射亮度值，就需要对获取的MODIS数据进行一个转换。用以下公式可得传感器所获得反射亮度值或辐射亮度值： (3)式中：offsets为偏移量，在HDF文件中，其名称为reflectance_offsets或radiance_offsets；scales为缩放比，在HDF文件中名称为rflectance_scales或radiance_scales。2. 1. 2云检测由于遥感数据在有云的情况下，不能反映地表真实情况。因此这里需要进行云检测。白天判断是否有云条件如下：其中，0.65和0.86分别为MODIS第一、第二波段反射率值，T12为MODIS第32波段亮度温度值。2 .1. 3归一化植被指数（NDVI）的计算植被的生长状况主要与水分有关，水分的供应程度便成了作物生长的关键因素，水分供应充足，植被生长良好，反之生长变差。因此植被的时空变化与土壤水分状况有一定的相关性，因此植被指数可以用于监测对作物生长不利的水分胁迫环境，达到对干旱程度的监测。目前最常用的是归一化植被指数，其定义为：(4)其中，1、2分别表示下垫面的红光波段和近红外波段的反射率，即MODIS数据的第一波段和第二波段的反射率值。2 .1. 4地表温度（LST）的反演方法现有的地表温度反演算法可以分为三大类：单通道算法、分裂窗算法、和多波段算法。其中单通道算法适合于只有一个热红外的数据。而多波段算法还在发展之中，目前没有一个简单可行的算法来进行地表温度反演。而分裂窗算法要使用两个相近的热红外波段遥感数据来进行温度反演。MODIS数据中31和32波段就符合分裂窗算法的要求。因此我们选择分裂窗算法。而在众多的分裂窗算法中，就简单和精度考虑，我们选择Qin et al.两因素地表反演模型。Qin et al.提出的两因素模型是根据星上亮度的线性组合来反演地表温度。其公式如下：(5)Ts为地表温度，T31和T32分别是MODIS第31和32波段的亮度温度（由其DN值计算得到），A0，A1和A2是分裂窗算法的参数，分别定义如下：(6)(7)(8)其中a31，b31，a32，b32是常量，由Qin et al算法可得，在地表温度050范围内，这些常量可以取a31=-64.60363，b31=0.440817，a32=-68.72575，b32=0.473453上述公式的中间参数计算如下： (9)(10)(11) (12) (13) (14)其中i是指MODIS的第31和32波段，分别为i=31或32；是视角为的大气透过率；是波段i的地表比辐射率。根据MODIS图像的DN值可以计算地表热辐射强度，其公式如下：(15)其中，Ii是MODIS第i波段的热辐射强度，NDi是第i波段的DN值，DRi和DRSi分别是第i波段的辐射常量，在MODIS的头文件中。根据普朗克函数可以求解星上亮度温度，计算公式如下： (16)式中：Ti是MODIS第i波段的亮度温度，即式(3)的T31和T32；Ii是MODIS第i波段的热辐射强度，由式(13)得出，i是第i波段的有效中心波长；其值可取31=11.03m、32=12.02m；C1和C2分别是第1是第2光谱常量；大气透过率i()是计算地表温度的基本参数，通常是通过大气水分含量来估计。在MODIS数据中，我们通过第2和19波段来反演大气水分含量，然后通过大气水分含量与大气透过率之间的函数估计大气透过率，其可能的大气水分含量用下式估计：(17)式中：w是大气水分含量；i表示MODIS第i波段的地面反射率。下表3-1为天顶视角为10度的大气透过率估计方程水分含量g/cm2大气透过率估计方程SEER2F0.4-2.0t31(10)=0.99513-0.08082w0.00440.9914804.4t32(10)=0.99377-0.11370w0.00550.99321028.72.0-4.0t31(10)=1.08692-0.12759w0.00250.999211553.0t32(10)=1.07900-0.15925w0.00080.9999173498.34.0-6.0t31(10)=1.07268-0.12571w0.00260.99919921.6t32(10)=0.93821-0.12613w0.00590.99551992.4表3-1 天顶角为10的大气都过率估计方程大气透过率还受遥感视角和大气剖面温度的影响，故要对大气透过率进行视角和温度校正。下表3-2为大气透过率的温度校正函数波段温度校正函数温度区间MODIS 31dt31(T)=0.08T318Kdt31(T)=-0.05+0.00325(T31-278)278T318Kdt31(T)=-0.05T318Kdt32(T)=-0.065+0.004(T32-278)278T318Kdt32(T)=-0.065T278K表3-2 大气透过率的温度校正函数大气透过率的遥感视角校正函数如下：dt31(q)=-0.00247+(2.365210-5)q2(18)dt32(q)=-0.00322+(3.096710-5)q2(19)天顶视角可以用下式简单估计：qVa*|D0Di|(20)其中Va是MODIS卫星高度的星下像元视角，Va=0.0812706；D0是星下像元所在的列号；Di是像元i所在的列号。经过改正后的大气透过率估计方法如下（Qin et al 2001a）：t31(q)=t31(10)+dt31(T)-dt31(q) (21)t32(q)=t32(10)+dt32(T)-dt32(q) (22)从式（11）知，要得出地表温度还需知道地表比辐射率。MODIS图像的地表比辐射率可以用下式估计（Qin et al . 2004）:ei=PvRveiv+(1-Pv)Rseisde(23)式中ei是MODIS图像第i(i=31,32)波段的地表比辐射率；eiv和eis分别是植被和裸土在第i波段的结表比辐射率，分别取e31v=0.98672，e32v=0.98990，e31s= 0.96767，e32s=0.97790。Pv是像元的植被覆盖率，通过植被指数估计(见下面)；de是热辐射相互作用校正，Rv和Rs分别是植被和裸土的辐射比率，建立它们与植被覆盖度之间的关系如下(Qin et al. 2004)：Rv=0.92762+0.07033Pv(24)Rs=0.99782+0.08362Pv(25)其中Pv为归一化植被指数，可由式（2）得出。我们根据Sobrino et al. (2004)的研究提出如下经验公式来估计de：当Pv=0或者Pv=1时，de最小，为de0.0 当0PvPv0.5时，de0.003796(1-Pv)当Pv=0.5时，de最大，为de0.001898(26)至此，由式（5）到式（26）可以算出地表温度。2 .1. 5植被供水指数的计算利用计算得出的植被指数NDVI和地表温度LST可以得出植被供水指数，公式如下：VSWI=NDVI/Ts (27)式中：Ts是地表温度。为了能和气象观测数据综合，这里需要把植被供水指数标准化，公式如下： (28)式中：SDI是标准化植被供水指数，范围为0-100，VSWId是最干旱旱的VSWI，VSWIw是最湿润的VSWI。2 .1. 6几何校正几何校正的处理流程一般包括4 个步骤: 1、控制点的选取和投影变换及输出范围的确定; 2、求输出图像空间到输入图像空间的逆变换函数; 3、逐个像素进行几何位置变换; 4、像素灰度值内插计算。2 .1. 7利用GDAL读取MODIS数据的方法利用软件编程读取MODIS数据是实现上述计算的前提，本项目采用开源库GDAL完成对MODIS数据的提取。这里介绍GDAL库读取MODIS数据的方法。详细代码见附录。GDAL所支持的一些栅格数据格式如IMG、TIFF，在GDAL的抽象数据模型中，每个文件仅有一个数据集。而一个通过HDF文件存储的MODIS数据却包含多个数据集。这些数据集放在抽象模型中的子数据域，而属性信息则存放在相应子数据在元数据中。因此是利用库GDAL提取HDF文件中的波段数据可以分为以下三步。第一步 打开HDF文件和指定的子数据集。HDF文件含有多个子数据集，因此要打开指定的子数据集，可以分为两步：第一步，打开文件，获取子数据名列表。第二步，根据第一步中获得的列表，选择并打开指定的数据集。代码如下：GDALAllRegister(); /注册已有的文件格式驱动const char* filename=”C:/test.hdf”;GDALDataset* fileDataset=(GDALDataset*)GDALOpen(filename, GA_ReadOnly)/以只读方式打开文件char *subList=GDALGetMetadata(GDALDatasetH)fileDataset,”SUBDATASETS”);/获取子数据列表此时获取的子数据集列表中每个列表项以下面形式表示：SUBDATASET_n_NAME=HDF4_SDS:subdataset_type:file_name:subdataset_index其中subdataset_type是预定义的HDF文件名，file_name是输入的文件名，subdataset_index是subdataset序号。访问子数据集时，子数据集的访问名字自能是“=”后面的内容，因此这里需要一个处理去掉前面的部分。代码如下：string subDataset1 = subList0; /以第一个数据集为例subDataset1Name = subDataset1Name.substr(subDataset1Name.find_first_of(“=”)+1);GDALDataset*subDataset=(GDALDataset*)GDALOpen(subDataset1Name.c_str(),GA_ReadOnly); /获取子数据集指针第二步 获取子数据集属性信息。数据集属性信息主要有数据集名字，数据大小，数据集中波段总数，波段数据位移量，波段数据缩放量等。在GDAL的抽象数据模型中，这些属性信息存放在子数据集得元数据项里。下面以读取数据辐射位移量（radiance_scales）为例说明。Const char* metadata=GDALGetMetadataItem(GDALDatasetH)subDataset,“radiance_scales”,NULL);/获取数据集元数据项”radiance_scales”的信息如果上述所获得信息是与波段有关的属性，而一般一个数据集具有多个波段，则上面所得的结果是所有波段的相应信息的集合，并且每个波段的属性之间以，隔开，这时，为以后使用，还需分离上述所得信息。代码如下：string metadataStr=metadata;vector attributeDataStr;/储存分离后的波段属性string tmpDataStr;/用于储存临时信息for (string:size_type i=0;iGetRasterBand(bandIndex); /bandIndex是波段顺序。GDAL库提供多个读取波段数据的函数，最常用的函数是GDALRasterBand:RasterIO()，这个函数能根据情况自动的执行数据转换，以及对数据窗口放大和缩小。对于波段数据，这个函数横向读取数据的效率远高于纵向读取的效率。因此本文读取波段数据时采用以行为单位读取。具体代码如下：int nXSize=poBand-GetXSize(); int nYSize=poBand-GetYSize();/获取波段数据的长宽double *pafScanline; /定义行缓冲区vector bandData; /定义储存波段数据容器for (int i=0;inYSize;+i)pafScanline = (int *) CPLMalloc(sizeof(double)*nXSize);bandData.push_back(pafScanline);pafScanline=NULL; /申请存放容器for (int j=0;jRasterIO( GF_Read, 0,j, nXSize, 1, bandDataj, nXSize, 1, GDT_Float32, 0, 0 ); /以行为单位读取数据通过以上步骤，可以获得MODIS数据中特定的波段数据，以及其相关的属性数据。完成读取工作后，可以通过GDALClose()函数关闭打开的子数据集和HDF文件。3.2 气象数据处理与降水空间插值法比较3.2.1 气象数据的预处理对于地面实际观测资料，我们现有的只有12个站点的数据，如果直接外推台站所在地及邻近地区的气象数据，精度就难以保证。考虑到遥感图像1000m的空间分辨率，就要求每个点上有相应的各种数据，因此要对点状分布的降雨量数据进行插值计算。数据插值就是：已知一组空间数据，可以是离散点的形式，也可以是分区数据的形式，插值就是从这些数据中找到一个函数关系式，使之最好地逼近已知的空间数据，并根据该关系式推求出区域范围内的其它任意点或任意区域的值。我们这里使用了反距离加权法、径向基函数法和普通克里金法三种方法对已知数据进行插值运算，通过比较最后的插值结果，选出最优的插值方法。首先将我们所得的站点坐标值导入到Excel里，将经纬度坐标都化成以度为单位的值，然后将得出的结果导入到Arcgis 9.3软件中，得出Shape文件，投影采用Albers Conical Equal Area，然后进行插值处理。3. 2. 2各种插值法的比较利用ARCGIS软件，使用反距离加权法、径向基函数法和普通克里金插值法实现对降水数据的插值处理，并分析结果等到最适合本项目的插值方法。a.反距离加权法(IDW)反距离权插值法是基于相似的原理：即两个物体离的近，它们就越相似，反之，离得越远则相似性越小。它以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样本赋予的权重越大。用表示待估点变量值，则有： (29)其中，为点的变量值；为对应权重系数。的计算公式为：, n为已知点的个数；为对于插值点与已知点之间距离的权重函数，其中最常用的一种是：,其中b 为合适常数，一般情况下b取2。反距离加权法的优点是简便易行，缺点是对权重函数的选择十分敏感，且受非均匀分布数据点的影响较大。插值结果如下所示：图3-2反距离加权法降雨量插值结果图b.径向基函数法(RBF)径向基函数插值法如同将一个软膜插入并经过各个已知样点，同时又使表面的总曲率最小。它属于精确插值方法。所谓精确插值方法就是指表面必须经过每一个已知样点。径向基函数包括五种不同的基本函数：平面样条函数，张力样条函数，规则样条函数，高次曲面函数和反高次曲面函数。选择何种基本函数意味着将以何种方式使径向基表面穿过一系列已知样点。我们采用了规则样条函数做径向基函数。插值结果如下所示：图3-3径向基函数法降雨量插值结果图c.普通克里金法(OK)一个待估点变量值的估计值就是其周围影响范围内的n个已知点变量观测值的线性组合，其数学表达式是： (30)其中为区域中位置的目标取值；为区域中点的目标取值。只需求出权重即可，权重可以由以下公式求得： (31)式中为变异函数；为均值。插值结果如下所示：图3-4普通克里金插值法空间降雨量插值结果图d. 插值结果比较与分析为了对不同插值方法的结果进行验证和对比，利用检验点的实测值与计算值之间的误差来计算和比较。误差以平均误差和相对均方根误差（RMSE）作为评价标准。其中平均误差总体反映估计值误差大小，相对均方误差则反映插值的精度。 (34)其中n为数据的个数，为误差值。计算的结果如下表：平均误差相对均方根误差IDW-7.44435.942OK-1.61214.136RBF-2.11618.402表3-4误差计算结果表1）平均误差由表4可知三种方法的平均误差都为负，说明估计值整体低于实测值普通克里金插值法的平均误差优于其它两种插值法，且反距离加权法的平均误差最大。2）相对均方根误差由表4可以看出：普通克里金插值法的相对均方根误差最小，因而此种方法的插值精度最高。其次是径向基函数插值法，精度最差的是反距离加权法。3.2.4 综合降水距平指数我们选用2009年5,6两个月共6旬的降雨量数据，并同时算出这六6旬的多年平均降雨量，运用普通克里金插值法插值，因而在ArcGIS 9.3中一共生成了12副图像，运用栅格计算工具并结合以下公式算出每旬的降雨距平指数。公式如下： (35)其中是降雨距平指数，越大越湿润，是该旬降雨量值，是该旬多年平均降雨量值，当时，取，因而我们可以求出 五幅矢量图。由于考虑到干旱是较长时间缺雨的气象现象，所以一旬以内的缺雨并不一定出现干旱，一旬以上时间的缺雨才会产生干旱问题，所以其评价最好以旬为单位，另外，不仅要考虑当旬，而且还要考虑前期供水量，这里我们考虑了2009年5、6月，总共6旬的降雨影响，从而计算出了综合降雨距平指数，公式如下： (36)是考虑降雨因素的综合降水距平干旱指数，其值范围是0-100，和是当旬的降雨距平指数及其权重。通过比较五降雨距平指数幅图，我们可以得出第三、第四、第五幅图效果较差，因而权重较小，第一、第二、第六幅图效果较好，应使用较大的权重。3. 3 指数耦合遥感监测土壤水分的方法具有宏观、高时效、经济等特点，而且它的监测点点俱到，便于对不同含水量区域面积的统计和分析。当然，由于气象卫星运行特点和遥感传感器本身性能的限制，该方法也存在一定的缺点和问题，例如，它在有云的情况下就不能使用，监测结果的精度也不如常规的广泛。而使用气象数据观测，虽然相对简单，但其性质都是以点代面来监测干旱的程度及范围，这样较难给出土壤不同含水量区域之间的分界线，对于相对较小区域间的土壤含水量的差异是难以反映出来的。遥感监测反映近实时的地表下垫面状态，降水因素反映历史降水量的多少，这两个因素都会影到干旱程度，并且能够优势互补，所以将两者结合起来应用于干旱监测。根据植被供水指数和降水距平指数耦合得出农业旱情指数：DI=B1*SDI+B2*MSRI (37)其中DI为农业旱情指数。0-100表示非常干旱到非常湿润；B1和B2为权重，B1=0.4，B2=0.6。得到农业旱情指数后，还需进行旱情级别的划分。首先进行数据标准化处理，公式如下：X=(DIi-DImin)/(DImax-DImin)*100% (38)式中DImax，DImin是得到DI结果中最大和最小值，X是标准值化后的值。这样处理后，就可以划分干旱等级了。1-5为重旱，5-15为中旱，15-25为轻旱，25-65为正常，65-100为湿润。根据DI给相应观测区域的地图上色，可以得到干旱监测合成图。3. 4 程序设计本系统采用主要采用C+语言编写，主要程序代码见附录1。该程序仅实现了数据运算部分。图像显示和图像制作功能由ESRI的ARCGIS软件实现。程序界面如下：图3-7系统界面图4 实验结果分析数据说明：遥感数据采用NASS公开下载HDF格式的MODIS 1B数据，搭载卫星为Aqua，数据时间范围为2009年8月10日到2009年8月30日，覆盖区域为四川省。气象数据采用分四川省12个降雨量观测站的地理坐标值及2000年到2010年每月3旬的降雨数据。结果图如下所示（1-重旱，2-中旱，3-轻旱，4-正常，5-湿润）：结果表明四川整体降雨量充足，但插值受单点的影响较大。图5-1 四川8月中旬降雨量插值干旱监测图能过较好的反应四川的干旱状态，但局部地区由于受云的影响，没有数据图5-2 四川8月中旬遥感干旱监测图结果能够很好的反应四川干旱状态。图5-3 四川8月中旬旱情监测图（指数耦合）从上面的结果来看，指数耦合的方法能正确反映出旱情分布的基本趋势。这种方法同时克服了气象数据和遥感数据观测的缺点，使得监测结果比较能反应现实旱情状况。5 不足与期望本项目能保证旱情基本趋势能正确监测，但同时还存在很多问题和不足：1. 编写的程序仅能实现数据处理，由于小组能力和知识有限，没有为系统添加数据显示和主题制图模块。这样该系统还需要和很多GIS软件结合使用。同时对降雨量数据要求较高，数据获取较难。2. 系统人机交互不强，因为使用GDAL库读取HDF文件的，数据读取速度不及HDF函数库快。同时系统输出结果单一，并且要实现对干旱的监测，需要大量的数据，不能大范围的推广使用。3. 由于降水数据数量过少，降雨量数据插值结果有可能不是很精确。还有四川地区云雾过多，也很有可能会造成遥感干旱监测数据的缺失。4. 在降雨量插值比较时，我们省略了很多气象上常用的插值方法。这样有可能会造成所选取的插值方法不是最优的。项目周期结束，并不意味这对项目相关知识的学习的结束。在以后的时间里，我们小组成员会逐渐完善干旱监测系统。主要有以下方面：1. GDAL库作为优秀的GIS开源库，其功能绝对不仅仅是读取HDF文件这样简单。在后续的学习中，需要更深层次的去认识GDAL库，学习其支持的主流栅格格式的组织方法。以便给系统添加多数据格式输出的功能，增加系统的灵活性。2. 关注干旱监测学术的动态，研究和改善现在的方法，不断完善系统。3. 学习GIS组件式开发，为系统加上数据显示和专题制图功能。6 项目感想2010年4月的时候，学校第五届SRTP项目申报开始了，几个同学凑一起决定去做点什么，于是就申请了这次项目，由曹云刚老师指导，项目名称是遥感与气象观测数据耦合的地表干旱状态监测。我们小组一共5个人，都是第一次真正接触项目，一开始很兴奋，但是只有一腔热情，完全不知道如何入手。然后我们开始找资料，图书馆，数据库，上网，从中积累了很多找资料的方法。然而，学校的资源有限，网络有时候使用也不是很方便，而且很多资料不知道去哪找，在这里还要感谢曹云刚老师，他给我们提供了很多宝贵的资料和数据，让我们对这个项目很快有了一个大概的认识，知道了要做什么，哪些问题需要解决。这个项目对我们目前所学的知识来说还是很有难度的，很多知识都没有学过，需要现学，同时需要熟练操作ArcGIS软件，而我们课上只学了很基础的，很多需要自学和练习。对任务进行分析之后，我们把小组成员分成了两组，从两个方向进行项目。随着项目的进行，慢慢的发现了我们的瓶颈，我们薄弱的编程能力，以及所需的关于GDAL库的信息的缺乏，所能找到的只有官方网站的资料，而且是全英文的，使用有点辛苦，寻找这方面的资料更是难上加难，我甚至去了国外的数据库找，还是相当的少，几乎没有。编程中遇到了很多的困难，很多很难克服，这里又要感谢曹云刚老师对我们的指导，最终克服了困难。回想这一年，有很多的收获，第一次真正意义上做了一次完整的项目，知道了项目怎么做，遇到问题了该怎么解决，而且小组同学在一起努力，一起解决问题，过程很愉