CN113051520B - 一种基于统计观测均权重粒子滤波数据同化方法 - Google Patents

Abstract

本发明提供一种基于统计观测均权重粒子滤波数据同化方法，获取模式积分初始背景场；判断是否达到统计观测开始时刻；在给定统计观测开始时刻，累计求取观测均值根据统计观测计算提议密度调整集合粒子；在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态；使用重采样方法，调整集合粒子维持粒子数稳定；计算均权重粒子滤波后的状态后验估计值。本发明利用统计观测代替传统均权重粒子滤波对于未来观测的依赖，通过提议密度在同化时刻之前调整集合粒子靠近统计观测信息，进一步改善均权重粒子滤波同化方法，使用该方法可以有效提高在非高斯，非线性模式下，数据同化的同化质量，可以更好的应用于当前实时数据同化应用领域。

Description

一种基于统计观测均权重粒子滤波数据同化方法

技术领域

本发明涉及一种均权重粒子滤波数据同化方法，尤其涉及一种基于统计观测均权重粒子滤波数据同化方法，属于大气与海洋数据同化领域。

背景技术

海洋动力学研究有两种方式，一种是使用数值模型进行研究，另一种是对大气与海洋进行直接观测。数值模式模拟主要用来反映海区特性，卫星等直接观测真实反映海洋观测特点。由于卫星遥感数据的大量获取，海洋环境的特殊性以及海洋观测在空间分布上存在不足。因此如何使用这些有限且不连续的数据、并将数据结果和数值模型合理结合起来与社会发展和海洋科技发展紧密相关。因此海洋数据同化的主要目的是将观测数据与理论数值模型结果相结合，吸收两者的优点，以期得到更接近实际的结果；同化的实质是将海洋模式与观测系统相结合以达到重建天气气候全方位历史图像和模式预报初始化的过程。除了模式是基础，还有两个重要部分，一是充分全面的观测资料收集和质量控制；二是资料同化算法设计。

数据同化方法主要分为两类：统计方法和顺序方法，统计方法以三维变分和四维变分方法为主，其主要特点是利用数值最优化算法求解目标函数的最小化问题。例如专利申请号 201911195801.7的专利申请一种判断降雨预报数据同化方案优劣的方法，该专利使用一种判断降雨预报数据同化方案优劣的方法，不仅可以提高降雨预报的精度，较大程度的简化数据同化过程,而且增加多源气象信息可以弥补不同气象数据的不足,充分发挥各类数据的优势，具有普遍适用性。例如专利申请号201910158554.7的专利申请代数多重网格三维变分数据同化技术，包括多重网格同化的基本框架，非结构网格同化方程解算方法和多重非结构网格构造的基本原理等。并消除了正交网格同化结果到非规则模式网格的插值过程,有效提高了同化精度和灵活度。针对非规则网格模式的特点开发的代数多重网格数据同化技术，将大幅改善近岸海区的同化效果。例如专利申请号201910430413.6的专利申请一种基于多源观测数据的水质模型粒子滤波同化方法。该专利构建二维水质模型；初始化粒子的状态变量和参数；生成粒子的边界条件；重采样获取新的粒子集合；计算二维水质模型的模拟状态变量和参数的最优估值；将粒子的参数从t时刻递推到t+1时刻；更新时刻，继续生成粒子的边界条件，直至所有时刻运行完成,实现对二维水质模型的粒子滤波同化。利用粒子滤波算法将水质多源观测数据合理地融入二维水质模型,动态更新二维水质模型参数,提高了二维水质模型的模拟精度和预测能力。

粒子滤波数据同化的主要的优势在于在保证集合粒子状态不变的情况下，根据观测调整集合粒子权重，保证了模式的动态平衡，使该方法更适合非线性非高斯的模式下进行数据同化，但该方法也存在自身的不足，在粒子滤波中，集合粒子根据观测信息不断调节自身权重，随着权重的选择，会出现单一粒子获得全部权重，从而出现粒子退化线性，单纯采用重采样方法虽然可以有效减弱粒子退化问题，但又回引入相应的粒子贫化问题。

针对粒子滤波中出现的粒子退化和粒子贫化问题，英国教授VanLeeuwen详细总结了高维模式数据同化中的粒子滤波算法，提出了提议密度的思想，该思想通过提议密度调整集合粒子靠近未来同化时刻观测信息，在保证采样粒子数稳定的前提下，确保集合中每一个粒子都可以贡献最大的后验概率密度，进而避免粒子退化和粒子贫化问题，该方法的特点是使用较少的粒子就可以达到传统方法需要较多粒子的同化效果。但是传统均权重粒子滤波方法需要依靠未来观测信息调整粒子向观测靠近，但是在实际应用中很难获得未来观测权重满足实时数据同化的要求。

发明内容

本发明的目的是为了使用统计观测代替未来观测引导集合粒子靠近同化时刻的观测信息，使该方法具有更好的实时应用价值而提供一种基于统计观测均权重粒子滤波数据同化方法。

本发明的目的是这样实现的：

步骤一：获取模式积分初始背景场

先将初始场引入到模式方程中，先积分模式方程，使模式方程达到混沌状态，避免模式方程波动问题的出现，随之模式方程的积分，将达到混沌状态的模式变量作为初始背景场，以背景场为积分起点，进行后续的模式状态积分

步骤二：判断是否达到统计观测开始时刻

根据统计观测阈值τ确定统计观测计算的起始时刻，选择合适的τ可以有效提高统计观测的可靠性，避免引入过多或较少的观测信息，可以更好的引导集合粒子向同化时刻观测靠近。

步骤三：在给定统计观测开始时刻，累计求取历史观测均值，计算提议密度调整集合粒子数。

当达到阈值统计观测开始时，对于需要同化的观测位置，根据该位置历史观测累加求取均值代替传统方法中的未来观测信息，进一步计算集合的最优提议密度，确定集合中每个粒子对于历史观测统计均值的后验概率密度，同时调整集合粒子靠近观测且更新粒子状态。

步骤四：在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态

根据使用历史观测统计观测均值计算和提议密度调整后的集合粒子，通过均权重粒子滤波中的权重公式确定集合中粒子的权重，保证集合中的粒子可以获得较为接近且最优的权重，进一步调整集合粒子状态。

步骤五：使用重采样方法，调整集合粒子维持集合粒子数稳定

使用重采样方法，主要是为了对于使用均权重粒子滤波方法后，集合粒子权重权重表现较差的粒子进行调整，对于权重较小的粒子进行剔除，维持集合粒子数的稳定。

步骤六：计算均权重粒子滤波后的状态后验估计值

对使用重采样方法后得到的集合粒子求取集合均值，将集合均值重新引入到模式方程中，在同化过程中更新模式积分，得到最终的同化分析结果。

步骤二具体为：判断是否达到统计观测开始时刻

根据统计观测阈值τ确定统计观测计算的起始时刻，τ的计算公式为：

其中t⁰为前一个需要同化的观测时刻，tⁿ为两个观测之间的观测间隔，t^j表示当前时刻；为了保证统计观测结果的可靠性，当需要同化的观测间隔较大时τ一般选择为0.8-0.9，当同化间隔较小时τ一般选择为0.5-0.8；

步骤三具体为：在给定统计观测开始时刻，累计求取历史观测均值，计算提议密度调整集合粒子数

步骤3.1：在给定统计观测开始时刻，求解统计观测均值；

使用统计的方法对于同化时刻开始前的观测进行统计求取均值，使用

表示需要同化的观测位置上观测间隔内的观测信息的时间序列，假设根据τ确定t^j为统计观测开始时刻，则统计观测均值可以表示为：

步骤3.2：根据统计观测均值，计算粒子概率密度；

在贝叶斯理论中，后验概率密度的求解需要依靠先验概率密度分布与似然概率求解，粒子滤波方法基于贝叶斯理论使用条件后验概率密度分布，其中对于统计观测

的后验概率密度分布可以表示为

其中

为提议密度，对于提议密度的求取需要先计算集合粒子的概率密度，使用模式方程，从n-1时刻的粒子状态

进行积分，集合粒子的概率密度表达式为：

表示从n-1时刻粒子状态到n时刻的概率密度，

表示n时刻第i个粒子的状态，

表示由n-1时刻模式方程积分结果，Q表示模式误差协方差表达式为：

M表示集合粒子总数，

表示集合粒子均值；通过该公式计算粒子的概率密度；

步骤3.3：通过集合粒子的概率密度求取粒子提议密度；

基于模式方程中上一步计算得出的先验密度可以进一步求得均权重粒子中的提议密度，由先验密度求取提议密度表达式为：

公式中

表示以观测y为目标的提议密度，Kⁿ表示松弛胁迫矩阵其表达式为 Kⁿ＝QH^T(HQH^T+R)^-1，Kⁿ与模式误差协方差Q和观测误差协方差R表达式为：

表示n时刻第i个粒子的状态，H表示观测算子一般取值为1，y表示观测信息；在基于统计观测均权重粒子滤波中，提议密度的计算不仅只局限于同化时刻的观测信息，提议密度可以使用同化时刻之前根据统计观测计算观测历史均值，目的是在同化时刻之前调整集合粒子位置靠近观测信息；

步骤3.4：选择最优提议密度，计算粒子权重；

寻找最优的提议密度假设

集合中每个粒子基于统计观测的权重表示为：

其中ω_i表示集合中第i个粒子的权重，

表示t^j时刻第i个粒子的状态，

表示t^j时刻的统计观测信息，

表示概率密度，

表示提议密度，

表示t^j时刻粒子状态对于统计观测的概率密度，最后得到每个粒子经过提议密度的权重；对于提议密度调整集合粒子向统计观测靠近，集合粒子状态可以表示为：

其中

为观测算子，先验松弛系数为B(τ)，

确保集合粒子向统计观测靠近，所以松弛系数可以表示为：

B(τ)＝bτQH^TR^-1

其中τ为统计观测开始阈值，b表示控制对观测的松弛程度的比例因子；

步骤四具体为：在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态

步骤4.1：在提议密度调整的集合粒子基础上，计算粒子的权重；

根据同化中观测信息计算提议密度调整集合中粒子位置，在同化时刻使粒子更加靠近观测信息，保证大多数的粒子在同化时刻的后验概率密度函数中获得相等的权重，提议密度使集合中每个粒子靠近历史统计观测信息，从而使其获得近乎相同的权重，每个粒子的权重表达式为：

其中

表示提议密度权重在同化间隔中的累乘结果，可以确定粒子的最小权重ω_i表示为：

其中

表示的是第i个粒子的先验权重，先验权重即提议密度计算过程中每个粒子权重；y表示观测向量；H表示观测投影算子，投影算子在简单模式中H＝1；x表示状态向量；上标T表示矩阵转置；Q表示模式误差协方差矩阵；R表示观测误差协方差矩阵，假设获得集合中粒子目标权重C_i，为保证在粒子集合中80％的粒子可以达到计算得到的权重，避免粒子退化和贫化问题的发生：

从而发现对于n时刻状态

集合中多数的粒子都可以保持相同的权重，对于有些没有达到均权重的粒子可以通过重采样方法进行调整；

步骤4.2：根据权重调整粒子状态；

在获得集合粒子权重后，根据权重调整粒子状态，则n时刻的状态可以表示为：

公式中y表示观测向量；H表示观测投影算子(H＝1)；x表示状态向量； K＝QH^T(HQH^T+R)^-1，Q是模式误差协方差，R是观测误差协方差，在均权重粒子滤波中集合粒子权重近似相等时，对于矢量α_i可以表示为：

式中

和

在这两个表达式中

C表示选择的目标权重，

表示在当前时刻粒子的相对权重；为了保证集合中粒子状态的随机效应加入随机项，得到最终的分析方程：

步骤六具体为：计算均权重粒子滤波后的状态后验估计值；

根据重采样后的集合粒子的状态，计算状态的后验估计集合均值：

将更新后的后验估计集合均值作为分析模型的初始值，重新带入模式积分方程中进行下一步预测和同化，在有可用观测的同化时间内重复上述步骤，得到最终的分析场，该分析场可以作为反映当前环境状态的数据场。

与现有技术相比，本发明的有益效果是：

(1)在保证均权重粒子滤波方法对于粒子滤波同化方法的改善能力的同时，对于解决观测间隔较大和未来观测难以实时获取的情况，保证使用较少的集合粒子可以达到较为稳定的同化结果。

(2)使用统计观测有效的代替传统方法对于未来观测的依赖，可以使该方法更好的应用于实时数据同化的应用需求，同时可以提高同化质量,使用统计观测的均权重粒子滤波方的同化结果均方根误差优于传统均权重粒子滤波方法，均方根对比结果如图所示。

附图说明

图1是统计观测均权重粒子滤波流程；

图2a-图 2 b是统计观测均权重粒子滤波与均权重粒子滤波方法均方根误差对比图；

图3是大气海洋数据同化流程图；

图4是传统均权重粒子滤波数据同化方法；

图5是基于统计观测均权重粒子滤波数据同化流程图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述。

为了更简单明确的描述基于统计观测均权重粒子滤波同化方法具体实施步骤，以简单的 Lorenz-63模式为例进行简单的说明。

步骤一：获取模式积分初始背景场

在Lorenz-63模式中先将初始场(0，1，0)输入模型方程积分100万步进行spin-up，使模式变量达到混沌状态，避免模式在同化积分过程中引入模式波动偏差，将达到混沌状态的模式变量作为模式后续积分的初始背景场，以此模式背景场为基础，通过模式方程获得粒子滤波集合粒子的初始值。可以有效提高模式方程求解得到的集合粒子状态的可靠性。

步骤二：判断是否达到统计观测开始时刻

在均权重粒子滤波数据同化中，使用统计历史观测计算提议密度，引导粒子更好的向需要同化的观测位置的历史统计观测信息靠近，提高采样过程中有效粒子数，摆脱对于未来观测的依赖，选择合适时刻开始计算统计观测具体计算方法如下：

对于统计历史观测的计算首先要确定统计观测开始阈值τ，它决定什么时刻开始进行观测统计，τ的选择与需要同化的观测信息的时间间隔有关，当τ过小时，将会将更多无用观测引入到统计当中，提高统计误差，当τ过大的时候，统计的观测信息会太少，无法正确的计算提议密度引导粒子靠近观测，τ的计算公式为：

其中t⁰为前一个需要同化的观测时刻，tⁿ为两个观测之间的观测间隔，t^j表示当前时刻。为了保证统计观测结果的可靠性，当需要同化的观测间隔较大时τ一般选择为0.8-0.9，当同化间隔较小时τ一般选择为0.5-0.8。

步骤三：在给定统计观测开始时刻，累计求取历史观测均值，计算提议密度调整集合粒子数

步骤3.1：在给定统计观测开始时刻，求解统计观测均值；

在数据同化的过程中，对于观测的获取是多种多样的，随着时间的推移观测信息源源不断的产生，在数据同化过程中仅考虑需要同化时刻的观测信息，但是同化时刻前该位置的历史观测信息往往被忽略，因此使用统计的方法对于同化时刻开始前的观测进行统计求取均值，使用

表示需要同化的观测位置上观测间隔内的观测信息的时间序列，假设根据τ确定t^j为统计观测开始时刻，则统计观测均值可以表示为：

步骤3.2：根据统计观测均值，计算粒子概率密度；

在贝叶斯理论中，后验概率密度的求解需要依靠先验概率密度分布与似然概率求解，粒子滤波方法基于贝叶斯理论使用条件后验概率密度分布，其中对于统计观测

的后验概率密度分布可以表示为

其中

为提议密度，对于提议密度的求取需要先计算集合粒子的概率密度，使用模式方程，从n-1时刻的粒子状态

进行积分，集合粒子的概率密度表达式为：

表示从n-1时刻粒子状态到n时刻的概率密度，

表示n时刻第i个粒子的状态，

表示由n-1时刻模式方程积分结果，Q表示模式误差协方差表达式为：

M表示集合粒子总数，

表示集合粒子均值。通过该公式计算粒子的概率密度。

步骤3.3：通过集合粒子的概率密度求取粒子提议密度；

基于模式方程中上一步计算得出的先验密度可以进一步求得均权重粒子中的提议密度，由先验密度求取提议密度表达式为：

公式中

表示以观测y为目标的提议密度，Kⁿ表示松弛胁迫矩阵其表达式为

Kⁿ与模式误差协方差Q和观测误差协方差R表达式为：

表示n时刻第i个粒子的状态，H表示观测算子一般取值为1，y表示观测信息。在基于统计观测均权重粒子滤波中，提议密度的计算不仅只局限于同化时刻的观测信息，提议密度可以使用同化时刻之前根据统计观测计算观测历史均值，目的是在同化时刻之前调整集合粒子位置靠近观测信息。

步骤3.4：选择最优提议密度，计算粒子权重；

提议密度的选择被认为是控制粒子与观测信息位置和粒子权重计算的重要标准，选择合适最优的提议密度可以保证采样中有效粒子的数量，同时保证粒子权重的可靠性，也是均权重粒子滤波中最重要的一部分，寻找最优的提议密度假设

集合中每个粒子基于统计观测的权重表示为：

其中ω_i表示集合中第i个粒子的权重，

表示t^j时刻第i个粒子的状态，

表示t^j时刻的统计观测信息，

表示概率密度，

表示提议密度，

表示t^j时刻粒子状态对于统计观测的概率密度，最后得到每个粒子经过提议密度的权重。对于提议密度调整集合粒子向统计观测靠近，集合粒子状态可以表示为：

其中

为观测算子，先验松弛系数为B(τ)，

确保集合粒子向统计观测靠近，所以松弛系数可以表示为：

B(τ)＝bτQH^TR^-1

其中τ为统计观测开始阈值，b表示控制对观测的松弛程度的比例因子。

步骤四：在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态

步骤4.1：在提议密度调整的集合粒子基础上，计算粒子的权重；

根据同化中观测信息计算提议密度调整集合中粒子位置，在同化时刻使粒子更加靠近观测信息，保证大多数的粒子在同化时刻的后验概率密度函数中获得相等的权重。提议密度使集合中每个粒子靠近历史统计观测信息，从而使其获得近乎相同的权重。每个粒子的权重表达式为：

其中

表示提议密度权重在同化间隔中的累乘结果，可以确定粒子的最小权重ω_i表示为：

其中

表示的是第i个粒子的先验权重，先验权重即提议密度计算过程中每个粒子权重；y 表示观测向量；H表示观测投影算子，投影算子在简单模式中H＝1；x表示状态向量；上标T表示矩阵转置；Q表示模式误差协方差矩阵；R表示观测误差协方差矩阵。假设获得集合中粒子目标权重C_i，为保证在粒子集合中80％的粒子可以达到计算得到的权重，避免粒子退化和贫化问题的发生：

从而可以发现对于n时刻状态

集合中多数的粒子都可以保持相同的权重，对于有些没有达到均权重的粒子可以通过重采样方法进行调整。

步骤4.2：根据权重调整粒子状态；

在获得集合粒子权重后，根据权重调整粒子状态，则n时刻的状态可以表示为：

公式中y表示观测向量；H表示观测投影算子(H＝1)；x表示状态向量； K＝QH^T(HQH^T+R)^-1，Q是模式误差协方差，R是观测误差协方差，在均权重粒子滤波中集合粒子权重近似相等时，对于矢量α_i可以表示为：

中

和

在这两个表达式中

C表示选择的目标权重，

表示在当前时刻粒子的相对权重。为了保证集合中粒子状态的随机效应加入随机项，得到最终的分析方程：

步骤五：使用重采样方法，调整集合粒子维持粒子数稳定；

重采样算法的基本思想是去除小权重的粒子，为保证集合中粒子总数的稳定，对于权重较大的粒子进行复制。均权重粒子滤波方法的本身保证大部分粒子得以保存，但是为了防止在均权重过程中，有极少数的粒子没有达到均等权重的要求，在均权重粒子滤波对粒子状态进行调整后，根据粒子权重使用重采样方法进行最后的状态调节，维持集合粒子数稳定。

步骤六：计算均权重粒子滤波后的状态后验估计值；

根据重采样后的集合粒子的状态，计算状态的后验估计集合均值：

将更新后的后验估计集合均值作为分析模型的初始值，重新带入模式积分方程中进行下一步预测和同化，在有可用观测的同化时间内重复上述步骤，得到最终的分析场，该分析场可以作为反映当前环境状态的数据场。

本发明公开了本发明的一种基于统计观测均权重粒子滤波数据同化技术，具体包括以下几个步骤：步骤一：获取模式积分初始背景场；步骤二：判断是否达到统计观测开始时刻；步骤三：在给定统计观测开始时刻，累计求取观测均值根据统计观测计算提议密度调整集合粒子；步骤四：在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态；步骤五：使用重采样方法，调整集合粒子维持粒子数稳定；步骤六：计算均权重粒子滤波后的状态后验估计值。本发明利用统计观测代替传统均权重粒子滤波对于未来观测的依赖，通过提议密度在同化时刻之前调整集合粒子靠近统计观测信息，进一步改善均权重粒子滤波同化方法，使用该方法可以有效提高在非高斯，非线性模式下，数据同化的同化质量，可以更好的应用于当前实时数据同化应用领域。

Claims (1)

1.一种基于统计观测均权重粒子滤波数据同化方法，其特征在于，包括如下步骤：

步骤一：获取模式积分初始背景场

先将初始场引入到模式方程中，先积分模式方程，使模式方程达到混沌状态，避免模式方程波动问题的出现，随之模式方程的积分，将达到混沌状态的模式变量作为初始背景场，以背景场为积分起点，进行后续的模式状态积分；

步骤二：判断是否达到统计观测开始时刻

根据统计观测阈值τ确定统计观测计算的起始时刻，τ的计算公式为：

其中t⁰为前一个需要同化的观测时刻，tⁿ为两个观测之间的观测间隔，t^j表示当前时刻；为了保证统计观测结果的可靠性，当需要同化的观测间隔较大时τ选择为0.8-0.9，当同化间隔较小时τ选择为0.5-0.8；

步骤三：在给定统计观测开始时刻，累计求取历史观测均值，计算提议密度调整集合粒子数

步骤3.1：在给定统计观测开始时刻，求解统计观测均值；

使用统计的方法对于同化时刻开始前的观测进行统计求取均值，使用

表示需要同化的观测位置上观测间隔内的观测信息的时间序列，假设根据τ确定t^j为统计观测开始时刻，则统计观测均值表示为：

步骤3.2：根据统计观测均值，计算粒子概率密度；

在贝叶斯理论中，后验概率密度的求解需要依靠先验概率密度分布与似然概率求解，粒子滤波方法基于贝叶斯理论使用条件后验概率密度分布，其中对于统计观测

的后验概率密度分布表示为

其中

为提议密度，对于提议密度的求取需要先计算集合粒子的概率密度，使用模式方程，从n-1时刻的粒子状态

进行积分，集合粒子的概率密度表达式为：

表示从n-1时刻粒子状态到n时刻的概率密度，

表示n时刻第i个粒子的状态，

表示由n-1时刻模式方程积分结果，Q表示模式误差协方差表达式为：

M表示集合粒子总数，

表示集合粒子均值；通过该公式计算粒子的概率密度；

步骤3.3：通过集合粒子的概率密度求取粒子提议密度；

基于模式方程中上一步计算得出的先验密度进一步求得均权重粒子中的提议密度，由先验密度求取提议密度表达式为：

公式中

表示以统计观测y为目标的提议密度，Kⁿ表示松弛胁迫矩阵其表达式为Kⁿ＝QH^T(HQH^T+R)^-1，Kⁿ与模式误差协方差Q和观测误差协方差R表达式为：

表示n时刻第i个粒子的状态，H表示观测算子取值为1，y表示观测；在基于统计观测均权重粒子滤波中，提议密度的计算不仅只局限于同化时刻的观测信息，提议密度使用同化时刻之前根据统计观测计算观测历史均值，目的是在同化时刻之前调整集合粒子位置靠近观测信息；

步骤3.4：选择最优提议密度，计算粒子权重；

寻找最优的提议密度假设

集合中每个粒子基于统计观测的权重表示为：

其中ω_i表示集合中第i个粒子的权重，

表示t^j时刻第i个粒子的状态，

表示t^j时刻的统计观测信息，

表示概率密度，

表示提议密度，

表示t^j时刻粒子状态对于统计观测的概率密度，最后得到每个粒子经过提议密度的权重；对于提议密度调整集合粒子向统计观测靠近，集合粒子状态表示为：

其中

为观测算子，先验松弛系数为B(τ)，

确保集合粒子向统计观测靠近，所以松弛系数表示为：

B(τ)＝bτQH^TR^-1

其中τ为统计观测开始阈值，b表示控制对观测的松弛程度的比例因子；

步骤四：在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态

步骤4.1：在提议密度调整的集合粒子基础上，计算粒子的权重；

根据同化中观测信息计算提议密度调整集合中粒子位置，在同化时刻使粒子更加靠近观测信息，保证大多数的粒子在同化时刻的后验概率密度函数中获得相等的权重，提议密度使集合中每个粒子靠近历史统计观测信息，从而使其获得近乎相同的权重，每个粒子的权重表达式为：

其中

第i个粒子的提议密度权重在同化间隔中的累乘结果，确定粒子的最小权重ω_i表示为：

其中

第i个粒子的提议密度权重在同化间隔中的累乘结果，先验权重即提议密度计算过程中每个粒子权重；y表示观测向量；H表示观测投影算子，投影算子在简单模式中H＝1；x表示状态向量；上标T表示矩阵转置；Q表示模式误差协方差；R表示观测误差协方差，假设获得集合中粒子目标权重C_i，为保证在粒子集合中80％的粒子达到计算得到的权重，避免粒子退化和贫化问题的发生：

从而发现对于n时刻状态

集合中多数的粒子都保持相同的权重，对于有些没有达到均权重的粒子通过重采样方法进行调整；

步骤4.2：根据权重调整粒子状态；

在获得集合粒子权重后，根据权重调整粒子状态，则n时刻的状态表示为：

公式中y表示观测向量；H表示观测投影算子，H＝1；x表示状态；K＝QH^T(HQH^T+R)^-1，Q是模式误差协方差，R是观测误差协方差，在均权重粒子滤波中集合粒子权重近似相等时，对于矢量α_i表示为：

式中

和

在这两个表达式中

C表示选择的目标权重，

第i个粒子的提议密度权重在同化间隔中的累乘结果；为了保证集合中粒子状态的随机效应加入随机项β，得到最终的分析方程：

步骤五：使用重采样方法，调整集合粒子维持粒子数稳定；

步骤六：计算均权重粒子滤波后的状态后验估计值；

根据重采样后的集合粒子的状态，计算状态的后验估计集合均值：

将更新后的后验估计集合均值作为分析模型的初始值，重新带入模式积分方程中进行下一步预测和同化，在有可用观测的同化时间内重复上述步骤一到步骤六，得到最终的分析场，该分析场作为反映当前环境状态的数据场。

Images