廖前瑜 ，任 超 ，冷 佩 ，段四波，韩晓静

（1. 桂林理工大学测绘地理信息学院，广西桂林541004；2. 中国农业科学院农业资源与农业区划研究所，北京100081）

0 引言

蒸散发是地表水分循环与能量平衡的重要组成分量，反映了大气与地球表面的水分和能量交换［1］。由于全球气候变化和人类对水资源的过度利用，导致了可利用水资源短缺和相关生态问题。因此，精确估算蒸散发在时空尺度上的变化，对农田灌溉管理、作物需水量估算、水资源优化配置等领域的研究具有重要意义。传统的蒸散发估算大多基于气象观测的单点计算，尽管能得到相对准确的蒸散发量，但无法实现非均匀下垫面的蒸散发研究，而遥感对于区域蒸散发估算具有无可替代的优势［2-3］。自20 世纪80 年代以来，众多学者围绕地表蒸散发遥感反演开展了大量工作，形成了许多切实可行的算法。这些算法大致可以分为基于能量平衡的方法［4-7］，基于水量平衡的方法［8］以及基于遥感地表温度—植被指数特征空间的方法［9-13］。其中，基于遥感地表温度—植被指数特征空间的方法近年来得到了最为广泛的应用，主要得益于其简单的原理和输入以及它能够避免用于计算水和热传递的空气动力学阻抗和表面阻抗的复杂参数化过程。应用地表温度—植被指数特征空间方法反演地表蒸散发的一个重要前提是要确定特征空间的干湿边。目前，确定地表温度—植被指数特征空间干湿边主要有基于传统的空间统计的方法［14-15］和基于能量平衡理论干湿边计算方法［13，16］。其中，基于传统的空间统计的方法通常具有一定主观性，且要求研究区气象条件一致，导致其在区域尺度上难以应用。基于能量平衡的方法则需要输入大量的参数，包括一些遥感难以获取的参数如风速、空气动力学阻抗及冠层阻抗等，中间参量过多会导致误差的累积和传递。尽管基于能量平衡估算干湿边的方法需要大量的输入参数，但因其具有良好的普适性和物理基础受到了研究学者们广泛的关注。

最近Leng 等［16］提出了逐像元确定地表温度—植被指数特征空间理论干湿边的思想，避免了传统特征空间方法在区域尺度应用的局限。由于在有云情况下无法获取有效的植被指数和地表温度参数，Leng 等［16］结合经典的Penman-Monteith 模型，发展了全天候地表蒸散发反演方法。然而，该研究仅在站点尺度进行了分析。在区域尺度上，个别站点观测的气象数据显然无法很好地描述气象因子的空间异质性。因此，文章在Leng等［16］研究的基础上，基于目前常用的卫星陆表遥感产品和较高空间分辨率的格网气象数据，深入探索区域尺度全天候蒸散发的反演。

1 数据

1.1 研究区与地面观测数据

黑河流域（97°E～102°E，38°N～43°N）面积约为 14.29 万 km2，位于我国西部甘肃省境内，流域气候主要受到中高纬度西风带环流和极地冷气团的影响。该流域从上游到中、下游，以水为纽带，形成了“冰雪／冻土—森林—草甸—人工/天然绿洲—荒漠—湖泊”的多元自然景观。流域内寒区和干旱区并存，地表异质性显著，尤其是山区冰冻圈和极端干旱的河流尾闾地区形成了鲜明对比。因此，准确获取流域内复杂环境状况下的地表蒸散发数据，对于深刻理解区域水循环和地气之间能量交换，具有重要意义。此外，该文选择黑河流域作为研究区的另一个重要原因是由中国科学院西北生态资源环境研究院主导开展一系列长期的生态水文实验，积累了丰富的水文气象观测数据，为研究工作的开展奠定了良好的基础。

该研究反演区域尺度全天候地表蒸散发，涉及的地面观测数据包括自动气象站气象观测数据和基于涡动的蒸散发观测数据。由中国科学院西北生态资源环境研究院主导的“黑河流域生态—水文过程综合遥感观测联合试验”于2012 年5 月底至9 月中旬在黑河流域中游地区布置了多套涡动相关仪和自动气象站，同步获取蒸散发与气象数据，为该文全天候蒸散的反演和验证提供了优质数据源。这些观测数据可从寒区旱区科学数据中心（http：//westdc.westgis.ac.cn/）免费下载。考虑到不同站点观测时期的差异，该文选取研究区内2012 年6 月1 日至9 月15 日的数据开展研究。根据获取的验证数据选择了4个具有代表性站点来探索区域地表全天候蒸散发反演，表1 给出了这4 个站点的编码名称、地表类型、高程、测量高度及经纬度。

图1 研究区域的地表高程Fig.1 DEM in the study area

表1 黑河生态水文遥感试验通量观测站点信息Table 1 Information of Heihe ecological hydrological remote sensing test flux observation sites

1.2 遥感数据

与地面观测数据对应，该文使用2012 年6 月1 日至9 月15 日的MODIS（Moderateresolution Imaging Spectroradiometer）陆表参数产品作为研究的遥感数据源。MODIS 是搭载在Terra 和Aqua 卫星上的一个重要传感器，具有36 个波段的观测通道，覆盖中国全境，具有较高的空间分辨率和时间分辨率。目前，基于MODIS 观测数据，已经发布了多种陆表和大气产品数据，在全球范围内得到了广泛应用。该文使用的MODIS 产品及其参数如表2 所示。

表2 MODIS 数据详情表Table 2 MODIS data detail

MODIS 产品数据来源于LAADS DAAC（Level-1 and Atmosphere Archive and Distribution System Distributed Active Archive Center）（https：//ladsweb.modaps.eosdis.nasa.gov/search/）网站。该文使用MRT 和MRTSwath 工具对不同级别的MODIS 数据进行重投影等预处理。研究中MOD021KM、MOD03 和MOD05_L2 产品被用来反演地表短波净辐射，其中MOD021KM 产品空间分辨率为250 m 和500 m 的大气顶部辐射亮度和反射率数据被整合成空间分辨率为1 km 的数据；MOD35_L2 云掩膜产品被用来判断MODIS 图像的每个像元是否受到云的干扰；MOD09A1、MOD11A1、MOD13A2 和MOD15A2 作为全天候地表蒸散发反演的参数输入，其中MOD09A1 和MOD15A2 在8 d 时间内被认为是不变的，MOD13A2 归一化植被指数通过hants 滤波后插值得到日尺度的植被指数［17］。

1.3 格网气象数据

研究所用的格网气象数据为中国气象局陆面数据同化系统CLDAS（CMA Land Data Assimilation System）的Verson 1.0 产品（简称CLDAS-v1.0），数据下载于中国气象数据网（http：//data.cma.cn/）。CLDAS-v1.0 产品覆盖范围为东亚区域（0°～60°N，70°E～150°E），空间分辨率为0.0625°，时间分辨率为1 h。研究中用到的CLDAS-v1.0 气象产品要素包括短波辐射、气温、风速和比湿。这些数据被重采样到MODIS 数据1 km空间分辨率。以2012 年6 月1 日为例，图2 显示了北京时间上午11：00 黑河流域的短波辐射、气温、风速和比湿。

2 方法

该文在Leng 等［16］建立的全天候地表蒸散发反演方法基础上，深入探索区域全天候地表蒸散发的反演。在该反演方法中，所有的像元根据是否存在有效的地表温度数据被划分为晴空和有云像元。其中，晴空像元的地表蒸散发利用逐像元地表温度—植被指数特征空间方法进行反演，有云像元的地表蒸散发则基于Penman-Monteith 模型计算。基于MODIS 陆表产品和CLDAS 格网气象数据，探索该全天候蒸散发反演方法在区域尺度上的可行性。图3 为区域尺度全天候蒸散发反演的技术流程图。

图2 2012 年6 月1 日CLDAS 产品：（a）短波辐射；（b）气温；（c）风速；（d）比湿Fig.2 CLDAS products on June 1，2012

图3 全天候蒸散发反演方法流程图Fig.3 Flow chart of the all-sky evapotranspiration retrieval method

2.1 LST-NDVI 特征空间估算晴天像元蒸散发

区别于传统基于空间信息的“地表温度—植被指数”特征空间，Leng 等［16］发展了独立于气象和下垫面条件的逐像元两阶段“地表温度—植被指数”特征空间新模式。图4 为示意图。该模式假设对于每一个晴空像元（如图1 中的像元C），都对应存在一个虚拟的特征空间，该虚拟特征空间由 4 个极端条件（Ts，max，Ts，min，Tv，max，Tv，min）对应的地表温度构成。其中，Ts，max和Tv，min两点构成的斜线把梯形特征空间划分为两个三角形，分别为由 Ts，max，Ts，min，Tv，min3 点构成的下三角形，由 Ts，max，Tv，max，Tv，min3 点构成的上三角形。在实际应用中，只要已知这4 个极端地表情形（即干燥裸土、受水分胁迫的全植被覆盖、供水良好的全植被覆盖和湿润裸土）对应的地表温度，便可确定像元C 的虚拟特征空间，从而能够根据像元C 在该虚拟特征空间中所处位置，估算蒸散发。该特征空间新模式的一个显著优点是在实际应用中能够避免传统基于空间信息确定干边和湿边时对气象和下垫面条件的苛刻要求，从而能够直接反演大范围蒸散发。

图4 LST-NDVI 特征空间示意图Fig.4 Schematic diagram of LST-NDVI feature space

LST-NDVI 特征空间假设NDVI 为某一值时，蒸散发随LST 增大而线性递增。当某像元的LST 和NDVI 落在下三角形中，土壤根层水分充足，植被蒸腾接近于潜在最大值，而表层土壤水分蒸散由像元处于三角形中的位置决定。在该情况下，土壤蒸发（ETs）和植被蒸腾（ETv）可以表示为：

式（1）～（2）中，Es，p和Ev，p分别为最大土壤蒸发量和最大植被蒸腾量。AB和AC为长度（图4）。

类似地，当某像元的LST 和NDVI 落在上三角形中，则表层土壤已干，土壤蒸发为0，植被开始受到胁迫，植被蒸腾量取决于该像元处于三角形中的位置。在该情况下，土壤蒸发（ETs）和植被蒸腾（ETv）可以表示为：

其中 Ts，max，Ts，min，Tv，max和 Tv，min分别是干燥裸土，饱和裸土，水分胁迫全植被和充分灌溉的全植被的理论地表温度（K），如图4 所示。Moran 等［18］基于能量平衡原理，假设4 种极端条件给出了4 个理论地表温度的计算公式：

式（7）中，Rn是地表净辐射（W/m2），rs是表面阻抗。ra是空气动力学阻抗，γ为干湿表常数（kPa/℃），Δ 为饱和水汽压随空气温度变化的斜率（kPa/℃），VPD为参考高度的空气饱和水汽压差（kPa），G是土壤热通量（W/m2），Cp为空气定压比热（J/（kg℃））。其中，rs和Rn参数的计算可参考 Leng 等［16］和 Todorovic 等［20］的文献。

式（5）中，Rn，s和Rn，v分别为土壤表面净辐射和植被表面净辐射（W/m2），G是土壤热通量（W/m2），Ta是空气温度（K）。Cv为空气定容比热（J/（℃m3）），γ为干湿表常数（kPa/℃），Δ 为饱和水汽压随空气温度变化的斜率（kPa/℃），VPD为参考高度的空气饱和水汽压差（kPa）。rcx是叶面气孔几乎完全关闭情况的冠层阻抗，设定为2 000 s/m；rcp是充足水分条件下的植被冠层阻抗，为100 与叶面积指数的比值；ra是空气动力学阻抗。

地表蒸散发可以表示为植被蒸腾和土壤蒸发的和：

2.2 Penman-Monteith 公式估算有云条件像元蒸散发

在有云条件下，无法有效获取LST 和NDVI 参数，因此使用经典的Penman-Monteith模型［19］估算地表蒸散发：

3 结果与分析

3.1 CLDAS 气象数据分析

该文使用所选4 个站点自动气象站实测的短波辐射、风速、气温、比湿对CLDAS 格网数据产品重采样为1 km 空间分辨率后的格网数据进行分析，结果如图5 所示。

图5 CLDAS 数据产品与实测数据散点图：（a）短波辐射；（b）风速；（c）气温；（d）比湿Fig.5 Comparison of CLDAS products and measured data

从图5 可以看出CLDAS 短波辐射、风速、气温、比湿数据误差较大，尤其是短波辐射和风速产品。这可能是由于CLDAS 气象数据在空间上受到气象站点疏密的影响，导致了在缺少气象站点的西北地区数据精度下降。Vinukollu 等［21］使用相同的一套驱动数据对比验证了P-M 模型［22］、P-T 模型［23］、表面能量平衡模型［5］，发现不同模型之间蒸散发反演的差异主要取决于净辐射数据的选择。考虑到短波辐射是地表蒸散发的能量来源，其精度高低不仅限制了地表蒸散发的反演结果，还是地表蒸散发的主要影响因子。基于此，该文采用Tang 等［24］的方法利用MODIS 数据估算短波净辐射转换为短波辐射来代替CLDAS 中的短波辐射，作为全天候蒸散发反演的输入参数。

为了验证MODIS 估算的短波净辐射在黑河地区的可靠性，该文选取了黑河流域4个站点的2012 年6 月1 日至9 月15 日的地面有效测量数据进行对比。其中，地面观测的短波净辐射采样频率为10 min，MODIS 数据反演的短波净辐射是分辨率为1 km 的瞬时值。另外，根据MOD35_L2 云产品判断每个像元是否会受到云的影响。MOD35_L2 中认为是clear 和probably clear 的像元，判断为晴空的像元；MOD35_L2 中认为是uncertain和cloudy 的像元，判断为有云的像元。图6 是该4 个站点的地表短波净辐射观测值和MODIS 估算的短波净辐射散点图，从图6 可以看出，晴空条件下，基于MODIS 反演的短波净辐射的均方根误差小于30 W/m2；有云时，均方根误差在70 W/m2左右。相比于利用CLDAS 数据获取的短波辐射精度，基于MODIS 反演的短波净辐射显然具有更高的精度。

图6 估算的地表短波净辐射与实测的地表短波净辐射散点图（a）1 号站；（b）7 号站；（c）10 号站；（d）17 号站Fig.6 Comparison of estimated net surface short-wave radiation and measured data

3.2 全天候区域地表蒸散发反演与分析

为了验证Leng 等［16］方法的可靠性，该文首先使用单点的实测气象数据和MODIS数据反演地表蒸散发，并利用黑河流域4 个站点的2012 年6 月1 日至9 月15 日的地面有效测量蒸散发数据进行对比。验证结果如图7 所示。使用单点的气象数据和MODIS数据反演的地表蒸散发RMSE 值在70 W/m2以内，具有较高的精度。这表明了利用Leng 等［16］的方法反演全天候地表蒸散发是可行的。在此基础之上，利用CLDAS 短波辐射、风速、比湿、气温的产品与MODIS 估算了黑河流域的全天候地表蒸散发，并将反演结果与地面实测数据进行了比较，4 个站点的平均偏差值为30.8 W/m2，平均均方根误差值为111.9 W/m2，误差相对较大。采用Tang 等［24］的方法转换得到短波辐射作为全天候蒸散发反演的输入参数后，精度有了明显的提高，其RMSE 在80 W/m2以内。

图7 不同驱动数据估算的蒸散发和实测的蒸散发散点图（a）1 号站、（b）7 号站、（c）10 号站、（d）17 号站Fig.7 Comparison of estimated evapotranspiration and measured data

基于能量平衡逐像元计算干湿边的方法具有良好的普适性和物理基础。然而大量的输入参数造成了误差的积累和传递。通过提高关键输入参数的精度，能够显著地提高蒸散发的反演精度。

4 结论

该文以黑河流域为研究区，基于遥感数据和较高空间分辨率的格网气象数据，实现了区域尺度地表蒸散发反演。在验证分析了CLDAS 数据精度的基础之上，利用MODIS数据估算的短波辐射代替CLDAS 短波辐射数据作为蒸散发反演模型的输入参数，最后利用2012 年6 月1 日至9 月15 日的CLDAS 数据和MODIS 数据反演了黑河流域的全天候地表蒸散发，并对结果进行精度检验。研究表明，基于MODIS 遥感数据估算得到精度较高的短波辐射替代CLDAS 中的短波辐射产品，能显著提高黑河地区蒸散发的反演精度。

虽然通过该文方法得到了精度较高的区域尺度全天候地表蒸散发，但研究仍然存在以下不足。首先是验证站点较少，且相对空间距离较近，不利于对地表蒸散发的全面评价；其次是由于蒸散发反演和验证过程中涉及不同尺度的数据源，该研究没有充分考虑尺度效应对反演结果可能造成的影响。