CN111783361B - 基于三重多层感知机的数值天气预报混合资料同化方法 - Google Patents

Abstract

本发明公开了基于三重多层感知机的数值天气预报混合资料同化方法，包括以下步骤：建立基于三维变分资料同化方法的第一感知机模型，并对第一感知机模型进行训练；建立基于集合卡尔曼滤波器资料同化方法的第二感知机模型，并对第二感知机模型进行训练；建立混合同化的第三感知机模型；对第三感知机模型进行训练；使用训练好的第一感知机模型、第二感知机模型和第三感知机模型来计算获得分析场数据。本发明方法利用大气运动的物理规律和大气状态变量的时间特性，对传统的资料同化方法进行模拟、优化和混合，充分考虑了大气变量的流依赖，以此引入大气数据本身所包含的时间特征，使得本发明方法的同化性能较传统方法得到有效提升。

Description

基于三重多层感知机的数值天气预报混合资料同化方法

技术领域

本发明属于数值天气预报技术领域，具体涉及基于三重多层感知机的数值天气预报混合资料同化方法。

背景技术

Vihelm Bjerknes(1904)把数值天气预报(NWP)定义为一个初边值问题。即大气的未来状态取决于其详细的初始状态，边界条件以及预测模型。资料同化(DA)可以被描述为一个过程，将观测和短期预报(也就是背景场或初猜测值)中的所有可用的数据资源结合起来，得到分析场，希望得到对实际大气状态的概率密度的最佳估计。即资料同化的目标是为数值天气预报产生最精确的初始条件。

在过去的三十年里，基于贝叶斯定理的两种主流资料同化方法蓬勃发展起来。一种是变分同化(Var)，另一种是基于集合卡尔曼滤波(EnKF)的集合同化。基本上，为了处理来自不同信息的不确定性，这些方法以相同的策略为指导，使用新的观测值来修正指定网格点的初猜测值，即通过最小化代价函数来组合背景场和观测值，代价函数由两个方面组成：一部分是对背景均值距离的惩罚项和另一个惩罚项是对观察值的距离。但是，具体的实现方式各不相同。首先，Var通过给出一个静态背景误差协方差矩阵(B)直接同化大量的观测数据，该矩阵是根据气候数据事先计算得到的。因此，静态B不能正确地表示真正的时空误差统计。Var分为3D Var和4D Var，前者利用单时刻数据，后者利用给定长度的时间窗口的数据估计初始状态。虽然在4D Var中，B隐式地随流依赖变化而演化，并且可以获得比3DVar更好的初始状态，但它需要昂贵的计算代价，并且需要建立线性模型和伴随性数值模型，这是很难建立的。其次，EnKF提供了集合驱动的B，随着系统的发展而变化，优于3D Var的B。此外，EnKF不需要线性模型和伴随性数值模型，这节省了巨大的计算开销。

尽管变分和集合同化方法在业务应用上取得了良好的效果，但是两者存在自身缺陷。为了最大限度地发挥这两种方法的优势，消除它们的不足，出现了一系列的技术来耦合两种方法，称为混合资料同化方法。集合-变分方法以背景误差协方差的流依赖性为出发点，结合变分方法的静态背景误差协方差和集合方法的动态背景误差协方差，重构背景误差协方差。实验证明，将两个背景误差协方差联合起来得到一个重新校准且更稳健的B是有效的。也就是说，混合资料同化方法与纯变分方法或纯集合方法更具有竞争性。然而，集合预报误差协方差和静态背景误差协方差不同的贡献比率直接影响分析场的质量。尽管NWP中心有能力在业务系统中实现混合资料同化方法，但在资料同化系统中要处理不确定性和自由度大的背景误差协方差矩阵(通常≥O(10⁷))需要很大的代价。因此，实际业务上使用经验比例来减轻巨大的计算成本，但牺牲了最佳的精度。

机器学习是一种数据驱动的方法，它在一定程度上有助于解决在地球系统科学中与时空相关的问题。其中，“物理”子模型的公式是指基于物理基础，采用经验参数，存在半经验性质，则这个子模型可以被机器学习模型所替代。自然地，将非线性神经网络与动力学模型相结合，形成混合神经-动力学模型，成为一类新的变分资料同化方法。在近似相同的资料同化效果下，神经网络在计算速度上有显著提升，多次实验检验，速度提升接近90倍，甚至更高。

然而，资料同化的最终目标是快速准确地获得最优分析场，为数值天气预报模式提供最优初始场。前述所提到的已拓展的神经网络资料同化方法，只是利用大量的气象资料，通过模拟已有的资料同化方法来加速同化过程，并没有提高同化效果。其主要原因是，研究人员仅仅关注了数据的统计特征，而忽略了大气状态本身的物理性质。

发明内容

有鉴于此，本发明目的在于提供基于三重多层感知机的数值天气预报混合资料同化方法，本发明方法利用大气运动的物理规律和大气状态变量的时间特性，对多种传统的资料同化方法进行模拟和优化，自适应地学习耦合权值矩阵，将多种资料同化方法进行耦合，使得本发明方法提出的基于神经网络的混合同化方法的同化性能较传统方法得到有效提升。

为达到该目的，本发明采用如下技术方案，基于三重多层感知机的数值天气预报混合资料同化方法，包括以下步骤：

步骤1，建立基于三维变分资料同化方法的第一感知机模型，并对第一感知机模型进行训练；

步骤2，建立基于集合卡尔曼滤波器资料同化方法的第二感知机模型，并对第二感知机模型进行训练；

步骤3，建立混合同化的第三感知机模型；

步骤4，根据第一感知机模型和第二感知机模型的输出对第三感知机模型进行训练；

步骤5，使用训练好的第一感知机模型、第二感知机模型和第三感知机模型来计算获得分析场数据。

所述的第一感知机模型、第二感知机模型和第三感知机模型均为前馈神经网络模型构成，第一感知机模型和第二感知机模型的输入均为背景场序列和观测值序列，输出分别为三维变分资料同化方法对应的分析场和集合卡尔曼滤波器资料同化方法对应的分析场，所述的第一感知机模型和第二感知机模型的输出作为第三感知机模型的输入，第三感知机模型的输出为优化后的分析场。

具体地，所述的第一感知机模型的训练过程包括以下步骤：

步骤101，准备第一训练数据集，包括背景场序列、观测值序列和对应的分析场；

所述第一训练数据集中的背景场序列和观测值序列长度均为6，即采用6次背景场和6次观测值进行资料同化，对应的分析场表示为：

其中

表示t时刻的分析场，

表示t-i时刻的背景场，

表示t-i时刻的观测值，3Dω_t-i和3Dv_t-i分别是

和

的权重，bias_3D表示偏差；

所述第一训练数据集中对应的分析场是通过三维变分资料同化方法计算所得；

步骤102，使用所述的训练数据集训练所述的第一感知机模型。

具体地，所述的第二感知机模型的训练过程包括以下步骤：

步骤201，准备第二训练数据集，包括背景场序列、观测值序列和对应的分析场；

所述第二训练数据集中的背景场序列和观测值序列长度均为6，即采用6次背景场和6次观测值进行资料同化，对应的分析场表示为：

其中

表示t时刻的分析场，

表示t-i时刻的背景场，

表示t-i时刻的观测值，Enω_t-i和Env_t-i分别是

和

的权重，bias_En表示偏差；

所述第二训练数据集中对应的分析场是通过集合卡尔曼滤波器资料同化方法计算所得；

步骤202，使用所述的训练数据集训练所述的第二感知机模型。

具体地，所述的三维变分资料同化的分析场计算公式为：

x^a＝x^b+[B^-1+H^TR^-1h]^-1H^TR^-1(y-h(x^b))

其中，x^a表示分析场，x^b表示背景场，B表示背景误差协方差矩阵，H是线性观测算子h的切线性算子，R表示观测误差协方差矩阵，y表示观测值；所述的三维变分资料同化的分析场计算公式通过适当变换获得目标函数的梯度信息，随后采用最优化方法迭代求解分析场计算公式的近似极小值。

更进一步地，所述的第一感知机模型、第二感知机模型和第三感知机模型通过数值实验测定隐藏层个数以及每个隐藏层中的神经元数目，在所述的多层感知机中以ReLU为激活函数得到非线性结果，其损失函数为均方根误差MSE。

更进一步地，所述的第一感知机模型、第二感知机模型训练时学习率η＝0.0005，当损失值达到10^-5时训练停止，所述的第三感知机采用无隐藏层结构，直接由输出层加权得到输出层，训练时权重从0.1到1，递增率为0.1，两个权重占比之和为1。

本发明方法是集合三维变分资料同化方法和集合卡尔曼滤波器资料同化方法的混合资料同化方法。在传统的分析周期的基础上，获取由背景场和观测值序列组成的训练数据集，分别利用第一感知机模型和第二感知机模型来优化两类分析场结果，此外还利用第三感知机模型来对第一感知机和第二感知机的输出进行混合优化，将三维变分资料同化方法和集合卡尔曼滤波器资料同化方法的分析场耦合起来。本发明方法充分考虑了大气变量的流依赖，以此引入大气数据本身所包含的时间特征，同时又综合利用了两种传统资料同化方法的分析场结果，最终使得在数值天气预报中的资料同化效果更佳。

附图说明

图1是本发明的流程示意图；

图2是本发明实施例的数据处理过程示意图；

图3是本发明实施例中多层感知机的结构示意图；

图4Lorenz-63模型经验混合同化方法的结果图；

图5是本发明实施例中Lorenz-63模型采用两种混合同化方法的结果图；

图6Lorenz-96模型经验混合同化方法的结果图；

图7是本发明实施例中Lorenz-96模型采用两种混合同化方法的结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例所用到的部分字母符号注释见表1。

表1本文关键符号注释

在数值天气预报的资料同化方法发展过程中，基于Bayes定理，变分方法在20世纪90年代以前已经是数值天气预报的主流业务方法。通常，3D-Var(三维变分方法)通过最小化代价函数J来估计数值天气预报模型的单个初始状态，找到分析场。代价函数J表示为：

其中：

J_b用于测量x和x^b之间的偏差，即使x基于背景误差协方差B拟合x^b。

J_o用于测量x和x^o之间的偏差，即使x基于观测误差协方差R拟合y^o，其中h为观测算子，它完成必要的内插并将模式变量转换为观测空间场值。

利用非线性共轭梯度法求出式(1)的最小值，即需满足式(2)的雅可比矢量式：

其中H是观测算子h的切线性算子，定义为：

根据式(2)，x^a是使函数J最小化的状态，并由式(4)解得：

[B^-1+H^TR^-1h]x^a＝B^-1x_b+H^TR^-1y (4)

则最优分析场为：

极小化代价函数J可以看作是一个反问题，由于误差协方差矩阵B、R的构造和求逆运算是三维变分同化的难点，且高维矩阵B、R的求逆是非常困难的，需要大量的计算资源和存储空间，因此在实际解决这个问题时，并不直接求解式(5)，而是通过适当变换获得目标函数的梯度信息，随后采用最优化方法(如牛顿法、拟牛顿法、BFGS法、LBFGS法等)迭代求解式(5)的近似极小值。

集合卡尔曼滤波器(EnKF)方法中，包括主要以下几步。

预报步：

分析步：

其中：

P^f是同时进行的N个资料同化系统的集合估计,即在完成t_i-1的分析集合和N个预测值

之后，P^f可以近似地由式(9)得到：

式(9)中的

表示样本平均值。

感知机是一个线性分类器，即对线性可分的数据集，通过感知机可以学习拟合到较好的线性方程。然而，单神经元的感知机模型表达能力不足以解决除简单线性回归外更复杂的问题。多层感知器又称为前馈神经网络(ANN)，具有全连接的结构，除了输入层外，每个节点都可以看作具有非线性激活函数的处理单元。如果每一个激活函数都是线性函数，则任意层中多层感知器的每个节点都可以简化为感知器。即多层感知器是一种经典的深度神经网络，它是一类通用的逼近器。一个多层感知器如果具有线性输出层和至少一层提供非线性变换的激活函数的隐藏层，给定网络各层合适的神经元数量，MLP可以实现高精度地拟合可测函数，即多层感知器的目标是拟合某个函数f(x)。

在数值天气预报的业务系统中，使用短期预报作为背景场(x^b)并结合新观测(y^o)的信息来获得分析场(x^a)被称为一个资料同化分析周期。定义资料同化的问题为，在时间为t时，根据式(10)找到耦合两个信息的最佳权重矩阵(ω)。

由此，如图1所示，本实施例提出基于三重多层感知机的数值天气预报混合资料同化方法，包括以下步骤：

步骤1，建立基于三维变分资料同化方法的第一感知机模型，并对第一感知机模型进行训练；

步骤2，建立基于集合卡尔曼滤波器资料同化方法的第二感知机模型，并对第二感知机模型进行训练；

步骤3，建立混合同化的第三感知机模型；

步骤4，根据第一感知机模型和第二感知机模型的输出对第三感知机模型进行训练；

步骤5，使用训练好的第一感知机模型、第二感知机模型和第三感知机模型来计算获得分析场数据。

所述的第一感知机模型、第二感知机模型和第三感知机模型均为前馈神经网络模型构成，第一感知机模型和第二感知机模型的输入均为背景场序列和观测值序列，输出分别为三维变分资料同化方法对应的分析场和集合卡尔曼滤波器资料同化方法对应的分析场，所述的第一感知机模型和第二感知机模型的输出作为第三感知机模型的输入，第三感知机模型的输出为优化后的分析场。所述的多层感知机的结构如图3所示；

具体地，所述的第一感知机模型的训练过程包括以下步骤：

步骤101，准备第一训练数据集，包括背景场序列、观测值序列和对应的分析场；

所述第一训练数据集中的背景场序列和观测值序列长度均为6，即采用6次背景场和6次观测值进行资料同化，对应的分析场表示为：

其中

表示t时刻的分析场，

表示t-i时刻的背景场，

表示t-i时刻的观测值，3Dω_t-i和3Dv_t-i分别是

和

的权重，bias_3D表示偏差；

所述第一训练数据集中对应的分析场是通过三维变分资料同化方法计算所得；

步骤102，使用所述的训练数据集训练所述的第一感知机模型。

具体地，所述的第二感知机模型的训练过程包括以下步骤：

步骤201，准备第二训练数据集，包括背景场序列、观测值序列和对应的分析场；

所述第二训练数据集中的背景场序列和观测值序列长度均为6，即采用6次背景场和6次观测值进行资料同化，对应的分析场表示为：

其中

表示t时刻的分析场，

表示t-i时刻的背景场，

表示t-i时刻的观测值，Enω_t-i和Env_t-i分别是

和

的权重，bias_En表示偏差；

所述第二训练数据集中对应的分析场是通过集合卡尔曼滤波器资料同化方法计算所得；如图2所示。

步骤202，使用所述的训练数据集训练所述的第二感知机模型。

具体地，所述的三维变分资料同化的分析场计算公式为：

x^a＝x^b+[B^-1+H^TR^-1h]^-1H^TR^-1(y-h(x^b))

其中，x^a表示分析场，x^b表示背景场，B表示背景误差协方差矩阵，H是线性观测算子h的切线性算子，R表示观测误差协方差矩阵，y表示观测值；所述的三维变分资料同化的分析场计算公式通过适当变换获得目标函数的梯度信息，随后采用最优化方法迭代求解分析场计算公式的近似极小值。

更进一步地，所述的第一感知机模型、第二感知机模型和第三感知机模型通过数值实验测定隐藏层个数以及每个隐藏层中的神经元数目，在所述的多层感知机中以ReLU为激活函数得到非线性结果，其损失函数为均方根误差MSE。

更进一步地，所述的第一感知机模型、第二感知机模型训练时学习率η＝0.0005，当损失值达到10^-5时训练停止，所述的第三感知机采用无隐藏层结构，直接由输出层加权得到输出层，训练时权重从0.1到1，递增率为0.1，两个权重占比之和为1。

为了验证本发明方法的有效性，本实施例利用两个经典的非线性混沌动力学系统作为验证模型，洛伦兹63模型和洛伦兹96模型。

洛伦兹63模型(Lorenz-63系统)是一个简化的大气对流数学模型，具有三个变量之间强非线性相互作用的优点。Lorenz-63系统的状态演化可用三个常微分方程表示：

其中产生一个混沌状态的标准参数设置为：

假设Lorenz-63系统是完美的，不存在模型误差，即模型误差协方差矩阵Q＝0；给定初始状态x₀＝1.508870,y₀＝-1.537121,z₀＝25.46091，并且Lorenz-63模型使用四阶Runge-Kutta时间差分格式进行时间积分，时间步长d_t＝0.01，在模拟积分步长111000下，可以通过积分公式(13)中的方程来获得真实场(x^t)。在实际大气模型中，变量的地面真值通常用来表示x^t。在真实状态中加入高斯分布的随机观测误差，以近似观测值。此外，我们还设置了线性观测算子h＝I和观测误差协方差矩阵R＝I，其中I是三阶单位矩阵。同样，在相同的初始条件下，采用传统的资料同化方法3D-Var(或EnKF)的同化系统分别完成111000次同化分析周期，生成分析场(x^a)。其中，为达到Lorenz-63系统的稳定状态和获得近似模拟大气，需进行前1000个积分时间步的spinning up；中间10万步用于准备Cache-MLP的训练集，最后1万步为测试数据。

此外，已知Lorenz-63混沌系统的初始状态，可以用NMC方法粗略地计算3D-Var的静态背景误差协方差B。

RMSE在数据资料同化中被广泛用于测量x^t和x^a之间的相似性，定义为：

其中i是第i个状态变量，j是第j个集合成员。

通过对第一感知机模型和第二感知机模型的训练，对训练数据集进行优化，为第三感知机模型做准备，并对优化过程进行重复，取得了明显的效果。为简化计算过程，选定3D-Var与EnKF分析场耦合的经验权重为一个固定的比例，即从0.1到1，递增率为0.1，两者占比之和为1，即比例(3DVar)+比例(EnKF)＝1。可知，当3D Var的x^a占比为10％，来自EnKF的x^a占比为0.9时，可以得到RMSE的最小平均值为0.1386，如图4。

但图5(a)和图5(b)相比，本发明方法通过神经网络学习两种方法之间耦合的最优权值矩阵，优于数理统计的最佳性能。得到平均RMSE为0.1243，性能提高了10.32％。

同样，将本发明方法在洛伦兹96模型(Lorenz-96系统)中进行验证。与Lorenze-63相比，Lorenz-96是一个更加复杂的非线性动力系统。其函数可定义为式(15)。其中，右边第一项为平流项，第二项为阻尼项。此外，F代表一个外部强迫常数，设置为8，这是一个已知的形成混沌行为的一般值。Lorenz-96是一个用于模拟大气变量时间演化的周期系统。

其中J是状态变量数。

基于数学经验统计方法，比例(3DVar)+比例(EnKF)＝1中，为了获得分析场的最佳性能，可以看出，当赋予3D-Var 30％的权重时，RMSE的值最小，值为0.2943，如图6所示。

然而，利用本发明方法，根据图7(a)和图7(b)，我们可以得到50次实验的较低的平均RMSE值为0.2530，其性能提升了14.03％。

由此可知，经验混合的资料同化方法与本发明方法的性能对比如表2所示。

表2经验混合资料同化法与本发明方法的性能比较

模型	经验混合资料同化法的RMSE	本发明方法的RMSE	性能提升(％)
				Lorenz-63	0.1368	0.1243	10.32
Lorenz-96	0.2943	0.2530	14.03

由发明内容和实施例内容可知，本发明方法利用大气运动的物理规律和大气状态变量的时间特性，对传统的资料同化方法进行模拟和优化，使得本发明方法的同化性能较传统方法得到有效提升。本发明方法是在传统的分析周期的基础上，获取由背景场和观测值组成的训练数据集，训练神经网络来优化分析场，充分考虑了大气变量的流依赖，神经网络模型的输出值不依赖于单一时间的数据，而是与前五个时间步长数据的有关，以此引入大气数据本身所包含的时间特征，使得在数值天气预报中的本发明方法资料同化效果更佳。

Claims (4)

1.基于三重多层感知机的数值天气预报混合资料同化方法，其特征在于，包括以下步骤：

步骤1，建立基于三维变分资料同化方法的第一感知机模型，并对第一感知机模型进行训练；

步骤2，建立基于集合卡尔曼滤波器资料同化方法的第二感知机模型，并对第二感知机模型进行训练；

步骤3，建立混合同化的第三感知机模型；

步骤4，根据第一感知机模型和第二感知机模型的输出对第三感知机模型进行训练；

步骤5，使用训练好的第一感知机模型、第二感知机模型和第三感知机模型来计算获得分析场数据；

所述的第一感知机模型、第二感知机模型和第三感知机模型均为前馈神经网络模型构成，第一感知机模型和第二感知机模型的输入均为背景场序列和观测值序列，输出分别为三维变分资料同化方法对应的分析场和集合卡尔曼滤波器资料同化方法对应的分析场，所述的第一感知机模型和第二感知机模型的输出作为第三感知机模型的输入，第三感知机模型的输出为优化后的分析场；

所述的第一感知机模型的训练过程包括以下步骤：

步骤101，准备第一训练数据集，包括背景场序列、观测值序列和对应的分析场；

所述第一训练数据集中的背景场序列和观测值序列长度均为6，即采用6次背景场和6次观测值进行资料同化，对应的分析场表示为：

其中，

表示t时刻的分析场，

表示t-i时刻的背景场，

表示t-i时刻的观测值，3Dω_t-i和3Dv_t-i分别是

和

的权重，bias_3D表示偏差；

所述第一训练数据集中对应的分析场是通过三维变分资料同化方法计算所得；

步骤102，使用所述的训练数据集训练所述的第一感知机模型；

所述的第二感知机模型的训练过程包括以下步骤：

步骤201，准备第二训练数据集，包括背景场序列、观测值序列和对应的分析场；

所述第二训练数据集中的背景场序列和观测值序列长度均为6，即采用6次背景场和6次观测值进行资料同化，对应的分析场表示为：

其中

表示t时刻的分析场，

表示t-i时刻的背景场，

表示t-i时刻的观测值，Enω_t-i和Env_t-i分别是

和

的权重，bias_En表示偏差；

所述第二训练数据集中对应的分析场是通过集合卡尔曼滤波器资料同化方法计算所得；

步骤202，使用所述的训练数据集训练所述的第二感知机模型。

2.根据权利要求1所述的数值天气预报混合资料同化方法，其特征在于，所述的三维变分资料同化的分析场计算公式为：

x^a＝x^b+[B^-1+H^TR^-1h]^-1H^TR^-1(y-h(x^b))

其中，x^a表示分析场，x^b表示背景场，B表示背景误差协方差矩阵，H是线性观测算子h的切线性算子，R表示观测误差协方差矩阵，y表示观测值；所述的三维变分资料同化的分析场计算公式通过变换获得目标函数的梯度信息，随后采用最优化方法迭代求解分析场计算公式的近似极小值。

3.根据权利要求1或2所述的数值天气预报混合资料同化方法，其特征在于，所述的第一感知机模型、第二感知机模型和第三感知机模型通过数值实验测定隐藏层个数以及每个隐藏层中的神经元数目，在所述的多层感知机中以ReLU为激活函数得到非线性结果，其损失函数为均方根误差MSE。

4.根据权利要求3所述的数值天气预报混合资料同化方法，其特征在于，所述的第一感知机模型、第二感知机模型训练时学习率η＝0.0005，当损失值达到10^-5时训练停止，所述的第三感知机采用无隐藏层结构，直接由输出层加权得到输出层，训练时权重从0.1到1，递增率为0.1，两个权重占比之和为1。

Images