CN111553459B - 一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法 - Google Patents

Abstract

本发明公开了一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法，所述估算方法通过将基于过程的蒸散发模型和基于温度的蒸散发模型进行集成来实现。本发明公开的基于遥感信息的灌溉地区逐日实际蒸散发估算方法，操作便捷，易于实现，能够对灌溉地区进行逐日、高精度的实际蒸散发估算，相较于S‑W模型，对逐日实际蒸散发的估算精度能够提高约50％。

Description

一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法

技术领域

本发明属于农业、遥感与水文测量技术领域，具体涉及一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法及系统。

背景技术

蒸散发(Evapotranspiration，ET)是地表水量平衡的重要分量，实际蒸散发在干旱和半干旱地区的灌溉系统设计和水资源管理中起着重要作用。由于干旱和半干旱地区降水稀少，灌溉成为作物生长的主要水量来源。例如，美国灌溉农业提供全国约50％的粮食供给，但农业灌溉用水量却占总淡水消耗量的60％。另外，由于地下水位下降和过渡灌溉导致的土地盐碱化问题，一直威胁着全球的粮食安全。因此，亟需对灌溉地区的水资源进行有效管理，而这些都建立在准确的估算逐日实际蒸散发的基础之上。

陆表蒸散发的估算通常有三种方法：站点观测、模型模拟和遥感方法。其中，利用蒸渗仪等进行地面站点观测是较为传统的方式，但由于成本昂贵，这种方法并不适合大尺度应用；依赖水文模型或陆面过程模型的模型模拟方法，由于缺少可靠的灌溉信息，水量平衡不能真实计算，往往会严重低估灌区的实际蒸散发。因此，基于遥感信息的蒸散发模型得到越来越广泛的关注。

遥感蒸散发模型包括两大类，一类是基于过程的模型，通常由大量的遥感可获取地表参量驱动，其可以估算连续的、日尺度实际蒸散发，但由于缺少可靠的地表表层土壤湿度和根区土壤湿度数据，以及较大的气象数据误差，容易造成对灌区蒸散发的严重低估；另一类是基于地表温度的遥感蒸散发模型，其在晴空日可以获取较准确的、非连续的日蒸散发，但受云影响较大，不能获取时间连续的逐日蒸散发。

因此，亟需提供一种能够对灌溉地区的实际蒸散发实现逐日、高精度估算的方法。

发明内容

为了克服上述问题，本发明人进行了锐意研究，设计出一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法，其首先利用基于温度的蒸散发模型获取晴空日、非逐日、但精度较高的日蒸散发；然后利用优化算法和获取的晴空日、高精度日蒸散发对基于过程的蒸散发模型进行优化；最后，利用优化后的模型高精度估算逐日实际蒸散发，从而完成了本发明。

具体来说，本发明的目的在于提供一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法，其中，所述估算方法通过将基于过程的蒸散发模型和基于温度的蒸散发模型进行集成来实现。

其中，所述基于过程的蒸散发模型为S-W双源模型；

所述基于温度的蒸散发模型为温度-植被指数模型。

其中，所述基于过程的蒸散发模型和基于温度的蒸散发模型通过粒子群优化算法集成。

其中，所述估算方法包括以下步骤：

步骤1，采用基于温度的蒸散发模型，获取灌溉地区晴空日的日蒸散发；

步骤2，在晴空日，对基于过程的蒸散发模型进行优化；

步骤3，估算灌溉地区逐日实际蒸散发。

其中，步骤1包括以下子步骤：

步骤1-1，选择晴空日，获得灌溉地区的蒸发比；

步骤1-2，根据获得的蒸发比，获取晴空日的日蒸散发。

其中，步骤2中，所述基于过程的蒸散发模型为S-W双源模型，

所述优化为对相关参数进行校正优化，相关参数为土壤阻抗和植被冠层阻抗。

其中，步骤2中，对土壤阻抗和植被冠层阻抗的优化校正包括以下步骤：

步骤2-1，设定原始估算的土壤阻抗与优化后土壤阻抗的关系式，并设定原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系式；

步骤2-2，获得步骤2-1所述关系式中的土壤阻抗变化系数和植被冠层阻抗变化系数。

其中，步骤2-1中，原始估算的土壤阻抗与优化后土壤阻抗的关系如下式(十六)所示：

原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系如下式(十七)所示：

其中，

为原始估算的土壤阻抗，

为优化后土壤阻抗，

为原始植被冠层阻抗，

为优化后植被冠层阻抗，a_s和b_s为土壤阻抗变化系数，a_c和b_c为植被冠层阻抗变化系数。

其中，在步骤2-2之后，还包括对土壤阻抗变化系数和植被冠层阻抗变化系数进行优化的步骤。

其中，步骤3中，所述灌溉地区逐日实际蒸散发通过基于过程的蒸散发模型估算获得；

其中，所述基于过程的蒸散发模型为优化后的模型。

本发明所具有的有益效果包括：

(1)本发明提供的基于遥感信息的灌溉地区逐日实际蒸散发估算方法，操作便捷，易于实现；

(2)本发明提供的基于遥感信息的灌溉地区逐日实际蒸散发估算方法，将基于过程的蒸散发模型和基于地表温度的蒸散发模型进行有效集成，能够对灌溉地区进行逐日、高精度的实际蒸散发估算；

(3)本发明提供的基于遥感信息的灌溉地区逐日实际蒸散发估算方法，相较于S-W模型，对逐日实际蒸散发的估算精度能够提高约50％；

(4)本发明提供的基于遥感信息的灌溉地区逐日实际蒸散发估算方法，估算效率和稳定性均较高，适合规模化使用。

附图说明

图1示出本发明实施例1中的黑河研究区的高程模型及研究区、站点分布图；

图2示出本发明实施例1、对比例1估算的日蒸散发与站点观测的时间序列比较图(2009～2011年)；

图3示出本发明实施例1、对比例2估算的日蒸散发与站点观测的时间序列比较图(2009～2011年)；

图4示出本发明实施例1、对比例2估算的日蒸散发与站点观测的数据比对图(2012年)；

图5示出本发明实施例1、对比例1和对比例2在黑河中游灌区的估算结果空间对比图。

具体实施方式

下面通过优选实施方式和实施例对本发明进一步详细说明。通过这些说明，本发明的特点和优点将变得更为清楚明确。

在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。

本发明人研究发现，基于过程的蒸散发模型和基于地表温度的蒸散发模型在连续性和精度方面具有非常好的互补性，将二者集成，将可以实现对逐日实际蒸散发的估算。

因此，本发明提供了一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法，所述方法通过将基于过程的蒸散发模型和基于温度的蒸散发模型进行集成来实现。

根据本发明一种优选的实施方式，所述基于过程的蒸散发模型为S-W双源模型。

其中，基于过程的蒸散发模型以相似理论为核心，以地-气间水热传输的物理过程为基础，在综合考虑土壤供水能力、植被吸水能力和大气需水能力的情况下，以太阳短波辐射和风能为主要能量来源，进行蒸散发估算完整推演的模型。该模型基于对地表表层和植被根区土壤水分经过土壤蒸发或植被蒸腾过程的精确表达，具有较强的逻辑性和过程严谨性。

基于过程的蒸散发模型，以Penman-Monteith模型和Shuttleworth and Wallace(S-W)双源模型为代表，其中，S-W双源模型，是为克服大叶模型在植被不完全覆盖情况下对地表水热交换过程不能有效刻画的问题而构建的，其在充分考虑土壤-植被系统内部水热交换过程的基础上，以串联方式组织蒸散发阻抗表达。

在本发明中，选择S-W双源模型的原因是该模型机理性强，对土壤-植被间的水热交互过程刻画细致明确。

在进一步优选的实施方式中，所述基于温度的蒸散发模型为温度-植被指数模型。

其中，基于温度的蒸散发模型以地表温度对地表蒸散发能力差异性表达为基础，通过地表温度来计算蒸散发。该模型主要包括两大类，一种是以地表能量平衡为基础，首先计算地表的热量交换量(感热通量)，然后通过能量平衡计算潜热通量，进而转换成蒸散发；另一种是以温度-植被指数为代表，通过遥感获取的地表温度和植被指数计算蒸发比，直接计算蒸散发。其中，由于蒸发吸热作用，地表温度在蒸散发估算中起着重要作用。

在本申请中，优选选择温度-植被指数模型，主要是由于该模型对温度反演误差不敏感，估算精度高，需要参数少。

根据本发明一种优选的实施方式，所述基于过程的蒸散发模型和基于温度的蒸散发模型通过粒子群优化算法集成。

其中，粒子群算法通过用大量粒子来模拟鸟群或鱼群的个体运动，通过空间寻优找到最佳目的地。每个粒子都具有速度和位置两个属性。其中，速度表示粒子移动的快慢，位置表示粒子移动的方向。每个粒子在搜索空间中独立运动，并搜寻个体极值，通过与群体中的其它粒子共享个体极值，找到群体极值，并在当前个体极值和群体极值的双重约束下来调整自己的位置和速度，并最终找到全局最优解。

本发明中优选采用粒子群算法对两个模型进行优化，有利于借助算法中粒子的记忆性传递信息，加快优化速度，提升优化效果。

根据本发明一种优选的实施方式，所述基于遥感信息的灌溉地区逐日实际蒸散发估算方法包括以下步骤：

步骤1，采用基于温度的蒸散发模型，获取灌溉地区晴空日的日蒸散发；

步骤2，在晴空日，对基于过程的蒸散发模型进行优化；

步骤3，估算逐日实际蒸散发。

以下进一步描述所述基于遥感信息的灌溉地区逐日实际蒸散发估算方法。

步骤1，采用基于温度的蒸散发模型，获取灌溉地区晴空日的日蒸散发。

其中，所述步骤1包括以下子步骤：

步骤1-1，选择晴空日，获得灌溉地区的蒸发比。

根据本发明一种优选的实施方式，在晴空日，采用温度-植被指数模型计算蒸散发。

具体地，所述温度-植被指数模型是指通过遥感获知的地表温度和植被指数计算蒸发比，直接计算蒸散发。采用该模型计算得到的蒸散发的估算精度较高。

其中，在干旱地区、地表相对平坦情况下，有足够多的干湿像元分布的地区，地表温度和植被指数可以构成典型的三角特征空间，其中，上边界称为干边，蒸散发为0，下边界称为湿边，蒸散发为潜在蒸散发。通过三角特征空间内任一点到湿边的距离除以干边到湿边的距离即为潜在蒸发的比例。

根据本发明一种优选的实施方式，在晴空日，利用遥感观测数据，采用下式(一)获得灌溉地区的瞬时蒸发比：

其中，EF为瞬时蒸发比，α为考虑空气动力学阻抗作用的综合参数，其取值范围为0～1.26；Δ为温度—饱和水汽压斜率(kPaK^-1)，Δ＝4098·(0.6108·exp(17.27·T/(T+237.3)))/(T+237.3)²，T为地表温度(℃)，来自于遥感数据；γ为干湿表常数(kPaK^-1)，由干湿球温度计测得。

优选地，上式(一)中的α按照文献“Tang,Ronglin,LiZhaoliang,and TangBohui.2010."An Application of the T-S-Vi Triangle Method with Enhanced EdgesDetermination for Evapotranspiration Estimation from Modis Data in and andSemi-Arid Regions:Implementation and Validation.".Remote Sensing OfEnvironment 114,no.3:540-51.”中所述方法获得。

在本发明中，对遥感数据的获取方法不做特别限定，可以采用现有技术中常见的方法获得。

本发明中，优选假定全天内的蒸散发不发生改变，将获得的瞬时蒸散发作为日尺度的蒸散发。

步骤1-2，根据获得的蒸发比，获取晴空日的日蒸散发。

根据本发明一种优选的实施方式，根据上述获得的蒸发比，结合日净辐射，采用下式(二)获得晴空日的日蒸散发：

其中，ET为日蒸散发，EF为瞬时蒸发比，R_n为日净辐射，λ为蒸发汽化潜热(Jkg^-1)，λ＝2501-2.361·T₀，T₀为气温(℃)，由再分析气象观测资料得到。

在本发明中，优选地，所述R_n采用下式(三)获得：

其中，

为下行短波辐射，由再分析气象观测资料得到(W m^-2)；

为下行长波辐射，由再分析气象观测资料得到(W m^-2)；α为反射率，来自于遥感数据；T_s为日均气温，由再分析气象观测资料得到(K)；δ为玻尔兹曼常数(5.67×10^-8，W m^-2 K-4)；ε为发射率，ε＝(1-FVC)·ε_s+FVC·ε_v，ε_s为土壤发射率，ε_v为植被发射率，FVC为植被覆盖度。

根据本发明一种优选的实施方式，所述晴空日为全年的晴空日。

步骤2，在晴空日，对基于过程的蒸散发模型进行优化。

根据本发明一种优选的实施方式，所述基于过程的蒸散发模型为S-W双源模型，其通过下式(四)获得总的蒸散发：

λ·ET＝C_c·PM_c+C_s·PM_s (四)

其中，λ为蒸发汽化潜热(Jkg-1)，ET为日蒸散发，PM_c、PM_s、C_c和C_s分别由下式(五)～式(八)获得：

C_C＝{1+R_CR_a/[R_S(R_C+R_a)]}^-1 (七)；

C_S＝{1+R_SR_a/[R_C(R_S+R_a)]}^-1 (八)；

其中，R_a、R_s、R_C、R_n,c、R_n,sR分别由下式(九)～式(十三)获得：

R_n,c＝FVC·R_n(十二)；

R_n,s＝(1-FVC)·R_n (十三)；

其中，R_n是日净辐射(Wm^-2)；R_n,s是土壤的净辐射(Wm^-2)；R_n,_c是植被部的净辐射(Wm^-2)；G是土壤热通量(Wm^-2)；FVC为植被覆盖度，FVC＝(NDVI-NDVI_min)/(NDVI_max-NDVI_min),NDVI是归一化植被指数，来自于遥感数据，NDVI_max和NDVI_min分别是图像覆盖区NDVI的最大值和最小值；

是冠层源区到参考高度处的空气动力学阻抗(sm^-1)；

是植被的整体边界层阻抗(sm^-1)；

是土壤到植被源区域间的空气动力学阻抗(sm^-1)；

是植被冠层阻抗(sm^-1)；

是土壤阻抗(sm^-1)；D是参考高度处饱和水汽压差(kpa)，由再分析气象观测资料得到；C_p是空气比热(Jkg^-1K^-1)；ρ是空气密度，(kgm^-3)；Δ是饱和水汽压斜率(kPaK^-1)；γ是干湿表常数(kPaK^-1)。

在本发明中，优选地，上述

和

可以按照文献“Shuttleworth,W.James,&Wallace,J.S.1985."Evaporation from Sparse Crops-an Energy CombinationTheory."[In English].Quarterly Journal Of the Royal MeteorologicalSociety111,no.469(1985):839-55.”中所述方法获得。

在进一步优选的实施方式中，所述S-W双源模型的优化为对相关参数进行校正优化，相关参数为土壤阻抗和植被冠层阻抗，其中，所述土壤阻抗

和植被冠层阻抗

分别由下式(十四)和式(十五)获得：

其中，θ为表层土壤湿度(cm³·cm^-3)，来自于遥感数据；θ_sat土壤饱和含水量(cm³·cm^-3)，来自于土壤质地数据。

其中，

LAI_eff＝LAI/(0.3LAI+1.2)；

f(T₀)＝1-((T₀-T_opt)/T_opt)²；

f(R_s,d)＝1-exp(-R_s,d/R_sd,opt)；

其中，r_min为气孔最小阻抗；LAI_eff为有效叶面积指数；LAI为叶面积指数，来自于遥感观测数据；f(T₀)为温度胁迫因子；T₀为气温(℃)，由再分析气象观测资料获得；T_opt为最优气温(℃)；R_s,d为短波辐射(Wm^-2)，由再分析气象观测资料获得；f(R_s,d)为短波辐射胁迫因子；R_sd,opt为最优短波(Wm^-2)；f(D₀)为水汽压胁迫因子；D₀为饱和水汽压差，由再分析气象观测资料获得；

为系数；f(θ₂)为土壤湿度胁迫因子；θ₂为根区土壤湿度(cm³·cm^-3)，来自于现有的陆面模式模拟产品；θ_wilt为植被凋萎土壤含水量(cm³·cm^-3)，来自土壤数据资料；θ_cr为田间持水量(cm³·cm^-3)，来自土壤数据资料。

其中，S-W双源模型是为了考虑土壤对稀疏植被覆盖下地表通量的贡献而提出的，该模型将底层土壤与植被冠层看作是上下叠加、彼此连续的湍流源，底层的水汽和热量汇合于植被冠层内部，并且只能通过顶层离开，整个冠层总热通量是各层热通量之和。其是在充分考虑土壤-植被系统内部水热交换过程的基础上，以串联方式组织蒸散发阻抗表达，而构建的模型。

该模型可以对蒸散发进行连续估算，对输入数据要求较高，土壤阻抗和植被冠层阻抗是估算的关键因素，但由于各种经验参数和气象遥感数据的不准确性，导致该模型估算的蒸散发会产生较大的误差，因此，本发明优选根据步骤1中由温度-植被指数模型获得的晴空日、高精度日蒸散发，对基于过程的蒸散发模型中土壤阻抗和植被冠层阻抗进行校正。

其中，对土壤阻抗和植被冠层阻抗的校正包括以下步骤：

步骤2-1，设定原始估算的土壤阻抗与优化后土壤阻抗的关系式，以及原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系式。

本发明人经过研究发现，S-W双源模型中通过式(十四)和式(十五)估算得到的原始土壤阻抗和原始植被冠层阻抗与优化后的土壤阻抗和植被冠层阻抗存在一定的关系，经过大量试验研究，综合估算的准确度、后续优化结果的稳定性及运算效率，得到原始估算的土壤阻抗与优化后土壤阻抗的关系如下式(十六)所示，原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系如下式(十七)所示：

其中，

为原始估算的土壤阻抗，

为优化后土壤阻抗，

为原始植被冠层阻抗，

为优化后植被冠层阻抗，a_s和b_s为土壤阻抗变化系数，a_c和b_c为植被冠层阻抗变化系数。

步骤2-2，获得步骤2-1所述关系式中的土壤阻抗变化系数和植被冠层阻抗变化系数。

本发明人研究发现，在晴空日，同时采用基于温度的“温度-植被指数”蒸散发模型和基于过程的S-W双源模型同时估算灌溉地区的日蒸散发，采用“温度-植被指数”蒸散发模型获得的高精度蒸散发对S-W模型中的原始土壤阻抗和植被冠层阻抗进行优化校正。

根据本发明一种优选的实施方式，在晴空日，采用粒子群优化算法获取优化后的土壤阻抗变化系数和植被冠层阻抗变化系数。

其中，采用粒子群优化算法有利于寻找全局最优解，获得最优的土壤阻抗变化系数和植被冠层阻抗变化系数。且该算法受原始估算过程中土壤阻抗和植被冠层阻抗误差的影响很小，有利于获得准确的阻抗变化系数。

根据本发明一种优选的实施方式，在步骤2-2之后，还包括以下步骤：

步骤2-3，重复步骤1、步骤2-1、步骤2-2四～六次，优选重复五次，以获得最优的土壤阻抗变化系数和植被冠层阻抗变化系数。

本发明人发现，在获取上述最优的阻抗变化系数的过程中，由于粒子群优化算法有一定概率落入局部最小解(本发明中小于10％)，因此需要重复步骤1、步骤2-1、步骤2-2。

本发明中经过大量研究发现，重复上述操作五次，能够获得最优的土壤阻抗变化系数和植被冠层阻抗变化系数。

步骤3，估算逐日实际蒸散发。

根据步骤2获得的最优土壤阻抗变化系数和植被冠层阻抗变化系数，结合式(十六)和式(十七)，能够得到优化后的土壤阻抗和植被冠层阻抗，在本发明中，该优化后的土壤阻抗和植被冠层阻抗可以作为真实土壤阻抗和植被冠层阻抗。

优选地，根据S-W双源模型，可以获得逐日、真实的土壤阻抗和植被冠层阻抗。

更优选地，根据获得的逐日、真实的土壤阻抗和植被冠层阻抗，结合遥感和气象资料，采用S-W模型的式(四)估算得到逐日、真实的蒸散发。

采用本发明中所述的基于遥感信息的灌溉地区逐日实际蒸散发估算方法，能够对日蒸散发进行逐日估算，且与S-W双源模型相比，能够将估算精度提高50％左右。

实施例

以下通过具体实例进一步描述本发明，不过这些实例仅仅是范例性的，并不对本发明的保护范围构成任何限制。

实施例1

本实施例以我国的黑河地区为研究区，其由三个较大的灌溉绿洲组成，如图1所示。其中，图1中的(a)为黑河流域数字高程模型图，图1中的(b)为研究区的土地覆盖图，图1中的(c)为研究区的验证站点分布图，(c)中的1～6及YK代表7个采用涡动观测仪进行数据观测的相关观测站。

在本实施例中，所述遥感数据来自MODIS的https://modis.gsfc.nasa.gov/。

在本实施例中，所述再分析气象观测资料来自国家青藏高原科学数据中心的https://data.tpdc.ac.cn/zh-hans/data/8028b944-daaa-4511-8769-965612652c49/。

按照下述步骤分别对上述研究区全年的逐日蒸散发进行估算：

步骤1，选择每个年度全年的晴空日，利用遥感观测的地表温度，使用温度-植被指数模型，采用下式(一)获得研究区的瞬时蒸发比：

其中，EF为瞬时蒸发比，α为考虑空气动力学阻抗作用的综合参数，按照文献“Tang,Ronglin,LiZhaoliang,and Tang Bohui.2010."An Application of the T-S-ViTriangle Method with Enhanced Edges Determination for EvapotranspirationEstimation from Modis Data in and and Semi-Arid Regions:Implementation andValidation.".Remote Sensing Of Environment 114,no.3:540-51.”中所述方法获得；Δ为温度—饱和水汽压斜率(kPaK^-1)，Δ＝4098·(0.6108·exp(17.27·T/(T+237.3)))/(T+237.3)²，T为地表温度(℃)，来自于遥感数据；γ为干湿表常数(kPaK^-1)，由干湿球温度计测得。

将瞬时蒸发比作为整日的蒸发比，结合气象数据计算的日净辐射，采用下式(二)获得晴空日的日蒸散发：

其中，ET为日蒸散发，EF为瞬时蒸发比，R_n为日净辐射，λ为蒸发汽化潜热(Jkg^-1)，λ＝2501-2.361·T₀，T₀为气温(℃)，由再分析气象观测资料得到。

R_n采用下式(三)获得：

其中，

下行短波辐射，由再分析气象观测资料得到(W m^-2)；

下行长波辐射，由再分析气象观测资料得到(W m^-2)；α为反射率，来自于遥感数据；T_s为日均气温，由再分析气象观测资料得到(K)；δ玻尔兹曼常数(5.67×10^-8W m^-2K^-4)；ε＝(1-FVC)·ε_s+FVC·ε_v，ε_s为土壤发射率，取值为0.96；ε_v为植被发射率，取值为0.98，FVC为植被覆盖度，FVC＝(NDVI-NDVI_min)/(NDVI_max-NDVI_min)，NDVI是归一化植被指数，来自于遥感数据，NDVI_max和NDVI_min分别是图像覆盖区NDVI的最大值和最小值，分别取0.9和0.1。

步骤2，

在晴空日，基于S-W双源模型，通过下式(四)获得总的蒸散发：

λ·ET＝C_c·PM_c+C_s·PM_s (四)

其中，λ为蒸发汽化潜热(Jkg^-1)，ET为日蒸散发。

PM_c、PM_s、C_c和C_s分别由下式(五)～式(八)获得：

C_C＝{1+R_CR_a/[R_S(R_C+R_a)]}^-1 (七)；

C_S＝{1+R_SR_a/[R_C(R_S+R_a)]}^-1 (八)；

其中，R_a、R_s和R_C分别由下式(九)～式(十一)获得：

R_n,c＝FVC·R_n (十二)；

R_n,s＝(1-FVC)·R_n (十三)；

其中，R_n是日净辐射(Wm^-2)；R_n,s是土壤的净辐射(Wm^-2)；R_n,c是植被部的净辐射(Wm^-2)；G是土壤热通量(Wm^-2)，本实施例为日尺度，取值为0；FVC为植被覆盖度，FVC＝(NDVI-NDVI_min)/(NDVI_max-NDVI_min)，NDVI是归一化植被指数，来自于遥感数据，NDVI_max和NDVI_min分别是图像覆盖区NDVI的最大值和最小值，分别取0.9和0.1；

是冠层源区到参考高度处的空气动力学阻抗(sm^-1)；

是植被的整体边界层阻抗(sm^-1)；

是土壤到植被源区域间的空气动力学阻抗(sm^-1)；

是植被冠层阻抗(sm^-1)；

是土壤阻抗(sm^-1)；D是参考高度处饱和水汽压差(kpa)，由再分析气象观测资料获得；C_p是空气比热，取值为1005(Jkg^-1K^-1)；ρ是空气密度，取值为1.29(kgm^-3)；Δ是饱和水汽压斜率(kPaK^-1)；γ是干湿表常数(kPaK^-1)。

其中，

和

按照文献“Shuttleworth,W.James,&Wallace,J.S.1985."Evaporation from Sparse Crops-an Energy Combination Theory."[In English].Quarterly Journal Of the Royal Meteorological Society 111,no.469(1985):839-55.”中所述方法获得。

所述土壤阻抗

和植被冠层阻抗

分别由下式(十四)和式(十五)获得：

其中，θ为表层土壤湿度(cm³·cm^-3)，来自于ESA CCI数据，https://www.esa-soilmoisture-cci.org/；θ_sat土壤饱和含水量(cm³·cm^-3)，来自于土壤质地数据，所述土壤质地数据为北京师范大学(Land-Atmosphere Interaction Research Group at BeijingNormal University,http://globalchange.bnu.edu.cn/)发布的面向陆面过程模型的中国土壤水文数据集(China Dataset of Soil Hydraulic Parameters UsingPedotransfer Functions for Land Surface Modeling)(Dai et al.,2013)。

其中，

LAI_eff＝LAI/(0.3LAI+1.2)；

f(T₀)＝1-((T₀-T_opt)/T_opt)²；

f(R_s,d)＝1-exp(-R_s,d/R_sd,opt)；

其中，r_min为气孔最小阻抗，农田取180，草地取110，森林为250；LAI_eff为有效叶面积指数；LAI为叶面积指数，来自于遥感数据；f(T₀)为温度胁迫因子；T₀为气温(℃)，由再分析气象观测资料获得；T_opt为最优气温(℃)，农田取25，草地取25，森林为20；R_s,d为短波辐射(Wm^-2)，由再分析气象观测资料获得；f(R_s,d)为短波辐射胁迫因子；R_sd,opt为最优短波(Wm^-2)，取值为500；f(D₀)为水汽压胁迫因子；D₀为饱和水汽压差，由再分析气象观测资料获得；K_D0为系数，农田取0.0238，草地取0.024，森林取0.037；f(θ₂)为土壤湿度胁迫因子；θ₂为根区土壤湿度(cm³·cm^-3)，来自于荷兰GLEAM的产品https://www.gleam.eu/；θ_wilt为植被凋萎土壤含水量(cm³·cm^-3)，来自土壤数据资料；θ_cr为田间持水量(cm³·cm^-3)，来自土壤数据资料，所述土壤数据资料为北京师范大学(Land-Atmosphere Interaction ResearchGroup at Beijing Normal University,http://globalchange.bnu.edu.cn/)发布的面向陆面过程模型的中国土壤水文数据集(China Dataset of Soil Hydraulic ParametersUsing Pedotransfer Functions for Land Surface Modeling)(Dai et al.,2013)。

设定上述S-W模型的原始估算的土壤阻抗与优化后土壤阻抗的关系如下式(十六)所示，原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系如下式(十七)所示：

其中，

为原始估算的土壤阻抗，

为优化后土壤阻抗，

为原始植被冠层阻抗，

为优化后植被冠层阻抗，a_s和b_s为土壤阻抗变化系数，a_c和b_c为植被冠层阻抗变化系数。

利用粒子群优化算法，结合步骤1得到的高精度的晴空日的蒸散发，获得a_s、b_s、a_c、b_c，重复步骤1和2共五次，获得最优的a_s为1.44、b_s为100、a_c为0.214、b_c为46。

步骤3，根据获得的上述最优阻抗变化系数，获得优化后的土壤阻抗和植被冠层阻抗，进而获得逐日、真实的土壤阻抗和植被冠层阻抗；

然后，结合遥感观测信息和气象资料，采用式(四)估算得到逐日、真实的蒸散发。

对比例

对比例1

本对比例仅采用温度-植被指数模型估算日蒸散发，采用式

(一)和式(二)进行估算。

对比例2

本对比例仅采用S-W双源模型估算日蒸散发，仅采用式(四)～式(十三)进行，不经过校正优化。

实验例

实验例1

将实施例1(本发明所述模型)和对比例1(温度-植被指数模型)中所述方法估算获得的2009～2011年的日蒸散发数据进行比对，其中，以盈科气象站2009～2011年的观测数据为对照组，所述观测数据采用涡动观测仪获得，上述比较结果如图2所示。

由图2可知，本发明所述模型(实施例1)与温度-植被指数模型(对比例1)估算结果与地面观测值(对照组)在数值和趋势上均接近，但本发明所述模型获取的结果为逐日的，而温度-植被指数模型(对比例1)获取结果为非连续的。

实验例2

将实施例1(本发明所述模型)和对比例2(原始S-W模型)中所述方法估算获得的2009～2011年的日蒸散发数据进行比对，其中，以盈科气象站2009～2011年的观测数据为对照组，所述观测数据采用涡动观测仪获得，上述比较结果如图3所示。

由图3可知，本发明所述模型与原始S-W模型获取结果虽都连续，但本发明所述模型和实测结果(对照组)在数值和趋势都优于原始S-W模型，原始S-W模型存在显著低估。

实验例3

将实施例1(本发明所述模型)和对比例2(原始S-W模型)所述方法估算得到的2012年的日蒸散发数据分别与黑河观测网的地面站点观测数据进行比对，结果分别如图4中的a)和b)所示，图中的EC_1～EC_6指的是图1中的六个涡动观测台站(1-6)，其中，黑河观测网的数据由2012年6个台站组成，数据采用涡动观测仪测定。

由图4中的a)和b)可知，本发明所述模型与6个观测台站数据吻合度(图4中的a))显著优于原始S-W模型(图4中的b))，原始S-W模型存在显著低估。

实验例4

将实施例1(本发明所述模型)、对比例1(温度-植被指数模型)和对比例2(原始S-W模型)在黑河中游灌区的估算结果进行空间比对，结果如图5所示。

由图5可知，本发明所述模型估算结果与温度-植被指数模型估算结果具有相似的空间分布和数量级，原始S-W模型在整个区域存在显著低估。

实验例5

将实施例1(本发明所述模型)和对比例2(原始S-W模型)所述方法对蒸散发的估测数据与在研究区7个验证站点

(EC_1～EC_6、YK)的观测数据进行比对，估测精度结果如表1所示。

表1

其中，R为相关系数，RMSE为均方根误差(mm·d^-1)，Bias为偏差(mm·d^-1)，MAD为平均绝对误差(mm·d^-1)。

由表1可知，与原始S-W模型的估算方法相比，本发明实施例1所述估算方法的估算精度提高约50％。

以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。

Claims (3)

1.一种基于遥感信息的灌溉地区逐日实际蒸散发估算方法，其特征在于，所述估算方法通过将基于过程的蒸散发模型和基于温度的蒸散发模型进行集成来实现；所述估算方法包括以下步骤：

步骤1，采用基于温度的蒸散发模型，获取灌溉地区晴空日的日蒸散发；

步骤2，在晴空日，对基于过程的蒸散发模型进行优化；

步骤3，估算灌溉地区逐日实际蒸散发；

步骤2中，根据步骤1获得的日蒸散发，对基于过程的蒸散发模型进行优化，所述优化为对相关参数进行校正优化，相关参数为土壤阻抗和植被冠层阻抗，

所述基于过程的蒸散发模型为S-W双源模型；

所述基于温度的蒸散发模型为温度-植被指数模型；

所述基于过程的蒸散发模型通过粒子群优化算法集成，对基于过程的蒸散发模型中土壤阻抗和植被冠层阻抗进行校正；

步骤3中，所述灌溉地区逐日实际蒸散发通过基于过程的蒸散发模型估算获得；

其中，所述基于过程的蒸散发模型为优化后的模型。

2.根据权利要求1所述的估算方法，其特征在于，步骤1包括以下子步骤：

步骤1-1，选择晴空日，获得灌溉地区的蒸发比；

步骤1-2，根据获得的蒸发比，获取晴空日的日蒸散发。

3.根据权利要求1所述的估算方法，其特征在于，步骤2中，对土壤阻抗和植被冠层阻抗的优化校正包括以下步骤：

步骤2-1，设定原始估算的土壤阻抗与优化后土壤阻抗的关系式，并设定原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系式；

步骤2-2，获得步骤2-1所述关系式中的土壤阻抗变化系数和植被冠层阻抗变化系数，采用粒子群优化算法获取优化后的土壤阻抗变化系数和植被冠层阻抗变化系数；

步骤2-1中，原始估算的土壤阻抗与优化后土壤阻抗的关系如下式（十六）所示：

（十六）；

原始估算的植被冠层阻抗与优化后植被冠层阻抗的关系如下式（十七）所示：

（十七）

其中，

为原始估算的土壤阻抗，

为优化后土壤阻抗，

为原始植被冠层阻抗，

为优化后植被冠层阻抗，a_s和b_s为土壤阻抗变化系数，a_c和b_c为植被冠层阻抗变化系数。

Images