SRTP結題報告

1、西南交通大學第五期大學生科研訓練計劃（SRTP）遙感與氣象觀測數據耦合的地表干旱狀態(tài)監(jiān)測結題報告2010年4月至2011年4月目錄目錄1 緒論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 不足與期望306 項目感想3_參考文獻321緒論1.1項目背景我國是一個旱災非常嚴重的國家，旱災給農業(yè)、農村和農民造成了

2、巨大的損 失。據統(tǒng)計，1950 1990年間，我國共有11年發(fā)生了重特大干旱，發(fā)生頻次 為26.83%,因干旱造成糧食損失占糧食總產量的 4.02%。而1991 2009年間, 我國共有8年發(fā)生重特大干旱，因干旱造成糧食損失占糧食總產量的6.09%。近年來，我國年年有干旱，平均不到3年發(fā)生一次重特大旱災，尤其經常發(fā)生區(qū) 域性特大旱災。同時我國旱災分布面積廣，過去，我國旱災高發(fā)的區(qū)域主要在干 旱缺水的北方地區(qū)，特別是西北地區(qū)。近幾年，在傳統(tǒng)的北方旱區(qū)旱情加重的同 時，南方和東部多雨區(qū)旱情也在擴展和加重，目前旱災范圍已遍及全國。與此同時，旱災影響范圍已由傳統(tǒng)的農業(yè)擴展到工業(yè)、 城市、生態(tài)等領域，工

3、農業(yè)爭水、 城鄉(xiāng)爭水、超采地下水和擠占生態(tài)用水現象越來越嚴重。因此實現對旱災的監(jiān)測，有利于實施抗旱措施，同時合理分配水資源，節(jié)約用水，對促進農業(yè)生產、保障 糧食安全和區(qū)域可持續(xù)發(fā)展具有重要的現實意義 1。傳統(tǒng)的干旱監(jiān)測采用氣象數據，氣象數據由氣象站觀測得到，通過建立區(qū)域 氣象站網同時輔以水文、社會和經濟等數據實施干旱監(jiān)測。這種方法應用十分廣 泛，在單點精度很高，但其性質是以點帶面來監(jiān)測，這樣很難給出不同干旱區(qū)域 的分界線，同時整體精度也相對較低，難以適應大范圍的干旱監(jiān)測。隨著衛(wèi)星遙感技術的迅速發(fā)展，國內外已在借助遙感手段進行大范圍的干旱 監(jiān)測方面開展了大量的研究和實際應用。 基于遙感手段的干旱

4、監(jiān)測能頻繁和持久 的提供地表特征的面狀信息，具有宏觀，高時效，經濟等特點，且點點俱到，適 應大范圍的干旱監(jiān)測。當然，由于氣象衛(wèi)星運行特點和遙感傳感器本身性能的限 制，該方法也存在一定的缺點和問題，比如目前相應的基于遙感數據的一些反演 模型不夠成熟，并且其受大氣影響較大，在很多有云的情況下，數據不能使用， 其監(jiān)測的結果在小范圍內也不如傳統(tǒng)的氣象監(jiān)測精確。1胡其峰.科學理性對待西南旱災減災需立足長遠.光明日報.2010.4針對單一遙感數據在進行地表干旱狀態(tài)監(jiān)測時的不足， 考慮降水對土壤水分 的階段性影響，我們應該把兩種方法結合在一起考慮，形成地表溫度、植被狀態(tài)、 空間降水等參數耦合的干旱監(jiān)測模型。

5、這樣能提高干旱監(jiān)測的準確性和及時性。1. 2項目介紹項目基于對地觀測MODIS數據與氣象站點實測資料，建立起多源數據耦合 下的地表干旱監(jiān)測方法與相關指標體系，在開源GDAL庫的基礎上，開發(fā)可用于大范圍區(qū)域地表干旱狀態(tài)監(jiān)測的應用系統(tǒng)。項目中涉及四個方面的內容：MODIS預處理算法程序設計、氣象觀測資料處 理、熱紅外地表溫度反演模型程序設計、多元數據耦合模型參數化方案及設計。 其中，MODIS預處理算法程序設計部分包括 HDF格式數據讀取、MODIS幾何校 正與輻射標定處理。氣象觀測資料處理包括氣象觀測站點降水資料讀取、空間插值與柵格化。熱紅外地表溫度反演程序設計主要實現利用分裂窗法的地表溫度反

6、演模型。耦合模型部分則包括了地表溫度參數與降水參數的耦合模型參數計算、 模型程序實現。本項目主要分為三個部分：遙感數據處理，氣象數據處理以及干旱監(jiān)測指數 的耦合。遙感數據處理：利用開源庫 GDAL提取MODIS 1B的數據，在此基礎上利用 MODIS數據中的相應波段反演地表溫度（LST）和歸一化植被指數（NDVI），得 到遙感干旱監(jiān)測因子作物供水指數（VSWI）。氣象數據處理：利用ArcGIS的二次開發(fā)平臺ArcGIS Engine實現氣象數據降 雨量的多種空間插值算法。在比較不同空間插值算法誤差的基礎上， 綜合多年數 據選取最佳的空間插值算法進行降雨量的插值運算， 并計算出氣象干旱監(jiān)測因子

7、綜合降水距平指數（MSRI）。干旱監(jiān)測指數的耦合：建立作物供水指數和綜合降水距平指數的耦合模型， 得到農業(yè)旱情監(jiān)測指數，并劃分干旱等級。項目流程圖如圖1所示:圖1-1項目流程圖1. 3開源庫GDAL與HDF文件格式介紹1.3. 1 GDAL庫介紹在MODIS數據的處理中，我們沒有選用 HDF Group開發(fā)的HDF函數庫， 而是選用了十分流行的開源庫 GDALo因為開源庫GDAL在處理柵格圖像格式有 很大優(yōu)勢，其不僅讀取數據較快，而且能夠把結果轉換成任一其支持的柵格格式。1998 年末，加拿大的 Frank Warmerdam 開始了 GDAL(Geospatial DataAbstracti

8、on Library)項目的編寫工作。該項目得到了許多個人和團隊的支持。GDAL是一個在X/MIT許可協(xié)議下的開源柵格空間數據轉換庫。 它利用抽象數據 模型來表達所支持的各種文件格式。它還有一系列命令行工具來進行數據轉換和 處理。OGR是GDAL項目的一個分支，功能和GDAL類似，只不過他提供對矢量 數據的支持。因此GDAL能同時提供柵格和矢量數據的操作。有很多著名的GIS類產品包括ESRI的ARCGIS Google Earth和跨平臺的GRASS GIS系統(tǒng)都使用了 GDAL/OGR開源庫的若干模塊。利用 GDAL/OGR庫，可以使基于Linux的地理 空間數據管理系統(tǒng)提供對矢量和柵格文件

9、數據的支持。GDAL提供對多種柵格數據的支持，包括 Arc/lnfo ASCII Grid(asc) , GeoTiff (tiff) , Erdas Imagine Images(img) , ASCII DEM(dem)等格式。提供讀取、寫入、 轉換、處理各種柵格數據格式(有些特定的格式對一些操作如寫入等不支持)操 作。GDAL使用抽象數據模型(Abstract Datamodel)來解析它所支持的數據格式， 抽象數據模型包括數據集(Dataset)，坐標系統(tǒng)，仿射地理坐標轉換(Affine GeoTransform)，大地控制點(GCPs),元數據(Metadata)，柵格波段(Rast

10、er Band)， 顏色表(ColorTable)，子數據集域(Subdatasets Domain)，圖像結構域 (lmage_StructureDomain)，XML 域(XML:Domains)。GDAL庫米用ANSI C和C+語言編寫，能通過所有c/c+編譯器的編譯。其 核心類框架如圖2所示圖1-2 GDAL核心類框架其中：GDALMajorObject類：帶有元數據的對象。GDALDataset類：通常是從一個文件中提取的相關柵格波段數據的集合，以及其相關元數據。GDALDriver類：文件格式驅動類。GDALDriverManager類：文件格式驅動管理類。GDALRasterBa

11、 nd類：柵格波段數據類。1.3. 2 HDF文件格式HDF (Hierarchical data format )格式是美國國家高級計算機應用中心為了滿足各種領域研究需求而研制的一種高效存儲和分發(fā)科學數據的層次數據格式，主要用來存儲由不同計算機平臺來產生的各種類型科學數據，便于在不同的計算機平臺上擴展。HDF被設計為：自述性：對于一個HDF文件里的每一個數據對象，有關于該數據的綜合信息。通用性：許多數據類型都可以被嵌入在一個 HDF文件里。靈活性：HDF允許用戶把相關的數據對象組合在一起， 放到一個分層結構中，向數據對象添加描述和標簽。它還允許用戶把科學數據放到多個HDF文件里。擴展性：HD

12、F極易容納將來新增加的數據模式，容易與其他標準格式兼容?？缙脚_性：HDF是一個與平臺無關的文件格式。HDF文件無需任何轉換就可 以在不同平臺上使用。目前HDF有兩個版本，即HDF4和HDF5，兩個版本的數據結構模型變化不 大，但其具體實現方法變化很大，因此這兩個版本是不兼容的。由于GDAL庫采用一個抽象的數據模型來支持遙感柵格圖像， 因此使用GDAL讀取這兩個版本的 HDF數據的方法基本相同。一個 HDF文件包括一個文件頭，至少一個數據描述 塊以及若干個數據元素（圖1-3）。HDF文件像一本帶目錄的書，數據描述塊可 以看成書的目錄，而書的每一章內容則是 HDF文件的數據元素。HDF文件圖1-3

13、 HDF文件格式本項目采用的MODIS 1B數據采用HDF文件格式儲存。圖1-4是通過HDFView軟件瀏覽器查看某MODIS 1km分辨率的文件列表示意圖IranA2D(39162 口3丘5.口口吁.20=02441 了幻&2了mcjfEV 1 KM EmissiveEV 1 KM Emissive Uncei?司EV_2 6 O_Ag g r 1 km_RefS 日波段：2*2040*1354£V_250_Aggr1 kmRefSBIJn cer find exesEV_2SO_Aggr"1 kmRefSB_Sampl彗 Um"EV_500i_Aggr

14、1kin_RefSB > 波段：5*2040*1354EV_5 0 O_Ag g r 1 km_RefS B_LJn ce r t_l ndexas&V_&0 O_Ag g1 fcm_RefSB_Sa mples_U sedHeight ensorZenith衛(wèi)星成像參數:408*271Senso rAzi m uthRangeSolarZenithSolarAzimuthgfl 日“圖1-4 HDF ( MODIS )文件列表示意圖2數據源及GDAL庫的配置2. 1數據源2. 1. 1遙感數據源中等分辨率成像光譜儀 MODIS數據采用HDF(Hierarchical d

15、ata format)數據格式存儲元數據，提供36個波段的地球綜合信息，分布在0.4T4卩m的電磁波 譜范圍內，對開展自然災害與生態(tài)環(huán)境監(jiān)測、 全球環(huán)境和氣候變化研究以及進行 全球變化的綜合性研究等有重要意義。MODIS儀器的地面分辨率為250m、500m和1000m，掃描寬度為2330km。在對地觀測過程中，每秒可同時獲得 6.1兆比 特的來自大氣、海洋和陸地表面信息，每日或每兩日可獲取一次全球觀測數據。本項目研究的數據來源于 NASA公開下載的MODIS 1B數據。本研究中用1000m分辨率的數據，主要在于1000m分辨率的數據包含了 36個波段全部數 據，而且其分辨率大小適中。所獲取的數

16、據是 2009年8月1日到2009年8月 30日白天數據，范圍覆蓋四川省。2. 1.2氣象數據源氣象數據主要有四川省12個降雨量觀測站2000年到2009年每月3旬的氣 象觀測數據，站點分布如下圖所示：圖例.*實驗數據測試數搦8T* 廿0先 W* 出E LUI1 O'O'E封豈 105*107* tf 0"E圖2-1氣象數據站點分布圖2. 1. 3其它數據源四川省省界來源于國家基礎地理信息中心（ Natio nal Geomatics Cen ter ofChi na）,另外還有中國區(qū)域的數字高程模型（DEM,分辨率為1km，來源于NASA）、 投影全部采用埃爾伯斯等

17、積圓錐投影 （Albers Conical Equal Area），這些基礎數 據用于生成旱情監(jiān)測系統(tǒng)的基礎數據庫。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安裝的地

18、方。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:Progr

19、am FilesMicrosoft Visual Studio .NET 2003Vc7binvcvars32.bat" 注冊vc的編譯環(huán)境。依次運行：nmake /f makefile.vcnmake /f makefile.vc installnmake /f makefile.vc devinstall如果編譯不成功，這時有可能是 GDAL代碼自帶的BUG。需要修改gdalfrmtslevellerlevellerdataset.cpp 文件 171 行：“ ?, kPI / 180.0, UNITLABEL DEGREE 將“修改為“” ？保存后重新編譯。e) 編譯好后，把&

20、quot;C:GDAL1.6bin"添加到系統(tǒng)環(huán)境變量中。并在項目預定義內容 中添加以下內容：#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")/編

21、譯時鏈接 GDAL 的 lib 庫這樣在項目中就能使用支持 HDF格式的GDAL庫了。3數據處理原理及實現3. 1遙感數據處理干旱是降水時空分布不均，長時間缺乏降水而形成的氣象現象。這樣就會造 成局部地區(qū)作物生長過程中因供水不足，阻礙作物的正常生長而發(fā)生的水量供應 不足的現象。一般來講，植被的生產狀況主要與水分有關，水分供應程度便成了作物生長 的關鍵因素，水分供應充足，植被生長良好，反之生長變差。植被指數的時空變 化與土壤水分狀況有一定的相關性，因此植被指數可以作為干旱監(jiān)測的因素。植 被指數是衛(wèi)星傳感器不同通道探測數據的線性或非線性組合，能夠反映綠色植物生長和分布的特征指數。本項目采用十分常用

22、的歸一化植被指數NDVI(Normalized Differe nee Vegetation In dex)，其定義為：NDVI=(CH1-CH2)/(CH1+CH2)(1)其中，CH1和CH2分別表示近紅外波段和紅光波段的反射率。地表溫度(Land Surface Temperature, LST)也可用于干旱監(jiān)測。地表溫度是控 制地球表面大多數物理、化學和生物過程的參數之一。影響地表溫度變化的因素 也比較多，比如地表濕度、氣溫、光照強度、地表材質(比如是草坪還是裸露土 地，還是水泥地面，或者是瀝青地面)等。綜合考慮上述的兩個干旱監(jiān)測因子，本項目采用植被供水指數VSWI作為遙感干旱監(jiān)測因子，

23、其定義為：VSWI=NDVI/Ts(2)VSWI的物理意義：當植被供水正常時，衛(wèi)星遙感的植被指數在一定的生長期內保持在一定的范圍，而衛(wèi)星遙感的植被冠層溫度也保持在一定的范圍，如果遇到干旱，植被供水不足，一方面植被的生長受到影響，衛(wèi)星遙感的植被指數將 降低，另一方面當植被受旱脅迫時，為減少水分損失，葉面的氣孔會部分關閉， 導致葉面溫度的升高，從而植被冠層溫度將升高。因此可以用植被供水指數來評 價旱情狀況2。遙感干旱監(jiān)測模塊流程圖如下所示：圖3-1遙感干旱監(jiān)測計算流程圖2. 1. 1遙感數據的預處理為了減小數據量的大小，MODIS把接收的數據轉轉換為整型數據存儲。 要獲 取真實的反射亮度值和輻射亮

24、度值，就需要對獲取的MODIS數據進行一個轉換。用以下公式可得傳感器所獲得反射亮度值或輻射亮度值：錯誤!未找到引用源。(3)式中：offsets為偏移量，在 HDF文件中，其名稱為reflectance_offsets或 radiance_offsets； scales為縮放比，在 HDF文件中名稱為 rflectance_scales或 radia nce_scales。2. 1.2云檢測由于遙感數據在有云的情況下，不能反映地表真實情況。因此這里需要進行 云檢測。白天判斷是否有云條件如下：(0.650.86 0.9) or (T|2 265K) or ( 0.650.86 0 7 and T

25、|2 285K)其中，錯誤!未找到引用源。0.65和錯誤!未找到引用源。0.86分別為MODIS第、第二波段反射率 值，T12為MODIS第32波段亮度溫度值。2 .1. 3歸一化植被指數(NDVI)的計算植被的生長狀況主要與水分有關，水分的供應程度便成了作物生長的關鍵因 素，水分供應充足，植被生長良好，反之生長變差。因此植被的時空變化與土壤 水分狀況有一定的相關性，因此植被指數可以用于監(jiān)測對作物生長不利的水分脅 迫環(huán)境，達到對干旱程度的監(jiān)測。目前最常用的是歸一化植被指數，其定義為：錯誤!未找到引用源。(4)其中，錯誤!未找到引用源。1、錯誤味找到引用源。2分別表示下墊面的紅光波段和近紅外波段

26、 的反射率，即MODIS數據的第一波段和第二波段的反射率值。2 .1.4地表溫度(LST)的反演方法現有的地表溫度反演算法可以分為三大類：單通道算法、分裂窗算法、和多 波段算法。其中單通道算法適合于只有一個熱紅外的數據。而多波段算法還在發(fā)展之中，目前沒有一個簡單可行的算法來進行地表溫度反演。而分裂窗算法要使用兩個相近的熱紅外波段遙感數據來進行溫度反演。MODIS數據中31和32波段就符合分裂窗算法的要求。因此我們選擇分裂窗算法。而在眾多的分裂窗算法 中，就簡單和精度考慮，我們選擇 Qin et al.兩因素地表反演模型。Qin et al.提出的兩因素模型是根據星上亮度的線性組合來反演地表溫度

27、。其公式如下：錯誤!未找到引用源Ts為地表溫度，T31和T32分別是MODIS第31和32波段的亮度溫度(由其DN值計算得到)，Ao，Ai和A是分裂窗算法的參數，分別定義如下：錯誤!未找到引用源。(6)錯誤!未找到引用源。(7)錯誤!未找到引用源。(8)其中a31，b31，a32，b32是常量，由Qin et al算法可得，在地表溫度 050 °C范圍內，這些常量可以取 a31=-64.60363，b31=0.440817 , a32=-68.72575，b32=0.473453上述公式的中間參數計算如下：錯誤!未找到引用源。(9)錯誤!未找到引用源。(10)錯誤!未找到引用源。(1

28、1)錯誤!未找到引用源。(12)錯誤!未找到引用源。(13)錯誤!未找到引用源。(14)其中i是指MODIS的第31和32波段，分別為i=31或32 ;錯誤!未找到引 用源。是視角為錯誤!未找到引用源。的大氣透過率；錯誤!未找到引用源。是波 段i的地表比輻射率。根據MODIS圖像的DN值可以計算地表熱輻射強度，其公式如下:錯誤!未找到引用源(15)其中，li是MODIS第i波段的熱輻射強度，NDi是第i波段的DN值，DRi和DRSi分別是第i波段的輻射常量，在 MODIS的頭文件中。根據普朗克函數可以求解星上亮度溫度，計算公式如下：ln(l _(16)式中：Ti是MODIS第i波段的亮度溫度，

29、即式(3)的T31和T32 ； li是MODIS第i波段的熱輻射強度，由式(13)得出，入i是第i波段的有效中心波長；其值可取入31=11.03卩m、入32=12.02卩m; C1和 C 分別是第1是第2光譜常量；大氣透過率T 9)是計算地表溫度的基本參數，通常是通過大氣水分含量來估 計。在MODIS數據中，我們通過第2和19波段來反演大氣水分含量，然后通 過大氣水分含量與大氣透過率之間的函數估計大氣透過率，其可能的大氣水分含 量用下式估計:(17)錯誤!未找到引用源式中：w是大氣水分含量；p表示MODIS第i波段的地面反射率。F表3-1為天頂視角為10度的大氣透過率估計方程水分含量g/cm2

30、大氣透過率估計方程SEER2F0.4-2.031(10)=0.99513-0.08082w0.00440.9914804.432(10)=0.99377-0.11370w0.00550.99321028.72.0-4.031(10)=1.08692-0.12759w0.00250.999211553.032(10)=1.07900-0.15925w0.00080.9999173498.34.0-6.03i(10)=1.07268-0.12571w0.00260.99919921.632(10)=0.93821-0.12613w0.00590.99551992.4表3-1天頂角為10 °

31、;的大氣都過率估計方程大氣透過率還受遙感視角和大氣剖面溫度的影響， 故要對大氣透過率進行視 角和溫度校正。F表3-2為大氣透過率的溫度校正函數波段溫度校正函數溫度區(qū)間MODIS 3131(T)=0.08T>318K31(T)=-0.05+0.00325(T 31-278)278<T<318K31(T)=-0.05T<278KMODIS3232(T)=0.095T>318K32(T)=-0.065+0.004(T 32-278)278<T<318K32(T)=-0.065T<278K表3-2大氣透過率的溫度校正函數大氣透過率的遙感視角校正函數如下:

32、31( )=-0.00247+(2.365210-5) 2(18)32( )=-0.00322+(3.096710-5) 2(19)天頂視角可以用下式簡單估計：=Va* D0- Di(20)其中Va是MODIS衛(wèi)星高度的星下像元視角，Va=0.0812706 °D0是星下像元所在的列號；Di是像元i所在的列號經過改正后的大氣透過率估計方法如下(Qin et al 2001a )：31()=31(10)+31 (T)-31()(21)32()=32(10)+32(T)-32()(22)從式(11 )知，要得出地表溫度還需知道地表比輻射率。MODIS圖像的地表比輻射率可以用下式估計(Qi

33、n et al . 2004):i = PvRv iv+(1_Pv)Rs is+ d(23)式中i是MODIS圖像第i(i=31,32)波段的地表比輻射率；iv和is分別是植被 和裸土在第i波段的結表比輻射率，分別取 31v=0.98672 , 32v=0.98990 , 31s= 0.96767, 32s=0.97790。 Pv是像元的植被覆蓋率，通過植被指數估計(見下面);d是熱輻射相互作用校正，Rv和Rs分別是植被和裸土的輻射比率，建立它們與 植被覆蓋度之間的關系如下(Qin et al. 2004):Rv=0.92762+0.07033P v(24)Rs=0.99782+0.08362

34、P v(25)其中Pv為歸一化植被指數，可由式(2)得出。我們根據Sobrino et al. (2004)的研究提出如下經驗公式來估計 d :當Pv=0或者Pv=1時，d最小，為d = 0.0當 0<Pv<0.5 時，d = 0.003796Pv當 1>Pv>0.5 時，d = 0.003796(1-Pv)當 Pv=0.5 時，d 最大，為 d = 0.001898(26)至此，由式(5)到式(26)可以算出地表溫度。2 .1. 5植被供水指數的計算利用計算得出的植被指數NDVI和地表溫度LST可以得出植被供水指數，公 式如下:VSWI=NDVI/Ts(27)式中：T

35、s是地表溫度。為了能和氣象觀測數據綜合，這里需要把植被供水指 數標準化，公式如下：錯誤!未找到引用源。(28)式中:SDI是標準化植被供水指數，范圍為0-100,VSWId是最干旱旱的VSWI, VSWIw是最濕潤的 VSWI。2 .1.6幾何校正幾何校正的處理流程一般包括4個步驟：1、控制點的選取和投影變換及輸出范圍的確定；2、求輸出圖像空間到輸入圖像空間的逆變換函數 ；3、逐個像素進行幾何位置變換；4、像素灰度值內插計算。2 .1. 7利用GDAL讀取MODIS數據的方法利用軟件編程讀取MODIS數據是實現上述計算的前提，本項目采用開源庫 GDAL完成對MODIS數據的提取。這里介紹 GD

36、AL庫讀取MODIS數據的方法。 詳細代碼見附錄。GDAL所支持的一些柵格數據格式如IMG、TIFF，在GDAL的抽象數據模型 中，每個文件僅有一個數據集。而一個通過 HDF文件存儲的MODIS數據卻包含 多個數據集。這些數據集放在抽象模型中的子數據域，而屬性信息則存放在相應 子數據在元數據中。因此是利用庫 GDAL提取HDF文件中的波段數據可以分為 以下三步。第一步打開HDF文件和指定的子數據集。HDF文件含有多個子數據集，因此要打開指定的子數據集，可以分為兩步： 第一步，打開文件，獲取子數據名列表。第二步，根據第一步中獲得的列表，選 擇并打開指定的數據集。代碼如下：GDALAllRegis

37、ter(); /注冊已有的文件格式驅動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:subdatase

38、t_index其中subdataset_type是預定義的HDF文件名，file_name是輸入的文件名, subdataset_index 是 subdataset 序號。訪問子數據集時，子數據集的訪問名字自能是=”后面的內容，因此這里需要一個處理去掉前面的部分。代碼如下:string subDataset1 = subList0; / subDataset1Name =以第一個數據集為例subDataset1Name.substr(subDataset1Name.find_first_of(“ =” )+1);GDALDataset*subDataset=(GDALDataset*)GDA

39、L0pen(subDataset1Name.c_str(),GA_ReadOnly); /獲取子數據集指針第二步獲取子數據集屬性信息數據集屬性信息主要有數據集名字，數據大小，數據集中波段總數，波段數 據位移量，波段數據縮放量等。在 GDAL的抽象數據模型中，這些屬性信息存放 在子數據集得元數據項里。下面以讀取數據輻射位移量(radiance_scales)為例 說明。Const char* metadata=GDALGetMetadataltem(GDALDatasetH)subDataset, “ radiance_scales ” ,NULL);/ 獲取數據集元數據項” radiance

40、scales ”的信息如果上述所獲得信息是與波段有關的屬性，而一般一個數據集具有多個波段，則上面所得的結果是所有波段的相應信息的集合，并且每個波段的屬性之間以，隔開，這時，為以后使用，還需分離上述所得信息。代碼如下：string metadataStr=metadata;儲存分離后的波段屬性vector<string> attributeDataStr;/string tmpDataStr;/用于儲存臨時信息for (string:size_type i=O;i<metadataStr.length();+i) if (metadataS tri=' , '

41、)attributeDataStr.push_back(tmpDataStr); /儲存tmpDataStr="”；continue;tmpDataStr=tmpDataSt葉metadataStri;遍歷全部字符，尋找'，'并刪除，并分割字符串第三步讀取波段數據。在獲取子數據集得基礎上，可以通過下面代碼獲取該子數據集里波段數據指 針。GDALRasterBand *poBand=subDataset ->GetRasterBand(bandlndex); /bandindex是波段順序。GDAL庫提供多個讀取波段數據的函數，最常用的函數是GDALRasterB

42、and:RasterlO()，這個函數能根據情況自動的執(zhí)行數據轉換，以及 對數據窗口放大和縮小。對于波段數據，這個函數橫向讀取數據的效率遠高于縱 向讀取的效率。因此本文讀取波段數據時采用以行為單位讀取。具體代碼如下：int nXSize=poBand->GetXSize();獲取波段數據的長寬/定義行緩沖區(qū)定義儲存波段數據容器int nYSize=poBand->GetYSize(); double *pafScanline; vector<double*> bandData; / for (int i=0;i<nYSize;+i)pafScanline = (i

43、nt *) CPLMalloc(sizeof(double)*nXSize); bandData.push_back(pafScanline);申請存放容器pafScanline=NULL; /for (int j=0;j<nYSize;+j) poBand->RasterIO( GF_Read, 0,j, nXSize, 1, bandDataj以行為單位讀取數據,nXSize, 1, GDT Float32, 0, 0 ); /通過以上步驟，可以獲得 MODIS數據中特定的波段數據，以及其相關的屬性數據。完成讀取工作后，可以通過 GDALClose()函數關閉打開的子數據集和H

44、DF文件3.2氣象數據處理與降水空間插值法比較3.2.1氣象數據的預處理對于地面實際觀測資料，我們現有的只有12個站點的數據，如果直接外推臺站所在地及鄰近地區(qū)的氣象數據，精度就難以保證?？紤]到遙感圖像1000m的空間分辨率，就要求每個點上有相應的各種數據，因此要對點狀分布的降雨量 數據進行插值計算。數據插值就是：已知一組空間數據，可以是離散點的形式， 也可以是分區(qū)數據的形式，插值就是從這些數據中找到一個函數關系式，使之最好地逼近已知的空間數據，并根據該關系式推求出區(qū)域范圍內的其它任意點或任 意區(qū)域的值。我們這里使用了反距離加權法、徑向基函數法和普通克里金法三種 方法對已知數據進行插值運算，通過

45、比較最后的插值結果，選出最優(yōu)的插值方法。首先將我們所得的站點坐標值導入到 Excel里，將經緯度坐標都化成以度為 單位的值，然后將得出的結果導入到 Arcgis 9.3軟件中，得出Shape文件，投影 采用Albers Conical Equal Area，然后進行插值處理。3. 2. 2各種插值法的比較利用ARCGIS件，使用反距離加權法、徑向基函數法和普通克里金插值法實現對降水數據的插值處理，并分析結果等到最適合本項目的插值方法。a. 反距離加權法(IDW)反距離權插值法是基于相似的原理：即兩個物體離的近，它們就越相似，反之，離得越遠則相似性越小。它以插值點與樣本點間的距離為權重進行加權平

46、均,(29)d bej離插值點越近的樣本賦予的權重越大。用e表示待估點變量值，則有：其中，jj1, ,n為點Xj,yj的變量值；j為對應權重系數。j的計j算公式為：f ( de j )nj)，n為已知點的個數；f(dej)為對于插值點(Xe,ye)與ej jj 1已知點xj，yj之間距離dej的權重函數，其中最常用的一種是:f (dej)其中b為合適常數，一般情況下b取21反距離加權法的優(yōu)點是簡便易行，缺點是對權重函數的選擇十分敏感，且受 非均勻分布數據點的影響較大。插值結果如下所示：圖3-2反距離加權法降雨量插值結果圖b. 徑向基函數法（RBF）徑向基函數插值法如同將一個軟膜插入并經過各個已

47、知樣點，同時又使表面的總曲率最小。它屬于精確插值方法。所謂精確插值方法就是指表面必須經過每 一個已知樣點。徑向基函數包括五種不同的基本函數：平面樣條函數，張力樣條函數，規(guī)則樣條函數，高次曲面函數和反高次曲面函數。 選擇何種基本函數意味著將以何種方式使徑向基表面穿過一系列已知樣點。我們采用了規(guī)則樣條函數做徑向基函數插值結果如下所示:J-J<_二】-壬 s ?s二 k一 主 二r I(30)圖3-3徑向基函數法降雨量插值結果圖c. 普通克里金法（0K）一個待估點變量值的估計值就是其周圍影響范圍內的n個已知點變量觀測值的線性組合，其數學表達式是:n(0)i ( i)i 1其中（J為區(qū)域 中位置

48、的目標取值；（0）為區(qū)域中點0的目標取值 只需求出權重即可，權重可以由以下公式求得：nj 1nj 1j1(i 1,2,n)廠(i),( j)-( i),( o)(31)式中為一變異函數；為均值。插值結果如下所示：W' 護CTE W 6' o'K 10: o' o*E L034lfl5: ii'o'k 1 o?: o' err:圖3-4普通克里金插值法空間降雨量插值結果圖d. 插值結果比較與分析為了對不同插值方法的結果進行驗證和對比， 利用檢驗點的實測值與計算值 之間的誤差來計算和比較。誤差以平均誤差和相對均方根誤差 （RMS）作為評價 標

49、準。其中平均誤差總體反映估計值誤差大小，相對均方誤差則反映插值的精度。100(35)RMSE 丄1n(34)其中n為數據的個數， 為誤差值計算的結果如下表:平均誤差相對均方根誤差IDW-7.44435.942OK-1.61214.136RBF-2.11618.402表3-4誤差計算結果表1）平均誤差由表4可知三種方法的平均誤差都為負，說明估計值整體低于 實測值普通克里金插值法的平均誤差優(yōu)于其它兩種插值法，且反距離加權法的平 均誤差最大。2）相對均方根誤差由表4可以看出：普通克里金插值法的相對均方根誤差 最小，因而此種方法的插值精度最高。 其次是徑向基函數插值法，精度最差的是 反距離加權法。3.

50、2.4 綜合降水距平指數我們選用2009年5,6兩個月共6旬的降雨量數據，并同時算出這六 6旬的 多年平均降雨量，運用普通克里金插值法插值，因而在ArcGIS 9.3中一共生成了 12副圖像，運用柵格計算工具并結合以下公式算出每旬的降雨距平指數。公 式如下：SRI2Rw其中SRI是降雨距平指數，SRI越大越濕潤，R是該旬降雨量值，Rw是該 旬多年平均降雨量值，當R 2Rw時，取SRI 100，因而我們可以求出SRI0, SRI1, SRI2, SRI3, SRI4, SRI5 五幅矢量圖。由于考慮到干旱是較長時間缺雨的氣象現象， 所以一旬以內的缺雨并不一定 出現干旱，一旬以上時間的缺雨才會產生

51、干旱問題，所以其評價最好以旬為單位， 另外，不僅要考慮當旬，而且還要考慮前期供水量，這里我們考慮了2009年5、6月，總共6旬的降雨影響，從而計算出了綜合降雨距平指數，公式如下：MSRI A0* SRI0 A1* SRI1 A5* SRI5 (36)MSRI是考慮降雨因素的綜合降水距平干旱指數，其值范圍是0-100，SRI0和A0是當旬的降雨距平指數及其權重。通過比較五降雨距平指數幅圖，我們可 以得出第三、第四、第五幅圖效果較差，因而權重較小，第一、第二、第六幅圖 效果較好，應使用較大的權重。3. 3指數耦合遙感監(jiān)測土壤水分的方法具有宏觀、高時效、經濟等特點，而且它的監(jiān)測點點俱到，便于對不同含

52、水量區(qū)域面積的統(tǒng)計和分析。當然，由于氣象衛(wèi)星運行特點和遙感傳感器本身性能的限制， 該方法也存在一定的缺點和問題，例如，它在有云的情況下就不能使用，監(jiān)測結果的精度也不如常規(guī)的廣泛。而使用氣象數據 觀測，雖然相對簡單，但其性質都是以點代面來監(jiān)測干旱的程度及范圍，這樣較難給出土壤不同含水量區(qū)域之間的分界線， 對于相對較小區(qū)域間的土壤含水量的 差異是難以反映出來的。遙感監(jiān)測反映近實時的地表下墊面狀態(tài)，降水因素反映歷史降水量的多少， 這兩個因素都會影到干旱程度，并且能夠優(yōu)勢互補，所以將兩者結合起來應用于 干旱監(jiān)測。根據植被供水指數和降水距平指數耦合得出農業(yè)旱情指數:DI=B1*SDI+B2*MSRI(3

53、7)其中DI為農業(yè)旱情指數。0-100表示非常干旱到非常濕潤；Bi和B2為權重, Bi=0.4, B2=0.6。得到農業(yè)旱情指數后，還需進行旱情級別的劃分。首先進行數據標準化處理， 公式如下：X=(Dli-Dlmin)/(DI max-Dlmin )*100%(38)式中DImax，DImin是得到DI結果中最大和最小值，X是標準值化后的值。這 樣處理后，就可以劃分干旱等級了。 1-5為重旱，5-15為中旱，15-25為輕旱， 25-65為正常，65-100為濕潤。根據DI給相應觀測區(qū)域的地圖上色，可以得到干旱監(jiān)測合成圖。3. 4程序設計本系統(tǒng)采用主要采用C+語言編寫，主要程序代碼見附錄1。該

54、程序僅實現 了數據運算部分。圖像顯示和圖像制作功能由 ESRI的ARCGIS軟件實現。程序 界面如下：圖3-7系統(tǒng)界面圖4實驗結果分析數據說明：遙感數據采用NASS公開下載HDF格式的MODIS 1B數據，搭載 衛(wèi)星為Aqua，數據時間范圍為2009年8月10日到2009年8月30日，覆蓋 區(qū)域為四川省。氣象數據采用分四川省12個降雨量觀測站的地理坐標值及 2000 年到2010年每月3旬的降雨數據。結果圖如下所示（1-重旱，2-中旱，3-輕旱，4-正常，5-濕潤）：四川百月中有降雨K S值干旱監(jiān)測團I1-圖5-1四川8月中旬降雨量插值干旱監(jiān)測圖冋II齊冃中自屋肢干旱監(jiān)測圖i- n三三J 數據

55、圖5-2四川8月中旬遙感干旱監(jiān)測圖能過較好的 反應四川的 干旱狀態(tài)，但局部地區(qū) 由于受云的 影響，沒有四川8月中旬旱情監(jiān)測團兩指數輻合、結果能夠很好的反應四川干旱狀態(tài)。圖5-3四川8月中旬旱情監(jiān)測圖（指數耦合）從上面的結果來看，指數耦合的方法能正確反映出旱情分布的基本趨勢。這 種方法同時克服了氣象數據和遙感數據觀測的缺點， 使得監(jiān)測結果比較能反應現 實旱情狀況。5不足與期望本項目能保證旱情基本趨勢能正確監(jiān)測，但同時還存在很多問題和不足：1. 編寫的程序僅能實現數據處理，由于小組能力和知識有限，沒有為系統(tǒng)添 加數據顯示和主題制圖模塊。這樣該系統(tǒng)還需要和很多GIS軟件結合使用。同時 對降雨量數據要求較高，數據獲取較難。2. 系統(tǒng)人機交互不強，因為使用 GDAL庫讀取HDF文件的，數據讀取速度 不及HDF函數庫快。同時系統(tǒng)輸出結果單一，并且要實現對干旱的監(jiān)測，需要 大量的數據，不能大范圍的推廣使用。3. 由于降水數據數量過少，降雨量數據插值結果有可能不是很精確。還有四 川地區(qū)云霧過多，也很有可能會造成遙感干旱監(jiān)測數據的缺失。4. 在降雨量插值比較時，我們省略了很多氣象上常用的插值方法。 這樣有可 能會造成所選取的插值方法不是最優(yōu)的。項目周期結束，并不意味這對項目相關知識的學習的結束。 在以后的