基于表面能量平衡模型的张承地区蒸散发研究

董祥旺，金晓媚，张绪财，殷秀兰，金爱芳，郎 捷，罗绪富，马靖宣

（1.中国地质大学（北京）水资源与环境学院，北京 100083；2.中国地质环境监测院，北京 100081）

蒸散发包括来自土壤、水域、冠层表面等各类地表的蒸发，也包括植物气孔的蒸腾[1]，准确的估算蒸散量对于水资源合理利用及生态环境保护具有极其重要的意义[2]。蒸散量估算方法大致可以分为实测法和模型法2类。实测法主要包括波文比法、涡度相关法、液流法、蒸发皿和蒸渗仪等方法，但是传统实测方法难以获得区域尺度的蒸散发[3]。随着遥感技术的兴起，出现许多与遥感相结合估算蒸散量的模型，主要分为与传统方法相结合的模型（如彭曼公式、互补相关理论模型、Priestley-Taylor模型）、经验统计模型、特征空间法、垂向能量平衡模型（如SEBAL、SEBS、TSEB）[4]。其中，Su[5−6]提出的表面能量平衡模型（Surface Energy Balance System, SEBS）已被证明能够在各种尺度上以可接受的精度估计湍流通量和蒸发比，众多学者[7−10]将该模型应用于不同空间和时间尺度上，均得到了较好的验证。

近年来，随着遥感云计算平台的发展，发布了众多地表蒸散量遥感产品，如MOD16A2（8 d, 500 m）、Gl-ASS-ET（8 d, 1 km）、GLEAM（1 d, 0.25°）等蒸散发遥感产品[11]，其中GLASS-ET，GLEAM等数据空间分辨率较低，不适用于空间面积较小的地区；而时空分辨率较高且应用广泛的MOD16A2数据产品在无植被覆盖区域进行了无值化处理，对于植被覆盖不是很高的地区也不能完全反映区域实际蒸散的空间分布。

张家口承德地区（简称“张承地区”）作为京津冀城市群生态安全的重要屏障，是生态环境关注的热点地区，2015年《全国生态功能区划（修编版）》将河北省的张家口、承德以及北京市北部城区等地区划为京津冀北部水源涵养重要区，并指出该地区存在水资源过度开发等生态环境问题。在京津冀协同发展的战略中，以往的蒸散研究多将京津冀作为一个整体进行探究，如于占江[12]等基于高桥公式对京津冀地区的实际蒸散开展了研究，发现京津冀地区的蒸散量呈现不显著上升趋势，但并未对其空间变化特征进行分析；少数学者对张家口部分地区的蒸散量曾进行过研究，如石嘉丽[13]等基于GLASS蒸散产品分析了河北坝上地区蒸散量的时空变化趋势，研究发现坝上地区的蒸散量呈略微下降趋势并且受土地利用变化影响显著。虽然前人的研究取得了一些重要成果，但是就张承地区而言，仍然缺少对该地区长时间序列的实际蒸散量空间变化规律以及影响因素的分析，不能完全反映该地区实际蒸散的时空分布特征。

因此，本文以张承地区为研究区，基于连续的MODIS产品数据和GLDAS气象数据，应用SEBS模型计算2001年1月——2020年12月逐月的区域蒸散量，与MOD16A2产品数据时间序列的变化趋势进行对比，通过野外实测数据在像元尺度上对蒸散量计算结果进行验证。在此基础上，对研究区20 a间的蒸散量时空变化特征及影响因素进行分析，所得结果期望为张承地区水资源的可持续利用以及生态环境保护提供参考。

1 研究区概况

张家口承德地区位于河北省西北部，北邻蒙古高原，南接华北平原，总面积约为7.6×104km2（图1）。研究区内以阴山山脉为界，分为坝上和坝下地区，在地势上呈现出西北高东南低的地势特点，研究区属于温带季风性气候，夏季高温多雨，冬季寒冷干燥，据张家口和承德2个气象站统计，多年平均降水量约为460 mm，年均气温9～12 °C，年日照率约为60%，日照充足。

图1 研究区及野外试验点位置示意图Fig.1 Location of the study area and the test site

2 研究方法与数据

2.1 SEBS模型

2.2 Sen+Mannkendall趋势分析

Sen趋势度可以很好地减少噪声的干扰，并判断趋势的上升或是下降[14]， Mannkendall检验是评估趋势显著性的非参数性检验，在水文趋势检验研究中得到了广泛的应用[15−16]，其优点是不需要测量值服从正态分布，也不要求趋势是线性的，并且不受缺失值和异常值的影响。采用上述2种方法结合可以增强方法的抗噪性，并在一定程度上提高检验结果的准确性[17]。

Sen趋势度计算公式为：

式中：xi、xj————时间序列数据。

β大于0表示时间序列呈现上升趋势，反之处于下降趋势。

MK检测过程如下：

假设H0：数据样本{xi,i=1,2,3,···,n}是独立且同分布的，H1：序列存在单调趋势。M-K的统计量S定义为：

式中：n——数据集的长度。

sgnθ定义为：

根据数据集的长度n值选取显著性检验统计量：当n<10时，使用统计量S进行双边趋势检验，在给定显著性水平α下，如果 |S|≥Sα/2则拒绝H0，认为原序列存在显著趋势；否则接受H0，认为序列趋势不显著[15]。当n≥10时，使用检验统计量Z进行趋势检验，Z值计算如下：

式中：m——序列中结（重复出现的数据组）的个数；

ti————结的宽度（第i个重复数据组中的重复数据个数）。

统计量Z值同样采用双边趋势检验，当|Z|>Z1−α/2时，则拒绝H0，认为原序列存在显著趋势；否则接受H0，认为序列趋势不显著。

2.3 相关性分析

（1）Pearson 相关性分析

Pearson相关性分析是分析2个变量间线性相关的方法，并且要求2个变量都呈正态分布，而且是随机变量。

（2）偏相关性分析

偏相关系数是衡量多个变量中某2个变量之间的线性相关程度的指标[18]，偏相关系数绝对值愈大，表明变量之间的线性相关程度愈高；反之愈低。样本的偏相关系数rij·h计算公式为：

式中：rij————变量xi和xj的简单相关系数；

rih————xi和xh的 简 单 相 关 系 数 ；

上述数据传输电路单元得到的电压数据和电流数据经过分压电阻R2和分压电阻R3流向通讯收发芯片U6，通讯收发芯片U6自带CAN总线通讯协议，在接收到单片机U1传输的电压数据和电流数据后对其进行通讯协议转化，转化后的电压数据和电流数据信号流向共模滤波电感L3，滤除掉信号中的干扰成分，并经过电阻R12和电阻R13的分压保护，经过瞬态抑制二极管Z1和瞬态抑制二极管Z2后流向保险F1和保险F2，最终通过接线端子J2和外部CAN总线相连，并通过CAN总线将测量得到的电压数据和电流数据上传至实时监测光伏组件运行状态的数据处理计算机，完成整个检测流程.

rjh————xj和xh的简单相关系数。

计算样本的检验统计量t，确定P值，做出推断结论，统计量t计算公式为：

式中：r——偏相关系数；

k——样本数；

q——偏相关阶数。

统计量服从k−q−2个自由度的分布。

2.4 研究数据

（1）MODIS 数据

（2）气象数据

由于研究区只有张家口承德2个气象站，所以本文采用2001——2020年间GLDAS数据中的气温、气压、风速、相对湿度、总降雨量5个波段，该数据的时间分辨率为1月、空间分辨率为25 km，GLDAS数据来源于GES DISC(https://disc.gsfc.nasa.gov/)。

（3）高程数据

本文选取高程SRTMDEM 90 m分辨率原始高程数据、SRTMSLOPE 90 m分辨率坡度数据、SETMASPECT 90 m分辨率坡向数据，数据来源为地理空间数据云（http://www.gscloud.cn/）。

（4）土地利用数据

为研究用地类型对蒸散发的影响，本文选取2000、2005、2010、2015、2020年5期中科院的中国土地利用遥感监测数据，空间分辨率为30 m，数据来源为资源环境与科学数据中心（https://www.resdc.cn/）。

为保证数据的一致性，所有数据均重采样为0.005°（约为500 m），地理坐标系为WGS84。

3 研究成果

3.1 蒸散量的时空分布特征

从张承地区2020年蒸散量（ET）的空间分布图（图2）可以看出，坝下地区的蒸散量明显高于坝上地区，其中承德地区的中部和西南部地区以及张家口的东南部地区蒸散量较高，年蒸散量主要在800～1 300 mm变化；而在阴山以北的坝上地区，包括承德围场北部靠近内蒙古的区域、张家口北部的康保和尚义以及沽源等地区的蒸散发相对较小，年蒸散量在60～700 mm之间变化。

图2 2020年研究区蒸散量空间分布图Fig.2 Spatial distribution map of evapotranspiration in the study area

通过对蒸散量统计结果进行分析可知，张承地区年蒸散量呈现略微波动上升趋势见图3，最大值为2013年的545 mm，最小值为2002年的348 mm；研究区内张家口承德2个地区在20 a间的蒸散变化趋势较为一致，且承德地区的蒸散量明显高于张家口地区，承德地区在2013年达到其最大值583 mm，张家口地区在2003年达到最大值508 mm。研究区多年平均的年内月蒸散量分布见图4，6——8月的蒸散量较高，12月至次年2月的蒸散量最小。

图3 张承地区2001——2020年蒸散量值Fig.3 Evapotranspiration values of Zhangjiakou and Chengde from 2001 to 2020

图4 研究区月均蒸散量值Fig.4 Average monthly evapotranspiration values in the study area

3.2 蒸散量结果验证

（1）与 MOD16A2 数据对比

由于MOD16A2数据在无植被覆盖区域进行了无值化处理，因此本文仅将计算结果与其做趋势上的对比。从图5可以看出，SEBS计算的月蒸散量值与MOD16A2数据在趋势上具有较强的一致性，确定系数为0.79，吻合较好，说明SEBS反演结果具有较高的可靠度。

图5 研究区2001——2020年SEBS与MOD16A2数据对比Fig.5 Comparison of the SEBS and MOD16A2 data from 2001 to 2020 in the study area

（2）与野外实测数据对比验证

本文利用自制蒸渗仪在研究区内选取了石庄屯（果园）、姚家庄（高植被覆盖）、王家楼（中等植被覆盖）3个地点进行了现场蒸发试验，对模型反演结果进行验证，试验期为2021年7月9——23日。

为保证验证结果的准确性，本文选取了2021年7月的MODIS数据以及GLDAS数据基于SEBS模型计算出7月的蒸散量与试验结果进行对比验证（表1），从结果可以看出，王家楼、姚家庄试验点的观测值与反演值的误差约为0.3 mm，石庄屯试验点的误差较小，仅为0.12 mm，相对误差为5.77%，其他2个试验点的误差均小于15%，总体来说SEBS模型反演的结果与现场实测值吻合相对较好。

表1 3个野外蒸发试验日蒸散发结果统计Table 1 Statistics of daily evapotranspiration results of three field evaporation tests

3.3 Sen+MannKendall的蒸散量变化趋势分析

为进一步探究研究区蒸散量的变化趋势，本文基于Sen+Mannkendall趋势检测方法对20 a研究区的蒸散量变化趋势进行分析。利用Sen斜率计算的β值作为变化斜率。根据经验判断，取0≤β≤10，蒸散量呈稳定趋势，由于研究序列长度为20，采用统计量Z进行检验，取显著水平α=0.05，Z1−α/2=Z0.975=1.96，按照表2将研究区蒸散量变化趋势分为5类，得到研究区蒸散变化趋势空间分布（图6）。

图6 研究区20 a实际蒸散量时空变化趋势图Fig.6 Trend of the actual evapotranspiration over the period of 20 years

表2 研究区蒸散量变化趋势分类标准Table 2 Classification criteria for the trend of evapotranspiration in the study area

对研究区蒸散量变化趋势进行统计（表3），从结果可以看出：蒸散量基本稳定的区域占总面积的75.41%；蒸散量显著增加的区域占据了研究区的5.12%，主要分布在张家口东部和南部地区以及承德的中部地区；蒸散量轻微降低的区域占研究区的18.35%，主要分布在张家口的张北北部、九连城、沽源北部、蔚县和怀来地区以及承德的围场北部地区、承德的南部地区；显著降低的区域占总面积的1.11%，主要分布在承德围场北部靠近内蒙古的地区，以及张家口安固里淖附近地区、怀来的部分地区。

表3 研究区蒸散发空间变化趋势面积统计Table 3 Area statistics of evapotranspiration spatial change trend in the study area

4 影响因素分析

4.1 气象因素

本文利用GLDAS 2001——2020年20 a间每月的气温降水数据与SEBS反演的月蒸散量，基于偏相关分析分别分析气温、降水与蒸散量的相关程度，并进行T检验，得到气温和降水与蒸散量的偏相关性检验空间分布图（图7），结果表明：研究区内气温与蒸散量具有较强的正相关性，特别是坝下地区，说明气温是蒸散量极其重要的影响因素，气温越高蒸散量越大；降水对于区域蒸散量也有一定的影响，其中承德东南部地区降水与蒸散量具有显著正相关性，张家口大部分地区与承德西部地区具有不显著正相关性。

图7 气温与蒸散量和降水与蒸散量的偏相关性检验空间分布图Fig.7 Spatial distribution plots for the partial correlation between air temperature and evapotranspiration (a) and for the partial correlation between precipitation and evapotranspiration (b)

4.2 植被因素

植被蒸腾是蒸发的一种方式，通常情况下，植被覆盖越大，植被指数越高，区域蒸散量越大，反之，植被覆盖越低或无植被区域，区域蒸散量越小。为探究研究区内植被对蒸散发的影响，本文利用MOD09A1的近红外波段和红外波段计算每月的归一化植被指数NDVI，然后逐像元的计算NDVI与每月的SEBS反演的蒸散量的相关系数，得到植被指数与蒸散量的相关系数空间分布（图8）。从图中可以看出，以阴山山脉为界，阴山以南大部分区域，包括张家口地区的尚义县南部、下花园区、怀安县、阳原县以及蔚县东南部，承德地区的丰宁、承德县、承德市、滦平、隆化、兴隆的大部分地区的蒸散发都与植被有较强的相关性；区域蒸散发与植被覆盖不明显的区域主要分布在阴山以北地区，包括康保、沽源、张北以及承德的围场北部，这些地区大部分属于坝上地区，气候较为寒冷，除了夏季植被发育较好外，其他季节的植被覆盖均较差，因此蒸散发的变化与植被的相关性不明显。

图8 研究区植被指数与蒸散量的 Pearson相关系数空间分布图Fig.8 Spatial distribution of Pearson correlation coefficients of vegetation index and evapotranspiration in the study area

4.3 人为因素

人类活动对区域蒸散发的分布也有较大的影响，这主要体现在用地类型的变化上。为探究用地类型对蒸散发的影响，本文根据 2000、2005、2010、2015、2020年5期的中科院土地利用空间分布图，并依据中科院土地利用一级分类标准分为耕地、林地、草地、水体、建设用地和未利用土地6类，统计得到研究区不同土地用地类型的面积（表4），结果显示，2020年不同用地类型面积从大到小依次为：林地>耕地>草地>建设用地>水体>未利用土地；研究区内耕地面积明显减少，主要转化为林地和草地；研究区建筑用地的面积增长了近一倍，使得城市及周边地区的蒸散量出现了明显的下降，特别是张家口、宣化以及怀来城市带地区；研究区的未利用土地面积很少，零星分布于坝上地区，主要分布于康保和沽源地区，在各类用地类型中占比最小。利用2020年中科院土地利用数据对研究区20 a蒸散量的平均值进行统计，结果显示，不同用地类型的蒸散量从大到小依次为：林地>水体>草地>耕地>建设用地>未利用土地。受城市热岛效应的影响[19]，一般情况下建设用地的表面温度高于未利用土地，因此建设用地的蒸散量略高于未利用土地。

表4 张承地区2000——2020年不同土地利用类型的面积以及年均蒸散量统计Table 4 Area statistics of different land use types from 2000 to 2020 and average annual evapotranspiration statistics of different land use types in the ZhangCheng area

5 结论

（1）SEBS反演结果与MOD16A2数据、野外现场实测数据吻合度较好，验证了SEBS模型在研究区反演的结果具有较高的可靠性。

（2）研究区内张家口承德地区2001——2020年的蒸散量变化趋势较为一致，均呈现略微上涨的趋势，其中张家口市的多年平均蒸散量为454 mm，承德市的多年平均蒸散量为506 mm。

（3）研究区蒸散量与气温具有较强的正相关性，温度越高，蒸散量越大，降雨对蒸散量的变化也有一定的影响，但是影响没有气温显著；植被与蒸散量的变化呈现较强的正相关性，植被指数越高，蒸散量越大；不同用地类型的蒸散量从大到小依次为：林地>水体>草地>耕地>建设用地>未利用土地。

（4）人类活动会改变下垫面以及用地类型情况，间接影响蒸散量的空间分布，而蒸散作为水均衡中的重要组成部分，其长序列的变化研究对于了解张承地区区域水循环，提高水资源的高效利用和管理以及生态环境的保护具有重要的意义。