基于GF-4和MODIS结合的蒸散发遥感反演方法和系统与流程

本发明涉及蒸散发反演技术领域，尤其是涉及一种基于gf-4和modis结合的蒸散发遥感反演方法和系统。

背景技术

蒸散发(evapotranspiration，et)是地-气系统中水分循环的重要环节，也是地表能量循环、水循环和碳循环中最难估算的分量。在干旱、半干旱地区，蒸散发量占到农田生态系统中总耗水量的80％以上，因此，定量计算蒸散发量对于实施例农业水文循环模拟、指导农业水管理具有重要意义。尤其在种植结构比较复杂和破碎的灌区作用更大。

遥感技术的发展使区域作物蒸散量的快速、准确、大面积估算成为可能，其改变了传统单点对地观测的理念，成为当今模拟大范围作物蒸散发的有效途径。随着对地观测技术的快速发展，可获得的遥感数据源越来越多，现有实施例仍普遍局限于基于单一传感器数据进行et的反演。但受云、雨、雾/霾、不同传感器参数及空间分辨率、重访周期、影像获取时复杂地表等不确定性因素影响，单一传感器在作物整个生育期过程中可获得的数据较少、时空连续性较差，造成其et反演结果在区域适应性方面距离精准农业监测业务化还存在较大差距。

技术实现要素：

有鉴于此，本发明的目的在于提供基于gf-4和modis结合的蒸散发遥感反演方法和系统，通过融合不同空间分辨率的多源传感器数据，可以反演得到高精度的et结果。

第一方面，本发明实施例提供了一种基于gf-4和modis结合的蒸散发遥感反演方法，包括：

获取目标区域的gf-4数据和modis数据，并根据所述gf-4数据和所述modis数据反演得到地表关键参数；

基于所述地表关键参数，利用地表能量平衡模型sebal得到均匀地表与大气间的热通量交换数据；

基于所述热通量交换数据，利用能量平衡模型tseb进行分解，得到植被冠层热通量交换数据和土壤热通量交换数据；

根据所述植被冠层热通量交换数据以及所述土壤热通量交换数据进行计算，得到瞬时蒸散发量；

将所述瞬时蒸散发量进行时间日尺度转换，得到所述目标区域的日蒸散发量。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，所述根据所述gf-4数据和所述modis数据反演得到地表关键参数的步骤，包括：

根据所述gf-4数据进行叶面积指数反演，得到叶面积指数；

根据所述modis数据反演得到地表温度、归一化植被指数、地表反照率以及地表比辐射率。

结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述根据所述gf-4数据进行叶面积指数反演，得到叶面积指数的步骤，包括：

基于所述gf-4数据，利用prosail辐射传输模型采用查找表方法进行叶面积指数反演，得到所述叶面积指数。

结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述根据所述植被冠层热通量交换数据以及所述土壤热通量交换数据进行计算，得到瞬时蒸散发量的步骤，包括：

获取植被冠层潜热通量的初始值；

根据能量平衡规则对所述植被冠层热通量交换数据、所述土壤热通量交换数据以及所述初始值进行迭代计算，得到瞬时蒸散发量；

结合第一方面的第三种可能的实施方式，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述获取植被冠层潜热通量的初始值的步骤，包括：

基于所述植被冠层热通量交换数据，根据priestley–taylor公式计算所述植被冠层潜热通量的初始值。

结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第五种可能的实施方式，其中，所述基于所述地表关键参数，利用地表能量平衡模型sebal得到均匀地表与大气间的热通量交换数据的步骤，包括：

基于所述地表关键参数，利用所述地表能量平衡模型sebal得到地表净辐射通量、感热通量和土壤热通量。

结合第一方面的第五种可能的实施方式，本发明实施例提供了第一方面的第六种可能的实施方式，其中，所述基于所述热通量交换数据，利用能量平衡模型tseb进行分解的步骤包括：

将所述地表净辐射通量分解为土壤净辐射通量和植被冠层净辐射通量；以及，将所述感热通量分解为土壤感热通量和植被冠层感热通量。

结合第一方面的第五种可能的实施方式，本发明实施例提供了第一方面的第七种可能的实施方式，其中，所述基于所述地表关键参数，利用地表能量平衡模型sebal得到感热通量的步骤，包括：

基于所述地表温度，并利用所述地表能量平衡模型sebal对空气动力学阻抗进行迭代计算，得到所述感热通量。

第二方面，本发明实施例还提供一种基于gf-4和modis结合的蒸散发遥感反演系统，包括：

数据获取模块，用于获取目标区域的gf-4数据和modis数据，并根据所述gf-4数据和所述modis数据反演得到地表关键参数；

sebal模型计算模块，用于基于所述地表关键参数，利用地表能量平衡模型sebal得到均匀地表与大气间的热通量交换数据；

tseb模型分解模块，用于基于所述热通量交换数据，利用能量平衡模型tseb进行分解，得到植被冠层热通量交换数据和土壤热通量交换数据；

tseb模型计算模块，用于根据所述植被冠层热通量交换数据以及所述土壤热通量交换数据进行计算，得到瞬时蒸散发量；

转换模块，用于将所述瞬时蒸散发量进行时间日尺度转换，得到所述目标区域的日蒸散发量。

结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，所述数据获取模块包括：

第一反演单元，用于根据所述gf-4数据进行叶面积指数反演，得到叶面积指数；

第二反演单元，用于根据所述modis数据反演得到地表温度、归一化植被指数、地表反照率以及地表比辐射率。

本发明实施例带来了以下有益效果：

本发明实施例提供了一种基于gf-4和modis结合的蒸散发遥感反演方法和系统，通过获取目标区域的gf-4数据和modis数据，并根据gf-4数据和modis数据反演得到地表关键参数；基于地表关键参数，利用地表能量平衡模型sebal得到均匀地表与大气间的热通量交换数据；基于热通量交换数据，利用能量平衡模型tseb进行分解，得到植被冠层热通量交换数据和土壤热通量交换数据；根据植被冠层热通量交换数据以及土壤热通量交换数据进行计算，得到瞬时蒸散发量；将瞬时蒸散发量进行时间日尺度转换，得到目标区域的日蒸散发量。结合gf-4数据和modis数据，通过融合不同空间分辨率的多源传感器数据，可以反演得到高精度的et结果。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的基于gf-4和modis结合的蒸散发遥感反演方法流程图；

图2为本发明另一实施例提供的基于gf-4和modis结合的蒸散发遥感反演方法流程图；

图3为本发明实施例提供的基于gf-4和modis结合的蒸散发遥感反演系统示意图；

图4为本发明实施例提供的电子设备示意图。

图标：10-数据获取模块；20-sebal模型计算模块；30-tseb模型分解模块；40-tseb模型计算模块；50-转换模块；1000-电子设备；500-处理器；501-存储器；502-总线；503-通信接口。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

目前，现有遥感反演实施例仍普遍局限于基于单一传感器数据进行et的反演。但受云、雨、雾/霾、不同传感器参数及空间分辨率、重访周期、影像获取时复杂地表等不确定性因素影响，单一传感器在作物整个生育期过程中可获得的数据较少、时空连续性较差，造成其et反演结果在区域适应性方面距离精准农业监测业务化还存在较大差距。基于此，本发明实施例提供的一种基于gf-4和modis结合的蒸散发遥感反演方法和系统，通过融合不同空间分辨率的多源传感器数据，可以反演得到高精度的et结果。

为便于对本实施例进行理解，首先对本发明实施例所公开的一种基于gf-4和modis结合的蒸散发遥感反演方法进行详细介绍。

高分四号(gf-4)卫星是地球同步轨道遥感卫星，搭载了一台可见光50米/中波红外400米分辨率、大于400公里幅宽的凝视相机，采用面阵凝视方式成像，具备可见光、多光谱和红外成像能力。

modis的全称为中分辨率成像光谱仪(moderate-resolutionimagingspectroradiometer)，是搭载在terra和aqua卫星上的一个重要的传感器，是卫星上唯一将实时观测数据通过x波段向全世界直接广播，并可以免费接收数据并无偿使用的星载仪器。在对地观测过程中，每秒可同时获得11兆比特的来自大气、海洋和陆地表面信息，每日或每两日可获取一次全球观测数据。

图1示出了本发明实施例提供的基于gf-4和modis结合的蒸散发遥感反演方法流程图。

如图1所示，本实施例提供了一种基于gf-4和modis结合的蒸散发遥感反演方法，可以应用于区域作物的蒸散发遥感反演，包括以下步骤：

步骤s101，获取目标区域的gf-4数据和modis数据，并根据gf-4数据和modis数据反演得到地表关键参数；

具体地，gf-4数据包括目标区域的gf-4卫星影像，可以获取地表反射率，modis数据包括反射率产品mod09a1和温度发射率产品mod11a2，可以得到另一种地表反射率和地表温度。

在蒸散发计算中，遥感反演的参数主要包括：地表温度，是地表-大气界面达到热量平衡的结果，决定了地表长波辐射能量的大小；地表反照率，决定了地表有效能量的大小；地表比辐射率，是地表温度反演的最基本参数；植被指数，是遥感监测地表植被分布植被盖度的指示参数。本实施例的地表关键参数还包括叶面积指数，是指植物叶片总面积与土地面积的比值。它与植被的密度、结构(单层或复层)、树木的生物学特性(分枝角、叶着生角、耐荫性等)和环境条件(光照、水分、土壤营养状况)有关，是表示植被利用光能状况和冠层结构的一个综合指标。

进一步地，步骤s101中根据gf-4数据和modis数据反演得到地表关键参数的步骤，包括：根据gf-4数据进行叶面积指数反演，得到叶面积指数；根据modis数据反演得到地表温度、归一化植被指数、地表反照率以及地表比辐射率。

(一)植被指数

植被指数是对植被盖度、生物量以及生长活力等具有指示意义的数值，是遥感监测地表植物生长与分布的方法之一。它通过对近红外与红光波段反射率的线性或非线性组合，不需要其他辅助资料，来消除影像中土壤光谱的影响，实现对植物信息状态的表达。归一化植被指数(ndvi)是植被生长状态以及植被覆盖度的最佳指示因子，可以消除部分与太阳高度角、卫星观测角、大气条件、云/阴影、地形和有关辐照度条件变化等的影响。

本实施例使用的数据是美国航空航天局nasa接收并经过处理得到的modis数据产品(hdf格式)。对于modis影像而言，ndvi是用近红外波段(第1波段)和红光波段(第2波段)数值之差α2-α1和这两个波段数值之和α1+α2的比值。如公式(1)：

对于陆地表面的主要覆盖物而言，雪、云、水在红光波段比近红外波段有较高的反射作用，因此其ndvi值为负值，裸岩在两波段有相似的反射作用，因而其ndvi值近于0，在有植被覆盖的情况下，ndvi为正值，且随着植被覆盖度的增大而增大。

植被覆盖度通常定义为植物群落总体绿叶的垂直投影面积与样方总面积之比。它是影响地表植被蒸腾和土壤水分蒸发损失过程的重要控制因子，是植被冠层形状、植被空间分布、叶子倾角、及重叠所形成的参量，与植被的光谱特征无关，但利用遥感手段反演植被覆盖度时，主要通过植被指数来估计，具体通过式(2)：

式中，ndvi是植被指数，ndviv和ndvis分别是茂密植被覆盖和完全裸土像元的ndvi值。对于modis数据，通常取ndviv＝0.9，ndvis＝0.15。因此，当ndvi>ndviv＝0.9时，fv＝1，表示该像元是一个茂密植被覆盖的地区，看不见裸露的土壤表面；否则，当ndvi<ndvis＝0.15时，fv＝0，表示该像元是一个完全裸露的地区，没有任何植被覆盖。

(二)地表反照率

地表反照率(albedo)是地表不同类型物体反射太阳辐射能量与入射到该物体的太阳总能量的比。表征地表不同下垫面对太阳辐射的反射能力，决定着能量在地-气间的分配比率，是影响地表能量收支平衡以及生态过程的最重要参数之一。albedo随着时空变化而变化，在一般情况下变化较缓慢，但是在气候环境突变(沙尘暴、暴风雪等)的情况下其值变化明显。地表反照率受下垫面状况、太阳天顶角、入射辐射的光谱分布等因子的影响，利用大气辐射传输模型，采用模拟方法对九种多波段卫星，包括avhhr、landsattm/etm+、modis、aster等建立了有关求取宽波段的方程组，精度约可达到0.02，本实施例可以采用梁顺林等(2000)提出的利用modis数据的窄波段反射率来推算宽波段的反照率算法，如公式(3)：

α＝0.160α1+0.291α2+0.243α3+0.116α4+0.112α5+0.081ε7-0.0015(3)

其中，α是宽波段的地表反照率，α1,2,3,4,5,7分别为经过大气纠正后的对应各波段的反射率。

(三)地表比辐射率

地表比辐射率(ε)，又称地表发射率，是物体与黑体在同温度、同波长下的辐射出射度的比值，是一个无量纲的量，ε的取值在0～1之间。比辐射率是地表物体向外辐射电磁波的能力，它受地表物体的组分、表面状态、温度、物体辐射能的波长、观测角度以及其他物理性质(介电常数、含水量等)都有密切关系。实施例表明，地表比辐射率存在0.01的误差，将导致地表温度的估算误差达到1-2℃。因此，估算地表比辐射率对地表温度的反演至关重要。

地球表面不同区域的地表结构虽然复杂，但从卫星象元尺度来看，其基本分为3中类型：水体、城镇和自然表面。城镇象元一般在影像中占得比例不大。而相反，自然表面通常占影像象元比例最大，是地表温度反演中需要考虑的重点。valor认为地表可以简单地看作是由不同比例的植被叶冠和裸土所混合组成的。因此本实施例，根据地表发射率的不同贡献，建立modis影像的地表比辐射率的估算公式，如公式(4)：

εi＝fvrvεiv+(1-fv)rsεis+dε(4)

其中，εi是modis影像第i波段(i＝31,32)的地表比辐射率；εiv和εis分别代表植被和裸土在第i波段的地表辐射率，一般取ε31v＝0.98672，ε32v＝0.98990，ε31s＝0.96767，ε32s＝0.97790；fv是植被覆盖度；dε是热辐射相互作用校正，由植被和裸土间的热辐射相互作用导致；rv和rs分别是植被和裸土的辐射比率，被定义为公式(5)和公式(6)：

rv＝bv(tv)/b(t)(5)

rs＝bs(ts)/b(t)(6)

式中，bv(tv)和bs(ts)分别是混合象元内植被和裸土的热辐射强度，是混合像元的热辐射强度。

通过模拟计算结果显示，rv和rs主要取决于温度变化和植被覆盖度，且植被覆盖度的影响更大，因此，建立rv、rs与fv的关系进行估算。如公式(7)和公式(8)：

rv＝0.92762+0.07033fv(7)

rs＝0.99782+0.08362fv(8)

最后对热辐射相互作用校正项dε进行估算，由于dε在植被与裸土各占一半时达到最大值，因此，可以根据如下经验公式来估计：

当fv＝0或fv＝1时，dε最小，dε＝0；

当0<fv<0.5时，dε＝0.003796fv；

当0.5<fv<1时，dε＝0.003796(1-fv)；

当fv＝0.5时，dε最大，dε＝0.001898；

用上述公式计算得到的εi若大于εiv，则取εi；若εi小于εis，则取εi＝εi。

进一步地，根据gf-4数据进行叶面积指数反演可以利用prosail辐射传输模型采用查找表方法进行叶面积指数反演，得到叶面积指数。

叶面积指数(leafareaindex，lai)反演的物理模型较多，但大都受制于有限数量的输入参数。prosail辐射传输模型以其对植被冠层光谱反射率模拟和lai、植被生理生化参数估算的稳定性和鲁棒性，成为应用较多的模型之一，得到了广泛的验证与应用。

prosail模型是由冠层反射率模型sail和叶片光学特性模型prospect耦合而成。prospect模型是从叶片尺度模拟叶片方向半球反射与透射的辐射传输模型。prospect是基于allen等的平板模型改进而来，可以表征新鲜叶片从400nm-2500nm的光学特性。模型的输入参数包括叶片结构参数(n)、叶绿素a+b的含量(cab)、等效水厚度(cw)、干物质含量(cm)等。sail模型是把植被当做浑浊介质，假设叶片方位角分布均匀，考虑任意倾角，冠层双向反射率作为观测角度的一个函数来描述植被冠层二向性反射率的辐射传输模型。sail是最早用于模拟冠层反射率的模型之一，随着发展模型又考虑了热点效应的影响。sail模型主要的输入参数有lai、平均叶倾角(alia)、热点参数(hot)、天空光散射比(skyl)、土壤亮度参数、土壤背景反射率、太阳天顶角、观测天顶角、相对方位角、叶片反射率和透射率。耦合后的prosail模型就可以把prospect模拟的叶片反射率和透射率作为sail的输入参数，来模拟可见光-近红外光谱的方向冠层反射率。

prosail模型输入参数范围及分布如表1。所有参数的输入范围及分布范围参考了以往的实施例。叶绿素含量、生物量、等效水厚度和平均叶倾角根据实测数据。天空光散射比由于对冠层反射率影响非常小，根据以往实施例，可以设定为0.1。土壤反射率光谱根据实地观测的值进行平均。太阳天顶角、观测天顶角、相对方位角可以直接从影像头文件中直接得到。

表1prosail模型输入参数范围及分布

查找表方法(lut)是解决模型反演问题最简单、有效的方法。克服了数值迭代优化造成的计算效率低下问题，也克服了人工神经网络方法黑箱操作造成的不易分析问题，适合于lai的反演。依据模型的输入参数范围，实施例首先对参数进行随机采样；然后结合遥感影像获取的天顶角、方位角等信息都输入到prosail模型中进行模拟窄波段反射率；利用不同传感器数据的光谱响应函数模拟不同传感器在可见光、近红外波段的宽波段反射率；生成叶面积指数与不同宽波段反射率构成的查找表。其中宽波段反射率可以利用以下公式(9)计算：

式中，ρs(λ)是模拟的传感器宽波段反射率，λminandλmax分别是波段波长的上、下界。ρs(λi)是窄波段第i波段的反射率，φ(λi)是窄波段第i波段的光谱响应函数。

为使查找表更有代表性、足够准确，查找表的大小级别应该足够大。因此，查找表的大小级别设置为200,000。为了避免观测误差或同物异谱、异物同谱造成的病态反演问题，反演策略应用最小化均方根误差代价函数方法。当代价函数值越小，模型模拟反射率值与遥感数据实际反射率越接近，则当代价函数值最小或较小时所对应的模拟lai值即被认为是反演结果。实施例选取的代价函数是相对均方根误差，参照公式(10)：

其中，m为波段数。

prosail模型模拟反射率与遥感反射率代价函数最小的条件下，对应输入参数的一组解。但由于观测误差和模型的不完整性，得到的解的结果并不一定是最优解、也可能不唯一。为克服这个问题，实施例根据相对均方根误差的值按从小到大排列参数组合，选取查找表中10％查找表个数的最小相对均方根误差值对应的反演参数组合的均值作为最终的lai估算结果。

步骤s102，基于地表关键参数，利用地表能量平衡模型sebal(surfaceenergybalancealgorithmforland，sebal)得到均匀地表与大气间的热通量交换数据；

进一步地，步骤s102包括：基于地表关键参数，利用地表能量平衡模型sebal得到地表净辐射通量、感热通量和土壤热通量。

具体地，基于sebal模型求取的相关通量的过程如下：

利用遥感手段来估算区域的蒸散发量，其基本思想是忽略植被光合作用耗能和水平方向上的能量输入，以能量平衡方程为基础建立不同的遥感蒸散模型，如sebal、地表能量平衡系统sebs(surfaceenergybalancesystem，sebs)模型等，他们都充分利用遥感反演的地表参数，结合少量的常规气象数据进行区域的地表蒸散反演，不同的遥感模型对气象要素以及实施例对象的特点不同，使其具有不同的实用性和反演精度。能量平衡方程可表示为公式(11)：

rn＝h+g+λet(11)

式中，rn为地表净辐射通量(w/m²)，是地表接收到的能量总和；h为地表与大气的感热通量(也称感热通量)(w/m²)；g为地表的土壤热通量(w/m²)，即下垫面与土壤间的热通量交换，用于下垫面地表升温；λet为潜热通量(w/m²)，即下垫面与大气间水汽的热交换，其中λ为水的汽化潜热，et为蒸散量。因此只要得到rn、g、h的值，就可以求出λet值，进而求算出蒸散发量et。

地表净辐射通量是地表能量、物质输送与交换的原动力，也是气候形成及变换的主要依据。地表净辐射通量可表述为地表短波和长波辐射的收入与支出的差值，可根据辐射平衡方程由入射能量减去出射能量求得，如公式(12)：

rn＝(1-α)rs↓+(rl↓-rl↑)-(1-ε)rl↓(12)

式中，rs↓为入射到地表的太阳短波辐射w/m²，rl↓为入射至地表的大气长波辐射w/m²，rl↑为反射至大气的长波辐射w/m²，α为地表反照率，ε为地表比辐射率，其中：

rl↓＝1.08(-lnτ)^0.265σta⁴

rl↑＝εσts⁴

式中，gsc是太阳常数(1367)；τ为大气透过率；θ为太阳天顶角；dr是以天文单位表示的日地距离，σ是斯蒂芬-玻尔兹曼常数(取5.67×10^-8)；t0是近地表气温的摄氏温度值；ta为空气温度的热力学温度值；ts为反演的地表温度；j是卫星影像获取日期在太阳历中的排列序号。

土壤热通量是指由传导作用而存储在土壤和植被中的能量，它对土壤蒸发、地表能量交换均有影响，相对于净辐射、潜热和感热通量比较小，但却是一重要的分量。土壤热通量一般难以利用遥感手段直接获取，但可以根据与地表温度、净辐射通量、地表反照率和归一化植被指数的经验统计关系进行估算求取，在计算过程中，sebal模型和tseb模型都有各自的利用以上参数而建立的估算公式，其中利用sebal模型进行土壤热通量求解的计算公式为(13)：

其中，c11表示卫星过境时间对g的影响，根据卫星过境时间进行适当的校正，过境时间在当地时间12点以前时c11取0.9；过境时间在当地时间12点到14点之间取1.0；在14点到16点之间取1.1。

进一步地，基于地表关键参数，利用地表能量平衡模型sebal得到感热通量的步骤，包括：基于地表温度，并利用地表能量平衡模型sebal对空气动力学阻抗进行迭代计算，得到感热通量。

具体地，感热通量，又称显热通量，是下垫面与大气间湍流形式的热交换，主要是下垫面由能量传导和大气对流作用而散失的能量。感热通量的计算相当复杂，受土地利用、土壤结构和组成、植被覆盖、地形等多种因素影响，一般通过如下公式(14)估算：

式中，ρ为空气密度(kg/m³)；cp为标准大气压下的空气比热(通常取1004j/kg·℃)；rah为空气动力学阻抗(s/m)，可用风速廓线理论计算；tr为空气动力学温度(k)，ta为大气平均作用温度。

通过公式(14)可知，感热通量与空气动力学阻抗、温度梯度密切相关，而利用遥感方法很难直接测得空气动力学温度和阻抗。为求得温度梯度，sebal模型引入两个极端的“锚点”(干点、湿点)作为边界条件，利用“湿点”的感热通量近似为零，蒸散量达到最大；而“干点”的蒸散量近似为零，感热通量达到最大的假设，进行温度梯度的计算。由于大气的稳定性对空气动力学阻抗影响明显，为消除地表热量引起的浮性效应，sebal模型引入monin-obukhov理论，通过复杂的循环迭代过程，得到稳定大气条件下的摩擦风、动力学阻抗等参数，而rah在校正过程又需要感热通量，因此采用迭代方法可求解得到稳定的感热通量h。

在感热通量计算中，干湿点选取的好坏直接影响最终的感热通量精度，利用温度-植被指数(ts-vi)的空间特征关系可进行干湿点的快速选取，一般认为水分供应充足、植被生长茂盛、温度很低、处于潜在蒸散水平的像元，可以是植被生长良好的完全覆盖的区域或开放的水体是湿(cold)点最优的候选区；而相对的没有植被覆盖的干燥的闲置农田或盐碱地，温度很高，蒸散量几乎为0的像元，是干点(“hot”)的最优候选区。且注意，在选取过程中，避免选取极端象元，因为sebal模型中温度梯度是建立在自然植被基础上的，若选择极低温度象元中的阴影或者云等“伪冷点象元”或者选择极高温度象元的城市水泥地等“伪干点象元”，则必然造成估算误差。

在感热通量的循环迭代计算中，需要首先计算动力粗糙度、摩擦风速、空气动力学阻抗；随后假设温度梯度dt与ts间存在线性关系，如公式(15)：

dt＝ats+b(15)

然后，根据干点、湿点温度梯度以及ts的数值，求取线性参数a和b，如公式(16)和公式(17)，进而求得整个影像各个象元温度梯度的分布。

计算求得参数a和b后，可以得到温度梯度dt，带入公式(18)可以求取感热通量h：

但在实际情况下，大气稳定性对空气动力学阻抗的影响明显，为消除大气的浮性效应，在此引入monin-obukhov相似理论，monin-obukhov长度的计算公式如(19)：

其中，u*是摩擦风速(m/s),k是vonkarman常数0.41，g是重力加速度(9.81m/s²)。l为monin-obukhov长度，它决定了大气的稳定状况，若l大于0，表示空气是稳定的；若l小于0，则空气是不稳定的；若l等于0，代表不需考虑大气的浮性效应。

如果l>0，则修正项ψm(200m)、ψm(0.1m)、ψm(2m)分别按照如下公式计算：

如果l<0则修正项ψm(200m)、ψm(0.1m)、ψm(2m)分别按照如下公式计算：

其中：

若l＝0，则ψm＝0，ψh＝0。将修正值带入u*和rah的修正公式(20)和公式(21)中，重新计算新的u*和rah：

其中，zom是动量表面粗糙度，利用新计算的u*和rah，根据公式(18)和公式(19)重新计算h和l，不断迭代计算，直到rah稳定为止(可以取rah变化幅度低于3％以下)。

步骤s103，基于热通量交换数据，利用能量平衡模型tseb(two-sourceenergybalance，tseb)进行分解，得到植被冠层热通量交换数据和土壤热通量交换数据；

进一步地，步骤s103包括：将地表净辐射通量分解为土壤净辐射通量和植被冠层净辐射通量；以及，将感热通量分解为土壤感热通量和植被冠层感热通量。

步骤s104，根据植被冠层热通量交换数据以及土壤热通量交换数据进行计算，得到瞬时蒸散发量；

进一步地，如图2所示，步骤s104包括：

步骤s1041，获取植被冠层潜热通量的初始值；

具体地，基于植被冠层热通量交换数据，根据priestley–taylor公式计算植被冠层潜热通量的初始值；

步骤s1042，根据能量平衡规则对植被冠层热通量交换数据、土壤热通量交换数据以及初始值进行迭代计算，得到瞬时蒸散发量。

具体地，基于tseb模型求取的相关通量的过程如下：

tseb模型虽也是同单层模型一样都以能量平衡方程为基础，但与单层模型明显不同，在tseb模型中，地表不再是单层模型中的单一、均匀下垫面，而是对土壤和植被进行了分解，分别基于各自的能量平衡方程进行求解，如公式(22)：

rns＝hs+lets+g

rnc＝hc+letc(22)

式中，rns和rnc分别为土壤净辐射通量和植被冠层净辐射通量；hs和hc分别为土壤感热通量和植被冠层感热通量；lets和letc分别为土壤潜热通量和植被冠层潜热通量。

而对于rns和rnc的计算，在tseb模型中可根据太阳高度角和叶面积指数，经验的估算植被冠层净辐射通量，再由总净辐射与植被冠层净辐射之差来求取土壤净辐射通量，如公式(23)：

rns＝rn-rnc(23)

其中，lai为叶面积指数，θ为太阳天顶角，k为是冠层衰减系数。

在土壤热通量的计算过程中，sebal和tseb模型都有各自的经验估算公式，都是根据土壤通量与净辐射的关系，建立估算公式，由于sebal模型中土壤热通量的计算在有植被覆盖时估算精度较高，当植被覆盖稀少、甚至裸土时，利用tseb的经验估算相对较准确，公式为：g＝0.15rn。

priestley-taylor公式可以估算湿润饱和下垫面的蒸散状况，在tseb模型中，根据priestley–taylor公式计算植被冠层潜热通量的初始计算值，如公式(24)：

式中，αpt为priestley-taylor常数(取1.26)，fg是绿叶面积占总叶面积指数的部分；δ为饱和水汽压-温度曲线在植被冠层温度处的斜率，s是干湿球常数(取66.0)。

在串联双层模型中，植被和土壤两层耦合，使得在混合高度处的空气温度等于空气动力学温度，整个感热通量的值是土壤感热通量和植被感热通量之和，把sebal模型计算得到的感热通量h分解为土壤感热通量hs和植被冠层的感热通量hc，如公式(25)：

h＝hc+hs(25)

综上，可知分解的感热通量和潜热通量是未知数，其余由sebal模型计算得出为已知数，利用迭代运算，可以得到：

hc＝rnc-letc

hs＝h-hc

lets＝rns-g-hs

根据实际的情况，估算土壤和植被的潜热通量的象元都应当是非负值，所以根据计算的结果，设所有土壤潜热通量为负的象元为0，再进行迭代计算：

hs＝rns-g-lets＝rns-g

hc＝h-hs

letc＝rnc-hc

当植被潜热通量象元均为非负值时停止计算，否则继续计算，知道所有的土壤和植被潜热通量象元均为非负值时停止，从而分别得到瞬时的土壤和植被的潜热通量值，两者之和即为最终的潜热通量。

步骤s105，将瞬时蒸散发量进行时间日尺度转换，得到目标区域的日蒸散发量。

具体地，利用遥感估算的地表热通量只能代表卫星过境时刻反演得到的地面瞬时值，而水文、生态、气候等领域对地表热通量的实施例在时间尺度上至少是日尺度，还有更长的月、年尺度，因此，在实际应用中需对蒸散估算的结果进行时间尺度的拓展。

天气晴朗时，农田等区域上的太阳辐射平衡各分量和蒸散发量的日变化曲线呈正弦变化，而相应的太阳辐射通量密度的日变化也存在正弦关系。因此，由瞬时蒸散速率推求日蒸散的计算公式可以表示为(26)：

式中，ne为日蒸散发时数，在实际日出后、日落前的一小时，由于蒸散量极少，可以忽略，因此实际的ne可以用日照时数减去2小时得到；t为从日出算起到卫星获取遥感数据时刻的时间间隔；由于从气象站仅得到每天的日照时数，没有确切的日出时刻，本实施例根据日照时数和卫星过境时刻进行经验估算，如公式(27)：

正弦函数法只需要知道日照时数以及日出时间即可以得到当天的日蒸散量，方法简便，因此本实施例利用正弦函数法进行日蒸散量拓展。

单层sebal模型物理概念清晰，对气象资料要求少，在地表辐射温度和空气动力学温度换算中利用循环迭代运算，避免了人为经验带来的误差；串联tseb模型也物理概念清晰，分离了土壤、植被冠层两湍流源，可以得到不同地表来源的蒸散发量。本实施例将两模型耦合，利用sebal模型求取均匀地表与大气间的热通量交换，然后考虑均匀地表的内部结构，利用求取的总热通量分解为土壤和植被冠层热通量。

本实施例既解决了单层sebal模型无法区分植被蒸腾和土壤蒸发的问题，又避免了串联tseb模型中土壤阻抗的复杂求解，且使tseb继承了sebal模型的稳健性，适合于et反演。本实施例基于的sebal和tseb耦合的et定量反演模型，首次实现了国产高分四号gf-4卫星和modis产品结合的区域作物et反演，与单纯modis产品的et反演结果相比，精度得到了明显的改善。

如图3所示，本实施例提供了一种基于gf-4和modis结合的蒸散发遥感反演系统，包括数据获取模块10、sebal模型计算模块20、tseb模型分解模块30、tseb模型计算模块40以及转换模块50：

数据获取模块10，用于获取目标区域的gf-4数据和modis数据，并根据gf-4数据和modis数据反演得到地表关键参数；

sebal模型计算模块20，用于基于地表关键参数，利用地表能量平衡模型sebal得到均匀地表与大气间的热通量交换数据；

tseb模型分解模块30，用于基于热通量交换数据，利用能量平衡模型tseb进行分解，得到植被冠层热通量交换数据和土壤热通量交换数据；

tseb模型计算模块40，用于根据植被冠层热通量交换数据以及土壤热通量交换数据进行计算，得到瞬时蒸散发量；

转换模块50，用于将瞬时蒸散发量进行时间日尺度转换，得到目标区域的日蒸散发量。

进一步地，数据获取模块10包括第一反演单元和第二反演单元：

第一反演单元，用于根据gf-4数据进行叶面积指数反演，得到叶面积指数；

第二反演单元，用于根据modis数据反演得到地表温度、归一化植被指数、地表反照率以及地表比辐射率。

本发明实施例提供的基于gf-4和modis结合的蒸散发遥感反演系统，与上述实施例提供的基于gf-4和modis结合的蒸散发遥感反演方法具有相同的技术特征，所以也能解决相同的技术问题，达到相同的技术效果。

本发明实施例还提供一种电子设备，包括存储器、处理器，存储器中存储有可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述实施例提供的基于gf-4和modis结合的蒸散发遥感反演方法的步骤。

本发明实施例还提供一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序，计算机程序被处理器运行时执行上述实施例的基于gf-4和modis结合的蒸散发遥感反演方法的步骤。

参见图4，本发明实施例还提供一种电子设备1000，包括：处理器500，存储器501，总线502和通信接口503，处理器500、通信接口503和存储器501通过总线502连接；存储器501用于存储程序；处理器500用于通过总线502调用存储在存储器501中的程序，执行上述实施例的基于gf-4和modis结合的蒸散发遥感反演方法。

其中，存储器501可能包含高速随机存取存储器(ram，randomaccessmemory)，也可能还包括非不稳定的存储器(non-volatilememory)，例如至少一个磁盘存储器。通过至少一个通信接口503(可以是有线或者无线)实现该系统网元与至少一个其他网元之间的通信连接，可以使用互联网，广域网，本地网，城域网等。

总线502可以是isa总线、pci总线或eisa总线等。所述总线可以分为地址总线、数据总线、控制总线等。为便于表示，图4中仅用一个双向箭头表示，但并不表示仅有一根总线或一种类型的总线。

其中，存储器501用于存储程序，处理器500在接收到执行指令后，执行所述程序，前述本发明实施例任一实施例揭示的流过程定义的装置所执行的方法可以应用于处理器500中，或者由处理器500实现。

处理器500可能是一种集成电路芯片，具有信号的处理能力。在实现过程中，上述方法的各步骤可以通过处理器500中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器500可以是通用处理器，包括中央处理器(centralprocessingunit，简称cpu)、网络处理器(networkprocessor，简称np)等；还可以是数字信号处理器(digitalsignalprocessing，简称dsp)、专用集成电路(applicationspecificintegratedcircuit，简称asic)、现成可编程门阵列(field-programmablegatearray，简称fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。结合本发明实施例所公开的方法的步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器501，处理器500读取存储器501中的信息，结合其硬件完成上述方法的步骤。

本发明实施例所提供的进行基于gf-4和modis结合的蒸散发遥感反演方法的计算机程序产品，包括存储了处理器可执行的非易失的程序代码的计算机可读存储介质，所述程序代码包括的指令可用于执行前面方法实施例中所述的方法，具体实现可参见方法实施例，在此不再赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。

在本申请所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个处理器可执行的非易失的计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read-onlymemory)、随机存取存储器(ram，randomaccessmemory)、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。