基于特征空间法的美国南部大平原蒸散发遥感反演

涂晨雨，官云兰，陈远芳

(东华理工大学测绘工程学院，330013，南昌)

0 引言

蒸散发(Evapotranspiration，ET)是指植被及地面向大气输送的水汽总通量，其过程需要消耗约60%的降水和地表净辐射，是维持全球陆-气水分和能量平衡的关键要素[1]。作为地表水分循环的主要组成部分，准确估算蒸散发对于农业灌溉与干旱监测、水资源管理、气候变化预估等研究至关重要，在水文学、生态学、气象学和农学等诸多领域备受关注[2]。

传统蒸散发估算方法以“点”的观测资料为基础，如水量平衡法、涡动相关法和大型蒸渗仪法等[3]。这些方法虽然观测精度较高，然而对于区域尺度而言，适用性却较差，不仅无法充分反映地表空间异质性，同时过于依赖地面数据的支持，在观测站点稀疏及资料匮乏地区应用较为困难[4]。遥感技术具有快速、周期性、宏观等特点，其迅速发展为区域或全球尺度的蒸散发估算提供了重要的驱动数据。在过去的几十年间，一系列基于遥感技术的蒸散发估算模型应运而生，主要分为以下三大类：经验统计模型、地表能量平衡模型和特征空间模型[5]。其中，特征空间模型是指地表温度、植被指数、反照率等遥感参数，在区域尺度上，构成的散点图呈三角形或梯形，进而以插值算法获得蒸散发[5]。该模型通常以地表温度和植被指数组成的特征空间为主，与其他模型相较，能有效地降低蒸散发估算对地面数据的依赖，提供了一个比较直接的反演地表通量的方法，具有较好的应用前景[6]。

对上述特征空间模型研究的过程中，许多学者发现模型估算蒸散发的准确度受干湿边界(三角形或梯形空间的上边界和下边界)限制。为了有效确定模型的干边及湿边，研究人员提出了经验回归法和理论计算法[7]。然而，前者受遥感数据空间尺度效应影响较大，而后者则需近地面气象和植被辅助，两者均有着不确定性[7]。因此，Zhu[8]等在提出的通用三角法中，将干边的确定从区域转化为像元尺度，同时仅需风速作为唯一辅助数据，有效地降低观察到的干边所涉及的主观性和不确定性。但该方法在极端干旱条件下的空气温度估算及湿边反演时，仍然受区域限制，具有经验性。

针对该方法存在的不足，本文在Szilagyi[9]等的研究成果的基础上，进一步给出了干湿边界相应的理论公式。同时，本文选定美国南部大平原(SGP)为研究区，利用MODIS(Moderate Resolution Imaging Spectroradiometer)遥感数据，构建时空二维地表温度-植被指数特征空间，实现2017年像元尺度SGP区域地表蒸散发的遥感连续估算。此外，还采用大气辐射测量(ARM)计划提供的站点数据对研究结果进行验证，评估该方法的适用性，并进一步分析研究区蒸散发时空变化规律。

1 研究区概况及数据来源

1.1 研究区概况

研究区位于美国大平原南部，经度范围为W95°18′～99°30′，纬度范围为N34°30′～38°30′，覆盖了约2/3的俄克拉荷马州以及堪萨斯州部分地区。图1展示了该地区多样的土地覆盖类型，其中以种植冬小麦、棉花和玉米等的农用地以及草地为主。整个研究区地势平坦，气候类型丰富，地表通量属性多样，因此被大气辐射测量(ARM)计划选定为第一个气侯试验场，拥有大量的连续测量大气辐射的观测站，可给予充足可靠的气象数据，为蒸散发模型研究及验证提供基础[10]。

图1 研究区土地覆盖类型及站点分布图

1.2 数据来源

本文使用数据包括遥感数据和气象站点数据。遥感数据主要来自于搭载在EOS-AM(上午轨道)和PM(下午轨道)系列卫星上的传感器MODIS。目前为止，共有44个MODIS数据产品发布，本文中涉及的MODIS产品为MOD03、MOD06_L2、MOD07_L2、MOD11A1、MOD13A2和MCD43A3，时间范围为2017年。上述产品在使用前均需通过预处理，提取研究所需有效数据，并重采样至相同的空间分辨率1 km。其中，MOD13A2提供的归一化植被指数(NDVI)用于估算研究区的植被覆盖度(fc)，该产品时间分辨率为16 d，不同于其他产品(1 d)，需进一步插值处理得到日尺度数据。气象观测数据来自于大气辐射测量计划(https://www.arm.gov/)的波文比能量平衡系统(EBBR)。该系统可提供气温、湿度和风速等气象数据，其中风速为模型所需唯一辅助数据。同时，该系统获得的显热通量(H)、潜热通量(LE)和地表净辐射，将用于验证研究过程中蒸发比(EF)、净辐射以及蒸散发的估算值。其中蒸发比的计算公式如下：

(1)

2 研究方法

本文采用遥感蒸散发估算模型为地表温度—植被指数特征空间，即在研究区内地表温度(TS)和植被指数(VI)构成的散点图，将呈三角/梯形空间，其中该空间形状的上下边界为干湿边界。通过确定特征空间模型构建的关键参数(干边和湿边)，进而获得修正温度植被干旱指数(MTVDI)求解逐像元蒸发比，在此基础上结合地表净辐射数据，估算研究区蒸散发。

2.1 干湿边求解

在水分胁迫条件下，地表温度将受地表能量影响，干边(TSmax)代表的是不同植被指数所对应的地表最大温度。根据地表能量平衡原理，并参考Zhu等[8]的研究，TSmax求解见公式(2)。湿边(TW)表示在水分饱和情况下的地表温度，本文选择湿球温度(TWb)来代替，该参数估算根据Szilagyi等[9]的研究成果表达，见公式(3)。

(2)

(3)

式中：εSS和CS分别为裸土地表发射率和纯裸土土壤热通量占地表净辐射的比例，取值为0.95和0.315；σ为Stefan-Boltzmann常数(5.67×10-8W/m2·K4)；cp为空气定压比热容(1 012 J/kg·K)；γ为湿度计常数(0.066 4 kPa/℃)；Ta(K)为近地表空气温度，根据Zhu等[11]提出的方法由MOD07_L2和MOD06_L2产品求得；Td(K)是由MOD07_L2产品提供的露点温度；Sd(W/m2)和εa分别表示下行太阳辐射和空气发射率，ras表示裸土的空气动力学阻抗，受地表粗糙度和风速等多影响；而Tdry表示在极端干旱条件下的空气温度，4个参数求解公式见式(4)～(9)。

(4)

(5)

(6)

(7)

式中：S0为大气层顶的太阳常数(1 367 W/m2)；β依据Bisht and Bras的研究[12]设为0.2；Lv为蒸发潜热，约为2.5×106J/kg；Rv是水汽的气体常数(461 J/kg·K，)；e0(hPa)是近地表水汽压，而太阳天顶角参数θ(rad)为MOD03产品提供。

(8)

(9)

式中：z(m)是参考高度，d(m)和z0m(m)是零平面位移高度和动量的表面粗糙度长度，对于裸土而言两者分别为0和0.005；k是冯卡曼常数0.41；u(m/s)是在参考高度测量的风速，采用所有站点位置测量仪器高度2 m处风速的平均值代替；e*表示对应温度的饱和水汽压。

2.2 遥感蒸散发估算

在忽略日尺度土壤热通量(G)的前提下，日蒸散发(ET)的估算公式如式(10)所示。

ET=Rn,day·EFday

(10)

式中：Rn,day为日尺度地表净辐射，是在瞬时地表净辐射(Rn)遥感反演的基础上，根据Rivas &Cormora[13]研究结果估算；EFday为日尺度蒸发比，可用卫星过境时刻的瞬时蒸发比近似代替[14]。

蒸发比依据Jiang & Islam[15]等研究采用Priestley-Taylor方程[16]进行估算，见公式(11)。

(11)

式中：△表示饱和水汽压随气温(Ta)变化的斜率；φ是无量纲参数，反映实际空气动力学和地表阻抗，是基于土壤湿度进行参数化获得，不同像元对应的φi会随着土壤水分有效性从干边向湿边增加由0到1.26线性变化。对于纯植被而言，φ值(φcanopy)恒等于1.26。而对于纯裸地的φ值(φsoil)，则需要采用能反映地表土壤湿度的修订型温度植被干旱指数(MTVDI)来计算。该指数为Zhu等[17]在Sandholt等[18]发现地表植被指数与土壤湿度之间存在变化关系，并据此提出的温度植被干旱指数(TVDI)基础上进一步改进获得，公式如下：

(12)

(13)

(14)

式中：Tsoil为任意像元通过温度分解求得的纯裸土地表温度；TS地表温度由MOD11A1产品提供；fc植被覆盖度可利用NDVI来估算，其中NDVImin和NDVImax分别设为0.05和0.86。

基于此，φsoil可表达为MTVDI的函数关系，见公式(15)，而混合像元的φ值以植被覆盖度(fc)为权重，通过线性插值的方法求得：

φsoil=1.26[1-exp(MTVDI-1)]

(15)

φ=φsoil(1-fc)+φcanopyfc

(16)

2.3 检验指标

本文采用4个统计指标来对蒸散发估算结果进行精度评估，分别为相关系数(r)、绝对偏差(MAE)、均方根误差(RMSE)和偏差(B)，公式如下：

(17)

(18)

(19)

(20)

3 结果与分析

3.1 精度验证

地表温度—植被指数特征空间模型主要估算的参数为蒸发比，需要地表净辐射数据结合获得蒸散发。因此，蒸散发估算结果的误差来源于两部分，一是蒸发比的估算；二是地表净辐射的估算。本文结合站点实测数据，分别对蒸发比、地表净辐射和蒸散发估算结果进行精度验证。

在站点尺度上，2017年本文估算的蒸发比精度如表1所示。总体而言，蒸发比存在高估现象，平均偏差为0.05，其中在E12和E40站点处蒸发比高估显著，偏差达到0.23，而在E32站点处偏差最小(B=0)。此外，蒸发比仅在E9、E11和E13站点存在低估，三者偏差分别为-0.14、-0.27和-0.10。各站点蒸发比估算值在0.10～0.84之间，与观测值比较，两者平均MAE、RMSE和r分别为0.13、0.15和0.69，在一定程度上达到了精度要求。以MAE、RMSE和r三者为蒸发比估算精度评价指标，从单个站点来看，E35站点处模型模拟效果总体而言最好(MAE=0.11、RMSE=0.14，r=0.76)，而在E11站点处模型模拟效果最差(MAE=0.16、RMSE=0.19)，但该站点的相关性却较高(r=0.81)。

表1 各站点观测日蒸发比与本文估算值比较

本文使用的地表净辐射数据是参考Bisht & Bras[12]和Rivas &Cormora[13]研究结果并结合MODIS数据估算获得。地表净辐射估算精度如表2所示，在各站点处，净辐射的估算结果总体略低于观测值，平均B为-0.03。从各指标来看，净辐射估算结果较好，总体精度较高。具体而言，净辐射MAE范围为19.26～24.53 W/m2，平均值21.54 W/m2，而RMSE范围则在24.63～31.03 W/m2，平均值达到27.26 W/m2，其中在E39站点处2个值均最高，而E35站点均最低。此外，净辐射估算值与观测值在站点尺度具有很好的正相关性，r在0.87～0.95之间，两者的一致性较好。

表2 各站点观测日净辐射与本文估算值比较

从蒸散发估算精度评价结果来看(表3)，本文蒸散发平均MAE和RMSE分别为0.82 mm和1.07 mm，平均偏差则为0.06。因此，可以说明蒸散发估算值与观测值相比存在一定程度的高估，同蒸发比估算结果相一致。除此之外，蒸散发同样在E11站点处模拟效果最差(MAE=1.36 mm，RMSE=1.78 mm)，而E34站点处模拟效果最好略优于E35站点。同时，结合站点尺度日蒸散比、净辐射及蒸散发三者估算值与实测值之间的散点图来分析(图2)，蒸散发估算值与观测值相关系数平均值为0.86，高于蒸发比(r=0.69)，略低于净辐射(r=0.91)。综合来看，本文估算的蒸散发结果与观测值一致性较好，基本满足精度要求，并且该结果受蒸发比估算精度影响较大，受净辐射影响影响相对较小。

图2 站点尺度日蒸散比、净辐射及蒸散发三者估算值与实测值的散点图

表3 各站点观测日蒸散发与本文估算值比较

3.2 蒸散发时空分析

2017年，SGP研究区蒸散发为978.98 mm，月际变化如图3所示，蒸散发变化呈显著单峰型。其中，1—6月蒸散发持续增加，于7月达到顶峰160.74 mm，之后便逐渐下降，于12月达到低谷，蒸散发为16.09 mm。其中，4—9月温度适宜，降水充足，为植被生长茂盛期，因而蒸发与蒸腾过程剧烈，约占全年蒸散发的79%。同时，这还与农作物种植周期有关，该地区农作物类型以冬小麦、玉米和棉花为主，冬小麦收割时段为5月下旬至6月上旬，而玉米及棉花约在9—10月采摘，从而在此之后蒸散量显著降低。统计2017年研究区四季蒸散发总量，不同季节之间存在较大差异，春夏秋冬四季的蒸散发均值分别为251.06 mm、454.88 mm、201.54 mm和71.50 mm。其中，夏季(6—8月)蒸散发最高，占全年总量的46.46%；

图3 2017年SGP研究区蒸散发月际变化

而春季(3—5月)蒸散发则略高于秋季(9—11月)；冬季(12—2月)蒸散发最低，仅有年蒸散发的7%。因而，夏季实际蒸散发对于年地表蒸散发总量具有较为明显的贡献。

通过统一空间分布图的显示范围，进一步分析研究区四季蒸散发空间分布情况(图4)。根据图4可以看出，研究区春夏秋冬四季蒸散发西北-东南分布特征明显，并且呈东南向西北递减趋势，其中蒸散发高值集中在林区，草原次之，而农业区蒸散发则最弱。由于研究区地处美国南部大平原，植被覆盖率较高，蒸散发总量受植被蒸腾影响较大，从而两极化相对较小。然而，考虑到受季节性气温、日照时长、降水量、相对湿度以及平均风速等年内分布不均的影响，研究区四季蒸散量空间分布仍具有一定差异性。对比春季和秋季蒸散发分布情况，两季蒸散发范围为21～300 mm，其中春季蒸散发高值区约占研究区1/5，而秋季蒸散发高值区相对较少，仅存在于东南地区。由于春季太阳辐射增强和地表温度升高，而秋季与之相反，从而植物在春季蒸腾作用较强，进而导致春季蒸散发略高于秋季。对比夏季和冬季，两季蒸散发差异显著，夏季气温较高、降水充足，各地区蒸散发均高于200 mm，而冬季与之相反，气温较低，降水稀少，蒸散发均低于200 mm，并且空间分布变化也较小。

图4 2017年SGP研究区四季蒸散发空间分布图

4 结论

针对特征空间法干湿边界确定时存在的主观性和不确定性，本文结合Zhu等[8]和Szilagyi等[9]的研究成果，给出相应的干湿边界理论公式，并利用地表温度—植被指数特征空间法，实现2017年美国南部大平原地区日尺度蒸散发估算。同时，采用站点观测数据对估算结果进行精度验证，并探析该研究区蒸散发的时空分布情况。

1)在站点尺度，日蒸散发估算值与观测值相关系数达到0.86，具有较好的一致性，两者平均MAE和RMSE分别为0.82 mm和1.07 mm，满足一定的精度要求。此外，本文蒸散发误差来源于蒸发比和地表净辐射的估算，其中蒸散发与蒸发比同观测值相比，均存在高估，而地表净辐射则略低于观测值。

2)2017年研究区蒸散发为978.98 mm，从月际变化来看，蒸散发呈现为单峰型分布，7月为峰顶，12月则为低谷，蒸散发主要集中在4—9月份，与植被生长周期相一致。从季节变化来看，研究区蒸散发均呈东南高西北低的分布趋势，并且蒸散发季节差异显著，蒸散发大小按季节排序依次为夏季>春季>秋季>冬季。