CN102737155A - 基于贝叶斯滤波的通用数据同化方法 - Google Patents

Abstract

一种基于贝叶斯滤波的通用数据同化方法：在预报步骤中将初始值集合输入分析模型获得预报集合值；在更新步骤中使用集合卡尔曼滤波计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每个预报集合；或采用粒子滤波利用集合预报值计算每个集合样本的重要性权重，利用归一化重要性权重计算有效粒子数，根据权重对集合进行重采样得到更新后的分析值和分析集合；或采用无迹卡尔曼滤波计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每个预报集合；再重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。使地球遥感观测信息和陆面过程模型信息有效融合，形成误差较小的陆面过程信息预报系统。

Description

基于贝叶斯滤波的通用数据同化方法

技术领域

本发明涉及地球系统科学信息处理领域，具体地说是涉及一种基于贝叶斯滤波的通用数据同化算法，使得地球遥感观测信息和陆面过程模型信息能够有效融合，从而能够形成误差较小的陆面过程信息(如土壤水分、土壤温度等)的预报系统。

背景技术

陆面数据同化的核心思想是在陆面过程模型的动力框架内，通过数据同化的算法融合不同来源和不同分辨率的直接与间接观测，将陆面过程模型和各种观测算子(如辐射传输模型)集成为不断地依靠观测而自动调整模型轨迹，并且减小误差的预报系统。

现代数据同化算法一般可分为连续同化和顺序同化两大类。其中，连续同化是指在一个同化窗口内，利用优化算法，通过迭代而不断调整模型初始场，最终将模型轨迹拟合到离散的观测点(多维的观测向量)上。连续同化的经典算法为变分法，常用的变分同化方法有3DVar和4DVar；顺序同化则一般利用滤波算法，在有观测的时刻，利用观测真值在误差加权的基础上对模型状态进行更新，从而获得模型状态的后验优化估计。顺序同化一般采用卡尔曼滤波和扩展卡尔曼滤波。近年来，以集合卡尔曼滤波为代表的非线性滤波方法也得到了广泛的应用，

尽管国内外学者对数据同化算法进行了大量的研究工作，但目前还存在以下不足：

1)无论是集合卡尔曼滤波还是经典的扩展卡尔曼滤波，都假设误差的先验概率分布为多维高斯分布，因此，这两种算法在解决非高斯问题的时候精度不够好；另外，需要着重指出的是服从高斯分布的变量在经过非线性系统的转换以后也会呈现非高斯分布的特征。

2)变分同化方法(例如3DVar和4DVar)一般需要求模型的伴随算子，但在陆面数据同化领域，由于大部分陆面过程模型和辐射传输模型都是非线性的，而且陆面过程模型广泛使用阈值而导致其不连续性(邱崇践，1997)，因此陆面过程模型的伴随算子在通常情况下很难得到(Reichle et al.，2002)。

3)在滤波同化方法中，扩展卡尔曼滤波作为非线性滤波领域最有名的算法(Anderson & Moore，1979)在数据同化领域得到了十分广泛的应用。但是，当这种算法应用于非线性特征明显的数据同化系统时，却遇到了很大的困难。首先EKF滤波需要对模型算子和观测算子做线性化变换，这种线性化对复杂的模型系统是非常困难的，而且，线性化只是原模型在局部的一阶近似，模型动力学的细节(高阶信息)都丢失了。另外在高维数据同化系统中对于EKF滤波，误差协方差矩阵的计算量和存储量都非常大，阻碍了EKF滤波的应用。

信号处理领域在发展非线性非高斯滤波算法方面起步早，知识积累非常丰富，有许多方法都可被陆面数据同化所借鉴。这些方法包括基于确定性采样的方法和基于随机采样的方法，其中，确定性采样的方法包括无迹卡尔曼滤波和中心差分卡尔曼滤波；随机采样的方法包括基于蒙特卡罗随机采样的采样重要性重抽样粒子滤波和无迹粒子滤波。这些新方法与集合卡尔曼滤波相比，都有一定的优势。但是，仍需要采用更精确的算法，使数据同化系统的精度不断提高。

发明内容

本发明的目的在于克服上述已有技术的不足，从贝叶斯(Bayes)估计的角度对现有非线性滤波方法进行抽象化，从理论和主要计算方法上对他们进行统一，建立基于贝叶斯滤波的通用数据同化算法。

本发明是这样实现的：

利用C++程序设计语言编程实现集合卡尔曼滤波、无迹卡尔曼滤波、和粒子滤波等贝叶斯滤波方法，由于贝叶斯滤波涉及大量矩阵运算，程序中的矩阵运算调用高性能的数学函数库BLAS(基本线性代数子例程)和LAPACK(线性代数程序包)来实现，并且利用OpenMP发展并行算法以提高系统的效率。编程实现以马其赛特旋转(Mersenne Twister)算法为核心的长周期的高性能随机数生成器，包括单变量正态分布、多变量正态分布、偏正态分布等各种概率分布函数的简单随机数生成算法；将地统计学的随机模拟方法与简单随机数生成算法相结合来构造集合卡尔曼滤波和粒子滤波算法所需的初始样本(图1)。利用函数指针的形式统一各个算法的接口，可以为以后的数据同化研究提供一个通用、快速、稳健、能够自动纠错、物理约束更强的同化系统算法平台。

本发明提供的基于贝叶斯滤波的通用数据同化方法具有三种算法，都是基于贝叶斯理论推导出来的。其中：

第一种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：

预报步骤：将初始值集合输入分析模型，获得预报集合值；

更新步骤：集合卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合；

重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

所述的基于贝叶斯滤波的通用数据同化方法，其中，预报步骤具体为：

$P_{x_{k|k-1}}=\frac{{A\left(A\right)}^T}{N-1}$

更新步骤依次为：

${\stackrel{\‾}{y}}_k=\frac{1}{N}\underset{i}{\overset{N}{\mathrm{\Σ}}}{\hat{y}}_k^i$

$P_{y_k}=\frac{{S\left(S\right)}^T}{N-1}$

$S=[{\hat{y}}_k^1-{\stackrel{\‾}{y}}_k,...,{\hat{y}}_k^N-{\stackrel{\‾}{y}}_k]$

$P_{{\mathrm{xy}}_k}=\frac{{A\left(S\right)}^T}{N-1}$

$K_k=P_{{\mathrm{xy}}_k}{[P_{y_k}+R_n]}^{-1}$

$P_{x_k}=\frac{{E\left(E\right)}^T}{N-1}$

$E=[x_k^1-{\hat{x}}_k,...,x_k^N-{\hat{x}}_k]$

${\hat{x}}_k=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^i.$

本发明提供第二种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：

预报步骤：将初始值集合输入分析模型，获得预报集合值；

更新步骤：粒子滤波利用集合预报值计算每个集合样本的重要性权重，利用归一化重要性权重计算有效粒子数，根据权重对集合进行重采样，减少权重小的样本，利用权重较大的样本进行补充，得到更新后的分析值和分析集合；粒子的权重的取值范围为0～1，权重越接近于0说明粒子的权重越小。

重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

所述的基于贝叶斯滤波的通用数据同化方法，其中，预报步骤依次为：

(1)初始化：从初始概率密度函数p(x₀)中采样得到样本

(2)预报：从转移概率密度函数

中采样得到样本

更新步骤依次为：

(3)计算重要性权重

(4)归一化重要性权重

(5)计算有效粒子数

(6)重采样：如果

小于预先给定的重采样阈值，则根据粒子的权重对其重采样，得到N个新的粒子，这些粒子的概率密度函数为p(x_k|y_k)

(7)重新设置每个粒子的权重：For i＝1，...，N，使

(8)计算滤波后的状态后验估计值：

本发明提供的第三种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：

预报步骤：将初始值集合输入分析模型，获得预报集合值；

更新步骤：无迹卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合；

重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

所述的基于贝叶斯滤波的通用数据同化方法，其中，预报步骤依次为：

(1)选择参数k，α和β的值。参数k(k≥0)确保协方差矩阵的半正定，一般情况下k＝0；参数α(0≤α≤1)控制点的散布；参数β(β≥0)，用来引入分布的高阶信息，对于高斯分布β＝2。

(2)构造一个新的状态向量χ：

${\mathrm{\χ}}_{k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k-1}& {\hat{x}}_{k-1}+\sqrt{(n_x+\mathrm{\λ})P_{x_{k-1}}}& {\hat{x}}_{k-1}-\sqrt{(n_x+\mathrm{\λ})P_{x_{k-1}}}]\end{array}\right]$

(3)预报：

${\mathrm{\χ}}_{k|k-1}^{*}=f({\mathrm{\χ}}_{k-1},v_{k-1})$

${\hat{x}}_{k|k-1}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{\mathrm{\χ}}_{i,k|k-1}^{*}$

$P_{x_{k|k-1}}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{k|k-1}){({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{k|k-1})}^T+R_v$

更新步骤依次为：

${\mathrm{\χ}}_{k|k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k|k-1}& {\hat{x}}_{k|k-1}+\sqrt{(n_x+\mathrm{\λ})P_{x_{k|k-1}}}& {\hat{x}}_{k|k-1}-\sqrt{(n_x+\mathrm{\λ})P_{x_{k|k-1}}}\end{array}\right]$

Y_k|k-1＝h(χ_k|k-1，n_k)

${\hat{y}}_{k|k-1}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{Y}_{i,k|k-1}$

$P_{y_k}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}(Y_{i,k|k-1}-{\hat{y}}_{k|k-1}){(Y_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$P_{{\mathrm{xy}}_k}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}-{\hat{x}}_{k|k-1}){(Y_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$K_k=P_{{\mathrm{xy}}_k}{[P_{y_k}+R_n]}^{-1}$

${\hat{x}}_k={\hat{x}}_{k|k-1}+K_k(y_k-{\hat{y}}_{k|k-1})$

$P_{x_k}=P_{x_{k|k-1}}+K_k{P}_{y_k}{K}_k^T$

本发明提供的基于贝叶斯滤波的通用数据同化方法，其中，随机生成100个均值为0.2，方差为0.0025的满足正太分布的随机数作为初值，输入分析模型；分析模型为Lorenz系统、陆面过程模型Noah、陆面过程模型VIC-3L和/或分布式水文模型GEOtop。

本发明具有如下效果：

本发明所实现的基于确定性采样方法的无迹卡尔曼滤波(UnscentedKalman Filter，UKF)(Julier&Uhlmann，2004)和基于蒙特卡罗随机采样的采样重要性重采样粒子滤波(Sampling Importance Resampling ParticleFilter，SIR-PF)(Arulampalam et al.，2002)。这些新方法与传统的集合卡尔曼滤波(EnsembleKalmanFilter，EnKF)相比，都有一定的优势。对于高斯分布，UKF可以捕捉到非线性函数的三阶精度的后验均值和方差，对于非高斯分布可以达到二阶精度；而且，与EnKF相比，UKF算法需要较少的样本就可以达到理想精度，所以它们的运算速度要比EnKF快得多。SIR-PF利用一系列的带有权重的粒子来近似地表达非线性系统状态的后验分布，因为SIR-PF可以近似完全地探索状态的后验空间，包含了完整的系统状态的后验信息，并且SIR-PF不需要对误差的概率分布做任何假设，误差可以是高斯也可以是非高斯的，所以它是一种有效的非线性滤波算法，从理论上来讲SIR-PF的精度比卡尔曼滤波的各种算法要高。

附图说明

图1是无迹卡尔曼滤波、集合卡尔曼滤波和粒子滤波三种算法在Lorenz系统数据同化中的应用；(a)是无迹卡尔曼滤波算法的结果；(b)是集合卡尔曼滤波算法的结果；(c)是粒子滤波算法的结果。

图2是集合卡尔曼滤波算法和粒子滤波算法在Noah模型土壤湿度同化中的应用；(a)中绿色线是集合卡尔曼滤波算法的结果、青色线是粒子滤波算法的结果、粉紫色是三维变分算法的结果、紫色是四维变分算法的结果；(b)中绿色线是集合卡尔曼滤波算法的结果、青色线是粒子滤波算法的结果、粉紫色是三维变分算法的结果、紫色是四维变分算法的结果、蓝色是土壤温度的真实值、红色是同化之前的模型模拟结果。

图3是集合卡尔曼滤波算法和粒子滤波算法在Noah模型土壤温度同化中的应用。(a)中绿色线是集合卡尔曼滤波算法的结果、青色线是粒子滤波算法的结果、粉紫色是三维变分算法的结果、紫色是四维变分算法的结果、蓝色是土壤温度的真实值、红色是同化之前的模型模拟结果；(b)中绿色线是集合卡尔曼滤波算法的结果、青色线是粒子滤波算法的结果、粉紫色是三维变分算法的结果、紫色是四维变分算法的结果、蓝色是土壤温度的真实值、红色是同化之前的模型模拟结果。

图4是集合卡尔曼滤波算法和粒子滤波算法在VIC-3L模型微波亮温同化中的应用。(a)中黑线是集合卡尔曼滤波的结果、粉紫色是粒子滤波的结果、红色是土壤湿度的真实值、蓝色是同化之前的模型模拟结果；(b)中黑线是集合卡尔曼滤波的结果、粉紫色是粒子滤波的结果、红色是土壤湿度的真实值、蓝色是同化之前的模型模拟结果。

图5是集合卡尔曼滤波算法在分布式水文模型GEOtop土壤湿度数据同化中的应用。(a)是GEOtop模型模拟的土壤湿度的空间分布；(b)利用遥感反演的土壤湿度的空间分布；(c)应用集合卡尔曼滤波算法同化以后的土壤湿度空间分布；(d)是在集合卡尔曼滤波中引入空间相关关系后的土壤湿度同化结果。

具体实施方式

首先在理论上论证贝叶斯理论是数据同化的基石：

贝叶斯理论为具有噪声的线性系统和非线性系统的顺序滤波提供了统一的方法论，因而为数据同化提供了更广义的理论基础。本发明运用数据同化的语言和标准的表达形式，从贝叶斯滤波的角度分析非线性系统中的数据同化。

1.1)数据同化与非线性动态系统

状态空间方法为描述动力学系统的状态估计问题提供了统一的框架，它分为状态预报模型和观测模型，二者在数据同化系统中又分别被称为模型算子和观测算子。

其中状态空间的非线性预报模型(即模型算子)表示为：

X^t(t_k)＝M_k(X^t(t_k-1))+η(t_k-1)

式中，X是系统的状态向量；上标t代表真值；η是独立同分布(i.i.d.)的模型误差；n是系统的维数；状态向量后的(t_k)或M的下标k表示时间；M是非线性模型，

它可能是依赖于时间的，因此使用了时间下标k，但一般在数据同化系统中，M不随时间变化。

观测模型表示为：

Y^o(t_k)＝H_k(X^t(t_k))+ε(t_k)

式中，t_k时刻的观测为Y^O；ε是独立同分布的观测误差；m是观测向量的维数；

在以上状态预报模型和观测模型中，误差的分布可以是任意形式的。

假定：(1)p(X^t(t₀))已知，即背景场的概率密度函数已知，并有背景场E(X^t(t₀))＝X^b，背景场的误差协方差矩阵Var(X^t(t₀))＝P(t₀)，或将t₀时刻的误差协方差矩阵记为背景场误差B。p(X)是模型状态的先验概率密度函数，可以看作我们根据以往的观测所积累起来的知识。

(2)模型误差和观测误差都是独立过程，而且二者相互独立，它们与初始状态也相互独立。它们的概率密度函数已知，即存在

p(η(t_k))，p(ε(t_k))

其均值为零，协方差矩阵分别用P(t_k)和R(t_k)表示。

(3)所有的概率密度函数都可以计算得到。

并且：

(1)从1到k时刻的所有观测为Y^o(t_1：k)≡{Y^o(t_i)，i＝1，...，k}；

(2)t_k时刻的真值X^t(t_k)≡X(t_k)；

(3)用概率形式表示的系统状态转换方程可表示为条件概率：p(X(t_k)|X(t_k-1))；

(4)用概率形式表示的观测方程p(Y^o(t_k)|X(t_k))，称为似然函数。它可以看作是当得到观测Y^o时，模型状态X的可能性是什么样的。似然函数建立了Y^o和X之间的依存关系，即观测算子的概率表达形式。

所谓贝叶斯滤波问题，就是在时刻k，根据时刻k之前所获得的所有观测信息，求得系统状态后验概率密度函数p(X(t_k)|Y^o(t_1：k))。

1.2)数据同化与贝叶斯估计

首先考虑一种最简单的情形，来说明贝叶斯估计与数据同化的关系。

不考虑模型状态随时间的演进，则根据贝叶斯定理，后验分布p(X|Y^o)可表示为

$p\left(X\right|Y^o)=\frac{p\left(Y^o\right|X)P\left(X\right)}{p\left(Y^o\right)}=\frac{p\left(Y^o\right|X)P\left(X\right)}{\∫p\left(Y^o\right|X)p\left(X\right)\mathrm{dX}}---\left(3\right)$

假设p(X)和p(Y^o|X)为期望为0且协方差矩阵分别为B和R的多维正态分布，展开公式(3)并对公式两边取对数，最大化后验概率密度函数，则可得到以下目标函数。

$J\left(X\right)=\frac{1}{2}{(X-X^b)}^T{B}^{-1}(X-X^b)+\frac{1}{2}{[Y^o-H\left(X\right)]}^T{R}^{-1}[Y^o-H\left(X\right)]---\left(4\right)$

它和从最佳线性无偏估计得到的3DVar目标函数完全相同，因此，可认为贝叶斯理论为数据同化提供了更广义的框架，也从更基础的数学理论上为我们揭示了数据同化的基本原理，即模型-作为一种重要的先验信息再加上观测，可提高估计的可信度，减小不确定性。数据同化实际上就是模型和观测信息的融合，可以概括为：

同化＝模型+观测

数据同化实质上就是在考虑模型误差和观测误差的基础上，融合模型动力学信息和观测信息的过程

以下对本发明作进一步详细的描述。由于本发明实现了多种算法的通用数据同化，虽然每种算法的具体实施步骤有所不同，从统一性上讲，都需要预报和更新两个步骤：

考虑模型状态随时间的演进，非线性系统的状态估计包括以下两个步骤。

(1)预报

假设已得到前一时刻的系统状态的后验概率密度函数p(X(t_k-1)|Y^o(t_1：k-1))，则根据Markov过程转移密度的Chapman-Kolmogorov方程，下一时刻状态预报的概率密度函数为

p(X(t_k)|Y^o(t_1：k-1))＝∫p(X(t_k)|X(t_k-1))p(X(t_k-1)|Y^o(t_1：k-1))dX(t_k-1)      (5)

其中，

p(X(t_k)|X(t_k-1))＝p[X(t_k)-M_k(X(t_k-1))]                (6)

因此，公式(5)可推导为

p(X(t_k)|Y^o(t_1：k-1))＝∫p[X(t_k)-M_k(X(t_k-1))]p(X(t_k-1)|Y^o(t_1：k-1))dX(t_k-1)      (7)

(2)更新

根据贝叶斯定理，系统状态的后验概率密度为

$p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k}\right))=\frac{p\left(Y^o\left(t_k\right)\right|X\left(t_k\right))p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k-1}\right))}{\∫p\left(Y^o\left(t_k\right)\right|X\left(t_k\right))p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k-1}\right))\mathrm{dX}\left(t_k\right)}---\left(8\right)$

其中，似然函数

p(Y^o(t_k)|X(t_k))＝p[Y^o(t_k)-H_k(X(t_k))]                 (9)

则，公式(8)可推导为：

$p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k}\right))=\frac{p[Y^o\left(t_k\right)-H_k\left(X\left(t_k\right)\right)]p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k-1}\right))}{\∫p[Y^o\left(t_k\right)-H_k\left(X\left(t_k\right)\right)]p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k-1}\right))\mathrm{dX}\left(t_k\right)}---\left(10\right)$

合并公式(7)和(10)得到

$p\left(X\left(t_k\right)\right|Y^o\left(t_{1:k}\right))=\frac{p[Y^o\left(t_k\right)-H_k\left(X\left(t_k\right)\right)]\∫p[X\left(t_k\right)-M_k\left(X\left(k_{k-1}\right)\right)]p\left(X\left(t_{k-1}\right)\right|Y^o\left(t_{1:k-1}\right))\mathrm{dX}\left(t_{k-1}\right)}{\∫\left\{p\right[Y^o\left(t_k\right)-H_k\left(X\left(t_k\right)\right)]\∫p[X\left(t_k\right)-M_k\left(X\left(t_{k-1}\right)\right)\left]p\left(X\left(t_{k-1}\right)\right|Y^o\left(t_{1:k-1}\right))\mathrm{dX}\left(t_{k-1}\right)\right\}\mathrm{dX}\left(t_k\right)}---\left(11\right)$

至此，公式(11)完整地表达出了顺序数据同化的贝叶斯(Bayes)递推滤波形式。它包含有3类信息：

1)系统动态演进的信息，即模型算子M_k；

2)观测信息，包括k时刻以前的所有观测Y^o(t_k)以及观测算子H_k；

3)误差信息，包括模型误差p[X(t_k)-M_k(X(t_k-1))]和观测误差p[Y^o(t_k)-H_k(X(t_k))]。公式(11)充分体现了Bayes理论对事物逐步求精的原理，针对数据同化来说，也就是在模型状态动态演进过程中，不断地融合新的观测信息，并在考虑模型和观测误差的基础上，得到模型状态的优化估计。

如果能够以解析形式给出公式(11)的解，则得到系统解。实际上，对于复杂的大系统，其概率分布函数的解析形式十分难以得到；而且即使知道它的解析解，求其统计量仍面临高维积分问题，计算十分困难。因此，在实际情况下，往往需要寻找公式(11)的次优解。基于随机采样-即蒙特卡罗(Monte Carlo)模拟的滤波方法是一种通用的数值近似算法，它可以表达任意形式的概率密度函数，对于任意非线性、非高斯的随机系统的滤波十分适用，也因此为数据同化系统中的集合预报和集合滤波提供了更广义的框架。在数据同化领域应用较多的基于蒙特卡罗(Monte Carlo)模拟的贝叶斯滤波方法是集合卡尔曼滤波，近年来，粒子滤波等其他非线性、非高斯滤波方法也越来越被重视。本发明分别在贝叶斯滤波的框架内介绍集合卡尔曼滤波、粒子滤波和无迹卡尔曼滤波的实现过程。

在介绍算法之前，首先简单介绍本文中用到的符号：u表示均值、N(u，∑)表示协方差矩阵为∑的高斯分布、和

分别表示状态x_k的分析值和预报值、

和

分别表示系统状态的分析误差协方差矩阵和预报误差协方差矩阵、

表示由

经过观测模型转换后的误差协方差矩阵、

表示状态和观测之间的交叉误差协方差矩阵、K_k表示卡尔曼增益矩阵、R_v和R_n分别表示模型和观测的误差协方差矩阵。

一、集合卡尔曼滤波的实现过程如下：

在预报部分，利用生成的初始状态向量的集合通过预报方程来获得预报场的集合，每个集合成员代表模型状态的一个具体实现，然后用预报集合计算预报误差协方差矩阵；在更新部分，EnKF利用观测向量和状态向量的误差协方差矩阵更新每个集合，得到分析场的集合。最后模型状态的后验估计值就是集合的均值。EnKF的计算过程描述如下：

a.预报：

$P_{x_{k|k-1}}=\frac{{A\left(A\right)}^T}{N-1}$

b.更新：

${\stackrel{\‾}{y}}_k=\frac{1}{N}\underset{i}{\overset{N}{\mathrm{\Σ}}}{\hat{y}}_k^i$

$P_{y_k}=\frac{{S\left(S\right)}^T}{N-1}$

$S=[{\hat{y}}_k^1-{\stackrel{\‾}{y}}_k,...,{\hat{y}}_k^N-{\stackrel{\‾}{y}}_k]$

$P_{{\mathrm{xy}}_k}=\frac{{A\left(S\right)}^T}{N-1}$

$K_k=P_{{\mathrm{xy}}_k}{[P_{y_k}+R_n]}^{-1}$

$P_{x_k}=\frac{{E\left(E\right)}^T}{N-1}$

$E=[x_k^1-{\hat{x}}_k,...,x_k^N-{\hat{x}}_k]$

${\hat{x}}_k=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^i$

其中N代表集合数，

代表k时刻预报场的第i个集合，代表k-1时刻分析场的第i个集合；

代表通过模型状态预测得到的第i个观测集合。

二、粒子滤波算法的实现过程如下：

a.预报：

(1)初始化：从初始概率密度函数p(x₀)中采样得到样本

(2)预报：从转移概率密度函数

中采样得到样本

b.更新：

(3)计算重要性权重

(4)归一化重要性权重

(5)计算有效粒子数

(6)重采样：如果小于预先给定的重采样阈值，则根据粒子的权重

对其重采样，得到N个新的粒子，这些粒子的概率密度函数为p(x_k|y_k)

(7)重新设置每个粒子的权重：For i＝1，...，N，使

(8)计算滤波后的状态后验估计值：

三、无迹卡尔曼滤波算法的实现过程如下：

UKF算法基于近似一个概率分布比近似一个任意的非线性函数更容易的观点，利用统计线性化技术来代替EKF中雅可比矩阵线性化的过程，通过一组仔细选择的确定性采样点来捕捉非线性系统的相关统计参数(如均值、方差等)。

这种确定性采样方法称为无迹转换(Unscented Transformation-UT)。UT方法的提出是为了计算一个随机变量经过非线性转换以后的统计量(如均值、方差)。n_x代表均值为

协方差为P_x的状态向量的维数，假设y是x的非线性函数，即：

y＝g(x)

为了利用UT计算y的均值和协方差，首先选取2n_x+1个带权重的点χ，表示为S＝{x_i，ω_i；i＝0，...，2n_x}，通过这些确定性选取的点，能够捕捉到向量x的真实均值和协方差。这些点的选取方式如下：

${\mathrm{\χ}}_0=\stackrel{\‾}{x}$ ${\mathrm{\ω}}_0^{\left(m\right)}=\frac{\mathrm{\λ}}{n_x+\mathrm{\λ}}i=0$

${\mathrm{\χ}}_i=\stackrel{\‾}{x}+{\left(\sqrt{(n_x+\mathrm{\λ})P_x}\right)}_i,i=1,...,n_x$ ${\mathrm{\ω}}_0^{\left(c\right)}=\frac{\mathrm{\λ}}{n_x+\mathrm{\λ}}+(1-{\mathrm{\α}}^2+\mathrm{\β}),i=0$

${\mathrm{\χ}}_i=\stackrel{\‾}{x}-{\left(\sqrt{(n_x+\mathrm{\λ})P_x}\right)}_i,i=n_x+1,....,{2n}_x$ ${\mathrm{\ω}}_i^{\left(m\right)}={\mathrm{\ω}}_i^{\left(c\right)}=\frac{1}{2(n_x+\mathrm{\λ})},i=1,...,{2n}_x$

其中λ＝α²(n_x+k)-n_x是一个尺度因子，选取的这些点称为Sigma点。权重的上标(m)和(c)代表这些权重分别被用于均值和协方差的计算。UT完整的实现过程如下：

1)选择参数k，α和β的值。参数k(k≥0)确保协方差矩阵的半正定，一般情况下k＝0；参数α(0≤α≤1)控制点的散布；参数β(β≥0)，用来引入分布的高阶信息，对于高斯分布β＝2。

2)计算2n_x+1个Sigma点及其权重

3)将每个Sigma点代入下面公式：

Y_i＝g(x_i)i＝0，...，2n_x

4)计算y的均值和协方差以及交叉协方差：

$\stackrel{\‾}{y}\≈\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{Y}_i$

$P_y\≈\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}(Y_i-\stackrel{\‾}{y}){(Y_i-\stackrel{\‾}{y})}^T$

$P_{\mathrm{xy}}\≈\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_i-\stackrel{\‾}{x}){(Y_i-\stackrel{\‾}{y})}^T$

基于上述UT转换，UKF的实现过程如下，假设噪声为加性噪声(Julier& Uhlmann，2004；van der Merwe，2004)：

根据上述Sigma点的选取方法，构造一个新的状态向量x：

${\mathrm{\χ}}_{k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k-1}& {\hat{x}}_{k-1}+\sqrt{(n_x+\mathrm{\λ})P_{x_{k-1}}}& {\hat{x}}_{k-1}-\sqrt{(n_x+\mathrm{\λ})P_{x_{k-1}}}\end{array}\right]$

a.预报：

${\mathrm{\χ}}_{k|k-1}^{*}=f({\mathrm{\χ}}_{k-1},v_{k-1})$

${\hat{x}}_{k|k-1}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{\mathrm{\χ}}_{i,k|k-1}^{*}$

$P_{x_{k|k-1}}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{k|k-1}){({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{k|k-1})}^T+R_v$

b.更新：

${\mathrm{\χ}}_{k|k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k|k-1}& {\hat{x}}_{k|k-1}+\sqrt{(n_x+\mathrm{\λ})P_{x_{k|k-1}}}& {\hat{x}}_{k|k-1}-\sqrt{(n_x+\mathrm{\λ})P_{x_{k|k-1}}}\end{array}\right]$

Y_k|k-1＝h(χ_k|k-1，n_k)

${\hat{y}}_{k|k-1}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{Y}_{i,k|k-1}$

$P_{y_k}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}(Y_{i,k|k-1}-{\hat{y}}_{k|k-1}){(Y_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$P_{{\mathrm{xy}}_k}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}-{\hat{x}}_{k|k-1}){(Y_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$K_k=P_{{\mathrm{xy}}_k}{[P_{y_k}+R_n]}^{-1}$

${\hat{x}}_k={\hat{x}}_{k|k-1}+K_k(y_k-{\hat{y}}_{k|k-1})$

$P_{x_k}=P_{x_{k|k-1}}+K_k{P}_{y_k}{K}_k^T$

四、数值试验

下面以陆面过程模型Noah中土壤湿度的数据同化为例来分别介绍两种方法的实现过程，同化中使用的集合数目为100：

a.预报：

(1)初始化：随机生成100个均值为0.2，方差为0.0025的满足正态分布的随机数作为土壤湿度的初值；

(2)预报：将随机生成的土壤湿度初值输入Noah模型，获得下一时刻的土壤湿度集合预报值；

b.更新：

·集合卡尔曼滤波：

利用土壤湿度集合预报值计算其预报误差协方差矩阵，然后通过预报误差协方差矩阵和观测误差协方差矩阵计算卡尔曼增益矩阵，最后通过土壤湿度观测值和卡尔曼增益矩阵更新每一个土壤湿度预报集合，得到更新后的土壤湿度分析值和分析集合。

·粒子滤波：

利用土壤湿度集合预报值和土壤湿度观测值计算每个集合样本的重要性权重，利用归一化重要性权重计算有效粒子数，根据权重对土壤湿度集合进行重采样，减少权重小的样本，利用权重较大的样本进行补充，得到更新后的土壤湿度分析值和分析集合。

·无迹卡尔曼滤波：

利用土壤湿度预报值计算和预报误差协方差矩阵，然后通过预报误差协方差矩阵和观测误差协方差矩阵计算卡尔曼增益矩阵，最后通过土壤湿度观测值和卡尔曼增益矩阵更新每一个土壤湿度预报集合，得到更新后的土壤湿度分析值和分析集合。

重新将更新后的土壤湿度分析集合作为Noah模型的初始值进行下一步土壤湿度预测和同化。

分别在Lorenz系统(图1)、陆面过程模型Noah(图2、3)、陆面过程模型VIC-3L(图4)和分布式水文模型GEOtop(图5)中检验了本发明所设计的通用贝叶斯滤波同化算法。

图1是三种滤波算法在Lorenz系统数据同化中的应用。在高度非线性的Lorenz系统中，无迹卡尔曼滤波、集合卡尔曼滤波和粒子滤波算法都可以很好地将观测信息融入模型中，改善了模型的运行轨迹。

图2和图3分别是集合卡尔曼滤波算法和粒子滤波算法在Noah模型土壤湿度和土壤温度同化中的应用。Noah模型土壤湿度和土壤温度同化试验的结果表明，集合卡尔曼滤波算法和粒子滤波算法可以充分利用土壤湿度和温度的观测信息来提高陆面过程模型的模拟精度，并通过模型土壤层之间的垂直相关信息，改善了下层土壤的模拟结果。并将同化结果与三维变分和四维变分算法做了对比。

图4是集合卡尔曼滤波算法和粒子滤波算法在VIC-3L模型微波亮温同化中的应用。通过集合卡尔曼滤波算法和粒子滤波算法在VIC-3L模型微波亮温中的一维同化试验，可以看到了两种同化算法在考虑垂向空间相关时，都可以将遥感观测到的表层信息向下传递，从而改善了下层土壤湿度的估计结果。

图5是集合卡尔曼滤波算法在分布式水文模型GEOtop土壤湿度数据同化中的应用。通过分布式水文模型GEOtop中的土壤湿度同化试验，可以发现集合卡尔曼滤波算法可以利用观测的水平空间相关信息，将观测信息传递到无观测区域，提高了无观测区域表层土壤湿度的估计结果。

Claims (8)

1.一种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：

预报步骤：将初始值集合输入分析模型，获得预报集合值；

更新步骤：集合卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合；

重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

2.根据权利要求1所述的基于贝叶斯滤波的通用数据同化方法，其中，预报步骤具体为：

$P_{x_{k|k-1}}=\frac{{A\left(A\right)}^T}{N-1}$

更新步骤依次为：

${\stackrel{\‾}{y}}_k=\frac{1}{N}\underset{i}{\overset{N}{\mathrm{\Σ}}}{\hat{y}}_k^i$

$P_{y_k}=\frac{{S\left(S\right)}^T}{N-1}$

$S=[{\hat{y}}_k^1-{\stackrel{\‾}{y}}_k,...,{\hat{y}}_k^N-{\stackrel{\‾}{y}}_k]$

$P_{{\mathrm{xy}}_k}=\frac{{A\left(S\right)}^T}{N-1}$

$K_k=P_{{\mathrm{xy}}_k}{[P_{y_k}+R_n]}^{-1}$

$P_{x_k}=\frac{{E\left(E\right)}^T}{N-1}$

$E=[x_k^1-{\hat{x}}_k,...,x_k^N-{\hat{x}}_k]$

${\hat{x}}_k=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^i.$

3.一种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：

预报步骤：将初始值集合输入分析模型，获得预报集合值；

更新步骤：粒子滤波利用集合预报值计算每个集合样本的重要性权重，利用归一化重要性权重计算有效粒子数，根据权重对集合进行重采样，减少权重小的样本，利用权重较大的样本进行补充，得到更新后的分析值和分析集合；粒子的权重的取值范围为0-1，权重越接近于0说明粒子的权重越小；

重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

4.根据权利要求3所述的基于贝叶斯滤波的通用数据同化方法，其中，预报步骤依次为：

(1)初始化：从初始概率密度函数p(x₀)中采样得到样本

(2)预报：从转移概率密度函数中采样得到样本

更新步骤依次为：

(3)计算重要性权重

(4)归一化重要性权重

(5)计算有效粒子数

(6)重采样：如果

小于预先给定的重采样阈值，则根据粒子的权重

对其重采样，得到N个新的粒子，这些粒子的概率密度函数为p(x_k|y_k)

(7)重新设置每个粒子的权重：For i＝1，...，N，使

(8)计算滤波后的状态后验估计值：

5.一种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：

预报步骤：将初始值集合输入分析模型，获得预报集合值；

更新步骤：无迹卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合；

重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

6.根据权利要求5所述的基于贝叶斯滤波的通用数据同化方法，其中，预报步骤依次为：

(1)选择参数k，α和β的值。参数k(k≥0)确保协方差矩阵的半正定，一般情况下k＝0；参数α(0≤α≤1)控制点的散布；参数β(β≥0)，用来引入分布的高阶信息，对于高斯分布β＝2。

(2)构造一个新的状态向量χ：

${\mathrm{\χ}}_{k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k-1}& {\hat{x}}_{k-1}+\sqrt{(n_x+\mathrm{\λ})P_{x_{k-1}}}& {\hat{x}}_{k-1}-\sqrt{(n_x+\mathrm{\λ})P_{x_{k-1}}}]\end{array}\right]$

(3)预报：

${\mathrm{\χ}}_{k|k-1}^{*}=f({\mathrm{\χ}}_{k-1},v_{k-1})$

${\hat{x}}_{k|k-1}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{\mathrm{\χ}}_{i,k|k-1}^{*}$

$P_{x_{k|k-1}}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{k|k-1}){({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{k|k-1})}^T+R_v$

更新步骤依次为：

${\mathrm{\χ}}_{k|k-1}=\left[\begin{array}{ccc}{\hat{x}}_{k|k-1}& {\hat{x}}_{k|k-1}+\sqrt{(n_x+\mathrm{\λ})P_{x_{k|k-1}}}& {\hat{x}}_{k|k-1}-\sqrt{(n_x+\mathrm{\λ})P_{x_{k|k-1}}}\end{array}\right]$

$Y_{k|k-1}=h({\mathrm{\χ}}_{k|k-1},n_k)$

${\hat{y}}_{k|k-1}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(m\right)}{Y}_{i,k|k-1}$

$P_{y_k}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}(Y_{i,k|k-1}-{\hat{y}}_{k|k-1}){(Y_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$P_{{\mathrm{xy}}_k}=\underset{i=0}{\overset{{2n}_x}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^{\left(c\right)}({\mathrm{\χ}}_{i,k|k-1}-{\hat{x}}_{k|k-1}){(Y_{i,k|k-1}-{\hat{y}}_{k|k-1})}^T$

$K_k=P_{{\mathrm{xy}}_k}{[P_{y_k}+R_n]}^{-1}$

${\hat{x}}_k={\hat{x}}_{k|k-1}+K_k(y_k-{\hat{y}}_{k|k-1})$

$P_{x_k}=P_{x_{k|k-1}}+K_k{P}_{y_k}{K}_k^T$

7.根据权利要求1、3或5所述的基于贝叶斯滤波的通用数据同化方法，其中，随机生成100个均值为0.2，方差为0.0025的满足正太分布的随机数作为初值，输入分析模型。

8.根据权利要求1、3或5所述的基于贝叶斯滤波的通用数据同化方法，其中，分析模型为Lorenz系统、陆面过程模型Noah、陆面过程模型VIC-3L和/或分布式水文模型GEOtop。

Images