基于海气耦合模式的均权重粒子滤波与集合卡尔曼滤波比较研究

李科，苑福利，刘厂

（1.海军研究院，天津 300061；2.哈尔滨工程大学 智能科学与工程学院，黑龙江 哈尔滨 150001）

数据同化算法已被广泛应用于大气和海洋环境的预报研究中，依靠模式状态的概率密度函数（沈哲奇等，2016），对大气和海洋信息进行准确估计。当前数据同化主要使用的模式都是非线性的（白杨等，2020）。针对模式的非线性问题，数据同化方法一般采用两种方案来解决（Chorin et al，2009），一种方法是使用高斯线性假设，如集合卡尔曼滤波数据同化方法（许立兵等，2020）。然而，实际应用的模式都是高维非线性的，针对该模式传统的线性化假设会产生截断误差，同时在海气耦合模式中也存在局限性（Evensen et al，1994）。其中最有代表性的是传统的EAKF 方法，使用蒙特卡洛方程积分控制的概率密度函数，进一步描述先验的假设和误差统计。如果模式是弱非线性时，EAKF 可以假设先验和似然密度分布都是高斯分布假设（Anderson et al，2001）。但如果系统是高度非线性的，这些迭代映射很难收敛，会影响协方差的准确性。因此，EAKF 在高度非线性的模式中具有一定的局限性（Freitas et al，2001）。针对该方法的局限性，近几年研究人员更多地将该方法与变分同化方法结合，构建集合变分方案，以En4Dvar为代表的混合变分方法（Tian et al，2018）可以较好地适应非高斯模式的要求，对于Kalman 的另一种改进方法是针对实际应用的需求构建局地化方案，其中Penny 等（2013）成功将LETKF 应用于在小尺度建模（COSMO）框架中。

另一种解决非线性问题的方法是粒子滤波方法（PF）。该方法基于统计和贝叶斯公式发展而来（Van Leeuwen et al，1996），同时结合先验和后验概率密度分布来更新粒子状态。该方法在同化过程中每个粒子将根据观测信息被赋予一个权重，该权重与观测的概率成比例，这是粒子滤波独有的特点（Van Leeuwen，2011）。但该方法也存在局限性，其中最为主要的就是“粒子退化”和“维数灾难”问题。针对粒子滤波自身存在的问题，对其改进研究也主要集中于两个方面：第一是采用混合方法，使用粒子滤波过滤状态空间的强非高斯部分，同时在高维状态使用卡尔曼滤波，该方法同时结合卡尔曼滤波与粒子滤波的优点，最有代表性的是Slivinskiet 等（2015）提出的混合粒子集合卡尔曼滤波。另一种方法是通过提议密度调整粒子防止粒子崩溃的发生，该方法的主要代表是Van Leeuwen（2011）提出了均权重粒子滤波（EWPF），该方法使用提议密度引导所有粒子靠近观测，同时提高粒子滤波的同化质量。而且Ades 等（2015） 和Chen 等（2020）分别将该方法进一步改进局地化方案并应用于高维模式中。

目前数据同化方法主要应用于耦合模式中。因此本文在一个简单的非线性耦合模型上比较两种算法的同化特点，进一步研究它们的同化特性。

1 试验方法

1.1 简单耦合模式介绍

在实际应用中，数据同化方法主要应用于耦合模式中，如CGCM 等，然而该动力学模式需要大量的计算资源支持，因此本文并没有使用复杂的CGCM 模式进行实验。Zhang 等（2013）在Lorenz 63 模型的基础上开发了简单耦合预测模型，该模型具有CGCM（Zhang et al，2012） 的基本特征，方程为：

该模型有五个变量，x1，x2和x3代表大气，w代表上层海洋，浊代表深层海洋。Zhang 等（2014）和Han 等（2014）详细介绍了参数的意义。大气模式误差为0.5，上层海洋模式误差为0.13，深海模式误差为0.01。

1.2 实验设计

本研究使用“孪生”实验框架设计实验，根据给定初始条件进行10 000 次迭代进行spin-up 扰动，然后积分10 000 次迭代产生真值。观测信息是由真值叠加白噪声近似，两个算法的集合成员数都是20。

实验框架设计参考了Zhang 等（2013）和Ades等（2015）的研究。该耦合模型初始条件为（0，1，0， 0，0），观测标准差，大气为2，上层海洋为0.5，由于深海缺乏观测，因此本实验假设深海没有观测值。根据Ades 等（2015）的实验设计，确定大气的同化观测间隔为45，上层海洋为2 356。实验积分100 000 次迭代，再根据最后的10 000 次迭代研究两种方法的同化特性。

均方根误差（root-mean-square-error，RMSE）是评价同化可靠性的重要标准，计算集合均值（x，y和z）与真值（xt，yt和zt）的差值（RMSET）：

本文采用两种均方根误差对同化结果进行评价。

1.3 集合调整卡尔曼滤波

在本文中，假设模式在时刻t 的状态为xt，观测状态为yt。基于模型方程产生观测值（Zhang et al，2011），其中观测误差假设为高斯：

1.4 均权重粒子滤波

对于数据同化的最优估计问题，可以近似理解为对于后验概率密度分布的估计。贝叶斯定理提供了一种在先验概率基础上，计算后验概率的方法。粒子滤波方法就是基于贝叶斯定理而来，它利用条件概率密度分布函数表示同化结果（Ades et al，2015）。根据观察，条件概率密度函数可以表示为（Van Leeuwen，2015）：

其中，p（x）为先验概率密度函数，p（y|x）为条件概率密度，p（x|y）为后验概率密度。

EWPF 由Van Leewen（2015）详细阐述，该方法使用提议密度保证集合中大多数粒子可以靠近观测信息，根据调整后的集合粒子确定后验概率密度：

其中子=（tj-t0）/（tn-t0）表示分析时间到观测时间的归一化时间尺度，t0和tn分别表示前一次观测时刻和后一次观测时刻。松弛项的函数B（子）是一个比例因子，控制对观测的松弛度。第二步是等权重调整，EWPF 利用松弛提议密度来调整靠近观测的粒子，保证集合中粒子的权重几乎相等，粒子权重可以表示为：

2 简单耦合模式实验结果

2.1 确定膨胀因子、松弛系数和时间尺度子

由于模态误差和采样误差的存在，对于EAKF方法中存在滤波发散的问题，使用膨胀因子避免滤波发散。而EWPF 使用时间尺度子和比例因子b控制对粒子趋向观测的松弛调整，从而调节提议密度调整粒子。因此，确定合适的膨胀因子和松弛系数可以有效地提高两种方法的同化质量。

根据上文介绍，从图1（a） 中可以看出当EAKF 中的大气膨胀因子为4 580.5 时，大气RMSE 最低。而从图1（b）可以得出，根据大气膨胀系数决定海洋膨胀系数，当EAKF 海洋膨胀系数为43.5 时，RMSE 最低。

图1 耦合模式中EAKF 中（a）大气膨胀因子和（b）海洋膨胀因子变化的RMSE 值

使用RMSE 评价EWPF 中的松弛比例因子对于同化结果的影响。由图2 可知，当EWPF 中的大气松弛参数为4，EWPF 中的海洋松弛参数为6时，大气和海洋同时RMSE 处于最低。

图2 在耦合模式中使用EWPF，松弛比例因子与RMSE 之间的关系

子的选取决定了在观测间隔中使用提议密度的时刻，针对该耦合模式中大气和海洋的不同观测频率，分别选择获得最优大气海洋状态RMSET的时间尺度子，绘制观测频率与时间尺度子的关系，如图3 所示。

图3 在耦合模式中使用EWPF 不同的观测频率与子的关系

由图3 可以看出，随着观测频率的增加，子的值逐渐降低。由公式可知，子越小表示开始使用提议密度调整粒子越早，造成该结果的主要原因是随着同化观测之间的间隔增加，为了保证粒子更好地靠近观测，需要在观测间隙更多地获取观测信息，更好地调整粒子状态，因此子的取值变小，是为了更好地调整粒子状态。

2.2 耦合模式比较实验结果

在海气耦合模式中，针对不同大气和海洋观测频率对于同化结果的影响，比较分析不同观测频率下大气和海洋状态的RMSET，结果如图4 所示。图4（a）和图4（b）给出了大气观测频率为20，海洋观测频率为1 500 时的RMSET，图4（c） 和图4（d）给出了大气观测频率为35，海洋观测频率为2 000 的结果，图4（e）和图4（f）给出了大气观测频率为45，海洋观测频率为2 356 的RMSET，由图中可以看出两种方法在不同大气观测频率下结果较为接近，随着观测频率增加，EWPF结果略优于EAKF 方法。对于上层海洋状态，EWPF 的同化结果在三种不同海洋观测条件下，其同化结果明显优于EAKF。在后续实验中将选择大气观测频率为45，海洋观测为2 356 继续实验，因为该观测频率可以更好地体现模式的非线性。

图4 两种算法的耦合模型根据不同观测频率的RMSET

在海气耦合模式中，图5（a） 描述了在海气耦合模式中选择在最后10 000 次迭代的统计结果，比较两种算法对于大气状态的估计，其中黑色菱形表示观测值，黑色点虚线表示真值。从图中可以看出EWPF（点划线）总优于EAKF（实线），而且EWPF 的轨迹更加趋近于观测值。图5（b） 给出了在60 000～100 000 模式迭代过程中两种算法对于上层海洋的状态估计，同样的EWPF 结果优于EAKF 方法。从两组实验中可以看出，EAKF 结果明显偏离真值，仅在同化时刻快速变化趋近真值，其原因是在实验中大气和海洋观测间隔选择较大，导致有效观测减少。而EWPF 可以使用提议密度来调整粒子，使其接近同化时刻前的观测，但EAKF 只在同化时刻使用观测值来调整当时的集合。因此，EAKF 状态估计在所有迭代中都有明显的突变，并在同化时迅速接近真值。为了进一步分析这一问题及其同化特征，图6 描述了两种同化方法的RMSET和RMSEO。

图5 在耦合模式下（a）给出大气x1 和（b）给出上层海洋两种方法的同化结果

根据模式状态统计结果，图6 给出了对于不同状态估计的RMSET和RMSEO，其中图6（a） 和（b）给出全时间序列的RMSET，EWPF（点划线）整体低于EAKF（实线），在同化时刻（黑点）EAKF 的RMSET明显下降，这一结果也印证了先前对EAKF 的描述，也证明了EAKF 状态估计与真值接近。图6（c）和（d）给出了两种同化方法在耦合模式下同化结果的RMSEO，可以清晰地看出EAKF 和EWPF 有相似的结果。EWPF 上层海洋低于EAKF，基于以上结果可以得出EWPF 结果更接近观测。

基于图6 对RMSE 描述，使用概率密度分布直方图进一步研究两种方法集合粒子的分布，结果如图7 所示。图7（a）和（c）给出了EAKF 中20个集合成员对于大气和海洋状态的概率分布，图7（b）和（d）给出了EWPF 中20 粒子在大气和海洋的概率密度分布，其中大气的实验选择在98 100时间步长的结果，海洋的实验选择在87 172 时间步长的结果。从图中可以得出，大部分EAKF 集合成员都分布在真值和观测值之间，而EWPF 粒子与观测值更接近。原因是EWPF 使用建议密度来调整粒子接近观测，而观测值对于最终确定粒子权重有着重要的意义。EAKF 的最终结果由观测增量决定，因此其估计值与真实值更接近。

图6 使用这两种算法在耦合模型进行不同实验条件下大气海洋对于真值和观测的RMSE 结果

图7 卡尔曼滤波和粒子滤波使用20 集合粒子数的概率密度分布图

在海气耦合模型中，假设深海没有可用的观测值，但通过公式（4）可知，深海的状态可以根据其他状态的观测来调整。图8 给出了两种同化方法在深海状态下，时间步长为92 000～100 000 时的RMSET。由于深海没有观测数据用于同化实验，因此最终结果仅由积分公式（1），通过同化调整后的大气和海洋状态来获得，从图中可以清晰地看出，EWPF 的RMSET明显低于EAKF 方法。说明在缺少观测的情况下，单纯依靠其他状态观测通过积分方程调整时，EWPF 方法可以更好地获得同化结果。

图8 深层海洋在时间序列92 000 到100 000 时间步长的RMSET

3 结论与展望

本文选择在非线性海气耦合模式上，比较分析EAKF 和EWPF 两种同化方法的同化特点。实验结果表明，在非线性耦合模式中，EWPF 模式状态估计明显优于EAKF，特别是在观测时间间隔较大的情况下，EWPF 可以有效避免EAKF 在同化时刻突然靠近模式真值，同化结果更加稳定。其次，基于两种不同的RMSE 和集合粒子分布直方图，可以看出EWPF 的同化结果更加接近观测值，而EAKF方法的同化结果更加接近状态真值。因此，EWPF同化结果的准确性与观测的质量有关。最后，在海气耦合中，假设深海中没有可用观测值，深海状态只受到积分方程中其他变量同化结果的影响，通过比较无观测深海状态RMSE 结果表明，EWPF 对于耦合模式有着更好的应用，对于缺少观测的状态，可以根据已有观测较好地调节。这是EWPF 的另一个特点，这些实验也证明了基于粒子滤波改进的EWPF 在该耦合模型中具有很好的应用潜力。

尽管EWPF 在简单模型上得到有效证明，但该方法缺少在中尺度耦合模式中与EAKF 方法的进一步验证，其特性和同化质量还需要更多的实验研究。