从太空卫星到社交媒体聊天，从信用卡交易业务到医疗保健记录，都得益于先进的传感技术、强大的计算平台和广泛的在线连接技术，日益增加的数据量极大地改变了现代处理和分析数据的计算和统计方法。与此同时，基于人工智能和机器学习的数据分析开始在推进科学发现和工程设计方面发挥出重要作用。过去，科学的进步是先提出假设，然后收集数据来证实或否定这些假设。然而，在大数据时代，虽然数据被收集了起来，却对其内部运作一无所知，导致形成了所谓的黑箱解决方案。这种观点根本不正确：海量数据可以使得在没有科学理解的情况下建立可操作的模型成为可能。理想情况下，将挖掘的模式信息和内在关系转化为可解释的理论和假设，从而推动科学认知的进步，这才是最重要的。因此，即使黑盒模型在特定的空间、时间和频谱情景下实现了某种程度的精确性能，但仍缺乏对底层演变过程作用机理的理解能力，也就不能可靠地成为后续科学和工程发展的基础。这一点尤其适用于当今地球物理流体动力学或热动力学领域。地球物理流体动力学，广义上说，指的是在地球和其他行星上自然发生的流动(例如熔岩流、海洋和行星大气流)的流体动力学。

在工程设计、控制和建模中，存在着各种各样的任务都需要求解基于偏微分方程(Partial Differential Equation, PDE)的多集合前向模型。例如，在顺序数据同化中，前向模型的多个集合被用来逼近协方差矩阵^[1]。随着所需精度要求的提高，时空数值模拟离散化的粒度也需要适当提升。这将导致计算工作量的显著增加，并可能成为设计和预测周期的外部循环中的瓶颈。例如，在计算流体力学(Computational Fluid Dynamics, CFD)中，湍流的三维模拟^[2]常被用于获得Navier-Stokes方程的精确解，但在形状优化等相关任务中却使用得较少^[3]。因此，建立一个计算代价比全阶模型(Full Order Models, FOM)更低的降阶模型(Reduced Order Model, ROM)，并提供足够精度水平的解决方案的研究得到了很多人的关注^[4-6]。

近年来，人们对开发用于物理系统的非侵入式ROM的研究兴趣与日俱增^[7-9]。“非侵入式”指的是只通过数据来构造降阶模型。非侵入式降阶模型特别适用于理想模型未知，但存在大量数据可用于模型恢复的系统。这种情况在地球物理流体中非常普遍。由于网格分辨率不足而导致的子网格尺度的近似、模型参数的不确定性或模型本身结构不正确等原因，使得地球物理系统模型可能是不完善的^[10]。然而，从遥感、卫星观测和地球物理流体实验测量中已经获得了大量的实测数据。因此，为了有效地预测地球物理过程，利用这些数据已经开展了多项研究^[11-12]。

在本研究中，主要使用线性降维技术并配合机器学习算法来演化隐空间^[13-14]。具体来说，就是利用本征正交分解(Proper Orthogonal Decomposition, POD)来识别全阶模型的隐空间，然后利用长短期记忆(Long Short-Term Memory, LSTM)神经网络来模拟隐空间的演化过程。这种替代建模技术的主要优点之一是它由纯数据驱动，因此特别适合地球物理流体的研究数据集。利用遥感和现场观测收集的档案数据建立替代模型，然后将替代模型用于完成预测任务。当获得新数据之后，可以使用迁移学习方法^[15]重新训练替代模型，从而提高预测精度。

除了ROM之外，神经网络在许多科学研究中的使用也显著增加^[16-17]。神经网络的主要挑战之一是其性能过于依赖其网络结构^[18]。此外，神经网络包含大量与问题相关的超参数。通常，人工神经网络是通过试错法来设计的，这一过程可能相当耗时。此外，要获得良好的性能，还需要设计者具有深刻的数据理解和丰富的神经网络领域知识。虽然有诸如“网格搜索”或“随机搜索”这样的方法来找到超参数的优化组合，但是随着搜索空间的维数增加，这些方法的可扩展性欠佳。因此，人们对神经网络结构和超参数搜索的自动化设计相关的内容越来越感兴趣，这将允许那些没有神经网络专业知识的用户可以将其应用于他们感兴趣的特定研究问题^[19]。本文采用遗传算法对LSTM神经网络的结构设计和超参数进行优化。将LSTM网络用于NOAA海表温度数据集的替代模型设计。

本文其余部分的结构如下。首先，在第1节描述了非侵入式降阶模型的建模方法和数据预处理。然后，在第2节中给出了所提出算法的细节。第3节讨论了实验结果及其分析。最后，在第4节中，对本文的结论和未来的研究方向进行了展望。

1 非侵入式降阶模型 1.1 本征正交分解

本文使用POD来提取表征上述非线性动力系统的主导模式。在不同的时刻收集了(动态网络的)数据快照 ${{\boldsymbol{u}}_1},{{\boldsymbol{u}}_2}, \cdots ,{{\boldsymbol{u}}_N} \in {{\bf{R}}^M}$ ，其中 $M$ 为空间自由度， $N$ 等于网格点总数，为数据快照的个数。利用本征正交分解，本文构造了一组标准正交基函数来优化描述系统的场变量。则快照的数据矩阵形式如式(1)所示。

$ {\boldsymbol{A}}=\left[ {{{\tilde {\boldsymbol{u}}}_1}\left| {{{\tilde {\boldsymbol{u}}}_2}} \right.\left| \cdots \right.\left| {{{\tilde {\boldsymbol{u}}}_N}} \right.} \right] \in {{\bf{R}}^{M \times N}} $

(1)

式(1)中，去均值(异常)场的计算公式如式(2)所示。

$ {\tilde {\boldsymbol{u}}_\iota }={{\boldsymbol{u}}_i} - \bar {\boldsymbol{u}},\;\;\bar {\boldsymbol{u}}=\frac{1}{N}\sum_{i=1}^N {{{\boldsymbol{u}}_i}} $

(2)

式(2)中， $\bar {\boldsymbol{u}}$ 是解域的时间均值。当构造出快照数据矩阵之后，本文使用奇异值分解(Singular Value Decomposition, SVD)计算矩阵 ${\boldsymbol{A}}$ 的左、右奇异向量。在矩阵形式中，SVD可以写成式(3)。

$ {\boldsymbol{A}}={\boldsymbol{W}}\sum {{{\boldsymbol{V}}^{\rm{T}}}} =\sum_{\kappa =1}^N {{\sigma _\kappa }{{\boldsymbol{w}}_\kappa }v_\kappa ^{\rm{T}}} $

(3)

式(3)中， ${\boldsymbol{W}} \in {{\bf{R}}^{M \times N}}$ ， ${\boldsymbol{V}} \in {{\bf{R}}^{M \times N}}$ 。 ${\boldsymbol{W}}$ 和 ${\boldsymbol{V}}$ 分别包含与 ${\boldsymbol{A}}{{\boldsymbol{A}}^{\rm{T}}}$ 和 ${{\boldsymbol{A}}^{\rm{T}}}{\boldsymbol{A}}$ 的特征向量相同的左奇异向量和右奇异向量。此外，奇异值的平方等于特征值，即 ${\lambda _\kappa }=\sigma _\kappa ^2$ 。向量 ${{\boldsymbol{w}}_\kappa }$ (也是 ${\boldsymbol{A}}{{\boldsymbol{A}}^{\rm{T}}}$ 的特征向量)是本征正交分解的基函数，本文将它们记为 ${{\boldsymbol{\phi}} _\kappa }$ 。POD的基函数是正交的(即 $\left\langle {{{\boldsymbol{\phi}} _i},{{\boldsymbol{\phi}} _j}} \right\rangle ={{\boldsymbol{\delta}} _{ij}}$ )，并且在 ${\rm{L}}2$ 范数下可被精确求解^[20-21]。则动力系统的状态可以用这些POD基函数来近似表征，如式(4)所示。

$ u(x,t)=\bar {\boldsymbol{u}} + \sum_{\kappa =1}^R {{a_\kappa }(t){{\boldsymbol{\phi}} _\kappa }(x)} $

(4)

式(4)中， $R$ 是使得 $R \ll N$ 的保留基函数的个数， ${a_\kappa }$ 为随时间变化的模型系数。POD的基函数最小化了场变量与其断面表示之间的均方误差。此外，它还使描述给定误差条件下的场变量所需的基函数的数目最小化。保留模式的数目通常是根据它们的场能含量来决定的。利用这些保留模式，可以构造本征正交分解基组 $\varPhi =\left\{ {{{\boldsymbol{\phi}} _\kappa }} \right\}_{\kappa =1}^R$ ，从而建立ROM。

1.2 ROM的无模型演化

特别是在地球物理动力系统中，ROM演化的无模型预测具有广阔的应用前景。对于各种各样的任务，如数据同化和地球物理流体的不确定性量化，需要运行前向模型的集合，伴随着巨大的计算代价。可以使用基于ROM的替代模型来取代前向模型的演化过程。对于许多地球物理动力系统，由于网格分辨率颗粒度较为粗糙，动力系统的力学描述是不可用的或不足的。然而，在过去的几十年里，从局部和卫星观测中获得了大量的数据，数据驱动的方法如递归神经网络(Recursive Neural Network, RNN)在时空混沌动力系统的无模型预测中被证明是成功的^[22-23]。

式(4)中的时变模态系数是通过在POD基础上的投影平均值减去场能求得的，如式(5)所示。

$ {a_\kappa }({t_i})=\left\langle {{{\tilde {\boldsymbol{u}}}_\iota },{{\boldsymbol{\phi}} _\kappa }} \right\rangle $

(5)

式(5)中，尖括号表示两个函数的内积。因此，如果能够对模态系数的演化进行预测，就可以重构未来的解域。一种流行的时间序列预测方法是递归神经网络，近来已被应用于物理系统非侵入式降阶模型的若干研究^{[13, 24-26]}。这是因为RNN可以通过系统的历史信息来进行将来的状态预测，以及通过可用的历史数据来训练神经网络(参数)。LSTM神经网络是解决消失梯度问题的RNN最成功的变体之一^[27]。本研究中，使用 ${{\boldsymbol{X}} _\kappa }$ 表示输入序列数据矩阵， ${{\boldsymbol{Y}}_\kappa }$ 表示输出序列数据矩阵。将输入训练矩阵 ${{\boldsymbol{X}} _\kappa }$ 的每个样本构造为 $\{ a_1^{(n)}, \cdots ,a_R^{(n)}; \cdots ; $ $ a_1^{(n - \sigma + 1)}, \cdots ,a_R^{(n - \sigma + 1)} \}$ ，输出序列数据矩阵 ${{\boldsymbol{Y}}_\kappa }$ 中相应的输出样本为 $\{ {a_1^{(n + 1)}, \cdots ,a_R^{(n + 1)}} \}$ 。参数 $\sigma $ 被称为回溯时间窗，它定义了LSTM模型进行训练和预测的递归时间步长。

1.3 数据预处理

无模型预测对于现实世界的实测数据尤其重要，因为实测数据是由多尺度过程产生的，这些过程不能用任何基于模型的方法精确地逼近。现实世界中的流体动力过程是由多个外部系统耦合控制的，而这种耦合可能是未知的或无法建模的。因此，现有模型对现实世界中的数据描述效果欠佳，非侵入式方法则非常适合这些复杂的流体。为此，本文研究了非侵入式降阶模型在NOAA最优差值海表温度数据集(第2版)上的应用。该数据集由一个分辨率为1°的网格上的每周平均海表温度数据快照组成，并结合卫星和本地测量生成。季节性波动使得该数据集的温度场具有很强的周期性结构。

沿纬度和经度以1°作为网格分辨率，数据集的每个快照的维数为180×360。利用掩模运算去除与陆地相对应的数据点，并建立了仅针对海洋表面的归一化数据的替代模型。该数据的时间跨度为1981年10月22日至2018年6月30日(即共1 914个快照)。

数据预处理完成后，利用去均值(异常)场的温度数据生成快照数据矩阵。ROM的模式数量是根据式(6)所示的相对信息含量(Relative Information Content, RIC)求得

$ {\rm{RIC}}(R)=\frac{{\displaystyle\sum_{\kappa =1}^R {\sigma _\kappa ^2} }}{{\displaystyle\sum_{\kappa =1}^N {\sigma _\kappa ^2} }} $

(6)

RIC表示可使用 $R$ 基函数恢复的总体数据的信息(方差)比例。图1显示了NOAA最优差值海表温度数据集的RIC百分比。本文将保留模式的数量固定为 $R=8$ ，它捕获了约92%的数据信息(方差)，并且这些模式足以捕获NOAA数据集中的季节趋势。可以看出，在8个模式之后，模式数量的增加使得总信息方差的增加非常小，这是由于这些模式主要负责捕捉小尺度的波动。

2 LSTM的遗传算法优化

神经网络的性能在很大程度上依赖于神经网络结构的设计和其他超参数的选择，例如激活函数、优化器、初始权重值等。类似于其他超参数优化的研究^[28-34]，本文使用不同类型的小型结构对整个LSTM架构进行编码。图2展示了作为LSTM主体网络组成构件的不同类型的小型体系结构。在这些小型结构中，每一个都采用了残差连接，因为它允许训练深度神经网络，而不会导致梯度消失的问题^[35-36]。最终的LSTM网络是通过对这些组件进行排序来构建的。LSTM存储单元数和模块层数是与使用遗传算法优化的神经网络结构设计相关的另外两个参数。需要优化的其他超参数包括优化器的类型、优化器的学习率、权值和偏差的初始化分布以及激活函数。

算法1列出了本文所提框架的伪代码。它首先初始化指定大小的填充。根据预定义的组件构造个体的LSTM子网络，并在此基础上分配其他超参数。在训练数据集上采用三重交叉验证对神经网络进行评估。一旦计算出种群中每个个体的适应度(即，本研究中的验证均方误差)，就会根据适应度值对个体进行排序。在这个种群当中，保留一定比例的精英个体(作为亲本)，通过联赛选择算法竞争产生下一代个体。利用交叉和变异算子从选择的双亲中产生子代个体。新种群是由老种群中保留下来的精英个体与子代个体结合而成的。因此，新种群是由前一次迭代的精英个体和新产生的子代组成的。增加遗传代数计数器，并在给定的指定遗传代数内重复上述过程。

算法1 本文提出的遗传算法

输入: 一组预定义的构件，种群规模，最大遗传代数N_g，

训练数据集，精英保留百分比。

输出: 找到的LSTM网络结构以及其他超参数。

1. $ { {P}_{0}\leftarrow }$ 以给定规模初始化种群；

2. $ {k\leftarrow 0} $ ；

3. 从 $ {k=0}$ 到最大遗传代数，循环以下操作：(for循环)

4. 计算种群当中每个个体的适应度 $ {{P}_{k} }$ ;

5. $ { {P}_{k}^{E}\leftarrow }$ 保留具有最大适应度值的精英群体;

6. $ { {P}_{k}^{T}\leftarrow }$ 使用联赛选择法挑选产生子代种群的亲本个体；

7. $ {{Q}_{k}\leftarrow} $ 使用本文提出的交叉和变异算子从所选亲本当中产生后代个体；

8. $ { {P}_{k+1}\leftarrow {P}_{k}^{E}\cup {Q}_{k} }$ ；

9. $ {k\leftarrow k+1 }$ ；

10. 结束循环。

11. 返回种群 $ { {P}_{k} }$ 产生的适应性最好的个体。

使用7个整数来定义种群中的每个个体。这些数字代表了LSTM网络的超参数。这些超参数的上、下限以及用于生成这些数字的分布类型参见表1。第1个参数是网络的组块(数量)，图2展示了4种不同类型的可用于设计LSTM网络的组块结构。第2和第3个参数与LSTM网络的深度(即隐藏层的数量)和宽度(即LSTM的cell单元的数量)有关。其他参数如激活函数、优化器、学习率、初始权重值等都与LSTM网络的训练有关。在评估过程中，针对每个个体对应的LSTM网络进行训练，并将验证均方误差(Mean Squared Error, MSE)作为该个体的适应度。种群根据适应度值按升序排序(MSE越低，适应度越好)。

下一步是为子代个体选择亲本。本文采用联赛选择算法来挑选亲本。算法2给出了联赛选择算法的选择过程。在联赛选择过程中，从种群中随机抽取少数个体，选出其中的最佳个体用于产生子代。

算法2 亲本选择

输入: 种群 $ { {P}_{k}} $ , 种群规模N, 精英个体数量M, 竞争者规模T。

输出： 后代种群 $ {{P}_{k}^{T} }$ 。

1. $ { {P}_{k}^{T}\leftarrow \varnothing} $ ；

2. 从 $ {i=0}$ 到 $ {N - M}$ ，循环以下操作：(for循环)

3. T $ { \leftarrow }$ 随机生成不重复的数字T;

4. $ { {P}_{k}^{T}\leftarrow {P}_{k}^{T}\cup {P}_{k}\left[\mathrm{min}\left(T\right)\right]; }$

5. 结束循环。

6. 返回选取的种群 $ { {P}_{k}^{T}} $ 。

算法3给出了生成子代的细节。子代种群的产生分为两个阶段。第一阶段是交叉运算(第2~13行)，这里使用均匀交叉算子。在均匀交叉运算中，从种群当中选择两个个体。然后，针对个体的每个基因生成一个随机数，如果该随机数大于交叉概率，则交换两个个体在位置 $j$ 处的对应基因。在突变过程中，如果随机数大于突变概率，则改变位置 $i$ 处的基因。

算法3 子代生成

输入: 包含适应度的被选取种群 $ { {P}_{k}^{T}} $ ，交叉概率 $ { {p}_{c}} $ ，变异概率 $ {{p}_{m} }$ 。

输出： 子代种群 $ { {Q}_{k} }$ 。

1. $ {{Q}_{k}\leftarrow \varnothing }$ ；

2. 当 $ { i < \left|{P}_{k}^{T}\right| }$ 时，执行以下循环操作：(while循环)

3. $ { {p}_{1},{p}_{2}\leftarrow }$ 从 $ {{P}_{k} }$ 中选择2个连续的个体；

4. 对于 $ {j=0}$ 到 $ {\left| {{p_1}} \right|}$ ，循环以下操作：(for循环)

5. $ { r\leftarrow }$ 随机生成一个介于 $ {\left[\mathrm{0,1}\right] }$ 之间的数字;

6. 如果 $ { r > {p}_{c} }$ ，则：

7. 交换 $ { {p}_{1}} $ 和 $ { {p}_{2}} $ 的基因；

8. 结束条件判断。

9. $ {j\leftarrow j+1} $ ；

10. 结束for循环。

11. $ { {Q}_{k}\leftarrow {Q}_{k}\cup {p}_{1}\cup {p}_{2}} $ ；

12. $ { i\leftarrow i+1 }$ ；

13. 结束while循环。

14. 对于 $ { {Q}_{k}} $ 中的每个个体p，循环以下操作：(for循环)

15. $ { r\leftarrow }$ 随机生成一个介于 $ {\left[\mathrm{0,1}\right] }$ 之间的数字;

16. 如果 $ {r < {p}_{m}} $ ，则：

17. $ { i\leftarrow }$ 在p中随机选择一个点;

18. 改变p的i点处的基因;

19. 结束条件判断。

20. 结束for循环。

21. 返回生成的子代种群 $ {Q}_{k} $ 。

3 实验与分析

在本节中，给出了遗传算法在优化神经网络结构设计和超参数方面的结果。然后，展示了优化后的LSTM网络在NOAA数据集仿真和预测中的性能。

NOAA数据集的时间跨度为1981年10月22日至2018年6月30日，对应1 914个样本快照。本文利用前1500个数据快照中随机选择的70%样本数据作为训练数据集。在LSTM网络的超参数优化过程中，对种群内的个体进行了100个epoch的训练(一个epoch即是将所有训练样本训练一次的过程，epoch数为100)，并采用三次交叉验证来避免过拟合。将LSTM网络的回溯时间窗口的步长设置为8，用于捕捉模态系数之间的时间相关性。每个个体的适应度是三折交叉验证数据集的MSE的平均值。一次训练所选取的样本数(Batch Size)固定为64。本实验的种群规模和遗传代数设置为20。图3展示了种群内个体适应度随遗传算法的遗传代数的演化轨迹。箱线图用于显示每一代的种群统计数据。图3还呈现了每一代种群数量的MSE的中位数和最小值。随着演化过程的进行，验证集上的均方误差也逐渐减小。矩形盒的高度表示种群在每一代的适应度的方差，可以看出，第一代之后的方差显著降低。MSE从第一代到第二代急剧下降，这可以归因于遗传算法开始时种群的随机初始化。对于本文中所研究的问题，3~10代(遗传代数)似乎就足以找到最佳的超参数集。

利用遗传算法找到最优网络结构和其他超参数后，对最优LSTM网络进行了1200个epoch(即1200次遍历整个训练数据集)的训练，批量大小(一次训练的样本数)为64。在训练后的网络部署过程中，采用自回归法对模态系数进行预测。给出了前8个时间步(等于LSTM的回溯时间窗口长度)模态系数的初始条件。该信息用于预测第9个时间步的模态系数。然后利用第2到第9个时间步的模态系数来预测第10个时间步的模态系数。重复此过程直到最后一个时间步，即第1914个时间步。由于本文自回归部署使用的是经过训练的LSTM网络，所以只需要设置LSTM回溯所对应的初始条件。在经过(初期的)几个时间步后，仅利用LSTM的预测结果作为未来(时间步)的模态系数预测值。图4展示了模态系数的真实值和预测值。可以看出，尤其对于负责捕捉大尺度波动和季节模式的前几个模态而言，它们的模态系数的真实值和预测值之间具有良好的一致性。

图5描绘了两个不同时刻的真实温度场和重建温度场。利用式(4)重建温度场，其中平均温度场是由前1500个数据快照计算所得。基于LSTM的ROM(降阶模型)能够准确地捕捉温度场中的大尺度模式，展现出数据驱动的ROM在地球物理流体分析中的应用潜力。图6展示了L2范数下的温度场真实值和预测值之平方误差。可以观察到，在最后的时间步，预测误差增加，这可能是由于模态系数预测不准确，以及POD基函数无法在外推时间域(即时间步超过1500)捕获高精度的空间模式。针对LSTM误差累积的一种补救方法是以非自回归的方式来训练和部署网络^[37]。

4 结论

本文提出了一种用于长短记忆神经网络结构自动搜索和超参数优化的遗传算法，可用于地球物理流体代理模型的建模任务。并且成功地证明了利用优化后的LSTM建立的代理模型能够预测海表温度场。为了便于更深层次神经网络的训练，本文采用了一种编码策略，其中LSTM网络使用了包含残差连接的更小的构件来设计。优化后的LSTM网络可以在足够长的时间内准确地预测海温场的季节变化，而不存在任何不稳定问题。

本文还观察到，真实模态系数和预测模态系数之间的差值在预测期间比训练期间的更大。这个问题可以通过使用非自回归部署或迁移学习的方法来解决，其中，LSTM网络在新数据可用时将被重新训练。在目前的研究中，本文假设每个隐藏层中的LSTM单元数目是恒定的。在未来的研究工作中，本文将消除这一约束，并对种群中的每个个体使用变长编码策略来优化LSTM网络。