CN101551459A - 基于遥感的区域蒸散量监测方法 - Google Patents

Abstract

本发明提供了一种基于遥感的区域蒸散量监测方法，包括如下的步骤：(1)获取特定区域的多时相多分辨率卫星遥感影像数据；(2)对地面的关键参量进行遥感反演；(3)通过遥感地气交换模型计算地气能量交换过程中每个像元的阻抗；计算出用于加热每个像元上方空气的显热通量，通过能量平衡余项式计算潜热能量，推求有效表面阻抗；(4)运用时间重建技术从瞬时间断的遥感观测中获取完整的植被生长期过程曲线；(5)在逐日尺度上建立表面阻抗模型以获取逐日地表阻抗；(6)反复运用P－M模型求得逐日蒸散量，进行时间累计和尺度下推，从而最终确定区域地表蒸散量。通过在多个区域的实测验证，证明本区域蒸散量监测方法行之有效。

Description

基于遥感的区域蒸散量监测方法

技术领域

本发明涉及一种基于遥感的区域蒸散量监测方法，尤其涉及一种利用遥感图像提供的可见光和热红外信息，结合近地面的气象信息，进行地表能量平衡计算，并通过适当的时间扩展，最终确定区域地表蒸散量的方法，属于遥感应用技术领域。

背景技术

当前，全球水资源日益匮乏，为了合理利用和分配水资源，迫切需要深入了解不同植被覆盖和土地利用条件下的蒸发耗水情况。蒸散量包括植被截流蒸发量、植被蒸腾量、土壤蒸发量和水面蒸发量，是区域水量平衡和能量平衡的主要成分，不仅在水循环和能量循环过程中具有极其重要的作用，而且是生态过程与水文过程的重要纽带。因此，准确及时地获取特定区域的地表蒸散量数据，在农业、水文、森林和生态等领域都能够发挥十分重要的作用。

目前，蒸散量的实地测定方法仍以传统的水文气象学方法为主，主要包括直接测定下垫面(大气圈以地球的水陆表面为其下界，称为大气层的下垫面)的失水速率和大气中的水汽传送速率。总的来说，由于下垫面不均匀的影响，有关测点的测值对大范围内的地表蒸散发状况代表性不强，且布设测点成本高，难以形成实用的观测网络。

相对于传统的水文气象学方法，通过遥感数据监测地表蒸散的方法具有空间上连续和时间动态变化的特点。遥感数据的多光谱信息能够提供与地表能量平衡过程和地表覆盖状况密切相关的参数。用遥感手段进行区域尺度下非均匀下垫面的蒸散发估算，已经成为遥感应用领域的重要研究方向。

目前，在区域尺度下估算地表蒸散主要有两种方法：一是使用遥感表面辐射温度结合气温以及系列阻抗公式求取显热，通过能量平衡余项表示蒸散；二是在联合国粮食及农业组织(FAO)推荐的Penman-Monteith模型(简称为P-M模型)的基础上，用地温-植被指数或其他方式来推算P-M模型中的表面阻抗，直接进行蒸散计算。一般认为，前者精度较高而后者因其参数较为简单而易于应用。受云或其它大气因素的影响，遥感在非晴好天气时无法获得可用的可见光和热红外数据，所获得的陆面参数实际上是不连续的，因此要充分利用遥感晴好数据和气象观测数据。但是，通过遥感数据监测地表蒸散毕竟是一种间接的蒸发测量方法，涉及的环节较多、数据量大，而且还存在一些实际的技术问题需要解决，例如如何获得空间上一致时间上连续的遥感数据，如何使遥感数据提供计算方法和模型所需要的计算参数，如何通过对中间过程的监控以及结合地面数据进行模型的标定以减少不确定误差等。因此，现有的蒸散量遥感定量化监测技术仍需要进一步加以改进。

发明内容

本发明提供了一种基于遥感的区域蒸散量监测方法。该方法利用多时相的可见光/热红外遥感数据及每日常规气象数据建立了区域尺度上陆面蒸散量监测的数据处理链，从而实现了区域蒸散量的运行性遥感监测。

为实现上述的目的，本发明采用下述的技术方案：

一种基于遥感的区域蒸散量监测方法，其特征在于包括如下的步骤：

(1)获取特定区域的多时相多分辨率卫星遥感影像数据；

(2)对地面的关键参量进行遥感反演；

(3)通过遥感地气交换模型计算地气能量交换过程中每个像元的阻抗，并结合气象要素计算每个像元的地气温差，计算出用于加热每个像元上方空气的显热通量，通过能量平衡余项式计算潜热能量，并由Penman-Monteith模型推求逐日有效的表面阻抗；

(4)运用时间重建技术从瞬时间断的遥感观测中获取完整的植被生长期过程曲线；

(5)在逐日尺度上建立表面阻抗模型以获取逐日的地表阻抗；

(6)反复运用Penman-Monteith模型求得逐日蒸散量，由此进行时间累计和尺度下推，最终确定区域地表蒸散量。

其中，所述步骤(1)中，将ASTER/TM和AVHRR/MODIS作为卫星遥感影像数据源。

所述步骤(2)中，所述关键参量至少包括地表反照率、净辐射量、地表温度、植被覆盖度、叶面指数和动量粗糙度。

在计算所述净辐射量时，直射辐射的系数为0.45，散射辐射的系数为0.1。

在计算所述动量粗糙度时，分别对植被、地形起伏、非植被覆盖表面下的几何粗糙度进行加权以便获得所述动量粗糙度。

所述步骤(3)中，使用叶面指数进行植被与裸地的分类。

所述步骤(4)中，在获取全年晴好日的叶面指数的基础上，通过时间序列谐波分析谐函数分解方法得到逐日的叶面指数影像。

所述步骤(4)中，在植被生长期过程曲线中首先舍弃云噪声影响的异值点。

所述步骤(5)中，首先在预定时间段中选择一幅具有代表性的高分辨率蒸散量图像；如果没有具有代表性的高分辨率蒸散量图像，就从植被覆盖情况相似的邻近时间段中选取；如果无法确定具有代表性的高分辨率蒸散量图像，就选择前后相邻的图像，按照时间距离进行加权，形成该时间段的具有代表性的高分辨率蒸散量图像。

本区域蒸散量监测方法是一种基于多尺度、多分辨率的蒸散量遥感定量化监测方法。与现有技术相比，具有以下优点：

1)充分利用地面和高空气象数据，描述了边界层地气交换的特征参数。在应用地表能量平衡模型时，提高了模型的精度，从而提高了下垫面表面阻抗的准确性。

2)将微观尺度的SEBAL与宏观尺度的SEBS模型相结合，发挥了模型各自的特点，在遥感数据源的时空分辨率上取得了良好的平衡：

3)遥感模型常常因为天气状况无法获取清晰的图像而造成数据缺失。本方法强化了遥感应用中的时间尺度扩展方法，将遥感反演得到的表面阻抗参数和气象参数输入到P-M模型，将逐日的气象数据与间断的遥感数据结合起来，最终获得逐日连续的蒸散产品，提供的遥感蒸散数据能同时反映出研究区的蒸散发时空分布规律，具有重大的应用价值。

附图说明

下面结合附图和具体实施方式对本发明作进一步的说明。

图1为本发明所提供的基于遥感的区域蒸散量监测方法的基本流程示意图；

图2(a)为2003年5月18日的叶面指数LAI直方图，图2(b)为2003年5月18日的植被指数NDVI直方图，图2(c)为2003年8月14日的叶面指数LA1直方图，图2(d)为2003年8月14日的植被指数NDVI直方图；

图3(a)为2002年5月23日的叶面指数LAI直方图，图3(b)为2002年5月23日的植被指数NDVI直方图，图3(c)为2003年5月31日的叶面指数LAI直方图，图3(d)为2003年5月31日的植被指数NDVI直方图。

具体实施方式

一般而言，实用化的地表蒸散量估算方法需要进行如下的工作：对地表参数进行遥感反演；进行地表阻抗模型的选取、优化；对蒸散估算结果在空间和时间尺度上的扩展等。这些工作在本发明中都有自己的特色。下面结合图1，对本发明所提供的基于遥感的区域蒸散量监测方法展开详细的说明。

一.获取特定区域的多时相多分辨率卫星遥感影像数据

在本发明中，推荐使用ASTER/TM(先进空间热辐射反射辐射计/美国陆地卫星专题制图仪)和AVHRR/MODIS(改进型甚高分辨率辐射仪/中分辨率航天成像光谱仪)作为卫星遥感影像数据源。这主要是考虑到所使用的遥感信息源需同时具有可见光、近红外和热红外波段来提供植被和热辐射的信息。AVHRR/MODIS数据具有多光谱特征，时间分辨率为1天，具有两个热红外波段。ASTER/TM的时间分辨率为26天和16天，也具有热红外波段。

尽管提供AVHRR/MODIS数据的NOAA气象卫星每日都过境，但是受天气和云遮掩的影响，使得可以用于蒸散量估算的MODIS景数受到极大限制。在云量多且云层厚的条件下，整景影像的信息都为云层屏蔽，失去了使用价值。而在实际应用中能够使用的影像数据是云覆盖范围少于60％～70％的数据。遥感数据的缺失和云区对蒸散量的精度有很大的影响。晴好天气多，时相分布均匀，蒸散量估算结果接近于时段内实际蒸散量的变化过程；晴好天气小，时相分布不均，蒸散量估算结果会偏大或偏小于实际蒸散量的变化过程。如果流域内某个区域长期被厚云覆盖(如海拔较高的林区)，缺失的影像时段间隔越长，那么相对于该时段的蒸散量误差会比较大。

二.对地面的关键参量进行遥感反演

地面的关键参量包括地表反照率、净辐射量、地表温度、植被覆盖度、叶面指数、动量粗糙度等。下面分别对这些产量的遥感反演过程展开详细的说明。

1.地表反照率：

r₀＝0.2183r₁+0.1349r₂+0.2838r₃+0.2536r₄+0.0645r₅+0.0327r₆+0.0123r₇

式中r1～r7为可见光～热红外(0.45～12.5)7个窄谱段的反射系数，而r0则为宽谱段的地表反射率。

2.净辐射量：

净辐射量决定着可供显热、潜热分配的有效能量，是直接影响蒸散量大小的重要参量。在运用现有的SEBAL模型中的净辐射公式时，发现其日净辐射结果与实际观测结果有较大出入，特别是在多云日：

R_n24＝(1-α)R_a24τ_sw-110τ_sw    (1)

式(1)中，R_n24为全天净辐射，α为地表反照率，R_a24为全天地外太阳辐射，τ_sw为短波透过率。为反映日净辐射受云的影响，在透过率中引入日照时数：

${\mathrm{\τ}}_{\mathrm{sw}}=a_s+b_s\frac{n}{N}---\left(2\right)$

式(2)中，n为日照时数，N为晴空最大日照时数，a和b为经验系数。在现有技术中，通常分别使用系数0.5和0.25代表直射辐射和散射辐射的情况。

在本发明中，利用全国范围内的辐射台站2002～2003年的逐日数据分冬季(12～2月)和春夏秋季(3～11)分月对a和b进行了回归。在分析了辐射台站实测数据后，对于华北地区，将这两个系数分别更正到0.45和0.1。使用更正后系数计算的净辐射与台站测值相对误差有了明显减小，与原经验系数相比，全年逐日净辐射的相对误差小于20％的天数从60％增加到了77％。

3.地表温度

在本发明中，使用分裂窗算法计算地表温度。该分裂窗算法在Mao，K.，Z.，Qin，J.Shi和P.Gong发表的论文《A Pratical Split-window algorithmfor retrieving land-surface temperature from MODIS data》(载于International Journal of Remote Sensing.2005，26(5)，3181-3204)中有详细的记载，在此就不详细说明了。

由于地温输入对于蒸散发结果非常敏感，为获得遥感质量控制信息，本发明中进一步使用了MODIS云检测小组2002年发表的Cloud Mask算法来生成云波段。关于Cloud Mask算法，可以参见Steve A.，S.，Kathleen，M.Paul等在《NASA MODIS cloud algorithm》发表的报告《DiscriminatingClear-Sky From Cloud With Modis Algorithm Theoretical Basis Document(Mod35)》。需要说明的是，在本发明的计算过程中使用了MODIS的10个波段数据，设定多个特征项进行阈值判定，并将原有的4级置信度(晴空，可能晴空，可能有云，有云)中的后三类都归为云区，将云污染可能引入的影响降到最低。

4.植被覆盖度

对于植被覆盖度fc，通过如下的公式获得：

$\mathrm{fc}={\left(\frac{\mathrm{ndvi}-{\mathrm{ndvi}}_{\min }}{{\mathrm{ndvi}}_{\max }-{\mathrm{ndvi}}_{\min }}\right)}^{1/2}---\left(3\right)$

其中NDVI为植被指数，ndvi_min、ndvi_max为区域内NDVI的最大值与最小值，其实质是纯土壤与纯植被的NDVI。在实际应用中，由于NDVI的易饱和性和影像因云造成的缺失，很难取得理想的参考值，同时也使得不同单景之间的植被覆盖度很难具有可比性，在时相差异较大的情况下尤其如此，因此设全年的ndvi_min和ndvi_max分别为：ndvi_min＝0.05，ndvi_max＝0.85，使得计算出的fc具有可比性。

5.叶面指数

叶面指数LAI通过如下的公式进行计算：

$\mathrm{LAI}=-\frac{\ln \left(\frac{0.69-\mathrm{SAVI}}{0.59}\right)}{0.91}---\left(4\right)$

其中SAVI是土壤调节植被指数，SAVI＝1.1(R_nir-R_red)/(0.1+R_nir+R_red)。

6.动量粗糙度

动量粗糙度z_0m对于能量平衡余项法是非常重要的参量。原有模型一般在平坦地表条件下建立，因此仅考虑植被高度对粗糙度的影响，或者根据土地利用分类来指定经验值。这种做法在地形起伏条件下的适用性较差。因此本发明使用植被、地形起伏、非植被覆盖表面的几何粗糙度等因素来表达区域的综合动量粗糙度，具体方法如下：

植被区的几何粗糙度z_0m ^v，是植被指数NDVI的函数：

$z_{0m}^v=0.001+{(0.5\·\frac{\mathrm{NDVI}}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }})}^{25}---\left(5\right)$

其中NDVI＝(ρ_NIR-ρ_RED)/(ρ_NIR+ρ_RED)，是由近红波段和红波段反射率构建的植被指数，NDVI_max为区域中植被指数的最大值，NDVI_min为最小值。

在地形起伏的情况下，考虑坡度因子对几何粗糙度的影响：

$z_{0m}^T=z_{0m}^v\·(1+\frac{\mathrm{slope}-a}{b})---\left(6\right)$

其中slope为从数字地形数据中获取的坡度文件，a与b是经验系数。

非植被覆盖表面下的地表粗糙度z_0m ^r，使用从雷达数据中反演得到的表面几何粗糙度表征其大小：

$\log \left(z_{0m}^r\right)=-1.221+0.0906\left({\mathrm{\σ}}_0\right)---\left(7\right)$

其中σ0是从Radarsat卫星影像中获取的后向散射系数，以分贝形式表示。

区域的动量粗糙度是以上三者的加权，如下式：

$z_{0m}=w_1\·z_{0m}^v+w_2\·z_{0m}^T+w_3\·z_{0m}^r---\left(8\right)$

三.通过遥感地气交换模型计算地气能量交换过程中每个像元的阻抗；结合气象要素计算出每个像元的地气温差，从而计算出用于加热每个像元上方空气的显热通量，通过能量平衡余项式计算潜热能量，并由Penman-Monteith模型推求逐日有效的表面阻抗(即水汽进入大气的难易程度)。

目前，常用的遥感地气交换模型主要有SEBAL模型和SEBS模型。其中SEBAL是由荷兰Water-Watch公司的Bastinnassen开发的基于遥感的陆面能量平衡模型，用于估算非均匀陆面的能量平衡。该模型尤其适合于微观尺度，其特色在于利用图像中极端干湿像元的通量已知的假设，对其他像元的通量进行初始化。SEBAL模型作为结合了陆面能量平衡和湍流交换原理的图像处理方法，在TM/ETM级影像上取得了较好的应用结果，但不适于低分辨率影像的处理。一方面，作为重要参数的图像锚点在低分辨率图像上难于选取；另一方面，由于缺少成熟的时间扩展模块的支持，图像幅数越多，在时段计算时所作的调整和重算工作量就越大，年际蒸散的变化受到影像时相分布、影像质量等诸因素的影响。

SEBS是由荷兰瓦赫宁根大学的Su等人提出的用遥感数据计算大气湍流通量的模型，已经在欧洲和亚洲等地得到应用。SEBS模型发展了一套用于陆面参数遥感反演的参数化方法，从遥感数据中获得反照率、发射率、地表温度和植被覆盖度；一个综合的计算热量传输粗糙度长度的模型；一个使用迭代方法同时计算摩擦速度、显热通量和稳定度长度的模型，以及用于计算蒸发比的表面能量平衡指数。该模型比较适合宏观尺度，其特点在于定义了表面能量平衡指数的概念，通过假设的极度干燥和湿润的边界条件来计算实际的蒸发比。该模型还考虑了不同情况下大气稳定度订正的不同方法，但由于遥感信号本身的限制，仅将下垫面分成植被与土壤两种组分不足以代表城区、水体、山区等情况。

在本发明中，计算显热通量及通过能量平衡余项式计算潜热通量的过程主要参考了SEBS(surface energy balance system)模型和SEBAL(Surface energy balance algorithm for land)模型，这两个模型分别应用于中低分辨率(数公里～250米)和高分辨率(30米或以下)的遥感影像。具体的算法在Bastiaanssen，W G和Menenti，M，Feddes在Journal ofHydrology，1998，212-213：p198-212上发表的论文《A remote sensingsurface energy balance algorithm for land(SEBAL).1.Formulation》以及SU Z在Hydrology and Earth System Sciences，2002，6(1)：p85-99上发表的论文《The Surface Energy Balance System(SEBS)for estimation ofturbulent heat fluxes》上有详细的记载，在此就不具体说明了。

由Penman-Monteith模型推求有效表面阻抗的公式如下：

$r_s=\frac{\mathrm{\Δ}(R_n-G)r_a+86400{\mathrm{\ρ}}_a{C}_p(e_s-e_a)}{\mathrm{\γ\λET}}-\frac{\mathrm{\Δ}{r}_a}{\mathrm{\γ}}-r_a---\left(9\right)$

式中：λET是参照腾发量，R_n是冠层表面净辐射，G是土壤热通量，ρ_a是空气密度，e_s是饱和水汽压，e_a是实际水汽压，Δ是饱和水汽压-气温曲线斜率，γ是湿度计常数，r_s是表面阻抗，r_a是空气动力学阻抗，简化为2米风速的函数。

在本发明所使用的遥感地气交换模型中，将微观尺度的SEBAL与宏观尺度的SEBS模型相结合，发挥了不同模型各自的特点，在遥感数据源的时空分辨率上取得了良好的平衡。具体而言，SEBAL模型在应用时需要对水体、植被与裸地进行分类，从而在不同的植被类别上应用不同的经验公式。由于NDVI的饱和性，其分类效果并不理想，特别是在不同时相的图像中，植被与裸地之间的分界阈值值域比较模糊，变化也较大。在本发明中，对于植被与裸地的分类使用叶面指数来完成。因为LAI由土壤调节植被指数SAVI计算得到，消除了土壤背景的影响，对植被与裸地区分能力增强，且LAI有其物理意义，可以通过实测手段对其进行订正，从实际处理效果和理论角度、及长期业务运行化角度均优于NDVI。

LAI对于植被信息中的低值部分(居民地)和高值部分(农作物)的对比度优于NDVI，NDVI反映的植被信息虽然趋势与LAI相同，但主要的信息被压缩在相对较小的值域中，这是由于没有消除土壤背景影响的缘故，可从图2(a)～图2(d)所示的直方图了解这一情况。

由图2(a)～图2(d)可以看到，由于NDVI的饱和性，植被部分的波峰出现在600(即0.6)左右，植被的信息被压缩在相对小较的区段中，而居民地的信息则并未形成波峰，分散在200-400地段之间。这与LAI的特性是根本不同的：在LAI直方图中，裸地和居民地信息主要处在500(即0.5)的左侧，并形成了波峰，被压缩在相对较小的值域中；而植被处于500的右侧，因植株生长状况的不同表现出分散分布的特征，有利于植被信息的进一步量化和提取。

对于邻近时相的LAI和NDVI图像，植被结构变化不大，植被生长状况也近似，对于遥感提取植被指数而言，关键是要表现出一定的稳定性。由图3(a)～图3(d)所示的直方图可知，邻近时相的LAI直方图的形状、分布、值域略有差异，但大体相似；NDVI图像则出现了一定的变化：31日相对于23日，由于作物的进一步生长，低值部分继续向高值比较移动，但LAI中居民地和裸地信息仍然能形成峰值，裸地与植被之间的过渡明显，而NDVI中的低值部分则迅速降低，过渡不明显，由于NDVI的饱和性，大量的植被像元集中在600左右的区段里，植被内部的差异不明显。

从以上结果来看，LAI相对于NDVI，无论是裸地与植被之间的区分，还是植被部分内部的量化，都具有较好的区分度，有利于遥感应用。

四.运用时间重建技术从瞬时间断的遥感观测中获取完整的植被生长期过程曲线

在获取全年晴好日的LAI的基础上，通过HANTS(Harmonic Analysisof Time Series，时间序列谐波分析)谐函数分解方法得到逐日LAI影像，其核心工作是最小二乘法拟合和傅里叶变换。谐函数分解方法对时间序列图像的要求不像其它的滤波方法严格，可以是不等间隔的、残缺时序遥感影像。它充分利用遥感图像存在时间性和空间性的特点，将其空间上的分布规律和时间上的变化规律联系起来共同分析，因此不但可以作为时序遥感数据去云重构的分析手段，对于定量反映植被动态变化也有重要作用。

在具体进行重建时，由于时序LAI的云噪声对LAI曲线的周期特征影响很大，是误差的主要来源，因此首先舍弃云噪声影响的异值点，对LAI特征点进行样条插值内插到每日，并使用大于5的张力系数来保留时间特征点，然后逐像元地对LAI进行HANTS平滑。

五.在逐日尺度上建立表面阻抗模型以获取逐日地表阻抗

在本方法中，为获得逐日连续的蒸散量，引入P-M模型，将有晴好图像日的蒸散结果作为“关键帧”，通过将关键帧的地表阻抗信息、逐日遥感信息、气象信息，建立地表阻抗模型反推逐日地表阻抗，以便利用逐日的气象数据基于P-M模型计算得到逐日的蒸散量数据。

以下是地表阻抗模型表达式：

${\mathrm{RS}}_{\mathrm{unc}}=\frac{{\mathrm{LAI}}_{\mathrm{clr}}\·{\mathrm{RS}}_{\mathrm{clr}}}{{\mathrm{LAI}}_{\mathrm{unc}}\×m\left(T_{\min }\right)\×m\left(\mathrm{VPD}\right)}---\left(10\right)$

其中RS_unc是无晴好图像日需要确定的阻抗项，而RS_clr是晴日阻抗项，LAI是对应日的叶面指数。对于影像缺失日的LAI，采用了HANTS谐函数分解方法来获得其时间序列，这种方法对于具有周期性的植被生长最具优势。

另外，式(10)中的m(T_min)和m(VPD)是两个用以表达极端温湿条件对植物气孔开闭的制约函数，定义如下：

$m\left(T\min \right)=\left\{\begin{array}{c}1.0\\ \frac{T\min -T{\min }_{\mathrm{close}}}{T{\min }_{\mathrm{open}}-T{\min }_{\mathrm{close}}}\left[\begin{array}{c}T\min \&GreaterEqual;{T\min }_{\mathrm{open}}\\ T{\min }_{\mathrm{close}}<T\min <{T\min }_{\mathrm{open}}\\ T\min \&le;{T\min }_{\mathrm{close}}\end{array}\right]\\ 0.1\end{array}\right.$

$m\left(\mathrm{VPD}\right)=\left\{\begin{array}{c}1.0\\ \frac{{\mathrm{VPD}}_{\mathrm{close}}-\mathrm{VPD}}{{\mathrm{VPD}}_{\mathrm{close}}-{\mathrm{VPD}}_{\mathrm{open}}}\left[\begin{array}{c}\mathrm{VPD}\&le;{\mathrm{VPD}}_{\mathrm{open}}\\ {\mathrm{VPD}}_{\mathrm{open}}<\mathrm{VPD}<{\mathrm{VPD}}_{\mathrm{close}}\\ \mathrm{VPD}\&GreaterEqual;{\mathrm{VPD}}_{\mathrm{close}}\end{array}\right]\\ 0.1\end{array}\right.$

其中的Open和Close表示的是在对应的温湿状态下植被气孔的临界状态。常数的取值可以由大量的文献调研中取得，如表1所示，并需要结合相应的土地利用图使用。

表1

六.反复运用P-M模型求得逐日蒸散量，由此进行时间累计和尺度下推，从而最终确定区域地表蒸散量

经过时间重建后，得到的日蒸散量是较低分辨率的。利用高分辨率的遥感数据得到的蒸散量是较高分辨率的，但只有TM卫星过境日的蒸散量结果。

为了使反演的蒸散量结果既有较高的空间分辨率又有较高的时间分辨率，需要实现不同分辨率蒸散量数据的融合。其实质是充分利用高分辨率遥感影像提供的蒸散量空间分布信息和低分辨率遥感影像提供的蒸散量时间变化信息，从而在更精细的时间空间尺度上反映蒸散量的时空变化。

为此，本发明在由分辨率为1公里的遥感影像估算的每半月蒸散数据的基础上，并参考作物物候情况。在每月选择一幅具有代表性的高分辨率蒸散量图像；如果没有这样的图像，就从植被覆盖情况相似的邻近月份中选取，如果无法确定具有代表性的高分辨率蒸散量图像，就选择前后相邻的图像，按照时间距离进行加权，形成本月的高分辨率蒸散量分布参考图像，即：

${\mathrm{ET}}_{30m}^r=c1\&times;{\mathrm{ET}}_{\mathrm{tm}}^a+c2\&times;{\mathrm{ET}}_{\mathrm{tm}}^b---\left(11\right)$

其中c1，c2为加权系数，设D1为a图像到本月参考点(每月15号)的距离，而D2为b图像到本月参考点(每月15号)的距离，则c1＝d2/(d1+d2)，c2＝d1/(d1+d2)。

对于高分辨率蒸散量分布图ET_30mr，要受到来自本月1km蒸散量图像所给出的蒸散量总量的限制，即得到用于时间外推的ET_30m ⁱ，这一ET_30m ⁱ即代表了本月蒸散的分布和量级：

${\mathrm{ET}}_{30m}^i=\frac{\stackrel{\&OverBar;}{{\mathrm{ET}}_{1\mathrm{km}}^i}}{\stackrel{\&OverBar;}{{\mathrm{ET}}_{30m}^r}}\&times;{\mathrm{ET}}_{30m}^r---\left(12\right)$

使用Penman-Monteith模型可从ET_30m ⁱ导出表面阻抗r_{s_30m} ⁱ。即：

$\mathrm{\&lambda;}{\mathrm{ET}}_{s\_30m}^i=\frac{\mathrm{\&Delta;}(R_n-G)+{\mathrm{\&rho;}}_a\mathrm{Cp}\frac{e_s-e_a}{r_a}}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+\frac{r_{s\_30m}^i}{r_a})}---\left(13\right)$

由高分辨率r_{s_30m} ⁱ图像，结合对应日的气象数据，使用Penman-Monteith模型计算逐日或时段平均蒸散量，并累积蒸散量到月或季，就可以最终确定区域地表蒸散量。

表2显示了利用本区域蒸散量监测方法在在2003年禹城暖温带半湿润旱作农田通量观测站的蒸渗仪数据基础上做的模拟计算结果与实测结果的比较情况，其中同步气象数据和蒸渗仪(Lysi)数据来自该地的同期观测。

月	ET模拟值	蒸渗仪(Lysi)
			1	8.57	17.16
2	18.95	25.42
			3	42.64	43.98
4	88.46	127.55
			5	131.64	202.71
6	87.84	77.02
			7	107.16	91.53
8	97.02	125.13
			9	77.35	86.54
10	41.56	49.25
			11	14.21	22.95
12	11.49	10.16
			小麦季	401.44	537.69
玉米季	325.46	341.71
			全年	726.89	879.40

表2

表2说明本区域蒸散量监测方法的结果精度可达到85％。除此之外，目前已运用多组地面实测数据对本区域蒸散量监测方法的计算结果进行了验证。例如对于小麦生长期内的逐日蒸散结果(禹城暖温带半湿润旱作农田通量观测站，2003年5月20～30日；北京昌平区小汤山国家精准农业研究示范基地，2004年5至7月，2005年5至6月)，逐日偏差在15％以内；对于县级的生长期内的总耗水量结果(河北馆陶，2005年观测作物耗水量；密云新河子，2002-2005观测耗水量)，偏差在10％以内；对于流域级的年水量平衡结果(2002-2005年海河水平衡数据，2002-2005年滦河山区、北京地区、永定河上、永定河到三家店区间、子牙河山区的水平衡数据)，偏差在8％以内。

以上对本发明所述的基于遥感的区域蒸散量监测方法进行了详细的说明。对本领域的一般技术人员而言，在不背离本发明实质精神的前提下对它所做的任何显而易见的改动，都将构成对本发明专利权的侵犯，将承担相应的法律责任。

Claims (9)

1.一种基于遥感的区域蒸散量监测方法，其特征在于包括如下的步骤：

(1)获取特定区域的多时相多分辨率卫星遥感影像数据；

(2)对地面的关键参量进行遥感反演；

(3)通过遥感地气交换模型计算地气能量交换过程中每个像元的阻抗，并结合气象要素计算每个像元的地气温差，计算出用于加热每个像元上方空气的显热通量，通过能量平衡余项式计算潜热能量，并由Penman-Monteith模型推求逐日有效的表面阻抗；

(4)运用时间重建技术从瞬时间断的遥感观测中获取完整的植被生长期过程曲线；

(5)在逐日尺度上建立表面阻抗模型以获取逐日的地表阻抗；

(6)反复运用Penman-Monteith模型求得逐日蒸散量，由此进行时间累计和尺度下推，最终确定区域地表蒸散量。

2.如权利要求1所述的区域蒸散量监测方法，其特征在于：

所述步骤(1)中，将ASTER/TM和AVHRR/MODIS作为卫星遥感影像数据源。

3.如权利要求1所述的区域蒸散量监测方法，其特征在于：

所述步骤(2)中，所述关键参量至少包括地表反照率、净辐射量、地表温度、植被覆盖度、叶面指数和动量粗糙度。

4.如权利要求3所述的区域蒸散量监测方法，其特征在于：

在计算所述净辐射量时，直射辐射的系数为0.45，散射辐射的系数为0.1。

5.如权利要求3所述的区域蒸散量监测方法，其特征在于：

在计算所述动量粗糙度时，分别对植被、地形起伏、非植被覆盖表面下的几何粗糙度进行加权，并在计算非植被覆盖表面下的几何粗糙度时引入雷达数据。

6.如权利要求1所述的区域蒸散量监测方法，其特征在于：

所述步骤(3)中，使用叶面指数进行植被与裸地的分类。

7.如权利要求1所述的区域蒸散量监测方法，其特征在于：

所述步骤(4)中，在获取全年晴好日的叶面指数的基础上，通过时间序列谐波分析谐函数分解方法得到逐日的叶面指数影像，并对叶面指数的特征点进行样条内值处理，使用大于5的张力系数来保留时间特征点。

8.如权利要求1所述的区域蒸散量监测方法，其特征在于：

所述步骤(4)中，在植被生长期过程曲线中首先舍弃云噪声影响的异值点。

9.如权利要求1所述的区域蒸散量监测方法，其特征在于：

所述步骤(5)中，首先在预定时间段中选择一幅具有代表性的高分辨率蒸散量图像；如果没有具有代表性的高分辨率蒸散量图像，就从植被覆盖情况相似的邻近时间段中选取；如果无法确定具有代表性的高分辨率蒸散量图像，就选择前后相邻的图像，按照时间距离进行加权，形成该时间段的具有代表性的高分辨率蒸散量图像。

Images