基于SEBAL模型的若尔盖高原蒸散量时空变化

张 亮，崔林林，苑 跃

(1.高原与盆地暴雨旱涝灾害四川省重点实验室，成都 610072；2.四川省气象灾害防御技术中心，成都 610072； 3.成都信息工程大学资源环境学院，成都 610225)

前 言

地表蒸散(Evapotranspiration, ET)是全球水循环的重要组成部分，蒸散量每年约占陆地降水量的60%[1]，并在区域水平衡中发挥着重要作用。定量研究不同时空尺度下蒸散的分布，对于全球气候变化、生态环境保护、水资源管理和农业生产都具有重要意义。

目前，蒸散的估算方法主要分为传统法和基于卫星遥感的方法。传统方法主要有蒸渗仪法、波文比法-能量平衡法、空气动力学法、涡度相关仪法等，这些方法能够在均匀的地表条件下提供较为精确的计算结果，但存在着理论模型复杂并且只能应用于小尺度区域等缺陷[2-3]。近年来随着卫星遥感技术的发展，基于卫星遥感数据的蒸散估算方法可对大范围区域的蒸散量进行计算[4～6]。基于地表能量平衡理论的SEBAL(Surface Energy Balance Algorithm for Land)模型是当前应用较多的模型之一，在蒸散的计算过程中具有以下优势：(1)最小程度的使用地面观测数据；(2)模型中近地表的空气温度的非强制性；(3)可通过识别干湿像元以及测定近地表空气温度差进行自校正[7～9]。研究表明，SEBAL模型蒸散计算结果在日尺度、季节尺度和年尺度上都具有较高的精度[9]。近年来国内外学者利用SEBAL模型对不同生态系统的地表蒸散量进行了大量研究，如：Opoku-Duah et al.利用ASTER和MODIS数据对西非热带草原地区的蒸散量进行估算，并将结果与Landsat ETM+图像中提取的ET进行了比较[10]；Teixeira et al.对SEBAL模型进行了进一步改进，并对巴西东北部旧金山河流域的蒸散量进行了估算[11]；程明瀚等对北京市日蒸散区域分布研究后发现，北京市蒸散整体呈现出西北山区高、东南平原低的趋势[12]；周妍妍等对疏勒河流域蒸散量的时空变化研究后发现，且其空间分布异质性十分明显，整体上从东南向西北呈逐渐减小趋势[13]；李宝富等估算了塔里木河干流区蒸散量，结果表明该区日蒸散量较高，其中水体蒸散最高，耕地和林地次之，居民和工业用地蒸散量最小[14]。

若尔盖高原地处青藏高原东北部，区内河网密集，是黄河上游重要水源补给区，被誉为“中国西部高原之肾”，同时也是全球气候影响的敏感区。近几十年，由于受人口增长、气候变化、过度放牧等影响，该区域的土地利用/覆被发生剧烈变化，沙丘面积逐步扩大，进而导致区域水文过程(蒸散、径流、截留等)和水循环过程的改变[15～17]。前人主要利用彭曼公式(Penman-Monteith formula)和MODIS的蒸散产品对若尔盖高原的潜在蒸散量进行计算，但得到的结果与站点观测值相差较大[18～20]。因此，本文利用SEBAL模型对近20年来若尔盖高原的实际蒸散量进行估算，并对其时空变化特征进行分析，研究结果可为本地区的气候变化、水资源管理和生态环境保护提供科学支撑。

1 研究地区与方法

1.1 研究区概况

若尔盖高原位于青藏高原东北部边缘，海拔介于3 400～3 900 m 之间，隶属四川省若尔盖县、红原县和阿坝县，以及甘肃省的玛曲县和碌曲县(图1)。若尔盖高原为大陆性高原气候，寒冷湿润，长冬无夏，霜冻期长，日温差大。研究区多年平均气温1℃左右，年日照时数均在2 000h以上，全年降水量 600～800mm，雨季分布在5～10月，约占全年降水量的 90%。本区植被类型以高山草甸和沼泽植被为主，区内河流多属于黄河水系，是黄河上游重要的水源补给区。

图1 研究区域位置图Fig.1 Location of research area

1.2 数据及方法

本研究所用若尔盖高原及周边299个国家及区域自动气象站逐日数据，来源于四川省气象局，包括气温(℃) 、平均风速(m/s) 、日照时数(h)、相对湿度(%) 和蒸散(mm)等，研究时段为2001～2017年，所有数据均经过质量控制，并通过回归插值法对缺测数据进行补充。遥感数据为2001～2017年EOS/MODIS数据，主要包括1 km分辨率的地表温度产品(MOD11A1)、1 km分辨率的植被指数产品(MOD13A2)、1 km分辨率的土地利用产品(MCD12Q1)、500 m分辨率的反照率产品(MCD 43A3)和1 km分辨率的反射率产品(MOD09GA)等，下载于美国国家航空和航天局的地球观测系统数据和信息系统(http://reverb.echo.nasa.gov/reverb/)。所有的输入数据均预处理为Albers等面积投影和1 km空间分辨率的数据集。

区域蒸散采用SEBAL模型进行计算，SEBAL模型基于地表能量平衡方程，通过卫星遥感与气象数据分别对各能量分量进行估算，最后得到区域的实际蒸散值。地表能量平衡方程为：

λET=Rn-G-H

(1)

公式(1)中λ为潜热蒸发系数，取值2.45×106 MJ/kg，λET是潜热通量(W/m2)，Rn是地表净辐射通量(W/m2)， G是土壤热通量(W/m2)，H是显热通量(W/m2)，Rn、G和H的计算方法介绍如下。

地表净辐射通量Rn根据公式(2)进行计算。

Rn=(1-α)RS↓-RL↑+εaRL↓

(2)

公式(2)中： RS↓是下行到达地表的太阳短波辐射(W/m2)； RL↓是下行的长波辐射(W/m2)； RL↑是上行的长波辐射(W/m2)；α为地表反照率，εα为地表比辐射率，具体的参数计算过程如下。

RS↓=GSCcosθdrτSW

(3)

(4)

(5)

(6)

ε0=1.0094+0.0047lnNVDI

(7)

τSW=0.75+2×10-5z

(8)

公式(3)～(8)中: GSC为太阳常数，1 367 W/m2;θ为太阳天顶角; dr为日地距离因子，可通过儒略日求得; τsw为大气单向透射率，根据研究区内气象观测站的海拔获得; z为高程; σ为5.67 × 10-8W/(m2·K4)；Ts为地表温度(K); Ts_cold为“冷”点地表温度(K)。α为地表反射率; αtoa为大气层顶反射率; αpath_radiance为大气层辐射值，NVDI为归一化植被指数。

土壤热通量G根据公式(9)进行计算。

(9)

显热通量H根据公式(10)进行计算。

(10)

(11)

(12)

公式(10)～(12)中：ρ为空气密度(kg/m3)；Cp为空气定压比热，取值1004J/(kg·K)；dT为地面Z1、Z2两个高度的温差，Z1和Z2是参考高度植被的零平面位移；rah为空气动力学阻抗(s/m)。u*是摩擦速度(m/s); k是卡门常数，取0.4; u200是200 m 的摩擦风速(m/s)；zom是地表动量传输的粗糙度长度(m)。其中，dT与rah为相互关联的2个未知量，直接求解比较困难，SEBAL 模型引入Monin-Obukhov循环迭代计算，直至感热通量稳定。

由于公式(1)～(12)计算的是卫星过境瞬时地面蒸散值(Λinst)。研究表明，一天内蒸发比(Λ24)随时间不发生变化，即Λ24=Λinst。因此利用蒸发比的这一特性对蒸散进行时间尺度的扩展，计算出一天的总蒸散量(式13和式14)。

(13)

ET24=86400Λ24(Rn24-G24)/λ

(14)

公式(13)(14)中：ET24为日蒸散量，Rn24和G24分别为日净辐射通量和日土壤热通量，其中日土壤热通量G24由于白天和夜晚保持能量收支平衡，可以忽略不计。日净辐射通量Rn24计算公式如下。

Rn 24= (1-α)Ra 24τSW- 110τSW

(15)

(16)

ωs=arccos(-tanφtanδ)

(17)

公式(15)～(17)中：φ为纬度，δ为太阳赤纬。

2 结果分析

2.1 SEBAL模型估算结果的对比验证

本文选取2015年的若尔盖、马尔康、红原、阿坝4个站点的蒸散实测数据(E-601B型蒸发器)与SEBAL模型估算出的蒸散量进行对比验证。如表1所示，SEBAL模型的估算结果与四个站点的实测值的相关性较高，相关系数在0.53～0.75之间，相对误差为25.05%～33.53%之间。由此可见，SEBAL模型在本研究区有较好的适用性。

表1 站点蒸散实测值与SEBAL模型估算结果对比Tab.1 The results comparation between SEBAL model and measured values by meteorological station

2.2 研究区蒸散量的时空变化特征

利用SEBAL模型估算研究区2001～2017年平均单位面积蒸散量，并统计出不同年份平均蒸散量的方差和变异系数于表2中。可以看出，17年间研究区平均单位面积蒸散量介于546.32～627.37mm，并且单位面积的蒸散量变异性较小，位于4.82%～6.67%之间。对于不同行政区域而言，玛曲县(5.0%)和红原县(4.78%)的单位面积蒸散量变异性最大、阿坝县居中(3.95%)、碌曲县(2.96%)和若尔盖县(2.68%)最小。拟合2001～2017年蒸散量的多年变化趋势后发现，研究区单位面积蒸散量总体上呈微弱波动上升趋势(图2)。

为了分析2001～2017年研究区蒸散量的空间分布格局，将单位面积蒸散量按照400～500mm、500～600mm、600～700mm和700～800mm划分为4个等级(图3)，从图3中可知，同一年份不同区域的蒸散量具有显著的差异性，且不同年份蒸散量的空间分布也具有差异性。总体上看，研究区单位面积蒸散量从西北向东南呈现出逐渐增加趋势。对于不同行政区域而言，17年平均单位蒸散量从大到小依次为红原县(622.11mm)、玛曲县(587.49mm)阿坝县(580.29mm)、碌曲县(572.70mm)、和若尔盖县(565.50 mm)。

表2 2001～2017年研究区单位面积蒸散量年均值统计Tab.2 Statistics of average annual evapotranspiration per unit area of Ruoergai Plateau

图2 若尔盖高原单位面积年均蒸散量时间变化趋势Fig.2 Change trend of average annual evapotranspiration per unit area of Ruoergai Plateau

图3 2001～2017年研究区蒸散量空间分布格局 Fig.3 Distribution pattern of evapotranspiration of Ruoergai Plateau from 2001 to 2017

2.3 研究区蒸散量时间变化趋势的空间分布特征

基于2001～2017年蒸散数据，逐像元计算蒸散量变化的斜率，并分为7个相应的蒸散变化等级：显著减少(slope<-10)、中等减少(-1010)(图4)。从图4中可以看出，研究区大部分都属于轻微增加等级，部分中等增加等级的斑块点缀其间，靠近东南边界有少量轻微减少等级的区域。轻微减少、轻微增加和中等增加区域所占的面积比例分别为2.11%、97.20%、0.68%，其余等级区域所占面积几乎为0。从行政区域上讲，研究区内五县均属于轻微增加趋势。

图4 研究区2001～2017蒸散量年际变化趋势Fig.4 Spatial distribution map of inter-annual evapotranspiration trend in Ruoergai Plateau from 2001 to 2017

2.4 不同植被类型下的蒸散量变化特征

本研究选择了研究区内三种主要的植被类型(湿地、草地和林地)作为分析对象，基于土地利用数据和蒸散数据，统计2001～2017年单位面积的平均蒸散量(表3)。从表3中可以看出，不同植被类型单位面积蒸散量具有差异，且同种植被在不同年份的单位面积蒸散量也变化较大，具体表现为湿地>草地>林地，17年单位面积蒸散量分别为598.17mm、587.46mm和549.83mm。分析三种植被类型单位面积蒸散量的年际变化趋势后发现(图5)，湿地、草地和林地单位面积蒸散量均呈现微弱增加趋势，且三者表现出同步变化规律。

表3 2001～2017年不同植被类型下单位面积蒸散量Tab.3 Statistics of average annual evapotranspiration per unit area for different vegetation types (mm)

图5 2001～2017年三种植被类型下 的单位面积蒸散量变化趋势Fig.5 Variation tendency of average annual evapotranspiration per unit area for different vegetation types from 2001 to 2017

3 讨 论

本文利用SEBAL模型计算得到2001～2017年若尔盖高原平均单位面积实际蒸散量介于546.32～627.37mm之间，小于刘佳等[18]和王建兵等[19]利用传统FAO56 P-M方法计算的若尔盖高原蒸散量，分别为1 200～1 600 mm和700～820 mm。这是由于FAO56 P-M法计算的为潜在蒸散量(ET0)，不能代表实际蒸散量(ET)。李志威等[17]将各类下垫面面积作为权重因子，利用FAO56 P-M 法加权计算若尔盖高原的实际蒸散量，结果表明1961～2011年若尔盖高原平均单位面积实际蒸散量介于452.6～535.2 mm之间。该结果略小于本文计算结果，但量级一致，造成这种差异的原因可能为研究方法及时空尺度不同。本文2001～2017年的多年平均蒸散量为587.06mm，与叶红等[20]利用改进的MODIS ET算法得到的2000～2014年多年平均蒸散量538.61mm类似，这可能是由于所用数据一致和研究年份相近的缘故。

研究区2001～2017年蒸散量整体上呈波动微弱增加趋势，其中2001～2008年蒸散量呈降低趋势，2009～2017年蒸散量呈波动增加趋势(图2)，这可能与研究区实施生态保护密切相关[21-23]。同时与叶红等[20]得到的结果类似：研究区蒸散量在2002、2008、2010和2013年处于低谷，且时间段整体和分段趋势一致。针对不同的土地利用类型，蒸散量从大到小依次为湿地>草地>林地，这也与文献[24-25]的研究结论相似。在模型参数设置方面，相比于传统FAO56 P-M方法，SEBAL模型能更好的对不同植被类型下的实际蒸散量进行计算。此外，本文结果从侧面印证了若尔盖高原湿地出现了明显的暖干化趋势。由于若尔盖高原是黄河上游重要的生态屏障，同时也是气候变化的敏感区和脆弱带，本文侧重于若尔盖高原实际蒸散量时空变化特征分析，后续工作将对影响若尔盖高原实际蒸散量的气象因子进行研究。

4 结 论

本文基于遥感和气象台站数据，利用SEBAL模型反演了若尔盖高原的蒸散量，并采用统计分析方法定量评价了区域蒸散的时空变化特征，得出主要结论如下。

4.1 SEBAL模型的反演结果与地面实测数据的相关系数介于0.53～0.75之间，计算误差小于35%，说明SEBAL模型在若尔盖高原具有较好的适用性。

4.2 2001～2017年研究区单位面积年均蒸散量介于546.32～627.37mm，变异性较小(<7%)。从时间上看，单位面积年均蒸散量呈波动上升趋势，其中轻微减少、轻微增加和中等增加等级区域所占的面积比例分别为2.11%、97.20%、0.68%，其余等级区域所占面积几乎为0%。从行政区域上讲，研究区内五县均属于轻微增加趋势。

4.3 从空间上看，研究区单位面积蒸散量呈现出从西北向东南逐渐增加趋势。不同行政区域而言，17年平均单位面积蒸散量从大到小依次为红原县(622.11mm)>玛曲县(587.49mm)>阿坝县(580.29mm)>碌曲县(572.70mm)>若尔盖县(565.50 mm)。

4.4 研究区不同植被类型(草地、林地和湿地)单位面积蒸散量均呈现微弱增加趋势，且三者表现出同步变化规律，17年平均单位面积蒸散量表现为湿地(598.17mm)>草地(587.46mm)>林地(549.83mm)。