深圳市前海桂湾公园土壤蒸发及植被蒸腾模拟

张珈玮, 陈 晗, 黄津辉, 李 晗, 蓝至清, 高俊杰

南开大学环境科学与工程学院/中加水与环境安全联合研发中心,天津 300071

蒸散发(Evapotranspiration,ET)是连接水文循环和地表能量平衡的重要组成部分,大约70%的降雨通过蒸散过程返回大气,地表的蒸散发过程消耗50%以上的净辐射能量[1]。之前关于蒸散发的研究多关注农田区域、草地区域和自然林地系统[2—4]。然而,对异质化程度较高的城市区域的蒸散发研究较少。城市化导致植被、土壤、水体等自然地表被不透水区域取代。不透水区域中建筑、道路等高蓄热体增加,降低了潜热通量,增加了显热通量的存储和传递,使得城市区域的气温升高,造成“城市热岛效应”。

城市林地作为城市地表环境的重要组成,其覆盖可占据城市总面积的30%—50%。城市林地可以发挥削弱热岛效应、促进低碳交通、土壤热能利用、吸纳有机垃圾及降低能耗的多重作用[5—7]。虽然城市林地在城市发展中扮演重要角色,但是由于不科学的灌溉措施,维护城市林地所耗费的水资源量极大,造成大量城市水资源浪费。在此背景下,有必要拆分城市林地蒸散发组分,准确估计植被蒸腾及土壤蒸发相对比例,获取真实植被耗水量信息,从而更科学的指导城市林地精准灌溉决策。深圳市前海桂湾公园位于高度开发的前海开发区中央,是深圳市前海地区重要的公共空间和集中绿地,准确估算深圳市前海桂湾公园的土壤蒸发和植被蒸腾,为园区实现智慧园林管控提供理论依据。

随着遥感技术的发展,结合高分辨率地物光谱信息与地面气象观测数据并基于能量平衡和湍流扩散理论,反演区域尺度蒸散发变化成为可能[8]。基于遥感的蒸散发模型大体可以分为单源模型和双源模型[9]。双源模型不仅能够区分土壤蒸发和植被蒸腾,模型精度也显著优于单源模型。Shuttleworth-Wallace(S-W)模型和双作物系数法(FAO dual-Kc)被认为是模拟区域土壤蒸发和植被蒸腾的标准模型。其中S-W模型是基于彭曼公式改进的双源模型,模型引入5个阻抗参数来描述土壤蒸发和植被蒸腾的湍流扩散过程[10],其性能在不同植被下垫面得到广泛验证。FAO dual-Kc模型是由联合国粮食和农业组织(FAO)开发的一种间接估算蒸散发的方法,并在FAO灌溉和排水第56号文件中进一步发展为可分别估算土壤蒸发和植被蒸腾的双作物系数法[11—12]。然而截至目前,基于遥感的城市蒸散发模型大多基于单源模型机制,如Faridatul等人[13]开发的城市地表能量平衡算法(USEBAL)模型、Wang等人[14]的经验城市(EU)模型和Qiu等人[15]的三温(3T)模型。这些模型无法单独估算城市区域土壤蒸发和植被蒸腾。关于目前可用的城市蒸散发模型,共同存在的问题是：1)尚未有研究者针对城市林地区域开发一个遥感双源蒸散发模型；2)适用于城市区域的遥感蒸散发模型分辨较低,不能很好的捕捉异质性程度较高的城市下垫面状况；3)城市区域的净辐射通量(Rn)和土壤热通量(G)没有合理的参数化方案。

MOD16(MODIS全球蒸散发产品)双源蒸散发模型基于经典的Penman-Monteith公式改进而来,使用互补相关理论结合相对湿度估算土壤蒸发。此外,模型假设植被冠层阻抗受到空气温度和水汽压差的约束,通过植被覆盖率对净辐射进行分配,以实现植被蒸腾的估算。目前MOD16算法已经广泛的应用于各类地表类型的土壤蒸发及植被蒸腾模拟。由于MOD16模型不需要输入遥感热红外波段即可实现区域蒸散发的模拟,因此MOD16模型被认为是可以实现城市蒸散发模拟的潜在模型。但是,该模型并不能直接应用在城市区域,由于存在以下局限性：1,MOD16模型应用1 km空间分辨率的MODIS的遥感数据,不能很好的反映异质化程度较高的城市下垫面类型；2,MOD16模型多用于自然区域下垫面植被蒸腾和土壤蒸发的模拟,而城市林地下垫面与自然下垫面具有截然不同的能量吸收和土壤热通量传递过程。因此,直接将MOD16模型应用于城市下垫面会带来较大误差。

针对现有城市蒸散发模型和MOD16模型在城市林地方面应用存在的问题,本研究的目标包含以下几点：1)对目前广泛使用的MOD16模型的物理机制进行改进,针对城市林地下垫面类型,提出了改进的净辐射通量和土壤热通量算法,改进后的模型可以更加准确的描述城市林地区域复杂下垫面的能量分配过程；2)基于10 m空间分辨率Sentinel- 2遥感影像,反演2020年深圳市前海桂湾公园20个无云日的土壤蒸发和植被蒸腾；3)基于S-W模型和FAO dual-Kc验证模型性能；4)识别所开发模型的关键敏感性参数并分析深圳市前海桂湾地区土壤蒸发和植被蒸腾空间分布格局。

1 研究区域与研究方法

1.1 研究区概况

本研究区域位于深圳市南山区前海桂湾公园(图1),其纬度为22°31′N—22°32′N,经度为113°53′E—113°54′E。深圳市属于亚热带海洋性气候,夏季高温多雨,其余季节气候温和,较为干燥,年平均气温在22.4℃左右。年平均降雨量为1935.8毫米,其中86%发生在雨季(4月—9月)。前海桂湾公园占地52万m2,是前海地区第一个水廊道公园,更是一带一路国际产能合作中引领片区产业发展的展示窗口。园内主要为灌木和常绿乔木组成的红树林湿地木本植物群落,绿地率达到72%。

图1 研究区域

1.2 研究方法及数据处理

由于该地区缺少地表ET涡动相关通量监测数据,本研究使用两类广泛使用的双源模型(S-W和FAO dual-Kc模型)验证改进的MOD16模型性能。选取一年四季中20个无云或少云日,计算得到卫星每日过境的瞬时蒸散量,采用蒸发比不变法将瞬时LE值转换为日值[16—17]。对于驱动S-W模型和FAO dual-Kc模型所需要的土壤水含量,每隔4—6 d,采用烘干法测定园区内6个采样点处深度1 m、间隔0.1 m土壤样品。因为深圳4、5、6月份为梅雨季节,且桂湾公园位于海湾附近,此时园区内水分供应充足,FAO dual-Kc模型中的土壤水分胁迫系数Ks取1。其余月份按照根区水分损耗Dr大于速效水RAW情况进行取值。FAO dual-Kc模型中植被基础系数(Kcb)和土壤蒸发系数(Ke)参照FAO灌溉和排水第56号文件[11]选取。

输入改进后的MOD16模型数据源包括Sentinel- 2遥感影像和气象数据因子。其中,Sentinel- 2影像下载于欧洲航天局(ESA)官网,下载地址为https://scihub.copernicus.eu 。Sentinel- 2遥感影像包含Level-0、Level- 1A、Level- 1B、Level- 1C数据。本研究中使用未经过大气校正的L1C级别数据,根据欧空局提供的预处理插件Sen2Cor(02.08.00版本),使用Windows命令配置调用插件,处理得到大气校正后的数据；将大气校正后的数据导入SNAP软件中重采样至10 m分辨率；之后将重采样数据在ENVI 5.3软件中进行波段融合,裁剪等处理。气象观测数据下载于中国气象数据网,包括日平均气压、日平均气温、日平均水汽压、日平均相对湿度、净辐射日总量、日平均风速等,原始下载地址为http://data.cma.cn/。为保证对比验证的可靠性,驱动S-W模型和FAO dual-Kc模型的地表气象监测数据来源与改进的MOD16模型的地表监测数据来源相同,均为中国气象数据网。

1.3 改进的MOD16模型构建

本研究基于MU等人[18]提出的MOD16双源模型,并针对城市林地区域的下垫面特点改进而来。在原始MOD16模型的基础上,本研究针对城市林地下垫面特点,提出了新的Rn和G的估算方法。MOD16模型关于土壤蒸发和植被蒸腾的参数化方案如等式(1)和(2)所示：

(1)

(2)

式中,LEv(W/m2)、LEs(W/m2)分别是植被潜热通量和土壤潜热通量；Δ(Pa/K)是饱和水汽压相对于温度的斜率；Rn(W/m2)是净辐射通量；G(W/m2)是土壤热通量；fc为植被覆盖度,具体计算方法参见Gutman等人[19]；ρ(kg/m3)是空气密度；Cp(J kg-1K-1)是空气定压比热容；VPD(hPa)饱和水汽压差是饱和水汽压和实际水汽压之间的差值；γ(Pa/K)是干湿表常数；ra(s/m)是从地表到参考高度处的空气阻抗；rv(s/m)是冠层气孔阻抗；rtot(s/m)空气动力学和水汽扩散阻抗；RH(%)是相对湿度；β是经验参数,根据Fisher等人研究[20]将其设置为1000。

改进的净辐射通量Rn参数化方案如(3—4)式：

Rn=(1-α)Ra↓-σεTa4+εRl↓

(3)

式中,Ra↓(W/m2)和Rl↓(W/m2)分别是太阳短波辐射和大气长波辐射；α是地表反照率；σ是斯忒藩-玻尔兹曼常数,通常取5.67×10-8W m-2K-4；Ta(K)是空气温度；ε是四种城市表面辐射率,依据不同的土地利用和土地覆盖情况,ε的计算方法如下[13]：

义7-6块直井长缝日均产液量为30.3t，日均产油量为18.8t；常规压裂日均产油量为12.9t；直井长缝压裂日均产油量为常规压裂日均产油量的1.92倍。

(4)

式中,NDVI是归一化植被指数。

其次,根据比尔-朗伯定律,将Rn分为土壤层净辐射通量(Rns,W/m2)和植被层净辐射通量(Rnv,W/m2),Rns和Rnv是叶面积指数(LAI)的函数[21—22],计算公式如下：

Rnv=Rn(1-exp(-kLAI))

(5)

Rns=Rnexp(-kLAI)

(6)

式中,k为消光系数,根据Impens和Lemur的研究[23]将其定为0.6；LAI为植被覆盖度(fc)的函数：

(7)

式中,kpar是经验系数,设为0.5[24—25]。

改进的土壤热通量估算方 法如(8)式所示：

(8)

式中,cg是非植被覆盖系数,取决于地表热容量和地表导热率[26]。

1.4 敏感度分析

为了评估不同输入参数对模型结果的影响程度,通过敏感性分析观察LEv和LEs的变化程度。输入参数的敏感度计算公式如下[27—28]：

(9)

式中,Si是模型对输入参数i的敏感度(%),LE±为变量增加或减少后模型输出的LEv(W/m2)和LEs(W/m2)。LE0是未改变变量时的LEv和LEs。根据输入改进的MOD16模型参数,重点关注的参数为气象参数和植被参数,设置(-20%,20%)区间内5%区间间隔。其中气象参数有气温(℃)、相对湿度(%)、水汽压差(hPa)和净辐射能(W/m2)；植被参数为叶面积指数(LAI)和植被冠层阻抗(rv)。考虑到研究区Sentinel- 2影像产品质量和各个参数的变化范围,选取2020年7月23的数据进行敏感性分析。

2 结果分析

2.1 模型模拟精度验证

改进的MOD16模型模拟LE与两类双源模型验证散点图如图2所示,相关统计指标见表1和表2。在20个测试日中,改进的MOD16与S-W和FAO dual-Kc两类双源模型模拟LE值验证的均方根误差(RMSE)分别为21.39 W/m2和20.41 W/m2,平均绝对误差(MAE)分别为18.81 W/m2和19.05 W/m2,相对偏差(Bias)分别为-16.71 W/m2和-12.05 W/m2。改进的MOD16与两类模型的验证性能相差较小,均在允许的误差范围之内。然而,从相关系数上来看,改进的MOD16模型模拟LE值与S-W模型得到的参考值的相关度为0.801,显著高于与FAO dual-Kc模型参考值的相关系数0.634。

图2 改进的MOD16模型和Shuttleworth-Wallace(S-W)模型、双作物系数法(FAO dual-Kc)计算的LE散点图

表1 改进的MOD16模型与Shuttleworth-Wallace(S-W)模型潜热通量LE差异

表2 改进的MOD16模型与双作物系数法(FAO dual-Kc)潜热通量LE差异

图3 改进的MOD16模型和Shuttleworth-Wallace (S-W)模型、双作物系数(FAO dual-Kc)模型计算的LEv/LE散点图

2.2 蒸散量时空变化特征

2.2.1研究区蒸散量时间变化特征

由图4可以看出,2020年整年,深圳市前海桂湾公园地区总蒸散月均值在85—165 W/m2之间。研究区月均总蒸散量在一年中呈现出单峰变化趋势。总体趋势为1—7月呈现上升趋势,8—12月呈现下降趋势。研究区总蒸散量月均最高值为162.06 W/m2,出现在夏季的七月份,占全年蒸散量的10.5%；总蒸散量月均最低值为87.98 W/m2,出现在冬季的一月份,占全年蒸散量的5.7%。

图4 2020年研究区净辐射量及潜热通量(LE)时间变化

2.2.2研究区蒸散量空间变化特征

选取2020年5月4日(春季)、9月6日(秋季)和11月25日(冬季)3个具有代表性的测试日(图5)。从图5可以看出：1)2020年研究区的土壤蒸发(LEs)和植被蒸腾(LEv)空间分布差异显著,土壤蒸散量高值(60—76 W/m2)分布区域植被蒸散量相对较低(22—58 W/m2),植被蒸散量高值(92—147 W/m2)分布区域土壤蒸散量相对较低。植被蒸散量高值多分布在前海桂湾沿岸,其沿岸植被覆盖率较高,以灌木和常绿乔木植被为主,低值出现在研究区西北和东北部。土壤蒸散量高值主要分布在研究区的西北、东北和东南部有较多裸地的区域。土壤蒸散量高值和植被蒸散量高值空间分布格局互补。2)春季、秋季和冬季蒸散量的空间分布一致。在不同季节中土壤蒸散量和植被蒸散量高值区域分布格局不变,没有随季节变化显示出空间分布差异。季节对研究区蒸散量空间分布几乎没有影响。3)LEv高值在蒸散量为零的不透水区域周边也有分布,主要因为在这些不透水区域周围进行了城市绿化,有较多树木和草地分布。总体来说,研究区夏季蒸散量达到最高值174.09 W/m2,秋冬蒸散量开始减少,LEv高值呈现零散分布,多集中在前海桂湾沿岸,LEs高值集在研究区的西北、东北和东南部。

图5 2020年不同季节植被蒸腾(LEv)和土壤蒸发(LEs)空间分布图

2.3 改进的MOD16模型敏感度分析

从表3、表4和图6可以看出,净辐射、水汽压差和植被覆盖度对植被蒸散量(LEv)显示出较强的正相关性,其中净辐射和植被覆盖度影响最大。随着净辐射和植被覆盖度5%、10%、15%到20%的变化,导致LEv增加3.07%、6.15%、9.22%和12.30%。温度显示出较弱的负相关性,随着温度5%、10%、15%到20%的变化,导致LEv减少1.29%、2.51%、3.66%和4.75%。而在本研究区域冠层阻抗的变化对LEv的影响甚微。因为研究区内存在较多的不透水表面,削弱了冠层阻抗这一参数的影响。对于土壤蒸散量(LEs),净辐射、水汽压差和湿度显示为正相关性,其中湿度影响最强。随着湿度5%、10%、15%到20%的变化,导致LEv增加6.03%、12.13%、18.29%和24.49%。温度和植被覆盖度对LEs呈现为负相关性,植被覆盖度对LEs有较强的抑制作用。当植被覆盖度分别增加5%、10%、15%和20%时,LEs减少3.09%、6.18%、9.28%和12.37%。随着植被覆盖的增加,地表接受到的太阳辐射能受到抑制,进而影响LEs。

表3 改进的MOD16模型估算LEv值对输入变量的敏感度

表4 改进的MOD16模型估算LEs值对输入变量的敏感度

图6 改进的MOD16模型估算LE对输入变量的敏感度

3 讨论

本研究通过改进能量通量,提出了适用于城市林地的双源模型。改进后的双源模型与S-W和FAO dual-Kc两类模型的统计学指标均在误差允许范围内。在总蒸散量(LE)的估算上,改进的MOD16模型和S-W的相关性显著高于和FAO dual-Kc模型的相关性。这可能是由于在FAO dual-Kc模型中,双源系数由植物生长周期确定,该系数与当地气候、植物类型与植被生长状况密切相关[29—30],直接采用广义作物系数给土壤蒸发和植被蒸腾模拟带来较大不确定性和误差。改进的MOD16模型能够很好的感知研究区蒸散量随季节改变的变化,总蒸散量最高值174.09 W/m2出现在夏季,最低值72.38 W/m2出现在冬季。这与Weng等人研究发现城市内潜热通量范围为96—195 W/m2的结果相一致[31]。原因是冬季研究区太阳辐射能低、气温低,加之此时不是深圳市的雨季降雨量少,不利于蒸散发；春季气温回升,植被开始生长,蒸散量开始升高；夏季时太阳辐射最强,气温升高,且是深圳市的雨季,为蒸散发提供了有利条件；转入秋季后,气温降低,植被进入枯萎期,蒸散量也开始降低。因此,建议在蒸散量较高的春夏季增加灌溉,蒸散量较低的秋冬季减少灌溉。由图4也可以看出蒸散量变化与净辐射和气温变化有较高的一致,这也证实了Faridatul等人认为太阳净辐射是影响城市区域蒸散发重要因素的结论[13]。

虽然改进的MOD16模型可以很好的估算城市林地这样小尺度区域的蒸散量,但是模型框架仍存在短板。阻抗参数作为MOD16模型关键参数,合理的参数化方法能够显著提升模型精度。本研究中,将应用于自然林地的阻抗参数直接用来估算城市林地的蒸散量。然而,城市林地的下垫面构成复杂、异质性程度高。因此,将均匀下垫面阻抗参数化方法直接应用于城市林地非均质表面,可能会给计算结果带来一定的误差。另外值得注意的是本研究区缺少涡动观测数据,不能将模型表现与实际观测结果进行验证。虽然针对观测数据不足,对改进的双源模型展开了敏感度分析,但是模型验证对于利用遥感数据进行蒸散量估算仍具有重要意义。针对以上问题,未来的研究还需要考虑以下两点1)构建植被、土壤和不透水面耦合关系下的非均质下垫面的阻抗参数；2)结合涡动观测数据率定模型参数,进一步提高模型精确度。

4 结论

本研究通过改进土壤热通量(G)和净辐射通量(Rn),提出了适合城市林地区域的改进的MOD16模型,并对深圳市前海桂湾地区2020年土壤蒸发和植被蒸腾的时空演变及其受输入参数的影响进行分析,得到以下结论：

1)在20个研究日中,改进的MOD16模型与参考模型验证结果总体精度较好。改进的MOD16与S-W和FAO dual-Kc两类双源模型模拟LE值验证的均方根误差(RMSE)分别为21.39 W/m2和20.41 W/m2,平均绝对误差(MAE)分别为18.81 W/m2和19.05 W/m2,R2分别为0.801和0.634；与两类模型估算LEv/LE值的RMSE分别为0.036和0.023,MAE分别为0.030和0.019,相关系数R2分别为0.539和0.575。该模型能够提供小区域高分辨率的蒸散量估算,有效适用于城市林地区域。

2)城市林地区域蒸散量具有明显的时间和空间分布差异。具体体现在：城市林地总蒸散量月均值在85—165 W/m2之间,春季到夏季呈现上升趋势,秋季到冬季呈现下降趋势；植被蒸散量高值在研究区域中零散分布,较多出现在前海桂湾公园,土壤蒸发与植被蒸腾展现出相反的分布特征。

3)敏感性分析表明,净辐射和植被覆盖度对LEv模拟结果的影响最大,其次是水汽压差和大气温度,而冠层阻抗对城市林地LEv模拟的影响很小；LEs的模拟结果受湿度和植被覆盖的影响最大。