基于植被指数－地表温度特征空间的伊河流域蒸散发量估算

周 峰，王文，王晓刚

（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］；刘朝顺等的研究同样表明受水分状况和植被覆盖度的影响，高蒸发量主要分布在水体、沼泽地区，而裸地及低植被覆盖度的自然植被区蒸发量较小［27］。而有关城镇蒸发特征的研究中，何延波等利用中分辨率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.