基于植被指數－地表溫度特征空間的伊河流域蒸散發量估算

周 峰，王文，王曉剛

（1.南京大學地理與海洋科學學院，江蘇 南京 210093；2.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210098）

0 引言

蒸散發（Evapotranspiration）包括土壤蒸發和植被蒸騰，是土壤－植被－大氣系統中水分和能量傳輸過程的重要組成部分。傳統的蒸散發量估算和測量方法（如波文比－能量平衡法、空氣動力學、互補相關理論和渦動相關法等）經過長時間的發展，形成了較為成熟的理論，并在單站點取得了較高的模擬精度［1］，但由單站點推算區域蒸發蒸騰量時，受到下墊面非均勻性的影響，一般很難取得準確結果。多時相、多光譜及具有傾斜角度的遙感資料，能客觀地反映下墊面的幾何結構和濕、熱狀況，從而使得區域蒸散發量的快速、準確估算成為可能。

基于遙感技術的地表蒸散發估算方法可分為經驗方法［2］、蒸發比方法［3］、能量平衡余項法［4，5］等［6，7］。其中，蒸發比方法主要通過確定潛熱通量和地表可利用能量的比值（即蒸發比），進而對蒸散發量進行估算，該方法簡單易行，同時降低了機理模型中對數據精度的要求，因而在實際中得到了廣泛應用［8，9］，其關鍵在于如何求解蒸發比。Price發現當研究區域的植被覆蓋度和土壤濕度條件變化范圍較大時，以遙感獲得的地表溫度和植被指數為橫縱坐標得到的散點圖呈近似三角形分布［10］，該特征空間反映了地表的干濕狀況，可用于蒸發比的估算，該方法也稱為三角形法［11］。三角形法只需常規的氣象資料和遙感資料，可操作性強，本文吸收該方法的優點，在植被指數平滑去噪處理和地表溫度高程訂正的基礎上，通過增加植被蒸騰受到的水分脅迫因子，提出改進的三角形算法，對黃河重要支流伊河流域的蒸散發量進行估算研究。

1 研究方法

1.1 三角形方法

三角形方法通過構建植被指數和地表溫度（VI－Ts）的特征空間，確定任意地表干濕狀況下的Priestlay－Taylor系數a值，求得蒸發比EF，并最終基于Priestlay－Taylor方程，利用地表可利用能量（Rn－G），進行潛熱通量（λE）的求解和實際蒸散發量（λE）的換算［8］。

式中：a為Priestlay－Taylor常數（1.26），Δ為飽和水汽壓與溫度關系曲線的斜率，γ為干濕表常數，Rn為地表凈輻射，G為土壤熱通量。

如圖1所示的三角形特征空間ABC中，A點表示完全裸露充分供水的區域，B點表示完全裸露極端干燥的區域，C點表示植被完全覆蓋且水分充分供應的地區。在濕邊AC上所有能量用于蒸發，因而a有最大值amax；在干邊BC上，a有最小值amin，且在完全裸土（B）處amin＝0，植被完全覆蓋（C）處amin＝1.26，并假定在干邊上amin隨植被覆蓋度（fveg）呈線性變化。

圖1 植被指數－地表溫度三角形特征空間Fig.1 The concept of a triangular VI－Ts diagram

在濕邊AC上，系數a具有最大值amax＝1.26，這樣利用不同植被指數類型下的amax和amin即可對空間任意坐標點（VI，Ts）下的a值（ai）進行線性插值：

結合式（1），蒸發比EF可表示為：

式中：Tsmax、Tsmin分別為各植被指數所對應的地表溫度的最大值和最小值。

1.2 方法的改進

如式（3）所示，三角形方法中EF由兩部分組成，即土壤部分的蒸發比（EFsoil）和植被部分的蒸發比（EFveg）。根據Priestlay－Taylor公式的定義，參數a取1.26表示“無平流條件下，凈輻射能量保持不變，非飽和陸面水分得到充分供應時的蒸散發量”，因而可以看出傳統三角形算法中植被蒸散發量處于潛在蒸發水平，并沒有考慮到蒸騰過程中的氣象因子和水分脅迫作用的影響。本文參考Nishida等提出的 MOD16算法［12］，利用 Penman－Monteith公式、Priestlay－Taylor公式并結合互補相關理論對式（3）中EFveg進行改進：

從而得到改進后的蒸發比計算公式：

式中：rc為冠層阻抗，ra為空氣動力學阻抗。入射的太陽輻射、水氣壓差、根區的土壤含水量、冠層的氣溫等因素影響植被冠層阻抗和空氣動力學阻抗，從而影響植被的蒸騰作用，與改進前（式3）的蒸發比計算相比，改進后的計算方法（式5）考慮了植被蒸騰受到的水分脅迫作用及氣象因子的影響。

2 參數的計算

2.1 地表通量

地表凈輻射Rn是地表能量、動量、水分傳輸和交換過程的主要來源，表現為地表輻射能量收支的差額，其計算公式如下：

式中：a為地表反照率，Rs↓為入射的短波輻射，Rl↑為地面向上的長波輻射項，Rl↓為大氣向下的長波輻射項，εs為地表比輻射率，σ為Stefan－Boltzman常數，εa為大氣比輻射率，Ts和Ta分別為衛星過境時間的瞬時地表溫度和大氣溫度。Rs為太陽短波輻射，本文采用Zillman（1972）晴空條件下的計算方法對其進行計算［13］，同時考慮研究區內地形的坡度和坡向對太陽輻射的影響，在其計算過程中對入射的太陽高度角進行了地形糾正［14］；大氣向下的長波輻射（Rl↓）及向上的長波輻射（Rl↑）通過Stefan－Boltzman公式進行計算；土壤熱通量（G）則采用SEBS模型參數化方案進行估算［15，16］。

2.2 蒸發比

利用增強型植被指數（EVI）和地表溫度構建特征空間，從而對EF進行估算，其中干濕邊確定采用最大最小值擬合的方法［8］，即：每隔0.01取不同EVI類型下對應的地表溫度的最大值和最小值，并進行最小二乘法線性擬合。改進后算法（式5）的關鍵在于冠層阻抗和空氣動力學阻抗參數的求解，其中冠層阻抗是葉面單個氣孔阻抗的并聯，參考Jarvis（1976）的經驗模式進行估算［17］：

式中：f1（Ta）、f2（PAR）、f3（VPD）、f4（Ψ）和 f5（CO2）分別為大氣溫度、光合有效輻射、飽和水汽壓差、葉水勢及CO2濃度影響因子，rcmin為冠層最小阻抗，rcuticle為與葉面表皮相關的冠層最大阻抗。大氣溫度和太陽有效輻射對冠層阻抗影響因子的計算參考文獻［12］，飽和水汽壓差計算參考文獻［18］。最小冠層阻抗的計算參考文獻［19］進行季節性訂正：

式中：LAI為葉面積指數，rcsto為最小氣孔阻抗，Sun確定冬小麥的最小氣孔阻抗為100s／m，不同植被最小氣孔阻抗可通過查詢表確定［20］。利用Thom模型對空氣動力學阻抗（ra）進行計算［21］，日蒸散發量則采用正弦函數法進行估算［9］。

3 研究區概況及數據處理

3.1 研究區及數據源

伊河流域發源于伏牛山北麓河南省欒川縣，流域面積6 029km2，全長268km，年平均降雨量為500～1 100mm，且年內分配不均。流域內的地形復雜，氣候垂直地帶性較明顯。

使用的數據源包括：1）衛星遙感數據，是搭載于EOS／Terra衛星上的中分辨率成像光譜儀（MODIS）獲取的數據，星下點空間分辨率為250～1 000m，原始數據經加工與處理后，共形成44個標準產品并免費向全球提供下載。產品數據經過了大氣輻射校正和幾何校正，可直接使用。本研究選擇2003－2004年20期無云或少云天氣條件下的MODIS產品數據作為主要信息源，包括地表反照率數據（MCD43A3，分辨率500m）、植被指數數據（MOD13A1，分辨率250m）、地表溫度和比輻射率（MOD11A1，分辨率1 000m）、土地覆蓋數據（MCD12Q1，分辨率500m）、地理定標數據（MOD03，分辨率1 000m），并利用 MRT軟件對其進行投影轉換、數據集提取及邊界裁剪。2）數字地面高程數據（SRTM，分辨率90m），用于生成流域邊界、地表溫度的地形訂正等。3）氣象數據，包括日最高（最低）氣溫、平均相對濕度、日照時數、平均風速等，空間尺度拓展上，利用氣溫、水汽壓、風速與高程的經驗關系將氣象站觀測值推算到整個研究區域［18］。4）欒川水文站E601蒸發皿日蒸發量數據，各種數據統一重采樣到500m空間分辨率。

3.2 數據的預處理

3.2.1 植被指數的平滑除噪 MODIS植被指數數據集雖然通過最大值合成處理，降低了云的影響，但受大氣及傳感器觀測視角等因素的影響，數據集仍存在大量的噪聲，作為特征空間構造中的核心變量，有必要對其進行重建。

Hird對多種除噪方法進行了對比分析［22］，發現4253H雙濾波法［23］具有較好的處理效果。本文采用該濾波法對EVI數據集（MOD13A1）進行預處理，該除噪過程不僅改善了EVI序列數據的異常高值和低值部分，并能保持數據的完整性。

3.2.2 地表溫度的地形訂正和分析 分析20期的地表溫度和高程數據，兩者存在明顯的負相關關系，而VI－Ts特征空間構建的原理，是基于水分蒸發和植被蒸騰對地表溫度的降低作用（蒸散發過程吸收熱量），因此在構造特征空間前，必須進行高程訂正，以去除高程對地表溫度的影響，使得訂正后的特征空間更真實地反映地表的干濕狀況。本文通過以下兩步對地表溫度進行訂正［24］：

式中：R 為干空氣氣體常數（287Jkg－1K－1），Cp 為大氣的比熱容（1 004Jkg－1K－1），Po為水平面平均大氣壓（101.3kpa），Z 為高程。

4 結果分析

4.1 模擬結果與參考作物蒸發量和蒸發皿數據對比

FAO－56Penman－Monteith方法計算的參照作物蒸散發量類似于表面開闊、高度一致、生長旺盛、完全覆蓋地表而不缺水的綠色草地的蒸散量［25］，而遙感反演的蒸散發考慮500m分辨率網格內平均地表狀態條件下的蒸散發，兩者相對比，一方面可以反映陸面地表不同下墊面的綜合蒸發能力和干濕狀況，同時也可用來對遙感估算的蒸散發量進行合理性檢驗。

圖2 遙感估算值與參考作物蒸發量的對比分析Fig.2 Comparison of the evapotranspiration estimation using RS and FAO－56 at Luanchuan Station

對比站點參考作物蒸散發量（FAO－56）和遙感估算的蒸散發量（MOT，圖2），兩者具有較高的相關性（r＝0.86）；用前30天降雨量（P30）反映季節性干濕狀況，對遙感估算結果做進一步分析，MOT方法反演的地表蒸散發量均小于FAO－56的結果（圖2a），且濕潤季節兩者差距較小，兩者間的比值（MOT／FAO－56）與地表的干濕狀況相一致（圖2b）。

利用蒸發皿觀測的水面蒸發量與遙感反演的蒸散發進行對比，用于反映流域不同季節的蒸發狀況（圖3）。結果表明，遙感估算的蒸散發和水面蒸發的對比呈現出季節性變化特征，夏季日平均氣溫（Ta＿mean）較高時兩者差別（Pan－MOT）較大，因為夏天地表可利用能量較多，蒸發過程受土壤水分影響較大，而冬天的蒸散發則主要受氣象條件的制約。水面蒸發反映了一定區域特定時段內蒸發潛能，其可視為實際蒸散發的上限，如果估算的陸面蒸散發高于蒸發皿蒸發，則結果不合理，從兩者的對比可知，遙感估算的蒸散發量多低于蒸發皿觀測值，受大氣和地表的快速波動以及云量的影響，個別日期（DOY，2004286）出現了高于水面蒸發的情形。

圖3 遙感估算值與蒸發皿觀測值的對比分析Fig.3 Comparison of the evapotranspiration estimation using RS and Panobservation data

4.2 不同地表覆蓋類型下蒸發量的差異分析

基于研究區2004年土地利用數據，選取具有代表性的4個典型月份的蒸散發估算結果，對比分析不同植被類型的差異性。圖4為改進的三角形算法在不同地表覆蓋類型下的均值和方差統計值，水體與其它地物類別相比具有最高的蒸發量和最小的內部方差；受氣象條件和土壤含水量空間差異性的影響，植被內部具有較大的方差，并在植被生長季節（DOY：2003104，2004253）更為明顯。

研究區內80%以上為農田和林地，夏季林地具有較高的蒸散發量（DOY：2004253），由于流域內種植冬小麥，林地主要以落葉林為主，冬春季節農田表現為較高的蒸散發量（DOY：2004343，2003104），其差異性主要與植被的生長狀況有關，灌叢和草地所占比重較少，夾雜在耕地和林地之間，其蒸發量也多介于兩者之間。由于研究區城鎮所占面積較小，分布較零散，夾雜著植被和水體等具有較高蒸散量的地物類型，因而呈現較高的蒸散發量。

很多研究者分析了不同下墊面條件下區域蒸散發量的差異。例如，郭曉寅對黑河流域的蒸散發空間分布特征進行研究時發現，水分條件較好的沼澤和人工綠洲蒸散量最高，其次是中覆蓋度的草地和林地，其它地區隨著植被覆蓋度的降低，蒸散發量也逐漸減少［26］；劉朝順等的研究同樣表明受水分狀況和植被覆蓋度的影響，高蒸發量主要分布在水體、沼澤地區，而裸地及低植被覆蓋度的自然植被區蒸發量較?。?7］。而有關城鎮蒸發特征的研究中，何延波等利用中分辨率MODIS數據（250～1 000m）分析發現，城鎮區具有較高的蒸發量［27］，但劉朝順等利用高分辨率TM影像進行分析時卻發現，城鎮區域表現出較低的蒸發量，因此蒸散發量的差異除了與地表覆蓋類型有關外，可能還與使用的數據源有關。MODIS中低分辨率尺度上，城鎮為水體、植被和裸地的混合體，實際的蒸發狀況多與三者之間的比例有關?？傊?，不同地表覆蓋類型下蒸散發差異與植被的季節性生長狀況有關，主要植被類型之間的對比分析與研究區的實際狀況相符。

圖4 典型日不同地表類型下的日蒸散發量特征Fig.4 Selected daily evapotranspiration estimation for different land cover

5 結論

本文通過構建植被指數和地表溫度特征空間對伊河流域蒸散發量進行估算，考慮到植被蒸騰過程受水熱環境因素的影響，引入植被冠層阻抗和空氣動力學阻抗對三角形算法進行改進，并在數據的預處理上進行了地表溫度的高程訂正和植被指數的平滑除噪。

結果表明：遙感估算的地表蒸散發量與地表的干濕狀況相一致，較好地體現了流域蒸散發量的季節性變化特征。不同地表覆蓋類型下，水體的蒸發量最高；受土壤含水量空間差異性的影響，植被內部具有較大的方差，其在生長季節更為顯著；耕地和林地蒸散發量的差異與其季節性生長狀況有關；城鎮夾雜著植被和水體等具有較高蒸散發量的地物類型，呈現出較高的蒸發量。估算結果合理地反映了不同地表類型之間蒸發量的差異性。

本文在特征空間構造上，干濕邊采用最大最小值擬合方法確定，當研究區普遍干旱和濕潤情況下，結果存在一定的不確定性，如何對干、濕邊進行合理的修正，需要進一步研究，另外方法的精度有待進一步的分析驗證。

［1］ 武夏寧，胡鐵松，王修貴，等.區域蒸散發估算測定方法綜述［J］.農業工程學報，2006，10（22）：257－262.

［2］ BELLA C，REBELLA C M，PARUELO J M.Evapotranspiration estimates using NOAA AVHRR imagery in Pampa region of Argentima［J］.International Journal of Remote Sensing，2000，21（4）：791－797.

［3］ SOBRINO J A，GóMEZ M，JIMéNEZ－MU?OZ J C，et al.A simple algorithm to estimate evapotranspiration from DAIS data：Application to the DAISEX campaigns［J］.Hydrology，2005（315）：117－125.

［4］ BASTIAANSSEN W G M，MENENTIA M.A remote sensing surface energy balance algorithm for land（SEBAL）［J］.Hydrology，1998（212）：198－212.

［5］ 倪猛，陳波，岳建華，等.洛河流域蒸散發遙感反演及其與各參數的相關性分析［J］.地理與地理信息科學，2007，23（6）：34－37.

［6］ MU Q Z，HEINSCH F A，ZHAO M S，et al.Development of a global evapotranspiration algorithm based on MODIS and global meteorology data［J］.Remote Sensing of Environment，2007（111）：519－536.

［7］ VENTURINI V，ISLAM S，RODRIGUE Z.Estimation of evaporative fraction and evapotranspiration from MODIS products using a complementary based model［J］.Remote Sensing of Environment，2008（112）：132－141.

［8］ JIANG L，ISLAM S.A methodology for estimation of surface evapotranspiration over large areas using remote sensing observations［J］.Geophysical Research Letters，1999，26（17）：2773－2776.

［9］ BATRA N，ISLAM S.Estimation and comparison of evapotranspiration from MODIS and AVHRR sensors for clear sky days over the Southern Great Plains［J］.Remote Sensing of Environment，2006（103）：1－15.

［10］ PRINCE J C.Using spatial context in satellite data to infer regional scale evapotranspiration［J］.IEEE Transactions on Geoscience and Remote Sensing，1990，28（5）：940－948.

［11］ CARLSON T.An overview of the"Triangle Method"for esti－mating surface evapotranspiration and soil moisture from satellite imagery［J］.Sensors，2007，7（8）：1612－1629.

［12］ NISHIDA K，NEMANI R R，RUNNING S W，et al.An operational remote sensing algorithm of land surface evaporation［J］.Geophysical Research，2003，108（D9）：4270.

［13］ BISHT G，VENTURINI V.Estimation of the net radiation using MODIS（Moderate Resolution Imaging Spectroradiometer）data for clear sky days［J］.Remote Sensing of Environment，2005，97（1）：52－67.

［14］ DUFFIE J A，BEXKMAN W A.Solar Engineering of Thermal Process［M］.New York：Wiley－Interscience，1980.

［15］ SU Z.The Surface Energy Balance System（SEBS）for estimation of turbulent heat fluxes［J］.Hydrology and Earth System Science，2002，6（1）：85－99.

［16］ 盧遠，華璀，韋燕飛.利用MODIS數據進行旱情動態監測研究［J］.地理與地理信息科學，2007，23（3）：55－58.

［17］ JARVIS P G.The interpretation of the variations in leaf water potential and stomatal conductance found in canopies in the field［J］.Philosophical Transactions of the Royal Society of London，1976，273：593－610.

［18］ MO X G，LIU S X，LIN Z H，et al.Simulating temporal and spatial variation of evapotranspiration over the Lushi basin［J］.Hydrology，2004，285：125－142.

［19］ SUN Z G，WANG Q X，OUYANG Z，et al.Evaluation of MOD16 algorithm using MODIS and ground observational data in winter wheat field in North China Plain［J］.Hydrological Processes，2007，21：1196－1206.

［20］ LDAS.Land Data Assimilation Systems［EB／OL］.http：／／ldas.gsfc.nasa.gov／nldas／web／web.veg.table.html.2012－01－31.

［21］ THOM A S.Momentum，Mass and Heat Exchange of Plant Communities［M］.London：Vegetation and the Atmosphere Academic Press，1972.

［22］ HIRD J N，MCDERMID G J.Noise reduction of NDVI time series：An empirical comparison of selected techniques［J］.Remote Sensing of Environment，2009，113：248－258.

［23］ VELLEMAN P.Definition and comparison of robust nonlinear data smoothing algorithms［J］.The American Statistical Association，1980，75：609－615.

［24］ HASSAN Q K，BOURQUE C P A，MENG F R，et al.A wetness index using terrain－corrected surface temperature and normalized difference vegetation index derived from standard MODIS products：An evaluation of its use in a humid［J］.Sensors，2007（7）：2028－2048.

［25］ 張可慧，劉劍鋒，劉芳圓，等.河北省潛在蒸發量計算與變化趨勢分析［J］.地理與地理信息科學，2010，26（6）：75－78.

［26］ 郭曉寅.黑河流域蒸散發分布的遙感研究［J］.自然科學進展，2005，15（10）：1266－1270.

［27］ 劉朝順，高志強，高煒.基于遙感的蒸散發及地表溫度對LUCC響應的研究［J］.農業工程學報，2007，23（8）：1－8.

［28］ 何延波，SU Z，JIA L，et al.遙感數據支持下不同地表覆蓋的區域蒸散［J］.應用生態學報，2007，18（2）：288－296.