CN105303040A - 时间连续的地表蒸散发数据的计算方法 - Google Patents

Abstract

本发明实施例提供了一种时间连续的地表蒸散发数据的计算方法。该方法主要包括：在非晴好日，使用SCE_UA算法利用代价函数计算出能量分配项的调整因子α和植被冠层阻抗项的调整因子β；根据所述能量分配项的调整因子α和植被冠层阻抗项的调整因子β优化非晴好日的地表蒸散发模拟值的计算方法，再利用优化后的计算方法计算出所述非晴好日的地表蒸散发数据。本发明致力于解决遥感估算地表蒸散发中的时间尺度扩展问题，结合同化方法和Penman-Monteith公式对日蒸散发进行日到月/年的时间尺度扩展，生产可供多领域使用的高精度时间连续区域蒸散发数据集，降低由于非晴好日观测数据带来的观测误差，促进遥感估算ET产品在水文循环、陆面能量相关领域的广泛应用。

Description

时间连续的地表蒸散发数据的计算方法

技术领域

本发明涉及地表水技术领域，尤其涉及一种时间连续的地表蒸散发数据的计算方法。

背景技术

蒸散发(Evapotranspiration，ET)是地表水量平衡的重要分量，其能量形式——潜热通量(LE)则是地表能量平衡的重要分量。潜热通量将下垫面贮存的热量通过湍流交换的方式传输给大气边界层中靠近下垫面的气层，是地球系统中重要的能量、物质交换过程。液态水的相变所需热量即为潜热，正是大气热量最大的单一热源，因此，对于潜热通量的研究在天气、气候动力学中具有重要的地位。地气间动量、水分、热量的交换过程，特别是潜热、显热通量的定量化对于气候变化模式、大气环流模式以及水文模式中边界条件的确定意义重大。

传统的地表蒸散发观测仪器观测到的蒸散发均限于点尺度或者网格尺度，数据缺乏区域代表性或难以反映地表蒸散发的空间分布。然而，蒸散发的估算和许多部门的现实问题息息相关，如南水北调等大型调水工程的设计、三峡水库的后效估计、科学用水合理灌溉方案的制定、农业旱灾的监测和预警、农作物产量的估算及预报等。在这些蒸散发的实际应用中，均需要区域尺度、长时间序列的蒸散量。由于仪器观测地表蒸散发通常选在平坦、均一的下垫面，很难代表所有类型的下垫面特征，因此仅通过内插法很难获得可靠的、能够反映下垫面异质性的区域地表蒸散发。遥感技术通过大面积测定地表辐射以及温度状况，可以提供与地表蒸散发变化相关的特征参量，这使得通过遥感估算蒸散发成为可能。目前，国内外学者已经建立了大量基于遥感数据或者与遥感数据相结合的区域蒸散发估算方法，通过与遥感像元尺度的地面观测蒸散发进行验证，认为在非均一下垫面上，遥感监测地表蒸散发具有优越性。

然而，由于陆面过程的复杂性，遥感估算蒸散发的精度受到了多种不确定性的影响，其中包括模型结构、模型输入数据、参数化方案以及空间和时间尺度方面的不确定性。过去的研究重点多集中于瞬时地表蒸散发的获取，对于瞬时蒸散发的时间尺度扩展方法尚未给予足够的重视。遥感瞬时蒸散发的时间扩展即将通过反演得到的瞬时蒸散发在时间尺度上进行扩展，得到小时、日、旬、月乃至年的蒸散发。一方面，遥感获得的仅是卫星过境时刻的瞬时地表信息，需要通过瞬时到日的时间扩展方法得到晴好日蒸散发；另一方面，在计算长时间序列的地表蒸散发时，由于受到云雨天气的影响，造成了大量数据的缺失，需要对晴好日的蒸散发进行时间尺度扩展，才能得到月尺度的蒸散发。

目前，对于日到月/年的时间尺度扩展方法的研究包括数学方法、地表阻抗方法和参考作物蒸发比方法等，然而这些方法均未考虑遥感观测自身的观测误差，因此具有一定局限性。

发明内容

本发明的实施例提供了一种时间连续的地表蒸散发数据的计算方法，以生成可供多领域使用的高精度时间连续区域蒸散发数据集。

为了实现上述目的，本发明采取了如下技术方案。

一种时间连续的地表蒸散发数据的计算方法，包括：

在非晴好日，使用SCE_UA算法利用代价函数计算出能量分配项的调整因子α和植被冠层阻抗项的调整因子β；

根据所述能量分配项的调整因子α和植被冠层阻抗项的调整因子β优化非晴好日的地表蒸散发模拟值的计算方法，再利用优化后的计算方法计算出所述非晴好日的地表蒸散发数据。

进一步地，所述的方法还包括：

晴好日的地表蒸散发数据用Penman-Monteith公式计算，所述Penman-Monteith公式的基本形式为：

$\mathrm{\λ}ET=\frac{\mathrm{\Δ}(R_n-G)+{\mathrm{\ρ}}_a{c}_p\frac{(e_s-e_a)}{r_a}}{\mathrm{\Δ}+\mathrm{\γ}(1+\frac{r_s}{r_a})}$

其中，ET是晴好日的地表蒸散发数据，λ是汽化潜热，Δ是饱和水汽压的平均斜率，R_n是净辐射，G是土壤热通量，ρ_a是空气密度，c_p是干空气比热容，r_s是冠层表面阻抗，r_a是空气动力学阻抗，e_a是实际水汽压，e_s是饱和水汽压，γ是干湿球温度计常数：

$\mathrm{\γ}=\frac{c_{p}P}{\mathrm{\ϵ}\mathrm{\λ}}=0.665\×{10}^{-3}P$

ε是水蒸气/干空气的分子比重，P是大气压，λ是汽化潜热

进一步地，所述的方法还包括：将晴好日的地表蒸散发遥感观测数据扩展到月/年，得到时间连续的地表蒸散发遥感观测数据。

进一步地，所述的在非晴好日，使用SCE_UA算法利用代价函数计算出能量分配项的调整因子α和植被冠层阻抗项的调整因子β，包括：

设定非晴好日的代价函数如下：

$J\left(x\right)=\sqrt{\frac{1}{2N}\underset{i=1}{\overset{N}{\Σ}}\[{(ET-{\mathrm{ET}}_{sim})}^2/Q+{(ET-{\mathrm{ET}}_{obs})}^2/R\]}$

J(x)是代价函数，ET是优化后的分析值，ET_obs是晴好日观测的地表蒸散发，ET_sim是晴好日的地表蒸散发，N是同化窗口内的遥感观测值个数，Q是模型误差，R是观测误差；

所述SCE_UA算法的步骤如下：

1、采样：在参数空间内随机采集s个样本点，并计算每个样本点的代价函数值；

2、排列样本点：按照每个样本点的函数值的大小对s个样本点进行升序排序，使第一个样本点的代价函数值最小，最后一个样本点的代价函数值最大；

3、拆分为复合体：将所有s个样本点按照排序平均分为p个复合体A¹…A^p，每个复合体包括m个样本点，其中，A^k＝{x^k _j/x^k _j＝X_{k+m(k- 1)}，j＝1，…，m}，k＝1，…，p；

4、复合体进化：根据竞争复合体进化算法进化每个复合体；

所述竞争复合体进化算法包括：

41、从复合体中随机选择q个点组成一个子复合体；

42、根据样本点的函数值，从子复合体中选择一个函数值最小的样本点，并计算剩余样本点的质心；

43、将所述函数值最小的样本点以所述剩余样本点的质心为中心投影，找到一个新的样本点，如果新样本点的函数值在合理范围内，则到第44步，否则，在合理范围内随机取一个样本点，到第45步；

44、如果所述新的样本点的函数值小于所述函数值最小的样本点的函数值，则用新样本点取代所述函数值最小的样本点，并到第47步，否则到第45步；

45、将计算质心和最差点之间的中点称为缩小点，计算所述缩小点的函数值，如果缩小点的函数值比所述函数值最小的样本点的函数值要小，则用缩小点取代所述函数值最小的样本点，并到第7步，否则到第46步；

46、在合理范围内随机取一个点，将所述函数值最小的样本点用这个点取代；

47、重复第42步-第46步r次，r是一个用户指定的参数；

48、重复第1步-第7步s次，s是在复合体混合前需要进化的次数；

5、混合复合体：将所有复合体中的样本点合为一个样本总体，将这个总体中的所有样本按照升序排序，重新分配样本点到p个复合体中，并重复3至5的过程；

6、检查迭代结果：如果迭代结果满足设定的迭代标准，则停止迭代，获取迭代停止时的α和β。

进一步地，所述的根据所述述能量分配项的调整因子α和植被冠层阻抗项的调整因子β优化非晴好日的地表蒸散发模拟值的计算方法，再利用优化后的计算方法计算出所述非晴好日的地表蒸散发数据，包括：

$\mathrm{\λ}ET=\frac{\mathrm{\Δ}\mathrm{\α}(R_n-G)+{\mathrm{\ρ}}_a{c}_p\frac{(e_s-e_a)}{r_a}}{\mathrm{\Δ}+\mathrm{\γ}(1+\frac{{\mathrm{\βr}}_s}{r_a})}$

其中，α是所述迭代停止时的能量分配项的调整因子，β是所述迭代停止时的植被冠层阻抗项的调整因子，ET是晴好日的地表蒸散发数据，ρ_a是空气密度，c_p是干空气比热容，r_s是冠层表面阻抗，r_a是空气动力学阻抗，e_a是实际水汽压，e_s是饱和水汽压，γ是干湿球温度计常数。

进一步地，在水体下垫面，水面蒸散发的变化用基于Penman公式的下式表示：

$ET=\frac{(\mathrm{\α}\frac{\mathrm{\Δ}}{\mathrm{\γ}}Rn+{\mathrm{\βE}}_a)}{\frac{\mathrm{\Δ}}{\mathrm{\γ}}+1}.$

由上述本发明的实施例提供的技术方案可以看出，本发明实施例致力于解决遥感估算地表蒸散发中的时间尺度扩展问题，结合同化方法和Penman-Monteith公式对日蒸散发进行日到月/年的时间尺度扩展，生产可供多领域使用的高精度时间连续区域蒸散发数据集，降低由于非晴好日观测数据带来的观测误差，促进遥感估算ET产品在水文循环、陆面能量相关领域的广泛应用。

本发明附加的方面和优点将在下面的描述中部分给出，这些将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

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

图1为本发明实施例一提供的一种时间连续的地表蒸散发数据的计算方法的处理流程图；

图2为2009年北京地表蒸散发的空间分布图(左：同化算法时间扩展结果；右：ET_MOD16，空白处为城市、水体或湿地下垫面)；

图3为北京市土地覆盖图。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的任一单元和全部组合。

本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

为便于对本发明实施例的理解，下面将结合附图以几个具体实施例为例做进一步的解释说明，且各个实施例并不构成对本发明实施例的限定。

实施例一

数据同化是指在考虑数据时空分布以及观测场和背景场误差的基础上，在数值模型的运行过程中融合新的观测数据的方法。本算法正是致力于解决遥感估算地表蒸散发中的时间尺度扩展问题，结合同化方法和Penman-Monteith公式对日蒸散发进行日到月/年的时间尺度扩展，生成1km分辨率的时间连续区域蒸散发数据集，促进遥感估算ET产品在水文循环、陆面能量相关领域的广泛应用。

本发明实施例构建了一个完整的遥感估算地表的蒸散发数据同化系统，主要包括：模型算子、观测算子、同化策略及数据同化算法等，选取Penman-Monteith模型作为同化系统的模型算子，主要用于模拟地表蒸散发的连续变化过程。选择SCE_UA数据同化算法，主要用于优化模型参数。采用同化晴好日遥感估算蒸散发数据优化非晴好日调整因子的同化策略。

本发明实施例提供的一种时间连续的地表蒸散发数据的计算方法的处理流程如图1所示，包括如下的处理步骤：

步骤S110、使用Penman-Monteith公式来计算晴好日地表蒸散发模拟值。

本发明实施例设计的具体算法内容如下：

1.Penman-Monteith公式

Penman-Monteith公式是同化系统中的模型算子，也是该算法的核心内容，具有坚实的物理基础，是目前公认的适应性强、计算精度高且可靠的植被蒸散发估算方法。公式中所涉及的大多数参数可以通过常规气象观测得到，在空气动力学阻抗和冠层表面阻抗已知的情况下，可用来直接计算作物蒸散发。Penman-Monteith公式的基本形式为：

$\mathrm{\λ}ET=\frac{\mathrm{\Δ}(R_n-G)+{\mathrm{\ρ}}_a{c}_p\frac{(e_s-e_a)}{r_a}}{\mathrm{\Δ}+\mathrm{\γ}(1+\frac{r_s}{r_a})}$

其中，ET是日蒸散发(mm/day)，Δ是饱和水汽压的平均斜率，R_n是净辐射(W/m²)，G是土壤热通量(W/m²)，λ是汽化潜热。

ρa是空气密度(kg/m³)，c_p是干空气比热容(C_p＝1.013×10-3MJ/kg·℃)，r_s是冠层表面阻抗(s/m)，r_a是空气动力学阻抗(s/m)，e_a是实际水汽压(kPa)，e_s是饱和水汽压(kPa)，γ是干湿球温度计常数：

$\mathrm{\γ}=\frac{c_{p}P}{\mathrm{\ϵ}\mathrm{\λ}}=0.665\×{10}^{-3}P$

其中，ε是水蒸气/干空气的分子比重(ε＝0.622)，λ是汽化潜热(λ＝2.45MJ/kg)，汽化潜热指在温度不变时，单位质量的水在汽化过程中所吸收的热量，λ是温度的函数，在高温时λ较小，低温时较大，但是在常温下λ的变化非常微小，为了简化计算过程，使用20℃时的汽化潜热进行计算。P是大气压(kPa)，如没有气压的观测数据，则

$P=101.3{\left(\frac{293-0.0065h}{293}\right)}^{5.26}$

其中，h是站点高程(m)。

Penman-Monteith公式是现有公式，本发明实例并没有对其作出改进，而是挑选了合适的参数化方案对公式中的各参数进行计算。

根据世界气象组织建议的饱和水汽压公式(戈夫-格雷奇公式)，在-49.9～100℃范围内：

$\begin{array}{c}{\mathrm{lge}}_s\left(T_a\right)=10.79574(1-T_{00}/(T_a+273.16\left)\right)-5.02800\lg \left(\right(T_a+273.16)/T_{00})\\ +1.50475\×{10}^{-4}\[1-{10}^{-8.2969\left(\right(T_a+273.16)/T_{00}-1)}\]+0.42873\×{10}^{-3}\[{10}^{4.76955(1-T_{00}/(T_a+273.16))}-1\]\\ +0.78614\end{array}$

其中，T₀₀＝273.16K，是水的三相点温度，T_a是大气温度(℃)，e_s(T_a)是在空气温度为T_a时的饱和水汽压(kPa)。Δ是饱和水汽压的平均斜率(kPa/℃)：

$\mathrm{\Δ}=\frac{4098\[0.6108e_s\left(T_a\right)\]}{{(T_a+237.3)}^2}$

其中，T_a是空气温度(℃)。

在计算日蒸散发时，假定为中性层结，日空气动力学阻抗r_a的定义为：

$r_a=\frac{\[\ln \left(\frac{z-d}{z_{0m}}\right)\]\[\ln \left(\frac{z-d}{z_{0h}}\right)\]}{k^{2}u}$

其中，z是风速观测高度(m)，d是零平面位移(m)，z_0m是空气动力学地表粗糙度(m)，von-Karman常数k＝0.41，u是风速(m/s)。z_0m是地表粗糙度，是对数风廓线公式中平均风速等于零时的高度。z_0m也是动力学粗糙度，是地表与大气之间进行动量交换的高度。与之类似，热力学粗糙度z_0h则是温度梯度外延到空气温度等于地表温度时的高度，由于近地面层中的热量交换主要依靠分子扩散的作用，而动量交换还包括气压梯度的作用，因此热力学粗糙度Z_0h要小于动力学粗糙度Z_0m。两者间的这一差异使用热量传输附加阻尼kB^-1＝ln(z_0m/z_0h)表达。kB^-1的影响因子包括气象条件、植被结构和下垫面状况等，在日、月、年的时间尺度上均具有一定差异，因此需要针对下垫面类型采用不同的kB^-1参数化方案。

(1)在有植被覆盖的下垫面：

${\mathrm{kB}}^{-1}=\frac{{\mathrm{kC}}_d}{4C_t\frac{u_{*}}{u}(1-e^{-n_{ec}}/2)}{f}_c^2+2f_c{f}_s\frac{k\frac{u_{*}{z}_{0m}}{uz}}{C_t^{*}}+{\mathrm{kB}}_s^{-1}{f}_s^2$

其中，C_d是叶片的拖曳系数，C_d＝0.2；C_t是叶片的热传输系数，在这里取C_t＝0.01；C_t ^*是土壤的热传输系数：

$C_t^{*}={\mathrm{Pr}}^{-2/3}{\left({\mathrm{Re}}^{*}\right)}^{-1/2}$

在这里，Pr是Prandtl数，取值为0.7，Re^*是土壤表面的粗糙雷诺数：

Re_*＝z_0su_*/ν

其中，v为动量分子粘滞系数：

ν＝1.327×10^-5(P₀/P)(T_a/T₀₀)^1.81

其中，P₀＝101.3kPa。

u^*是摩擦风速(m/s)：

u^*＝ku/ln((z-d)/z_0s)

其中，z是风速的观测高度(m)，u是风速(m/s)。

nec是冠层内风速廓线削弱系数：

$n_{ec}=\frac{C_{d}LAI}{2u_{*}^2/u^2}$

f_c是地表的植被覆盖度，是绿叶的垂直投影面积与土壤总面积之比，LAI是叶面积指数。

f_s＝1-f_c，是土壤覆盖度。当f_c＝0时相当于的裸土kB^-1参数化方案：

${\mathrm{kB}}_s^{-1}=2.46{\left({\mathrm{Re}}_{*}\right)}^{1/4}-\ln 7.4$

(2)居民地或交通用地等城市下垫面属于粗糙刚性下垫面，使用Brutsaert的参数化方案：

kB^-1＝2.46(Re^*)^0.25-2

在水体下垫面，采用Brutsaert的水面kB^-1参数化方案：

$z_{oh}=0.169\exp (-1.53u_{*}^{1/4})$

Re^*是粗糙雷诺数：

Re^*＝u^*z_0m/ν

其中，u^*是摩擦风速(m/s)：

u^*＝ku/ln((z-d)/z_0m)

其中，z是风速的观测高度(m)，u是风速(m/s)。

在遥感区域应用中，采用Liu-2007模型中的z_0m参数化方案估算植被覆盖区域的z_0m：

$z_{0m}=\exp (a\frac{NDVI}{\mathrm{\α}}+b)$

其中，NDVI是归一化植被指数、α是地表宽波段反照率，a和b是通过观测资料拟合的参数，a＝0.26，b＝-2.21，再根据z_0m推出零平面位移d：

d＝4.9z_0m

在其他下垫面上，z_0m取值方案如表1所示。

表1典型下垫面z_0m取值方案

Penman-Monteith公式中的冠层表面阻抗r_s采用Landsberg&Gower的参数化方案，该方案也被用于生产MODIS全球地表蒸散发产品：

C_s＝c_Lm(T_min)m(VPD)

C_c＝C_sLAI

r_s＝1/C_c

其中，c_L是单位叶面积的平均潜在气孔导度(m/s)，m(T_min)是日最低温度限制因子，m(VPD)是水汽压差VPD(VPD＝e_s-e_a)的限制因子，LAI是叶面积指数，通过对LAI积分将气孔导度升尺度到冠层导度。以上的Penman-Monteith公式参数化方案主要用于植被以及城市下垫面，在水体下垫面，则使用计算一定气象条件下开阔水面蒸发量的Penman公式作为模型算子：

${\mathrm{ET}}_0=\frac{(\frac{\mathrm{\Δ}}{\mathrm{\γ}}Rn+E_a)}{\frac{\mathrm{\Δ}}{\mathrm{\γ}}+1}$

其中，E_a是水面附近的空气干燥力(mm)：

E_a＝0.35(0.5+u₂/100)(e_s-e_a)

其中，u₂是2米处的风速(miles/day，1mile＝1609m)。

步骤120、在非晴好日，使用SCE_UA算法利用代价函数计算出能量分配项的调整因子α和植被冠层阻抗项的调整因子β。

利用数据同化方法对晴好日的地表蒸散发数据进行时间尺度扩展，该时间尺度扩展是将晴好日的地表蒸散发数据扩展到月/年，从而得到了全年完整的蒸散发时间序列。然后对全年的时间序列进行月平均或年平均，就可以得到月/年的地表蒸散发。

本发明实施例采用的数据同化策略在能量平衡项Rn-G和植被阻抗项r_s/r_a前分别加入调整因子α和β，通过优化这两个因子来改进Penman-Monteith公式的预报。

数据同化方法是根据“最优”统计的方法将各种已知信息(模型预报和观测)融合在一起，从而得到能够最优地描述给定分辨率下的地表状态的一个过程。在本系统的同化策略中，进行同化的目的是通过调整参数尽量减小模拟值和观测值之间的差异，因此，在非晴好日，可以建立以下的代价函数：

$J\left(x\right)=\sqrt{\frac{1}{2N}\underset{i=1}{\overset{N}{\Σ}}\[{(ET-{\mathrm{ET}}_{sim})}^2/Q+{(ET-{\mathrm{ET}}_{obs})}^2/R\]}$

其中，J(x)是代价函数，ET是优化后的分析值，ET_obs是晴好日观测的地表蒸散发，ET_sim是晴好日Penman-Monteith公式模拟的地表蒸散发，N是同化窗口内的遥感观测值个数。Q是模型误差，R是观测误差，属于先验知识。

同化窗口的大小主要根据经验进行选择，即需要用户提前进行设置，如下垫面为农田，同化窗口不应大于农作物的物候期。在北京市的示例中，同化窗口为8天，窗口内的遥感观测值个数根据天气状况各不相同。

在代价函数中，没有明显的时间参数，是因为每一个非晴好日的时候会使用SCE_UA算法计算代价函数达到最小值时的α和β，而代价函数只和同化窗口内的晴好日观测值有关(即N)。

通过对晴好日卫星过境时刻瞬时蒸散发的瞬时到日时间尺度扩展，可得到晴好日的遥感估算蒸散发数据集。在非晴好日，同化系统则自动调用SCE_UA算法，通过迭代搜索同化窗口内代价函数J(x)得最小值时的最优解，最终得到同化策略中的模型参数α和β，再代入Penman-Monteith公式得到优化后的非晴好日地表蒸散发，从而对遥感观测的地表蒸散发数据集进行插补，生产时间连续的遥感估算地表蒸散发数据集。

本发明中采用的SCE_UA算法是一种全局最优策略，对参数初值的设定并不敏感。SCE_UA的算法步骤如下：

1、采样：在参数空间内随机采集s个样本点，并计算每个样本点的代价函数值。在缺少先验知识的情况下，认为样本点符合均匀分布；

2、排列样本点：按照函数值的大小对s个样本点进行升序排序，使第一个点的函数值最小，最后一个点的函数值最大；

3、拆分为复合体：将所有s个样本点按照排序平均分为p个复合体A¹…A^p，每个复合体包括m个样本点，其中，

A^k＝(x^k _j/x^k _j＝X_k+m(k-1)，j＝1，…，m}，k＝1，…，p；

4、复合体进化：根据竞争复合体进化算法(theCompetitiveComplexEvolutionalgorithm，CCE)进化每个复合体；

5、混合复合体：将所有复合体中的样本点合为一个样本总体，将这个样本总体中的所有样本按照升序排序，重新分配样本点到p个复合体中，并重复(3)至(5)的过程；

6、检查迭代结果：如果迭代结果满足设定的迭代标准，则停止迭代，获取迭代停止时的α和β。

SCE_UA算法的关键组成部分是第四步中的CCE算法，主要步骤如下：

41、从复合体中随机选择q个点组成一个子复合体；

42、根据样本点的函数值，从子复合体中选择一个函数值最小的样本点，并计算剩余样本点的质心；

43、将所述函数值最小的样本点以所述剩余样本点的质心为中心投影，找到一个新的样本点，事先根据已有的观测数据和先验知识得到样本点的合理范围，如果新样本点的函数值在合理范围内，则到第44步，否则，在合理范围内随机取一个样本点，到第45步；

44、如果所述新的样本点的函数值小于所述函数值最小的样本点的函数值，则用新样本点取代所述函数值最小的样本点，并到第47步，否则到第45步；

45、将计算质心和最差点之间的中点称为缩小点，计算所述缩小点的函数值，如果缩小点的函数值比所述函数值最小的样本点的函数值要小，则用缩小点取代所述函数值最小的样本点，并到第7步，否则到第46步；

46、在合理范围内随机取一个点，将所述函数值最小的样本点用这个点取代；

47、重复第42步-第46步r次，r是一个用户指定的参数；

48、重复第1步-第7步s次，s是在复合体混合前需要进化的次数；

步骤130、根据所述能量分配项的调整因子α和植被冠层阻抗项的调整因子β优化晴好日的地表蒸散发数据的计算方法，得到所述非晴好日的地表蒸散发数据。

我们假设在非晴好日，地表蒸散发的变化可以用基于Penman-Monteith公式的下式表示：

$\mathrm{\λ}ET=\frac{\mathrm{\Δ}\mathrm{\α}(R_n-G)+{\mathrm{\ρ}}_a{c}_p\frac{(e_s-e_a)}{r_a}}{\mathrm{\Δ}+\mathrm{\γ}(1+\frac{{\mathrm{\βr}}_s}{r_a})}$

其中，α是所述迭代停止时的能量分配项的调整因子，β是所述迭代停止时的植被冠层阻抗项的调整因子，ET是晴好日的地表蒸散发数据，ρ_a是空气密度，c_p是干空气比热容，r_s是冠层表面阻抗，r_a是空气动力学阻抗，e_a是实际水汽压，e_s是饱和水汽压，γ是干湿球温度计常数。将迭代停止时的调整因子带入上式后，即可改变非晴好日的地表蒸散发估算值。

在水体下垫面，水面蒸散发的变化用基于Penman公式的下式表示，同样将迭代停止时的调整因子带入，得到优化后的非晴好日的水面蒸散发：

$ET=\frac{(\mathrm{\α}\frac{\mathrm{\Δ}}{\mathrm{\γ}}Rn+{\mathrm{\βE}}_a)}{\frac{\mathrm{\Δ}}{\mathrm{\γ}}+1}$

数据同化方法是根据“最优”统计的方法将各种已知信息(模型预报和观测)融合在一起，从而得到能够最优地描述给定分辨率下的地表状态的一个过程。不同的数据同化方法使用不同的标准来定义所谓的“最优”。在同化方法中，优化后的地表蒸散发的分析增量X^a-X⁰是观测增量ET_obs-ET_sim的加权平均，不同同化方法则是在不同的数学假定条件下通过算法推导这个后验权重。

实施例二：

北京是中国的首都，也是全国的政治、经济、文化中心，然而，随着城市用水矛盾日益加深，水资源已成为制约城市繁荣与发展的瓶颈。本算法通过SCE_UA方法将Liu-2007模型反演的晴好日蒸散发扩展到月，得到年蒸散发的空间分布。从图2中的同化方法估算结果(ET_SCEUA)可以看出，北京市城八区的年蒸散发明显低于周围的远郊区县，这主要是由于城市下垫面多为水泥等人造表面，年蒸散发较小。而周边区县中，北部的怀柔、密云年蒸散发高于南部的大兴、房山区。这种地表蒸散发分布特征与北京市的城市规划有关，北京的南部大兴、房山区开发的房地产项目较多，在城八区住房接近饱和的状况下，不少买房者选择购买南部的房源。从土地利用图(图3)中可以看出，除了城八区基本以居民区为主，在北部，居民地主要分布于延庆、怀柔县城，而在南部，居民地在大兴、房山区全境皆有分布。在主城区以及远郊区县的县城，由于城市中降水基本通过径流汇入下水道，人造地表的蒸发量很小。而在以耕地、林地为主的远郊区县，地表植被覆盖度较高，灌溉也进一步增加了水分来源，因此蒸散发较大。北京北部的密云水库具有宽阔的水面，应是北京地区地表蒸散发的高值区。

与图2中的MOD16A3产品估算结果相比，同化方法的年蒸散发的空间分布趋势与MOD16A3产品基本一致。MOD16A3产品估算结果整体偏低，这可能与MOD16采用的气象数据分辨率为1°×1.25°，反照率数据分辨率为0.05°有关。

对同化系统的结果和MOD16产品在不同的土地利用类型上进行了统计，表2说明在同化算法时间尺度扩展结果中，年蒸散发较大的土地利用类型依次为水体、耕地、园地、牧草地和未利用地。MOD16产品不包括城市下垫面和水体下垫面的蒸散发，其他土地利用类型中，年蒸散发较大的土地利用类型依次为林地、牧草地、耕地、园地，由于城区以人造水泥地面为主，因此这一结果是合理的。水体的蒸发量较大，在时间尺度扩展结果中是蒸散发第一的土地利用类型。耕地也是蒸散发量较大的土地利用类型，部分是因为北京属于半干旱半湿润地区，2009年较为干旱，年降雨量比2008年少30％，比多年平均少23％，作物主要依靠灌溉，因此蒸散发量较大。从地表蒸散发的空间分布来看，平原地区耕地的蒸散发明显大于山区林地的蒸散发。北京地区的耕地基本位于山前的平原地区，年降雨量在450mm左右，而山区的年降雨量在350mm左右，而考虑到耕地受到灌溉的影响，耕地的蒸散发能力实际上高于林地。

表22009年北京各土地利用/覆盖类型的年蒸散发

表32009年北京区域平均年蒸散发

根据水量平衡公式可以得到北京地区2009年的区域平均蒸散发为387mm，同化算法的估算结果大于水量平衡计算的年蒸散发，区域平均年蒸散发与水量平衡计算的蒸散发相比，偏高约46mm。

综上所述，本发明实施例致力于解决遥感估算地表蒸散发中的时间尺度扩展问题，结合同化方法和Penman-Monteith公式对日蒸散发进行日到月/年的时间尺度扩展，生成1km分辨率的时间连续区域蒸散发数据集，促进遥感估算ET产品在水文循环、陆面能量相关领域的广泛应用。

本发明的有益效果为有效利用已有的晴好日遥感估算蒸散发，生产可供多领域使用的高精度时间连续区域蒸散发数据集，降低由于非晴好日观测数据带来的观测误差。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本发明可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置或系统实施例而言，由于其基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。以上所描述的装置及系统实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (6)

1.一种时间连续的地表蒸散发数据的计算方法，其特征在于，包括：

在非晴好日，使用SCE_UA算法利用代价函数计算出能量分配项的调整因子α和植被冠层阻抗项的调整因子β；

根据所述能量分配项的调整因子α和植被冠层阻抗项的调整因子β优化非晴好日的地表蒸散发模拟值的计算方法，再利用优化后的计算方法计算出所述非晴好日的地表蒸散发数据。

2.根据权利要求1所述的时间连续的地表蒸散发数据的计算方法，其特征在于，所述的方法还包括：

晴好日的地表蒸散发数据用Penman-Monteith公式计算，所述Penman-Monteith公式的基本形式为：

其中，ET是晴好日的地表蒸散发数据，λ是汽化潜热，Δ是饱和水汽压的平均斜率，R_n是净辐射，G是土壤热通量，ρ_a是空气密度，c_p是干空气比热容，r_s是冠层表面阻抗，r_a是空气动力学阻抗，e_a是实际水汽压，e_s是饱和水汽压，γ是干湿球温度计常数：

ε是水蒸气/干空气的分子比重，P是大气压，λ是汽化潜热。

3.根据权利要求2所述的时间连续的地表蒸散发数据的计算方法，其特征在于，所述的方法还包括：将晴好日的地表蒸散发遥感观测数据扩展到月/年，得到时间连续的地表蒸散发遥感观测数据。

4.根据权利要求1或2或3所述的时间连续的地表蒸散发数据的计算方法，其特征在于，所述的在非晴好日，使用SCE_UA算法利用代价函数计算出能量分配项的调整因子α和植被冠层阻抗项的调整因子β，包括：

设定非晴好日的代价函数如下：

J(x)是代价函数，ET是优化后的分析值，ET_obs是晴好日观测的地表蒸散发，ET_sim是晴好日的地表蒸散发，N是同化窗口内的遥感观测值个数，Q是模型误差，R是观测误差；

所述SCE_UA算法的步骤如下：

（1）、采样：在参数空间内随机采集s个样本点，并计算每个样本点的代价函数值；

（2）、排列样本点：按照每个样本点的函数值的大小对s个样本点进行升序排序，使第一个样本点的代价函数值最小，最后一个样本点的代价函数值最大；

（3）、拆分为复合体：将所有s个样本点按照排序平均分为p个复合体A¹…A^p，每个复合体包括m个样本点，其中，A^k＝{x^k _j/x^k _j＝x_k+m(k-1)，j＝1，…，m}，k＝1，…，p；

（4）、复合体进化：根据竞争复合体进化算法进化每个复合体；

所述竞争复合体进化算法包括：

（41）、从复合体中随机选择q个点组成一个子复合体；

（42）、根据样本点的函数值，从子复合体中选择一个函数值最小的样本点，并计算剩余样本点的质心；

（43）、将所述函数值最小的样本点以所述剩余样本点的质心为中心投影，找到一个新的样本点，如果新样本点的函数值在合理范围内，则到第44步，否则，在合理范围内随机取一个样本点，到第45步；

（44）、如果所述新的样本点的函数值小于所述函数值最小的样本点的函数值，则用新样本点取代所述函数值最小的样本点，并到第47步，否则到第45步；

（45）、将计算质心和最差点之间的中点称为缩小点，计算所述缩小点的函数值，如果缩小点的函数值比所述函数值最小的样本点的函数值要小，则用缩小点取代所述函数值最小的样本点，并到第7步，否则到第46步；

（46）、在合理范围内随机取一个点，将所述函数值最小的样本点用这个点取代；

（47）、重复第42步-第46步r次，r是一个用户指定的参数；

（48）、重复第1步-第7步s次，s是在复合体混合前需要进化的次数；

（5）、混合复合体：将所有复合体中的样本点合为一个样本总体，将这个总体中的所有样本按照升序排序，重新分配样本点到p个复合体中，并重复3至5的过程；

（6）、检查迭代结果：如果迭代结果满足设定的迭代标准，则停止迭代，获取迭代停止时的α和β。

5.根据权利要求4所述的时间连续的地表蒸散发数据的计算方法，其特征在于，所述的根据所述述能量分配项的调整因子α和植被冠层阻抗项的调整因子β优化非晴好日的地表蒸散发模拟值的计算方法，再利用优化后的计算方法计算出所述非晴好日的地表蒸散发数据，包括：

其中，a是所述迭代停止时的能量分配项的调整因子，β是所述迭代停止时的植被冠层阻抗项的调整因子，ET是晴好日的地表蒸散发数据，ρ_a是空气密度，c_p是干空气比热容，r_s是冠层表面阻抗，r_a是空气动力学阻抗，e_a是实际水汽压，e_s是饱和水汽压，γ是干湿球温度计常数。

6.根据权利要求5所述的时间连续的地表蒸散发数据的计算方法，所述的方法还包括：

在水体下垫面，水面蒸散发的变化用基于Penman公式的下式表示：