考虑坡面地形的蒸散发遥感估算及空间分布研究

尹 剑，欧照凡

(东北农业大学水利与土木工程学院，黑龙江 哈尔滨 150050)

蒸散发(Evapotranspiration，ET)是流域水循环过程和能量过程的重要组成部分[1]，准确估算实际ET，对合理调配水资源、应对地区旱情、指导灌溉具有重要的理论和应用价值[2]。遥感技术是估算区域ET的有效方法[3]。多年来，ET遥感模型由传统方法的空间扩展、经验统计公式，到基于能量平衡的单层、双层模型，取得了长足的发展[4]。其中，双层模型物理机理更完善，较好地描述了土壤-植被-大气间的水热交换[5]。张仁华等[6]结合实验发展了一种易操作的双层ET遥感模型，减小了温度分解困难和阻抗计算量大的问题，近年来算法不断改进，并取得了丰富的研究经验[6-9]。目前，ET遥感模型主要针对平坦地区，较少考虑地形因素。坡度、坡向等微地形直接影响太阳辐射和地表温度的计算，进而影响ET精度。傅抱璞[10]利用地形起伏度权重，提出了估算山区ET的方法。Thomton等[11]将高程和坡度融入山区ET的估算方法之中。Allen等[12]考虑坡度和坡向的影响，提出了METRIC模型。高永年等[13]考虑了地表净辐射的地形效应，构建了复杂地形ET模型。当前地形对于ET的影响多为非遥感手段。而遥感模型中的地表参数，特别是关键的太阳辐射，受地形遮蔽影响较大，在ET反演中不容忽视[14]。本研究针对目前遥感ET模型估算的局限性，利用Landsat 8数据，基于张仁华等的双层模型，考虑地形效应对辐射的影响，进行ET遥感反演，并进行验证和空间分布分析，以期为区域生态水文研究提供借鉴。

1 研究方法

遥感反演思路为：①构建考虑地形特征的辐射计算方法获得太阳辐射通量；②基于Landsat8波段特征构建地表特征参数的计算方法；③将太阳辐射和地表特征参数作为输入数据，利用双层ET模型[6-10]计算水热通量；④借助蒸发比不变法[15]和日净辐射累加值[13]获得日ET。

1.1 考虑地形的太阳辐射计算

受到地形因素控制，起伏地形下坡面太阳辐射由直接辐射Is、反射辐射Ir和散射辐射Id构成。到达任意坡面的太阳辐射强度为：

(1)

式中：i为入射角；是坡度β、坡向α、纬度、赤纬和时角的函数；Ib为垂直到达地球水平表面的太阳辐射，取决于大气透过率τ和大气上界太阳辐射强度；r为地表植被反射率；τd和τr分别为散射和反射辐射透过系数，根据与τ的经验关系求得[16]。

总辐射It计算有2种情况：无地形遮蔽的栅格It(0)和有地形遮蔽的栅格It(1)：

(2)

地形遮蔽度取决于周围地形的相互遮蔽，将最大水平地平线角与同一时刻太阳高度角进行比较，如果大于太阳高度角，则处于阴影范围内，太阳直接辐射为零，否则为正常。对于起伏地形中任意山地网格M，在给定方位角和判断遮蔽范围的半径后，对其沿周边地形沿圆周进行积分，获得各方位角上积分线Li的遮蔽度Vi：

Vi=1-sinai

(3)

式中：ai为Li方向周围山地对M的遮蔽角，基于DEM求解获得；M的遮蔽度V为Vi的均值[17]。

1.2 地表特征参数反演

植被覆盖度、地表反照率的估算方法与landsat 5/7的TM/ETM波段方法类似。地表比辐射率基于改进的Sobrino混合模型[20]，结合NDVI确定的地表分类来确定。地表温度基于Landsat 8-TIRS第10波段数据，借鉴改进型单窗算法[18]反演。

1.3 水热通量和日蒸散发

水热通量基于双层遥感模型获得，详见文献[6-9]。核心算法包括分解地表温度的像元排序对比法和分解地表可供能量的分层能量切割算法。像原排序对比法是根据“地表温度-植被覆盖度”的散点图呈梯形的特点，结合干、湿点定标场的参数测定值，推算4个极端像元的真实温度，从而确定梯形框架，再利用线性内插方法得出地表温度与土壤温度和植被温度的关系，实现组分温度的分解。分层能量切割算法的本质是求算土壤和植被波文比，再利用波文比能量平衡法切割地表可利用能量，进而获得土壤蒸发和植被蒸腾，从而得出潜热通量。日尺度ET是在获得考虑地形因素的地表日净辐射量[13]后，采用较为常用的蒸发比率不变法[15]估算。

1.4 案例研究区概况

研究选择北京和张家口交界处的试验地块(见图1)，处于海河流域北支流域片区，地处太行山和燕山山脉交会处。区域范围大致为北纬40°00′～40°30′，东经115°40′～116°20′，海拔区间为25～1 513 m。主要有山地和平原2种地貌，山地面积占65%。山地主要为林地和灌丛，平原多为农田和建设用地。该区域是官厅水库进京渠道和北京市内流域的上游，是重要的水源涵养区和农产品供应区，构建合理的ET估算方法对监测生态水文过程实现节水灌溉具有重要作用。

图1 研究区图Fig.1 Map of the study area

2 模型验证

研究反演了2013年5月12日、2013年6月13日、2013年7月31日、2013年9月1日、2013年10月3日、2013年11月4日、2013年11月20日，2014年4月29日、2014年10月6日，2015年4月16日10个晴好日的ET。采用研究区怀来地表通量站(见图1)的实测数据验证感热、潜热。怀来地表通量站的大孔径闪烁仪和涡度相关系统测量了2013-2015年每天30 min步长的感热和潜热通量[23-25]，其中2014年4月29日潜热通量缺测，2013年7月31日、2013年10月3日有局部云层干扰，不参与验证。

由表1发现，感热和潜热的相对误差均在15%以内，平均相对误差分别为4.21%和8.36%。图2表明感热通量和潜热通量估算值与实测值决定系数R2分别达0.82和0.98，相关关系较强。感热通量相对误差较低，潜热通量与实测值的相关性更强。

表1 水热通量估算值与实测值对比 W/m2

图2 水热通量遥感估算值与实测值对比Fig.2 Comparison of estimated and measured water and heat fluxes

2013年在北京昌平军都山的缓坡山地架设了波文比实验站，耦合得出日ET。怀来站采用蒸渗仪对2013年的日尺度ET进行了监测[25]。图3为日ET的验证情况。昌平实验点处，决定系数R2为0.93，平均相对误差为8.12%，均方差为0.35 mm/d。怀来站处R2为0.92，平均相对误差为11.09%，虽然样本较小，仍表现出较高的精度。

图3 日蒸散发遥感估算值与实测值对比Fig.3 Comparison of estimated and measured daily

3 结果分析

3.1 太阳总辐射与坡度的关系

太阳总辐射与微地形关系密切，冬季太阳高度角减小，地形效应显著。图4为2013年11月20日的太阳总辐射与坡度的散点图。由图4可以看出，坡度为0～10°的平原地区，地形起伏小，总辐射值相近，散点较为聚集；随着坡度增大，像元散点逐渐稀疏，散点图上下包络线处的像元明显比同一坡度的其他区域更为密集，如图4中聚类1和聚类2所示。聚类1处太阳总辐射与坡度成正相关，主要为向南阳坡；聚类2处太阳总辐射与坡度成负相关，主要为阴坡；聚类3对应太阳总辐射为0～100 W/m2的像元，受地形遮蔽影响较大，无法接收太阳直接辐射，总辐射很低。聚类3出现在坡度30°以上的区域，表明只有当地形起伏效应明显时才体现出显著的遮蔽作用。

图4 太阳总辐射与坡度散点图(2013-11-20)Fig.4 The scatter plot of total solar radiation and slope (2013-11-20)

3.2 ET的空间分布

选择2013年9月1日(2013244)、2013年11月20日(2013324)与2015年4月16日(2015106)作为典型代表日进行分析。图5和图6为通量和日ET的空间及频率分布图。

图5 潜热通量、感热通量、日蒸散发的空间分布图Fig.5 Spatial distribution of latent heat flux, sensible heat flux and daily ET

2015年4月16日感热通量均值为221.35 W/m2，主要集中于180～300 W/m2，20～40 W/m2时有小型波峰，此处为地形遮蔽像元。潜热通量均值为233.17 W/m2，主要集中在140～300 W/m2；日ET均值为3.55 mm，主要集中于1～5 mm，在2～4 mm像元最多。潜热高值区对应水体和阳坡山地，平原和山地阴坡潜热较低，山区地形效应显著，阳坡阴坡差异明显。2013年9月1日感热通量均值为212.19 W/m2，潜热均值为357.53 W/m2，日ET均值为4.82 mm，3者的频率分布没有明显的峰值，高频区主要为山区林地和水体。该时段由于植被覆盖率高，在经历7-8月雨季后，土壤水分供应充足，山地潜热通量大。平原区的潜热通量低于山区，但是农田中夏玉米尚未成熟，蒸散吸收了较多辐射，潜热较高，感热通量相对较低。建设用地为区域内低值区。2013年11月20日的感热均值为109.98 W/m2，潜热通量均值为102.97 W/m2，日ET均值为1.01 mm，甚至出现0值。3个代表日内，日ET的空间分布与潜热类似，感热通量的空间分布与潜热相反。11月20日的植被覆盖率低，潜热、感热和日ET相对较小，频率分布更集中。受植被蒸腾的影响，9月1日的空间差异较另外2 d小。

图6 潜热通量、感热通量、日蒸散发的频率分布图Fig. 6 Frequency distribution of latent heat flux, sensible heat flux and daily evapotranspiration注：黑细线为2013年9月1日；灰色线为2013年11月20日；黑粗线为2015年4月16。

将研究区分为山地(山区林地灌丛的综合体)、城镇绿地、水体、建设用地、裸地和农田6种土地利用类型。图7为不同土地利用类型潜热、感热和日ET的最高、最低和均值统计。

3个代表日潜热通量均值最高的是水体和山地。9月1日由于山地植被覆盖率高，植被蒸腾作用大，山地潜热与水体十分接近，此时农田和城镇绿地的潜热通量均值也较高。11月20日整体潜热都较低，除水体外，其他地类均值相差不大，该时期植被覆盖率低，空气相对湿度较低，夜间气温已到零下会发生凝结，所以ET差异不明显。水体表面存在冰冻，山区林地存在日ET量超过水体极值的局部区域，是该区域积雪经日照迅速融化蒸发或升华的结果。4月16日山地植被出芽不久，农田冬小麦刚开始返青，潜热通量仍处于较低水平。建设用地的潜热通量很低，裸地潜热受干湿影响较大。日ET和潜热呈现类似的分布，3个代表日水体的均值都最高，建设用地和裸地最低。不同代表日各地类的感热通量大小规律相近，水体的感热通量最小，其次是山地，建设用地的感热通量均值最大。从潜热变化幅度看，城镇用地与山地最大，水体与裸土的变化幅度最小，耕地和城镇绿地介于中间。

图7 不同土地利用类型潜热通量、感热通量和日ET统计Fig.7 Latent heat flux, sensible heat flux and daily ET statistics of different land use types

3.3 日蒸散发与地形参数的相关分析

研究区属于半山半平原地形，山区面积占比达65%。为了考察坡度与日ET的关系，以山地为研究对象，以5°为间隔长度，将坡度划分为11级(大于55°为一级)，统计各坡度等级对应像元的日ET均值。以20°为间隔长度，将坡向划分为18级，统计各坡向等级对应像元的日ET均值。日ET与坡度、坡向关系见图8。

图8 日蒸散发与坡度、坡向的关系Fig.8 Relationship between daily ET and slope and aspect

从日ET与坡度关系图可以看出，2013年9月1日ET与坡度拟合曲线更近似于抛物线，R2为0.96。2013年11月20日和2015年4月16日的ET～坡度拟合曲线为线性关系，R2分别为0.92和0.93。在坡度小于35°时，随着坡度上升，ET呈现出较为明显增加趋势，当坡度大于35°时，2013年11月20日和2015年4月16日仍呈现增加趋势，而2013年9月1日则呈现出缓慢下降趋势。3个代表日的ET与坡向的关系均呈反抛物线关系，在0～180°范围内呈负相关，在180°～360°范围内呈正相关。在关系图中可以看出，不论是基于坡度还是坡向2015年4月16日的ET变异程度都是3个代表日中最大的，与该时段植被分布不均匀有关。11月20日植被覆盖率低，因此坡度方面ET变化不明显，而坡向相对明显。9月1日植被覆盖率高，气温也相对高、蒸腾作用大因此坡向和坡度方面的ET变异不大。

4 结 语

研究对双层ET遥感模型增加了考虑地形的辐射计算模块，改进了地表参数计算方法，结合Landsat 8波段特点实现ET反演。经验证，该方法具有较高的精度，感热通量、潜热通量和日ET的平均相对误差分别为4.21%、8.36%和8.12%。

通过分析ET的空间、频率分布以及不同土地利用类型的ET，获得了不同代表日内的ET分布特征；结合地形发现，坡度、坡向对ET的分布有显著影响，在坡度小于35°时，ET与坡度有显著的正相关关系，当坡度大于35°时，受植被覆盖率影响，各季节代表日的ET呈现不同的变化趋势。ET与坡向同样存在相关关系，趋势线程反抛物线。