CN104573393B

CN104573393B - 一种基于贝叶斯理论的土壤水分站点数据升尺度方法

Info

Publication number: CN104573393B
Application number: CN201510041606.4A
Authority: CN
Inventors: 高胜国; 潘耀忠; 朱忠礼; 朱秀芳
Original assignee: Beijing Normal University
Current assignee: Beijing Normal University
Priority date: 2015-01-28
Filing date: 2015-01-28
Publication date: 2017-10-13
Anticipated expiration: 2035-01-28
Also published as: CN104573393A

Abstract

本发明公开了一种基于贝叶斯理论的土壤水分站点数据升尺度方法，包括，基于升尺度区域中稀疏站点观测数据估计目标变量的先验概率密度分布pre_pdf；通过建立SM和MODIS ATI之间的非线性回归关系，将MODIS ATI反演为SM，并且估计土壤水分非线性回归的估计置信区间以及概率分布作为概率形式的软数据；通过贝叶斯公式融合目标变量的先验分布以及来自于MODIS ATI的概率形式的辅助信息，获取目标变量的后验概率密度分布函数post_pdf；通过后验概率密度函数的最大化求取概率最大时目标变量的取值。本发明有效减少土壤水分遥感产品和地面站点验证数据之间尺度差异而引入的不确定性，亦可应用到其他地表参数的升尺度应用当中。

Description

一种基于贝叶斯理论的土壤水分站点数据升尺度方法

技术领域

本发明涉及导航遥感领域，具体是一种基于贝叶斯理论的土壤水分站点数据升尺度方法。

背景技术

尺度问题是多源土壤水分数据被充分利用的一个巨大阻碍，尺度不匹配问题不仅存在于土壤水分的多途径观测过程中，而且存在于土壤水分模拟以及数据同化等各个相关领域。低分辨率土壤水分产品(SMAP为大约100km²，SMOS为大约1600km²，等)的地面验证工作也受到尺度因素的限制，这种制约来源于卫星传感器分辨率和地基仪器的点观测之间的尺度差异。另外，严格意义上讲点观测的尺度(support area小于1m²)也同样不匹配于模型格网的尺度(support area大于10m²)，不同尺度的土壤水分观测以及模拟结果之间的尺度转化对于现在以及未来星载传感器土壤水分产品的验证以及数据同化技术的推广应用至关重要。所以，发展尺度转换的策略和方法就显得尤为急迫，当然，尺度问题不仅仅存在于土壤水分相关的领域，同样是整个地球科学领域的关键问题。本发明主要以土壤水分为实施例，解决针对微波遥感土壤水分产品地面验证的升尺度问题。到目前为止，虽然升尺度策略已经有多种，但升尺度方面的研究仍处于初级阶段，还没有一种通用的升尺度方法可以满足应用的需求，而且现阶段大部分的大尺度微波遥感土壤水分产品的验证工作仍然是基于地面多点观测的平均值来进行，所以仍然有必要发展新的升尺度方法。另外，随着多源多尺度遥感数据的持续增加，发展新的升尺度方法时需要重视海量遥感数据中的宝贵信息。

贝叶斯理论在融合遥感信息中可以发挥重要作用。目前，还很少有人基于贝叶斯数据融合理论进行升尺度方法的构建。本发明的目的是将土壤水分的稀疏站点观测数据升尺度到微波遥感的像元尺度。

发明内容

本发明在贝叶斯理论框架下，通过融合多源信息，尤其是具有不确定性的信息(软数据)，构建了一种非线性的站点观测数据升尺度方法。本发明中涉及的多源数据包括土壤水分无线传感器网络数据和MODIS的表观热惯量(ATI)。

本发明采用如下技术方案：一种基于贝叶斯理论的土壤水分站点数据升尺度方法，包括如下步骤：

步骤一、基于升尺度区域中稀疏站点观测数据估计目标变量(区域尺度土壤水分)的先验概率密度分布(pre_pdf)；

步骤二、通过建立SM和MODIS ATI之间的非线性回归关系，将MODIS ATI反演为SM，并且估计土壤水分非线性回归的估计置信区间以及概率分布作为概率形式的软数据；

步骤三、通过贝叶斯公式融合目标变量的先验分布以及来自于MODIS ATI的概率形式的辅助信息，获取目标变量的后验概率密度分布函数(post_pdf)；

步骤四、通过验概率密度函数的最大化求取概率最大时目标变量的取值，即通过融合MODIS ATI信息将站点土壤水分观测升尺度后的土壤水分。

优选的，上述步骤一具体是利用目标尺度(例如，微波遥感像元足迹)范围内的稀疏站点观测数据估计目标尺度下目标变量的先验概率密度分布，在样本数据服从正态分布的假设下，可以利用基于样本数据的正态总体均值的区间估计，取95％置信区间内目标变量的概率分布来表示目标变量的先验概率密度分布，估计方法具体为正态整体均值的区间估计，估计公式为：

其中，f_G(μ)为总体的先验概率密度分布；为样本均值，S²为样本方差，T_A服从自由度为n-1时的学生t分布。而T_A的分布只依赖于样本个数n，并不依赖于总体均值μ和方差σ²；是T_A的概率密度分布。

优选的，上述步骤二中选择的土壤水分(SM)和表观热惯量(ATI)之间的非线性回归关系表达式为：

其中，a，b，c为系数。

优选的，上述步骤二中利用的辅助数据并非直接利用，而是首先将辅助数据表达为目标变量的软数据，例如本发明实施例中将表观热惯量表达为土壤水分的概率分布形式软数据。

优选的，上述步骤二中的非线性回归置信区间是基于渐进线理论进行估计，并且取95％的置信区间以及在此区间目标变量的概率密度分布来表达为目标变量的软数据；非线性回归的估计置信区间的估计利用公式为：

其中，F_S和f_S分别为SM₀的概率分布函数和概率密度函数；为在自由度n-p情况下t分布的概率密度函数；SM₀表示回归估计置信区间中土壤水分的可能取值；表示对应于指定表观热惯量ATI₀的土壤水分回归估计值；p表示模型参数个数；表示模型回归估计值误差的标准差。

优选的，所述步骤三还包括，需要构建块与块之间的空间协方差C_BB′，空间协方差C_BB′可以表示为：

其中，n是块B中子块的个数；n′是块B′中子块的个数；在这里，子块之间的协方差刚好是直接由协方差函数估计得到的协方差。

优选的，上述步骤三中目标变量后验概率的获取方法是基于基于贝叶斯原理，其中先验分布为基于站点数据估计的目标变量的先验分布，辅助信息为通过遥感数据获取的概率分布软数据。

与现有技术相比，本发明的优点：本发明基于贝叶斯理论的土壤水分站点数据升尺度方法，提供的升尺度方法可以有效减少微波遥感土壤水分产品地面验证工作中由于土壤水分遥感产品和地面站点验证数据之间尺度差异而引入的不确定性，同时本发明亦可应用到其他地表参数的升尺度应用当中。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1升尺度研究区概况及土壤水分无线传感器网络分布示意图；

图2贝叶斯框架下构建的升尺度方法示意图；

图3升尺度方法案例数据流程图；

图4升尺度范围内站点数据和辅助数据分布及编号示意图；

图5块与块之间协方差示意图；

图6大尺度范围内平均MODIS ATI时间序列以及土壤水分站点数据多点平均时间序列；

图7土壤水分和MODIS ATI之间非线性回归以及置信区间估计示意图；

图8MODIS ATI和六天平均MODIS ATI(前后各三天)相关性示意图；

图9试验协方差及理论协防差函数示意图；

图10基于M7，M16and M21三个站点的升尺度及验证结果；

图11基于M5，M13and M19三个站点的升尺度及验证结果；

图12三个站点所有随机组合情况下，贝叶斯升尺度方法和站点平均方法升尺度精度对比图；

图13贝叶斯方法和站点平均方法升尺度精度对利用站点数量的敏感性分析示意图；

图14是本发明实施例中贝叶斯理论框架下升尺度方法概念示意图.

具体实施方式

以下结合附图对本发明的实施例进行详细说明，应当理解，此处所描述的实施例仅用于说明和解释本发明，并不用于限定本发明。

研究区及数据

本发明实施例的研究区位于青藏高原的中部。如图1所示。研究区中生物种类较少，土壤水分变化幅度较大，并且季节性的冻融循环。在研究区存在一个密集的观测网络来监测土壤水分和土壤温度的变化，观测网络包括三个尺度上的嵌套网络(1.0，0.3，0.1degree)，每个网络节点共包括四层的观测(0-5，10，20，and 40cm)。观测网络在2010-2011年的土壤水分数据已经在国际土壤水分观测网络网站(International SoilMoisture Network)上发布。该网络的土壤水分和土壤温度的观测进行过严密的定标工作。在发明实施例中，我们使用了2010年8月1日-2011年12月31日之间大尺度网络和中尺度网络的0-5cm的土壤水分观测数据。这两个尺度的网络观测节点的分布如图1中所示。

MODIS标准产品由美国陆地过程分布式活动档案中心(Land ProcessesDistributed Active Archive Center，LPDAAC)发布，该中心位于美国地质调查局的地球资源观测和科学中心(Earth Resources Observation and Science Center(EROS)，U.S.Geological Survey(USGS)，https://lpdaac.usgs.gov)。在本发明实施例中，我们利用了V5.0的地表反照率和温度产品，分别为MCD43B3(Terra+Aqua 16-day albedo)、MOD11C1(Terra daily land surface temperature)和MYD11C1(Aqua daily landsurface temperature)，分别率均为1km。范围为100km*100km，基本上是大尺度网络的大概覆盖范围。

目标尺度土壤水分先验分布估计

在需要升尺度的区域中，一般情况下地面的土壤水分站点观测会相对稀疏。而像元尺度范围内稀疏的几个点仍然对区域土壤水分“真值”有一定的先验指示作用，首先，我们需要根据样本点来对区域土壤水分“真值”的先验概率密度分布进行估计。假设土壤水分观测点{SM₁，SM₂，...，SM_n}是均值为μ方差为σ²正态分布总体的随机样本，那么：

其中，为样本均值，S²为样本方差，根据概率论和数理统计原理，我们可以构建统计量：

T_A服从自由度为n-1时的学生t分布。而T_A的分布只以来于样本个数n，并不依赖于总体均值μ和方差σ²，当然总体均值和方差也是不可观测得到的。这样，我们可以计算得到总体的先验概率密度分布f_G(μ)为：

在本发明中，我们将像元足迹范围内少数的土壤水分站点观测定义为升尺度方法中的一般数据(general data)，计算得到的整体均值的概率密度分布f_G(μ)作为大尺度土壤水分的先验信息使用。

MODIS ATI的估计

土壤热惯量反映着土壤对温度变化的阻抗能力，它随着土壤含水量的变化而变化，因此它可以被用作土壤含水量大小的指标。正因如此，在稀疏或低矮植被覆盖下，我们可以将MODIS反演的地表表观惯量(ATI)用作土壤水分升尺度过程中的辅助数据。表观热惯量(ATI)定义为：

其中：

ATI是表观热惯量[K-¹]，α₀是表观反照率[-]，C是太阳校正系数[-]，A是地表温度日变化范围[K]，是地球纬度[rad]，δ太阳赤纬角[rad]。太阳赤纬角δ通过以下模型计算得到：

δ＝0.006918-0.399912cos(Γ)+0.070257sin(Γ)

-0.006758cos(2Γ)+0.000907sin(2Γ) (7)

-0.002697(3Γ)+0.00148sin(3Γ)

其中：

Γ是日角度变化量[rad]，n_d是天数[-]。另外，这里使用MODIS短波(0.3-5μm)反照率产品代替地表反照率α₀。日温度变化范围A也需要通过MODIS温度产品推算，在本研究中，日温度变化通过余弦函数来模拟：

其中，T(t_i)是在t_i时刻的地表温度[K]，是日平均温度[K]，A是温度日变化范围[K]，ω是地球自转角速度[rad s^-1]，t_i是一天中的时间[s]，ψ是相位[rad]，可以容易计算得到ψ。那么，剩下的两个未知变量为日温度变化范围A和日平均温度如果每个像元每天可以有两次或更多次的温度测量值(n，正常情况下MODIS每天有四次)，可以根据最小二乘法计算得到，经演算，A和可以表示为：

当n＝2时，日温度变化范围可以表示为：

这样，当一天中最少有两次的MODIS温度时就可以计算得到当天的地表表观热惯量。不过当只有两次温度时，建议两次温度中一次是白天温度，一次是晚上温度，如果两次都是白天或是晚上，会导致公式(12)的分子分母都很小。由于MODIS地表温度数据会由于云覆盖等影响而缺失，所以在某些天会由于MODIS地表温度数据的缺失而影响到ATI的获取。另外，在地表为浓密植被覆盖的情况下，MODIS ATI更多反映的是植被冠层的信息，基本不能反映地表水分信息，此时，ATI也不建议作为辅助信息加入到升尺度过程中。

基于ATI的土壤水分反演及置信区间估计

根据公式(5)，MODSIATI是和日温度变化幅度成正比的。由于水的比热容比较大，在土壤质地和土壤吸收热量均相同的条件下，土壤含水量越高，温度变化越小，那么ATI也就越大，这样，地表ATI可以作为地表土壤含水量多少的指标。在此部分中，我们首先需要通过建立SM和ATI的经验关系来根据表观热惯量反演土壤水分。在本研究中我们使用对数模型来描述SM和ATI之间的关系：

在模型中，表示地面土壤水分多点观测的平均值，表示对应于所有的土壤水分地面观测站点的ATI像元平均值，a，b和c是回归系数。

实际上，根据回归关系来反演土壤水分会有一定的不确定性。在升尺度过程中，我们以不确定性数据来对待基于非线性回归得到的土壤水分数据。可行的方案是我们可以根据非线性拟合的估计置信区间将将估计得到的土壤水分表达为概率形式的软数据(softdata)。下面基于经典的渐近线理论进行非线性拟合估计置信区间的估计。这里我们定义对应于任何表观热惯量ATI_i，土壤水分回归估计值表示为

其中，是系数向量。模型的具体形式可以参考公式(8)。这里，定义随机变量T_B：

其中：

这里，SM₀表示回归估计置信区间中土壤水分的可能取值；表示对应于指定表观热惯量ATI₀的土壤水分回归估计值；p表示模型参数个数；表示模型回归估计值误差的标准差；表示大小为n*p的雅克比矩阵，在雅克比矩阵的第j列，包括了n个渐近线理论告诉我们，公式(15)中T_B分布的极限(当n趋向于无穷大)是一个标准正态分布N(0，1)。然而，实际情况下，n是一个比较小的定植，不一定为正态分布而且未知。不过，当n不是无穷大时，我们可以利用学生t分布t_n-p来较为精确得逼近真实的分布，这样，SM₀的概率分布可以由学生t分布表示为：

这里，F_S和f_S分别为SM₀的概率分布函数和概率密度函数；为在自由度n-p情况下t分布的概率密度函数。把T₁作为SM₀的函数，在公式(19)两边对SM₀求导，可以得到由表示的f_S：

上面介绍了如何将MODIS ATI预处理成为概率分布形式的土壤水分软数据。在具体升尺度过程中，由于每个微波像元足迹范围内对应太多MODIS ATI像元，具体如图4所示，如果表示为多个土壤水分软数据，那么在以下步骤中的贝叶斯条件化涉及的数值积分过程中运算量就会过大，所以本研究中，将30km*30km区域等分为四块，以每一块的MODIS ATI的像元平均值来求取对应的土壤水分软数据，则相当于每个升尺度范围内有最多四个土壤水分软数据，具体方案如图4中所示。

土壤水分空间协方差的估计

在贝叶斯条件化过程中，如果事件之间非独立，那么需要考虑事件之间的相关性。在图4中，我们定义区域“真值”为事件A1，四个土壤水分软数据依次为事件B1，B2，B3，B4。根据地理学第一定律，事件{A1，B1，B2，B3，B4}在空间上彼此相关，这里我们利用协方差函数(covariogram)来描述事件之间的空间相关性。由于土壤水分具有较大的时空变异性，所以我们需要在每一个升尺度日期估计土壤水分空间协方差函数。一般情况下，区域土壤水分站点观测较稀疏，不能以此构建稳定的协方差函数，所以本研究中我们在大尺度土壤水分观测网络覆盖范围(100km*100km)内利用MODIS ATI的土壤水分回归估计数据来建立每天的试验协方差(experimental covariance)，然后利用块金模型和指数模型嵌套的模型来模拟理论协方差函数。

确定理论协方差函数之后，我们可以得到在相关性范围内任意距离上的协方差，严格意义上这个协方差反映的是1km*1km的区域之间的空间协方差。而这里的事件B1，B2，B3，B4均代表了15km*15km的范围，事件A1代表了30km*30km的范围，这样事件之间的协方差应该是块与块(block to block)之间的协方差。在本研究中，我们以块与块之间所有空间协方差的平均值来代表块与块之间的协方差。图4-5给出了基于算术平均方法计算块与块之间协方差的示意图，块与块之间协方差C_BB′可以表示为：

贝叶斯条件化

贝叶斯条件化其实就是贝叶斯公式在本发明中的应用。根据贝叶斯公式，目标变量的后验概率密度函数(post_pdf)可以由先验概率密度函数pre_pdf和概率形式的辅助数据表示为：

其中，A＝∫_RdF_S(X_B)f_G(X_B)是归一化系数；x_k表示目标变量，即升尺度后的土壤水分；F_S表示X_B的概率分布函数；f_G表示一般数据的概率密度函数；R是辅助数据向量X_B的变化范围；x_k|X_B表示在已知知识(辅助数据)X_B后x_k可能发生的值。更详细的介绍及推理请参考Christakos and Li(1998)and Christakos(2000)的贝叶斯最大熵理论，这里的先验概率分布为步骤一中基于站点数据估计的先验概率分布，概率形式的辅助数据位步骤二中基于MODIS ATI估计的土壤水分软数据。

得到目标变量的后验概率密度函数f_K(x_k)之后，我们可以通过最大化后验概率密度函数求取目标变量的估计最大化后验概率密度即求解以下方程：

这里，即为由稀疏地面观测升尺度后的30km*30km的区域土壤水分，升尺度过程中融合了土壤水分地面站点观测以及来自更高分辨率遥感数据的辅助信息。

结果与分析

(1)SM和ATI的回归关系和置信区间估计

大尺度范围内土壤水分站点观测一般较为稀疏，为了和实际情况保持一致，在本发明实施例中我们选择青藏高原大尺度的土壤水分观测网络表示较稀疏的土壤水分实际观测点，而更为密集分布的中尺度土壤水分网络观测作为升尺度方法的验证数据。在本发明中我们主要在中尺度范围内测试我们构建的升尺度方法。我们需要构建SM和ATI之间的经验关系，而在实际情况下遥感足迹范围内可能只有3-5个土壤水分观测点，仅根据这几个点建立的经验关系可信度不高，一般在遥感像元周围更大的范围内会有更多的土壤水分观测，我们考虑用升尺度范围周边更多的土壤水分站点数据来共同建立更稳定的SM和ATI的经验关系。在本发明实施例中，虽然在30km*30km范围内有中尺度网络的20个节点，但是我们仍然考虑利用大尺度土壤水分的观测数据来建立此经验关系，这样会尽量和实际情况保持一致，且可以减小中尺度网络数据的使用程度，增加最后验证过程的客观性。这样，我们选择了在2010年8月1日-2011年12月31日期间大尺度网络中土壤水分时间序列较为完整的28个站点，这28个站点的分布如图1中的large network_1。图6中显示了这28个站点土壤水分平均值的时间序列，同时，图6中包括了大尺度网络范围内MODIS ATI所有像元平均值的时间序列。在上文中提高，由于云覆盖等的影响，空间上某些区域MODIS ATI是缺失的，图6为可获取的ATI的统计结果。从SM和ATI的时间变化趋势可以初步看出，在大尺度范围内SM和ATI之间存在着一定的相关关系，这里，MODIS ATI可以作为土壤水分升尺度有效的辅助数据。

图7是站点土壤水分和MODIS ATI之间的非线性回归结果。需要说明的是，我们需要SM和MODIS ATI之间的稳定的经验关系，我们更关心的是较大尺度SM和ATI之间的关系，而不是站点土壤水分和像元平均ATI之间的关系，所以，本实施例中利用多个站点土壤水分平均值和每个站点对应的MODIS ATI像元平均值建立的此经验关系。为了减小不确定性，同时我们对此数据集作了如下预处理，首先我们提取了所有土壤水分站点观测以及对应的MODIS AIT像元值，数据集中SM和AIT一一对应，然后我们对ATI数据进行从小到大排序，当然，根据一一对应关系，同时调整了SM数据的排列顺序，最后，按照ATI从小到大的顺序，对SM和ATI分别依次每80个数据取平均值。图7为预处理之后数据的回归结果。图中绿色圆圈为SM平均值和ATI平均值的散点图；红色的曲线为公式(13)的回归曲线。具体回归参数如图7图中所示，此非线性回归的R²为0.95。图中蓝色的虚曲线为回归曲线在95％置信水平下的估计置信区间，同时根据本节标题4中介绍的方法我们可以估计出估计置信区间内土壤水分的概率密度分布。图7中右子图显示了在MODIS ATI为0.08时对应的土壤水分的估计置信区间以及在置信区间内的概率密度分布，这样，我们可以将由所有辅助数据(MODIS ATI)反演为土壤水分，并表达为我们所需要的概率形式的软数据。

Covariogram估计结果

在试验期间的每一天，我们将根据MODIS ATI回归反演得到的土壤水分数据来构建协方差函数。由于在某些天，大部分区域MODIS ATI可能会缺失进而无法继续反演土壤水分，这将可能由于样本量不足或是样本分布不合理导致无法构建合理的协方差函数。为此，我们对MODIS ATI数据在时间序列上进行了插补，当然，为减小升尺度不确定性，此插补程序只应用在了构建协方差函数过程中，升尺度中的软数据只利用未查补的原始数据。这里的时间序列插补方法基于我们的一个统计结论，我们定义：

ATI(i)_mean＝[ATI(i-3)+…+ATI(i-1)+ATI(i+1)+…+ATI(i+3)]/6 (24)

其中，i表示指定的某一天，ATI(i)_mean为ATI六天的均值(指定天的前三天和后三天)。图8显示了在大尺度网络区域研究时期内ATI(i)和ATI(i)_mean之间的统计相关关系。相关系数为0.78，回归直线和1∶1线基本重合。基于此统计关系对MODIS ATI进行时间序列插补之后，可以在试验时期的每一天可以根据插补后ATI反演的土壤水分建立稳定的协方差函数。图9给出了验证期间随机选取三天的协方差函数。结合本节标题5中块与块之间协方差函数的计算方法，可以估计各事件之间的空间协方差。

三个点的随机情况1

在任意一指定区域，比如在微波遥感像元足迹区域内，如果存在有土壤水分稀疏站点观测，并且MODIS ATI可以作为辅助数据，那么我们可以用本研究中构建的升尺度方法将稀疏站点观测升尺度到区域尺度的土壤水分。但是如何验证升尺度的结果是一个比较棘手的问题，因为我们没有区域尺度的土壤水分地面“真值”。本发明实施例中，我们按如下方案解决这个问题。首先，我们相信一个原则，那就是如果指定区域中土壤水分观测站点越多，那么站点观测平均值对区域土壤水分“真值”代表性越好。而在我们的研究区中，中尺度网络区域和微波遥感像元区域的尺度大小差不多，并且中尺度网络区域中包括了高达20个土壤水分观测站点，而在微波像元尺度这样的土壤水分观测的分布已经相当密集，所以根据以上假设，此区域中这20个站点的平均土壤水分已经可以非常接近区域“真值”，这意味着如果只利用少数站点的升尺度结果和20个站点平均土壤水分越接近，那么升尺度精度就越高。这样，虽然没有地面“真值”，我们可以在中尺度网络区域内，通过20个土壤水分观测站点的平均值间接得对我们的升尺度方法进行验证。一般情况下，在像元足迹范围内的站点观测只有3-5个，这里首先尝试只有三个站点的情况。我们选择在2011年6月5日-2011年9月15日时间内在中尺度网络范围内进行升尺度方法的测试和验证。在中尺度网络中共有20个土壤水分观测点，为保证合理，我们需要从20个站点中随机选择三个站点作升尺度使用，20个站点的平均值作相应的验证数据使用。本发明实施例中，随机列出了三个站点随机组合的两种情况。需要特别说明的是作为升尺度软数据使用的MODIS ATI并没有进行时间序列插补，以空间上可获取的数据为准。在验证过程中，我们同时加入了三个站点的平均值作为升尺度方法的参考，因为目前站点平均方法在遥感产品验证中仍然是使用最多的方法。

首先，我们给出了三个站点的第一种随机情况下升尺度方法的表现。随机选择的三个站点为M7，M16and M21，他们在中尺度网络中的位置可以参考图1。在2011年6月5日-2011年9月15日时期内，我们以日为单位进行升尺度，图10显示了升尺度及验证结果。图10-a中，红色的曲线表示中尺度网络的20个站点的平均值STAM(20)，即升尺度的验证数据，我们发现区域土壤水分随着时间在一直减小，研究时期的末期有所升高；绿色的圆圈表示站点平均方法的升尺度结果，即用于升尺度的三个站点的平均值STAM(3)，和STAM(20)相比，在验证时期的初期，STAM(3)要比STAM(20)大，而在后期要小于STAM(20)，STAM(3)表现出大值区偏大小值区偏小的整体特点；蓝色的五角星是由回归方法反演得到的土壤水分像元平均值(ATI-SM)，从整个时间序列来讲，ATI-SM可以正确反映土壤水分变化趋势，但是相邻几天的ATI-SM往往会剧烈的上下抖动，这可能是由于在卫星过境日大气或是气候原因导致的不确定性，这也是光学遥感自身的弊端；图10-a中紫色的五角星为本研究构建的升尺度方法(UPS-SM)的升尺度结果，整体来讲，UPS-SM和STAM(20)在时间序列上匹配的很好，尤其是在开始的二十天以及最后的二十天，UPS-SM比STAM(3)表现出明显的优势。图10-b包括了STAM(3)vs STAM(20)和UPS-SM vs STAM(20)的散点图。两个散点集均紧紧分布在1∶1线附近，STAM(3)的RMSE为0.022m³/m³，而UPS-SM的RMSE可以降低到0.015m³/m³，这可以说明UPS-SM的升尺度结果可以更接近STAM(20)，即UPS-SM可以比STAM(3)对区域“真值”更具代表性。

三个点的随机情况2

在这一部分，我们给出了三个站点组合的另外一种随机情况。随机选取的三个站点为M5，M13，M19，三个站点的位置同样可以参考图1。和图10-a相同，图11-a显示了升尺度结果验证数据的时间序列。验证数据同样为20个站点的平均值STAM(20)，且此图中符号代表的意义和图10-a中相同。在图11-a中，和STAM(20)相比，在550天-580天三个站点的均值可以很好的代表STAM(20)，在590天-610天之间，STAM(3)要明显小于STAM(20)。然而，UPS-SM可以在整个验证时期内基本和STAM(20)保持一致。图11-b显示了验证结果的散点图。STAM(3)的RMSE为0.035m³/m³，这要比上一个组合情况稍差一些，而同样情况下，UPS-SM的RMSE仍然可以达到0.018m³/m³，这说明在此情况下，UPS-SM的升尺度效果同样可以明显优于STAM(3)。另外，在图11-a中，前半段时期，STAM(3)表现较好，但是在后半时期表现却变差，同时考虑在图10-a中STAM(3)出现高值区高估低值区低估的现象，我们认为三个站点的平均值在试验时期的时间稳定性可能会发生变化，所以在本试验区基于时间稳定性分析的升尺度方法并不合适。相比之下，本研究中的UPS-SM方法在三个站点的情况下，通过融合MODIS ATI反映的湿度信息，可以提供稳定的升尺度估计，可以很好代表土壤水分区域地面“真值”。

三个点所有情况的讨论

在上述内容中，我们列举了三个站点随机组合的两种情况，两种情况下UPS-SM均和STAM(20)很好的吻合。然而，二十个站点中三个站点的组合情况有1140种，而且每个站点的空间代表性也不一致，这样就会导致不同情况下升尺度的精度也会有差异，就像本研究中列出两种情况估计精度会不同，并且上千情况中只列举两种情况并不具有说服力。那么在本部分中，我们考虑了三个站点的所有组合情况并且给出了UPS-SM和站点平均两种升尺度方法的验证结果。图12是两种方法在所有情况下RMSE的比较，x轴表示UPS-SM方法在1140种情况下的RMSE，y轴表示站点平均方法在1140种情况下的RMSE。从散点图中，我们可以发现站点平均方法的最大RMSE接近0.10m³/m³，而UPS-SM方法的最大RMSE小于0.05m³/m³。从1140种情况的平均精度(Mrmse)来讲，站点平均方法的Mrmse是0.03033m³/m³，而UPS-SM方法的只有0.02017m³/m³。根据以上的所有分析，我们可以确定在随机情况1和随机情况2中得出的结论是合理的，即UPS-SM可以完成较理想的升尺度估计，并且可以优于站点平均方法。从散点图分布上我们还可以发现大部分散点分布在1∶1线的上方，而分布在1∶1线下方的散点均在UPS-SM的RMSE较小的区域，而且此时站点平均方法的RMSE更小。这说明在某些情况下，站定平均方法的升尺度结果也是可以很好的代表区域“真值”的，而此时由于在UPS-SM方法中融入了MODIS ATI为辅助数据，UPS-SM的升尺度精度反而不及常用的站点平均方法。在图10和图11的分析中我们也发现通过回归方法反演的土壤水分在局部呈现上下波动的现象，我们可以确定在ATI反演的土壤水分中也会存在不确定性，而当随机选择的三个土壤水分站点均具有很好的空间代表性的时候，即简单的站点平均方法可以提供很好的升尺度效果，而融合MODIS ATI的UPS-SM会在升尺度过程中额外引入部分不确定性，这样就会导致在某些情况下站点平均方法要优于UPS-SM。但是，这种情况只会发生在三个随机站点的空间代表性均很好的情况下，而且虽然UPS-SM的精度不及站点平均方法，但是此时UPS-SM的RMSE也都小于0.02m³/m³。在站点平均方法的RMSE大于0.02m³/m³的所有情况下，UPS-SM都表现出了明显的优势。总的来说，本研究中的UPS-SM在三个站点的情况下可以提供非常不错的升尺度估计，并且要整体上明显优于常用的站点平均方法。

升尺度精度对站点个数的敏感性

在以上部分，我们讨论了UPS-SM方法在只有三个土壤水分地面观测站点的情况下的表现，而地面站点数量对升尺度精度的影响也是我们需要考虑的问题。在这部分，我们将就升尺度精度对地面观测站点数量的敏感性进行讨论，作为比较，讨论中同时也包括了站点平均方法。图13显示了UPS-SM和站点平均方法对土壤水分站点数量的敏感性分析结果，即将升尺度精度表示为地面观测站点数量的函数。讨论同样里基于中尺度土壤水分网络进行，图中的纵轴表示在指定的站点个数下所有站点组合情况升尺度精度的平均值(MRMSE)，和图12中的Mrmse意义相同，图中横轴表示升尺度中涉及的站点个数。图中红色圆圈表示UPS-SM方法的分析结果，蓝色方框表示站点平均方法的分析结果。总的来讲，两种方法的MRMSE均随着地面观测站点数量的增加而减小。有趣的是两种方法的MRMSE在地面站点数量为九个的时候发生相交，在九个点之前，UPS-SM的MRMSE要小于站点平均方法的MRMSE，而且站点越少，UPS-SM方法比站点平均方法的优越性越大，这可能主要得益于在UPS-SM中融入了MODIS ATI反映的土壤湿度信息；当站点数量大于九个的时候，站点平均方法的MRMSE仍然会随着站点数量的增加而减小，这是因为土壤水分观测站点的总数量一定，参与升尺度的站点数越多，则和验证数据(20个站点平均值)越接近，可以预见当参与升尺度的站点个数接近20，站点平均方法的MRMSE可以趋向于0；而UPS-SM的MRMSE在参与升尺度的站点数量大于9个以后却逐渐趋于平缓，我们认为这个时候UPS-SM方法的不确定性是由MODIS ATI引入的不确定性主导，我们在之前分析过在升尺度过程中融入MODIS ATI的同时也会引入一定的不确定性，这样，UPS-SM的MRMSE会随着站点数量的增加趋于稳定。所以，通过以上讨论，我们可以知道，在本发明实施例研究区，当可参与升尺度的土壤水分地面站点的数量小于九个的时候，UPS-SM方法可以提供非常精确的升尺度估计，当在地面站点数量为九个的时候，UPS-SM方法和站点平均方法的升尺度精度差不多，而当地面站点数量大于九个的时候，UPS-SM方法的精度上升空间有限，站点平均方法的精度仍然会随着站点数量的增加不断提高，所以，当站点数量增加到一定程度，站点平均方法反而会提供更高精度的升尺度估。

结论

本发明实施例中，在贝叶斯理论框架下，通过融合MODIS ATI，我们构建的土壤水分稀疏站点观测数据的升尺度方法(UPS-SM)可以成功的将此处稀疏站点土壤水分站点观测升尺度为区域尺度土壤水分。在升尺度精度对地面土壤水分观测站点数量的敏感性分析中，同样得到了有趣的结论，即随着升尺度区域中土壤水分观测站点的数量的增加，UPS-SM和站点平均方法的RMSE都在减小，当站点个数小于九个的时候，UPS-SM方法可以为区域尺度土壤水分提供可信的估计，其精度要比常用的站点平均升尺度方法更好，而且观测点的数量越少，UPS-SM比站点平均方法的优势越大，当然此研究中UPS-SM最少需要三个地面站点；当站点个数等于9个的时候，UPS-SM方法和站点平均方法给出的精度基本相同；当站点个数大于9个的时候，UPS-SM方法的精度变化趋于平缓，而站点平均方法的RMSE可以随着站点个数的增加继续减小。

综上所述，本发明中构建的UPS-SM升尺度方法可以完成将土壤水分站点观测成功升尺度到区域尺度土壤水分，尤其在站点土壤水分比较稀疏时UPS-SM可以表现出更大的优势。同时，在未来的工作中，UPS-SM方法可以在更多的区域及下垫面类型尝试。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种基于贝叶斯理论的土壤水分站点数据升尺度方法，其特征在于，包括以下步骤：

步骤一、基于升尺度区域中稀疏站点观测数据估计目标变量的先验概率密度分布pre_pdf；

步骤二、通过建立土壤水分SM和MODIS ATI之间的非线性回归关系，将MODIS ATI反演为SM，并且估计土壤水分非线性回归的估计置信区间以及概率分布作为概率形式的软数据；所述非线性回归的估计置信区间的估计利用公式为：

其中，f_s为SM₀的概率密度函数；为在自由度n-p情况下t分布的概率密度函数；SM₀表示回归估计置信区间中土壤水分的可能取值；表示对应于指定表观热惯量ATI₀的土壤水分回归估计值；p表示模型参数个数；表示模型回归估计值误差的标准差，的公式为：

其中，是系数向量，表示大小为n*p的雅可比矩阵，在雅可比矩阵的第j列，包括了n个j＝1，…，p；

步骤三、通过贝叶斯公式融合目标变量的先验分布以及来自于MODIS ATI的概率形式的辅助信息，获取目标变量的后验概率密度分布函数post_pdf；

步骤四、通过后验概率密度函数的最大化求取概率最大时目标变量的取值，即通过融合MODIS ATI信息将站点土壤水分观测升尺度后的土壤水分。

2.根据权利要求1所述的一种基于贝叶斯理论的土壤水分站点数据升尺度方法，其特征在于：所述步骤一中目标变量是指由区域稀疏站点观测数据估计的升尺度后的区域尺度土壤水分，估计方法具体为正态整体均值的区间估计，根据概率论和数理统计原理，构建统计量：

先验概率密度分布f_G(μ)为：

其中，f_G(μ)为总体的先验概率密度分布；为样本均值，S²为样本方差，T_A服从自由度为n-1时的学生t分布，而T_A的分布只依赖于样本个数n，并不依赖于总体均值μ和方差σ²；是T_A的概率密度分布。

3.根据权利要求1所述的一种基于贝叶斯理论的土壤水分站点数据升尺度方法，其特征在于：所述步骤二中选择的土壤水分SM和表观热惯量ATI之间的非线性回归关系表达式为：

其中，a，b，c为系数。

4.根据权利要求1所述的一种基于贝叶斯理论的土壤水分站点数据升尺度方法，其特征在于：所述步骤三还包括，需要构建块与块之间的空间协方差C_BB′，空间协方差C_BB′可以表示为：

5.根据权利要求1所述的一种基于贝叶斯理论的土壤水分站点数据升尺度方法，其特征在于：所述步骤三中目标变量的后验概率的获取方法是基于贝叶斯原理，其中先验分布为基于站点数据估计的目标变量的先验分布，辅助信息为通过遥感数据获取的概率分布软数据。