CN110503142B - 基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法 - Google Patents

Abstract

本发明提供了一种基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法，针对具有季节变化规律的地学参数的遥感产品存在的时空不连续、时空分辨率单一的问题，在传统的时间序列统计模型中增加时空随机效应模型组分，构造时空序列加法过程模型，嵌入到层次贝叶斯融合模型中，利用多源遥感产品信息的互补性，进行产品信息融合，得到时空完整的、多尺度、高精度、空间信息量丰富的融合遥感产品。本发明可以满足地学参数真实的时空过程模拟、参数的多尺度无缝转换、遥感产品及动态过程模拟参数的不确定表达、地学参数局部时空结构模拟，并且实现快速高效的海量遥感数据计算，得到高精度、局部空间模式丰富的多尺度时空完整的融合数据。

Description

基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯 融合方法

技术领域

本发明属于地学遥感技术领域，具体涉及一种基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法。

背景技术

用定量方法反演地表参数是当前遥感数据最重要的应用。遥感产品的生产者一直致力于研究的问题是通过改进传感器性能、反演策略及反演方法来提高产品的精度，但是遥感观测数据获得的都是瞬时值，而地表过程具有时空连续性，且绝大多数遥感观测和地表过程参数变量之间的关系是隐含的，导致反演的难度较大，即使是微小的误差也可能导致“病态”反演。

遥感产品的应用者根据应用环境的不同，提出了多种统计方法进行单遥感产品时空插值或多源遥感产品融合插值来提高现有遥感产品的时空完整性、精度及空间结构信息。如滤波方法、数据同化方法、多元数据融合方法等。滤波方法缺乏对地学参数动态过程模型的利用，且方法不能同时将历史的先验知识、数据的不确定性、时空异质性同时纳入方法中，噪声和真实信息取舍中存在较大的不确定性。虽然数据同化方法可以将地学观测数据与动态过程模型相结合，但动态过程模型通常选用物理模型，虽然物理模型具有明确的物理机理，能精确地模拟地学参数随时间的演变过程，但基于物理机理的过程模型往往需要大量的地学参数观测作为驱动，而这些参数的可获取性是有限的，为数据同化带来了一定的难度。目前应用的时空融合方法中对地学参数的时空过程模型利用较少，且不能同时满足融合过程中需要的尺度无缝转换、数据及参数的不确定性定量表达。

例如在公开号为CN107341513A，专利名称为基于稳健的固定阶数滤波模型的多源海洋表面温度遥感产品融合方法中虽然可以满足融合过程中需要的尺度无缝转换，产品的空间完整性融合结果也可以实现全覆盖，但是稳健固定阶数滤波模型模拟值作为先验均值嵌入层次贝叶斯，也就是稳健固定阶数滤波模型是单独运行的，没有将该过程模型直接嵌入层次贝叶斯模型，过程模型中的参数为固定值，而不是随机变量，因此参数的不确定性不能科学表达。

发明内容

本发明所要解决的技术问题在于针对上述现有技术的不足，提供一种基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法，利用多源遥感产品的互补性解决现有遥感产品的时空完整性、精度及空间结构问题。

为解决上述技术问题，本发明采用的技术方案是：

基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法，针对具有季节变化规律的地学参数的遥感产品存在的时空不连续、时空分辨率单一的问题，在传统的时间序列统计模型中增加时空随机效应模型组分，构造时空序列加法过程模型，嵌入到层次贝叶斯融合模型中，利用多源遥感产品信息的互补性，进行产品信息融合，得到时空完整的、多尺度、高精度、空间信息量丰富的融合遥感产品；具体包括如下步骤：

S1，构建时空序列加法过程模型，描述宏观趋势特征及局部空间结构特征。在不同时刻存在不同程度的局部空间变异以及在不同的空间位置存在不同程度的局部时间变异，但相邻空间、时间点上又具有明显的相关性，基于独立同分布假设的统计方法很难对其进行科学的分析，因此，所述构建时空序列加法过程模型在传统的时间序列分析方法的基础上引入了时空相关性分析，即在传统的时间序列统计模型基础上，构建时空随机效应模型并与时间序列结构相结合，采用如下公式：

Y(s,t)＝M(s)+S(t,β(s))+R(s,t)+γ(s,t)

Y(s,t)表示潜在的时空过程，M(s)为时间均值空间趋势项，表示每一像元位置上的时间均值空间趋势，用来描述空间特性，M(s)＝(M(1),...,M(s_i))，随空间位置的变化而改变，不随时间变化。S(t；β(s))为季节项，表示大尺度上的具有空间位置特定参数β(s)的季节特征模型，β(s)＝(β(1),...,β(s_i))。R(s,t)为时空随机效应项，用来模拟时空序列去除时间均值及季节均值后的局部时空变异，γ(s,t)是服从正态分布，且均值为0的随机变量，用来描述模型分解的噪声，s为二维欧几里得空间中的位置，s＝{s₁,...s_n}，t表示时间，t＝{t₁,t₂,...}；

S2，构建层次贝叶斯融合模型，将时空随机效应加法过程模型作为时空过程的模拟，嵌入到层次贝叶斯的框架下；

地学参数的时空过程作为一个潜在的过程，无法实际观测，但可以通过分析有限的观测数据，构建数据驱动的经验时空统计模型进行模拟。层次贝叶斯模型可以很方便地集成时空统计模型进行地学参数的后验推断，具体地说是将地学参数的时空过程看作一个时空随机场，根据层次贝叶斯融合理论、地表参数所表现的季节变化特征及遥感观测数据探索性分析的基础，建立经验的时空序列加法过程模型，充分利用各种先验知识对模型中的参数及超参数赋先验分布，利用多源遥感观测数据更新这些先验分布，从而估计地学参数的时空过程均值的后验分布，即

[时空过程|遥感观测数据]∝[遥感观测数据|时空过程，参数][时空过程|参数][参数]

其中，[|]表示条件分布，[时空过程|遥感观测数据]表示地学参数时空过程均值的后验分布，[遥感观测数据|时空过程，参数]表示遥感观测数据模型，该数据模型与传统的统计模型不同，是定义观测数据的误差模型，遥感观测数据被视为有观测误差的数据，[时空过程|参数]表示在给定潜在时空过程模型各种参数基础上的时空过程均值的条件分布，时空过程是无法观测到的，所获得的遥感观测只是在此过程支配下的不同时刻、不同空间位置上的瞬时值，但该过程在层次贝叶斯框架下可以通过模型进行模拟，即构建的时空序列加法过程模型，[参数]表示参数模型，包括数据模型、时空序列加法过程模型中所有的随机变量，这些变量在参数模型阶段需要根据先验知识确定他们的先验分布。

作为本方案的进一步描述，所述构建时空序列加法过程模型中每一时间点的每一像元位置上的时空序列加法过程条件依赖时间均值空间趋势项、季节项和时空随机效应项三个子过程及相应的参数θ₂，

其中，

作为本方案的进一步描述，所述时间均值空间趋势项、季节项和时空随机效应项三个子过程条件依赖于参数集合θ₃＝{θ_M,θ_β,θ_R}，即：

[M,S,R|θ₃]＝[M|θ_M][S|θ_β][R|θ_R]

通过计算时间均值的方法探索分析每个像元的时间均值的空间趋势并可视化表达，有季节规律演变的地学参数通常表现出大尺度的纬向地带性，用简单的线性趋势面模拟，协变量选择截距、纬度、经度，所述时间均值空间趋势项针对每一个像元的表达式为M(s)＝X(·)α，其中，α是系数，X(·)＝(X₁(·),...,X_p(·))是p维的协变量，为了体现时间均值随空间位置的变化，将空间位置(经纬度)作为协变量，即：M(s)＝μ+α₁*s_{i_lat}+β₁*s_{i_lon}，其中，μ为像元的时间均值，α₁与β₁为系数，s_{i_lat}与s_{i_lon}分别表示像元的纬度和经度，北半球为正值，南半球为负值，其作用是对像元时间均值的空间变化微量调整；

将遥感观测移除时间均值空间趋势，剩余残差部分，按纬向地带性规律，分纬度探索季节规律，构建余弦波普曲线季节项，余弦曲线可以复合从一个时间周期到另一个时间周期的平滑变化，同时还能保持季节特性，以纬度、经度为自变量，构建线性结构函数，以使其随着空间位置的变化而变化，所述季节项采用余弦函数表达，

其中，

δ₁与δ₂分别为余弦、正弦系数，按纬度分区为基本单元，不同单元系数不同，

f＝1/t,t＝1,2,...,46，为了区分南北半球不同的季节特性，显示季节特性的空间结构，构造余弦、正弦系数为纬度和经度两个协变量的线性结构如下：

δ₁(s)＝a₀+a₁*s_{i_lat}+a₂*s_{i_lon}

δ₂(s)＝b₀+b₁*s_{i_lat}+b₂*s_{i_lon}

其中，a为系数，a₀，a₁，a₂，b₀，b₁，b₂均为独立同分布的高斯随机变量；

将移除时间均值空间趋势的剩余残差，继续移除季节特征的空间趋势，探索此次剩余的残差部分的自相关关系即时间依赖关系，并根据自相关关系确定时间依赖步长，通常取相关系数达到0.8以上的时间步长建立模型。除了要探索残差的时间依赖关系外，还需要探索残差的局部空间结构，不能简单地把残差看作是均值为0的随机变量，然后根据残差特征构建局部空间结构模型-时空随机效应项，所述时空随机效应项R(s,t)应用马尔科夫随机场理论，在一阶马尔科夫时间演变模型的基础上，增加空间8-邻域加权模型，模拟空间局部变异，即一阶矢量自回归基础上的空间8-邻域加权时空随机效应项：

其中，R(s,t+1)为t+1时刻的随机效应项，H_t为传播矩阵，

为t+1时刻的随机误差项，s_i是索引为i的像元，s_j为索引为i的像元的8-邻域像元，i，j均为整数，

是均值为0，方差为

的随机误差项，时间随机效应项子过程的参数集合为

W_ij为临界权重若两像元相邻，则W_ij＝1，否则W_ij＝0，

表示空间变异组分，R(s,t)随周围8-邻域像元的加权均值的改变而变化。

作为本方案的更进一步描述，所述构建层次贝叶斯融合模型，首先需要定义数据模型，还需要构建多遥感数据的无缝尺度转换模型，数据模型结构表示为[Z(s,t)|Y(s,t),θ₁]，其中，θ₁为数据模型阶段的参数，Z(s,t)表示空间不连续的遥感时空序列观测数据。

作为本方案的更进一步描述，所述时空序列加法过程模型和数据模型中的所有随机变量的先验分布通过参数模型定义，即[θ₁,θ₂,θ₃|θ₄]，其中θ₄为超参数变量，所述时空序列加法过程模型和数据模型的参数集合条件依赖超参数变量θ₄，

[θ₁,θ₂,θ₃|θ₄]＝[θ₁|θ₄(1)][θ₂|θ₄(2)][θ₃|θ₄(3)]

针对所述时空序列加法过程模型的每一个子过程M(s)，S(t；β(s)，R(s,t)，相应的超参数变量为θ₄(3)＝{θ₄(M),θ₄(β),θ₄(R)}，条件分布概念模型为θ₃|θ₄(3)＝[θ_M|θ₄(M)][θ_β|θ₄(β)][θ_R|θ₄(R)]。

作为本方案的进一步描述，所述时空序列加法过程模型、数据模型和参数模型中，经度、纬度为确定性参数，除经度、纬度之外的其他随机变量参数均为服从正态分布的随机变量，

作为本方案的更进一步描述，所述每一个模型阶段的方差均采用共轭先验分布，即逆伽马分布，

其中，IG表示逆伽马分布，A用来标识不同的模型阶段，q_A为形状参数，r_A为尺度参数。

本发明与现有技术相比具有以下有益效果：

1、本发明将传统的时间序列统计模型进行扩展，增加时空随机效应模型组分，构建时空序列加法过程模型，真实准确地模拟具有季节变化规律的地学参数的时空演变过程。

2、本发明将时空序列加法过程模型嵌入层次贝叶斯模型中，地学参数的时空过程模型条件依赖多子过程、多随机变量复合的模式，可以实现两种不同空间尺度遥感产品的无缝集成。

3、本发明中的层次贝叶斯模型可以将复杂的空间结构转换为一系列条件分布的乘积，通过马氏链蒙特卡罗(MCMC)方法，实现海量遥感数据的高效计算。

4、本发明可以实现遥感产品及动态过程模拟参数的不确定表达，并得到融合有多遥感产品信息的高精度、局部空间模式丰富的多尺度时空完整的融合数据。

附图说明

图1是本发明的AMSR-E海洋表面温度遥感产品每个像元的时间均值空间趋势示意图。

图2是本发明的AMSR-E海洋表面温度遥感产品纬度2°子区去除时间均值空间趋势后的季节特征示意图，其中，黑色加粗线条为北半球的季节特征，浅灰色线条为南半球季节特征。

图3是本发明的AMSR-E海洋表面温度遥感产品移除时间均值与季节均值时空趋势后的残差抽样空间自相关函数示意图。

图4是本发明的MODIS、AMSR-E和融合数据的空间完整性评价示意图。

图5是本发明的MODIS、AMSR-E和融合数据的局部方差比较示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明。

本发明包括基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法，针对具有季节变化规律的地学参数的遥感产品存在的时空不连续、时空分辨率单一的问题，在传统的时间序列统计模型中增加时空随机效应模型组分，构造时空序列加法过程模型，嵌入到层次贝叶斯融合模型中，利用多源遥感产品信息的互补性，进行产品信息融合，得到时空完整的、多尺度、高精度、空间信息量丰富的融合遥感产品；具体包括如下步骤：

S1，构建时空序列加法过程模型，描述宏观趋势特征及局部空间结构特征。在不同时刻存在不同程度的局部空间变异以及在不同的空间位置存在不同程度的局部时间变异，但相邻空间、时间点上又具有明显的相关性，基于独立同分布假设的统计方法很难对其进行科学的分析，因此，构建时空序列加法过程模型在传统的时间序列分析方法的基础上引入了时空相关性分析，即在传统的时间序列统计模型基础上，构建时空随机效应模型并与时间序列结构相结合，采用如下公式：

Y(s,t)＝M(s)+S(t,β(s))+R(s,t)+γ(s,t)

Y(s,t)表示潜在的时空过程，M(s)为时间均值空间趋势项，表示每一像元位置上的时间均值空间趋势，用来描述空间特性，M(s)＝(M(1),...,M(s_i))，随空间位置的变化而改变，不随时间变化。S(t；β(s))为季节项，表示大尺度上的具有空间位置特定参数β(s)的季节特征模型，β(s)＝(β(1),...,β(s_i))。R(s,t)为时空随机效应项，用来模拟时空序列去除时间均值及季节均值后的局部时空变异，γ(s,t)是服从正态分布，且均值为0的随机变量，用来描述模型分解的噪声，s为二维欧几里得空间中的位置，s＝{s₁,...s_n}，t表示时间，t＝{t₁,t₂,...}；构建时空序列加法过程模型中每一时间点的每一像元位置上的时空序列加法过程条件依赖时间均值空间趋势项、季节项和时空随机效应项三个子过程及相应的参数θ₂，

其中，

时间均值空间趋势项、季节项和时空随机效应项三个子过程条件依赖于参数集合θ₃＝{θ_M,θ_β,θ_R}，即：

[M,S,R|θ₃]＝[M|θ_M][S|θ_β][R|θ_R]

(1)通过计算时间均值的方法探索分析每个像元的时间均值的空间趋势并可视化表达，有季节规律演变的地学参数通常表现出大尺度的纬向地带性，用简单的线性趋势面模拟，协变量选择截距、纬度、经度，所述时间均值空间趋势项针对每一个像元的表达式为M(s)＝X(·)α，其中，α是系数，X(·)＝(X₁(·),...,X_p(·))是p维的协变量，为了体现时间均值随空间位置的变化，将空间位置(经纬度)作为协变量，即：M(s)＝μ+α₁*s_{i_lat}+β₁*s_{i_lon}，其中，μ为像元的时间均值，α₁与β₁为系数，s_{i_lat}与s_{i_lon}分别表示像元的纬度和经度，北半球为正值，南半球为负值，其作用是对像元时间均值的空间变化微量调整；

(2)将遥感观测移除时间均值的空间趋势，剩余残差部分，按纬向地带性规律，分纬度探索季节规律，构建余弦波普曲线季节项，余弦曲线可以复合从一个时间周期到另一个时间周期的平滑变化，同事还能保持季节特性，以纬度、经度为自变量，构建线性结构函数，以使其随着空间位置的变化而变化，季节项采用余弦函数表达，

其中，

δ₁与δ₂分别为余弦、正弦系数，按纬度分区为基本单元，不同单元系数不同，

f＝1/t,t＝1,2,...,46，为了区分南北半球不同的季节特性，显示季节特性的空间结构，构造余弦、正弦系数为纬度和经度两个协变量的线性结构如下：

δ₁(s)＝a₀+a₁*s_{i_lat}+a₂*s_{i_lon}

δ₂(s)＝b₀+b₁*s_{i_lat}+b₂*s_{i_lon}

其中，a为系数，a₀，a₁，a₂，b₀，b₁，b₂均为独立同分布的高斯随机变量；

(3)将移除时间均值空间趋势的剩余残差，继续移除季节特征的空间趋势，探索此次剩余的残差部分的自相关关系即时间依赖关系，并根据自相关关系确定时间依赖步长，通常取相关系数达到0.8以上的时间步长建立模型。除了要探索残差的时间依赖关系外，还需要探索残差的空间局部结构，不能简单地把残差看作是均值为0的随机变量，然后根据残差特征构建局部空间结构模型-时空随机效应项，时空随机效应项R(s,t)应用马尔科夫随机场理论，在一阶马尔科夫时间演变模型的基础上，增加空间8-邻域加权模型，模拟空间局部变异，即一阶矢量自回归基础上的空间8-邻域加权时空随机效应项：

其中，R(s,t+1)为t+1时刻的随机效应项，H_t为传播矩阵，

为t+1时刻的随机误差项，s_i是索引为i的像元，s_j为索引为i的像元的8-邻域像元，i，j均为整数，

是均值为0，方差为

的随机误差项，时间随机效应项子过程的参数集合为

W_ij为临界权重若两像元相邻，则W_ij＝1，否则W_ij＝0，

表示空间变异组分，R(s,t)随周围8-邻域像元的加权均值的改变而变化。

S2，构建层次贝叶斯融合模型，首先需要定义数据模型，还需要构建多遥感数据的无缝尺度转换模型，数据模型结构表示为[Z(s,t)|Y(s,t),θ₁]，其中，θ₁为数据模型阶段的参数，Z(s,t)表示空间不连续的遥感时空序列观测数据。

时空序列加法过程模型和数据模型中的所有随机变量的先验分布通过参数模型定义，即[θ₁,θ₂,θ₃|θ₄]，其中θ₄为超参数变量，时空序列加法过程模型和数据模型的参数集合条件依赖超参数变量θ₄，

[θ₁,θ₂,θ₃|θ₄]＝[θ₁|θ₄(1)][θ₂|θ₄(2)][θ₃|θ₄(3)]

针对时空序列加法过程模型的每一个子过程M(s)，S(t；β(s)，R(s,t)，相应的超参数变量为θ₄(3)＝{θ₄(M),θ₄(β),θ₄(R)}，条件分布概念模型为θ₃|θ₄(3)＝[θ_M|θ₄(M)][θ_β|θ₄(β)][θ_R|θ₄(R)]。

时空序列加法过程模型、数据模型和参数模型中，经度、纬度为确定性参数，除经度、纬度之外的其他随机变量参数均为服从正态分布的随机变量，

每一个模型阶段的方差均采用共轭先验分布，即逆伽马分布，

其中，IG表示逆伽马分布，A用来标识不同的模型阶段，q_A为形状参数，r_A为尺度参数。

实施例1

本算法已经成功应用于4km的MODIS和25km的AMSR-E海洋表面温度遥感产品的融合。

S1，构建时空序列加法过程模型的探索性数据分析，以AMSR-E海洋表面温度遥感产品数据为例，结果如下：

图1为每个像元的时间均值空间趋势示意图，图中明显的表现出从赤道向南北两个半球存在明显的大尺度空间趋势，即从低纬度到高纬度沿纬度SST呈条带状降低趋势，基本表现出一个简单的线性趋势面，用截距、纬度、经度为协变量的线性结构和的结构进行模拟。

图2为像元移除时间均值空间趋势后，每2°纬度的季节特征的空间分布图，从图中可以看出，将每一像元移除空间趋势均值后，从北半球到南半球以2°纬度为统计单元，海洋表面温度表现出了明显的正弦余弦曲线特征的季节特性。因此，去除了时间均值的空间趋势残差，可以用余弦波普曲线的季节模型模拟，表达季节特征。

图3为一出时间均值与季节均值时空趋势后的残差抽样自相关函数图，可以明显看出在移除季节特征和时间均值等空间趋势后的残差项之间的自相关关系，即很强的时间依赖关系。步长为3的自相关函数的系数达到0.6以上，而步长为1的自相关系数达到0.8以上。

此外，移除空间趋势均值与季节均值后的残差除了像图3所示的表现强烈的时间依赖关系外，每周的残差分布还表现出了明显的空间结构。因此不能简单地把残差项视为0均值的随机变量，而应构建局部空间模型模拟其局部空间结构，本发明中引入时空随机效应模型描述这种局部空间变异。

综上，海洋表面温度时空过程模型为时间均值趋势项、余弦波普季节特征项、时空随机效应项三个时空结构的和复合而成。在传统的时间序列统计模型的基础上，进行扩展，增加时空随机效应项，构建时空序列加法过程模型，更加符合SST时空演变的复杂特征。

S2，层次贝叶斯时空融合，经MODIS和AMSR-E海洋表面温度的遥感产品验证结果如下：数据模型，

空间多层次尺度转换模型，

其中，

为观测误差方差，在融合阶段，只需要给定每个像元潜在海洋表面温度时空序列加法过程Y条件下的观测数据Z，即遥感观测的海洋表面温度的误差分布，为了计算效率，在每一个子区，又以经度2°为计算单元进行并行计算。

时空序列加法过程模型和参数模型采用公式均与上述公式相同。

融合结果评价：以漂流浮标实测数据为参考，分别评价了融合数据、MODIS数据和AMSR-E的空间完整性、整体精度、局部精度和局部空间结构保持能力。

空间完整性评价：经计算，年平均有效观测率AMSR-E为87.53％，MODIS为80.38％，融合数据为100％，如图4所示，由此可见融合后的SST实现了海洋像元的全覆盖。

精度评价：以漂流浮标SST为参考值，计算三种数据的平均偏差

误差标准差(STD)、均方根误差(RMSE^a)、相关系数R^a及回归关系。表1融合数据、MODIS数据及AMSR-E数据的整体精度验证与比较

当MODIS为有效观测时，不论AMSR-E是否为有效观测，融合数据的平均偏差、误差标准差、均方根误差均比MODIS小，说明融合数据的精度均高于MODIS。

表2MODIS为有效观测时融合SST、MODIS SST局部精度验证与比较

当MODIS为无效观测时，融合数据局部精度评价结果，绝对平均偏差均在0.5℃范围内，符合全球海洋数据同化试验融合绝对偏差在0.5℃范围内的要求。

表3MODIS为无效观测时融合数据局部精度验证

局部空间结构保持能力评价：采用局部方差定量地评价融合数据局部空间结构保持能力的大小，如图5所示，融合后的数据局部细节信息量明显丰富于AMSR-E，接近于MODIS，局部方差的年平均值，MODIS为0.2409℃，AMSR-E为0.0562℃，融合数据为0.2320℃，说明融合后的数据保持了细分辨率的MODIS的局部空间细节。

以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制。凡是根据发明技术实质对以上实施例所作的任何简单修改、变更以及等效变化，均仍属于本发明技术方案的保护范围内。

Claims (1)

1.基于时空随机效应加法过程模型的多源遥感产品层次贝叶斯融合方法，其特征在于，包括如下步骤：

S1，构建时空序列加法过程模型，即在传统的时间序列统计模型基础上，构建时空随机效应模型并与时间序列结构相结合，采用如下公式：

Y(s,t)＝M(s)+S(t；β(s))+R(s,t)+γ(s,t)

Y(s,t)表示潜在的时空过程，M(s)为时间均值空间趋势项，表示每一像元位置上的时间均值空间趋势，M(s)＝(M(1),...,M(s_i))，随空间位置的变化而改变，不随时间变化；S(t；β(s))为季节项，表示大尺度上的具有空间位置特定参数β(s)的季节特征模型，β(s)＝(β(1),...,β(s_i))；R(s,t)为时空随机效应项；γ(s,t)是服从正态分布，且均值为0的随机变量，s为二维欧几里得空间中的位置，s＝{s₁,...s_n}，t表示时间，t＝{t₁,t₂,...}；

所述构建时空序列加法过程模型包括时间均值空间趋势项、季节项和时空随机效应项三个子过程及相应的参数θ₂，

其中，

所述时间均值空间趋势项、季节项和时空随机效应项三个子过程条件依赖于参数集合θ₃＝{θ_M,θ_β,θ_R}，即：

[M,S,R|θ₃]＝[M|θ_M][S|θ_β][R|θ_R]

所述时间均值空间趋势项M(s)＝X(·)α，其中，α是系数，X(·)＝(X₁(·),...,X_p(·))是p维的协变量，为了体现时间均值随空间位置的变化，将经纬度作为协变量，即：M(s)＝μ+α₁*s_{i_lat}+β₁*s_{i_lon}，其中，μ为像元的时间均值，α₁与β₁为系数，s_{i_lat}与s_{i_lon}分别表示像元的纬度和经度，北半球为正值，南半球为负值；

所述季节项采用余弦函数表达，

其中，

δ₁与δ₂分别为余弦、正弦系数，

f＝1/t,t＝1,2,...,46，为了区分南北半球不同的季节特性，构造余弦、正弦系数为纬度和经度两个协变量的线性结构如下：

δ₁(s)＝a₀+a₁*s_{i_lat}+a₂*s_{i_lon}

δ₂(s)＝b₀+b₁*s_{i_lat}+b₂*s_{i_lon}

其中，a为系数，a₀，a₁，a₂，b₀，b₁，b₂均为独立同分布的高斯随机变量；

所述时空随机效应项R(s,t)应用马尔科夫随机场理论，在一阶马尔科夫时间演变模型的基础上，增加空间8-邻域加权模型，模拟空间局部变异，即一阶矢量自回归基础上的空间8-邻域加权时空随机效应项：

其中，R(s,t+1)为t+1时刻的随机效应项，H_t为传播矩阵，

为t+1时刻的随机误差项，s_i是索引为i的像元，s_j为索引为i的像元的8-邻域像元，i，j均为整数，

是均值为0，方差为

的随机误差项，时间随机效应项子过程的参数集合为

W_ij为临界权重若两像元相邻，则W_ij＝1，否则W_ij＝0，

表示空间变异组分，R(s,t)随周围8-邻域像元的加权均值的改变而变化；

S2，构建层次贝叶斯融合模型，将时空随机效应加法过程模型作为时空过程的模拟，嵌入到层次贝叶斯的框架下；

所述构建层次贝叶斯融合模型，首先需要定义数据模型，数据模型结构表示为[Z(s,t)|Y(s,t),θ₁]，其中，θ₁为数据模型阶段的参数，Z(s,t)表示空间不连续的遥感时空序列观测数据；

所述时空序列加法过程模型和数据模型中的所有随机变量的先验分布通过参数模型定义，即[θ₁,θ₂,θ₃|θ₄]，其中θ₄为超参数变量，所述时空序列加法过程模型和数据模型的参数集合条件依赖超参数变量θ₄，

[θ₁,θ₂,θ₃|θ₄]＝[θ₁|θ₄(1)][θ₂|θ₄(2)][θ₃|θ₄(3)]

针对所述时空序列加法过程模型的每一个子过程M(s)，S(t；β(s)，R(s,t)，相应的超参数变量为θ₄(3)＝{θ₄(M),θ₄(β),θ₄(R)}，条件分布概念模型为θ₃|θ₄(3)＝[θ_M|θ₄(M)][θ_β|θ₄(β)][θ_R|θ₄(R)]；

所述时空序列加法过程模型、数据模型和参数模型中，经度、纬度为确定性参数，除经度、纬度之外的其他随机变量参数均为服从正态分布的随机变量，

每一个模型阶段的方差均采用共轭先验分布，即逆伽马分布，

其中，IG表示逆伽马分布，A用来标识不同的模型阶段，q_A为形状参数，r_A为尺度参数。

Images