一种面向参数和状态变量时空异步更新的数据同化方法

技术领域

本发明涉及地理信息科学技术领域，尤其涉及一种面向参数和状态变量时空异步更新的数据同化方法。

背景技术

模型与观测是现代科学中两种基本的研究手段，它们共同促进了人类对于物理世界的认知。但同时这两者都存在不确定性，导致人类认知存在偏差。数据同化是融合模型与观测、减少认知不确定性、提高物理世界模拟与预报研究精度的方法论。

数据同化是在模型的动力框架内，融合不同来源和不同分辨率(如时间、空间、光谱等)的直接与间接观测，从而增强系统的可预报性和可观测性。数据同化起源于20世纪50年代气象预报中对于客观分析的需求。在数学上主要借助于现代估计理论、控制论和优化方法。现代数据同化方法出现于1990年代之后，可被分为连续数据同化和顺序数据同化两大类。其中。顺序数据同化方法是指在系统运行过程中，当有观测的时刻，利用观测信息在观测和模型误差分别加权的基础上对模型状态进行更新，从而获得模型状态的后验优化估计；状态更新后，模型利用新的状态重新初始化，继续向前积分，直到获得新的观测信息。

参数作为数学物理模型中的重要元素，其作用是控制和调节模型的行为与轨迹，使其更接近真实世界的动力过程。合理准确的参数可以提高模型的稳定性、加快收敛速度、改进精度，从而帮助人类更好地理解科学问题的本质和特征。因此，如何在动力系统的模拟过程中对参数进行修改和优化，不仅是模型校准和适应性研究的主要内容，也是数据同化研究的重要方向。

然而，现有的数据同化参数估计研究，或者将参数视为模型状态扩展向量中的附属要素，在状态更新的同时对参数进行更新，弱化了外来信息对参数的影响，并且强制使参数与状态的更新周期一致，这与真实场景不符(因为参数的时空异质性往往弱于状态变量，其变异周期更长)，给同化系统带来更多不确定性；或者没有对状态变量进行估计，无法充分发挥数据同化融合观测信息以改进模型预报能力的关键作用。

发明内容

本发明所要解决的技术问题是提供一种可以实现与真实场景一致的参数估计和状态变量实时更新的面向参数和状态变量时空异步更新的数据同化方法。

为解决上述问题，本发明所述的一种面向参数和状态变量时空异步更新的数据同化方法，包括以下步骤：

(1)利用不同时刻、不同空间及不同物理含义的观测信息，基于贝叶斯最大后验概率估计理论，确定模型参数与状态变量及其时空变异特征，确定其估计的时间周期和空间距离，建立外部观测与参数、状态的关系；

(2)基于参数估计算法，对模型的参数进行估计；

(3)基于非线性滤波算法，优化模型状态变量；

(4)依据模型的时空变异特征，将模型的运行周期划分为参数估计和数据同化两个窗口，并将参数估计算法和改进状态更新集成到数据同化方法框架中，采取不同的计算流程更新参数和状态变量；

(5)运行同化系统，依据估计周期在不同的时域周期上分别对状态和参数进行最优估计，并在系统运行节点开放接口以实现参数对模型结构、状态对模型运行轨迹的调整。

所述步骤(2)中参数估计算法是指改进后的马尔可夫链-蒙特卡洛采样方法，该方法基于一种典型的马尔可夫链-蒙特卡洛采样即Metropolis-Hastings采样，在贝叶斯推断的概念下，对参数θ进行估计，等价于在给定观测集Ω的情况下，求后验概率P(θ|Ω)；P(θ|Ω)∝P(Ω|θ)P(θ)；其中：P(θ)表示参数先验分布，P(Ω|θ)为似然函数，表示参数为θ时观测集为Ω的概率。

所述步骤(3)中非线性滤波算法是指改进后的集合卡尔曼滤波算法；

系统状态分析值为

其中：t代表时间；x

本发明与现有技术相比具有以下优点：

1、本发明聚焦于发展一种更合理、高效的数据同化方法，用于动力系统中参数最优估计和状态变量更新，其估计和优化思路既参考外来观测信息，也分别依据参数和状态的时空变异特征，从而构建一种呈现双循、交叉式的参数估计和状态更新策略，可以有效地估计与真实场景一致的参数，并实现状态变量的实时更新，从而面向非线性动力系统提供一种预报能力更强的数据同化方法。

2、本发明构建了一种基于贝叶斯最大后验概率分布估计的方法，依据外部观测信息和时空异质性特征，合理划分参数估计和状态更新的区间，并采用一系列措施(如方差定时更新、动态调整参数估计区间、数据的预处理等)显著提高算法运行效率，从而提供更合理和高效的数据同化应用方案。

3、本发明修正了传统数据同化方法不能面向真实场景为模型参数和状态进行差异化估计的问题，将改进后的马尔可夫链-蒙特卡洛采样方法集成到数据同化算法中，实现了更合理和有效的参数估计和状态更新策略。

4、本发明适用于在非线性动力系统中，融合观测信息以改进随时空变化的参数及状态变量的估计研究。

附图说明

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

图1为本发明在非线性动力系统中实施面向参数和状态变量时空异步更新的数据同化方法流程图。

图2为基于改进后的马尔可夫链-蒙特卡洛采样实施参数估计的流程图。

图3为基于改进后的集合卡尔曼滤波实施状态变量优化的流程图。

图4为基于本发明方法的实施例1：基于接触网络模型的状态和参数异步更新结果图。

图5为基于本发明方法的实施例2：基于传染病模型的状态和参数异步更新结果图。

具体实施方式

本发明是一种在非线性动力系统中，基于贝叶斯最大后验估计理论，发展了非线性滤波结合马尔可夫链-蒙特卡洛采样的算法，在模型运行的同时，参考时空变异特征与观测信息，分别对参数和状态变量进行优化估计，进而提高数据同化系统的预报精度。

如图1所示，一种面向参数和状态变量时空异步更新的数据同化方法，包括以下步骤：

(1)利用不同时刻、不同空间及不同物理含义的观测信息，基于贝叶斯最大后验概率估计理论，确定模型参数与状态变量及其时空变异特征，确定其估计的时间周期和空间距离，建立外部观测与参数、状态的关系，以获得更符合真实场景的同化系统预报结果。

具体过程：设定模型参数θ与状态变量的初始值，设定参数估计的时间周期和空间距离，根据先验知识设定每个参数的估计范围。根据与参数相关的外部观测数据Ω，设定参数估计的时空窗口。

(2)启动模型，当运行至参数估计的窗口时，基于参数估计算法，采用如图2所示策略对模型的参数进行估计。

参数估计算法是指改进后的马尔可夫链-蒙特卡洛采样方法，该方法基于一种典型的马尔可夫链-蒙特卡洛采样即Metropolis-Hastings采样，在贝叶斯推断的概念下，对参数θ进行估计，等价于在给定观测集Ω的情况下，求后验概率P(θ|Ω)。这里P(θ|Ω)∝P(Ω|θ)P(θ)；其中：P(θ)表示参数先验分布，P(Ω|θ)为似然函数，表示参数为θ时观测集为Ω的概率。P(Ω|θ)计算需要借助马尔可夫链-蒙特卡洛采样，根据提议分布q(·)生成近似后验样本。

具体计算步骤为：

①构建n维参数空间{[a

②迭代步k＝0，参数空间中随机选择一个n维参数向量θ

③迭代步k＝k+1，根据先验的提议分布q(·)对n维参数向量θ

④随机选择接受率ω∈(0，1)；

⑤根据似然函数P(Ω|θ

i将参数θ

ii计算似然函数

iii计算

iV如果ω≤α，接受采样θ

⑥k＝k+1，判断迭代步数是否达到预先设定的最高值，如是则停止迭代，否则跳到步骤②；

⑦基于所有接受的采样集合及其分布特征，计算均值、中位数或众数，得到参数的最优估计值。

在应用时，设定提议分布分别为均匀分布和正态分布，对应的最高迭代步数分别为s和r，运行步骤①～⑥，其中正态分布的协方差矩阵由均匀分布的采样结果给出，其对角线值为每个参数中所有被接受的采样的估计方差。

(3)基于非线性滤波算法，优化模型状态变量。

当模型运行至数据同化的窗口时，采用如图3所示策略实现参数估计。图中：t代表在时刻，θ代表参数，x

其中：非线性滤波算法是指基于一个典型的动力方程M和观测方程H系统的集合卡尔曼滤波算法；

一个典型的动力方程M和观测方程H系统为

其中：t代表时间；x

因此，基于该系统的集合卡尔曼滤波算法为：

①t＝0，扰动初始状态场x

②利用动力方程计算t+1时刻的状态预报值x

③如果t+1时刻有可用观测y，则计算状态分析值

其中：K＝P

没有观测则

④计算t+1时刻的分析误差

⑤令

(4)依据模型的时空变异特征，将模型的运行周期划分为参数估计和数据同化两个窗口，并将参数估计算法和改进状态更新集成到数据同化方法框架中，采取不同的计算流程更新参数和状态变量。

(5)运行同化系统，依据估计周期在不同的时域周期上分别对状态和参数进行最优估计，并在系统运行节点开放接口以实现参数对模型结构、状态对模型运行轨迹的调整。

本发明所采用的改进后的马尔可夫链-蒙特卡洛采样方法是一种有记忆的计算方法，仅在设定的时刻后再更新正态分布的协方差矩阵，可以避免迭代过程中反复运行动力模型引起的运算耗时问题，有效提高了提议分布的方差计算效率。同时根据采样的统计分布特征决定使用均值、中位数还是众数作为最优参数，避免了因低效采样(如有效采样量过低，采样代表性不足等)而导致的参数误估问题。另外，还可以根据实时采样结果动态调整参数估计区间，可有效解决在高维空间采样效率不高的问题。

本发明所采用的非线性滤波算法是在标准的集合卡尔曼滤波方法流程上设计了动态的误差赋值方法，如令

实施例1使用接触网络模型为模型算子(50万个网络节点)，集合卡尔曼滤波为同化算法(集合成员为20个，集合扰动比率为10％)，Metropolis-Hastings采样为参数估计算法(均匀分布采样迭代次数s＝3000；正态分布采样迭代次数r＝1500)，观测误差和模型误差比率为1：2。结果如图4所示。

由图4可以发现，采用本发明方法，其结果要明显优于只使用参数估计或直接模拟的结果(利用最后3个观测验证，平均精度提高了86.19％)。

实施例2使用传染病模型为模型算子，集合卡尔曼滤波为同化算法(集合成员为100个，集合扰动比率为40％)，Metropolis-Hastings采样为参数估计算法(均匀分布采样迭代次数s＝30000；正态分布采样迭代次数r＝5000)，观测误差和模型误差比率为1：10。随机选择30个时间点做100次实验，结果如图5所示。

由图5可以发现，采用本发明方法，进行预报，与仅使用参数估计的预报相比，决定系数从0.62上升到0.98，均方根误差RMSE从221.29降低到30.57，因此显著提升了预报精度。