CN103810387A - 基于modis数据的地表蒸散发全遥感反演方法及系统 - Google Patents

Abstract

本发明涉及一种基于MODIS数据的地表蒸散发全遥感反演方法及系统，其步骤包括：1）确定研究区域，获取研究区域相关MODIS数据产品并进行预处理；2）利用经预处理的MODIS数据产品构建空间变化信息散点图，并进行像元尺度的地表蒸发比EF的全遥感反演；3）根据相关MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；4）根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度的地表蒸散发LE的全遥感反演。本发明仅需要MODIS数据作为输入，克服了目前蒸散发遥感反演一般需要较多地面观测数据辅助的问题，可用于无资料或资料缺乏地区的蒸散发遥感反演。

Description

基于MODIS数据的地表蒸散发全遥感反演方法及系统

技术领域

本发明涉及地表蒸散发反演方法及系统，尤其涉及一种基于MODIS数据的地表蒸散发全遥感反演方法及系统。

背景技术

陆面蒸散发是地表能量平衡和水量平衡的重要组成部分，直接影响土壤–植被–大气系统中的水、热传输。在农学、水文学、气象学等研究领域，如田间尺度的农田灌溉和耕作制度管理、区域尺度的地表水和地下水预测和估算、全球尺度的大气环流和气候变化模拟等方面，蒸散发都发挥着重要的作用。几乎所有水问题的解决，都离不开蒸散发的研究。

传统的蒸散发估算方法，大多基于站点观测的气象数据，并依赖于不同下垫面的先验知识，如水文、气象等领域广泛使用的Penman-Monteith方程。早期的此类蒸散发模型大多只适用于水分充足的下垫面，并且通常只能得到“点”尺度上的蒸散发，无法满足目前日益增长的学科交叉研究的需要，限制了在无资料或资料缺乏地区的应用。近年来，随着各种新兴技术的发展以及蒸散发估算理论的不断完善，基于遥感技术的蒸散发估算已成为水文、气象、农业等领域的重点热门问题。目前，基于遥感数据的蒸散发估算模型主要依靠地表能量平衡原理实现。根据土壤和植被组分是否分开考虑，基于地表能量平衡原理的蒸散发估算模型主要分为一源和二源能量平衡模型。这两类能量平衡模型具有如下缺点：

一源能量平衡模型需要一定的地面观测数据作为输入，对输入数据的精度较为敏感，而且没有考虑地表温度反演中方向性对蒸散发估算的影响。具体来说，一源模型一般需要像元尺度的近地面气温、风速、湿度、植被高度等气象和植被数据作为输入，而目前的遥感技术无法直接反演获得这些参数。通过空间插值得到的区域地面观测气象数据在非均一地表上会产生较大的误差。此外，由于空气动力学温度难以进行测量，通常利用遥感地表温度代替空气动力学温度进行地表温度梯度的计算。对于全植被覆盖地表，空气动力学温度和辐射地表温度之间的较小差异可能不会造成较大的通量估算误差。然而，对于部分植被覆盖地区，该温差可能会达到10°C以上，从而造成非稳定大气条件下显热通量的高估。由于动量和能量传输机制的区别，当利用遥感地表温度进行近地表温差和显热通量的估算时，考虑到能量传输比动量传输具有更高的阻抗，需要在空气动力学阻抗的基础上增加一项剩余阻抗。显然，这种对空气动力学阻抗的校正无法在所有的情况下开展，造成一源模型的应用受到严重的限制。另外，传感器定标、大气校正以及地表比辐射率的确定均会对依赖于绝对地表温度或地气温差的一源能量平衡估算方法产生致命的影响。虽然一源地表能量平衡模型可以提供可靠的地表显热通量和潜热通量估算，但是它们一般需要田间校正，因而很难在广泛变化的下垫面条件下应用。

二源能量平衡模型虽然在模型结构、模型参数化方面较一源能量模型具有一定的优势（如降低了对输入数据的敏感性、考虑了卫星观测几何的影响等），但仍需要大量的近地面观测气象和植被数据作为输入，限制了它在无资料或资料缺乏地区的应用。

为克服一源和二源能量平衡遥感蒸散发反演模型遇到的常规地面气象和植被数据在卫星像元尺度上难以获取的问题，我们利用MODIS遥感数据发明一种地表蒸散发的全遥感反演方法。

发明内容

本发明所要解决的技术问题是针对现有技术的不足，提供一种基于MODIS数据的地表蒸散发全遥感反演方法及系统。

本发明解决上述技术问题的技术方案如下：一种基于MODIS数据的地表蒸散发全遥感反演方法，包括如下步骤：

步骤1：确定研究区域，获取研究区域相关MODIS数据产品并进行预处理；

步骤2：利用经预处理的数据产品构建空间变化信息散点图，并进行像元尺度的地表蒸发比EF的全遥感反演；

步骤3：利用经预处理的数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；

步骤4：根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度的地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF。

本发明的有益效果是：本发明避免了在蒸发比遥感反演时一般需要对地表温度进行准确大气校正的要求；本发明也避免了在利用现有一源能量平衡模型进行蒸散发和蒸发比遥感反演时需要对空气动力学阻抗进行较为复杂参数化的要求，降低了在使用遥感地表温度代替空气动力学温度对蒸散发和蒸发比估算所产生的不确定性；本发明中地表蒸散发反演模型的输入除了由MODIS遥感传感器反演的地表参数外，不需要其他近地面辅助观测数据，解决了过去蒸散发遥感反演需要依赖较多地面气象和植被数据的问题，首次实现了区域地表蒸散发的全遥感反演，使无资料地区和资料缺乏地区地表蒸散发的估算成为可能，在水文、农业、气候变化等研究领域展现了广阔的应用前景。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，步骤1中所述相关MODIS数据产品包括地表温度/比辐射率产品、归一化植被指数产品、定标辐射产品、几何定位产品和大气可降水产品。

进一步，步骤1中对研究区域内相关MODIS数据产品进行预处理的过程如下：

步骤1.1：在表层相对土壤水分和植被覆盖度全范围变化的条件下，确定研究区域；

步骤1.2：对该研究区域内相关MODIS数据产品进行投影转换和剪切处理，进而将处理后的MODIS数据产品数值转换为实际的物理量；

步骤1.3：利用地表温度/比辐射率产品中的质量控制文件剔除研究区内受云污染的像元数据；

步骤1.4：利用几何定位产品中的数字高程数据，设定高程阈值，剔除受地形冷却效应影响的像元数据。

进一步，步骤2中利用经预处理的数据产品构建空间变化信息散点图，进行像元尺度的地表蒸发比的全遥感反演的过程为：

步骤2.1：利用该研究区域的MODIS数据产品中的归一化植被指数NDVI估算该研究区域每个像元的地表植被覆盖度F_r，计算公式为

$F_r={\left(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^2$

其中，NDVI为每个像元的归一化植被指数，NDVI_min和NDVI_max分别为对应裸土和全植被覆盖的归一化植被指数；

步骤2.2：利用该研究区域的MODIS数据产品中每个像元的地表温度T和估算得到的植被覆盖度F_r，构建其以植被覆盖度F_r为横坐标，以地表温度T为纵坐标的空间变化信息散点图；

步骤2.3：通过鲁棒的线性回归方程回归空间变化信息散点图中地表表层相对土壤水分为0时不同植被覆盖度像元的地表温度T_max,_i；同时，设定地表表层相对土壤水分为1时不同植被覆盖度的像元具有相同的地表温度T_min,_i，

T_max，i=T_smax+F_r(T_smin-T_smax)

T_min，i=T_smin

其中，T_smax为在植被覆盖度和地表表层相对土壤水分均为0时的地表温度，T_smin指在植被覆盖度为1时的地表温度；

步骤2.4：利用线性内插法估算空间变化信息散点图中不同土壤水分和植被覆盖度条件下每个像元的蒸发比指示因子

计算公式为，

其中，T_max,i和T_min,i分别指在植被覆盖度为i时，地表表层相对土壤水分为0和1时的地表温度，T_si指在植被覆盖度为i时某一像元所对应的地表温度，

指在植被覆盖度为i时，地表表层相对土壤水分为1时的蒸发比指示因子，

指在植被覆盖度为i时，地表表层相对土壤水分为0时的蒸发比指示因子；

步骤2.5：利用地表温度T计算空间变化信息散点图内每个像元的温度-饱和水汽压曲线斜率Δ，计算公式为

$\mathrm{\Δ}=4098\×\frac{0.6108\×\exp \left[\frac{17.25\×(T-273.15)}{(T-273.15)+237.3}\right]}{{[(T-273.15)+237.3]}^2}$

步骤2.6：利用温度-饱和水汽压曲线斜率Δ和蒸发比指示因子计算像元尺度的蒸发比EF，计算公式为

其中，γ为湿度计常数，为0.066kpa/℃。

进一步，步骤3中根据相关MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演的过程为：

步骤3.1：利用该研究区域MODIS数据产品中的太阳天顶角产品、大气可降水产品和大气顶部窄波段反射率产品估算地表短波净辐射R_sw，计算公式为，

$R_{\mathrm{sw}}=\frac{E_0\cos {\mathrm{\θ}}_s}{D^2}(\mathrm{\α}\′-\mathrm{\β}\′r)$

其中，

$\begin{array}{c}\mathrm{\α}\′=1-\frac{a_1}{\mathrm{\μ}}-\frac{a_2}{{\mathrm{\μ}}^x}-\frac{(1-\exp (-\mathrm{\μ}))(a_3+a_4{w}^y)}{\mathrm{\μ}}\\ \mathrm{\β}\′=1+a_5+a_{6}1\mathrm{n\μ}+a_7{w}^z\\ r=b_0+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{b}_i{\mathrm{\ρ}}_i\end{array}$

其中，E₀为大气层顶的太阳辐照度；θ_s为从MODIS几何定位产品得到的太阳天顶角；D为日地距离；r为大气层顶处宽波段反照率；μ为太阳天顶角的余弦函数；a1～a7、x、y和z均为代表各种地表类型的常数；w为从MODIS大气可降水产品得到的可降水数据；b0～b7为依赖于观测天顶角和太阳天顶角的系数；ρ_i为从MODIS定标辐射产品得到的大气层顶窄波段的反射率；

步骤3.2：利用该研究区域MODIS数据产品中热红外波段的辐射亮度产品和地表温度/比辐射率产品估算地表长波净辐射R_lw，计算公式为

L_d＝c₀+c₁×M₂₉+c₂×M₃₄+c₃×M₃₃+c₄×M₃₆+c₅×M₂₈+c₆×M₃₁

R_lw＝ε_sL_d-5.67×10^-8ε_sT_s ⁴

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε₃₂ ²

其中，L_d为大气下行长波辐射；c_i（i=0～6）为依赖于观测天顶角和地表高程的系数；M为从MODIS定标辐射产品提取的热红外波段辐射亮度数据，下标中的数字表示MODIS传感器热红外波段；ε_s为地表比辐射率；ε₃₁和ε₃₂为从MODIS地表温度/比辐射率产品提取到的MODIS31和32波段的地表比辐射率；

步骤3.3：利用估算得到的地表短波净辐射和地表长波净辐射计算全遥感的地表净辐射R_n，计算公式为

R_n＝R_sw+R_lw

其中，R_sw为地表短波净辐射，R_lw为地表长波净辐射；

步骤3.4：利用地表净辐射R_n和植被覆盖度F_r估算部分植被覆盖的土壤热通量G，计算公式为

Γ＝Γ_v+(1-F_r)(Γ_s-Γ_v)

$\mathrm{\Γ}=\frac{G}{R_n}$

G＝[Γ_v+(1-F_r)(Γ_s-Γ_v)]R_n

其中，Γ为部分植被覆盖情况下土壤热通量与地表净辐射的比值，Γ_v和Γ_s分别为全植被覆盖和干燥裸土对应的土壤热通量与地表净辐射的比值，Γ_v和Γ_s为已知量。

本发明解决上述技术问题的技术方案如下：一种实现基于MODIS数据的地表蒸散发全遥感反演方法的系统，包括数据预处理模块、地表蒸发比反演模块、地表可利用能量反演模块和地表蒸散发反演模块；

所述数据预处理模块，其用于对获取的研究区域相关MODIS数据产品进行预处理；

所述地表蒸发比反演模块，其用于利用经预处理的MODIS数据产品构建空间变化信息散点图，实现像元尺度地表蒸发比EF的全遥感反演；

所述地表可利用能量反演模块，其用于利用经预处理的MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；

所述地表蒸散发反演模块，其用于根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述地表蒸发比反演模块包括地表植被覆盖度计算单元、空间变化信息散点图构建单元、极端地表温度计算单元、蒸发比指示因子计算单元和蒸发比计算单元；

所述地表植被覆盖度计算单元，其用于利用该研究区域的MODIS数据产品中的归一化植被指数NDVI估算该研究区域的地表植被覆盖度F_r，计算公式为

$F_r={\left(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^2$

其中，NDVI_min和NDVI_max分别为对应裸土和全植被覆盖的归一化植被指数；

空间变化信息散点图构建单元，其用于利用该研究区域的MODIS数据产品中的地表温度和估算得到的植被覆盖度F_r，构建空间变化信息散点图，其以植被覆盖度F_r为横坐标，以地表温度为纵坐标；

所述极端地表温度计算单元，其用于计算空间变化信息散点图内地表表层土壤水分分别为0和1时不同植被覆盖度像元所对应的地表温度T_max,i和T_min,i，

T_max,_i＝T_smax+F_r(T_smin-T_smax)

T_min,_i＝T_smin

其中，T_smax为在植被覆盖度和地表表层相对土壤水分均为0时的地表温度，T_smin指在植被覆盖度为1时的地表温度；

蒸发比指示因子计算单元，其用于利用线性内插法估算空间变化信息散点图中不同土壤水分和植被覆盖度条件下每个像元的蒸发比指示因子计算公式为，

其中，T_max,i和T_min,i分别指在植被覆盖度为i时，地表表层相对土壤水分为0和1时的地表温度，T_si指在植被覆盖度为i时某一像元所对应的地表温度，

指在植被覆盖度为i时，地表表层相对土壤水分为1时的蒸发比指示因子，

指在植被覆盖度为i时，地表表层相对土壤水分为0时的蒸发比指示因子；

所述蒸发比计算单元，其用于利用蒸发比指示因子，实现像元尺度蒸发比EF的全遥感参数化

其中，

为蒸发比指示因子，γ为湿度计常数，为0.066kpa/℃，Δ为空间变化信息散点图内每个像元的温度-饱和水汽压曲线斜率，计算公式为

$\mathrm{\Δ}=4098\×\frac{0.6108\×\exp \left[\frac{17.25\×(T-273.15)}{(T-273.15)+237.3}\right]}{{[(T-273.15)+237.3]}^2}$

其中，T为地表温度。

进一步，所述地表可利用能量反演模块包括地表短波净辐射计算单元、地表长波净辐射计算单元、地表净辐射计算单元和土壤热通量计算单元；

所述地表短波净辐射计算单元，其用于利用该研究区域MODIS数据产品中的太阳天顶角产品、大气可降水产品和大气顶部窄波段反射率产品估算地表短波净辐射R_sw，计算公式为，

$R_{\mathrm{sw}}=\frac{E_0\cos {\mathrm{\θ}}_s}{D^2}(\mathrm{\α}\′-\mathrm{\β}\′r)$

其中，

$\begin{array}{c}\mathrm{\α}\′=1-\frac{a_1}{\mathrm{\μ}}-\frac{a_2}{{\mathrm{\μ}}^x}-\frac{(1-\exp (-\mathrm{\μ}))(a_3+a_4{w}^y)}{\mathrm{\μ}}\\ \mathrm{\β}\′=1+a_5+a_{6}1\mathrm{n\μ}+a_7{w}^z\\ r=b_0+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{b}_i{\mathrm{\ρ}}_i\end{array}$

其中，E₀为大气层顶太阳辐照度；θ_s为从MODIS几何定位产品得到的太阳天顶角；D为日地距离；r为大气层顶处宽波段反照率；μ为太阳天顶角的余弦函数；a1～a7、x、y和z均为代表各种地表类型的常数；w为从MODIS大气可降水产品得到的可降水数据；b0～b7为依赖于观测天顶角和太阳天顶角的系数；ρ_i为从MODIS定标辐射产品得到的大气层顶窄波段的反射率；

所述地表长波净辐射计算单元，其用于利用该研究区域MODIS数据产品中热红外波段的辐射亮度产品和地表温度/比辐射率产品估算地表长波净辐射R_lw，计算公式为

L_d＝c₀+c₁×M₂₉+c₂×M₃₄+c₃×M₃₃+c₄×M₃₆+c₅×M₂₈+c₆×M₃₁

R_lw＝ε_sL_d-5.67×10^-8ε_sT_s ⁴

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε₃₂ ²

其中，L_d为下行长波辐射；c_i（i=0～6）为依赖于观测天顶角和地表高程的系数；M为从MODIS定标辐射产品提取的热红外波段辐射亮度数据，下标中的数字表示MODIS传感器热红外波段；ε_s为地表比辐射率；ε₃₁和ε₃₂为从MODIS地表温度/比辐射率产品提取到的MODIS31和32波段的地表比辐射率；

所述地表净辐射计算单元，其用于利用估算得到的地表短波净辐射和地表长波净辐射计算全遥感的地表净辐射R_n，计算公式为

R_n＝R_sw+R_lw

其中，R_sw为地表短波净辐射，R_lw为地表长波净辐射；

所述土壤热通量计算单元，其用于利用地表净辐射R_n和植被覆盖度F_r估算土壤热通量G，计算公式为

Γ＝Γ_v+(1-F_r)(Γ_s-Γ_v)

$\mathrm{\Γ}=\frac{G}{R_n}$

G＝[Γ_v+(1-F_r)(Γ_s-Γ_v)]R_n

其中，Γ为部分植被覆盖情况下土壤热通量与地表净辐射的比值，Γ_v和Γ_s分别为全植被覆盖和干燥裸土对应的土壤热通量与地表净辐射的比值，Γ_v和Γ_s为已知量。

附图说明

图1为本发明所述一种基于MODIS数据的地表蒸散发全遥感反演方法流程图；

图2为本发明所述步骤1的具体实现流程图；

图3为本发明所述计算蒸散发比EF的流程图；

图4为本发明所述构建的空间变化信息散点图的示意图；

图5为本发明所述计算地表净辐射R_n和土壤热通量G的流程图；

图6为本发明所述计算地表蒸散发LE的流程图；

图7为本发明所述一种基于MODIS数据的地表蒸散发全遥感反演系统框图；

图8为本发明所述所述地表蒸发比反演模块结构框图；

图9为本发明所述地表可利用能量反演模块结构框图。

附图中，各标号所代表的部件列表如下：

1、数据预处理模块，2、地表蒸发比反演模块，3、地表可利用能量反演模块，4、地表蒸散发反演模块,201、地表植被覆盖度计算单元，202、空间变化信息散点图构建单元，203、极端地表温度计算单元，204、蒸发比指示因子计算单元，205、蒸发比计算单元，301、地表短波净辐射计算单元，302、地表长波净辐射计算单元，303、地表净辐射计算单元，304、土壤热通量计算单元。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

如图1所示，一种基于MODIS数据的地表蒸散发全遥感反演方法，包括如下步骤：

步骤1：确定研究区域，获取研究区域相关MODIS数据产品并进行预处理；

步骤2：利用经预处理的MODIS数据产品构建空间变化信息散点图，并进行像元尺度的地表蒸发比EF的全遥感反演；

步骤3：利用经预处理的MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；

步骤4：根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度的地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF。

其中，步骤1中所述相关MODIS数据产品包括地表温度/比辐射率产品、归一化植被指数产品、定标辐射产品、几何定位产品和大气可降水产品。

在LP DAAC（Land Processes Distributed Active Archive Center，http://lpdaac.usgs.gov/）下载获取MODIS地表温度/比辐射率产品（MOD11_L2或MYD11_L2）和归一化植被指数产品NDVI（MOD13A2或MYD13A2）产品，在LAADS（Level1and Atmosphere Archive and Distribution System，http://ladsweb.nascom.nasa.gov）下载定标辐射产品（MOD021KM或MYD021KM）、几何定位产品（MOD03或MYD03）和大气可降水产品（MOD05或MYD05）。在获取MODIS数据产品后，需要利用Modis Reprojection Tool和MODIS Swath Reprojection Tool对数据进行一些预处理以便能够实现对研究区内所有数据的投影转换、剪切处理，并需要将经过处理后的MODIS产品数值转换为实际的物理量。由于云会显著影响卫星反演得到的地表温度，因此，需要剔除研究区内受云影响的像元。最后，需要利用数字高程产品，根据设定的高程上下限，剔除研究区内高程与像元平均高程相差较大的像元。

如图2所示，步骤1中对研究区域内相关MODIS数据产品进行预处理的过程如下：

步骤1.1：在表层相对土壤水分和植被覆盖度全范围变化的条件下，确定研究区域；

步骤1.2：对该研究区域相关MODIS数据产品进行投影转换和剪切处理，进而将处理后的MODIS数据产品数值转换为实际的物理量；

步骤1.3：利用地表温度/比辐射率产品中的质量控制文件剔除研究区内受云污染的像元数据；

步骤1.4：利用几何定位产品中的数字高程数据，设定高程阈值，剔除受地形冷却效应影响的像元数据。

如图3所示，步骤2中利用经预处理的数据产品构建空间变化信息散点图，进行像元尺度的地表蒸发比的全遥感反演的过程为：

步骤2.1：利用该研究区域的MODIS数据产品中的归一化植被指数NDVI估算该研究区域每个像元的地表植被覆盖度F_r，计算公式为

$F_r={\left(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^2$

其中，NDVI为每个像元的归一化植被指数，NDVI_min和NDVI_max分别为对应裸土和全植被覆盖的归一化植被指数，本实施例中NDVI_min=0.2，NDVI_max=0.86；

步骤2.2：利用该研究区域的MODIS数据产品中每个像元的地表温度T和估算得到的植被覆盖度F_r，构建以植被覆盖度F_r为横坐标，以地表温度T为纵坐标的空间变化信息散点图；

步骤2.3：通过鲁棒的线性回归方程回归空间变化信息散点图中地表表层相对土壤水分分别为0和1时不同植被覆盖度像元的地表温度T_max,i和T_min,i，

T_max,i＝T_smax+F_r(T_smin-T_smax)

T_min,i＝T_smin

其中，T_smax为在植被覆盖度和地表表层相对土壤水分均为0时的地表温度，T_smin指在植被覆盖度为1时的地表温度；

步骤2.4：利用线性内插法估算空间变化信息散点图中不同土壤水分和植被覆盖度条件下每个像元的蒸发比指示因子计算公式为，

其中，T_max,i和T_min,i分别指在植被覆盖度为i时，地表表层相对土壤水分分别为0和1时的地表温度，T_si指在植被覆盖度为i时某一像元所对应的地表温度，

指在植被覆盖度为i时，地表表层相对土壤水分为1时的蒸发比指示因子，

指在植被覆盖度为i时，地表表层相对土壤水分为0时的蒸发比指示因子；

步骤2.5：利用地表温度T计算空间变化信息散点图内每个像元的温度-饱和水汽压曲线斜率Δ，计算公式为

$\mathrm{\Δ}=4098\×\frac{0.6108\×\exp \left[\frac{17.25\×(T-273.15)}{(T-273.15)+237.3}\right]}{{[(T-273.15)+237.3]}^2}$

步骤2.6：利用温度-饱和水汽压曲线斜率Δ和蒸发比指示因子计算像元尺度的蒸发比EF，计算公式为

其中，γ为湿度计常数，为0.066kpa/℃。

其中，步骤2.3中确定空间变化信息散点图内地表表层土壤水分分别为0和1时，不同植被覆盖度条件所对应的地表温度（即干边和湿边）。

在由步骤2.3构建的二维散点图空间中，由于步骤1.3和1.4中剔除了一些受云污染和受地形冷却效应影响的像元，这就造成了一些植被覆盖度条件下对应的最高和最低地表温度并非极端表层土壤水分条件下的温度，即形成了包括真实干点像元和伪干点像元的干边和包括真实湿点像元和伪湿点像元的湿边，如图4所示，地表表层相对土壤水分为0的是干边，地表表层相对土壤水分为1的是湿边。这些所谓的伪干点和伪湿点会显著影响一般的线性回归方法结果。在拟合不同植被覆盖度条件下对应的最高地表温度的时候，需要利用鲁棒性强的线性回归方法消除这些“伪”干点的影响。

在利用步骤2.3确定干湿限之后，需按如下方式开展空间变化信息散点图内每个像元的蒸发比指示因子

参数的估算。首先，设置空间变化信息散点图内干燥裸土（F_r=0和干边上的最大地表温度处）的值为0（即在位置（F_r=0,T_smax）处=0），全植被覆盖处和湿边的

为1.26（即在位置（F_r=1,T_smin）处

=1.26）。其次，利用线性内插法估算空间变化信息散点图内像元i（F_r，T_si）的

值：

1）假定干边上的

值随着植被覆盖度F_r在

和

之间线性变化以确定干边上的值（即

）；由于湿边上的

一般波动很小，所以假定湿边上的值为一常数，即

进行湿边

的确定；

2）在给定的F_r件下，通过假定值随着地表温度的降低而在

和之间线性增加，确定像元i的

值

对于空间变化信息散点图内任一像元(F_r，T_si)的

值：

其中，

T_max,i＝T_smax+F_r(T_smin-T_smax)

T_min,i＝T_smin

如图5所示，步骤3中根据相关MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演的过程为：

步骤3.1：利用该研究区域MODIS数据产品中的太阳天顶角产品、大气可降水产品和大气顶部窄波段反射率产品估算地表短波净辐射R_sw，计算公式为，

$R_{\mathrm{sw}}=\frac{E_0\cos {\mathrm{\θ}}_s}{D^2}(\mathrm{\α}\′-\mathrm{\β}\′r)$

其中，

$\begin{array}{c}\mathrm{\α}\′=1-\frac{a_1}{\mathrm{\μ}}-\frac{a_2}{{\mathrm{\μ}}^x}-\frac{(1-\exp (-\mathrm{\μ}))(a_3+a_4{w}^y)}{\mathrm{\μ}}\\ \mathrm{\β}\′=1+a_5+a_{6}1\mathrm{n\μ}+a_7{w}^z\\ r=b_0+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{b}_i{\mathrm{\ρ}}_i\end{array}$

其中，E₀为大气层顶的太阳辐照度；θ_s为从MODIS几何定位产品得到的太阳天顶角；D为日地距离；r为大气层顶处宽波段反照率；μ为太阳天顶角的余弦函数；a1～a7、x、y和z均为代表各种地表类型的常数（陆地、海洋、冰/雪等）；w为从MODIS大气可降水产品得到的可降水数据；b0～b7为依赖于观测天顶角和太阳天顶角的系数；ρ_i为从MODIS定标辐射产品得到的大气层顶窄波段（波段1-7）的反射率；

步骤3.2：利用该研究区域MODIS数据产品中热红外波段的辐射亮度产品和地表温度/比辐射率产品估算地表长波净辐射R_lw，计算公式为

L_d＝c₀+c₁×M₂₉+c₂×M₃₄+c₃×M₃₃+c₄×M₃₆+c₅×M₂₈+c₆×M₃₁

R_lw＝ε_sL_d-5.67×10^-8ε_sT_s ⁴

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε₃₂ ²

其中，L_d为大气下行长波辐射；c_i（i=0～6）为依赖于观测天顶角和地表高程的系数；M为从MODIS定标辐射产品提取的热红外波段辐射亮度数据，本实施例中分别为波段28、29、31、33、34和36，下标中的数字表示MODIS传感器热红外波段；ε_s为地表比辐射率；ε₃₁和ε₃₂为从MODIS地表温度/比辐射率产品提取到的MODIS31和32波段的地表比辐射率；

步骤3.3：利用估算得到的地表短波净辐射和地表长波净辐射计算全遥感的地表净辐射R_n，计算公式为

R_n＝R_sw+R_lw

其中，R_sw为地表短波净辐射，R_lw为地表长波净辐射；

步骤3.4：利用地表净辐射R_n和植被覆盖度F_r估算部分植被覆盖的土壤热通量G，计算公式为

Γ＝Γ_v+(1-F_r)(Γ_s-Γ_v)

$\mathrm{\Γ}=\frac{G}{R_n}$

G＝[Γ_v+(1-F_r)(Γ_s-Γ_v)]R_n

其中，Γ为部分植被覆盖情况下土壤热通量与地表净辐射的比值，Γ_v和Γ_s分别为全植被覆盖和干燥裸土对应的土壤热通量与地表净辐射的比值，Γ_v和Γ_s为已知量。

如图6所示，为本发明所述计算地表蒸散发LE的流程图，根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度的地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF

如图7所示，一种实现基于MODIS数据的地表蒸散发全遥感反演方法的系统，包括数据预处理模块1、地表蒸发比反演模块2、地表可利用能量反演模块3和地表蒸散发反演模块4；

所述数据预处理模块1，其用于对获取的研究区域相关MODIS数据产品进行预处理；

所述地表蒸发比反演模块2，其用于利用经预处理的MODIS数据产品构建空间变化信息散点图，实现像元尺度地表蒸发比EF的全遥感反演；

所述地表可利用能量反演模块3，其用于利用经预处理的MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；

所述地表蒸散发反演模块4，其用于根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF。

如图8所示，所述地表蒸发比反演模块包括地表植被覆盖度计算单元201、空间变化信息散点图构建单元202、极端地表温度计算单元203、蒸发比指示因子计算单元204和蒸发比计算单元205；

所述地表植被覆盖度计算单元201，其用于利用该研究区域的MODIS数据产品中的归一化植被指数NDVI估算该研究区域的地表植被覆盖度F_r，计算公式为

$F_r={\left(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^2$

其中，NDVI_min和NDVI_max分别为对应裸土和全植被覆盖的归一化植被指数；

空间变化信息散点图构建单元202，其用于利用该研究区域的MODIS数据产品中的地表温度和估算得到的植被覆盖度F_r，构建空间变化信息散点图，其以植被覆盖度F_r为横坐标，以地表温度为纵坐标；

所述极端地表温度计算单元203，其用于计算空间变化信息散点图内地表表层土壤水分分别为0和1时不同植被覆盖度像元所对应的地表温度T_max,i和T_min,i，

T_max,i＝T_smax+F_r(T_smin-T_smax)

T_min,i＝T_smin

其中，T_smax为在植被覆盖度和地表表层相对土壤水分均为0时的地表温度，T_smin指在植被覆盖度为1时的地表温度；

蒸发比指示因子计算单元204，其用于利用线性内插法估算空间变化信息散点图中不同土壤水分和植被覆盖度条件下每个像元的蒸发比指示因子

计算公式为，

其中，T_max,i和T_min,i分别指在植被覆盖度为i时，地表表层相对土壤水分为0和1时的地表温度，T_si指在植被覆盖度为i时某一像元所对应的地表温度，

指在植被覆盖度为i时，地表表层相对土壤水分为1时的蒸发比指示因子，

指在植被覆盖度为i时，地表表层相对土壤水分为0时的蒸发比指示因子；

所述蒸发比计算单元205，其用于利用蒸发比指示因子，实现像元尺度蒸发比EF的全遥感参数化

其中，

为蒸发比指示因子，γ为湿度计常数，为0.066kpa/℃，Δ为空间变化信息散点图内每个像元的温度-饱和水汽压曲线斜率，计算公式为

$\mathrm{\Δ}=4098\×\frac{0.6108\×\exp \left[\frac{17.25\×(T-273.15)}{(T-273.15)+237.3}\right]}{{[(T-273.15)+237.3]}^2}$

其中，T为地表温度。

如图9所示，所述地表可利用能量反演模块包括地表短波净辐射计算单元301、地表长波净辐射计算单元302、地表净辐射计算单元303和土壤热通量计算单元304；

所述地表短波净辐射计算单元301，其用于利用该研究区域MODIS数据产品中的太阳天顶角产品、大气可降水产品和大气顶部窄波段反射率产品估算地表短波净辐射R_sw，计算公式为，

$R_{\mathrm{sw}}=\frac{E_0\cos {\mathrm{\θ}}_s}{D^2}(\mathrm{\α}\′-\mathrm{\β}\′r)$

其中，

$\begin{array}{c}\mathrm{\α}\′=1-\frac{a_1}{\mathrm{\μ}}-\frac{a_2}{{\mathrm{\μ}}^x}-\frac{(1-\exp (-\mathrm{\μ}))(a_3+a_4{w}^y)}{\mathrm{\μ}}\\ \mathrm{\β}\′=1+a_5+a_{6}1\mathrm{n\μ}+a_7{w}^z\\ r=b_0+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{b}_i{\mathrm{\ρ}}_i\end{array}$

其中，E₀为大气层顶太阳辐照度；θ_s为从MODIS几何定位产品得到的太阳天顶角；D为日地距离；r为大气层顶处宽波段反照率；μ为太阳天顶角的余弦函数；a1～a7、x、y和z均为代表各种地表类型的常数；w为从MODIS大气可降水产品得到的可降水数据；b0～b7为依赖于观测天顶角和太阳天顶角的系数；ρ_i为从MODIS定标辐射产品得到的大气层顶窄波段的反射率；

所述地表长波净辐射计算单元302，其用于利用该研究区域MODIS数据产品中热红外波段的辐射亮度产品和地表温度/比辐射率产品估算地表长波净辐射R_lw，计算公式为

L_d＝c₀+c₁×M₂₉+c₂×M₃₄+c₃×M₃₃+c₄×M₃₆+c₅×M₂₈+c₆×M₃₁

R_lw＝ε_sL_d-5.67×10^-8ε_sT_s ⁴

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε₃₂ ²

其中，L_d为下行长波辐射；c_i（i=0～6）为依赖于观测天顶角和地表高程的系数；M为从MODIS定标辐射产品提取的热红外波段辐射亮度数据，下标中的数字表示MODIS传感器热红外波段；ε_s为地表比辐射率；ε₃₁和ε₃₂为从MODIS地表温度/比辐射率产品提取到的MODIS31和32波段的地表比辐射率；

所述地表净辐射计算单元303，其用于利用估算得到的地表短波净辐射和地表长波净辐射计算全遥感的地表净辐射R_n，计算公式为

R_n＝R_sw+R_lw

其中，R_sw为地表短波净辐射，R_lw为地表长波净辐射；

所述土壤热通量计算单元304，其用于利用地表净辐射R_n和植被覆盖度F_r估算土壤热通量G，计算公式为

Γ＝Γ_v+(1-F_r)(Γ_s-Γ_v)

$\mathrm{\Γ}=\frac{G}{R_n}$

G＝[Γ_v+(1-F_r)(Γ_s-Γ_v)]R_n

其中，Γ为部分植被覆盖情况下土壤热通量与地表净辐射的比值，Γ_v和Γ_s分别为全植被覆盖和干燥裸土对应的土壤热通量与地表净辐射的比值，Γ_v和Γ_s为已知量。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于MODIS数据的地表蒸散发全遥感反演方法，其特征在于，包括如下步骤：

步骤1：确定研究区域，获取研究区域相关MODIS数据产品并进行预处理；

步骤2：利用经预处理的MODIS数据产品构建空间变化信息散点图，并进行像元尺度的地表蒸发比EF的全遥感反演；

步骤3：利用经预处理的数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；

步骤4：根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度的地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF。

2.根据权利要求1所述一种基于MODIS数据的地表蒸散发全遥感反演方法，其特征在于，步骤1中所述相关MODIS数据产品包括地表温度/比辐射率产品、归一化植被指数产品、定标辐射产品、几何定位产品和大气可降水产品。

3.根据权利要求1所述一种基于MODIS数据的地表蒸散发全遥感反演方法，其特征在于，步骤1中对研究区域内相关MODIS数据产品进行预处理的过程如下：

步骤1.1：在表层相对土壤水分和植被覆盖度全范围变化的条件下，确定研究区域；

步骤1.2：对该研究区域相关MODIS数据产品进行投影转换和剪切处理，进而将处理后的MODIS数据产品数值转换为实际的物理量；

步骤1.3：利用地表温度/比辐射率产品中的质量控制文件剔除研究区内受云污染的像元数据；

步骤1.4：利用几何定位产品中的数字高程数据，设定高程阈值，剔除受地形冷却效应影响的像元数据。

4.根据权利要求2所述一种基于MODIS数据的地表蒸散发全遥感反演方法，其特征在于，步骤2中利用经预处理的MODIS数据产品构建空间变化信息散点图，进行像元尺度的地表蒸发比的全遥感反演的过程为：

步骤2.1：利用该研究区域的MODIS数据产品中的归一化植被指数NDVI估算该研究区域每个像元的地表植被覆盖度F_r，计算公式为

$F_r={\left(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^2$

其中，NDVI为每个像元的归一化植被指数，NDVI_min和NDVI_max分别为对应裸土和全植被覆盖的归一化植被指数；

步骤2.2：利用该研究区域的MODIS数据产品中每个像元的地表温度T和估算得到的植被覆盖度F_r，构建以植被覆盖度F_r为横坐标，以地表温度T为纵坐标的空间变化信息散点图；

步骤2.3：通过鲁棒的线性回归方程回归空间变化信息散点图中地表表层相对土壤水分为0时不同植被覆盖度像元的地表温度T_max,_i；同时，设定地表表层相对土壤水分为1时不同植被覆盖度的像元具有相同的地表温度T_min,_i，

T_max,_i＝T_smax+F_r(T_smin-T_smax)

T_min,i＝T_smin

其中，T_smax为在植被覆盖度和地表表层相对土壤水分均为0时的地表温度，T_smin指在植被覆盖度为1时的地表温度；

步骤2.4：利用线性内插法估算空间变化信息散点图中不同土壤水分和植被覆盖度条件下每个像元的蒸发比指示因子

计算公式为，

其中，T_max,i和T_min,i分别指在植被覆盖度为i时，地表表层相对土壤水分为0和1时的地表温度，T_si指在植被覆盖度为i时某一像元所对应的地表温度，

指在植被覆盖度为i时，地表表层相对土壤水分为1时的蒸发比指示因子，

指在植被覆盖度为i时，地表表层相对土壤水分为0时的蒸发比指示因子；

步骤2.5：利用地表温度T计算空间变化信息散点图内每个像元的温度-饱和水汽压曲线斜率Δ，计算公式为

$\mathrm{\Δ}=4098\×\frac{0.6108\×\exp \left[\frac{17.25\×(T-273.15)}{(T-273.15)+237.3}\right]}{{[(T-273.15)+237.3]}^2}$

步骤2.6：利用温度-饱和水汽压曲线斜率Δ和蒸发比指示因子

计算像元尺度的蒸发比EF，计算公式为

其中，γ为湿度计常数，为0.066kpa/℃。

5.根据权利要求1所述一种基于MODIS数据的地表蒸散发全遥感反演方法，其特征在于，步骤3中根据相关MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演的过程为：

步骤3.1：利用该研究区域MODIS数据产品中的太阳天顶角产品、大气可降水产品和大气顶部窄波段反射率产品估算地表短波净辐射R_sw，计算公式为，

$R_{\mathrm{sw}}=\frac{E_0\cos {\mathrm{\θ}}_s}{D^2}(\mathrm{\α}\′-\mathrm{\β}\′r)$

其中，

$\begin{array}{c}\mathrm{\α}\′=1-\frac{a_1}{\mathrm{\μ}}-\frac{a_2}{{\mathrm{\μ}}^x}-\frac{(1-\exp (-\mathrm{\μ}))(a_3+a_4{w}^y)}{\mathrm{\μ}}\\ \mathrm{\β}\′=1+a_5+a_{6}1\mathrm{n\μ}+a_7{w}^z\\ r=b_0+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{b}_i{\mathrm{\ρ}}_i\end{array}$

其中，E₀为大气层顶的太阳辐照度；θ_s为从MODIS几何定位产品得到的太阳天顶角；D为日地距离；r为大气层顶处宽波段反照率；μ为太阳天顶角的余弦函数；a1～a7、x、y和z均为代表各种地表类型的常数；w为从MODIS大气可降水产品得到的可降水数据；b0～b7为依赖于观测天顶角和太阳天顶角的系数；ρ_i为从MODIS定标辐射产品得到的大气层顶窄波段的反射率；

步骤3.2：利用该研究区域MODIS数据产品中热红外波段的辐射亮度产品和地表温度/比辐射率产品估算地表长波净辐射R_lw，计算公式为

L_d＝c₀+c₁×M₂₉+c₂×M₃₄+c₃×M₃₃+c₄×M₃₆+c₅×M₂₈+c₆×M₃₁

R_lw＝ε_sL_d-5.67×10^-8ε_sT_s ⁴

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε₃₂ ²

其中，L_d为大气下行长波辐射；c_i（i=0～6）为依赖于观测天顶角和地表高程的系数；M为从MODIS定标辐射产品提取的热红外波段辐射亮度数据，下标中的数字表示MODIS传感器热红外波段；ε_s为地表比辐射率；ε₃₁和ε₃₂为从MODIS地表温度/比辐射率产品提取到的MODIS31和32波段的地表比辐射率；

步骤3.3：利用估算得到的地表短波净辐射和地表长波净辐射计算全遥感的地表净辐射R_n，计算公式为

R_n＝R_sw+R_lw

其中，R_sw为地表短波净辐射，R_lw为地表长波净辐射；

步骤3.4：利用地表净辐射R_n和植被覆盖度F_r估算部分植被覆盖的土壤热通量G，计算公式为

Γ＝Γ_v+(1-F_r)(Γ_s-Γ_v)

$\mathrm{\Γ}=\frac{G}{R_n}$

G＝[Γ_v+(1-F_r)(Γ_s-Γ_v)]R_n

其中，Γ为部分植被覆盖情况下土壤热通量与地表净辐射的比值，Γ_v和Γ_s分别为全植被覆盖和干燥裸土对应的土壤热通量与地表净辐射的比值，Γ_v和Γ_s为已知量。

6.一种实现权利要求1-5任一项基于MODIS数据的地表蒸散发全遥感反演方法的系统，其特征在于，包括数据预处理模块、地表蒸发比反演模块、地表可利用能量反演模块和地表蒸散发反演模块；

所述数据预处理模块，其用于对获取的研究区域相关MODIS数据产品进行预处理；

所述地表蒸发比反演模块，其用于利用经预处理的MODIS数据产品构建空间变化信息散点图，实现像元尺度地表蒸发比EF的全遥感反演；

所述地表可利用能量反演模块，其用于利用经预处理的MODIS数据产品进行像元尺度的地表净辐射R_n与土壤热通量G的全遥感反演；

所述地表蒸散发反演模块，其用于根据地表蒸发比EF、地表净辐射R_n和土壤热通量G，进行像元尺度地表蒸散发LE的全遥感反演，计算公式为LE＝(R_n-G)*EF。

7.根据权利要求6所述一种基于MODIS数据的地表蒸散发全遥感反演系统，其特征在于，所述地表蒸发比反演模块包括地表植被覆盖度计算单元、空间变化信息散点图构建单元、极端地表温度计算单元、蒸发比指示因子计算单元和蒸发比计算单元；

所述地表植被覆盖度计算单元，其用于利用该研究区域的MODIS数据产品中的归一化植被指数NDVI估算该研究区域的地表植被覆盖度F_r，计算公式为

$F_r={\left(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}\right)}^2$

其中，NDVI_min和NDVI_max分别为对应裸土和全植被覆盖的归一化植被指数；

空间变化信息散点图构建单元，其用于利用该研究区域的MODIS数据产品中的地表温度和估算得到的植被覆盖度F_r，构建空间变化信息散点图，其以植被覆盖度F_r为横坐标，以地表温度为纵坐标；

所述极端地表温度计算单元，其用于计算空间变化信息散点图内地表表层土壤水分分别为0和1时不同植被覆盖度像元所对应的地表温度T_max,i和T_min,i，

T_max,i＝T_smax+F_r(T_smin-T_smax)

T_min,i＝T_smin

其中，T_smax为在植被覆盖度和地表表层相对土壤水分均为0时的地表温度，T_smin指在植被覆盖度为1时的地表温度；

蒸发比指示因子计算单元，其用于利用线性内插法估算空间变化信息散点图中不同土壤水分和植被覆盖度条件下每个像元的蒸发比指示因子

计算公式为，

其中，T_max,i和T_min,i分别指在植被覆盖度为i时，地表表层相对土壤水分为0和1时的地表温度，T_si指在植被覆盖度为i时某一像元所对应的地表温度，

指在植被覆盖度为i时，地表表层相对土壤水分为1时的蒸发比指示因子，

指在植被覆盖度为i时，地表表层相对土壤水分为0时的蒸发比指示因子；

所述蒸发比计算单元，其用于利用蒸发比指示因子，实现像元尺度蒸发比EF的全遥感参数化

其中，为蒸发比指示因子，γ为湿度计常数，为0.066kpa/℃，Δ为空间变化信息散点图内每个像元的温度-饱和水汽压曲线斜率，计算公式为

$\mathrm{\Δ}=4098\×\frac{0.6108\×\exp \left[\frac{17.25\×(T-273.15)}{(T-273.15)+237.3}\right]}{{[(T-273.15)+237.3]}^2}$

其中，T为地表温度。

8.根据权利要求6所述一种基于MODIS数据的地表蒸散发全遥感反演系统，其特征在于，所述地表可利用能量反演模块包括地表短波净辐射计算单元、地表长波净辐射计算单元、地表净辐射计算单元和土壤热通量计算单元；

所述地表短波净辐射计算单元，其用于利用该研究区域MODIS数据产品中的太阳天顶角产品、大气可降水产品和大气顶部窄波段反射率产品估算地表短波净辐射R_sw，计算公式为，

$R_{\mathrm{sw}}=\frac{E_0\cos {\mathrm{\θ}}_s}{D^2}(\mathrm{\α}\′-\mathrm{\β}\′r)$

其中

$\begin{array}{c}\mathrm{\α}\′=1-\frac{a_1}{\mathrm{\μ}}-\frac{a_2}{{\mathrm{\μ}}^x}-\frac{(1-\exp (-\mathrm{\μ}))(a_3+a_4{w}^y)}{\mathrm{\μ}}\\ \mathrm{\β}\′=1+a_5+a_{6}1\mathrm{n\μ}+a_7{w}^z\\ r=b_0+\underset{i=1}{\overset{7}{\mathrm{\Σ}}}{b}_i{\mathrm{\ρ}}_i\end{array}$

其中，E₀为大气层顶太阳辐照度；θ_s为从MODIS几何定位产品得到的太阳天顶角；D为日地距离；r为大气层顶处宽波段反照率；μ为太阳天顶角的余弦函数；a1～a7、x、y和z均为代表各种地表类型的常数；w为从MODIS大气可降水产品得到的可降水数据；b0～b7为依赖于观测天顶角和太阳天顶角的系数；ρ_i为从MODIS定标辐射产品得到的大气层顶窄波段的反射率；

所述地表长波净辐射计算单元，其用于利用该研究区域MODIS数据产品中热红外波段的辐射亮度产品和地表温度/比辐射率产品估算地表长波净辐射R_lw，计算公式为

L_d＝c₀+c₁×M₂₉+c₂×M₃₄+c₃×M₃₃+c₄×M₃₆+c₅×M₂₈+c₆×M₃₁

R_lw＝ε_sL_d-5.67×10^-8ε_sT_s ⁴

ε_s＝0.273+1.778ε₃₁-1.807ε₃₁ε₃₂-1.037ε₃₂+1.774ε₃₂ ²

其中，L_d为下行长波辐射；c_i（i=0～6）为依赖于观测天顶角和地表高程的系数；M为从MODIS定标辐射产品提取的热红外波段辐射亮度数据，下标中的数字表示MODIS传感器热红外波段；ε_s为地表比辐射率；ε₃₁和ε₃₂为从MODIS地表温度/比辐射率产品提取到的MODIS31和32波段的地表比辐射率；

所述地表净辐射计算单元，其用于利用估算得到的地表短波净辐射和地表长波净辐射计算全遥感的地表净辐射R_n，计算公式为

R_n＝R_sw+R_lw

其中，R_sw为地表短波净辐射，R_lw为地表长波净辐射；

所述土壤热通量计算单元，其用于利用地表净辐射R_n和植被覆盖度F_r估算土壤热通量G，计算公式为

Γ＝Γ_v+(1-F_r)(Γ_s-Γ_v)

$\mathrm{\Γ}=\frac{G}{R_n}$

G＝[Γ_v+(1-F_r)(Γ_s-Γ_v)]R_n

其中，Γ为部分植被覆盖情况下土壤热通量与地表净辐射的比值，Γ_v和Γ_s分别为全植被覆盖和干燥裸土对应的土壤热通量与地表净辐射的比值，Γ_v和Γ_s为已知量。

Images