基于Lorenz-96模型的数据同化方法比较研究

Research on Data Assimilation Methods Based on the Lorenz-96 Model

摘要: 目的：数据同化将观测数据和基于模型的方法有机结合，以实现更精准的预测和状态估计，在多个领域都发挥着积极的作用。近年来，数据同化在核工业领域的重要性不断上升。为了分析现有的同化算法在不同数据场景下的表现，本文对集成卡尔曼滤波(EnKF)、三维变分(3D-Var)和四维变分(4D-Var)数据同化方法进行了详细的对比分析。方法：为了验证上述同化方法在不同数据场景下的同化效果，本文利用Lorenz-96模型生成仿真数据，并添加不同误差水平的噪声，同时获取具有不同观测间隔的数据。利用不同的数据集对同化方法进行分析，通过统计分析时空均方根误差来评估同化效果。结果：实验结果表明：EnKF在不同观测噪声水平和观测间隔的数据下展现出卓越的同化效果，适用于实际复杂系统的数据同化。3D-Var由于仅在当前数据点进行同化，因此在同化速度上较为迅速；而4D-Var方法则对给定窗口内的数据进行同化，导致同化时间相对于3D-Var较长。

Abstract: Introduction: Data assimilation seamlessly combines observational data with model-based methods to achieve more accurate predictions and state estimations, playing a positive role in multiple domains. In recent years, the significance of data assimilation in the nuclear industry has been steadily increasing. In order to assess the effectiveness of assimilation algorithms across various data scenarios, this paper conducts a comprehensive comparative analysis of data assimilation methods, including Ensemble Kalman Filter (EnKF), 3D-Var, and 4D-Var. Method: To validate the assimilation effectiveness of the methods above in various data scenarios, this study uses the Lorenz-96 model to generate simulated data. Different noise levels are added, and data with varying observation intervals are obtained. The assimilation methods are applied to different datasets, and the assimilation performance is evaluated through statistical analysis of the root mean square error in space and time. Result: The experimental results indicate that EnKF demonstrates excellent assimilation performance across different levels of observational noise and observation intervals, making it suitable for data assimilation in practical complex systems. Due to its assimilation being confined to the current data points, 3D-Var exhibits a faster assimilation speed. On the other hand, 4D-Var assimilates data within a given window, resulting in a relatively longer assimilation time than 3D-Var.

1. 引言

作为一种低碳、高效的清洁能源，核能在应对全球气候变化中发挥着积极的作用[1]。随着人们对环境的日益重视，核电的发展将成为未来能源布局的重要组成部分。随着科技的进步，核电行业也迎来了数字化和智能化的转型[2] [3]。在工程科学中，数据同化(Data Assimilation, DA) [4]是一种强大的工具，它将观测数据和基于模型的方法结合起来，更好的进行预测和状态估计。数据同化在天气预测[5] [6]，气象学[7] [8]，海洋学[9]，遥感和水文模型[10]等均有应用。在核工程中的第一次实际应用始于60年代，旨在充分利用当时稀缺的实验数据，发展核反应堆设计概念，并改进以问题为导向的核数据库[11]-[14]。

由于核电厂反应堆运行期间无法直接测量设备的状态，如堆芯反应状态、堆芯燃耗分布等[4] [15]，这些状态量对核电厂的安全性和经济性至关重要。因此，在反应堆工程中广泛应用数值模拟。然而，由于组件制造偏差、流量分布不均等因素，数值模拟与核电厂机组实际状态之间存在偏差。数据同化[4]方法能综合利用数值模拟和实际观测值的信息，实现对核电站真实状态的最有效估计。

近年来，数据同化在核工程中也逐渐发挥积极作用。例如，Garcia等[16]的工作赋予了数据同化在处理更广泛的安全嵌入式智能和集成状态感知的功能。Lu等[3]提到以数据同化为基础的数字孪生侧重于推进核电站作为人–网络–物理系统的概念。Argaud等[4]将数据同化方法应用于核反应堆堆芯，该方法被用来研究如何集成来自多个仪器的信息以及评估仪器损失的影响。同时，通过处理大量数据得到最佳参数估计，并反向改善核反应堆堆芯的仿真模拟精度。Gong等[17]将改进的三维变分(3D-Var)方法在核反应堆的各种场景进行了测试，集成到标准代码包中并应用于实际的核反应堆观测数据的同化。Cheng等[18]提出了一种综合降阶模型和机器学习的方法，创建了一个数字孪生体，在机器学习提取的隐变量上利用变分方法进行数据同化，用于在核反应堆运行阶段预测功率分布。Guo等[15]提出了一种用于PWR燃耗分布的数据同化方法。采用了3D-Var算法进行燃耗分布校准，并采用人工神经网络(ANN)算法建立了燃耗分布与功率分布之间的复杂关系。数值结果表明，该方法使得燃耗分布误差显著减小，功率分布的相对误差最大值可以从9.53%减小到5.11%。

随着计算能力以及计算物理学的进步，例如卡尔曼滤波方法(Kalman)、3D-Var以及四维变分(4D-Var) [15] [19]-[22]等数据同化方法的可靠性逐渐提高。卡尔曼滤波及其变体[19] [22]是最常见的顺序同化方法。集成卡尔曼滤波[19] (Ensemble Kalman Filter, EnKF)常用来估计系统状态，特别是在非线性和高维系统。其基本原理是通过集成多个系统状态的样本，对系统的状态进行估计。EnKF克服了传统卡尔曼滤波在非线性系统中的局限性，提供了一种更适用于实际复杂系统的数据同化方法。例如，Zheng等[23]结合随机游走扩散模型采用EnKF同化核素释放率进行核事故评估。Tang等[24]利用EnKF和高斯多烟团扩散模型对核事故现场实时释放量快速反演。

3D-Var方法可以从最优控制或贝叶斯分析的角度推导出[8]，其基本思想是通过最小化模型模拟结果与观测值之间的差异，来调整模型的初始条件。具体而言，3D-Var [17]方法将观测数据与模型的先验场结合，通过调整先验场，使得模型模拟结果更好地拟合观测数据。这一过程可以通过最小化一个代表观测与模型差异的代价函数来实现。Gong等[17]提出了一种基于逆距离的拟合项，用于3D-Var数据同化。通过对核反应堆模拟中的各种情景进行了一系列数值测试，观察到与其他现有数据同化方法相比，所提出的方法表现出有效性。此外，该方法已经集成到CORCA-3D [25]代码包中，用于同化核反应堆的观测数据。

4D-Var方法[22] [26]是一种高级的数据同化方法，用于改进数值模型的预测能力。其原理是通过最小化一个代价函数，将模型模拟结果与观测在时间窗口内的时刻进行比较。这一时间窗口可以包含过去的观测数据和模型状态，使得4D-Var不仅调整初始条件，还调整了模型在整个时间窗口内的演化过程。通过迭代地调整模型状态和模型参数。代价函数包括了观测与模型的差异以及先验信息的惩罚项。然而，上述的几种同化方法在不同数据场景下的性能比较还有待进一步的研究。

本文主要利用Lorenz-96模型[20]生成仿真数据来对EnKF、3D-Var和4D-Var方法，在观测误差水平和观测间隔因素上进行数据同化性能对比分析。本文的结构如下：在第1节介绍了3个数据同化方法的原理，第2节介绍了Lorenz-96模型及评价数据同化方法的指标，第3节为不同方法的实验，第4节为结论。

2. 数据同化方法

数据同化是将模型预测与观测数据结合的过程；如图1所示，数据同化方法是连接观测数据和数学物理模型的桥梁。数据同化[15] [27]作为一种数据分析技术，将观测系统与状态系统结合起来，通过将观测系统的测量值与状态系统的背景信息结合起来，用于估计无法获取的状态系统真实值。

2.1. EnKF方法

EnKF数据同化方法起源于1969年由Epstein提出的随机动力预测理论，是一种用于数值天气预报、海洋学等领域的数据同化方法，通过整合观测数据来提高数值模型的准确性，目前在各类数据同化领域中得到了长足发展。EnKF的基本思想是构建一个预报状态的集合，用$x_f^{\left(i\right)}$ 表示(其中$i=1,2,\cdots ,N$ , $N$ 是集合的大小)。该集合旨在表示系统的不确定性。集合的均值$\left({\overline{x}}_f\right)$ 被视为最佳估计的总体均值，而集合的样本协方差矩阵被用作最佳估计的预报误差协方差的估计。

来源和权限：数据同化研究团队，https://aics.riken.jp/jp/index.html

Figure 1. Data assimilation “bridges” observation data and models

图1. 数据同化(DA)“连接”观测数据和模型

数学上，预报集合的均值和协方差矩阵分别为：

${\overline{x}}_f=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{x}_f^{\left(i\right)}}$ (1)

$P_f=\frac{1}{N-1}{\displaystyle {\sum }_{i=1}^N\left(x_f^{\left(i\right)}-{\overline{x}}_f\right){\left(x_f^{\left(i\right)}-{\overline{x}}_f\right)}^{\text{T}}}$(2)

其中，${\overline{x}}_f$ 表示预报集合的均值，$P_f$ 是预报误差协方差矩阵，$x_f^{\left(i\right)}$ 表示集合的各个成员。该方法的关键是集合的均值作为总体均值的估计，集合的协方差矩阵提供了集合成员在均值周围的不确定性或分散的估计。然后使用这些估计更新带有观测数据的模型预报，使预报更接近系统真实状态。EnKF特别适用于状态空间为高维且演化是非线性的问题。

2.2. 3D-Var方法

在3D-Var方法中，代价函数被用来描述模拟值和真实值之间的差异，数据同化问题被转化为解决代价函数。代价函数可以表示为方程(3)：

$J\left(x\right)={\left(x-x_b\right)}^{\text{T}}{B}^{-1}\left(x-x_b\right)+{\left(Hx-y_0\right)}^{\text{T}}{R}^{-1}\left(Hx-y_0\right)$ (3)

公式中的数学符号的物理含义如表1所示。

Table 1. The physical meaning of math symbols in cost function of 3D-Var

表1. 3D-Var代价函数中数学符号的物理含义

编号	数学符号	物理含义
1	J(x)	代价函数
2	x	状态向量
3	xb	状态向量的背景值
4	y0	观测向量
5	H	观测算子
6	B	背景误差的协方差矩阵
7	R	观测误差的协方差矩阵

根据协方差矩阵的定义，矩阵被定义为方程(4)：

$B=E\left[\left(x_b-x_t\right)\cdot {\left(x_b-x_t\right)}^{\text{T}}\right]$ (4)

其中，$E[\cdot ]$ 表示数学期望；$x_t$ 是由真实状态构成的向量。然而，实际情况中无法直接测量到真实值，因此无法通过方程(4)直接获得矩阵$B$ 。在本研究中，采用了二阶自回归(the second order autoregressive, SOAR)函数[28]来近似协方差矩阵$B$ ，可以表示为以下方程：

$B=B_{0}C$ (5)

$C_{ij}=\left(1+\frac{r_{ij}}{L}\right){\text{e}}^{-\frac{r_{ij}}{L}}$ (6)

其中，$B_0$ 是一个适当的标量系数，由具体问题确定；$C$ 是相关矩阵；$r_{ij}$ 是径向距离，下标$i$ 和$j$ 表示不同组件；$L$ 是相关长度。协方差矩阵$R$ 由对角矩阵近似，可以通过观测向量求解出。

2.3. 4D-Var方法

3D-Var方法可以被称为静止案例，因为观测、背景和分析都对应于一个固定的时间点。换言之，最小化方程(3)的优化问题仅在状态空间中发生。作为一个扩展，4D-Var数据同化旨在解决空间和时间的优化问题，提供了一个非静止的框架。特别是，模型的动力学被纳入优化问题中，以将不同时间点联系起来。4D-Var的代价函数可以写成如下形式：

$J\left(u\left(t_0\right)\right)={\displaystyle {\sum }_{t_k\in \Gamma }\frac{1}{2}{\left(w\left(t_k\right)-h\left(u\left(t_k\right)\right)\right)}^{\text{T}}{R}^{-1}\left(t_k\right)\left(w\left(t_k\right)-h\left(u\left(t_k\right)\right)\right)}$ (7)

其中，$w\left(t_k\right)$ 是时间$t_k$ 时的测量值，而$\text{Γ}$ 定义了可用观测的时间点集。请注意，这个代价函数的参数是初始条件$u\left(t_0\right)$ 。即4D-Var方法的目的是评估一个初始状态估计，如果在时间上演变，它将产生一条轨迹，该轨迹尽可能接近收集到的测量值(根据干扰噪声的协方差矩阵的逆加权)。

4D-Var通过一阶近似[29]的方式进行求解。从初始条件的先验猜测开始，通过在时间上向前求解模型，计算出基础轨迹，直到对应于最后一个观测点的最终时间(即，$k=0,1,2,\cdots ,ON$ )。然后，在最终时间评估特征值。之后，一直沿着时间点迭代，直到结果收敛。

3. 不同数据同化方法的敏感性实验

3.1. 实验仿真模型及数据

Lorenz-96模型是1996年Lorenz和Saltzman在研究流体有限振幅对流时提出的非线性模型[30]。作为一种强非线性模型，经常被用来验证数据同化系统的性能。Lorenz-96模型[20]有40个状态变量，可以在研究中进行不同的组合和取舍，其表达式如方程(7)所示。其中，右侧第一项是平流项，而第二项表示阻尼。此外，$F$ 代表外部迫使常数，设定为8，这是一个已知会导致混沌行为的整体值[31]。Lorenz-96是一个周期系统，常用于模拟大气随时间演变。

$\frac{\text{d}{x}_i}{\text{d}t}=\left(x_{i+1}-x_{i-2}\right)x_{i-1}+F$ (7)

其中，$i=1,2,\cdots ,N$ 是循环标号，在本文中取$N=40$ 。如图2所示，Lorenz-96模型生成了40个状态变量的数据，图中展示了数据的热力图。横轴表示状态变量$x_i$ ，纵轴表示时间$t$ 。模型使用四阶Runge-Kutta方案进行积分，时间步长为$h=0.025$ 。真实模型的解算子被表示为${ℳ}_{L96}$ 。该模型具有上述给定的$N$ 和$F$ 值，并处于混沌状态。

Figure 2. Lorenz-96 model simulation data

图2. Lorenz-96模型仿真数据

3.2. 评价指标

鉴于真实模型${ℳ}_{L96}$ 的解算子和真实场$x_{1:K}$ ，我们使用分析结果的时空均方根误差(${\text{RMSE}}_a$ , Root Mean Square Error - analysis)来评估数据同化过程结束时的同化准确性，${\text{RMSE}}_a$ 计算公式如下：

${\text{RMSE}}_a=\sqrt{\frac{1}{m\left(K-k_0\right)}{\displaystyle {\sum }_{k=k_0}^K{\displaystyle {\sum }_{n=0}^{m-1}{\left(x_{n,k}-x_{m,k}^a\right)}^2}}}$ (8)

其中，$x_k=\left[x_{1,k},x_{2,k},\cdots ,x_{m-1,k}\right]$ ，而$k_0>0$ 是评估中考虑的第一个值。本文中，我们选择$k_0=100$ 以消除初始状态的影响。

4. 不同数据同化方法的敏感性实验

为了研究EnKF、3D-Var和4D_Var数据同化方法的同化效果，选择观测数据噪声水平和观测频率等作为关键因子，研究关键因子对数据同化方法的同化效果的影响。如图3所示，三个同化方法对Lorenz-96模型产生的数据进行同化，图中展示的是仿真数据中的第1、8维状态参数被三个同化方法同化的结果。其中，黑色曲线表示Lorenz-96模型真实的状态值，红色点表示观测数据点，该数据是由真实值加上噪声得到的。蓝色、橙色和绿色的曲线分别表示3D-Var、4D-Var和EnKF方法的同化结果。

4.1. 观测误差对同化效果的影响

本实验中，为了研究观测数据质量对同化效果的影响，观测误差的方差取值从0.6到1.9变化。其中，EnKF的集合数为24，4D-Var方法中的窗口系数为2，该值与观测频率相乘得到4D-Var方法的同化窗口长度。观测误差的变化对三个数据同化方法的影响如图4所示，蓝色、橙色和绿色的线分别表示3D-Var、4D-Var和EnKF方法的同化效果趋势。曲线颜色表示的同化方法在下面的图中保持一致。

Figure 3. Assimilation effects of three methods on Lorenz-96 data

图3. 三种同化方法对数据的同化效果

Figure 4. The impact of observation errors on assimilation effectiveness

图4. 观测误差对同化效果的影响

从图4可以看出三个同化方法具有相似的变化趋势，随着观测误差的变大同化结果与真实值之间的${\text{RMSE}}_a$ 也随之变大。三个同化方法中，EnKF的同化效果最好，3D-Var较4D-Var方法没有考虑到时间信息，因此效果稍差。

4.2. 观测频率对同化效果的影响

在本实验中，观测误差值固定为1，观测的频率从1到7变化。其中，观测频率表示对仿真数据点的观测间隔，如观测频率为2时，即表示对仿真数据隔一个点进行观测。三个同化方法的同化效果随观测频率的变化如图5所示。

Figure 5. The impact of observation frequency on assimilation effectiveness

图5. 观测频率对同化效果的影响

从图5可以看出，随着观测频率的增大，三个同化方法的${\text{RMSE}}_a$ 值均有轻微的上升趋势。其中EnKF方法在三个同化方法中${\text{RMSE}}_a$ 最小，同化效果最好。但随着观测频率的增加，EnKF方法的${\text{RMSE}}_a$ 的增长较其他两个方法更快，说明EnKF方法容易受到观测频率的影响。

4.3. 同化方法的时间效率对比

本实验中，我们探究了数据的观测频率、观测误差对同化方法的计算时间的影响。如图6所示，图中展示了三个同化方法的同化时间随着观测频率的变化情况，纵轴表示数据同化时间，单位为秒。

Figuure 6. The impact of observation frequency on assimilation time

图6. 观测频率对同化时间的影响

从图6中可以看出，3D-Var和EnKF方法随着观测频率的增加，同化时间呈现出递减的趋势，而4D-Var方法的同化时间较前两者高出很多，且随着观测频率的增加，同化时间出现上升的趋势。这可能是4D-Var方法在同化过程中考虑了时间信息，同化过程在一个同化窗口中进行的，随着观测频率的增加同化窗口也随着变长，使得同化的运算量增加。而3D-Var和EnKF方法只在当前观测点和当前观测点前后进行同化，随着观测频率的增加同化的次数变少使得同化时间减少。

Figure 7. The impact of observation errors on assimilation time

图7. 观测误差对同化时间的影响

观测误差对同化时间的影响如图7所示，从图中可以看出，三个同化方法在不同观测误差数据下的同化时间基本上没有太大的变化，说明观测误差对同化时间的影响很有限，更多的是影响同化的精度。其中，4D-Var方法的同化时间还出现了小的波动，这说明4D-Var方法在同化过程中，寻找最优解中受到观测噪声的影响而表现出的不稳定。

5. 结论

通过Lorenz-96模型进行数据仿真，添加不同误差水平的噪声生成了多个仿真数据集。同时，获取了具有不同观测间隔的观测数据。利用这些数据集对EnKF、3D-Var和4D-Var数据同化方法在观测数据误差和观测间隔上的同化效果进行对比分析。实验结果显示：

1) EnKF在不同观测噪声水平和观测间隔的数据下同化的综合效果优秀，适用于实际复杂系统的数据同化。

2) 3D-Var由于只考虑当前数据点的同化，在同化速度上较快。

3) 观测数据噪声水平对三个同化方法同化时间影响不大，主要影响同化的精度。

NOTES

^*通讯作者。