基于Lorenz-96 的数据同化基础算法数值实验系统设计与实现
2020-12-16王际朝臧绍东阮宗利吴国丽
王际朝，臧绍东，阮宗利，吴国丽
（中国石油大学（华东） 理学院，山东 青岛 266580）
数据同化（data assimilation）亦称资料同化，根据一定的标准和方法，将不同时空、不同分辨率、不同手段获取的观测数据与动力模式有机结合，建立数据与模型相互协调的优化关系，从而给出当前模式状态的变量分析场。作为一种技术，数据同化已广泛应用于大气[1]、海洋[2]、陆面[3]、石油工程[4]等诸多领域。尤其是当前各领域的观测数据飞速积累，数值模拟技术也日趋成熟，数据同化作为连接观测和模式的桥梁，其作用愈发重要。
国内很多高校，尤其是涉及大气、海洋、石油等专业的院校，面向高年级本科生或研究生开设数据同化课程。数据同化课程的教学内容是一系列同化方法或者算法的理论讲解和上机实验。此课程的理论基础是估计理论和控制理论[5]，教师在教学过程中易偏重理论讲解，大量的复杂公式使学生感到抽象而缺乏兴趣，很难发挥学生学习的积极性和主动性。算法是数据同化课程的本质，在完成算法详细讲解的同时，教师应适当利用计算机可视化技术，使得教学内容更形象和生动。
本课题组查阅了大量资料，发现鲜有对此课程进行教学研究的文献，这可能是由于数据同化课程的“科研” 属性较强。本文基于Lorenz-96 模型，设计并实现了数据同化课程基础算法数值实验系统，旨在为数据同化教学过程中基础算法的讲解提供一个可视化的载体，从而加深学生对算法的理解程度，培养学生的动手能力，提高学生的学习兴趣[6]，进而提升教学质量。
1 系统设计
考虑到数据同化各种方法的理论基础、出现时间、是否业务化应用主流、是否当前研究前沿等，本文将数据同化课程中涉及的基础算法分为5 类：早期同化方法、变分同化算法、卡尔曼滤波算法、粒子滤波算法和混合方法。所有算法都在统一的数学模型下进行程序设计，以方便比较各种方法的同化效果。
1.1 算法模块
早期同化方法是指出现时间较早、实现简单的方法。变分同化算法和卡尔曼滤波算法是当前业务化应用的主流方法。粒子滤波算法由于不需要高斯分布假设，近些年越来越受到关注。混合方法是把前述的算法进行混合，充分利用不同算法的优点，也是当前同化方法的一种趋势[7]。由于方法名字本身的中文翻译没有统一的标准，本文同时给出各同化方法的英文。
早期同化方法包括：多项式插值（polynomial interpolation）、逐步订正法（successive corrections method）、松弛逼近法（nudging method）和最优插值（optimal interpolation，OI）法。
变分算法包括： 三维变分（ three-dimensional variational，3D-VAR）和四维变分（four-dimensional variational，4D-VAR）算法。
卡尔曼滤波算法包括：卡尔曼滤波（Kalman filter，KF）、扩展卡尔曼滤波（extended Kalman filter，EKF）、集合卡尔曼滤波（ensemble Kalman filter，EnKF）、集合调整卡尔曼滤波（ensemble adjustment Kalman filter，EAKF）、集合转换卡尔曼滤波（ensemble transform Kalman filter，ETKF）和局地化集合转换卡尔曼滤波（local ensemble transform Kalman filter，LETKF）算法。
粒子滤波算法包括：基本粒子滤波（particle filter，PF）、重采样粒子滤波（resampling particle filter，RPF）、隐式粒子滤波（implicit particle filter，IPF）、等权重粒子滤波（equivalent-weights particle filter，EWPF）和映射粒子滤波（mapping particle filter，MPF）算法。
混 合 方 法 包 括 ： 4DVAR+EnKF 、 PF+EnKF 、4DVAR+KF 和EnVarPF 方法。
该实验系统的主要功能模块如图1 所示。
图1 数据同化基础算法实验系统功能模块
数据同化的算法仍在不断改进和发展中，部分新的算法或者非常复杂的算法并未包括在本实验系统中。本文设计的模块基本涵盖了数据同化课程的基础算法，同时兼括了部分近几年发展起来的同化算法。
1.2 数学模型
数据同化算法的设计与实现需要数学模型作为基础。经过多年数据同化课程教学经验的积累，本文选择Lorenz-96 模型作为所有算法的模型基础，即本文的所有算法均基于此数学模型进行编程和设计。Lorenz-96 模型是高度非线性的，且对初值极端敏感，是验证各种同化算法的理想模型[8]。此模型是Lorenz等于1996 年从动力模式中推导出的简化的动力模型，被广泛地用来验证各类同化方法。其动力学方程为
其中：F为外强迫项，n为变量个数。Lorenz-96 系统的线性耗散项使得总能量迫项F能增大或降低总能量，二次平流项对贡献，具有保存总能量的性质。
2 系统实现
Matlab 软件为数据同化课程的实验教学提供了强有力的设计平台[9]。 GUI 是指图形用户界面（graphical user interface，GUI），由窗口、按钮、文字、图形、菜单等对象构成的用户视窗，是实现人机交互的有效工具[10-11]。本文利用Matlab GUI，设计了数据同化基础算法数值实验系统。
2.1 实验系统界面
Matlab GUI 形成包含GUI 组件属性的扩展名为“.fig” 的文件和保存程序代码的扩展名为“.m” 的文件[12]。GUI 程序主要包括各控件的初始化函数、输出函数、布局以及回调函数。实验系统的主界面采用经典的菜单模式，主菜单上的5 个部分对应5 个功能模块，点击可以弹出各子模块对应的二级子菜单。单击每个子菜单会调出各子模块对应的界面。实验系统主界面如图2 所示。
图2 数据同化基础算法数值实验系统主界面
点击主菜单的下拉二级菜单，会调出各个算法的执行模块。由于算法众多，本文优选出集合卡尔曼滤波算法进行实验系统的功能解析，并给出算法的原理和主要程序的说明。
2.2 实验系统功能解析
集合卡尔曼滤波（EnKF）是由Geir Evensen 利用蒙特卡洛思想提出的，随后涌现了基于此方法思想的各种集合卡尔曼滤波类同化方法。集合卡尔曼滤波在大量的理论和实际业务化数据同化问题上已经被证明是非常有效的。
2.2.1 算法流程
EnKF 算法中，计算预测误差协方差矩阵是通过集合成员完成的，因此能够巧妙地避开切线性算子的求解问题，对非线性系统有着更好的适用性。这种处理方式不仅解决了KF 方法应对非线性系统的局限，同时保留了误差协方差矩阵随时间演化的特点。由于EnKF 方法中对于除了卡尔曼增益矩阵的计算外其他所有在集合成员之间的操作都是互相独立的，因此很容易运用计算机进行并行计算，这是EnKF 方法得到推广的重要原因。EnKF 算法流程如下：
（4）分析步和预测步循环进行。
2.2.2 程序实现
在EnKF 算法中，数据同化分析步是整个算法的核心。同化分析步作为一个函数 “EnKF_analysis”，供实验系统调用。用户只需输入预测状态矩阵、观测向量、观测算子、观测误差协方差矩阵和集合数，即可得到同化分析值等。其Matlab 程序如下：
Function [xa,mean_xf,mean_xa] =
EnKF_analysis(xf0, y, H, R,n)
%% EnKF 的分析步
%% 输入:
%% xf0 -预报状态矩阵
%% y - 观测向量
%% H - 观测算子
%% R - 观测误差协方差矩阵
% %n -集合数
%% 输出:
% xa -分析的集合.
% mean_xf - 预报集合的均值.
% mean_xa - 分析集合的均值.
%%----------------------------------------------%%
%预报集合均值
mean_xf = mean(xf0,2);
% 预报误差协方差矩阵
Pf = ((xf0 - repmat(mean_xf,1,n)) * (xf0 -repmat(mean_xf,1,n))') / (n-1);
% 观测数
p = size(y,1);
% 观测集合随机误差
u = 0 + sqrtm(R) * randn(p,n);
% 观测集合
y_ensemble = repmat(y,1,n) + u;
% 观测误差协方差矩阵
Ru = (u * u') / (n-1);
% 集合卡尔曼增益
Ku = Pf * H' * inv((H * Pf * H' + Ru));
% 分析集合
xa = xf0 + Ku * (y_ensemble - H * xf0);
% 分析均值
mean_xa = mean(xa,2);
end
2.2.3 功能演示
在主界面下点击菜单“卡尔曼滤波算法”，在弹出的下拉菜单中点击 “Ensemble Kalman Filter (EnKF)”（如图3 所示），实验系统弹出EnKF 算法子模块的界面。
图3 EnKF 二级菜单
在弹出的EnKF 算法界面中，输入变量数、集合数、同化间隔、观测算子等初始化参数后，点击“Run_EnKF” 按钮，如图4 所示，实验系统自动给出同化效果的结果图形。效果图展示了EnKF 同化观测数据后与真值比较的结果，其中黑色虚线表示 “真值”，红色实线表示EnKF 的分析值，蓝色方框表示输入参数规定的观测数据。
图4 第1 个变量的同化效果
图4 展示的是第1 个变量的同化效果。在不改变其他输入参数的情况下，实验系统亦可展示其他变量的同化效果。图5 给出的是第19 个变量的同化效果。从效果图4 和5 中可以看出，EnKF 算法具有较好的同化效果。根据计算结果，亦可分析均方根误差、相关系数等误差统计指标值。
图5 第19 个变量的同化效果
此EnKF 子模块可对输入参数进行任意符合规则的设置，提高了算法的人机交互性能。由于EnKF 算法输入参数较多，有些参数是作为子程序进行输入的，比如观测算子H的构造，其本身就是一个子程序，输入和输出要符合实验系统的要求。这里给出“参数”Generate_H 的Matlab 程序。实验系统的用户可在后台修改这些参数，以满足个性化的需求。
function [H] = Generate_H(p,n)%%产生观测算子矩阵
H = zeros(p,n); %预分配空间
if rem(n,p) == 0; %如果能整除
j = 1:n/p: n;
for i = 1: p
H(i,j(i)) = 1;
end else
sd=1; %随机种子
rng(sd);
k_randn = randperm(n,p);
[k_randn] = sort(k_randn,'ascend');
j = k_randn; for i = 1 : p
H(i,j(i)) = 1;
end
end
end
除了可以展示前述功能外，EnKF 子模块也可修改各输入参数以查看算法的同化效果，分析各参数的敏感性。甚至，用户可以直接修改Maltab 的“.m” 文件来实现定制的个性化算法，以满足教学和科研的需要。本文仅给出了EnKF 子模块的功能展示，根据算法原理的不同，其他子模块的界面、输入参数、输出结果的展示方式等略有不同。培养学生的能力是大学教学的核心[13]，本系统的所有参数、代码均可修改，可培养学生的自主学习和动手能力，为数据同化课程的研究性实验教学提供了很好的实验平台。
3 结语
基于Lorenz-96 数学模型，利用Matlab GUI，设计和实现了数据同化课程基础算法数值实验系统，并以EnKF 算法为例展示了实验系统的部分功能。系统界面友好，交互简单，内容丰富，流程清晰，演示和扩展方便。师生输入算法所需参数即可获得可视化实验结果，并可对比分析不同条件下、不同算法的数据同化效果，使得原来算法的枯燥理论推导过程可在“所见即所得” 的环境下变得更加生动和形象。本实验系统还应用到了线上教学活动中，丰富了教学手段，取得了理想的教学效果。
数据同化课程基础算法数值实验系统为课程理论教学和实验教学的有机结合提供了抓手，为实验教学提供了强有力的载体，丰富了教学手段，可用在数据同化课程的的辅助性教学中以提升教学质量。本系统亦可为学生提供一个同化算法的数值实验平台，有益于加深学生对算法的理解和提高学生的动手能力。
