CN102034027A - 流域尺度土壤湿度遥感数据同化方法 - Google Patents

Abstract

本发明提供了流域尺度土壤湿度遥感数据同化方法。其步骤为：首先改进流域产流计算模块，发展了一个适合同化遥感土壤湿度信息的具有土壤水动力学过程描述的分布式水文模型；然后引入来源于信息科学的粒子滤波顺序数据同化的方法，在分布式水文过程数值模拟的动态运行过程中不断融合、同化新的遥感观测数据，顺序同化中得到更新的流域土壤湿度同化数据，将其反馈至分布式水文模型平台中，逐步估计流域土壤水分含量在时间和空间上的分布格局。实践证明本发明不但能够为水文、生态、环境和农业研究提供高精度的具有物理一致性的流域土壤湿度数据集，亦为利用遥感反演得到的土壤表层土壤湿度数据进行四维数据同化处理，提高模型的精度奠定了基础。

Description

流域尺度土壤湿度遥感数据同化方法

技术领域

本发明涉及一种遥感数据同化方法，更具体的说是流域尺度土壤湿度遥感数据同化方法。

背景技术

土壤湿度状况既是气候和水资源系统演化的最终结果，同时自身对气候和土壤植被大气作用系统施加强烈反馈，深刻影响其未来变化。在水文水资源领域，获取准确的土壤水分空间分布，可作为水文模型的洪水模拟的土壤含水量初始值，提高洪水预报的精度。同时经过校正的连续的精确土壤水分数据，对于长期水文过程模拟，径流预报和预测水资源的变化有潜在应用价值。发展一个适合同化遥感土壤湿度信息的具有土壤水动力学过程描述的分布式水文模型为模型算子，遥感反演表层土壤含水量为观测算子，结合合适的同化算法发展的四维数据(4DDA)同化系统所获取的具有物理一致性的高分辨率时、空土壤湿度数据集，在大气科学、生态水文、环境、农业等领域均有广泛的应用和需求。

流域尺度和田间尺度土壤水分状况的研究有利于加深人们对水文过程的理解并为最终掌握地表水和地下水资源的复杂交互作用提供帮助。在水文水资源领域，准确的土壤水分空间分布，可作为分布式水文模型洪水模拟的所需的土壤湿度变量的初值，提高洪水预报的准确性；另外，将经过校正的连续准确的土壤水分数据代入水文模型，对于提高长期水文过程模拟和径流预报的精度以及准确估算地表水资源量有潜在价值。因此论文以流域尺度的土壤湿度变量作为主要研究对象。

流域土壤湿度的研究离不开模型和观测两种手段。首先，水文模型是计算和分析土壤含水量的有效工具，Western et al.强调了水文模型在模拟土壤水分相关的时空过程中具有重要意义。近年来，随着计算机科学和信息技术的发展，RS和GIS技术支持下的分布式流域水文模型成为研究的热点。RS和GIS的DEM处理技术提供了分布式水文模型所需的不同时空分辨率的水文气象参数和下垫面植被、土壤资料以及流域地形参数，如坡度、坡向、水流路径和流域边界等参数；结合RS、GIS、DEM等空间信息技术建立的分布式流域水文模型，可在更小的水文模拟单元上，精确地模拟复杂气候条件和地表状况下的水文循环过程，因而分布式水文模型可以提供流域中多种水文变量的空间输出(如ESSI模型可得到土壤含水量、径流深、蒸散发等23个水文过程变量)。这与主要通过概念参数和土壤湿度变量来模拟流域出水口流量的集总式水文模型有显著区别。正因为以上特点，分布式水文模型不仅能够为相关的水文生态研究如流域产、输沙、营养物输移、污染物扩散以及人类活动对水循环的影响等提供先进的计算和模拟平台，而且在洪水预报，如与气象预报模式耦合延长洪水灾害的预警时间方面有着重要的用途。但是，分布式水文模型还存在一定的不足，例如，由于受到实验室模型的动力学方案在流域上应用还有待改进，如何建立物理基础更强的分布式水文模型，获取更加客观真实的土壤水分含量等水文参数的时空格局是分布式模型研究的重要课题；另外，由于受到模型输入数据误差和模型自身结构误差的影响，分布式水文模型的模拟结果还存在的一定的不确定性，如何量化并减少水文模型预测的不确定性，也是水文科学当前研究的前沿热点。利用遥感数据通过数据同化方法降低模型的模拟误差积累，是获取高精度和高分辨率土壤湿度数据的很有前景的方法。

其次，土壤水分可以通过观测直接获取，土壤水分直接监测方法大致可分为地面调查和机载/卫星遥感两种。基于地面调查的监测方法根据获取数据方式和手段不同，又可以分为土钻取土称重法、烘干法、中子仪法、电阻法、TDR法，这些方法采用空间上定点观测，采样速度慢、人力物力财力消耗大，难以满足大面积连续动态监测的需要，受到其空间采样密度的限制和大气水分能量过程随机扰动的影响，土壤湿度的实际空间格局也无法准确观测。与地面调查相比，遥感技术具有大面积同步观测、时效性、经济性的特点，是解决土壤湿度观测的有效途径。随着遥感技术向高空间分辨率、高光谱分辨率方向发展，遥感技术在土壤水分大面积动态监测中将发挥日益重要的作用。

可用于监测土壤湿度状况波段包括可见光与近红外、热红外波段和微波波段，其中可见光与近红外波段数据主要通过上述的归一化植被指数法、温度状态指数法、TVDI法计算反演得到土壤湿度；热红外波段数据可采用热惯量法获取湿度信息；微波数据则主要通过亮度温度与土壤湿度的经验统计关系、基于物理的微波辐射传输模型反演的方法得到土壤湿度。与可见光与红外波段相比，利用微波数据反演土壤湿度具有坚实的物理基础，其应用远较可见光与红外波段的数据广泛，反演可靠性和精度也更高。用微波数据反演土壤湿度的优势主要体现在以下几点；1)微波的物理特性。微波遥感土壤水分反演具有坚实的物理基础：地物微波比辐射率主要决定于目标物的介电常数，而土壤介电常数主要决定于土壤的水分含量。在微波波段，水的介电常数大约为80，干土仅为3，随着水分含量的变化，土壤的比辐射率从湿土的0.6(30％体积土壤湿度)到干土的0.9(9％体积土壤湿度)之间变化，它们之间具有较大的反差，因此微波图像对水分十分敏感。2)地区适应性。在我国南方，由于气候湿润，多云多雨，限制了光学遥感数据的及时获取，影响了调查任务的按时完成。我国南方多云多雨地区的土壤湿度调查急需具有全天候对地观测能力的高分辨率雷达遥感数据源。3)由于土壤湿度变化多发生在阴雨天气，在阴雨天气这一重要的湿度特征变化的时间窗口，微波数据是唯一可行的土壤湿度调查方法。微波遥感观测土壤湿度也有一定的局限性，由于受到植被的影响和微波波长的限制，目前星载微波遥感能够探测得到的土壤湿度仅为地面表层(仅几厘米)的水分含量。

目前卫星观测得到的土壤湿度仅限于上层(2-10厘米)的土壤，无法将其直接用于水文模型中，因为水文模型的模拟需要的是整个土壤剖面的水分含量。通过同化系统集成卫星观测数据于分布式水文模型中，可以提供准确的物理一致的土壤剖面水分含量，有利于改进流域水循环过程的模拟。

构建适合于土壤湿度模拟的，具有土壤水过程动力学描述的分布式水文模型，在流域强化测试和验证的基础上，将其发展成可有效同化遥感土壤湿度信息的分布式水文模型同化平台，同时利用现有的对地观测卫星的资料对水文模型中的关键变量，土壤湿度进行四维数据同化处理。通过连续水文模拟和同化实验，获取高精度高分辨率的土壤湿度时空数据集。

模型和观测误差统计特征的确定是数据同化的关键，误差特征只有为数据同化方法准确把握后，才能得到正确有效的同化结果，对状态作出最优估计，数据同化算法理论研究和应用一直是同化研究的重要组成部分。当前大气、海洋和陆面数据同化中用到的数据同化方法主要有最优插值法、四维变分法、集合卡尔曼滤波和粒子滤波等。这些算法主要可顺序同化和变分同化两类，顺序同化根据前一时刻的模式状态生成当前时刻模式状态的预报值，在观测更新阶段，引入观测数据，利用最小方差估计方法对模式状态进行重新分析。而变分方法在同化窗口内建立目标函数，并通过调整模型预报值以极小化与观测数据之间的差值达到最优化状态的目的。以下分别对各主要方法的研究作一介绍，

(1)最优插值法。最优插值又称统计插值法，该方法以线性最小平方估计理论为基础，通过给定最佳线性无偏估计方程中的权矩阵，建立气象场的统计协方差函数。其优点是实施简单，需要的计算量相对较小，可处理精度不同的各种观测数据，考虑了观测和预报之间的线性相关及观测要素固有的统计结构。由于易于应用，可操作性好，最优插值法目前仍然为大多数气象业务化同化系统所采用。该方法主要缺点是当模型同化采用不同类型和来源的观测数据集时，会使分析场产生伪噪声，无法确保大尺度和小尺度分析间的一致。

(2)四维变分法。它利用变分法的思想，将同化问题转化成对一个极值方程的求解，通过全局的调整模型模拟使其与同化期所有的观测值的差值和最小。变分数据同化的优点是不需要计算模式预报误差的方差阵，具有较高的运行效率。近年来在一些陆面数据同化研究中也采用了该方法。但变分同化方法也存在一些明显的缺陷，如该方法必需的预报模型切线性伴随方程的求解十分耗时且容易出错，甚至在理论上是不存在的。尤其是当需要对预报模型进行改进时，其伴随方程必须随之进行更复杂的改进和测试，这些特点制约了它在陆面模型同化中的应用。

(3)集合卡尔曼滤波(EnKF)。集合卡尔曼滤波算法是顺序同化方法中的一种，是扩展卡尔曼(EKF)滤波方法的一种Monte Carlo近似。该算法的最大特点是利用MonteCarlo集合预报方法估算状态的预报误差协方差，通过模型状态集合代表随机动态预报中的概率密度函数，其中集合的均值代表最优的状态估计，集合的分布代表误差的方差。EnKF算法不需要模型算子的切线性模式和伴随模式，克服了扩展Kalman滤波要求线性化的模型算子和在一些复杂非线性模型同化中不稳定的缺点，与四维变分(4DVAR)和EKF同化方法相比，基于EnKF的数据同化易于实现，可移植性强，被广泛地应用到大气、海洋和陆面数据同化研究中。

(4)粒子滤波。粒子滤波也称为顺序蒙特卡罗SMC，属于顺序同化方法的一种。该算法的优点是可以得到误差后验分布的全部信息，准确计算其概率分布的一些统计特征，如均值、方差、峰度等，并把握模型系统的非线性，相对于EnKF，粒子滤波求取了整个后验分布，而EnKF只求取后验分布的均值和方差。因此它是一种能处理非线性非高斯问题的滤波方法，在通信理论，信号处理和目标跟踪领域，有大量的应用，粒子滤波有几个很重要也必须好的问题，一是粒子的产生，常用的是序贯重要性采样(SIS)，序贯重要性重采样(SIR)等；其次是怎么解决好粒子退化和粒子的多样性问题。

发明内容

本发明的目的是提供流域尺度土壤湿度遥感数据同化方法，构建适合于土壤湿度模拟的，具有土壤水过程动力学描述的分布式水文模型，在流域水文过程模拟和验证的基础上，将其发展成可有效同化遥感土壤湿度信息的分布式水文模型同化平台，利用现有对地观测卫星的资料对水文模型中的关键变量，表层土壤湿度进行四维数据同化处理。以获取较高时间精度和空间分辨率的流域土壤湿度时空同化数据集。

本发明的目的是通过以下步骤来实现的：

流域尺度土壤湿度遥感数据同化方法，步骤为：

A)数据准备；

B)土壤湿度反演；

C)分布式水文模型同化平台构建；

D)流域尺度土壤湿度同化。

步骤A)中的数据准备，包括：

1)水文、气象数据；

2)DEM(数字高程模型)数据；

3)植被、土壤参数库，流域土地利用/覆被和土壤类型这些下垫面参数数据，并根据流域所在地区的“土壤志”收集了研究区各土壤亚类的典型土壤剖面数据，建立了模型所需的土壤参数空间数据库；

4)卫星影像数据，包括主动微波ENVISAT-ASAR遥感数据和MODIS影像数据。

5)逐日径流和土壤湿度实测数据

步骤B)实现的具体流程见图2。

可见光遥感：温度植被干旱指数TVDI可由图像数据派生的植被指数和地表温度计算求得，由以下公式给出：

$\mathrm{TVDI}=\frac{\mathrm{LST}-{\mathrm{LST}}_{\min }}{{\mathrm{LST}}_{\max }-{\mathrm{LST}}_{\min }}$

LST为地表温度；LST_min为相同NDVI值的最小地表温度，对应的是LST-NDVI特征空间的湿边；LST_max为相同NDVI值的最大地表温度，对应的是LST-NDVI特征空间的干边。

ASAR湿度反演：

微波散射过程

入射的主动微波与自然地表之间存在着复杂的交互作用，雷达的后向散射系数除了与微波系统的物理参数有关之外，还与地物的介电特征和几何特征密切相关。公式列出了主要的影响后向散射的主要因子。理解微波和地物之间的相互作用的原理，将有利于从主动微波数据中反演得到地表参数。

σ⁰＝f(pp，θ，λ，∈r，r)

其中pp表示微波数据的极化特征，θ是入射角，λ是入射波长，∈_r是介电常数，r是地表的粗糙度，实际中较难以测定，通常由均方根高度k和相关长度l来表示自然地表粗糙度。

目前已经有多种地表散射理论模型被提出，用于模拟散射系数与各物理参数之间的相互作用，包括几何光学模型、物理光学模型、SPM(Small Perturbtatoin Model)模型和IEM(Integral Equation Model)模型。这些模型都将均方根高度、相关长度和介电常数作为模型的输入变量，来计算模拟散射作用。由于地表参数的复杂性，很多模型都对地表粗糙度的输入范围做了较强限定，因此在实际应用中的土壤特征往往超出了模型的适用范围。针对这一情况，Fung et al.提出了IEM模型，这个模型可以适用于具有更宽范围粗糙度的地表，并且在研究中被证明可以较好的模拟自然地表的后向散射系数参数。

IEM模型适用于裸土地表的后向散射模拟，而在植被影响显著的区域，则需要结合植被散射模型，对植被散射的影响加以考虑

ASAR数据预处理：

辐射校正

根据雷达方程，后向散射系数可由辐射强度值(Intensity)计算得到。对于ENVISATASAR数据，影像DN值和后向散射系数之间的关系式如下

${\mathrm{\σ}}^0=\frac{{\mathrm{DN}}^2}{k}\sin \mathrm{\θ}$

其中σ⁰是线性单位的雷达后向散射系数，值域在0～1之间，k是外部校正参数，可以从ASAR影像的头文件中得到，θ是空间分布的入射角。分析数据时常用以分贝为单位的σ⁰，因此需要对上式得到以线性单位表示的σ⁰进行对数转换，

σ⁰(dB)＝10[log(σ⁰)]

几何校正

需要对获取的属于一级产品的ASAR影像进行几何校正，这主要分为两个步骤，首先是几何粗校正，其次是与数字高程DEM影像数据匹配的几何精确校正。

ESA的雷达处理影像工具NEST软件提供了ASAR影像的几何校正功能，在利用该工具对ASAR数据进行粗校正之后，进一步通过与DEM的比对，以DEM为基准影像，进行图像对图像的精确校准，如此循环后将所有的ASAR数据统一的精确订正到DEM的坐标系下。几何校正是必须的，这一步骤保证了遥感反演结果用于与模型对比和模型同化的空间精度。

光斑滤波

雷达影像通常都有不同程度的光斑(speckle)噪声，这是由于雷达系统的图像分辨率通常大于信号的波长和独立散射体的大小，而象元内不同散射体的散射之间的相长干涉和相消干涉导致了后向散射在图像上的随机分布。光斑噪声在图像上看起来类似于椒盐状分布，使图像的可解译性变差。有效的降低噪声对于计算雷达影像的后向散射系数十分关键，目前去除噪声的方法主要有多视方法(multi-look)和空间滤波方法(Specklefiltering)。

利用Loew et al.提出的针对ASAR WS数据的湿度估算模型，根据研究区特性和实测土壤湿度数据，对估算模型的相关参数做了一定的修正，得到了适合于半干旱研究区沂河流域的ASAR土壤湿度估算的模型参数。反演了沂河流域连续的ASAR WS影像数据的土壤湿度，总体结果满意。

详细的ASAR数据湿度估算方法如下：

散射系数归一化

Loew et al.指出不同覆盖的地表后向散射与入射角存在着线性相关关系。IEM辐射传输模型的数值模拟试验也验证了这一关系的存在，地表后向散射与入射角关系可以用下式表达，

σ⁰＝pθ+q

此处，θ表示入射角，p、q是与土地覆盖有关的经验参数，具体数值可由Loew etal.中查到。尽管每一景ASAR WS模式影像的入射角各不相同，同一景影像中的不同区域入射角变化范围也较大，可通过下面公式，将所有影像的后向散射系数归一化至参考的入射角θ₀，此处θ₀为ERS的入射角度值23度，具体的归一化方法如下，

${\mathrm{\σ}}_{\mathrm{nor}}^0=\frac{{\mathrm{p\θ}}_0+q}{\mathrm{p\θ}+q}{\mathrm{\σ}}^0$

其中，σ⁰为图像原始后向散射系数，

为归一化后的图像后向散射反演参数

众多研究表明，后向散射系数与地表介电常数的相关关系具有如下的正相关关系，这一经验关系主要由设在德国的20个观测试验站的经验数据库得出，并被用于ERS和中尺度的雷达数据的土壤湿度反演，利用该方法得到的表层土壤湿度的反演的精度(RMSE)在0.04m3/m3到0.07m3/m3之间。本发明中利用上一步归一化后处理得到的ASAR后向散射，结合式2-6可以得到土壤的介电常数，最终反演得到表层的土壤湿度值，

∈_r＝a+bσ⁰+c(σ⁰)²

其中∈_r为介电常数，a、b、c为经验参数，与土地覆盖类型相关，由于受到研究区差异，根据不同的土地覆盖类型，对此公式的参数做一定的修正，其具体值见下表。

本发明中使用的ASAR遥感土壤湿度估算参数

介电模型

由于介电常数与土壤湿度具有较好的相关性，因此通过上述方法获取得到的介电常数结果，可以利用介电模型将其转化为需要的土壤湿度信息。

目前已发展了很多的理论和半经验模型用于描述自然地表条件下介电常数与土壤湿度的关系。其中在遥感领域最为常用的是Hallikainen et al.提出的介电常数模型，该模型给出的土壤介电常数与土壤湿度之间的函数关系如下，图2-6给出了由该模型模拟得到的介电常数和土壤湿度相关关系(砂粒比重51.5％，粘粒比重为13.4％)，

∈_r＝(a₀+a₁S+a₂C)+(b₀+b₁S+b₂C)m_v+(c₀+c₁S+c₂C)m_v ²

其中，m_v是需要计算的体积含水量，S和C是土壤的砂粒和粘粒的重量百分比，即土壤质地参数，可以由我们建立的流域土壤物理参数库得到，a_i，b_i和c_i是模型的经验常数，已由文献给出。

基于上述模型，利用获取的流域土壤介电常数图像和已有的土壤质地参数，通过迭代计算方法(误差阈值为0.01)逐栅格推求得到最终的ASAR遥感数据的表层土壤湿度估算结果。

步骤C)简化了流域水文产流模拟，在此基础上对纵向水文过程进行较为细的建模，主要描述了土壤水分传输、土壤水与地下水的交互以及根系的水分提取作用。其中涉及到的主要过程变量包括：冠层截留，地表产流，土壤水分传输，潜在蒸散发，植物蒸腾和土壤表面蒸发，地下水与不饱和土壤水交互，地表汇流，地下汇流。

建模思路：

模型将流域看成由不同厚度的土壤柱体组成，而土壤层下方为饱和地下水。土壤柱体的厚度是由流域实际的土壤调查资料给定，主要根据不同土壤类型的典型剖面实地调研资料确定，因此不同土壤类型的柱体厚度并不相同。这一特点决定流域土壤蓄水容量的空间异质性，通过格点土壤柱体的蓄满产流与流域汇流计算，模拟得到流域总出口流量。

土壤水分传输和地表径流的产生过程主要发生和决定于根区的土壤层，因此该土壤体的描述至关重要。模型中，将流域的土壤柱体按照相同的分层标准划分为不超过10层的土壤层，最大土壤深度为3.43米，因而格点上的土壤柱体由于土壤厚度的差别具有不同的土壤层数。土壤各层的土壤质地参数来自于土壤调查资料所建立的土壤参数库，同样也是根据不同土壤类型的典型剖面的实地调研资料得到。在获取土壤分层参数的基础上，利用至地面的降水量和产流差得到土壤表面的平均下渗量，基于不饱和土壤水动力方程，模拟各层土壤的水分状态和层间的纵向水分传输，模型还同时考虑了植被根系吸水作用和土壤表面蒸发对土壤湿度的影响。

模型的具体构建框架如附图6所示，其中k1...k3表示土壤水力传导度，θ1...θ3表示土壤水分含量，Φ1...Φ3表示土壤水势，Q1...Q3表示土壤水流，dz1...dz3表示土壤层厚度，图中给出的水文过程均包含于新模型之中

水文过程概化：

林冠截留量计算公式如下：

$I_{\mathrm{int}}=C_m\·(1-e^{-K{\mathrm{PC}}_m})+e_w\·T$

其中C_m为冠层截留能力，k为降雨截留系数，P为降雨量(mm)，T为降雨历时(hr)，e_w为林冠蒸发率，因此e_w·T表示树冠蒸发以及树干蒸发之和。

冠层截留能力C_m采用Von Hoyningen-Huene所提出的方法计算：

C_m＝C_p·(0.935+0.498·LAI-0.00575·LAI²)

其中Cp为植被盖度，LAI为格点单元平均的逐日叶面积指数

降雨截留系数k是植被叶面积指数的函数，由下式给出：

k＝0.046·LAI

林冠蒸发率e_w与降雨初期的干燥程度和风速等因素相关，在无实测数据验证时，通常取以下经验参数：亚热带地区取0.031～0.047，温带湿润、半湿润地区取0.063～0.093，温带干旱、半干旱地区取0.170～0.200。

水量平衡计算公式如下：

P-Es-Et-G＝Δw+Δc+R

其中，P是降水(包括降雨和降雪)，Es是土壤表面蒸发，Et是植物蒸腾，G是地下水出流，Δw是土壤水的增量，Δc是冠层截留，R是地表径流。

土壤水分传输中涉及到不饱和土壤水通量定律Buckingham-Darcy方程和连续方程，Buckingham-Darcy方程如下：

$q=-K\left(h\right)(\frac{\∂h}{\∂z}+1)$

其中q是不饱和水流通量，h是土壤水基质势，K(h)是水力传导度，z是垂直方向的深度。

连续方程如下：

$\frac{\∂\mathrm{\θ}}{\∂t}+\frac{\∂q}{\∂z}=0$

基于垂向流的假设，公式中的q由Buckingham-Darcy方程替换，求得解即为一维Richard’s方程表达式：

$\frac{\∂\mathrm{\θ}}{\∂t}=\frac{\∂}{\∂z}\left(K\left(h\right)\frac{\∂h}{\∂z}\right)+\frac{\∂K\left(h\right)}{\∂z}$

其中t是时间，θ是土壤体积含水量。

根系土壤水分提取，此土壤水含量计算中需考虑根系吸水作用，本模型将其作为S项添加到Richard’s方程中。本发明使用以下方法考虑不同土壤层的根系水分吸收作用：

$\frac{\∂\mathrm{\θ}}{\∂t}=\frac{\∂}{\∂z}\left(K\left(h\right)\frac{\∂h}{\∂z}\right)+\frac{\∂K\left(h\right)}{\∂z}-S$

其中S表示土壤根系吸水或者土壤表面蒸发，由以下公式描述

$S_i=E\frac{r_i{w}_i}{\mathrm{\Σ}{r}_i{w}_i}$

其中i是模拟土壤层的序号(i≥2)，对于第一层土壤(i＝1)，S_i为表层土壤表面蒸发，E_veg是植物的蒸腾，r_i是土壤层中根系分布参数，与植被类型有关，w_i是土壤水饱和系数，当湿度处于凋萎点时等于0，完全饱和时为1，由下式给出：

其中

是凋萎点的水势，等于-1.5×10⁵，

是饱和土壤水势，

是土壤水基质势。

潜在蒸散发，当流域存在蒸散发观测站时，潜在蒸散发计算由观测数据插值得到；当流域无蒸散发的观测时，利用气象观测数据根据Penman-Monteith公式计算得出，Penman-Monteith公式由下式给出：

$E=\frac{1}{\mathrm{\λ}}\left[\frac{\mathrm{\ΔA}+{\mathrm{\ρ}}_a{c}_{p}D/r_a}{\mathrm{\Δ}+\mathrm{\γ}(1+r_s/r_a)}\right]$

其中E是蒸发量，λ是水汽潜热，Δ是饱和水汽压温度曲线的斜率，A是净辐射强度，ρ_a是湿空气的密度，c_p是常压下的空气汽化潜热，D是水汽压与饱和水汽压的差值，r_a是空气动力学阻力系数，γ是干湿计常数，r_s是作物叶面的阻力系数。

植物蒸腾，蒸腾量计算如下：

$E_{\mathrm{veg}}=E\frac{\mathrm{LAI}}{3}\left(\frac{1+p}{2}\right)$ 0＜LAI≤3.0

$E_{\mathrm{veg}}=E\left(\frac{1+p}{2}\right)$ LAI＞3.0

其中E_veg是植物的蒸腾量，E是潜在蒸散发量，p是格点上的植被覆盖率，LAI是格点上的逐日叶面积指数。

植被覆盖率p通过逐日叶面积指数LAI来计算，公式如下：

p＝1-e^-α·LAI

式中：α为地表太阳辐射的消光系数，缺省值为0.45。

土壤表明蒸发，分别考虑土壤水分蒸发的两个主要阶段，即当土壤水分供给充足，蒸发量较高且相对稳定；当土壤水分降低，土表蒸发随着土壤水分含量呈非线性下降：

$E_{\mathrm{soi}}=E\frac{\mathrm{dz}}{\mathrm{dz}+\exp (2.374-0.00713\·\mathrm{dz})}\exp \left(\frac{2.5\·(\mathrm{\θ}-f_c)}{f_c-w_p}\right)$ θ＜f_c

$E_{\mathrm{soi}}=E\frac{\mathrm{dz}}{\mathrm{dz}+\exp (2.374-0.00713\·\mathrm{dz})}$ θ≥f_c

其中E_soi是土壤表面蒸发量，E是潜在蒸散发，dz是表层土壤的厚度(mm)，θ是表层土壤的水分含量(体积含水率)，f_c是表层土壤的田间持水量(体积含水率)，w_p是表层土壤的凋萎系数(体积含水率)。

地下水与不饱和土壤水交互，为在流域模拟中考虑交互作用，模型主要作了以下三个假设，一是假设深层土壤下的含水层中存在具有一定水深的饱和地下水库，二是土壤层与地下水库之间的水力传导度随着深度增加呈指数下降，三是流域空间上的地下水库深是相同的。地下水的补给计算由以下公式给出：

其中Q是地下水补给量(mm/s)，当水流向下流向地下含水层时为正值，反之为负值，

是地下水位的深度，

是底层土壤水基质势，z_bot是底层土壤的节点深度，K_a是含水层的水力传导度，由下式给出：

$K_a=\frac{{\∫}_{z_{\mathrm{bot}}}^{z_{\▿}}{k}_{\mathrm{bot}}{e}^{-f\left({z-z}_{\mathrm{bot}}\right)}\mathrm{dz}}{z_{\▿}-z_{\mathrm{bot}}}=\frac{k_{\mathrm{bot}}(1-e^{-f(z_{\▿}-z_{\mathrm{bot}})})}{f(z_{\▿}-z_{\mathrm{bot}})}$

其中f是延时因子，用于水力传导度的计算。

地下汇流，地下水运动是渗流水力学的问题，Boussinesq根据连续方程和达西定律，建立了饱和带中的一维水平运动的微分方程，并对方程线性化得到了扩散方程，其解求得为：

Q_t＝Q₀e^-αt

式中，Q₀是地下径流初值，α是地下水的消退系数，是固定的经验值，模型中设为0.984，Q_t为t时刻的地下径流出流。

步骤C)完成水文过程概化后，需要进行流域水文过程模拟与验证。

对结果的精度验证采用如下五个统计指标作为模型精度的评价标准：平均偏差(MD)、均方根误差(RMSE)、确定性系数(NE)、相关系数(R²)和平均绝对偏差(MAE)。平均偏差MD反映了水量平衡模拟误差，均方根误差RMSE反映了模拟值和观测值序列之间误差的标准差，确定性系数NE反映了模型对观测变异性的描述能力，相关系数R²反映了径流模拟和观测序列间相关性，MAE反映了模拟的绝对值偏差，同时使用这些评价指数可以综合反映径流或土壤湿度模拟的精度以及与观测值之间的相关紧密程度。五个统计参数的表达式具体如下：

$M_D=n^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(O_i-M_i)$

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(O_i-M_i)}^2}{n}}$

$R^2=\sqrt{\frac{{\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(M_i-\stackrel{\‾}{M})(O_i-\stackrel{\‾}{O})\right]}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(M_i-\stackrel{\‾}{M})}^2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(O_i-\stackrel{\‾}{O})}^2}}$

$\mathrm{NE}=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(O_i-M_i)}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(O_i-\stackrel{\‾}{O})}^2}$

$\mathrm{MAE}=n^{-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}|O_i-M_i|$

式中，M_i为模拟径流量值，O_i为水文站观测的径流量值，

为观测径流量的均值，为模拟径流量的均值，n为时间序列的长度，需注意的是NE系数类似于R²，当模拟完全吻合实测时NE系数等于1，但与相关系数不同的是NE系数值可以小于0。由于MAE和M_D指标的意义较为接近，而M_D适合于反映水量平衡的整体模拟效果，因此在日径流模拟精度评价中，主要采用前四个统计指标。MAE指标将在土壤湿度的模拟精度评价中采用。

步骤D)构建基于分布式水文模型和粒子滤波同化算法的流域土壤湿度遥感数据同化方案，主要包括模型算子、观测算子、驱动气象数据以及同化算法四个部分。

粒子滤波同化算法

在Bayes框架下，最优滤波就是基于所有量测信息构造状态的后验概率分布函数(Posterior Distribution Function，PDF)，系统状态的各种估计值，如均值、协方差等都可从PDF获得。然而对大量实际中出现的非线性非高斯模型，无法求得PDF的解析解。为了解决该问题，出现了各种逼近方法，如扩展Kalman滤波(Extened Kalman Filter，EKF)，高斯混合滤波(Gaussian Mixture Filter，GMF)、交互式多模型(Interacting MultipleModel，IMM)等次优滤波方法，但算法的稳定性、精度等问题很难满足实际要求。

1954年，Hammersley和Morton提出了一种被称为“序贯重要性采样(SIS)”的蒙特卡洛(Monte Carlo methods，MC)方法，通过离散的随机测度逼近概率分布，并被应用到物理和工程领域，但由于高度的计算复杂性和退化问题，很长的一段时间内SIS算法没有太大的进展。为了解决退化问题，Gordon提出了重采样(Resamping)概念，由此实现了第一个可操作的Monte Carlo滤波，这些MC滤波在不同的领域的名称有Bootstrap，适者生存(Survival of The Fittest)，凝聚算法(Condensation)，序贯MonteCarlo法(Sequential Monte Carlo)等，现通称为粒子滤波。

粒子滤波是用于非线性非高斯状态空间模型最优估计问题的算法，其具体思想的数学描述如下：对于平稳的随机过程，假定k-1时刻系统的后验概率密度为p(x_k-1|z_1:k-1)，依据一定原则选取n个随机样本点，k时刻获得观测信息后，经过状态和时间更新过程，n个粒子的后验概率密度可近似为p(x_k|z_1:k)，达到最优贝叶斯估计的效果。PF算法的核心是基于MC方法，通过一组在状态空间中传播的粒子群对状态概率密度函数进行近似，获得状态的最小方差，因此PF也称顺序MC方法。尽管算法中的概率分布只是真实分布的一种近似，但由于非参数化的特点，PF可较精确地表达基于观测量和控制量的后验概率分布。对于线性高斯状态空间模型，最优滤波方法是Kalman滤波。

粒子滤波技术在非线性、非高斯系统表现出来的优越性，决定了它的应用范围非常广泛，目前国际上，PF算法已成功应用于信号处理、自动控制、金融、无线通讯等领域，有许多会议和讨论组都将粒子滤波作为专题进行深入讨论和学术交流。PF算法摆脱了解决非线性滤波问题时随机量必须满足高斯分布的制约，能表达比高斯模型更广泛的分布，也对变量参数的非线性特性有更强的建模能力，这也是PF优于EnKF最重要的一方面。由于流域水文系统及土壤水分运动具有高度复杂性和非线性特征，因此基于分布式水文模型的土壤湿度同化研究必须使用PF算法，以保证同化试验的合理性，这也是采用粒子滤波作为同化算子的原因。

动态系统的状态空间方程可表述为

x_k＝f_k(x_k-1，v_k-1)

z_k＝h_k(x_k，u_k)

x_k表示系统状态，本研究中为水文模型模拟的表层土壤湿度，z_k表示观测，为由MODIS和ASAR遥感数据估算得到的表层土壤湿度，u_k，v_k为独立的观测噪声和系统噪声。假设x_k服从一阶Markov过程，k-1时刻的后验分布函数p(x_k-1|z_1：k-1)已知，时刻k的模型状态PDF可以通过下式求得，

p(x_k|z_1:k-1)＝∫p(x_k|x_k-1)p(x_k-1|z_1:k-1)dx_k-1

p(x_k|x_k-1)由系统方程和已知的系统噪声v_k所定义，在k时刻观测量z_k即遥感估算结果获取后，通过贝叶斯方程对系统的状态进行更新，

$p\left(x_k\right|z_{1:k})=\frac{p\left(z_k\right|x_k)p\left(x_k\right|z_{1:k-1})}{p\left(z_k\right|z_{1:k-1})}$

其中归一化常数p(z_k|z_1:k-1)可由下式给出，

p(z_k|z_1:k-1)＝∫p(z_k|x_k)p(x_k|z_1:k-1)dx_k

式中p(z_k|x_k)由观测方程和已知的观测噪声u_k定义。以上两式构成最优的Bayes解，但其解析解只对有限的模型成立，EKF、GMF和IMM等逼近方法的估计精度有限，粒子滤波基于Monte Carlo思想，采用重要性采样的方法(Importance Sampling，IS)从重要性分布函数q(x_0:k|z_1:k-1)中独立抽取样本

模型状态的PDF逼近为

$p\left(x_{0:k}\right|z_{1:k})=\underset{i=1}{\overset{N_s}{\mathrm{\Σ}}}{w}_k^i\mathrm{\δ}(x_{0:k}-x_{0:k}^i)$

其中权重需归一化处理，即满足

为递推估计，选取重要性函数

q(x_0:k|z_1:k)＝q(x_0:k-1|z_1:k-1)q(x_k|x_0:k-1，z_1:k)

因此从q(x_k|x_0:k-1，z_1:k)中抽取样本

其权重

计算方法如下，本研究中，上式的右项即是模型模拟和遥感估算土壤湿度结果的误差比率，

$w_k^i=w_{k-1}^i\frac{p\left(z_k\right|x_k^i)p\left(x_k^i\right|x_{k-1}^i)}{q\left(x_k^i\right|x_{k-1}^i,z_k)}$

以上公式组成了顺序IS(Sequential IS，SIS)方法，针对迭代之后出现的权退化现象，Gordon提出对样本重新采样，繁殖重要性权重高的粒子，淘汰权重低的粒子从而抑制退化现象，最常用的重采样方法是SIR(Sampling Importance Resampling)。一般粒子滤波就由SIS和重采样法构成。

粒子滤波算法的计算流程如下，1.对模型的模拟得到的土壤湿度进行扰动，产生粒子群，即从q(x_k|x_0:k-1，z_1:k )中随机抽取N个有限样本，N为粒子个数；2.在遥感数据的获取时刻，使用遥感估算得到的土壤湿度结果，基于

计算对应粒子的权重；3.对粒子的权重作归一化处理；4.重采样，利用重要性采样SIR法对粒子权重进行重新采样；5.在权重采样后，结合模型状态的PDF逼近对模型状态进行最后的更新计算，得到同化后的土壤湿度结果。

本发明与现有技术相比，其创新主要体现如下：

(1)通过土壤参数库建立水文模型来统一刻画产流和土壤水分传输等水文过程，模型利用流域土壤类型差异描述了流域蓄水容量的空间异质性即实际空间格局，并由此参数动态计算流域栅格的缺水深获得流域产流量，克服传统模型需较多调节参数的缺点；

(2)在水文模型框架下对土壤剖面水分传输过程进行较为详细的物理参数化方案描述，解决了水文模型无法直接作为遥感土壤湿度同化研究模型算子的问题，可以得到合理的径流和具有稳定精度的表层土壤湿度模拟结果。

(3)基于粒子滤波的流域表层土壤湿度数据同化方案可以有效的融合遥感湿度信息，并提高湿度模拟精度；

(4)在空间分布方面，同化结果较好的融合了遥感估算和模型模拟结果的空间特征，与模型模拟相比土壤湿度空间差异性更加显著并表现出更为细致合理的空间纹理。

有益效果：

1)针对现有水文模型中的土壤水模拟这一薄弱环节，解决了水文模型无法给出表层土壤湿度量化信息，而无法直接作为遥感土壤湿度同化研究的模型算子的问题。利用土壤水动力学方法结合蓄满产流原理，构建了一个新的能够有效融合微波遥感信息、具有一定物理基础的分布式流域水文模型，经过典型半干旱半湿润区沂河流域的水文模拟检验，结果表明，逐日径流模拟效果较好、表层土壤湿度模拟精度具有稳定性，可作为流域土壤湿度遥感数据同化的模型算子。

2)建立了基于分布式水文模型的流域尺度土壤湿度同化方案，利用对地观测卫星TERRA-MODIS和ENVISAT-ASAR的土壤湿度估算结果对水文模型中的表层土壤湿度变量进行了四维数据同化。通过在半干旱半湿润区沂河流域连续水文模拟和粒子滤波同化试验的研究表明，同化结果有效融合了遥感湿度的空间格局信息，并提高了水文模型的土壤湿度模拟精度，该方案可有效的用于时空分布的流域尺度土壤湿度同化数据集的获取。

说明书附图

图1为本发明的流程示意图；

图2为遥感土壤湿度同化的观测算子构建流程；

图3为研究区内不同格点(以实测土壤湿度观测点为例)的理想试验同化结果；

图4为表层土壤湿度模拟和同化结果的实测验证(DOY：2006148)；

图5为不同日期的沂河流域LST-NDVI空间；

图6为水文同化模型的构建思路图。

具体实施方式

以下通过实施例对本发明作进一步描述：

实施例

为得到更加准确的流域表层土壤湿度数据，融合模型模拟和遥感数据反演二者的优势，以沂河流域为研究区，建立了流域土壤湿度同化方案。在利用已建立的分布式流域水文模拟平台对研究区2006年的土壤水分时空动态数值模拟的基础上，结合粒子滤波顺序同化算法，考虑了模拟和遥感反演的误差，同化主动微波ENVISAT-ASAR和MODIS可见光影像反演得到的地表土壤水分信息。顺序同化中得到更新的流域土壤湿度同化数据，将其反馈至分布式水文模型平台中，并通过递推模拟最终得到具有物理一致性的流域表层土壤湿度同化时空数据集。实测验证结果表明，同化结果精度要高于模拟结果的精度，在空间分布上，同化后的表层土壤湿度较好的融合了遥感估算和模型模拟结果的空间特征，与模型模拟值相比土壤湿度的空间异质性增强，表达更为合理

土壤湿度反演包括可见光遥感反演和主动微波遥感反演。

可见光遥感反演：

使用了2006年全年可以获取的无云覆盖MODIS影像反演沂河流域的表层土壤湿度，并利用地面观测数据对反演结果进行验证。首先，分别建立了流域四个季度的LST-NDVI空间以及相关干边和湿边关系，用于TVDI计算。其次利用流域内分布的站点观测数据建立了TVDI和实测表层湿度的线性关系，并据此关系从2006年的MODIS连续影像序列中反演得到了流域表层土壤湿度。最后我们分析了MODIS估算得到土壤湿度空间格局，并通过4个固定段观测站点的实测数据对反演得到的2006年连续土壤湿度估算结果进行验证。

在对连续的MODIS影像进行预处理后，利用热红外和其他相关波段基于劈窗算法反演了地表土壤温度，并同时计算了NDVI指数，得到研究区的LST-NDVI空间散点图，见附图5。

土壤湿度反演结果，为了检测流域内表层土壤湿度空间格局的逐日变化，选择了时间上连续的表层土壤湿度估算结果，其时间范围由2006年DOY260日到DOY267日，由于受到云的影响，DOY 260日到DOY 263日没有进行反演计算。可以较为明显的看出，由于局部降水作用流域表层土壤由干旱逐渐湿润再由湿润转化为干旱的过程，而MODIS估算结果也合理描述了该过程中的土壤湿度空间差异性，可以看出流域内的土壤湿度增高和降低幅值在空间上并不同步，其变率受到降水和土壤特性因素较为显著的影响。

主动微波遥感反演：

利用2006年获取的14景研究区ASAR WS影像反演得到了流域表层土壤湿度结果，得出了2006年DOY148日和DOY164日的反演结果。可以看出，土壤湿度值的空间分布并没有明显的地域特征，具有一定的随机性。与DOY14日的土壤湿度结果相比，在DOY164日研究区大部分范围内的土壤湿度有一定升高，可以看出，ASAR的后向散射系数可较好的描述和再现表层土壤湿度变化的空间差异性。

采用ASAR主动微波数据和MODIS影像，分别反演了沂河流域表层土壤湿度序列，并基于这一遥感估算结果构建了流域土壤湿度同化的观测算子。将流域内土壤湿度实际观测数据用于遥感估算结果的精度验证，给出了观测算子误差的定量评估，为进一步与模型数值模拟结果对比以及土壤湿度的顺序数据同化研究打下基础。土壤湿度反演的主要工作和结论分以下两点阐述，

(1)基于改进的ERS归一化方法，利用获取的ENVISAT ASAR影像反演了研究区的表层土壤湿度时间序列，并利用已有的土壤湿度观测资料对反演结果进行了验证，给出了量化评估结果，结果表明，ASAR影像估算的均方根估算误差(RMSE)为0.063m³/m³，平均偏差(MAE)为0.053m³/m³。这一估算精度较为合理，达到较成熟的ERS土壤湿度估算精度范围(0.03m³/m³-0.07m³/m³)。另一方面，从空间结果看，遥感估算结果空间分布格局总体合理，DOY148日土壤湿度估算结果的直方图分布总体类似于高斯分布，流域土壤湿度平均值为0.237m3/m3，大部分数值位于0.200m³/m³和0.330m³/m³之间，与实际的湿度观测经验相符，结果合理。

(2)基于TVDI方法，利用2006年可获取的MODIS影像反演了沂河流域的表层土壤湿度，实际观测资料的验证表明，MODIS估算结果与实测数据具有良好的正相关性。从四个观测站点的验证结果来看，均方根误差水平较低，数值范围在0.049m³/m³到0.102m³/m³之间，平均误差值为0.072m³/m³，总体精度令人满意。为了检测流域内表层土壤湿度空间格局的逐日变化，利用反演结果给出了连续多日的表层土壤湿度分布，可以看出，MODIS的土壤湿度估算结果可较清晰地反映局部降水作用下流域表层土壤的干湿转化过程，并合理描述了土壤湿度的空间异质性。

总的来说，ASAR和MODIS的土壤湿度反演结果总体满意，可有效反映表层土壤湿度变化的空间差异性。本章根据研究区实测数据对ASAR土壤湿度估算模型的经验参数做了一定的修正，给出了适合研究区的ASAR湿度反演参数以及MODIS影像TVDI反演方法的经验参数。结合遥感估算方法和水文模型的模拟方法研究流域土壤湿度，对进一步揭示土壤水分变化规律有重要意义。而未来完善的流域观测系统将有助于遥感土壤湿度反演的进一步发展。

在水文模型概化的基础上，构建分布式水文模型同化平台。在流域气象驱动参数、土壤参数和土地利用等输入数据的支持下，将模型应用于江口站控制流域(湿润区流域)和沂河流域(半干旱半湿润区流域)的水文过程模拟，使用流域出口水文站的实际观测逐日流量数据对模型的模拟精度进行检验和评价，并与常见水文模型进行了比较验证，最后从多个统计评价指标对模型的模拟结果进行了检验和分析。为了验证模型的流域表层土壤湿度的模拟能力，利用沂河流域获取的2006年多个站点土壤湿度观测资料，对10cm和20cm深度的土壤湿度模拟结果进行了验证分析。

逐日径流模拟及评价：

为检验模型的模拟能力，首先选择汉江上游的江口流域进行模型的逐日径流模拟测试，其中1981年～1983年为模型模拟的率定期，1984年～1985年作为模型验证期。

需要说明的是ESSI模型在不同的分辨率下模拟结果有一定的差异，总体上模拟分辨率较高时，模型的模拟结果较好，如240m栅格情形下的模拟结果要好于1000m的模拟结果。

模型与ESSI模拟结果进行比较表明，新模型对流量过程的拟合逼近程度要明显高于ESSI模型1000m栅格分辨率的模拟结果，同时模型模拟的确定性系数也量化反映了二者在模拟效果上的差别。

为进一步检验模型的模拟能力，选择沂河临沂以上流域进行模型的逐日径流模拟测试，其中2001～2003年为模型模拟的率定期，2006～2007年作为模型验证期。模拟结果所使用的评价指标同上，与上节江口流域模拟类似，新模型的模拟效果和拟合逼近程度要显著好于ESSI模型1000m栅格分辨率的模拟结果。

土壤湿度模拟及评价：

为了验证模型的流域表层土壤湿度模拟能力，利用在沂河流域获取的2006年多个站点土壤湿度观测资料，对2006年模型模拟得到的10cm和20cm深度的流域土壤湿度结果进行了验证分析。

表层土壤(10cm)湿度模拟：

表层土壤湿度模拟值与实测湿度数据具有较好的相关性，尤其是在固定地段的观测点。其中平邑固定地段模拟值的变化趋势和数值大小均与观测值序列十分接近，表明土壤水分模型的纵向模拟有一定可靠性，在数据精度满足条件下，可达到较高的模拟精度，这也为同化试验生成具有物理一致性的连续土壤湿度数据集保证了基本精度和可靠性(验证精度见下表)。总体来看，相对于实测数据，模型模拟值存在一定程度的低估趋势，这主要是由于对土壤层水分的提取项计算过大造成土壤层的失水过快所造成，而模型中与蒸腾相关的根系吸收水分作用模拟偏大是首要因素，因此植物蒸腾量以及分层根系水分提取的计算仍需进一步改进。在某些土壤中间层的水分低估趋势较为严重，表明格点上剖面根系分布参数不合理。通过同化独立的遥感土壤湿度信息将有助于改善模拟低估的情况。

土壤水分(10cm)模拟误差统计结果

上层土壤(20cm)湿度模拟：

上层土壤湿度模拟结果总体上与表层土壤湿度模拟结果显示了较高的相似度。固定地段的实际观测数据与模拟值的趋势拟合较好，虽然在整体上有一定偏低估计的趋势，该层土壤湿度的低估原因与10cm土壤层基本一致。对于非固定地段，除了在沂水丘陵和平原站点，部分时间的模拟值接近于0，与实测数据偏差较大之外，其他站点的模拟数值与观测数据较为接近(验证精度见下表)。沂水丘陵和平原站点的部分模拟值过低而接近零的原因主要在于格点上纵向根系分布参数赋值不均，过于集中在20cm土壤层，致使无雨期土壤湿度模拟值下降较快，并达到不合理低点。

虽然存在一定的不足，但从模拟验证总体结果看，可以预见，从模型数值模拟的角度入手研究，获取一定精度的表层土壤湿度信息是可行的。进一步完善和改进水文模型的不同类型输入数据和下垫面参数的精度，将有助于流域表层土壤湿度动态模拟的研究。

土壤水分(20cm)模拟误差统计结果

总的来说，流域内所有可以利用的观测验证数据分布于8个不同站点(沂水麦田，平邑麦田，临沂麦田，平邑平原，沂水洼地，沂水平原，沂水山地，沂水丘陵)。总的来说，模型模拟结果精度评价的一个重要结论就是各点均方根误差RMSE稳定在0.080m³/m³左右。所有观测点的RMSE误差值范围在0.040～0.130m³/m³之间，平均RMSE误差是0.083m³/m³，这一结果表明通过水文模型得到的表层湿度模拟值具有可靠性。不仅对于固定点，对于其他非固定观测点，模拟结果也表现出相似的误差大小。

将遥感观测数据反演得到的表层土壤湿度作为理想资料，利用粒子滤波进行连续的同化试验，验证粒子滤波用于复杂非线性分布式水文模型遥感数据同化的可行性。理想试验以沂河流域为研究区，模型模拟及同化时段为2006年，模拟空间分辨率为1000m。试验中，土壤同化深度为表层10cm，观测数据为ASAR和MODIS、数据的湿度反演结果，模型的湿度模拟时间步长为30分钟。模型只在遥感数据可以获取的当天进行同化计算和更新，同化后得到的土壤湿度作为流域水文状态变量当天的准确估计值，代入模型中用于下一时刻的模型模拟。结果如附图3。

附图3给出了2006年全年所有同化时刻通过模型、遥感估算以及同化得到的表层10cm土壤湿度结果。图中横坐标代表同化的天数，按照遥感影像获取的日期顺序排列，并以0和0.60为阈值剔除了遥感资料中的少数异常点引起的土壤湿度极值。可以明显看出，8个格点的同化效果均十分显著。从趋势上看，土壤湿度的同化结果与理想“真值”吻合程度很高，而主要的差异部分在于，遥感观测的理想值发生突变的位置即土壤湿度值过高超出正常的范围时，同化结果较接近模型模拟值，显示了算法滤波的特性，这一特性较好滤除遥感噪声影响，维持表层土壤湿度同化结果的相对稳定。以上结果说明PF用于表层土壤湿度同化试验有较好的可靠性。

为了进一步对比验证同化后湿度的空间分布格局，从给出的同化结果可以知道，首先，同化后的土壤湿度空间格局较为显著的融合了遥感信息和模型模拟信息，同化结果中可见遥感估算中未有的湿度信息如模拟结果中出现的河道，但同化结果仍表现出与遥感估算结果较一致的空间分布。其次，同化结果在结合模拟信息的同时，改善了模型模拟结果片状分布的缺点，通过融合遥感土壤湿度信息，具有更高程度的土壤湿度空间差异以及更精细的纹理特征。总体上来说，同化结果在空间分布上较为趋近于遥感估算结果，尤其是ASAR影像的湿度估算结果。以上分析结果进一步支持了流域尺度土壤湿度空间同化的可行性以及粒子滤波同化方案的适用性。

为了更直观的了解同化后的土壤湿度结果，尤其是同化对模型模拟值的改进效果，下表给出了各测点的同化以及模拟误差统计结果。

模拟与同化结果的总体统计评价(m3/m3)

这里列出流域内8个观测站点2006年全年表层土壤湿度同化结果的精度评价以及和模型模拟精度的对比，可以看出，与模型模拟结果相比，同化后的土壤湿度结果有较为明显的改进，8个观测站点中6个站点的均方根误差减小，6个站点平均的RMSE误差降低0.020m³/m³，所有站点总体平均误差降低约0.010m³/m³。其中沂水平原站和沂水麦田站的改进较为明显，分别达到0.041m³/m³和0.029m³/m³，这主要是由于同化的土壤湿度值中所融合的遥感反演信息缓解了模型模拟的低估趋势。

获取流域尺度表层土壤湿度的模型数值模拟法和遥感估算法各有利弊，二者精度和误差结构亦不相同。模型数值模拟结果具有物理一致性好和时间连续的优点，但空间精度受到输入数据精度的影响较大，而遥感估算方法在土壤湿度空间格局上有较好的捕获信息的能力。为融合两种方法的优势，有必要采用数据同化方法改进土壤湿度的计算。

基于此，在对模型和遥感估算误差的一定先验认识的基础上，以水文模型为平台建立了流域尺度土壤湿度同化框架，针对水文模型的高度复杂性和非线性特征，以非线性贝叶斯估计的粒子滤波方法作为同化算子，观测算子是通过ASAR和MODIS遥感数据估算得到的土壤湿度。具体的同化策略如下，在表层土壤湿度数值模拟过程中，当遥感观测数据获取时，采用粒子滤波同化ASAR或MODIS估算得到土壤湿度结果，并将同化后的结果代回分布式水文模型继续计算，纠正模型数值模拟的运行轨迹，如此迭代循环最终得到同化后的流域土壤湿度时空数据集。

同化试验选择沂河流域为研究区，2006年全年为模拟同化时段，在流域气象驱动数据、土壤参数、地表覆盖等数据支持下，利用分布式水文模型共同化遥感影像149景，其中ASAR影像14景，MODIS影像135景。在同化试验前，首先对遥感估算和模型模拟结果作了比较分析，后利用同化平台分别进行了同化理想试验、一般性同化数值试验和敏感性试验，最后利用所有可以获取的站点实测数据对同化结果进行了验证分析。

结果表明，同化后的表层土壤湿度结果相对于模型模拟结果有明显改善。在空间分布方面，同化结果较好的融合了遥感估算和模型模拟结果的空间特征，与模型模拟相比土壤湿度空间差异性更加显著并表现出更为细致合理的空间纹理。进一步通过流域空间内8个观测站点连续时间资料的验证，与模拟结果相比，6个站点的均方根误差减小，6个站点平均RMSE误差降低为0.020m³/m³，所有站点总体平均误差降低0.010m³/m³。在表层土壤湿度空间误差方面，同化结果精度要高于模拟结果，以ASAR遥感影像和观测数据同时可获取的2006年DOY148日为例，流域内模型模拟RMSE误差为0.081m³/m³，同化后RMSE误差降低至0.045m³/m³。

与遥感估算的土壤湿度值相比，虽然模型误差较为稳定，但模型模拟值有较明显低估趋势而且波动性相对较小，同化后的土壤湿度结果能较好纠正这一低估趋势，使土壤湿度变化更加合理。理想试验结果同时表明，粒子滤波算法在表层土壤湿度同化中具有很好的同化能力，同化结果与理想资料的拟合趋势一致，精度较高，并能够去除遥感估算中过高的异常噪声。从同化敏感性试验来看，同化土壤湿度验证结果对于观测均方差参数并不敏感，在0.030～0.090m³/m³范围内对同化结果精度影响不大，同化算子相对稳定。

Claims (8)

1.一种流域尺度土壤湿度遥感数据同化方法，其步骤为：

A)数据准备；

B)土壤湿度反演；

C)分布式水文模型同化平台构建；

D)流域尺度土壤湿度同化。

2.根据权利要求1所述的流域尺度土壤湿度遥感数据同化方法，其特征在于步骤A)中，所述数据准备包括水文、气象数据，DEM数据，植被、土壤参数库，流域土地利用/覆被和土壤类型下垫面参数数据，卫星影像数据以及逐日径流和土壤湿度实测数据。

3.根据权利要求2所述的流域尺度土壤湿度遥感数据同化方法，其特征在于步骤B)中使用MODIS遥感和主动微波ASAR遥感构建同化的观测算子。

4.根据权利要求3所述的流域尺度土壤湿度遥感数据同化方法，MODIS遥感使用的温度植被干旱指数TVDI由以下公式给出：

LST为地表温度；LST_min为相同NDVI值的最小地表温度，对应的是LST-NDVI特征空间的湿边；LST_max为相同NDVI值的最大地表温度，对应的是LST-NDVI特征空间的干边。

5.根据权利要求3或4所述的流域尺度土壤湿度遥感数据同化方法，ASAR遥感反演需要的使用以下方法和公式：

ASAR数据预处理

辐射校正

对于ENVISAT ASAR数据，影像DN值和后向散射系数之间的关系式如下

其中σ⁰是线性单位的雷达后向散射系数，值域在0～1之间，k是外部校正参数，可以从ASAR影像的头文件中得到，θ是空间分布的入射角；分析数据时常用以分贝为单位的σ⁰，因此需要对上式得到以线性单位表示的σ⁰进行对数转换，

σ⁰(dB)＝10[log(σ⁰)] 

几何校正

需要对获取的属于一级产品的ASAR影像进行几何校正，这主要分为两个步骤，首先是几何粗校正，其次是与数字高程DEM影像数据匹配的几何精确校正；

ESA的雷达处理影像工具NEST软件提供了ASAR影像的几何校正功能，在利用该工具对ASAR数据进行粗校正之后，进一步通过与DEM的比对，以DEM为基准影像，进行图像对图像的精确校准，如此循环后将所有的ASAR数据统一的精确订正到DEM的坐标系下；

光斑滤波

光斑噪声在图像上看起来类似于椒盐状分布，使图像的可解译性变差；有效的降低噪声对于计算雷达影像的后向散射系数十分关键，目前去除噪声的方法主要有多视方法和空间滤波方法；

详细的ASAR数据湿度估算方法如下：

散射系数归一化

IEM辐射传输模型的数值模拟试验也验证了这一关系的存在，地表后向散射与入射角关系可以用下式表达，

σ⁰＝pθ+q

此处，θ表示入射角，p、q是与土地覆盖有关的经验参数，具体数值可由Loew etal.中查到；尽管每一景ASAR WS模式影像的入射角各不相同，同一景影像中的不同区域入射角变化范围也较大，可通过下面公式，将所有影像的后向散射系数归一化至参考的入射角θ₀，此处θ₀为ERS的入射角度值23度，具体的归一化方法如下，

其中，σ⁰为图像原始后向散射系数， 

为归一化后的图像后向散射系数；

反演参数

中利用上一步归一化后处理得到的ASAR后向散射，结合下式可以得到土壤的介电常数，最终反演得到表层的土壤湿度值，

∈_r＝a+bσ⁰+c(σ⁰)²

其中∈_r为介电常数，a、b、c为经验参数，与土地覆盖类型相关，介电模型 

∈_r＝(a₀+a₁S+a₂C)+(b₀+b₁S+b₂C)m_v+(c₀+c₁S+c₂C)m_v ²

其中，m_v是需要计算的体积含水量，S和C是土壤的砂粒和粘粒的重量百分比，即土壤质地参数，可以由我们建立的流域土壤物理参数库得到，a_i，b_i和c_i是模型的经验常数，已由文献给出；

基于上述模型，利用获取的流域土壤介电常数图像和已有的土壤质地参数，通过迭代计算方法(误差阈值为0.01)逐栅格推求得到最终的ASAR遥感数据的表层土壤湿度估算结果。

6.根据权利要求3或4所述的流域尺度土壤湿度遥感数据同化方法，水文过程概化方法如下：

林冠截留量计算公式如下：

其中Cm为冠层截留能力，k为降雨截留系数，P为降雨量(mm)，T为降雨历时(hr)，ew为林冠蒸发率，因此ew·T表示树冠蒸发以及树干蒸发之和；

冠层截留能力Cm采用Von Hoyningen-Huene所提出的方法计算：

C_m＝C_p·(0.935+0.498·LAI-0.00575·LAI²)

其中Cp为植被盖度，LAI为格点单元平均的逐日叶面积指数

降雨截留系数k是植被叶面积指数的函数，由下式给出：

k＝0.046·LAI

林冠蒸发率e_w与降雨初期的干燥程度和风速等因素相关，在无实测数据验证时，通常取以下经验参数：亚热带地区取0.031～0.047，温带湿润、半湿润地区取0.063～0.093，温带干旱、半干旱地区取0.170～0.200；

水量平衡计算公式如下：

P-Es-Et-G＝Δw+Δc+R

其中，P是降水(包括降雨和降雪)，Es是土壤表面蒸发，Et是植物蒸腾，G是地下水出流，Δw是土壤水的增量，Δc是冠层截留，R是地表径流；

土壤水分传输中涉及到不饱和土壤水通量定律Buckingham-Darcy方程和连续方 程，Buckingham-Darcy方程如下：

其中q是不饱和水流通量，h是土壤水基质势，K(h)是水力传导度，z是垂直方向的深度；

连续方程如下：

基于垂向流的假设，公式中的q由Buckingham-Darcy方程替换，求得解即为一维Richard’s方程表达式：

其中t是时间，θ是土壤体积含水量；

根系土壤水分提取，此土壤水含量计算中需考虑根系吸水作用，本模型将其作为S项添加到Richard’s方程中；本发明使用以下方法考虑不同土壤层的根系水分吸收作用：

其中S表示土壤根系吸水或者土壤表面蒸发，由以下公式描述

其中i是模拟土壤层的序号(i≥2)，对于第一层土壤(i＝1)，Si为表层土壤表面蒸发，Eveg是植物的蒸腾，ri是土壤层中根系分布参数，与植被类型有关，wi是土壤水饱和系数，当湿度处于凋萎点时等于0，完全饱和时为1，由下式给出：

其中 

是凋萎点的水势，等于-1.5×105， 

是饱和土壤水势， 

是土壤水基质势；

潜在蒸散发，当流域存在蒸散发观测站时，潜在蒸散发计算由观测数据插值得到；当流域无蒸散发的观测时，利用气象观测数据根据Penman-Monteith公式计算得出， Penman-Monteith公式由下式给出：

其中E是蒸发量，λ是水汽潜热，Δ是饱和水汽压温度曲线的斜率，A是净辐射强度，ρa是湿空气的密度，cp是常压下的空气汽化潜热，D是水汽压与饱和水汽压的差值，ra 是空气动力学阻力系数，γ是干湿计常数，rs是作物叶面的阻力系数；

植物蒸腾，蒸腾量计算如下：

0＜LAI≤3.0

LAI＞3.0

其中Eveg是植物的蒸腾量，E是潜在蒸散发量，p是格点上的植被覆盖率，LAI是格点上的逐日叶面积指数；

植被覆盖率p通过逐日叶面积指数LAI来计算，公式如下：

p＝1-e^-α·LAI

式中：α为地表太阳辐射的消光系数，缺省值为0.45；

土壤表明蒸发，分别考虑土壤水分蒸发的两个主要阶段，即当土壤水分供给充足，蒸发量较高且相对稳定；当土壤水分降低，土表蒸发随着土壤水分含量呈非线性下降：

θ＜f_c

θ≥f_c

其中Esoi是土壤表面蒸发量，E是潜在蒸散发，dz是表层土壤的厚度(mm)，θ是表层土壤的水分含量(体积含水率)，fc是表层土壤的田间持水量(体积含水率)，wp是表层土壤的凋萎系数(体积含水率)；

地下水与不饱和土壤水交互，地下水的补给计算由以下公式给出：

其中Q是地下水补给量(mm/s)，当水流向下流向地下含水层时为正值，反之为负值， 

是地下水位的深度， 

是底层土壤水基质势，z bot是底层土壤的节点深度，K a是含水层的水力传导度，由下式给出：

其中f是延时因子，用于水力传导度的计算；

地下汇流，地下水运动是渗流水力学的问题，Boussinesq根据连续方程和达西定律，建立了饱和带中的一维水平运动的微分方程，并对方程线性化得到了扩散方程，其解求得为：

Q_t＝Q₀e^-αt

式中，Q0是地下径流初值，α是地下水的消退系数，是固定的经验值，模型中设为0.984，Qt为t时刻的地下径流出流。

7.根据权利要求1和6所述的流域尺度土壤湿度遥感数据同化方法，需要进行流域水文过程模拟与验证；

对结果的精度验证采用如下五个统计指标作为模型精度的评价标准：平均偏差(MD)、均方根误差(RMSE)、确定性系数(NE)、相关系数(R²)和平均绝对偏差(MAE)；五个统计参数的表达式具体如下：

式中，M_i为模拟径流量值，O_i为水文站观测的径流量值， 

为观测径流量的均值， 

为模拟径流量的均值，n为时间序列的长度，需注意的是NE系数类似于R²，当模拟完全吻合实测时NE系数等于1，但与相关系数不同的是NE系数值可以小于0；由于MAE和MD指标的意义较为接近，而MD适合于反映水量平衡的整体模拟效果，因此在日径流模拟精度评价中，主要采用前四个统计指标；MAE指标将在土壤湿度的模拟精度评价中采用。

8.根据权利要求1所述的流域尺度土壤湿度遥感数据同化方法，粒子滤波算法的计算流程如下，a对模型的模拟得到的土壤湿度进行扰动，产生粒子群，即从q(x_k|x_0:k-1，z_1:k)中随机抽取N个有限样本，N为粒子个数；b在遥感数据的获取时刻，使用遥感估算得到的土壤湿度结果，基于 

计算对应粒子的权重；c对粒子的权重作归一化处理；d重采样，利用重要性采样SIR法对粒子权重进行重新采样；e在权重采样后，结合模型状态的PDF逼近对模型状态进行最后的更新计算，得到同化后的土壤湿度结果。 

Images