CN114417553A - 一种基于遥感的陆表蒸散计算方法 - Google Patents

Abstract

本发明属于地理学与生态学领域，涉及一种基于遥感的陆表蒸散计算方法。陆表蒸散由植被蒸腾、土壤蒸发与拦截蒸发组成。本发明主要对植被蒸腾中的温度和土壤蒸发中的土壤水分进行了优化，主要包括两部分：使用物理模型计算叶片温度代替原方法中环境温度进行植被蒸腾计算；使用不同土壤层的土壤含水量来重新定义土壤蒸发过程中的土壤水分限制。本发明的方法简单实用，所涉及的数据均可从遥感平台、中国气象数据网、涡旋通量塔以及相关文献免费获取。本发明属于机理模型，外推时无需参数校正，可以更低的成本快速地对陆表蒸散进行可靠的估算。

Description

一种基于遥感的陆表蒸散计算方法

技术领域

本发明属于地理学与生态学领域，涉及一种陆表蒸散的计算方法。具体使用为通过遥感技术获得相关MODIS与FLDAS数据、通过国家气象中心获取站点相对湿度数据、通过中国通量塔及相关文献获取蒸散数据，并基于现有的Priestley-Taylor模型框架，优化其土壤水分子模型和温度子模型，将优化后的模型用于陆表蒸散的估算。

背景技术

蒸散(Evapotranspiration，E)是陆表水热循环的重要组成部分，由植被蒸腾、土壤蒸发与拦截蒸发三部分组成。在最适环境下的最大蒸散量叫做潜在蒸散；在受到气候、植被与土地覆盖等实际状况限制下的蒸散量叫做实际蒸散。超过50％的太阳辐射与60％的年总降水都是通过实际蒸散返回到大气中。

蒸散与可利用能的关系可以表示为：

LE＝R_n-G-H

其中，L表示蒸发潜热(MJ/kg)，E表示蒸散(mm)，LE表示潜热通量(W/m²)，R_n表示净辐射(W/m²)，G表示土壤热通量(W/m²)，H表示显热通量(W/m²)。

涡旋通量塔是观测蒸散最好的方式，但是全球仅有400个通量塔，记录时间短，费用高，且很难从点尺度扩展到面尺度。想要获取空间上连续的蒸散，一种方法是通过加大通量塔密度并插值，但是这种方法成本高昂，费时费力，且只能获取建立通量塔以后的数据，不具有较大的可行性；另一种方法是通过反距离权重与克里金、样条函数等方法将气象数据进行插值，然后基于气象数据的栅格数据与Penman-Monteith等模型得到蒸散，但是插值的前提是地形平坦，空间异质性较小，蒸散所需的各个要素都有足够多且连续的站点记录数据。实际上，在一些发展中国家和地区，气象站数量较少，时间跨度较短，插值效果受到很大的限制，特别是在青藏高原地区，地形与气候变化极大，并且在青藏高原西北部缺少气象站，插值效果会更差。因此，长时间序列高时空分辨率的免费遥感数据便成了估算高时空分辨率蒸散的重要数据来源。基于遥感的蒸散估算包括统计模型与机理模型。

统计模型可以分为两类：(1)回归模型：回归模型是通过现有的通量塔蒸散数据与温度、降水等气象数据、NDVI等遥感指数数据建立回归关系，并将该回归关系应用于全区进行蒸散估算。回归模型结构简单，通过简单的栅格运算即可获得蒸散，但是回归模型的参数与结构是基于已有实测数据拟合的，在运用到新的地区或者新的气候带时，往往需要进行参数的重新校正。(2)机器学习模型：利用随机森林、人工神经网络、与支持向量机等方法，通过建立气象因素、地形因素与蒸散的非线性关系进行模型训练，并将训练好的模型外推到其他区域。但是机器学习模型忽略了内在的物理机理，很难说清楚蒸散与驱动因素之间的交互作用，同时也受限于数据的质、量和数据的典型性。Chen等(2014)在评估8个模型的效果时，其中有三个机器学习模型(人工神经网络、支持向量机和决策树)，模拟精度最差。

机理模型包括表面能量平衡模型(Surface energy balance algorithm forland，SEBAL)、Penman-Monteith(PM)模型与Priestley-Taylor(PT)模型等。地表能量平衡模型在计算过程中涉及到一些微气候参数只能通过经验获取，如动力学粗糙度、热力学粗糙度等参数，因此，此方法发展受到限制；在此基础上，Penman将地表能量平衡模型与空气动力学原理结起来提出湿润下垫面的蒸散发计算方法，Monteith引入了表面抗阻的概念，提出了Penman-Monteith模型。理论上来说，PM模型比PT模型机理更加清楚，但是当利用遥感进行全区域计算时，很难量化PM模型中的气孔导度与表面导度等参数。相对PM模型，PT模型需要的参数较少且能实现近似的效果，因此PT模型得到了更广泛的应用，特别是结构简单，机理清晰的Priestley-Taylor Jet Propulsion Lab(PT-JPL)模型(Fisher等，2008；Yao等，2013；Yao等，2015)。

包括PT-JPL在内的所有模型都存在一个相同的问题：即所有的模型都是基于陆表温度进行计算。但根据实际过程，真正决定植被蒸腾的应该是叶温(leaf temperature，T_l)而非气温(air temperature，T_a)(Leuning等，1995；Tuzet等，2003)，并且一般叶温与气温也并不相等(Blonder和Michaletz，2018；Michaletz等，2016；Michaletz等，2015)，此外，气温与叶温的关系也并不固定，叶温有时候会比气温高7-10℃(Korner等，2002；Gates等，1964；Lange，1959)，有时甚至高出20℃(Doughty和Goulden，2008)。叶温随气温的调节其实是植物自适应的结果，自适应将会有利于植物的生长、生存与生产(Blonder和Michaletz，2018)。因此，使用气温代替叶温的做法会错误的评估温度对蒸腾的影响，特别是在气温与叶温相差较大的地方。

除了叶温，土壤水分对蒸散有非常重要的影响(Detto等，2006；Krishnan等，2006；Wang和Liang，2008)。多项研究建议使用遥感反演的高空间异质性的土壤水分来估算土壤蒸发(Gokmen等，2012；Jin等，2011；Miralles等，2011；Yao等，2015)。但是之前基于卫星的土壤水分只考虑了约2-5cm的表层土壤水分，忽略了更深层的土壤水分对土壤蒸发的影响(Anderson等，2007；Yao等，2015)。近几年，新兴的不同土壤层的土壤水分数据为估算土壤蒸发提供了直接有效的数据源。水汽压偏差能够直接描述空气的干燥程度，同一水分状况下，温度越高，VPD越大。

发明内容

本发明要解决的技术问题在于提供一种精度更高的基于遥感的陆表蒸散计算方法。

本发明对机理模型PT-JPL中的温度子模型与土壤水分子模型进行优化，并将该模型运用于青藏高原潜在蒸散与实际蒸散的计算，将实际蒸散与青藏高原涡旋通量塔的观测数据(真实值)对比做精度验证。并与FLDAS、MOD16A2等模型精度与本方法进行对比，通过均方根误差、偏差与决定系数三个指标来验证本发明优化模型的效果。

本发明的技术方案为：

一种基于遥感的陆表蒸散计算方法，实际蒸散与植被表面蒸腾E_s、土壤表面蒸发E_s及拦截蒸发E_i有关，由下面模型计算：

LE＝L(E_c+E_s+E_i)

其中，L表示蒸发潜热(MJ/kg)，E为蒸散(mm)，蒸发潜热与蒸散的乘积LE表示潜热通量(W/m²)；E_c为干燥的植被表面蒸腾(mm)，E_s为干燥的土壤表面蒸发(mm)，E_i表示拦截蒸发(mm)；LE_c，LE_s与LE_i分别表示各部分的潜热通量；α_PT为无量纲系数，一般取值1.26；Δ为饱和水汽压对温度的曲线斜率(kPa/℃)；γ是干湿表常数(kPa/℃)；R_nc为植被表面的净辐射(W/m²)；R_ns为土壤表面净辐射(W/m²)；G为土壤热通量(W/m²)。

由以上公式，可得潜在蒸散(Potential Evapotranspiration,PE)与实际蒸散(Estimated Evapotranspiration,E)分别为：

PE在k_c、k_i与k_s下转为E，k_c、k_i与k_s分别表示植被潜在蒸腾到实际蒸腾、植被表面潜拦截蒸发到实际拦截蒸发、土壤拦截潜在蒸发到实际蒸发的限制因子，计算方法分别为：

k_c＝(1-f_wet)f_gf_Tf_M

k_s＝f_wet+(1-f_wet)f_sm

k_i＝f_wet

其中，f_wet表示拦截表面占比；f_g、f_T、f_M分别表示植物表面受到的绿度、温度以及水分的胁迫系数；f_sm为土壤表面的水分胁迫系数。以上五个指标的取值范围都是0-1，当数值为1时，表示不存在环境胁迫；数值越小，表示环境胁迫越大。

各因子的计算方法如下：

(1)温度胁迫系数(f_T)

改进后f_T的计算方式如下：

其中，T_l表示叶片温度；T_opt表示最适温度；f_T的数值范围为0-1。

(2)拦截表面占比(f_wet)

式中R_h表示相对湿度，计算方法如下：

其中，p为大气压，q为绝对湿度，T_a为空气温度。

(3)叶温(T_l)

其中，T_a为空气温度；R_nc为落在植被冠层的净辐射；c_P为空气定压比热；δ_b为叶片边界导度，

其中，μ代表风速，d与叶片宽度有关。

(4)植被绿度(f_g)

植被绿度表示植被蒸腾中，被绿色植被吸收的部分，植被绿度的计算方法为：

其中，f_g表示植被绿度；f_APAR表示活体植被吸收的陆表太阳辐射占比；f_IPAR表示由植被拦截的所有太阳辐射占比。各个因子的具体表达式为：

f_APAR＝1.36SAVI-0.048

f_IPAR＝NDVI-0.05

SAVI＝1.5×(ρ_NIR-ρ_Red)/(ρ_NIR+ρ_Red+0.5)

DNVI＝(ρ_NIR-ρ_Red)/(ρ_NIR+ρ_Red)

其中，SAVI为土壤校正的植被指数，ρ_NIR与ρ_Red分别为近红外波段和红色波段的反射率。

(5)植被水分胁迫系数(f_M)

植被水分限制的表达式为：

其中，

表示该像元绿色植被吸收的太阳辐射占比的最大值。

(6)土壤水分胁迫系数(f_sm)

本方法优化的土壤水分胁迫为：

其中，SM表示土壤水分含量，SM_max与SM_min分别表示土壤水分含量的最大值与最小值。

(7)干湿表常数(γ)

其中，C_p为空气定压比热，表示在一定气压下，单位体积的空气温度每升高1℃所需的能量；P为大气压，与海拔有关；λ为蒸发潜热；ε为水汽分子量与干空气分子量之比。

(8)植被表面(R_nc)和土壤表面(R_ns)所能吸收的太阳辐射分别为：

R_nc＝FVC×R_n

R_ns＝(1-FVC)×R_n-G

其中，FVC表示植被覆盖度，R_n为植被净辐射，G表示土壤热通量。

附图说明

图1，遥感模拟与气象站实测相对湿度的对比验证。

图2，潜热通量的验证，不同模型之间真实值与模拟值的对比：a)Fisher；b)MOD16A2；c)TerraClimate；d)FLDAS；e)PT-Hybrid；f)本文。

具体实施方式

本发明以青藏高原作为研究对象。

1研究区概况

青藏高原具有“世界屋脊”的美称，位于中国的西南部，经度26°03’N～40°00’N，纬度73°06’E～105°08’E，平均海拔4000米以上，覆盖面积2572.4×10³km²，是长江、黄河、澜沧江以及诸多河流的发源地，中国地表径流的37％都来自青藏高原。该区域重要的水源涵养功能使得蒸散的评估对于准确认识整个生态系统的水热循环十分重要。

2数据与数据预处理

2.1气象数据

气象数据来自中国气象数据网(http://www.nmic.cn/)的中国地面气候资料日值数据集(V3.0)，该数据集收集了全国近700个气象站点的蒸发、气温、相对湿度等8类信息。由于研究需要，本文将2000-2019青藏高原区域的65个气象站的相对湿度日值数据由MATLAB求取算术平均值得到月值数据。

2.2遥感数据

遥感数据包括MODIS原始数据与FLDAS、TerraClimate再分析数据集。MODIS是搭载于Terra卫星(10:30)与Aqua卫星(1:30)上的光学传感器，传感器上带有36个波段，波长范围自0.415微米到14.235微米不等，空间分辨率自250米至1000米不等(Running等，1994；Wang等，2011)。MODIS提供了NDVI(MOD13Q1)、陆表温度(MOD11A2)、叶面积指数(MOD15A2)等多种数据产品，已经成为估算蒸散(Chen等，2014；Hao等，2019；Yao等，2015)的重要数据源。FLDAS数据集(NcNally等，2017)旨在帮助数据少的发展中国家进行粮食安全评估。FLDAS使用Noah 3.6.1版地面模型，其CHIRPS-6每小时降水量已使用NASA陆地地面数据工具包进行了缩减，是NASA土地信息系统(NASA Land Information System，LIS；http://lis.gsfc.nasa.gov/)的一部分，提供了1982至今的蒸散、下行太阳辐射、土壤热通量、不同土壤层的水含量等多种月值数据产品。TerraClimate数据集(Abatzoglou等，2018)提供了全球陆表的最大空气温度、最小空气温度与10m处风速等。

数据产品及相关描述如表1所示。MOD13Q1的近红外波段NIR与红色波段Red可计算得到植被绿度所需的SAVI与NDVI，增强植被指数EVI可用来求解植被覆盖度。MOD15A2的光合有效辐射分量FPAR与FLDAS的净辐射可用来计算植被吸收的太阳辐射总量，进而可以计算植被的水分限制因子。SRTMDEM高程用来求解大气压，进而求解干湿表常数，同时也配合FLDAS含湿量、TerraClimate空气温度来求解相对湿度。TerraClimate的风速用来求解叶温。FLDAS的土壤水分数据用来求解土壤蒸发过程中的水分限制。MOD16A2的蒸散、TerraClimate的蒸发潜热与FLDAS的蒸散用来进行模型精度的对比验证。

为了对比FLDAS土壤水分含量的精度，本发明还收集了一套降尺度数据集(SoilMoisture in China dataset，后文简称CN_SM)与一套欧洲航天局的全球数据产品(ESACCI SM v06.1，后文简称ESA_SM)。数据源与数据单位如表所示。

表1数据源

3通量数据

使用地表观测的涡旋通量塔数据(Ground-measurements of eddy covariance,ECOR)验证蒸散算法的精确性。蒸散数据来自中国通量塔与相关文献。中国通量塔的数据下载自相关网站(www.chinaflux.org/)，数据包括海北灌丛、海北湿地与当雄三个观测站，每个观测站提供了半小时观测记录、日观测记录、月观测记录与年观测记录。唐古拉、玛曲与那曲的观测数据来自相关文献，具体信息与相关文献如表2所示。

表2中国通量塔与相关研究的ECOR地面观测

其中，ECOR提供的LE数据存在能量不平衡的问题，即太阳辐射与土壤热通量的差不等于潜热通量与显热通量的和。即：

R_n-G≠LE+H

因此，在使用LE之前需要进行一个波文比的校正从而将观测的LE转化为实际的LE(Twine等，2000)，校正方法如下：

其中，LE_r表示波文比校正后的潜热通量，LE和H分别表示观测的潜热通量和显热通量。R_n表示陆表净辐射，G表示土壤热通量。为了长时间序列研究的需要，最终将通量塔观测数据以及相关文献获取的数据全部转换为月尺度。

4模型结构

经过Priestley和Taylor在开放水域和干旱区的实验研究，确定PT模型的基本结构为：

L为蒸发潜热，取值2.45MJkg^-1；E为陆表蒸散；α_PT为无量纲系数，取值1.25-1.31，一般取1.26；Δ为饱和水汽压对温度的曲线斜率；γ是干湿表常数；R_n为植被表面的净辐射。

由于PT模型计算的是潜在蒸散，潜在蒸散转换为实际蒸散时，并未考虑到温度、水分与植被绿度等环境因子对LE的影响，因此Fisher等(2008)提出PT-JPL模型，该模型认为实际蒸散E与冠层表面蒸腾(E_c)、土壤表面蒸发(E_s)及拦截损失(E_i)有关，具体计算方法如下：

LE＝L(E_c+E_s+E_i)

其中，L表示蒸发潜热(MJ/kg)；E为蒸散(mm)；LE为蒸发潜热与蒸散的乘积，表示潜热通量(W/m²)；E_c为冠层表面蒸腾(mm)；E_s为土壤表面蒸发(mm)；E_i表示拦截蒸发(mm)；LE_c，LE_s与LE_i分别表示各部分的潜热通量；α_PT为无量纲系数，一般取值1.26；Δ_l，Δ_a与Δ_s分别表示饱和水汽压对叶片温度、气温及土壤温度的曲线斜率(kPa/℃)；γ是干湿表常数(kPa/℃)；R_nc为冠层表面的净辐射(W/m²)；R_nc为土壤表面净辐射(W/m²)；G为土壤热通量(W/m²)。

由以上公式，可得潜在蒸散(Potential Evapotranspiration，PE)与实际蒸散(Estimated Evapotranspiration，E)分别为：

PE在k_c、k_i与k_s下转为E，k_c、k_i与k_s分别表示植被潜在蒸腾到实际蒸腾、植被表面潜拦截蒸发到实际拦截蒸发、土壤拦截潜在蒸发到实际蒸发的限制因子，计算方法分别为：

k_c＝(1-f_wet)f_gf_Tf_M

k_s＝f_wet+(1-f_wet)f_sm

k_i＝f_wet

其中，f_wet表示拦截表面占比；f_gf_Tf_M分别表示植物表面受到的绿度、温度以及水分的限制系数；f_sm为土壤表面的水分限制系数。以上五个指标的取值范围都是0-1，当数值为1时，表示不存在环境胁迫；数值越小，表示环境胁迫越大。

以上参数所设涉及到的具体模型于下方分别计算：

(1)温度胁迫系数(f_T)

模型本身f_T的计算方法为：

其中，T_a为空气温度；T_opt表示最适温度。由于影响植被蒸腾的实际上是植被表面的温度而不是空气温度，因此，本发明使用植被表面的温度来替代空气温度来计算来计算温度胁迫系数。改进后f_T的计算方式如下：

其中，T_l表示叶片温度；T_opt表示最适温度；f_T的数值范围为0-1。

(2)拦截表面占比(f_wet)

拦截表面占比与空气的相对湿度有关：

其中，R_h表示相对湿度。相对湿度是计算水气压差的重要数据源，但目前许多遥感产品只提供了绝对湿度。因此，有必要将绝对湿度与相对湿度做一个转换。采用https://earthscience.stackexchange.com/questions/2360/how-do-i-convert-specific-humidity-to-relative-humidity的推导方程：

其中，w为实际水汽压，w_s为饱和水汽压，p为大气压，q为绝对湿度，T_a为空气温度。

(3)叶温(T_l)

其中，T_l表示叶片温度，T_opt表示最适温度。f_T的数值范围为0-1。

叶温影响植被蒸腾与植被光合，植被蒸腾又反过来调节叶温。叶温的变化与空气温度、陆表太阳辐射等因素非线性相关，本文方法采用Dong等(2017)进行叶温与气温的转换：

其中，T_a为空气温度；R_nc为落在植被冠层的净辐射；c_P为空气定压比热；δ_b为叶片边界导度，

其中，μ代表风速，d与叶片宽度有关。

(4)植被绿度(f_g)

植被绿度表示植被蒸腾中，被绿色植被吸收的部分。植被绿度的计算方法为：

其中，f_g表示植被绿度，f_APAR表示绿色植被吸收的陆表太阳辐射占比，f_IPAR表示由植被拦截的所有太阳辐射占比。各个因子的具体表达式为：

f_APAR＝1.36SAVI-0.048

f_IPAR＝NDVI-0.05

SAVI＝1.5×(ρ_NIR-ρ_Red)/(ρ_NIR+ρ_Red+0.5)

NDVI＝(ρ_NIR-ρ_Red)/(ρ_NIR+ρ_Red)

其中，SAVI为土壤校正的植被指数，ρ_NIR与ρ_Red分别为近红外波段和红色波段的反射率。

(5)植被水分胁迫系数(f_M)

植被水分限制的表达式为：

其中，

表示该像元绿色植被吸收的太阳辐射占比的最大值。

(6)土壤水分胁迫系数(f_sm)

原方法中f_sm定义土壤水的胁迫系数为：

f_sm＝R_h ^VPD

其中，R_h表示相对湿度，VPD为水气压偏差。

R_h和VPD都能反应空气的干燥程度，可以近似替代但不能直接反应土壤的水分胁迫状况。为了提高评估结果的可靠性，我们直接使用近年来发展的遥感土壤水数据量化土壤水分胁迫。土壤水、VPD与土壤水分胁迫(水分胁迫系数)存在以下关系：土壤水含量越高，水分胁迫越小(水分胁迫系数越大)；VPD越高，对土壤水分的需求越大，土壤水分胁迫越小(水分胁迫系数也越大)。因此，土壤水分胁迫与VPD都与土壤水分胁迫系数成正比。在两者的协同作用下，本方法定义土壤水分胁迫方程为：

其中，SM表示土壤水分含量，SM_max与SM_min分别表示土壤水分含量的最大值与最小值。

由于研究对象是草地，根深较浅，因此土壤水含量的取值由0-10cm与10-40cm获取。0-40cm土壤层水含量按加权平均表示为：

SM＝0.25×SM1+0.75×SM2

其中，SM1为0-10厘米的表层土壤水，SM2为10-40厘米的土壤水。

SM_max与SM_min分别表示土壤水分含量的最大值与最小值。本文计算20年土壤水分的最大值作为SM_max，土壤水分最小值作为SM_min。

(7)干湿表常数(γ)

干湿表常数的表达式为：

其中，C_p为空气定压比热，表示在一定气压下，单位体积的空气温度每升高1℃所需的能量；P为大气压，与海拔有关；λ为蒸发潜热；ε为水汽分子量与干空气分子量之比。

(8)植被表面(R_nc)和土壤表面(R_ns)所能吸收的太阳辐射分别为

R_nc＝FVC×R_n

R_ns＝(1-FVC)×R_n-G

其中，FVC表示植被覆盖度，R_n为植被净辐射，G表示土壤热通量。

基于EVI计算植被覆盖度的表达式为：

G＝a_g(1-FVC)R_n

其中，EVI表示增强植被指数，EVI_max与EVI_min分别表示该像元增强植被指数的最大值与最小值；K、L、C₁、和C₂分别是2.5、1、6.0和7.5；a_g为常数0.18，FVC为植被覆盖度，R_n为净辐射。

5数据分析与评价方法

本发明以均方根误差、偏差及R²三种统计分析方法对蒸散的模拟效果进行评价。

均方根误差(RootMean Square Error，RMSE)也称为标准误差，能够解释模拟值多么集中的围绕着观测值。它是观测值与模拟值偏差的平方和与观测值个数比值的平方根。

其中，RMSE表示均方根误差，y_i为第i个观测值，

为第i个模拟值，N为观测值个数。

偏差(Bias)描述了观测值与模拟值的偏差程度，反映了模型本身的拟合能力。

其中，Bias表示偏差，

表示模拟值，y_i表示观测值，N为观测值个数。偏差的单位与观测值一致。

R²也称为拟合优度，表征观测值中可有模拟值解释的百分比。

其中，R²范围为0-1，数值越大，表示可解释的百分比越高。

表示模拟值，y_i表示观测值，N为观测值个数，

表示所有观测值的平均值。

6结果验证与分析

6.1相对湿度验证

本节主要是对相对湿度的推算方程的精度进行验证。验证数据使用国家气象数据中心相对湿度的站点数据。数据包括2000-2019年14400条数据记录。对比分析结果如图1所示，均方根误差RMSE为8.55％，决定系数R²为0.401，偏差为2.80％。由此可以证明，由绝对湿度到相对湿度的转换是可靠的。

6.2土壤水分验证

结果如表3所示，我们计算了FLDAS数据集(FLDAS_SM)、SMC数据集(SMC_SM)和ESACCI SM v06.1(ESA_SM)数据集中土壤水分与观测数据的相关性，并通过加权平均法得到不同模型的平均相关性的值。在对比前，我们将日值ESA_SM转换为月值，但仍存在许多数据空缺。因此，为了结果的客观性，我们只使用它们的交集进行对比验证。每个站点的空缺值不同，每个站点的采样数(n)与其他站点也不完全相同。另外需要注意的是，在清水河站点，它只有一个有效记录，因此这里不做准确性验证。三者的比较结果中，FLDAS_SM数据集在贵德站的精度最高，ESA_SM在玉树的精度最差。综合来看，FLDAS_SM、SMC_SM、ESA_SM对应的平均系数分别为0.46、0.34、0.40。

表3三种土壤水分产品与观测数据的相关性对比

6.3不同蒸散模型的对比

为了验证本发明模型的效果，MOD16(Mu等，2007；Mu等，2011)，PT-JPL(Fisher等，2008)，PT-Hybrid(Yao等，2015)与FLDAS(McNally等，2017)都被用来作为对比。本节对模型的RMSE,Bias和R²做了统计分析。

各个遥感模型与通量塔观测潜热通量的统计结果如表4与图2所示。表4中对比了唐古拉、玛曲、那曲、海北灌丛、海北湿地与当雄六个站点的偏差与均方根误差。图2对比了六个站点的R²。

从表4来看，当雄的模拟效果最好，偏差为2.69W/m²，均方根误差为10.15W/m²；其次分别是海北湿地与海北灌丛；玛曲、那曲与唐古拉的模拟效果较差，因为这三个站点的数据来自相关文献，没有基于站点的显热通量、土壤热通量等数据，因此没有进行波文比校正。从模型之间的对比来看，本文的模拟结果最好，FLDAS的效果较好，MOD16与PT-Hybrid的效果最差。由于FLDAS本身就是基于非洲农业发展起来的(McNally等，2017)，因此适用于干旱区与半干旱区。MOD16是基于PM的全球通用模型，涉及气孔导度等难以量化的参数，最终表现在局地的模拟效果较差。PT-Hybrid在更换气象数据时，模型的参数往往需要进行重新校正，因此在本文中的数据进行模拟时，其效果较差。

从图2可知，PT-JPL模型的总体R²为0.727，MOD16的总体R²为0.547，TerraClimate的总体R²为0.878，PT-Hybrid的总体R²为0.118，本文的总体R²为0.903。从总体的分布趋势来看，PT-JPL模型，MOD16与PT-Hybrid都有较为明显的低估趋势。FLDAS的模拟值与观测值趋势较为一致，模拟效果与本文模拟结果最为接近，只有极个别值FLDAS比本文更好。但是由于FLDAS更适合于干旱与半干旱区，因此，在海北湿地站，该模型会有一定的低估。本文直接量化了土壤水含量，在海北湿地依然具有较好的模拟效果。

表4 PT-JPL、MOD(MOD16A2)、TC(TerraClimate)、FLDAS、PT^H(PT-Hybrid)与本发明的均方根误差、决定系数的对比

本发明优化了一种蒸散的计算方法，将植被蒸腾中的气温替换成叶温，将土壤蒸发中的水分限制因子替换为归一化的土壤水分。本发明具有以下优点：

1、数据易获取，蒸散的计算全部依赖于免费的遥感数据、气象数据与通量塔数据，这些数据均可通过相关网站免费下载获得。

2、计算结果更加可靠。本发明使用不同土壤层的土壤含水量来量化土壤水分对土壤蒸发的影响；利用机理模型计算叶温与气温的非线性关系，并利用叶温来计算植被蒸腾中的温度限制。这种计算方式的理论基础和结果精度都更加可靠。

3、本发明采用遥感的方法计算蒸散，能够实时的反映不同时空尺度的潜在蒸散与实际蒸散现状，能够综合反映一个地区的水分状况，有利于准确地量化水循环过程。

Claims (1)

1.一种基于遥感的陆表蒸散计算方法，其特征在于，实际蒸散(E)与冠层表面蒸腾(E_c)、土壤表面蒸发(E_s)及拦截蒸发(E_i)有关，具体计算方法如下：

LE＝L(E_c+E_s+E_i)

其中，L表示蒸发潜热(MJ/kg)；E为蒸散(mm)；LE为蒸发潜热与蒸散的乘积，表示潜热通量(W/m²)；E_c为冠层表面蒸腾(mm)；E_s为土壤表面蒸发(mm)；E_i表示拦截蒸发(mm)；LE_c，LE_s与LE_i分别表示各部分的潜热通量；α_PT为无量纲系数，一般取值1.26；Δ_l，Δ_a与Δ_s分别表示饱和水汽压对叶片温度、气温及土壤温度的曲线斜率(kPa/℃)；γ是干湿表常数(kPa/℃)；R_nc为冠层表面的净辐射(W/m²)；R_ns为土壤表面净辐射(W/m²)；G为土壤热通量(W/m²)；

由以上公式，可得潜在蒸散(PE)与实际蒸散(E)分别为：

PE在k_c、k_i与k_s下转为E，k_c、k_i与k_s分别表示植被潜在蒸腾到实际蒸腾、植被表面潜拦截蒸发到实际拦截蒸发、土壤拦截潜在蒸发到实际蒸发的限制因子，计算方法分别为：

k_c＝(1-f_wet)f_gf_Tf_M

k_s＝f_wet+(1-f_wet)f_sm

k_i＝f_wet

其中，f_wet表示拦截表面占比；f_gf_Tf_M分别表示植物表面受到的绿度、温度以及水分的胁迫系数；f_sm为土壤表面的水分胁迫系数；以上五个指标的取值范围都是0-1，当数值为1时，表示不存在环境胁迫；数值越小，表示环境胁迫越大；

各因子的计算和改进方法如下：

(1)温度胁迫系数(f_T)

其中，T_l表示叶片温度；T_opt表示最适温度；f_T的数值范围为0-1。

(2)拦截表面占比(f_wet)

其中，R_h表示相对湿度；

(3)叶温(T_l)

叶温(T_l)表达式为：

其中，T_a为空气温度；R_nc为落在植被冠层的净辐射；c_P为空气定压比热；δ_b为叶片边界导度，

其中，μ代表风速，d与叶片宽度有关；

(4)植被绿度(f_g)

植被绿度表示植被蒸腾中，被绿色植被吸收的部分，植被绿度的计算方法为：

其中，f_g表示植被绿度；f_APAR表示活体植被吸收的陆表太阳辐射占比；f_IPAR表示由植被拦截的所有太阳辐射占比；

(5)植被水分胁迫系数(f_M)

植被水分限制的表达式为：

其中，

表示该像元绿色植被吸收的太阳辐射占比的最大值；

(6)土壤水分胁迫系数(f_sm)

其中，SM表示土壤水分含量，SM_max与SM_min分别表示土壤水分含量的最大值与最小值；VPD表示水汽压偏差；

(7)干湿表常数(γ)

其中，C_p为空气定压比热，表示在一定气压下，单位体积的空气温度每升高1℃所需的能量；P为大气压，与海拔有关；λ为蒸发潜热；ε为水汽分子量与干空气分子量之比；

(8)太阳辐射：植被表面(R_nc)和土壤表面(R_ns)

R_nc＝FVC×R_n

R_ns＝(1-FVC)×R_n-G

其中，FVC表示植被覆盖度，R_n为植被净辐射，G表示土壤热通量。

Images