Search
E-mail
RSS

Online submission Editor-in-chief Peer review Staff Login

Journal of Tropical Oceanography >

2024 , Vol. 43 >Issue 1: 56 - 63

DOI: https://doi.org/10.11978/2023052

Marine Hydrology

Construction and ideal experimental verification of hybrid data assimilation method based on particle filter and 3Dvar

  • YAO Changkun ,
  • WEI Kun
Expand
  • College of Mathematical Sciences, Harbin Engineering University, Harbin 150001, China
YAO Changkun. email:

Copy editor: SUN Cuici

Received date: 2023-04-25

  Revised date: 2023-06-08

  Online published: 2023-06-26

Fold

Abstract

In this paper, a new hybrid data assimilation method is designed based on particle filter and 3Dvar. The new method generates an ensemble deviation with background error information through an optimal estimation of particle filter, thus providing flow-dependent background error covariance for 3Dvar. Particle degeneracy has always been the main obstacle of particle filtering in data assimilation field. In order to make the hybrid method work better, an improved residual resampling method is proposed to solve the problem of particle degeneracy. By sampling particles in the normal distribution, the lack of particle diversity caused by degeneracy is solved. Data assimilation experiments were tested on the ideal lorenz-63 model. The results show that the new method is better than the ETKF-3Dvar method when the model error is large, and as the model error increases, the new method is also better than the traditional data assimilation method. In the comparison experiment with hierarchical resampling and general residual resampling, the improved residual resampling method can ensure the stability of the assimilation results within a given time window, while the other two methods have a large deviation in the assimilation results.

Key words： hybrid data assimilation; particle filter; 3Dvar; residual resampling

Cite this article

YAO Changkun , WEI Kun . Construction and ideal experimental verification of hybrid data assimilation method based on particle filter and 3Dvar[J]. Journal of Tropical Oceanography, 2024 , 43(1) : 56 -63 . DOI: 10.11978/2023052

大气-海洋系统十分复杂, 各种过程在不同的时间和空间维度上具有高度的异质性(李新 等, 2007)。所以, 如何更高效地利用观测数据成为当前全球应对环境变化研究的关键课题, 数据同化(data assimilation, DA)技术因此得到了广泛的运用和长足的发展(Mclaughlin, 1995; Zhou et al, 2020)。通过数据同化技术充分利用现有数据和信息产品, 有助于模拟预测全球环境变化和开展研究, 其中海洋数据同化可为预防海洋灾害、保证航行安全、海洋生态治理等提供有效帮助(Kim et al, 2020)。就本质来说, 数据同化的目的是通过恰当的数据同化方法融合不同尺度和来源的直接或者间接的观测数据, 并利用这些数据对模型结果进行调整, 从而为下一阶段的预测提供良好的初始条件, 达到增加模型精度的目的。
目前主流的数据同化方法主要分为三大类, 分别是用动力模型作为约束条件构造变分问题, 并用伴随模式去迭代求解的变分方法, 例如三维变分(three-dimensional variational method, 3Dvar) (Dimet et al, 2016)和四维变分(four-dimensional variational method, 4Dvar)(Zeltmann et al, 2020), 和基于估计理论将蒙特卡洛(Leith et al, 1974; Wang et al, 2021)方法融入到卡尔曼滤波中诞生的集合卡尔曼滤波方法(ensemble Kalman filter method, EnKF) (Evensen et al, 1994; Solberger et al, 2020), 以及两种方法的混合方法(hybrid method)。尤其是混合方法, 如何将变分方法和集合滤波方法更好地结合, 在近些年中吸引了越来越多的关注。Hamill等(2000)给出了一种混合集合滤波方法和变分方法的原型方案, 把通过EnKF计算出的动态背景误差协方差加权代入3Dvar, 从而将流依赖(flow-dependent)的背景误差引入到3Dvar中。并指出, 如果样本足够大, 那么根据集合计算出的背景误差协方差可以完全代替3Dvar的静态背景误差协方差。Hunt等 (2004)将EnKF扩展到了四维同化形式, 设计了四维集合卡尔曼滤波方法(four-dimensional ensemble Kalman filter method, 4DEnKF), 以便EnKF能够处理同化窗口内的所有观测数据。Zupanski等 (2005)提出了最大似然集合滤波(maximum-likelihood ensemble filter, MLEF), 该方法使用集合样本计算目标函数的Hessian矩阵和梯度, 从而允许在变分方法中使用通过集合样本得到的背景误差协方差, 同时改进了非线性算子的处理方法。Wang等(2007)等应用了和Hamill类似的方案, 基于现有的WRF模型, 混合了3Dvar和集合变换卡尔曼滤波(ensemble transform Kalman filter, ETKF), 设计了一种ETKF-3Dvar方法, 并证明其性能大多数时候优于原型方案EnKF-3Dvar。Liu等(2008)设计了一种基于集合的四维变分方法(ensemble four-dimensional variational, En4Dvar), 该方法使用由集合预测结果构建的流依赖背景误差协方差, 并通过4Dvar优化得到最优估计。与标准4Dvar相比, 这种设计的一个巨大优势在于可以在公式化和实现中避免切线线性模型和伴随模型。Goodliff 等(2015)系统地比较了多种混合数据同化方法(例如ETKF-4Dvar和En4Dvar等)和两种传统方法(4Dvar和ETKF)在系统非线性较强时的性能。
以上提到的混合方法大多都是建立在变分方法和集合卡尔曼滤波方法上的, 虽然可以融合它们各自的优点, 但不能改变其本身需要满足的一些假设条件。集合卡尔曼滤波通过集合的均值和协方差来对真实值和模型的不确定性进行近似, 这其中就暗含了对模型和观测误差的高斯假设, 在误差非高斯变化的情况下不满足最优估计(Anderson, 2001)。并且, 卡尔曼滤波本身容易发生滤波发散问题(Kelly et al, 2015), 表现为随着同化时间的增加, 分析对预测的依赖程度越来越大, 观测起的作用越来越小, 从而使模型误差不断累积, 最终使同化完全失效(Lorenc et al, 2003)。
本文在Hamill给出的原型框架下(Hamill et al, 2001), 将粒子滤波(particle filter, PF)和三维变分相结合, 用粒子滤波的最优估计结果作为期望求解集合扰动, 进而得到流依赖的动态背景误差协方差。由于粒子滤波是通过一堆粒子近似后验概率密度函数, 所以理论上在粒子数量足够时, 可以完美近似任意分布, 不需要满足高斯假设(Leng et al, 2013)。同时, 由于脱离了卡尔曼滤波理论, 也不存在滤波发散问题, 取而代之的是粒子退化问题。在基本的粒子滤波中, 更新步不改变粒子的状态, 而是只更新粒子的权重, 离观测越近的粒子权重越大, 但随着观测时间的推进, 会出现一个或几个粒子的权重远大于其他粒子的情况, 导致粒子失去多样性, 从而使下一步更新失效。虽然重采样(resampling)可以缓解这种问题(冯驰 等, 2010), 但传统的重采样方法不能稳定地抑制粒子退化, 所以本文提出了一种改进的残差重采样方法, 不同于传统残差重采样方法直接复制选中粒子, 而是用选中粒子当前的状态作为均值构建正态分布, 在正态分布中采样, 实验证明可以有效避免粒子退化问题。
本文的结构安排如下: 第一部分介绍本文设计的混合数据同化方法和辅助新方法的改进残差重采样方法; 第二部分为在理想Lorenz-63模型下的数据同化实验和结果分析; 第三部分为总结。

1 基于粒子滤波和三维变分的混合数据同化方法的构建

1.1 3Dvar和PF

传统三维变分方法的目标函数为:
$\begin{align} & J\left( x \right)={{J}_{\text{b}}}\left( x \right)+{{J}_{\text{o}}}\left( x \right)= \\ & \frac{1}{2}\left[ {{\left( x-{{x}_{\text{b}}} \right)}^{\text{T}}}{{B}^{-1}}\left( x-{{x}_{\text{b}}} \right)+{{\left( y-Hx \right)}^{\text{T}}}{{R}^{-1}}\left( y-Hx \right) \right] \\ \end{align}$
其中${{J}_{\text{b}}}\left( x \right)$与${{J}_{\text{o}}}\left( x \right)$ 分别表示背景项和观测项, $x$是状态向量, ${{x}_{\text{b}}}$是背景向量, $B$是静态背景误差协方差矩阵, $y$是观测向量, $H$是观测算子, 主要作用是将模式状态向量映射或插值到观测向量, $R$是观测误差协方差矩阵。目标函数的梯度为:
$\nabla J\left( x \right)={{B}^{-1}}\left( x-{{x}_{\text{b}}} \right)+{{H}^{\text{T}}}{{R}^{-1}}\left( Hx-y \right)$
使目标函数达到最小值的状态量即为最优估计状态量, 等价于使目标函数梯度等于0的状态量。
粒子滤波脱胎于贝叶斯滤波, 基于蒙特卡洛方法和贝叶斯估计理论构建一个粒子集合来近似后验概率密度函数(王法胜 等, 2014)。由$N$个粒子构成的集合$\left\{ x_{k}^{i},w_{k}^{i} \right\}_{i=1}^{N}$近似后验概率密度函数$p\left( {{x}_{k}}|{{y}_{k}} \right)$, 其中$\left\{ x_{k}^{i},i=1,\ldots,N \right\}$表示$k$时刻的所有粒子, $\left\{ w_{k}^{i},i=1,\ldots,N \right\}$表示各粒子的权重, 且$\sum\limits_{i=1}^{N}{w_{k}^{i}=1}$。$k$时刻的后验概率密度可以近似为：
$p\left( {{x}_{k}}|{{y}_{k}} \right)=\sum\limits_{i=1}^{N}{w_{k}^{i}\delta \left( {{x}_{k}}-x_{k}^{i} \right)}$
其中$\delta \left( \centerdot \right)$为狄拉克函数, 且有
$\int_{-\infty }^{+\infty }{f}\left( x \right)\delta \left( x-a \right)dx=f\left( a \right)$
所以由(4)式, $k$时刻粒子滤波的估计值为：
$\begin{align} & {{{\hat{x}}}_{k}}=\int_{-\infty }^{+\infty }{{{x}_{k}}}p\left( {{x}_{k}}|{{y}_{k}} \right)dx \\ & =\int_{-\infty }^{+\infty }{{{x}_{k}}\sum\limits_{i=1}^{N}{w_{k}^{i}\delta \left( {{x}_{k}}-x_{k}^{i} \right)}}dx \\ & =\sum\limits_{i=1}^{n}{w_{k}^{i}}x_{k}^{i} \end{align}$
其中权重为：
$w_{k}^{i}=w_{k-1}^{i}\frac{p\left( {{y}_{k}}|x_{k}^{i} \right)p\left( x_{k}^{i}|x_{k-1}^{i} \right)}{q\left( x_{k}^{i}|x_{k-1}^{i},{{y}_{k}} \right)}$
$q\left( \centerdot \right)$表示重要性函数, 权重$w$的推导过程这里不过多赘述, 详见参考文献(王法胜 等, 2014)。
一般为了方便计算, 会直接选$q=p\left( {{x}_{k}}|{{x}_{k-1}} \right)$, 这样上式(6)权重就变为了
$w_{k}^{i}=w_{k-1}^{i}p\left( {{y}_{k}}|x_{k}^{i} \right)$
本文应用的粒子滤波器方案实现步骤如下:
1) 初始化, 在$k=0$时刻, 对由以往模型运行获取的初值${{x}_{0}}$加满足模型误差的扰动, 生成具有$N$个粒子的初始粒子集合, 并把所有粒子权重设为${1}/{N}\;$。
2) 重要性采样, 从重要性函数中抽取$N$个粒子, 由于选取$q=p\left( {{x}_{k}}|{{x}_{k-1}} \right)$, 所以这步等价于预测步, 通过给定模型预测所有粒子下一时间步的状态, 组成新的粒子集合$\left\{ x_{k}^{i} \right\}_{i=1}^{N}$。
3) 更新步, 根据式(7)更新粒子权重, 并归一化$w_{k}^{i}={w_{k}^{i}}/{\sum\limits_{i=1}^{N}{w_{k}^{i}}}\;$。
4) 根据有效粒子数${{\hat{N}}_{\text{eff}}}={1}/{\sum\limits_{i=1}^{N}{{{\left( w_{k}^{i} \right)}^{2}}}}\;$, 判断是否需要重采样, 如果小于给定阈值则执行步骤5。
5) 进行重采样。
6) 根据式(5)计算此时刻估计值。

1.2 新混合同化方法

1.2.1 融合粒子滤波和三维变分方法

新方法在原型混合方案基础上, 融合粒子滤波和三维变分方法, 引入粒子滤波的一些机制到三维变分方法中, 使得三维变分方法具有流依赖的背景误差协方差。有别于传统三维变分的目标函数, 新方法的目标函数变为如下形式:
$\begin{align} & J\left( x \right)=\frac{1}{2}\left[ {{\left( x-{{x}_{\text{b}}} \right)}^{\text{T}}}{{\left( \beta B+\left( 1-\beta \right)P \right)}^{-1}}\left( x-{{x}_{\text{b}}} \right) \right.\text{+} \\ & \left. {{\left( y-Hx \right)}^{\text{T}}}{{R}^{-1}}\left( y-Hx \right) \right] \\ \end{align}$
其中新增的$P$为由粒子滤波提供的动态背景误差协方差, 并通过加权的方式控制静态与动态背景误差协方差的比例, $\beta $为加权系数, 取值在0到1之间。
一般来说, 真实的背景误差协方差可以表示为${{P}^{t}}=\overline{\left( x-{{x}^{t}} \right){{\left( x-{{x}^{t}} \right)}^{T}}}$, 其中横线代表期望, $t$代表真实。但是在实际中, 几乎不可能得知真实状态${{x}^{t}}$, 所以对背景误差协方差只能采取估计的手段。在传统的集合卡尔曼滤波体系下, 通过用先验集合的集合平均${{\bar{x}}_{k}}$代替${{x}^{t}}$来近似得到$P$。在新方法中, 我们用由(5)式计算出的后验状态估计${{\hat{x}}_{k}}$代替${{x}^{t}}$, 用粒子集合的样本方差代替真实状态的期望。真实背景误差协方差可以近似地表示为：
$P=\frac{{{X}_{k}}X_{k}^{T}}{N-1}$
其中
${{X}_{k}}=\left[ x_{k}^{1}-{{{\hat{x}}}_{k}},\cdots,x_{k}^{N}-{{{\hat{x}}}_{k}} \right]$
$\beta $的选取与粒子数量有关, 当粒子数足够时, 样本方差足以模拟模型运行中误差的动态演化过程, 可以选择小一点的$\beta $。反之, 粒子数不足时可以适当选取较大的$\beta $。当$\beta \text{=1}$时, 新方法与传统三维变分等效。
理论上来说, 通过粒子滤波得到的${{\hat{x}}_{k}}$为从后验概率密度函数中得到的状态近似, 性能上会优于集合卡尔曼滤波类方法中通过集合平均得到的先验估计。由于EnKF-3Dvar方案中, 更新步完全依赖于卡尔曼增益的选取(吴新荣 等, 2011), 而卡尔曼增益和$P$有关, $P$的准确与否直接决定了混合同化方案的性能优劣。另外, 建立在卡尔曼滤波理论基础上的传统混合方法对模型参数的准确性要求较高, 模型参数与真实参数相差较大就有可能在同化时间推进的过程中因为误差不断堆积发生滤波发散问题, 这会直接导致之后的更新步失效。而新方法因为应用了粒子滤波代替集合卡尔曼滤波, 更新步只对粒子权重进行更新, 模型参数的不准确性直接被包含在了更新步中, 不确定性大的粒子会直接被重采样筛选掉, 减小了参数不准确给同化带来的影响。为了方便称呼, 以后称新方法为PF-3Dvar。

1.2.2 改进的残差重采样

粒子退化一直是阻碍粒子滤波方法应用的主要问题之一。主要成因有二: 第一, 受计算机算力限制, 实际应用中粒子数不能取过多, 粒子数不足就容易引发权重集中在少数粒子上。所以, 直接增加粒子数也是解决粒子退化的手段之一; 第二, 在更新步中, 粒子权重通过观测满足的概率密度函数进行更新, 常见的观测误差大多满足正态分布, 而正态分布的概率密度函数是指数型的, 有时会导致粒子权重之间相差上千倍。
除了暴力地增加粒子数外, 重采样是解决粒子退化问题最重要的手段(冯驰 等, 2010)。基本思想是去除权重小的粒子, 复制权重大的粒子。一般分层重采样利用分层统计的思想, 把采样空间分成多个部分, 每个不同的部分称为层, 然后在每个分层中进行简单随机的采样, 详见参考文献(Carpenter et al, 1999; Veeramalla et al, 2020)。一般的残差重采样步骤如下:
1) 对粒子$x_{k}^{i}$, 复制$n_{k}^{i}=\left\lfloor Nw_{k}^{i} \right\rfloor $个, 其中$\left\lfloor \centerdot \right\rfloor $为向下 取整运算符, 组成新的粒子集合$\left\{ x_{k}^{j} \right\}_{j=1}^{{{N}'}}$, 其中${N}'=\sum\limits_{i=1}^{N}{n_{k}^{i}}$, 显然$N\ge {N}'$。
2) 计算粒子残差$M=N-{N}'\ge 0$, 对新的粒子集合$\left\{ x_{k}^{j} \right\}_{j=1}^{{{N}'}}$进行随机采样$M$次, 得到$\left\{ x_{k}^{q} \right\}_{q=1}^{M}$, 合并两个粒子集合, 之后把所有粒子权重设为${1}/{N}\;$。
一般的残差重采样有两个问题: 第一, 不能保证输出的粒子数等于输入的粒子数, 还需要从粒子集合中随机采样填补缺失的粒子; 第二, 只是简单地复制权重大的粒子, 剔除权重小的粒子, 导致粒子多样性损失。针对这两个问题, 本文提出了一种改进的残差重采样方法。
对问题一, 改进残差重采样不再直接计算每个粒子生成的子代个数, 而是分别计算前$m$个粒子生成的子代总数${{N}^{m}}=\left\lfloor \sum\limits_{i=1}^{m}{Nw_{k}^{i}} \right\rfloor $, 再通过${{n}^{m}}={{N}^{m}}-{{N}^{m-1}}$计算每个粒子的子代个数, 其中$m=1,....,N$, 并规定${{N}^{0}}=0$, 这样就可以保证最后输出的粒子个数不变。对问题二, 不再单纯复制选中粒子, 而是针对每个粒子构建一个正态分布, 在保留原粒子的同时, 从分布中随机采样${{n}^{m}}-1$个新粒子, 这样就一定程度上保证了粒子的多样性。对于如何构建正态分布, 均值可以直接采用选中粒子的状态量, 而协方差的选取仍需要一些考虑。本文中一共出现了三个协方差, 分别是静态背景误差协方差$B$, 动态背景误差协方差$P$和观测误差协方差$R$, 但显然它们与粒子的正态分布应服从的协方差关系不大。本文提供了一种协方差的选取思路: 首先找一个合适的矩阵$D$作为一个单位协方差矩阵(可以参考$RPB$, 也可以凭经验选取), 并建立有效粒子数${{\hat{N}}_{\text{eff}}}$和选定矩阵之间的联系, 在有效粒子数较小时, 粒子退化严重, 可以适当选取大的协方差来增加粒子多样性, 例如$2\times D$; 反之, 如果有效粒子数很接近选定的开始重采样阈值, 可以选取小一些的协方差, 例如$0.5\times D$(图1)。这样就不会因为协方差过小而使粒子失去多样性, 或过大引入不必要的采样误差干扰同化结果。
图1 改进残差重采样方法示意图

其中圆点表示粒子, 圆点的大小表示粒子权重

Fig.1 Diagram of improved residual resampling method.

The black point represents the particle, and the size of the point represents the weight of the particle

2 利用理想实验对PF-3Dvar的验证

2.1 模型、观测、参数配置和性能指标

低阶模型的检验对于理解数据同化方法的运行逻辑是必不可少的, 因为现实中使用的预测模型往往太过复杂, 不方便从中直观地分析问题。这里, 我们选用Lorenz-63模型, 它是一组三个耦合的非线性三变量微分方程(Anderson et al, 1999)。
$\left\{ \begin{align} & \frac{d{{x}_{1}}}{dt}=-a{{x}_{1}}+a{{x}_{2}} \\ & \frac{dx{}_{2}}{dt}=-{{x}_{1}}{{x}_{3}}+b{{x}_{1}}-{{x}_{2}} \\ & \frac{d{{x}_{3}}}{dt}={{x}_{1}}{{x}_{2}}-c{{x}_{3}} \\ \end{align} \right.$
其中参数为$a=10,b=28,c=2.6$, 用四阶龙格库塔法求解, 时间步长取0.01。该模型的计算成本低, 可以进行大样本量重复实验, 并且它的连续形式是混沌的, 离散形式对初始条件有很大的敏感性。
设原始参数是真实参数, 取模型以(0, 1, 0)为真实初值运行$1\times {{10}^{4}}$个时间步的结果为真值。真值上加满足均值为0, 观测误差协方差$R$的扰动, 周期性地生成观测。观测间隔选为30, 即每隔30个时间步设置一个观测。在真实初值上加均值为0, 协方差为$B$的扰动生成初始粒子集合, 为了方便比较, 实验中所有集合滤波均采用相同的初始集合。并且, 由于粒子滤波本身对粒子数量十分依赖, 为了消除粒子数不足对实验结果产生的影响, 实验中涉及的集合滤波的集合样本个数均取为50。由于在粒子数充足的情况下, 可以选取较小的加权系数, 所以实验中$\beta $选为0.2。实验中PF-3Dvar除了2.3节外均采用本文设计的改进残差重采样方法, 开始进行重采样的有效粒子数阈值为40, 并且由于模型较为简单, 改进方法中用于构建正态分布的协方差直接采用$B$。所有实验中, 静态背景误差协方差$B$、观测误差协方差$R$和观测算子$H$采用如下形式:
$B=\left[ \begin{matrix} 8 & 4 & 4 \\ 4 & 8 & 4 \\ 4 & 4 & 8 \\ \end{matrix} \right]R=\left[ \begin{matrix} 2 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \\ \end{matrix} \right]H=\left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right]$
性能指标采用整体均方根误差(root mean squared error, RMSE) (计算舍弃前2000个时间步, 避免同化早期初始状态对轨迹的影响)
$\frac{1}{n}\sum\limits_{j=1}^{n}{{{\left[ \frac{1}{L}\sum\limits_{k=1}^{L}{{{\left( x_{k}^{j}-x_{k}^{j,t} \right)}^{2}}} \right]}^{1/2}}}$
和时刻均方根误差：
${{\left[ \frac{1}{n}\sum\limits_{i=1}^{n}{{{\left( {{x}_{i}}-x_{i}^{t} \right)}^{2}}} \right]}^{1/2}}$
其中上标$t$代表真实, $L$代表时间步数, $n$代表模型变量个数。前者用于评估同化方法在整个同化过程中的性能, 后者用于评估某一时刻的性能。

2.2 模型误差对同化结果的影响

本节通过在模型参数上加扰动模拟参数不准确对同化结果带来的影响。首先对比同为混合方法的ETKF-3Dvar和PF-3Dvar。实验中ETKF-3Dvar采用了和Wang等(2007)相同的方法防止滤波发散, 即在分析扰动上乘一个膨胀因子来缓解对分析扰动的低估。经测试最佳膨胀因子在1.8左右, 但在模型误差较大时, ETKF-3Dvar仍会产生局部滤波发散。之后对比PF-3Dvar和其他传统数据同化方法。选用三种不同的传统数据同化方法, 分别是确定性方法(deterministic methods)中极具代表性的集合平方根滤波(EnSRF, ensemble square-root filter)(Whitaker et al, 2002), 扰动观测集合卡尔曼滤波(PO-EnKF, perturbed observation EnKF) (Burgers et al, 1998)和一般的三维变分方法。为了确保统一性和公平性, 防止由于滤波发散使对照失去意义, 实验中EnSRF和PO-EnKF均采用和ETKF-3Dvar相同的方法防止发散, 并采用相同的膨胀因子。并比较四种方法在表1给出的六组扰动依次增大的模型参数下的表现。
表1 六组扰动依次增大的模型参数

Tab.1 Six groups of model parameters with successively increasing perturbations

参数 1 2 3 4 5 6
$a$ 10.1 10.6 11.1 11.6 12.1 12.6
$b$ 28.1 28.6 29.1 29.6 30.1 30.6
$c$ 2.7 3.2 3.7 4.2 4.7 5.2
首先让每个参数都增大0.1, 即a=10.1, b=28.1, c=2.7。如图2 所示, 在集合数充足, 模型误差较小的时候, ETKF-3Dvar和PF-3Dvar均有良好的效果, 其中ETKF-3Dvar和PF-3Dvar的整体均方根误差分别为1.3574和1.4257, 仅相差0.0683。PF-3Dvar整体误差略高于ETKF-3Dvar的原因可能为在模型误差较低时, 改进残差重采样方法中引入了不必要的采样误差。
图2 两种混合同化方法在模型误差较小时, 8500到9000步上的同化结果对比图

Fig. 2 Comparison of assimilation results of the two hybrid assimilation methods at 8500 to 9000 steps when the model error is relatively small

继续增大模型参数, 让$c$再增加1, 即a=10.1, b=28.1, c=3.7。如图3所示, 随着模型误差的增大, PF-3Dvar的分析场更接近真实场且更平滑, ETKF-3Dvar在8700步和8850步到8900步存在明显偏差。ETKF-3Dvar和PF-3Dvar的整体均方根误差分别为3.0324和2.3588。因为此时ETKF-3Dvar在同化过程中由于误差累积(Bishop et al, 2001; Prasad et al, 2016), 对动态背景误差协方差$P$的估计出现一定偏差, 不能充分反映误差的动态传播, 从而产生局部发散, 而PF-3Dvar摆脱了卡尔曼滤波的更新方式, 仍可以较准确地估计动态背景误差协方差。并且, 如图4所示, PF-3Dvar在每个时刻的均方根误差相比同框架下的ETKF-3Dvar更为稳定, 在同化过程中的抗干扰能力更强。
图3 两种混合同化方法在模型误差较大时, 8500到9000步上的同化结果对比图

Fig. 3 Comparison of assimilation results of the two hybrid assimilation methods at 8500 to 9000 steps when the model error is relatively large

图4 两种混合同化方法在模型误差较大时, 8500到9000步上的时刻均方根误差(RMSE)图

Fig. 4 The error of moment root mean square at 8500 to 9000 steps when the model error of the two hybrid assimilation methods is large

如图5所示, 参数组1、2由于施加扰动较小, 模型误差不明显, 改进重采样引入的采样误差对同化结果影响较大, 所以PF-3Dvar的整体误差略大于EnSRF和PO-EnKF, 但仍小于3Dvar。参数组3处PO-EnKF和PF-3Dvar整体误差基本持平, 均小于EnSRF和3Dvar。参数组4—6部分随着模型误差的进一步增大, PF-3Dvar的整体误差显著小于其他三种传统数据同化方法。
图5 PF-3Dvar和其他三种同化方法在六组不同模型参数下的整体均方根误差(RMSE)图

Fig. 5 Global root mean square error figure of PF-3Dvar and three other assimilation methods under six different model parameters

2.3 不同种重采样对PF-3Dvar方案同化结果的影响

在控制其他条件一致的情况下, 选用分层重采样和一般残差重采样与改进残差重采样进行比较。三种方法对应整体均方根误差分别为4.9410、4.9429和1.5957。如图6所示, 可以看出应用其他两种重采样方法的轨迹高度重合, 这是因为粒子退化导致所有粒子的状态都被替换为了一个权重几乎为1粒子的状态, 导致之后的更新失效。在5050步和5250步处, 另外两种重采样方法的轨迹略优于改进方法, 但在5300步后则逐渐偏离真实状态, 这是因为粒子退化虽然会导致之后的更新失效, 但对当前步的更新没有影响, 甚至有利于当前步对状态的估计, 但由于此时的粒子已经失去了多样性, 所以对任何小的扰动都十分敏感。改进残差重采样方法虽然在计算的过程中也会出现权重集中在少数几个粒子上的情况, 但因为生成的粒子不是单纯的复制, 而是在正态分布中采样, 所以仍然可以在一定程度上保证粒子的多样性。
图6 三种不同重采样方法下的PF-3Dvar同化结果对比图

Fig. 6 Comparison of assimilation results of PF-3Dvar under three different resampling methods

3 总结

本文基于Hamill提出的混合数据同化框架设计了一种新的混合数据同化方法。通过混合粒子滤波和三维变分方法, 将传统三维变分中固定不变的静态协方差和由粒子滤波提供的流依赖的动态协方差相结合估计背景误差协方差, 再通过三维变分的目标函数得到对状态的最优估计。研究发现, 在模型误差较小且集合数足够时, PF-3Dvar和同为混合方法的ETKF-3Dvar(Wang et al, 2007)相比同化效果近似, 但在模型参数不准确, 即模型误差较大时, PF-3Dvar性能优于ETKF-3Dvar, 并且具有更强的稳定性。和传统数据同化方法的比较中, 在模型误差较小时, PF-3Dvar可能因为改进残差重采样引入的采样误差使同化效果略逊于其他集合滤波, 但在模型误差超过一定大小之后, PF-3Dvar的同化效果显著优于传统数据同化方法。另外, 本文针对如何解决粒子滤波中容易发生的粒子退化问题, 提出了一种改进的残差重采样方法, 通过在一般残差重采样的复制过程中引入正态分布, 达到增大粒子多样性, 防止粒子退化的目的。实验表明, 在与其他重采样方法对比时, 改进的残差重采样方法可以有效防止粒子退化。
与此同时, 我们也注意到一些问题。粒子滤波对粒子数有一定要求, 在粒子数不足(Zhang et al, 2023)时不能取得良好的效果(当然, 这也是很多基于集合滤波混合方法的通病), 但在实际问题中, 可以取到的最大粒子数往往受到限制; 改进残差重采样方法虽然对保存粒子多样性有显著效果, 但在模型误差与观测误差本就不大的情况下, 从正态分布中采样引入的采样误差就变得十分明显, 会对同化结果有一定干扰作用; PF-3Dvar在简单模型下有较好的表现, 但在现实应用中, 状态变量维数的增加、模型复杂程度等对同化效果影响较大, 方法的局限性依然明显; 评价同化方法好坏, 除了方法精度外, 最重要的就是运行速度, 整体上PF-3Dvar的运行速度略慢于其他同化方法, 推测原因是频繁地重采样拖慢了方法的运行速度。
综上所述, 使用改进残差重采样方法的PF-3Dvar虽然仍有许多问题, 但其在模型不确定性大时能获得相较其他方法更准确和更稳定的同化结果, 这在实际使用中的优势也是十分明显的。数据同化的应用场景主要是各种复杂的陆面-大气-海洋模型, 不确定性的巨大可想而知, 所以我们有理由相信, 随着PF-3Dvar方法的不断完善, 其在条件合适的情况下将有广阔的应用前景。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
冯驰, 赵娜, 王萌, 2010. 一种改进残差重采样算法的研究[J]. 哈尔滨工程大学学报, 31(1): 120-124.

FENG CHI, ZHAO NA, WANG MENG, 2010. Improving the residual resampling algorithm[J]. Journal of Harbin Engineering University, 31(1): 120-124 (in Chinese with English abstract).

[2]
李新, 黄春林, 车涛, 2007. 中国陆面数据同化系统研究的进展与前瞻[J]. 自然科学进展, 17 (2): 163-173.

LI XIN, HUANG CHUNLIN, CHE TAO, et al, 2007. Progress and prospect of land surface data assimilation system in China[J]. Progress in Natural Science, (2): 163-173 (in Chinese with English abstract).

[3]
王法胜, 鲁明羽, 赵清杰, 等, 2014. 粒子滤波算法[J]. 计算机学报, 37(8): 1679-1694.

WANG FASHENG, LU MINGYU, ZHAO QINGJIE, et al, 2014. Particle filtering algorithm[J]. Chinese Journal of Computers, 37(8): 1679-1694 (in Chinese with English abstract).

[4]
吴新荣, 韩桂军, 李冬, 等, 2011. 集合滤波和三维变分混合数据同化方法研究[J]. 热带海洋学报, 30(6): 24-30.

DOI

WU XINRONG, HAN GUIJUN, LI DONG, et al, 2011. A hybrid ensemble filter and 3D variational analysis scheme[J]. Journal of Tropical Oceanography, 30(6): 24-30 (in Chinese with English abstract).

[5]
ANDERSON J L, 2001. An ensemble adjustment Kalman filter for data assimilation[J]. Monthly Weather Review, 129(12): 2884-2903.

DOI

[6]
ANDERSON J L, ANDERSON S L, 1999. A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts[J]. Monthly Weather Review, 127(12): 2741-2758.

DOI

[7]
BISHOP C H, ETHERTON B J, MAJUMDAR S J, 2001. Adaptive sampling with the ensemble transform Kalman filter. part I: theoretical aspects[J]. Monthly Weather Review, 129(3): 420-436.

DOI

[8]
BURGERS G, VAN LEEUWEN P J, EVENSEN G, 1998. Analysis scheme in the ensemble Kalman filter[J]. Monthly Weather Review, 126(6): 1719-1724.

DOI

[9]
CARPENTER J, CLIFFORD P, FEARNHEAD P, 1999. Improved particle filter for nonlinear problems[J]. IEE Proceedings-Radar Sonar and Navigation, 146(1): 2-7.

DOI

[10]
EVENSEN G, 1994. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[J]. Journal of Geophysical Research: Oceans, 99(C5): 10143-10162.

[11]
GOODLIFF M, AMEZCUA J, VAN LEEUWEN P J, 2015. Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[J]. Tellus A: Dynamic Meteorology and Oceanography, 67(1): 26928.

DOI

[12]
HAMILL T M, SNYDER C, 2000. A hybrid ensemble Kalman filter-3D variational analysis scheme[J]. Monthly Weather Review, 128(8): 2905-2919.

DOI

[13]
HAMILL T M, WHITAKER J S, SNYDER C, 2001. Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter[J]. Monthly Weather Review, 129(11): 2776-2790.

DOI

[14]
HUNT B R, KALNAY E, KOSTELICH E J, et al, 2004. Four-dimensional ensemble Kalman filtering[J]. Tellus A: Dynamic Meteorology and Oceanography, 56(4): 273-277.

DOI

[15]
KELLY D, MAJDA A J, TONG X T, 2015. Concrete ensemble Kalman filters with rigorous catastrophic filter divergence[J]. Proceedings of the National Academy of Sciences of the United States of America, 112(34): 10589-10594.

DOI PMID

[16]
KIM J, YOO J, DO K, 2020. Wave data assimilation to modify wind forcing using an ensemble Kalman filter[J]. Ocean Science Journal, 55(2): 231-247.

DOI

[17]
DIMET F-X L, TALAGRAND O, 2016. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[J]. Tellus A: Dynamic Meteorology and Oceanography, 38(2): 97-110.

DOI

[18]
LEITH C E, 1974. Theoretical skill of Monte Carlo forecasts[J]. Monthly Weather Review, 102(6): 409-418.

DOI

[19]
LENG HONGZE, SONG JUNQIANG, 2013. Hybrid three-dimensional variation and particle filtering for nonlinear systems[J]. Chinese Physics B, 22(3): 030505.

DOI

[20]
LIU CHENGSI, XIAO QINGNONG, WANG BIN, 2008. An ensemble-based four-dimensional variational data assimilation scheme. part I: technical formulation and preliminary test[J]. Monthly Weather Review, 136(9): 3363-3373.

[21]
LORENC A C, 2003. The potential of the ensemble Kalman filter for NWP—a comparison with 4D‐Var[J]. Quarterly Journal of the Royal Meteorological Society, 129(595): 3183-3203.

[22]
MCLAUGHLIN D, 1995. Recent developments in hydrologic data assimilation[J]. Reviews of Geophysics, 33(S2): 977-984.

[23]
PRASAD V S, JOHNY J, SODHI J S, et al, 2016. Impact of EnVar hybrid assimilation using EnKF ensembles[C]. SPIE Asia-Pacific Remote Sensing, 98820J-98820J-10.

[24]
SOLBERGER M, SPÅNBERG E, 2020. Estimating a dynamic factor model in EViews Using the Kalman filter and smoother[J]. Computational Economics, 55(1): 875-900.

[25]
VEERAMALLA S K, TALARI V K H R, 2020. Multiple dipole source localization of EEG measurements using particle filter with partial stratified resampling[J]. Biomedical engineering letters, 10(2): 205-215.

[26]
WANG JIAN, GAO XIANG, CAO RUNAN, et al, 2021. A multilevel Monte Carlo method for performing time-variant reliability analysis[J]. IEEE Access, 9: 31773-31781

[27]
WANG XUGUANG, BARKER D M, SNYDER C, et al, 2007. A hybrid ETKF-3DVAR data assimilation scheme for the WRF Model. Part I: observing system simulation experiment[J]. Monthly Weather Review, 136(12): 5116-5131.

[28]
WHITAKER J S, HAMILL M, 2002. Ensemble data assimilation without perturbed observations[J]. Monthly Weather Review, 130(7): 1913-1924.

[29]
ZELTMANN S, MINOR A, OPHUS C, 2020. Denoising of sparse three- and four-dimensional hyperspectral electron microscopy data using a total variational method[J]. Microscopy and Microanalysis, 26(S2): 1724-1726.

[30]
ZHANG QINSHENG, TAGHVAEI A, CHEN YONGXIN, 2023. An optimal control approach to particle filtering[J]. Automatica, 151: 110894.

[31]
ZHOU WEI, LI JINGHUI, XU FANGHUA, et al, 2021. The impact of ocean data assimilation on seasonal predictions based on the National Climate Center climate system model[J]. Acta Oceanologica Sinica, 40(5): 58-70.

[32]
ZUPANSKI M, 2005. Maximum likelihood ensemble filter: theoretical aspects[J]. Monthly Weather Review, 133(6): 1710-1726.

Options
Outlines

/

〈 〉