基于IDL的MODIS影像地表蒸散發(fā)參數(shù)反演系統(tǒng)

鄧世贊,張友靜,2,張子衡,謝麗軍,王文種

(1.河海大學(xué)水文水資源學(xué)院,江蘇南京 210098;2.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室,江蘇南京 210098)

地表蒸散發(fā)參數(shù)是涉及多個圈層水熱交換與平衡的重要參數(shù).利用氣象數(shù)據(jù),采用彭曼 蒙特斯模型、彭曼 組合模型和基于太陽輻射的日蒸散發(fā)模型等,計算點尺度的蒸散發(fā),已有許多研究[1-2],但由于這些研究局限于離散的點觀測與估算,存在插值外延精度低、大范圍高密度觀測成本大等缺陷.隨著多時相、多分辨率衛(wèi)星遙感數(shù)據(jù)的應(yīng)用,通過遙感信息提取來估算地表蒸散發(fā)方法的應(yīng)用越來越廣.

在地表蒸散發(fā)遙感反演方法中,SEBAL(地表能量平衡)模型[3-4]是目前應(yīng)用較多的模型之一.但蒸散發(fā)計算中所需要的地表參數(shù)眾多,流程復(fù)雜,需要遙感專業(yè)人員進行重復(fù)性的演算,不利于長時間序列的蒸散發(fā)反演.Wang J等[5]改進了SEBAL模型,并利用C++語言對農(nóng)作物——山核桃的消費性用水量進行了計算,但模型參數(shù)局限于產(chǎn)品數(shù)據(jù),無法從原始影像中反演得到.吳炳方等[6]利用能量平衡余項式方法和Penman-Monteith模型相結(jié)合的方法開發(fā)了區(qū)域蒸散發(fā)遙感監(jiān)測系統(tǒng),該系統(tǒng)利用逐日氣象數(shù)據(jù)與遙感反演參數(shù),能動態(tài)反映區(qū)域蒸散發(fā)的時空變化規(guī)律,但其簡化了地表與大氣由于溫度差引起的熱量傳輸過程,使反演的蒸散發(fā)量具有一定的不確定性.

針對以上問題,本文利用MODIS遙感數(shù)據(jù)并結(jié)合地面氣象觀測數(shù)據(jù),應(yīng)用SEBAL模型,通過ENVI/IDL二次開發(fā)語言實現(xiàn)了日蒸散發(fā)量的批處理,并給出了該方法的應(yīng)用實例.

1 SEBAL(地表能量平衡)模型

SEBAL模型的優(yōu)點是可充分利用由遙感圖像反演得到的地表參數(shù),結(jié)合少量的常規(guī)氣象資料就能得到區(qū)域地表蒸散發(fā)量.該模型建立的基礎(chǔ)是能量平衡方程[7],即

式中:R n——地表凈輻射通量;G——土壤熱通量,即下墊面土壤中的熱交換量;H——下墊面到大氣的顯熱通量(也稱感熱通量);λ——水的汽化潛熱;E——蒸散發(fā)量;λ E——下墊面到大氣的潛熱通量.計算得到各通量后,再利用能量平衡方程就可以反演日蒸散發(fā)量.

1.1 地表凈輻射通量

地表凈輻射通量是地表接收到的太陽輻射和大氣長波輻射減去地表面反射的太陽輻射和發(fā)射的長波輻射后得到的輻射差值.它是地表能量、動量、水分輸送與交換過程中的主要能量來源.地表輻射平衡方程為[8]

式中:Q——太陽總輻射;a——反照率;εα——地表反射率;εs——地表比輻射率;σ——斯蒂芬-波爾茲曼常數(shù);Tα——空氣溫度;Ts——地表溫度.

1.2 土壤熱通量

土壤熱通量是指用于土壤熱交換的那部分能量,表征土壤表層和深層的熱量傳遞.對于植被覆蓋區(qū)域,土壤熱通量與地表凈輻射通量的比值為[9]

式中b為歸一化植被指數(shù).

不同地物類型,土壤熱交換不同:非植被-裸土,G=0.2Rn;水體,G=0.9Rn-40.

1.3 顯熱通量

顯熱通量表征下墊面與大氣間湍流形式的熱交換,又稱感熱通量,其表達式為[10]

式中:ρ——空氣密度;Cp——空氣定壓質(zhì)量熱容;ra——空氣動力學(xué)阻抗.

1.4 日蒸散發(fā)量

由于遙感所獲得的是瞬時蒸散發(fā)量,需推算到日蒸散發(fā)量.本文采用謝賢群的正弦公式[11],考慮到1 d中不可避免地會有云的出現(xiàn),而氣象站觀測到的日照時數(shù)可以反映1d中云的出現(xiàn)時長,因而在進行時間尺度擴展的過程中,加入日照時數(shù)來進行校正[12].

式中:a,b——系數(shù),一般分別取0.25和0.50;n——氣象站觀測日照時數(shù);NE——理論日照時數(shù),可以通過地理緯度計算得出;t——影像接收時間;Et——瞬時蒸散發(fā)量;Ed——日蒸散發(fā)量.

2 日蒸散發(fā)量算法的IDL實現(xiàn)

2.1 氣象數(shù)據(jù)插值

SEBAL模型蒸散發(fā)反演過程中部分地表參數(shù)的計算必須依靠氣象站實測數(shù)據(jù),而由于氣象站數(shù)據(jù)往往是離散的,必須通過插值方法獲取該區(qū)域任意點的數(shù)據(jù).本系統(tǒng)通過IDL語言自動交互讀取氣象站矢量文件,利用griddata函數(shù)中反距離插值方法實現(xiàn)氣象數(shù)據(jù)的插值.

2.2 地表參數(shù)反演

SEBAL模型涉及地表參數(shù)較多,其中大氣透過率、地表反照率、地表比輻射率等大多數(shù)參數(shù)通過經(jīng)驗?zāi)Ｐ颓蟪鯷13].本文以地表溫度為例,采用劈窗算法,簡述其IDL的實現(xiàn)過程.算法的實現(xiàn)主要利用地表輻射率和大氣透過率2個因子.在計算大氣透過率時不僅要利用MODIS的B31和B32波段對其進行溫度校正,還要結(jié)合MODIS原始影像中Sensor Zenith數(shù)據(jù)集對傳感器進行視角糾正,具體實現(xiàn)流程見圖1.

為提高系統(tǒng)運行效率,將眾地表參數(shù)合成一個文件,每個地表參數(shù)相當于其中1個波段,需要該數(shù)據(jù)時可直接讀取該波段數(shù)據(jù).對本模塊而言,用戶可以通過界面選擇不同的參數(shù)和計算方法,并可以方便地按不同渲染方式輸出地表參數(shù)影像,如圖2所示.

圖1 地表溫度反演算法流程Fig.1 Flow chart of inversion of surface temperature

2.3 SEBAL模型各通量IDL的實現(xiàn)

圖2 地表參數(shù)查詢窗口Fig.2 Interface for querying surface parameters

2.3.1 凈輻射通量

考慮到地形起伏導(dǎo)致的陽坡和陰坡太陽輻射不同,將對太陽輻射進行地形校正.本系統(tǒng)考慮了DEM及所反映的地形特征參數(shù)如坡度、坡向的太陽天頂角校正,校正時只需輸入影像獲取時間、DEM、經(jīng)緯度波段等初始條件.

2.3.2 土壤熱通量

計算土壤熱通量時,將地表分為水體、植被和裸土3類,并利用IDL提供的where判別函數(shù)來判斷地物類型.計算某類型地表土壤熱通量時,首先計算歸一化植被指數(shù)NDVI(以INDV表示),然后分別查找INDV＞0.25-INDV＜0和 INDV=0～0.25之間所包含的像元并獲取各個像元所在的位置,最后再計算這些像元的土壤熱通量.裸土地表土壤熱通量IDL實現(xiàn)算法為

2.3.3 顯熱通量

顯熱通量的計算最為復(fù)雜.SEBAL模型引入了“干點”和“濕點”概念,利用“干點”處蒸散發(fā)量基本為0,而“濕點”處蒸散發(fā)量達到最大的假設(shè),結(jié)合近地表面溫度差的線性函數(shù),并考慮大氣穩(wěn)定度,通過Monin-Obukhov方法[14]求得不同狀態(tài)條件下的摩擦風(fēng)速、動力學(xué)粗糙度和阻抗等參數(shù)[15],從而計算獲得顯熱通量H.本文通過IDL語言并以顯熱通量的變化小于1%作為循環(huán)終止條件,實現(xiàn)了整個循環(huán)迭代過程.

在程序?qū)崿F(xiàn)過程中,“干點”和“濕點”根據(jù)地表溫度參數(shù)來選取(圖3).為保證“干點”和“濕點”選取的合理性,應(yīng)考慮研究區(qū)局部云層覆蓋導(dǎo)致的地表溫度降低而引起的誤差.利用云層的高反射低溫度特點,通過多波段合成方法將云層檢測出來,獲取其所在像元位置,并進行掩膜處理.

3 SEBAL模型日蒸散發(fā)量反演

黃河三花間(三門峽—花園口區(qū)間)是我國中部地區(qū)主要的半干旱區(qū),該系統(tǒng)的研究區(qū)包括陜西、山西、河南3省的部分地區(qū).該區(qū)域地表類型復(fù)雜,春季易干旱,夏季多暴雨.本研究采用的遙感影像是Terra衛(wèi)星的2002年3月8日、4月18日和6月12日3期MODIS影像.

3.1 日蒸散發(fā)量IDL批處理結(jié)果

充分利用ENVI/IDL二次開發(fā)函數(shù),實現(xiàn)了日蒸散發(fā)量反演的批處理.運算過程中,僅需要人工選擇預(yù)處理后的MODIS影像、對應(yīng)日期的氣象shapefile文件以及DEM數(shù)據(jù)和初始條件,經(jīng)系統(tǒng)運算就可得到日蒸散發(fā)量分布,如圖4所示.

圖3 顯熱通量計算顯示窗口Fig.3 Interface for calculating sensib le heat flux

圖4 黃河三花間日蒸散發(fā)量分布(單位:mm)Fig.4 Distribution of daily evapotranspiration in Sanmenxia-Huayuankou watershed of Yellow River(unit:mm)

分析圖4可知,日蒸散發(fā)量高值區(qū)分布在研究區(qū)西南熊耳山地帶以及北部植被覆蓋度較高的山區(qū)及有灌溉的農(nóng)業(yè)區(qū),一般在4mm以上,而植被覆蓋相對較少的裸土和城鎮(zhèn),蒸發(fā)量相對較小,一般在1.0～2.5mm之間.不同下墊面具有不同的地表特征和水熱性狀,有植被覆蓋區(qū)明顯大于周圍無植被區(qū)和少植被區(qū).日蒸散發(fā)量的空間分布基本符合實際.

圖5 日蒸散發(fā)量與各地表參數(shù)相關(guān)分析顯示窗口(2002-04-18)Fig.5 Interface for corrlation between daily evapotranspiration and surface parameters(April 18,2002)

3.2 日蒸散發(fā)量與地表參數(shù)相關(guān)分析

IDL語言在數(shù)據(jù)處理、統(tǒng)計分析與圖形顯示方面具有強大的優(yōu)勢,可以使開發(fā)者通過較少的命令完成大量的數(shù)據(jù)預(yù)處理、變換及統(tǒng)計分析等工作.通過對日蒸散發(fā)量與地表參數(shù)的統(tǒng)計計算和相關(guān)分析(圖5),可以知道各參數(shù)對蒸散發(fā)量的影響程度.

分析結(jié)果表明:日蒸散發(fā)量與地表溫度相關(guān)性最好,整體呈線性負相關(guān)關(guān)系,相關(guān)系數(shù)R達到0.95,這是由于地表溫度高的區(qū)域多為植被覆蓋相對比較少的地區(qū);植被覆蓋度高的地區(qū)由于植被蒸騰作用造成冠層溫度降低;日蒸散發(fā)量與地表溫度差和植被指數(shù)的相關(guān)系數(shù)分別為0.63和0.61.

4 參數(shù)計算結(jié)果存取

為了將日蒸散發(fā)量反演過程中計算獲得的地表參數(shù)和其他參數(shù)進行影像匹配,可以將其結(jié)果保存為ENVI標準格式或指定格式,并定義好地理坐標投影.應(yīng)用時,通過自定義map-info函數(shù)讀取校正后的影像坐標并賦予參數(shù)影像就可以實現(xiàn)了.具體代碼如下:

5 結(jié) 語

本文基于SEBAL模型,在IDL環(huán)境下編程建立了地表蒸散發(fā)反演系統(tǒng).該系統(tǒng)界面友好,操作簡單,運算速度快,可進行批量、快速的日蒸散發(fā)量計算,統(tǒng)計分析功能可用于精度評價、蒸散發(fā)與地表參數(shù)回歸分析等.由于該系統(tǒng)具有平臺無關(guān)性特點,因而可將其作為獨立模塊加入水文分析等專業(yè)系統(tǒng),從而可有效解決通用遙感軟件缺乏專業(yè)信息提取功能的問題.

[1]PETER W.A discussion on and alternative to the Penman-Monteith equation[J].Agricultural Water M anagement,2009,96(4):711-721.

[2]CHAUHAN S,SHRIVASTAVA R K.Performance evaluation of reference evapotranspiration estimation using climate based methods and artificial neural networks[J].Water Resources Management,2009,23(5):825-837.

[3]TIMMER MANS W J,KUSTAS W P,ANDERSON M C,et al.An intercomparison of the surface energy balance algorithm for land(SEBAL)and the two-source energy balance(TSEB)modeling schemes[J].Remote Sensing of Environment,2007,108(4):369-384.

[4]BASTIAANSSEN W G M.SEBAL-based sensible and latent heat fluxes in the irrigated Gediz Basin[J].Journal of Hydrology,2000,229(1/2):87-100.

[5]WANG J,SAMMIS T W,MEIER C A,et al.A modified SEBAL model for spatially estimating pecan consum ptive water use for Las Cruces,New Mexico[C]//American Meteorological Society 15th Conference on Applied Climatology.Savannah,GA:[s.n.],2005:paper number:7.13.

[6]吳炳方,熊雋,閆娜娜,等.基于遙感的區(qū)域蒸散量監(jiān)測方法:ETWatch[J].水科學(xué)進展,2008,19(5):671-678.(WU Bing-fang,XIONG Jun,YAN Na-na,et al.ETWatch for monitoring regional evapotranspiration with remote sensing[J].Advances in Water Science,2008,19(5):671-678.(in Chinese))

[7]楊永民,馮兆東,周劍.基于SEBS模型的黑河流域蒸散發(fā)[J].蘭州大學(xué)學(xué)報:自然科學(xué)版,2008,44(5):1-3.(YANG Yongmin,FENG Zhao-dong,ZHOU Jian.Evapotranspiration in heihe river basin based on SEBS model[J].Journal of Lanzhou University:Natural Science,2008,44(5):1-3.(in Chinese))

[8]趙英時.遙感應(yīng)用分析原理與方法[M].北京:科學(xué)出版社,2003:432-436.

[9]李守波,趙傳燕.基于能量平衡的關(guān)川河流域蒸散發(fā)的遙感反演[J].遙感技術(shù)與應(yīng)用,2006,21(6):521-525.(LI Shou-bo,ZHAO Chuan-yan.Estimating evapotranspiration based on energy balance in guanchuan river basin using remote sensing[J].Remote Sensing Technology and App lication,2006,21(6):521-525.(in Chinese))

[10]韓惠.基于遙感技術(shù)的祖厲河流域土地利用/土地覆蓋變化與蒸散發(fā)研究[D].蘭州:蘭州大學(xué),2006.

[11]謝賢群.遙感瞬時作物表面溫度估算農(nóng)田全日蒸散總量[J].遙感學(xué)報,1991,6(4):253-260.(XIEXian-qun.Estimation of daily evapo-transpiration from one time-of-day remotely sensedcanopy temperature[J].Journal of Remote Sensing,1991,6(4):253-260.(in Chinese))

[12]姜紅.基于MODIS影像的新疆奇臺縣區(qū)域蒸散發(fā)量的研究[D].烏魯木齊:新疆大學(xué),2007.

[13]鮑平勇.半干旱區(qū)域日蒸散發(fā)估算的遙感研究[D].南京:河海大學(xué),2007.

[14]KOLOSKOV G,MUKHAMEJANOV K H,TANTON T W.Monin-Obukhov length as a cornerstone of the SEBAL calculations of evapotranspiration[J].Journal of Hydrology,2007,335(1/2):170-179.

[15]唐榮林.基于SEBAL模型的海河流域蒸散發(fā)的時空反演研究[D].武漢:武漢大學(xué),2007.