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

Abstract

本发明提供一种基于统计观测局地化均权重粒子滤波数据同化方法，获取模式积分初始背景场；判断是否达到统计观测开始时刻，累加观测求取统计观测均值；根据统计观测均值计算提议密度调整集合粒子；在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态；使用重采样方法，调整集合粒子维持粒子数稳定，更新观测对应位置粒子状态；使用局地化函数，确定同化观测对应位置周围的粒子权重；根据局地化权重更新粒子权重，更新周围粒子状态；计算统计观测局地化均权重粒子滤波的状态后验估计值。本发明可以有效提高非高斯网格化模式的数据同化质量，可以更好的应用于实时数据同化在网格化复杂模式中，提高同化质量。

Description

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

技术领域

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

背景技术

海洋动力学研究有两种方式，一种是使用数值模型进行研究，另一种是对大气与海洋进行直接观测。数值模式模拟主要用来反映海区特性，卫星等直接观测真实反映海洋观测特点。由于卫星遥感数据的大量获取，海洋环境的特殊性以及海洋观测在空间分布上存在不足。数据同化是一种能够有机结合数值模式和观测这两种海洋学研究基本手段的研究方法。数据同化是指在考虑数据时空分布以及观测场和背景场误差的基础上，在动力学模型的动态运行过程中不断融入新的观测数据的方法。通过在模型中不断融入新的观测数据，可以逐渐校正模型模拟预测的轨迹，使之更加接近真实的轨迹，提高模型模拟预测精度。数据同化的主要目的是将观测数据与理论模型结果相结合，吸收两者的优点，以期得到更接近实际的结果；目前数据同化方法已被广泛应用于大气、海洋和陆面等领域，为海洋模式状态的预报提供更加准确的初始场并优化海洋模式参数，以提高海洋模式的气候预报能力。

数据同化算法作为数据同化的核心，主要依赖于准确的观测数据和合理的数值模型。按同化算法与模型之间的关联性，数据同化算法分为连续数据同化算法和序贯数据同化算法两大类。例如专利申请号201910038258.3的专利申请基于最优观测时间窗口的耦合数据同化与参数优化方法，该专利使用基于最优观测时间窗口的耦合数据同化与参数优化方法，属于耦合气候模式系统的数据同化、参数优化与数值预报技术领域。最大程度上提取有效的观测信息以拟合耦合模式状态的特征变率并忽略模式内部参数的时变特征并引入时间窗口内的时间平均系数,实现对模式参数的更加精确的估计与优化,强化耦合模式的大气与海洋的数值预报能力。例如专利申请号202010013183.6的专利申请卫星数据同化在垂直方向的适应性局地化方法及集合卡曼滤波天气同化预报方法，该专利使用卫星数据同化在垂直方向的适应性局地化方法及集合卡曼滤波天气同化预报方法。适应性局地化方法根据集合卡曼滤波同化系统中给出的任意观测资料和模式变量，计算出观测资料和模式变量的相关系数；接着利用分组后的相关系数估计该观测资料和模式变量的原始局地化函数；根据相关系数的廓线估计出卫星观测的位置，将所得到的适应性局地化参数用于区域模式中预报台风,预报结果与没有使用本发明的预报结果相比,相对于观测的误差明显减小,同时使用本发明还明显改进了台风快速增强阶段的预报。例如专利申请号201910430413.6的专利申请一种基于多源观测数据的水质模型粒子滤波同化方法。该专利构建二维水质模型；初始化粒子的状态变量和参数；生成粒子的边界条件；重采样获取新的粒子集合；计算二维水质模型的模拟状态变量和参数的最优估值；将粒子的参数从t时刻递推到t+1时刻；更新时刻，继续生成粒子的边界条件，直至所有时刻运行完成,实现对二维水质模型的粒子滤波同化。利用粒子滤波算法将水质多源观测数据合理地融入二维水质模型,动态更新二维水质模型参数,提高了二维水质模型的模拟精度和预测能力。

粒子滤波算法是一种集合数据同化方法，由于粒子滤波算法不受模型状态量和误差高斯分布假设的约束，适用于任意非线性非高斯动态系统。它还采用蒙特卡罗采样方法来近似状态量的后验概率密度分布，能更好地表现非线性系统的变化信息。粒子滤波算法简单且容易实现，同时与当前主流的卡尔曼滤波系列算法相比，粒子滤波算法中没有复杂的矩阵转置和求逆等运算，因此计算效率更高。相比卡尔曼滤波系列算法直接对粒子的状态值进行更新，粒子滤波算法在对粒子进行更新时，只更新了粒子的权重，粒子实际表示的状态值保持不变，这样能避免在更新过程中出现粒子状态值超出其物理取值范围的情况。

对于粒子滤波的改进，英国教授VanLeeuwen提出的均权重粒子滤波方法，使用提议密度思想有效改善传统粒子滤波中出现的粒子退化和粒子贫化问题，使用较少的集合粒子就可以达到传统方法较多粒子的同化效果。同时基于统计观测的均权重粒子滤波方法有效改善其对于未来观测信息的依赖，使其更好的应用于实时数据同化领域并且有效提升了同化质量。但是均权重粒子滤波方法缺少相应的局地化方案，使该方法很难适应复杂的网格化高纬度模式当中。

发明内容

本发明的目的是为了为均权重粒子滤波提供了了合适的局地化方案，解决了该方法在复杂网格化模式应用的限制，使该方法具有了更好的实际应用的潜力与价值而提供一种基于统计观测局地化均权重粒子滤波数据同化方法。

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

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

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

步骤二：判断是否达到统计观测开始时刻，累加观测求取统计观测均值

根据给定的τ值确定统计观测计算的起始时刻，选择合适的τ可以有效提高统计观测的可靠性，可以更好的引导集合粒子向同化时刻的历史观测靠近。当统计观测开始时，对于同化观测对应位置的历史观测累加求取均值，求解均值可以有效避免对应观测的历史信息出现突然跳变，使用统计均值代替传统方法中的未来观测信息。

步骤三：根据统计观测均值计算提议密度调整集合粒子在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态

根据统计同化时刻观测位置的观测历史均值计算观测对应位置集合粒子的提议密度，根据集合中粒子的提议密度选择其中的最优提议密度，确定集合中每个粒子对于统计观测的后验概率密度，同时调整集合粒子靠近统计历史观测均值，同时根据提议密度更新粒子状态。

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

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

步骤五：使用重采样方法，调整集合粒子维持粒子数稳定，更新观测对应位置粒子状态

使用重采样方法对观测所对应位置的集合粒子进行调整，主要是为了保证对于均权重方法中权重表现较差的集合粒子进行调整，对于权重较小的粒子进行提出，维持集合粒子数的稳定。

步骤六：使用局地化函数，确定同化观测位置周围的粒子权重

在对同化时刻观测对应位置集合粒子进行状态更新后，使用局地化函数继续调整该观测影响半径内的集合粒子状态。对于局地化函数的选择，参考局地化粒子滤波方法中的局地化方案，在计算过程中使用局地化函数描述区域内集合粒子权重与给定观测之间的位置关系，根据已有的同化观测，参考局地化参数确定调整的集合粒子的权重。

步骤七：根据局地化权重更新粒子权重，更新周围粒子状态；

在根据同化观测确定影响的集合粒子的权重后，需要依靠此权重调整对应集合粒子的状态，在状态调整前需要对集合粒子的权重进行归一化处理，根据归一化权重，配合采样粒子与先验粒子的线性组合，对于实现局地化后的集合粒子状态进行更新，使观测附近集合粒子可以根据非对应点观测进行状态调整。

步骤八：计算统计观测局地化均权重粒子滤波的状态后验估计值

在使用局地化方案调整观测周围集合粒子状态后，将全部集合均值重新引入到模式方程中，在同化过程中更新模式积分，得到最终的同化分析结果。

步骤二具体为：判断是否达到统计观测开始时刻，累加观测求取统计观测均值；

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

其中t⁰为前一个需要同化的观测时刻，tⁿ为两个观测之间的观测间隔，t^j表示当前时刻；为了保证统计观测结果的可靠性，当需要同化的观测间隔较大时τ一般选择为0.8-0.9，当同化间隔较小时τ一般选择为0.5-0.8；使用统计的方法对于同化时刻开始前的观测进行统计求取均值，使用

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

步骤三具体为：根据统计观测均值计算提议密度调整集合粒子在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态；

步骤3.1：根据统计观测均值，计算观测对应位置粒子概率密度；

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

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

其中

为提议密度，对于提议密度的求取需要先计算集合粒子的概率密度，使用模式方程，从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可以表示为：

中

和

在这两个表达式中

表示选择的目标权重，

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

步骤六具体为：使用局地化函数，确定同化观测对应位置周围的粒子权重；

在使用均权重粒子滤波方法和重采样方法更新同化观测对应位置的集合粒子状态后，使用局地化函数继续调整该观测影响半径r内的集合粒子状态；局地化函数的选择，参考局地化粒子滤波方法中的局地化方案，在计算过程中使用局地化函数描述区域内集合粒子权重与给定观测位置之间的关系，粒子滤波中权重表示观测的似然概率，在计算中使用局地化算子l[y,x_i,r]，该算子是用来描述观测y和集合粒子状态x_i的相对位置系数，对应的粒子权重表示为：

ω_i,j＝p(y|x_i,j)l[y,x_j,r]+[1-l[y,x_j,r]]

其中局地化算子l[y,x_i,r]主要应用于局部分析，判断集合粒子与观测信息的相对位置信息；当观测y和集合粒子x_i重合时，该函数取最大值为1，当两者之间的距离超过给定的影响半径r则函数值为0，代表该观测对于此集合粒子状态没有调整作用；最后集合粒子的矢量权重表示为：

其中

为观测误差协方差，h_j为测量算子；由此可以看出集合粒子的权重不止于观测信息有关，也与观测和粒子的相对位置有关；

步骤七具体为：根据局地化权重更新粒子权重，更新周围粒子状态；

在获得同化时刻观测信息影响半径内集合粒子的权重后，需要依靠此权重调整对应集合粒子的状态，在状态调整前需要对集合粒子的权重进行归一化处理，确保集合中权重的和为1，归一化权重公式可以表示为：

根据归一化的权重，重新调整集合中粒子状态，对于后验的粒子状态可以表示为：

其中

表示后验均值，k_n是第n个采样粒子，其中向量r₁和r₂可以将新粒子形成采样粒子与先验粒子进行线性组合，最后实现局地化粒子状态的后验更新，其中向量r₁和r₂的计算公式可以表示为：

r_2,j＝c_jr_1,j

公式中

是直到y_i为止的所有观测的误差方差，在这个过程中忽略集合粒子状态之间的后验相关性，但是会通过相应的采样步骤提供对应的相关性，其中对于局地化算子当l[x_j,y_j,r]→1时，c_j→0，并且有

因为当l[x_j,y_j,r]＝1时后验方差近似等于采样粒子方差，同理可得

这就使后验粒子获得采样粒子的状态，当l[x_j,y_j,r]→0时，c_j→∞，

因为当l[x_j,y_j,r]＝0时后验方差等于先验方差，同理可得

该采样方法提供了一种调整粒子以适合观测附近的一般贝叶斯后验解的方法.因为每个采样的粒子都与一个先验粒子组合，所以后验集合包含集合粒子唯一的模型状态，这样就避免了观测值同化期间集合方差的崩溃；

步骤八具体为：计算统计观测局地化均权重粒子滤波的状态后验估计值；

根据观测信息跟新观测对应位置与观测影响半径内所有粒子的状态，计算状态的后验估计集合均值：

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

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

(1)在统计观测均权重粒子滤波中引入局地化方案，在稀疏观测时可以有效提高同化质量，统计观测局地化均权重粒子滤波同化结果均方根误差优于局地化粒子滤波方法和传统均权重粒子滤波的局地化改进；

(2)在引入局地化方案后，统计观测均权重粒子滤波可以更好的应用于复杂的网格化中、高维度模式中，使粒子滤波方法有更好的实际应用潜力。

附图说明

图1是统计观测局地化均权重粒子滤波过程；

图2是统计观测局地化均权重粒子滤波的均方根误差对比图；

图3是传统均权重粒子滤波数据同化流程图；

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

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

具体实施方式

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

为了更简单明确的描述基于统计观测局地化权重粒子滤波同化方法具体实施步骤，以简单的Lorenz-96模式为例进行简单的说明，该模式可以较好的反应局地化方法对于同化结果的影响。

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

选择Lorenz-96模式是因为该模式有40个状态变量，可以更好的应用局地化方案，在Lorenz-96模式中先将模式初始起点输入模型方程积分100万步进行spin-up，使模式变量达到混沌状态，避免模式在同化积分过程中引入模式波动偏差，将达到混沌状态的模式变量作为模式的初始背景场，以此模式背景场为基础，通过模式方程获得粒子滤波集合粒子的初始值，可以有效提高模式方程求解得到的集合粒子状态的可靠性。

步骤二：判断是否达到统计观测开始时刻，累加观测求取统计观测均值

在统计历史观测均权重粒子滤波数据同化中，使用统计观测代替传统方法的未来观测信息，引导集合粒子靠近同化时刻观测的历史统计结果，提高采样过程中有效粒子数，选择合适的阈值τ判断合适开始统计观测。τ的选择与需要同化的观测信息的时间间隔有关，当τ过小时，将会将更多无用观测引入到统计当中，提高统计误差，当τ过大的时候，统计的观测信息会太少，无法正确的计算提议密度引导粒子靠近观测，τ的计算公式为：

其中t⁰为前一个需要同化的观测时刻，tⁿ为两个观测之间的观测间隔，t^j表示当前时刻。为了保证统计观测结果的可靠性，当需要同化的观测间隔较大时τ一般选择为0.8-0.9，当同化间隔较小时τ一般选择为0.5-0.8。使用统计的方法对于同化时刻开始前的观测进行统计求取均值，使用

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

步骤三：根据统计观测均值计算提议密度调整集合粒子在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态

步骤3.1：根据统计观测均值，计算观测对应位置粒子概率密度

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

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

其中

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

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

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

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

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

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

步骤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可以表示为：

中

和

在这两个表达式中

表示选择的目标权重，

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

步骤五：使用重采样方法，调整集合粒子维持粒子数稳定，更新观测对应位置粒子状态

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

步骤六：使用局地化函数，确定同化观测对应位置周围的粒子权重

在使用均权重粒子滤波方法和重采样方法更新同化观测对应位置的集合粒子状态后，使用局地化函数继续调整该观测影响半径r内的集合粒子状态。局地化函数的选择，参考局地化粒子滤波方法中的局地化方案，在计算过程中使用局地化函数描述区域内集合粒子权重与给定观测位置之间的关系，粒子滤波中权重表示观测的似然概率，在计算中使用局地化算子l[y,x_i,r]，该算子是用来描述观测y和集合粒子状态x_i的相对位置系数，对应的粒子权重表示为：

ω_i,j＝p(y|x_i,j)l[y,x_j,r]+[1-l[y,x_j,r]]

其中局地化算子l[y,x_i,r]主要应用于局部分析，判断集合粒子与观测信息的相对位置信息。当观测y和集合粒子x_i重合时，该函数取最大值为1，当两者之间的距离超过给定的影响半径r则函数值为0，代表该观测对于此集合粒子状态没有调整作用。最后集合粒子的矢量权重表示为：

其中

为观测误差协方差，h_j为测量算子。由此可以看出集合粒子的权重不止于观测信息有关，也与观测和粒子的相对位置有关。

步骤七：根据局地化权重更新粒子权重，更新周围粒子状态；

在获得同化时刻观测信息影响半径内集合粒子的权重后，需要依靠此权重调整对应集合粒子的状态，在状态调整前需要对集合粒子的权重进行归一化处理，确保集合中权重的和为1，归一化权重公式可以表示为：

根据归一化的权重，重新调整集合中粒子状态，对于后验的粒子状态可以表示为：

其中

表示后验均值，k_n是第n个采样粒子，其中向量r₁和r₂可以将新粒子形成采样粒子与先验粒子进行线性组合，最后实现局地化粒子状态的后验更新。其中向量r₁和r₂的计算公式可以表示为：

r_2,j＝c_jr_1,j

公式中

是直到y_i为止的所有观测的误差方差，在这个过程中忽略集合粒子状态之间的后验相关性，但是会通过相应的采样步骤提供对应的相关性，其中对于局地化算子当l[x_j,y_j,r]→1时，c_j→0，并且有

因为当l[x_j,y_j,r]＝1时后验方差近似等于采样粒子方差，同理可得

这就使后验粒子获得采样粒子的状态。当l[x_j,y_j,r]→0时，c_j→∞，

因为当l[x_j,y_j,r]＝0时后验方差等于先验方差，同理可得

该采样方法提供了一种调整粒子以适合观测附近的一般贝叶斯后验解的方法.因为每个采样的粒子都与一个先验粒子组合，所以后验集合包含集合粒子唯一的模型状态，这样就避免了观测值同化期间集合方差的崩溃。

步骤八：计算统计观测局地化均权重粒子滤波的状态后验估计值

根据观测信息跟新观测对应位置与观测影响半径内所有粒子的状态，计算状态的后验估计集合均值：

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

本发明提出一种基于统计观测局地化均权重粒子滤波数据同化技术。与传统的统计观测均权重粒子滤波数据同化技术相比，本发明的显著特征在于：针对统计观测均权重粒子滤波方法提出适应其应用的局地化方案，有效解决该方法对于复杂高纬网格化模式应用的限制，在同化时刻，采用统计观测方法更新观测对应位置的集合粒子状态，然后根据同化观测的位置信息和局地化参数，使用该观测调整观测位置周围，对于观测敏感可用集合粒子的权重，再根据权重调节观测附近集合粒子状态，提高观测的利用率，同时提高同化质量。本专利提出的方法可以有效改善传统统计观测均权重粒子滤波方法在复杂网格化模式和稀疏观测条件下的同化能力，可以有效提高同化质量，同时改善统计均权重粒子滤波方法的应用前景。

Claims (1)

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

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

先将模式初始场引入到模式积分方程中，先积分模式方程，使模式方程达到混沌状态，将达到混沌状态的模式变量作为初始背景场，以背景场为积分起点，进行后续的模式状态积分；

步骤二：判断是否达到统计观测开始时刻，累加观测求取统计观测均值；

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

其中t⁰为前一个需要同化的观测时刻，tⁿ为两个观测之间的观测间隔，t^j表示当前时刻；为了保证统计观测结果的可靠性，当需要同化的观测间隔较大时τ一般选择为0.8-0.9，当同化间隔较小时τ一般选择为0.5-0.8；使用统计的方法对于同化时刻开始前的观测进行统计求取均值，使用

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

步骤三：根据统计观测均值计算提议密度调整集合粒子在给定同化时刻，使用均权重方法计算粒子权重，调整粒子状态；

步骤3.1：根据统计观测均值，计算观测对应位置粒子概率密度；

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

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

其中

为提议密度，对于提议密度的求取需要先计算集合粒子的概率密度，使用模式方程，从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可以表示为：

中

和

在这两个表达式中

表示选择的目标权重，

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

步骤五：使用重采样方法，调整集合粒子维持粒子数稳定，更新观测对应位置粒子状态；

步骤六：使用局地化函数，确定同化观测对应位置周围的粒子权重；

在使用均权重粒子滤波方法和重采样方法更新同化观测对应位置的集合粒子状态后，使用局地化函数继续调整该观测影响半径r内的集合粒子状态；局地化函数的选择，参考局地化粒子滤波方法中的局地化方案，在计算过程中使用局地化函数描述区域内集合粒子权重与给定观测位置之间的关系，粒子滤波中权重表示观测的似然概率，在计算中使用局地化算子l[y,x_i,r]，该算子是用来描述观测y和集合粒子状态x_i的相对位置系数，对应的粒子权重表示为：

ω_i,j＝p(y|x_i,j)l[y,x_j,r]+[1-l[y,x_j,r]]

其中局地化算子l[y,x_i,r]主要应用于局部分析，判断集合粒子与观测信息的相对位置信息；当观测y和集合粒子x_i重合时，该函数取最大值为1，当两者之间的距离超过给定的影响半径r则函数值为0，代表该观测对于此集合粒子状态没有调整作用；最后集合粒子的矢量权重表示为：

其中

为观测误差协方差，h_j为测量算子；由此可以看出集合粒子的权重不止于观测信息有关，也与观测和粒子的相对位置有关；

步骤七：根据局地化权重更新粒子权重，更新周围粒子状态；

在获得同化时刻观测信息影响半径内集合粒子的权重后，需要依靠此权重调整对应集合粒子的状态，在状态调整前需要对集合粒子的权重进行归一化处理，确保集合中权重的和为1，归一化权重公式可以表示为：

根据归一化的权重，重新调整集合中粒子状态，对于后验的粒子状态可以表示为：

其中

表示后验均值，k_n是第n个采样粒子，其中向量r₁和r₂可以将新粒子形成采样粒子与先验粒子进行线性组合，最后实现局地化粒子状态的后验更新，其中向量r₁和r₂的计算公式可以表示为：

r_2,j＝c_jr_1,j

公式中

是直到y_i为止的所有观测的误差方差，在这个过程中忽略集合粒子状态之间的后验相关性，但是会通过相应的采样步骤提供对应的相关性，其中对于局地化算子当l[x_j,y_j,r]→1时，c_j→0，并且有

因为当l[x_j,y_j,r]＝1时后验方差近似等于采样粒子方差，同理可得

这就使后验粒子获得采样粒子的状态，当l[x_j,y_j,r]→0时，c_j→∞，

因为当l[x_j,y_j,r]＝0时后验方差等于先验方差，同理可得

该采样方法提供了一种调整粒子以适合观测附近的一般贝叶斯后验解的方法.因为每个采样的粒子都与一个先验粒子组合，所以后验集合包含集合粒子唯一的模型状态，这样就避免了观测值同化期间集合方差的崩溃；

步骤八：计算统计观测局地化均权重粒子滤波的状态后验估计值；

根据观测信息跟新观测对应位置与观测影响半径内所有粒子的状态，计算状态的后验估计集合均值：

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

Images