植被水分遙感監(jiān)測(cè)處理系統(tǒng)的建立及應(yīng)用

(中國科學(xué)院 遙感應(yīng)用研究所 北京 100101)

摘 要：在分析了遙感監(jiān)測(cè)植被水分可行性的基礎(chǔ)上，闡述了以MODIS HDF為數(shù)據(jù)源進(jìn)行植被水分遙感監(jiān)測(cè)的處理流程；根據(jù)處理流程，利用IDL結(jié)合ENVI中提供的庫函數(shù)開發(fā)出了植被水分遙感監(jiān)測(cè)處理系統(tǒng)。介紹了系統(tǒng)實(shí)現(xiàn)中的關(guān)鍵技術(shù)和核心功能實(shí)現(xiàn)，系統(tǒng)在黃淮海農(nóng)作物水分監(jiān)測(cè)中得到了成功的應(yīng)用，為農(nóng)作物水分監(jiān)測(cè)的業(yè)務(wù)化處理發(fā)揮了核心作用。

關(guān)鍵詞：植被水分；遙感；接口定義語言；ENVI

中圖分類號(hào)：TP391文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-3695(2009)05-1819-03

Establishment and application of vegetation water monitoring system

XIONG Jinguo，WANG Litao，WANG Shixin，ZHOU Yi

(Institute of Remote Sensing Application Chinese Academy of Sciences Beijing 100101 China)

Abstract:Based on feasibility analysis of monitoring vegetation water using remote sensing，stated a flow chart of monitoring vegetation water with MODIS HDF data.According to this，developed vegetation water monitoring system utilizing IDL 6.0 combined with the ENVI routine library.Demonstrated some key techniques in the process of implementing the system.The system has been successfully applied in monitoring crop water content of Huanghuai region.

Key words:vegetation water;remote sensing;IDL(interactive data langage);ENVI

0 引言

植被水是植被體內(nèi)部含有的游離態(tài)或化合態(tài)的水，是植物進(jìn)行光合作用的主要原料之一，對(duì)植被的生長具有重要意義。通過對(duì)植被水的監(jiān)測(cè)可以反映出植被體的生長狀況。另外，植被水研究對(duì)于森林或草原火災(zāi)、旱災(zāi)以及生態(tài)安全監(jiān)測(cè)等方面也有著重大的意義。對(duì)植被水的估算有多種野外實(shí)測(cè)的手段，但是野外實(shí)測(cè)獲取的植被水信息只反映了采樣點(diǎn)周邊較小范圍內(nèi)的情況，并且只是量測(cè)當(dāng)時(shí)的較短時(shí)間內(nèi)的狀況。野外實(shí)測(cè)要獲取較大范圍內(nèi)的植被水動(dòng)態(tài)情況，就需要布控非常多的采樣點(diǎn)，需要長時(shí)間連續(xù)的野外作業(yè)，這些在經(jīng)濟(jì)上和技術(shù)上實(shí)現(xiàn)起來代價(jià)都是很高的。

遙感可以很好地滿足植被水監(jiān)測(cè)的空間廣泛性和時(shí)間連續(xù)性的要求。不同分辨率的遙感影像可以提供不同區(qū)域范圍上的影像信息，對(duì)于大范圍動(dòng)態(tài)的植被水反演監(jiān)測(cè)工作而言，利用遙感手段是一個(gè)很好的選擇。

遙感監(jiān)測(cè)植被水分涉及到遙感數(shù)據(jù)的獲取、預(yù)處理和植被水分反演模型計(jì)算以及水分圖發(fā)布等一系列環(huán)節(jié)，處理過程復(fù)雜，如果僅用通用的遙感圖像處理軟件工作量較大，耗用時(shí)間長。本文利用IDL結(jié)合ENVI中提供的一些功能函數(shù)，根據(jù)遙感監(jiān)測(cè)農(nóng)作物水分業(yè)務(wù)化運(yùn)行的需要，在總結(jié)、優(yōu)化和集成前人研究成果的基礎(chǔ)上，開發(fā)出了農(nóng)作物水分遙感監(jiān)測(cè)的應(yīng)用軟件，實(shí)現(xiàn)用遙感監(jiān)測(cè)農(nóng)作物水分的快速處理，為農(nóng)作物水分監(jiān)測(cè)的業(yè)務(wù)化運(yùn)行奠定了良好的基礎(chǔ)。

1 農(nóng)作物水分監(jiān)測(cè)的業(yè)務(wù)流程

1.1 農(nóng)作物水分遙感監(jiān)測(cè)所用數(shù)據(jù)源

本文采用的數(shù)據(jù)源是MODIS傳感器接收到的HDF格式的數(shù)據(jù)。MODIS是當(dāng)前世界上新一代圖譜合一的光學(xué)遙感儀器，具有36個(gè)光譜通道，分布在0.4~14 μm的電磁波譜內(nèi)。MODIS儀器的地面分辨率分別為250、500和1 000 m。MODIS有如下幾個(gè)特點(diǎn)^[1]：a)涉及的光譜范圍廣(36個(gè)波段)、譜帶窄，可見光—近紅外通道除659和2 100 nm外，譜帶寬度為10~35 nm，有很多大氣糾正的特征波段，便于大氣參數(shù)的反演。b)空間分辨率ch1、ch2為250 m；ch3~ch7為500 m。對(duì)于農(nóng)作物水分監(jiān)測(cè)通過將前兩個(gè)波段重采樣成500 m分辨率，以方便進(jìn)行波段運(yùn)算，構(gòu)造水分反演模型。c)TERRA和AQUA衛(wèi)星都是太陽同步極軌衛(wèi)星，TERRA在地方時(shí)上午過境，AQUA在地方時(shí)下午過境。TERRA與AQUA上的MODIS數(shù)據(jù)在時(shí)間更新頻率上相配合，加上晚間過境數(shù)據(jù)，對(duì)于接收MODIS數(shù)據(jù)來說，可以得到每天最少兩次白天和兩次黑夜更新數(shù)據(jù)。這樣的數(shù)據(jù)更新頻率，對(duì)實(shí)時(shí)地球觀測(cè)和應(yīng)急處理有較大的實(shí)用價(jià)值。d)NASA對(duì)MODIS數(shù)據(jù)實(shí)行全球免費(fèi)接收的政策，提供了大范圍的農(nóng)作物水分監(jiān)測(cè)數(shù)據(jù)上的支持。

1.2 HDF 1B/1A數(shù)據(jù)的預(yù)處理

本系統(tǒng)采用的數(shù)據(jù)是從http://ladsweb.nascom.nasa.gov/data/search.html 下載的Lever 1B產(chǎn)品，以及中國科學(xué)院遙感應(yīng)用所接收站下載的沒有進(jìn)行定標(biāo)的Lever 1A數(shù)據(jù)。針對(duì)不同級(jí)別的數(shù)據(jù)要進(jìn)行的處理有所不同，但一般都包括以下幾個(gè)步驟：反射率定標(biāo)、幾何糾正、投影轉(zhuǎn)換、云水掩模處理。對(duì)于要進(jìn)行反射率計(jì)算的Lever 1A級(jí)數(shù)據(jù)，R=reflectance_scaleB×(DN-reflectance_offsetB)/cos(A)。其中：DN為圖像像素記錄下的亮度值；reflectance_scaleB為校正偏移系數(shù)；reflectance _offsetB為校正偏移系數(shù)；A為太陽天頂角。它們的值可在相應(yīng)HDF數(shù)據(jù)集中的屬性中找到。幾何糾正可以利用HDF內(nèi)建的經(jīng)緯度波段信息，根據(jù)圖像的分辨率來采集控制點(diǎn)，然后利用控制點(diǎn)對(duì)MODIS 500 m的七個(gè)波段進(jìn)行幾何糾正。在進(jìn)行農(nóng)作物水分監(jiān)測(cè)時(shí)，需要將投影轉(zhuǎn)換成適合中國區(qū)域特點(diǎn)的Albers圓錐等面積投影，其投影參數(shù)是中央經(jīng)線為105°，標(biāo)準(zhǔn)平行緯度分別為25°、47°。對(duì)于MODIS圖像來說，容易受天氣的影響，需要掩模處理有云的區(qū)域，通常的做法是將紅波段反射率值大于0.18的認(rèn)為是云，水體的判斷以NDVI 值小于0為判斷依據(jù)^[2]。NDVI=(Rnir-Rred)/(Rnir+ Rred)，即近紅波段（NIR）和紅波段（RED）的值的差除以兩者之間的和。其中NIR對(duì)應(yīng)MODIS的第二波段，RED對(duì)應(yīng)MODIS的第一波段。水的反射波譜曲線異于土壤和植被的波譜曲線，其中在NIR處的反射率小于RED處的反射率，這就導(dǎo)致了計(jì)算出的水體NDVI值為負(fù)值。由于水體的該特征較為明顯，可以作出判斷：所有NDVI值小于0的區(qū)域都認(rèn)為是水體。預(yù)處理流程如圖1所示。

1.3 農(nóng)作物水分監(jiān)測(cè)指標(biāo)的選取

在農(nóng)作物水分遙感監(jiān)測(cè)的應(yīng)用中，通常采用等效水層厚度(equivalent water thickness，EWT)來表征植被含水量的多少。EWT被定義為農(nóng)作物水分含量(鮮重與干重的差值)與葉面面積的比值^[3]，即EWT=(FW-DW)/A。其中，F(xiàn)W、DW、A分別為植被的鮮重、干重和葉面面積，均可通過野外測(cè)量得到。

在葉片光譜特征中，可見光譜段受葉片葉綠素含量的控制、近紅外譜段受葉內(nèi)細(xì)胞結(jié)構(gòu)的控制、短波紅外譜段受葉細(xì)胞內(nèi)水分含量的控制。可見光中綠色波段0.52~0.59 μm對(duì)區(qū)分植被類型敏感，紅光波段0.63~0.69 μm對(duì)植被覆蓋度、植物生長狀況敏感。對(duì)于復(fù)雜的農(nóng)作物遙感，僅用個(gè)別波段或多個(gè)單波段數(shù)據(jù)分析對(duì)比來提取農(nóng)作物信息是相當(dāng)局限的，因而往往選用多光譜遙感數(shù)據(jù)經(jīng)分析運(yùn)算（加、減、乘、除等線性或非線性組合方式）產(chǎn)生某些對(duì)農(nóng)作物長勢(shì)、生物量等有一定指示意義的數(shù)值，即所謂的植被指數(shù)^[4]。針對(duì)植被的光譜曲線和MODIS通道的波段設(shè)置，很多研究者建立了對(duì)葉片水分反應(yīng)敏感的植被指數(shù)，并建立這些植被指數(shù)與葉片水分間的經(jīng)驗(yàn)關(guān)系，從而達(dá)到反演農(nóng)作物水分反演的目的。

常用的植被指數(shù)有水分脅迫指數(shù)（moisture stress index，MSI）^[5]、歸一化差異水分指數(shù)（normalized difference water index，NDWI）^[6]、全球植被水分指數(shù)（global vegetation moisture index，GVMI）^[7]、短波紅外水分脅迫指數(shù)（shortwave infrared water stress index，SIWSI）^[8]、地表含水量指數(shù)（surface water capacity index，SWCI）^[9]、短波角度歸一化指數(shù)（shortwave angle normalized index，SANI）^[10]等。通過建立采樣點(diǎn)處植被指數(shù)與實(shí)測(cè)含水量之間的回歸曲線得到植被水分的回歸方程，然后利用回歸方程即可反演出圖像覆蓋區(qū)的植被水分。

1.4 農(nóng)作物水分遙感監(jiān)測(cè)流程

通過以上分析可知，要進(jìn)行農(nóng)作物水分監(jiān)測(cè)需要進(jìn)行模型構(gòu)建、模型比較、模型選擇以及影像數(shù)據(jù)處理等幾個(gè)步驟。主要工作流程如圖2所示。

2 農(nóng)作物水分監(jiān)測(cè)處理系統(tǒng)的總體設(shè)計(jì)

本文選擇具有高級(jí)圖像處理能力的IDL，結(jié)合ENVI遙感圖像處理軟件中提供的函數(shù)庫來進(jìn)行系統(tǒng)的開發(fā)。IDL作為面向矩陣、語法簡單的第四代可視化語言，是進(jìn)行科學(xué)數(shù)據(jù)的可視化和分析以及跨平臺(tái)應(yīng)用的最佳選擇。IDL面向矩陣，可以靈活方便地分析任何有格式或無格式的數(shù)據(jù)類型，支持通用文本及圖像數(shù)據(jù)，并且支持在NASA、TPT、NOAA等機(jī)構(gòu)中大量使用的HDF、CDF及netCDF等科學(xué)數(shù)據(jù)格式。ENVI作為美國RSI公司的旗艦產(chǎn)品，它是由遙感領(lǐng)域的科學(xué)家采用交互式數(shù)據(jù)語言IDL開發(fā)的一套功能強(qiáng)大的遙感圖像處理軟件，是處理分析并顯示多光譜數(shù)據(jù)、高光譜數(shù)據(jù)和雷達(dá)數(shù)據(jù)的高級(jí)工具。

根據(jù)圖2的業(yè)務(wù)流程，可以將系統(tǒng)抽象為三個(gè)主要的模塊：a)HDF數(shù)據(jù)讀入顯示模塊，包括數(shù)據(jù)集屬性和HDF全局屬性的讀入以及數(shù)據(jù)集本身的顯示；HDF與TIFF和IMG格式的相互轉(zhuǎn)換。b)HDF圖像處理模塊，主要包括輻射糾正、幾何糾正、投影轉(zhuǎn)換、圖像鑲嵌、掩模處理、圖像裁剪。c)模型計(jì)算模塊，根據(jù)反演得到的模型對(duì)預(yù)處理后的HDF文件進(jìn)行植被水分的反演。系統(tǒng)總體結(jié)構(gòu)如圖3所示。

該系統(tǒng)的主界面如圖4所示。

3 系統(tǒng)核心功能實(shí)現(xiàn)及其關(guān)鍵技術(shù)

3.1 IDL對(duì)HDF文件的讀寫處理

IDL中有很多函數(shù)可以直接讀取HDF數(shù)據(jù)信息。一般對(duì)HDF文件中一個(gè)數(shù)據(jù)集文件的讀取步驟如下：用函數(shù)hdf_sd_start打開文件，返回文件標(biāo)志符hdfid，然后調(diào)用函數(shù)hdf_sd_nametoindex獲取某一名稱的數(shù)據(jù)集索引。對(duì)MODIS Lever 1B數(shù)據(jù)，不同分辨率數(shù)據(jù)對(duì)應(yīng)的數(shù)據(jù)集名稱分別為“EV_250_RefSB”“EV_500_RefSB”“EV_1KM_RefSB”。得到索引后，調(diào)用函數(shù)hdf_sd_select返回該變量的標(biāo)志符；然后再用此標(biāo)志符作為參數(shù)，通過函數(shù)hdf_sd_getdata得到變量的內(nèi)容；最后用函數(shù)hdf_sd_endaccess關(guān)閉變量通道，用hdf_sd_end關(guān)閉打開的文件。

sd_id=hdf_sd_start(hdf_filename，/read)

dataset_index=hdf_sd_nametoindex(sd_id dataset_name)

dataset_id=hdf_sd_select(sd_id，dataset_index)

hdf_sd_getdata，dataset_id，data /noreverse

hdf_sd_endaccess，dataset_id

hdf_sd_end，sd_id

3.2 MODIS IB輻射定標(biāo)計(jì)算

將DN值轉(zhuǎn)換為反射率定標(biāo)涉及到讀取HDF數(shù)據(jù)集中的屬性信息，主要是reflectance_scales、reflectance_offsets和solarZenith三個(gè)屬性的讀取。讀取屬性信息除了打開數(shù)據(jù)集等通用操作外，還要用到函數(shù)hdf_sd_attrfind獲得屬性的索引，再通過hdf_sd_attrinfo讀取屬性數(shù)據(jù)，通過太陽天頂角、反射率比例和位移參數(shù)對(duì)數(shù)據(jù)集波段進(jìn)行計(jì)算。波段計(jì)算可以利用ENVI 中的math_doit函數(shù)進(jìn)行，此函數(shù)的關(guān)鍵參數(shù)是波段計(jì)算表達(dá)式的確定^[11]。模塊運(yùn)行界面如圖5所示。

3.3 圖像云水掩模處理

云水掩模處理的原理在數(shù)據(jù)預(yù)處理部分已經(jīng)闡述。通常的云檢測(cè)是通過對(duì)紅波段進(jìn)行閾值分割處理，水體判斷是通過NDVI指數(shù)設(shè)置某一閾值進(jìn)行分割。在具體實(shí)現(xiàn)中用envi_get_data函數(shù)獲得特定波段數(shù)組，用make_array數(shù)組產(chǎn)生一個(gè)與該波段數(shù)組維數(shù)相同的初值為1的數(shù)組；根據(jù)設(shè)定的閾值檢索出滿足條件的像元并將數(shù)組的值設(shè)為0；然后用該二值數(shù)組對(duì)HDF中的波段進(jìn)行掩模處理，具體可以用ENVI中的envi_mask_apply_doit進(jìn)行計(jì)算。主要代碼如下：

t_fid=[fid，fid]

pos1=[0，1]

exp=′(b2-b1)/(b2+b1)′

bnames=[′歸一化植被指數(shù)′]

envi_doit，′math_doit′ $

fid=t_fid，pos=pos1，dims=dims，$

exp=exp，r_fid_r=ndvi_fid，out_bname=bnames，out_name=ndvi

ndvi=envi_get_data(fid=ndvi_fid，dims=dims，pos=0)

maskarray=make_array(ns，nl，/float，value=1)

maskarray[where(ndvi LT 0)]=0

envi_enter_data，maskarray，bnames=bnames，def_stretch=def_stretch，file_type=file_type，r_fid=fid_water

envi_mask_apply_doit，fid=fid，pos=pos，dims=dims，m_fid=fid_water，m_pos=m_pos，value=value，r_fid=watermasked_fid，out_name=out_name，in_memory=0

envi_file_mng，fid=fid_water，/remove

3.4 用shape文件對(duì)圖像進(jìn)行裁剪

可以用某一區(qū)域的行政邊界矢量文件對(duì)柵格圖像進(jìn)行裁剪，實(shí)現(xiàn)過程關(guān)鍵是對(duì)矢量文件進(jìn)行柵格化處理。首先讀取shape文件中的每一個(gè)封閉區(qū)域要素的坐標(biāo)，通過envi_convert_file_coordinates將地理坐標(biāo)轉(zhuǎn)換為圖像的文件坐標(biāo)，即數(shù)組的行列號(hào)，通過這個(gè)要素邊界數(shù)組用函數(shù)envi_create_roi、envi_define_roi建立一個(gè)感興趣區(qū)，再用envi_get_roi獲得感興趣區(qū)內(nèi)所有像元的文件坐標(biāo)，從而建立掩模二值矩陣。用此掩模矩陣對(duì)原始柵格圖像進(jìn)行掩模處理，就完成了裁剪。

3.5 圖像鑲嵌功能實(shí)現(xiàn)

圖像鑲嵌的功能能夠?qū)Χ喾鵐ODIS圖像進(jìn)行拼接，從而滿足大區(qū)域范圍的監(jiān)測(cè)。在ENVI函數(shù)庫中提供了MOSAIC_DOIT函數(shù)，需要確定的參數(shù)有鑲嵌后圖像的范圍，各個(gè)圖像左上角X0、Y0的值以及圖像的投影信息map_info。首先用envi_file_query查詢每個(gè)圖像的維數(shù)，可以得到每幅圖像的四個(gè)角點(diǎn)的文件坐標(biāo)，通過函數(shù)envi_convert_file_coordinates將四個(gè)角點(diǎn)坐標(biāo)轉(zhuǎn)換為地理坐標(biāo)；然后通過循環(huán)比較得到所有圖像中x坐標(biāo)和y坐標(biāo)的最小值和最大值以及最小值和最大值的差值deltaX和deltaY；再用deltaX、deltaY除以xsize和ysize后就可以得到圖像的大小。鑲嵌后圖像的投影map_info可通過以下幾個(gè)函數(shù)確定：

proj=envi_get_projection(fid=fids[0])

map_info = envi_map_info_create(proj=proj mc=[0，0，west，north] ps=out_ps)

temp = bytarr(10，10)

envi_enter_data temp map_info=map_info /no_realize r_fid=tmp_fid

在得到鑲嵌后具有地圖投影信息圖像的tmp_fid后，就可以將每幅圖像的左上角坐標(biāo)的地理坐標(biāo)轉(zhuǎn)換為新的文件坐標(biāo)，也即X0[i]、Y0[i]的值。在確定了這些關(guān)鍵參數(shù)后，用mosaic_doit函數(shù)進(jìn)行圖像的鑲嵌。

3.6 觀測(cè)點(diǎn)反射率值批提取

在進(jìn)行植被水分模型的構(gòu)建時(shí)，是利用野外觀測(cè)臺(tái)點(diǎn)的實(shí)測(cè)值與圖像中對(duì)應(yīng)像元的反射率進(jìn)行統(tǒng)計(jì)回歸而得來的。本功能實(shí)現(xiàn)了快速將野外諸多樣點(diǎn)對(duì)應(yīng)像元在每個(gè)波段的反射率值的快速提取并保存成文件。具體實(shí)現(xiàn)如下：首先將觀測(cè)點(diǎn)的經(jīng)緯度坐標(biāo)保存成指定格式的文本文件，通過IDL提供的OPENR，READF函數(shù)將坐標(biāo)讀入到內(nèi)存數(shù)組中；循環(huán)執(zhí)行數(shù)組中的每一個(gè)樣點(diǎn)，用函數(shù)envi_convert_file_coordinates，fid，xf，yf，valuelong，valuelat算出該地理坐標(biāo)對(duì)應(yīng)的文件坐標(biāo)xf、yf 再用函數(shù)data=envi_get_data(fid=fid，dims=dims，pos=i)取出每一個(gè)波段對(duì)應(yīng)的數(shù)據(jù)數(shù)組，通過數(shù)組索引xf和yf就可以得到該點(diǎn)處的反射率值。