CN107341513A - 基于稳健的固定阶数滤波模型的多源海洋表面温度遥感产品融合方法 - Google Patents

Abstract

本发明公开了一种基于稳健的固定阶数滤波模型的多源海洋表面温度(SST)遥感产品融合方法，综合SST遥感产品时空融合过程中的尺度转换问题、产品的不确定性表达问题及SST时空过程复杂的空间结构和时间结构，充分利用遥感产品在空间分辨率、时空完整性、精度特征等方面的互补特性，采用基于稳健的固定阶数滤波过程模型的层次贝叶斯框架融合红外SST和微波SST遥感产品，得到高精度的、空间完整的、局部空间模式丰富的细尺度融合数据，本发明对SST空间趋势的模拟更加合理，实现了多源遥感产品间的无缝时空尺度转换、遥感产品不确定性的定量化，适用于遥感大数据量的高效计算。

Description

基于稳健的固定阶数滤波模型的多源海洋表面温度遥感产品 融合方法

技术领域

本发明涉及海洋遥感技术领域，具体涉及一种基于稳健的固定阶数滤波模型的多源海洋表面温度遥感产品融合方法。

背景技术

海表面温度是全球海气耦合系统中的重要参数之一，在海洋和大气之间的能量、水汽等交换过程中起着最基本的作用，已经被广泛应用于上层海洋过程、气候变化检测、海气热量交换、海洋大气数值模拟与预报及海洋大气同化模型等的研究中。

目前获得SST的手段主要包括传统的现场观测和现代的遥感观测。传统的现场观测主要通过船舶、海上固定的和漂流的浮标及沿岸站点近海观测平台等常规观测系统观测，精度高，但空间覆盖范围有一定的限制，且时空分布不连续。因此基于空间上离散分布的浮标等传统的海面温度观测，很难满足多尺度海洋过程模拟等研究的需要。虽然卫星SST以其较高的时空连续性、近实时的观测特点被广泛应用，高精度、高时空分辨率、时空完整的卫星SST数据集越来越成为局部尺度下海气耦合系统研究中不可缺少的数据源。但目前单一传感器反演SST产品受到传感器性能、卫星轨道的变化、大气环境、海洋表面环境、反演算法等因素的影响，单一传感器的SST产品往往存在不同程度的不确定性，还不能满足遥感产品的业务化生产及科学研究的需要，尤其是局部尺度下的科研与生产的需要，SST遥感产品仍需改进。主要表现在以下几个方面：目前单一传感器的SST遥感产品存在时空不完整、分辨率难以满足多尺度模型建模的需要；单一传感器的SST产品普遍存在缺值像元，时空不完整；由于反演策略的差异，不同传感器反演得到的SST产品之间缺少物理上的一致性。

近年来，针对单一传感器定量遥感产品的改进，一方面，产品发布者通过改进反演算法提高产品精度，陆续发布改进版本的产品。另一方面，研究者从不同的角度进行了广泛的研究，提出了一系列定量遥感产品改进算法，主要有滤波方法、数据同化方法、多源遥感产品时空融合方法等。

滤波方法主要是通过时间滤波和空间插值得到时空完整的、符合时空演变规律的参数时空数据集，广泛应用于NDVI、LAI、气溶胶等遥感产品的时间序列重建。基于时间滤波的遥感产品时间序列重建算法中，主要存在以下问题：缺乏对遥感参数动态过程模型的利用；历史的先验知识利用不足；在时间序列重建时没有对数据的不确定性进行定义和量化；没有考虑空间异质性；大多数方法在计算的过程中，时间窗口的确定受主观经验的影响，在噪声和真实信息取舍中存在较大的不确定性，因此滤波后的产品精度受到一定的影响。虽然空间滤波方法充分考虑了空间异质性，仍存在以下问题：空间滤波根据像元邻域空间关系进行插值，没有考虑数据间的时间动态过程的依赖关系，且由于最优无偏估计的条件使得空间插值的方差比原始数据的方差低，具有明显的平滑效应，局部空间细节信息不易保持；受观测数据密度的影响，在大片缺失数据区域空间插值误差较大；受协方差求逆的限制，很难实现大数据量的高效计算，虽然FRK和FRF方法通过多分辨率小波基函数实现空间降维，解决了传统地统计方法中大数据量协方差求逆的计算问题，计算效率较高，但仍未解决时空插值后的数据过于平滑的问题，且这两种方法对SST这种不规则区域的时空插值，多分变率小波基函数计算不稳定，需要改进。

数据同化方法在考虑数据时空分布的基础上，将遥感时间数据与动态过程模型相结合，在数值模型的动态运行过程中融合新的观测数据，广泛应用于大气科学、海洋科学和陆面过程参数的估计中。虽然基于动态过程模型的数据同化技术将遥感观测与动态过程模型(如作物生长模型、生态过程模型、水文动态模型、大气数值模型、海洋动态模型、陆面过程模型等)相结合，具有明确的物理机理，有效利用不同时空分辨率的观测数据，将观测信息通过时间演变和物理属性规律的一致性约束而融入模型中，可以精确地模拟过程随时间的演变规律，因此得到的同化产品具有较高的精度。但在同化的过程中，基于物理机理的过程模型往往需要大量的参数观测值作为驱动，而通过遥感获得的参数是有限的，这为数据同化带来了一定的难度。且目前SST产品的常用同化方法，均未涉及时空过程模型的构建，也没有提供无缝的尺度转换模型。

多源遥感产品的时空融合方法。目前的融合方法，主要包括多分辨率树(MRT)方法、经验正交函数法(EOF)、数据插值经验正交函数法(DIEOF)、小波分析法、最小二乘法、人工神经网络技术、贝叶斯最大熵方法(BME)等，但均对参数的时空过程模型利用不够。MRT模型可以填充缺失值，降低误差，但对于尺度转换问题需要设置严格的假设条件；EOF方法的优点之一就是它只根据可获得的数据计算必要的信息，没有主观的参数需要估计，该方法计算量小。但该方法在时间序列重建中，如果图像的有效像元太少，或同一个像元位置的时间序列存在较多的缺失值，都将影响重建序列的精度，甚至影响全部插值的结果。SST可能存在大范围空间上或时间上或时空上连续的缺值像元，因此EOF插值方法难以提高SST插值精度。小波分析法虽然能够保持高分辨率数据的细节特征，但仅仅依靠该方法，不能获得全覆盖的SST时空数据，还需要用其他方法进行融合后缺值数据的插补。最小二乘法虽然不需要人为定义参数的假设条件，降低了产品的不确定性，但误差的时空变化需要进一步细化。BME方法在贝叶斯框架下显示地表达了各种信息、知识的不确定性，利用时空协方差结构很好地表达了数据的时空变异，充分考虑了多源知识的不确定性，不仅可以融合多传感器产品、点位实测数据，还可以有效考虑时间变化规律的信息，融合后的数据很好地保持高分辨率遥感产品的局部细节，具有较丰富的信息量。但该方法借助时空协方差模型实现时空时差，存在较大的计算速率问题，且难以与时空过程模型相结合。

发明内容

为解决上述问题，本发明提供了一种基于稳健的固定阶数滤波模型的多源海洋表面温度遥感系统融合方法，利用多源遥感产品的互补性解决单源海洋表面温度遥感产品的空间完整性、精度及空间模式问题。

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

基于稳健的固定阶数滤波模型的多源海洋表面温度(SST)遥感产品融合方法，综合SST时空融合过程中的尺度转换问题、产品的不确定性表达问题及SST时空过程复杂的空间结构和时间结构，充分利用遥感产品在空间分辨率、时空完整性、精度特征等方面的互补特性，采用基于稳健的固定阶数滤波过程模型的层次贝叶斯框架融合红外SST遥感产品和微波SST遥感产品，得到高精度的、空间完整的、局部空间模式丰富的细尺度融合数据，具体包括如下步骤：

潜在的真值SST时空过程在该时空域上的有限的卫星SST观测过程Z(s；t)与其潜在的真值过程Y(s；t)之间的关系表示为：

Z(s；t)＝Y(s；t)+ε(s；t) (1)

其中，{ε(s；t)：s∈D，t∈{1，2，...}}为时空域上的高斯白噪声过程，均值为0，方差是待估参数。在实际中，获得的是在时刻t，空间{s_1，t，...，s_n，t}上的n_t维的有效观测于此相对应，则 所有时刻的观测为条件依赖过程的高斯分布为：

S1、建立稳健的固定阶数滤波模型

潜在的真值时空过程Y(s；t)为由大尺度宏观趋势μ_t(s)、小尺度空间变异过程v(s；t)和0均值的高斯白噪声ξ(s；t)三部分组成的线性结构；

Y(s；t)＝μ(s；t)+v(s；t)+ξ(s；t) (3)

其中，μ(s；t)是一个确定性的时空均值函数，描述SST时空过程的大尺度宏观趋势，μ(s；t)＝X_t(·)′β_t，X_t(·)≡(X_1，t(·)，...，X_p，t(·))′为协变量，β_t≡(β_1，t，...，β_p，t)为未知系数；SST小尺度局部空间变异过程v(s；t)是一个时空随机过程，用时空随机效应模型模拟；ξ(s；t)为模型分解误差，描述随机变异，同地统计学中的块金效应，用0均值方差为的时空高斯白噪声过程模拟；

小尺度局部空间变异v(s；t)，在任一固定的时刻t，为0均值的空间随机效应模型：

v(s；t)＝S_t(s)′η_t (4)

其中，S_t(·)≡(S_1，t(·)，...，S_r，t(·))′是一组r维的细尺度下小波基函数，η_t≡(η_1，t，...，η_r，t)′是一个0均值的高斯随机过程，协方差为K_t，矩阵大小为r×r，即用r维的空间基函数S_t(·)描述任一时刻t的局部空间变异；公式(4)中，如果t＞1，且通过统计模型模拟随机变量η_t的时间依赖关系，通过一阶矢量自回归模型模拟(公式5)，在时空随机效应模型中，细尺度下的小波基函数的个数可以随着时间变化而变化，也可以是随时间不变的，如果随时间而变，则r≡max{r_t}，v(s；t)＝S(s)′η_t；

η_t+1＝H_t+1η_t+ζ_t+1；t＝1，2，.... (5)

其中，H_t+1是r×r的传播矩阵(一阶矢量自回归矩阵)，ζ_t+1是独立于η_t的r维高斯奇异矢量，均值为0，方差var(ζ_t+1)≡U_t+1；

交叉协方差K_t1，t2定义为：

K_t1，t2≡cov(η_t1，η_t2)＝K_t1(H_t2H_t2-1…H_t1+1)′ (6)

相应的，K_t+1＝H_t+1K_tH′_t+1+U_t+1′

根据公式6，则时间步长为1的交叉协方差L_t+1为(公式7)：

L_t+1≡K_t，t+1＝K_tH′_t+1 (7)

综上所述，根据公式(1)，用层次模型表示潜在的真值时空过程Y(s；t)为：

η₁～N(0，K₁) (8)

η_t+1|η₁，...，η_t～N(H_t+1η_t，U_t+1)，t＝1，2，... (9)

Y(s；t)的边际分布为：

Y(s；t)协方差结构为：

其中，I(·)为指示矩阵；

根据公式(12)，可以定义出潜在的真值时空过程Y(s；t)和观测数据Z(r；u)之间的协方差结构：

cov(Y(s；t)，Z(r；u))＝cov(Y(s；t)，Y(r；u)) (13)

S2、构建层次贝叶斯融合模型

将稳健的固定阶数滤波模型作为一个SST时空过程的模拟，嵌入层次贝叶斯的框架下，表示为如下结构：

μ_a(s_a；t)＝μ(s_a；t)+v(s_a；t)+ξ(s_a；t) (17)

μ(s；t)＝X_t(·)′β_t (18)

v(s；t)＝S_t(s)′η_t (19)

其中，μ_a(s_a；t)是确定的，I是单位阵；

在层次贝叶斯框架下，SST时空演变过程μ_a(s_a；t)的先验均值的确定采用稳健的固定阶数滤波模型估计的SST时空模拟值，条件依赖子过程μ(s_a；t)、v(s_a；t)，随机变量参数为两个卫星观测数据的方差，采用逆伽马分布，即：

其中，IG表示逆伽马分布，A用来标识不同的卫星观测；形状参数q_A与尺度参数r_A均为0.1。

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

在层次贝叶斯的融合框架下，可以实现：

1)基于稳健的固定阶数滤波的SST时空过程模型的嵌入。稳健的固定阶数滤波过程模型利用空间降维实现海量数据的高校计算，并且该模型在SST空间趋势模拟过程中既考虑了像元间的空间关系，又结合了像元间的时间依赖关系，因此对SST空间趋势的模拟更加合理。

2)多源遥感产品间的无缝时空尺度转换。本算法通过层次贝叶斯框架下的数据模型构建，以条件概率分布的形式实现两种不同空间尺度遥感产品的无缝集成。

3)遥感产品不确定性的定量化；

4)遥感大数据量的高效计算。

附图说明

图1融合SST、MODIS SST、AMSR-E SST空间完整性评价。

图2融合SST、MODIS SST、AMSR-E SST局部方差比较。

具体实施方式

为了使本发明的目的及优点更加清楚明白，以下结合实施例对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明实施例提供了一种基于稳健的固定阶数滤波模型的多源海洋表面温度遥感产品融合方法，综合SST时空融合过程中的尺度变换问题、遥感产品的不确定性表达问题及SST时空过程复杂的空间结构和时间结构，充分利用遥感产品在空间分辨率、时空完整性、精度特征等方面的互补特性，采用基于稳健的固定阶数滤波过程模型的层次贝叶斯框架融合红外SST遥感产品和微波SST遥感产品，得到高精度的、空间完整的、局部空间模式丰富的细尺度融合数据。具体包括如下步骤：

潜在的真值SST时空过程在该时空域上的有限的卫星SST观测过程Z(s；t)与其潜在的真值过程Y(s；t)之间的关系表示为：

Z(s；t)＝Y(s；t)+ε(s；t) (1)

其中，{ε(s；t)：s∈D，t∈{1，2，...}}为时空域上的高斯白噪声过程，均值为0，方差是待估参数。在实际中，获得的是在时刻t，空间{s_1，t，...，s_n，t}上的n_t维的有效观测于此相对应，则 所有时刻的观测为条件依赖过程的高斯分布为：

S1、建立稳健的固定阶数滤波模型

潜在的真值时空过程Y(s；t)为由大尺度宏观趋势μ_t(s)、小尺度空间变异过程v(s；t)和0均值的高斯白噪声ξ(s；t)三部分组成的线性结构；

Y(s；t)＝μ(s；t)+v(s；t)+ξ(s；t) (3)

其中，μ(s；t)是一个确定性的时空均值函数，描述SST时空过程的大尺度宏观趋势，μ(s；t)＝X_t(·)′β_t，X_t(·)≡(X_1，t(·)，...，X_p，t(·))′为协变量，β_t≡(β_1，t，...，β_p，t)为未知系数；SST小尺度局部空间变异过程v(s；t)是一个时空随机过程，用时空随机效应模型模拟；ξ(s；t)为模型分解误差，描述随机变异，同地统计学中的块金效应，用0均值方差为的时空高斯白噪声过程模拟；时空随机效应子模型基于卡尔曼滤波及多分辨率小波基函数，由于SST数据除了有大量的无效观测外，还有大量的陆地像元不需要参与计算，为了在计算过程中使得矩阵满秩，本算法空间基函数的选择时增加两个限定参数——非零值的海洋小波基函数个数与非零值的小波基函数个数的比值、观测像元非零值的小波基函数个数与非零值的海洋小波基函数个数的比值，从而保证了固定阶数滤波是稳健的。

小尺度局部空间变异v(s；t)，在任一固定的时刻t，为0均值的空间随机效应模型：

v(s；t)＝S_t(s)′η_t (4)

其中，S_t(·)≡(S_1，t(·)，...，S_r，t(·))′是一组r维的细尺度下小波基函数，η_t≡(η_1，t，...，η_r，t)′是一个0均值的高斯随机过程，协方差为K_t，矩阵大小为r×r，即用r维的空间基函数S_t(·)描述任一时刻t的局部空间变异；公式(4)中，如果t＞1，且通过统计模型模拟随机变量η_t的时间依赖关系，通过一阶矢量自回归模型模拟(公式5)，在时空随机效应模型中，细尺度下的小波基函数的个数可以随着时间变化而变化，也可以是随时间不变的，如果随时间而变，则r≡max{r_t}，v(s；t)＝S(s)′η_t；

η_t+1＝H_t+1η_t+ζ_t+1；t＝1，2，.... (5)

其中，H_t+1是r×r的传播矩阵(一阶矢量自回归矩阵)，ζ_t+1是独立于η_t的r维高斯奇异矢量，均值为0，方差var(ζ_t+1)≡U_t+1；

交叉协方差K_t1，2定义为：

K_t1，t2≡cov(η_t1，η_t2)＝K_t1(H_t2H_t2-1…H_t1+1)′ (6)

相应的，K_t+1＝H_t+1K_tH′_t+1+U_t+1′

根据公式6，则时间步长为1的交叉协方差L_t+1为(公式7)：

L_t+1≡K_t，t+1＝K_tH′_t+1 (7)

综上所述，根据公式(1)，用层次模型表示潜在的真值时空过程Y(s；t)为：

η₁～N(0，K₁) (8)

η_t+1|η₁，...，η_t～N(H_t+1η_t，U_t+1)，t＝1，2，... (9)

Y(s；t)的边际分布为：

Y(s；t)协方差结构为：

其中，I(·)为指示矩阵；

根据公式(12)，可以定义出潜在的真值时空过程Y(s；t)和观测数据Z(r；u)之间的协方差结构：

cov(Y(s；t)，Z(r；u))＝cov(Y(s；t)，Y(r；u)) (13)

S2、构建层次贝叶斯融合模型

将稳健的固定阶数滤波模型作为一个SST时空过程的模拟，嵌入层次贝叶斯的框架下，表示为如下结构：

μ_a(s_a；t)＝μ(s_a；t)+v(s_a；t)+ξ(s_a；t) (17)

μ(s；t)＝X_t(·)′β_t (18)

v(s；t)＝S_t(s)′η_t (19)

其中，μ_a(s_a；t)是确定的，I是单位阵；

在层次贝叶斯框架下，SST时空演变过程μ_a(s_a；t)的先验均值的确定采用稳健的固定阶数滤波模型估计的SST时空模拟值，条件依赖子过程μ(s_a；t)、v(s_a；t)，随机变量参数为两个卫星观测数据的方差，采用逆伽马分布，即：

其中，IG表示逆伽马分布，A用来标识不同的卫星观测；形状参数q_A与尺度参数r_A均为0.1。

实施例

本算法已经成功应用于MODIS map SST产品和AMSR-E map SST数据的融合。

MODIS map SST产品：空间分辨率：4km；时间分辨率：8天合成；

AMSR-E map SST产品：空间分辨率：25km；时间分辨率：每天数据；

融合得到的数据产品：空间分辨率：4km；时间分辨率：8天合成；

MODIS map SST产品、AMSR-E map SST产品和融合得到的数据产品的空间完整性比较：

年平均有效性AMSR-E SST为87.53％，MODIS SST为80.38％，基于稳健的固定阶数滤波过程模型的层次贝叶斯融合SST为100％，实现了海洋像元的全覆盖(图1)。

融合SST的精度更接近于MODIS SST，但比MODIS SST高，比AMSR-E SST低，绝对平均偏差比MODIS低0.2205℃，比AMSR-E高0.0952℃；误差标准差比MODIS低0.0098℃，比AMSR-E高0.2104℃；均方根误差比MODIS低0.0855℃，比AMSR-E高0.2233℃(表1)。

表1 融合SST、MODIS SST、AMSR-E SST整体精度验证与比较

当MODIS有效观测时，不论AMSR-E是否为有效观测，融合SST的平均偏差明显低于MODIS SST，在AMSR-E为有效观测区域，绝对平均偏差降低0.2668℃，在AMSR-E为无效观测区域，绝对平均偏差降低0.3028，误差标准差、均方根误差比MODIS SST小，降低的幅度比平均偏差小，误差标准差分别降低0.0619℃、0.0772℃，均方根误差分别降低0.1479℃、0.1353℃。说明融合SST的精度均高于MODIS SST(表2)。

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

当MODIS为无效观测时，尽管无论在AMSR-E是否为有效观测时，融合SST存在偏差相对较大的点，在AMSR-E为无效观测的海域与其为有效观测的海域，平均偏差小了0.0566℃，但两种情况下的绝对平均偏差均在0.5℃范围内，符合GODAE融合SST绝对偏差在0.5℃范围内的要求；而误差标准差与均方根误差在AMSR-E为有效观测的区域小于AMSR-E为无效观测的区域，分别小0.1603℃、0.1290℃。总体来看，当AMSR-E SST为有效观测时，融合SST的精度高于AMSR-E SST为无效观测时的融合SST的精度(表3)。

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

局部方差的年平均值，MODIS为0.2409℃，AMSR-E为0.0562℃，融合SST为0.2134℃。图2显示了三种数据局部方差随时间的变化。从图中可以看出，融合SST局部方差明显高于AMSR-E SST，接近MODIS SST，说明融合SST保持了4km空间分辨率的MODIS的局部空间细节信息(图2)，改善了25km空间分辨率的AMSR-E局部空间细节信息。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.基于稳健的固定阶数滤波模型的多源海洋表面温度(SST)遥感产品融合方法，其特征在于，综合SST时空融合过程中的尺度转换、多源遥感产品的不确定性表达及SST时空过程复杂的空间结构和时间结构，充分利用遥感产品在空间分辨率、时空完整性、精度特征等方面的互补特性，采用基于稳健的固定阶数滤波过程模型的层次贝叶斯框架融合红外SST遥感产品和微波SST遥感产品，得到高精度的、空间完整的、局部空间模式丰富的细尺度融合数据。

2.如权利要求1所述的基于稳健的固定阶数滤波模型的多源海洋表面温度遥感产品融合方法，其特征在于，包括如下步骤：

潜在的真值SST时空过程在该时空域上的有限的卫星SST观测过程Z(s；t)与其潜在的真值过程Y(s；t)之间的关系表示为：

Z(s；t)＝Y(s；t)+ε(s；t) (1)

其中，{ε(s；t)：s∈D，t∈{1，2，...}}为时空域上的高斯白噪声过程，均值为0，方差 是待估参数。在实际中，获得的是在时刻t，空间{s_l，t，...，s_n，t}上的n_t维的有效观测于此相对应，则 所有时刻的观测为条件依赖过程的高斯分布为：

<mrow> <msubsup> <mi>Z</mi> <mi>t</mi> <mo>*</mo> </msubsup> <mo>|</mo> <msubsup> <mi>Y</mi> <mi>t</mi> <mo>*</mo> </msubsup> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <msubsup> <mi>Y</mi> <mi>t</mi> <mo>*</mo> </msubsup> <mo>,</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>&amp;epsiv;</mi> <mn>2</mn> </msubsup> <msubsup> <mi>V</mi> <mi>t</mi> <mo>*</mo> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

S1、构建稳健的固定阶数滤波模型

潜在的真值SST时空过程Y(s；t)为由大尺度宏观趋势过程μ_t(s)、小尺度空间变异过程v(s；t)和0均值的高斯白噪声ξ(s；t)三部分组成的线性结构；

Y(s；t)＝μ(s；t)+v(s；t)+ξ(s；t) (3)

其中，μ(s；t)是一个确定性的时空均值函数，描述SST时空过程的大尺度宏观趋势，μ(s；t)＝X_t(·)′β_t，其中X_t(·)≡(X_l，t(·)，...，X_p，t(·))′为协变量，β_t≡(β_l，t，...，β_p，t)为未知系数；SST小尺度局部空间变异过程v(s；t)是一个时空随机过程，用时空随机效应模型模拟；ξ(s；t)为模型分解误差，描述随机变异，同地统计学中的块金效应，用0均值、方差为的时空高斯白噪声过程模拟；

小尺度局部空间变异v(s；t)，在任一固定的时刻t，为0均值的空间随机效应模型：

v(s；t)＝S_t(s)′η_t (4)

其中，S_t(·)≡(A_l，t(·)，...，S_r，t(·))′是一组r维的细尺度下小波基函数，η_t≡(η_l，t，...，η_r，t)′是一个0均值的高斯随机过程，协方差为K_t，矩阵大小为r×r，即用r维的空间基函数S_t(·)描述任一时刻t的局部空间变异；公式(2)中，如果t＞1，通过统计模型模拟随机变量η_t的时间依赖关系，本算法采用一阶矢量自回归模型模拟(公式5)，在时空随机效应模型中，细尺度下的小波基函数的个数可以随着时间变化而变化，也可以是随时间不变的，如果随时间而变，则r≡max{r_t}，v(s；t)＝S(s)′η_t；

η_t+1＝H_t+1η_t+ζ_t+1；t＝1，2，.... (5)

其中，H_t+1是r×r的传播矩阵(一阶矢量自回归矩阵)，ζ_t+1是独立于η_t的r维高斯奇异矢量，均值为0，方差var(ζ_t+1)≡U_t+1；

交叉协方差K_t1，t2定义为：

K_t1，t2≡cov(η_t1，η_t2)＝K_t1(H_t2H_t2-1…H_t1+1)′ (6)

相应的，K_t+1＝H_t+1K_tH′_t+1+U_t+1′

根据公式6，则时间步长为1的交叉协方差L_t+1为(公式7)：

L_t+1≡K_t，t+1＝K_tH′_t+1 (7)

综上所述，根据公式(1)，用层次模型表示潜在的真值时空过程Y(s；t)为：

η₁～N(0，K₁) (8)

η_t+1|η₁，...，η_t～N(H_t+1η_t，U_t+1)，t＝1，2，... (9)

<mrow> <mi>Y</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>;</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>|</mo> <msub> <mi>&amp;eta;</mi> <mn>1</mn> </msub> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msub> <mi>&amp;eta;</mi> <mi>t</mi> </msub> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <mi>&amp;mu;</mi> <mo>(</mo> <mrow> <mi>s</mi> <mo>;</mo> <mi>t</mi> </mrow> <mo>)</mo> <mo>+</mo> <msub> <mi>S</mi> <mi>t</mi> </msub> <mo>(</mo> <mi>s</mi> <mo>)</mo> <msub> <mi>&amp;eta;</mi> <mi>t</mi> </msub> <mo>,</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>&amp;xi;</mi> <mn>2</mn> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>

Y(s；t)的边际分布为：

<mrow> <mi>Y</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>;</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <mi>&amp;mu;</mi> <mo>(</mo> <mrow> <mi>s</mi> <mo>;</mo> <mi>t</mi> </mrow> <mo>)</mo> <mo>,</mo> <msub> <mi>S</mi> <mi>t</mi> </msub> <mo>(</mo> <mi>s</mi> <mo>)</mo> <msub> <mi>K</mi> <mi>t</mi> </msub> <msub> <mi>S</mi> <mi>t</mi> </msub> <mo>(</mo> <mi>s</mi> <mo>)</mo> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>&amp;xi;</mi> <mn>2</mn> </msubsup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>

Y(s；t)协方差结构为：

<mrow> <mi>cov</mi> <mrow> <mo>(</mo> <mi>Y</mi> <mo>(</mo> <mrow> <mi>s</mi> <mo>;</mo> <mi>t</mi> </mrow> <mo>)</mo> <mo>,</mo> <mi>Y</mi> <mo>(</mo> <mrow> <mi>r</mi> <mo>;</mo> <mi>u</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>S</mi> <mi>t</mi> </msub> <mrow> <mo>(</mo> <mi>s</mi> <mo>)</mo> </mrow> <msub> <mi>K</mi> <mrow> <mi>t</mi> <mo>,</mo> <mi>u</mi> </mrow> </msub> <msub> <mi>S</mi> <mi>u</mi> </msub> <mrow> <mo>(</mo> <mi>r</mi> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>&amp;xi;</mi> <mn>2</mn> </msubsup> <mi>I</mi> <mrow> <mo>(</mo> <mi>s</mi> <mo>=</mo> <mi>r</mi> <mo>,</mo> <mi>t</mi> <mo>=</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>

其中，I(·)为指示矩阵；

根据公式(12)，可以定义出潜在的真值时空过程Y(s；t)和观测数据Z(r；u)之间的协方差结构：

cov(Y(s；t)，Z(r；u))＝cov(Y(s；t)，Y(r；u)) (13)

S2、构建层次贝叶斯融合模型

将稳健的固定阶数滤波模型作为一个SST时空过程的模拟，嵌入层次贝叶斯的框架下，表示为如下结构：

<mrow> <msub> <mi>Z</mi> <mrow> <mi>M</mi> <mi>O</mi> <mi>D</mi> <mi>I</mi> <mi>S</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>s</mi> <mi>m</mi> </msub> <mo>;</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;mu;</mi> <mi>m</mi> </msub> <mo>(</mo> <mrow> <msub> <mi>s</mi> <mi>m</mi> </msub> <mo>;</mo> <mi>t</mi> </mrow> <mo>)</mo> <mo>,</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>m</mi> <mn>2</mn> </msubsup> <mi>I</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>Z</mi> <mrow> <mi>A</mi> <mi>M</mi> <mi>S</mi> <mi>R</mi> <mo>-</mo> <mi>E</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>s</mi> <mi>a</mi> </msub> <mo>;</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;mu;</mi> <mi>a</mi> </msub> <mo>(</mo> <mrow> <msub> <mi>s</mi> <mi>a</mi> </msub> <mo>;</mo> <mi>t</mi> </mrow> <mo>)</mo> <mo>,</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>a</mi> <mn>2</mn> </msubsup> <mi>I</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>&amp;mu;</mi> <mi>m</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>s</mi> <mi>m</mi> </msub> <mo>;</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>~</mo> <mi>N</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;mu;</mi> <mi>a</mi> </msub> <mo>(</mo> <mrow> <msub> <mi>s</mi> <mi>a</mi> </msub> <mo>;</mo> <mi>t</mi> </mrow> <mo>)</mo> <mo>,</mo> <msubsup> <mi>&amp;sigma;</mi> <mrow> <mi>a</mi> <mi>m</mi> </mrow> <mn>2</mn> </msubsup> <mi>I</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>

μ_a(s_a；t)＝μ(s_a；t)+v(s_a；t)+ξ(s_a；t) (17)

μ(s；t)＝X_t(·)′β_t (18)

v(s；t)＝S_t(s)′η_t (19)

其中，μ_a(s_a；t)是确定的，I是单位阵；

在层次贝叶斯框架下，SST时空演变过程μ_a(s_a；t)的先验均值的确定采用稳健的固定阶数滤波模型估计的SST时空模拟值，条件依赖子过程μ(s_a；t)、v(s_a；t)，随机变量参数为两个卫星观测数据的方差，采用逆伽马分布，即：

<mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>A</mi> <mn>2</mn> </msubsup> <mo>~</mo> <mi>I</mi> <mi>G</mi> <mrow> <mo>(</mo> <msub> <mi>q</mi> <mi>A</mi> </msub> <mo>,</mo> <msub> <mi>r</mi> <mi>A</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>

其中，IG表示逆伽马分布，A用来标识不同的卫星观测；形状参数q_A与尺度参数r_A均为0.1。