通过卷积神经网络增强的多层集成卡尔曼滤波算法

《Computers & Fluids》：Multi-level ensemble Kalman Filter algorithms enhanced by Convolutional Neural Networks

　　汤姆·穆西（Tom Moussie）| 保罗·埃兰特（Paolo Errante）| 马塞洛·梅尔迪（Marcello Meldi）

摘要
本研究提出了一种使用卷积神经网络（CNN）提升数据同化策略的方法。具体而言，通过CNN工具改进了集合卡尔曼滤波器的多层次算法，旨在减少不同模型预测结果之间的差异。该方法通过分析雷诺数Re=1000、马赫数Ma=0.5时NACA 0012型翼型的流动特性来评估其效果。根据攻角α的不同，可以观察到流动的非稳态特征。结果表明，使用基于数据同化（DA）过程训练的CNN工具可以显著提高低保真度模型的准确性，同时计算成本的增加很小。此外，CNN工具的使用还加快了数据同化算法的收敛速度，从而在计算资源需求方面带来了显著的优势。

1. 引言
开发新的数值技术和算法是实现多个科学领域实际应用中系统性技术进步的关键步骤。这一点在流体力学研究中尤为明显，因为该学科涉及许多研究案例，如非线性多尺度现象的研究及其在运输工程、医学研究、能量收集和城市环境中的应用。因此，提高这一学科预测工具的准确性将为广泛应用带来巨大的技术进步。对于像计算流体动力学（CFD）求解器这样的尺度解析工具来说，实现这一目标面临双重挑战。一方面，高雷诺数下的流动包含众多动态活跃尺度，所需的计算资源量巨大，而这些情况在工程中非常常见，且流动行为具有强烈的非稳态特性。可以考虑使用近似方法来减轻这种计算负担[2]。这在湍流建模中尤为有用。这些方法在雷诺平均纳维-斯托克斯（RANS）建模[3]和大涡模拟（LES）[4]等研究中得到了广泛研究，虽然具有相对较高的准确性，但对具体测试案例的敏感性极高。另一个问题是，即使所采用的数值模型具有完美的准确性，由于缺乏初始/边界条件的知识，也可能无法与实际情况进行比较[5]。由于高雷诺数流动的多尺度特性，瞬时值的微小变化可能导致观察到的涡流结构发生完全不同。因此，使用实验数据等参考资料对算法的数值参数进行校准是一项复杂的任务。
在过去二十年里，人们已经采用数据驱动的方法对这些问题进行了广泛研究。文献中报道的研究包括不确定性量化（UQ）和不确定性传播（UP）[6]、[7]、数据同化（DA）[8]、[9]、[10]以及机器学习（ML）[11]、[12]的应用。这些方法为改进湍流建模[13]和确定合适的边界条件[14]提供了指导。在所研究的方法中，DA技术显示出广泛的应用潜力。它们与UQ方法相辅相成，并且当用于校准模型的参考数据在空间和时间上分布稀疏时，可以与ML方法高效结合[15]。对于大多数工业应用来说，这尤其成立，因为这些应用中往往无法在关键位置部署传感器。在流体力学研究中，变分DA方法已被广泛应用于稳态/集合平均配置的情况，因为存在解吸引子的特性可以利用这些方法的精度[16]、[17]。对于非稳态、时间分辨的流动，基于贝叶斯理论的集合方法最近也被应用于三维问题[18]、[19]，其中特别包括了基于集合卡尔曼滤波器（EnKF）[20]的方法。这类方法的一个有前景的特点是能够将预测模型与从实际应用中获得的稀疏数据相结合，从而利用数字孪生体来控制物理结构[21]。然而，由于上述计算成本的缘故，这种复杂的应用在流体力学中难以实现。此外，进行DA计算所需的资源通常也比CFD求解器的计算资源更多。因此，目前的应用主要限于周转时间较长的物理流动现象，例如天体流动。
为了在这一领域实现技术进步，人们提出了几种在保持相同准确性的同时减少DA计算资源的策略。包括定位和膨胀在内的策略已被证明是有效的方法[9]。对于像EnKF这样的集合方法，另一种可能性是采用多层次技术。这些方法通常结合少量高精度数值模型的模拟结果和大量低精度运行，以优化准确性、计算成本和鲁棒性。文献中报道的提案中，Popov等人[22]开发了一种多保真度EnKF（MFEnKF），它使用异构模型结果并在控制变量算法中加以整合。Moldovan等人[23]则利用数值求解器的多网格能力获得不同网格分辨率的流动表示。在该算法中，他们将传统的DA分析阶段（称为外循环）与第二次DA更新相结合，其中使用高精度模型的数据来校准低精度求解器中包含的修正项，后者称为内循环。
在本文中，提出了一种使用ML技术改进多层次集合DA方法的新技术。具体来说，利用DA应用的结果来训练低精度模型的修正项，以提高其精度。这个修正项是通过卷积神经网络（CNN）[24]获得的，CNN接收DA算法提供的瞬时和完整的流动快照作为输入。这一策略的一个有趣特点是，由于DA工具能够在广泛的参数空间内生成完整的数据集，因此CNN工具的收敛速度更快且更稳健。这种分析预示了未来的应用场景，即将从物理应用中的有限传感器获取的流数据输入到多保真度DA算法中。这种架构将同时利用即时可用的数据改进数值模型，并通过闭环为物理/数字孪生体架构提供有效的控制策略。本文通过分析雷诺数Re=1000、马赫数Ma=0.5时NACA 0012型翼型的流动特性来研究所提出的算法。在这种情况下，可以观察到流动的非稳态特征，因此需要一个能够考虑流动时间变化的复杂CNN修正项。

文章结构如下：第2节介绍CFD求解器和实验案例的数值细节；第3节介绍数据驱动方法；第4节讨论了先进的EnKF技术和多保真度DA研究；第5节阐述了DA-ML方法的原理以及CNN工具对多保真度算法的改进；最后第6节总结了研究结果并探讨了未来发展方向。

2. 数值工具和实验案例
2.1. 控制方程和数值求解器
本研究使用的数值求解器基于可压缩牛顿流体的有限体积离散化[1]方法来求解纳维-斯托克斯方程。所考虑的方程组如下：
(1) ?ρ?t + ??(ρu) = 0
(2) ?ρu?t + ??ρu?u = ??p + ??μ??u + ??u? ? ??u?δij
(3) ?ρE?t + ??ρE + pu = ??μ??u + ??u? ? ??u?δij?u ? κ?T
其中u表示速度场，p表示压力，ρ表示密度，E表示总能量，T表示温度，μ表示动力粘度，κ表示热导率，δij表示克罗内克符号。方程(1)、(2)和(3)分别代表质量守恒、动量守恒和能量守恒。
数值模拟是通过使用开源工具OpenFOAM对动力方程进行离散化来实现的。该库包含多种基于有限体积（FV）离散化的求解器。不同求解器之间的差异取决于待研究物理流动特性（如湍流、可压缩性效应、热传递等），这些特性通常需要特定的算法来提高效率。本研究中使用的是rhoPimpleFoam求解器，该求解器基于Pimple算法[25]、[26]，适用于非稳态、可压缩流动的分析。OpenFOAM提供了多种时间和空间导数离散化的方案，其中使用了二阶向后时间推进方案。关于空间离散化，对平流项和扩散项都采用了二阶中心格式。

2.2. NACA 0012型翼型绕流的流动
本研究的实验案例是NACA 0012型翼型绕流的二维流动[27]。如图1所示，这种流动具有简单的几何形状，允许使用较少的网格元素进行快速的CFD分析。然而，可以清晰地研究边界层的发展、分离和失速等现象。现在描述该实验案例：翼型沿坐标系的x方向排列，y方向为法线方向。系统原点位于翼型的前缘。数值域的大小为x×y=[?4c,14.5c]×[?0.48c,0.48c]，如图1所示。在入口处设置了一个恒定的速度幅度‖u∞‖，使得雷诺数Re=ρ‖u∞‖c/μ=1000，上游马赫数Ma∞=‖u∞‖/a∞=0.5，其中a∞是声速。a∞定义为a∞=γrT∞，T∞是入口温度，r是特定气体常数，γ是热容比。在这些条件下，流动是层流的，并表现出二维特性。通过改变入口处施加的速度场来控制翼型相对于上游流动的迎角α（见图1）。α的变化会影响流动特性。在域的顶部和底部表面施加周期性边界条件。翼型之间的间距选择为s=0.84c，该值根据文献中的先前研究[27]确定。由于采用了FV离散化方法，OpenFOAM即使在二维计算中也需要三维网格。在这种情况下，网格在展向z方向上只有一个元素，并在侧面施加无通量边界条件。在出口处施加一个非反射的、质量守恒的边界条件。在出口边界施加静态压力狄利克雷条件。OpenFOAM代码的性能通过分析NACA 0012型翼型在开放域中的流动进行了初步评估，结果与文献中报告的结果一致[28]。
下载：高分辨率图像（240KB）
下载：全尺寸图像
图1. 计算域和边界条件
下载：高分辨率图像（420KB）
下载：全尺寸图像
图2. (a) 二维数值域和攻角的变化。(b) OpenFOAM工具snappyHexMesh生成的网格示例

2.3. 网格收敛性和初步模拟
为了评估结果对网格分辨率的敏感性，进行了初步模拟。为此使用了三种网格：参考网格、细网格和粗网格。这些网格是用OpenFOAM内置的snappyHexMesh工具创建的，网格详细信息见表1。图2(b)展示了翼型周围的可视化结果，图3展示了尾流区域内网格分辨率的定性表示。对于参考网格和细网格，在靠近翼型的区域生成了多层元素。如图3(b)所示，细网格在壁法向方向上设置了一层元素，层厚为3?10?3c。图3(c)中显示的参考网格是使用围绕剖面的十层网格创建的。最靠近壁面的层的厚度为7?10^-5c。随着远离剖面，层的厚度以1.5的膨胀比率增加。使用之前介绍的三个网格进行了初步模拟。对于不同的攻角α=0,8,12,20,25°，根据方程(4)改变入口处的速度边界条件，进行了多次运行。使用参考网格进行的模拟采用的时间步长值为δt=10^-5。这个选择允许最大局部Courant数低于2。另一方面，使用细网格和粗网格的模拟采用的时间步长为δt=10^-4。这个选择分别提供了大约3.5和1的最大局部Courant数。所有模拟都是在总物理时间T=25t=4250tA内进行的，其中t代表物理模拟时间，tA=c/||u∞||是特征平流时间。初始条件在整个域内是均匀的，因此ux=u∞,x和uy=u∞,y，见方程(4)。这里需要提醒的是，入口处的速度被设置为Ma=0.5。结果包括升力系数CL、阻力系数CD和斯特劳哈尔数（Strouhal number），在表2中展示。斯特劳哈尔数St的计算公式为[28]：(5)St=f?c||u∞||，其中f是从升力系数CL的时间演变中通过快速傅里叶变换（FFT）分析得到的特征频率。统计是基于[5,25]t时间窗口内的CL和CD信号计算得出的。为了排除初始条件的影响，初始瞬态时期[0,5]t被排除在统计分析之外。模拟结果显示，在α≤20°时观察到稳定流动配置，而在α=25°时观察到非稳定特征。级联配置的 solidity s 负责延迟向非稳定状态的过渡，这在较小攻角时观察到。对于稳定状态，使用精细网格和细网格对时间平均量CLˉ和CDˉ的预测之间的差异很小。然而，将参考网格与粗网格的结果进行比较时，差异最大可达16%。表1列出了用于分析NACA 0012级联配置的计算网格的详细信息。

使用α=25°的攻角进行的模拟显示，尾流特征是涡脱落。图4展示了在这个攻角下的下游流动模式的定性示意图。虽然三个网格使用的尾流看起来相似，但参考网格的St=1.11，细网格的St=1.08，粗网格的St=0.99。通过对升力和阻力系数标准差的时域分析，确认了网格分辨率会影响流动预测。对于升力系数，使用参考网格的模拟得到的值为CL′=0.0075，而细网格得到CL′=0.0071，粗网格得到CL′=0.009。对于阻力系数CD′，参考网格得到CD′=0.0037，细网格得到CD′=0.0026，粗网格得到CD′=0.0039。尽管这些差异明显，但它们的幅度低于表2中观察到的时间平均值的差异。这可以在图5中报告的13个脱落周期内的CL和CD的时间历史中看到。

使用α=25°的攻角进行的模拟表明，对于α≤20°观察到稳定流动配置，而对于α=25°观察到非稳定特征。级联配置的实体度（solidity s）导致了向非稳定状态过渡的延迟，这在单个剖面的较小攻角时观察到。对于稳定状态，使用精细网格和细网格对时间平均量CLˉ和CDˉ的预测之间的差异很小。相比之下，将参考网格与粗网格的结果进行比较时，差异较大，最大达到16%。表1提供了用于分析NACA 0012级联配置的计算网格的详细信息。

图3. 尾流区域网格分辨率的定性表示。(a) 粗网格；(b) 细网格；(c) 参考网格。

使用攻角α=25°进行的模拟显示，尾流特征是涡脱落。图4展示了在这个攻角下的下游流动模式的定性示意图。虽然三个网格使用的尾流看起来相似，但参考网格的St=1.11，细网格的St=1.08，粗网格的St=0.99。通过对升力和阻力系数标准差的时域分析，确认了网格分辨率会影响流动预测。对于升力系数，使用参考网格的模拟得到的值为CL′=0.0075，而细网格得到CL′=0.0071，粗网格得到CL′=0.009。对于阻力系数CD′，参考网格得到CD′=0.0037，细网格得到CD′=0.0026，粗网格得到CD′=0.0039。虽然这些差异明显，但它们的幅度低于表2中观察到的时间平均值的差异。这可以在图5中报告的13个脱落周期内的CL和CD的时间历史中看到。

使用三种网格对五个攻角进行的初步模拟结果。报告了时间平均升力系数CLˉ、时间平均阻力系数CDˉ和斯特劳哈尔数（Strouhal number）。

现在研究了α=25°下的时间平均流场。图6显示了归一化时间平均速度幅度‖uˉ‖/a∞的等值线，流线表示由于边界层分离而形成的再循环区域。虽然参考网格和细网格获得的结果在定性上是相似的，但可以看到使用粗网格得到的流动特征有显著不同。

使用α=25°的攻角进行的模拟显示，尾流特征是涡脱落。图4展示了在这个攻角下的下游流动模式的定性示意图。(a) 粗网格；(b) 细网格；(c) 参考网格。

使用直接指标进行了定量分析，以测量再循环区域的大小。Δxr/c表示沿流向的最大尺寸，而Δyr/c表示法向的最大尺寸。这些量是在满足uxˉ=0的条件下计算的。结果见表3。与参考网格相比，细网格模拟的Δxr/c和Δyr/c的差异分别为8%和4%。最后，使用粗网格得到的流动解在xsep/c处的差异达到27%和7%。总之，数值结果对网格分辨率敏感。这一特性将利用数据同化（Data Assimilation）技术加以利用。内循环通常针对校正项C的优化，通过调整一些可用的自由系数ψ来减少使用粗网格和细网格进行预测时的准确性差异。使用之前为EnKF引入的相同参数，MGEnKF算法的结构包括以下操作：

1. 预测：在粗网格（C）上计算的Ne集合成员使用粗网格模型MC进行时间推进：(13) xi,kCf = MCxCi,k?1 + aθi,k?1 + CxCi,k?1 + aψi,k?1。同样，细网格（F）上的解xf也通过相应的细网格模型MF进行时间推进：(14) xfkf = MFxFk?1。

2. 从细网格到粗网格的投影：使用投影算子Φr?将之前在细网格上获得的预测结果投影到粗网格空间：(15)xC?kf = Φr?xfkf。

3. 内循环：对粗网格模型进行数据辅助（DA）优化。使用从场xC?kf中采样的替代观测值来优化自由系数ψi,k?1。

4. 外循环：在粗网格上进行EnKF计算。通过卡尔曼增益KC使用可用的外部观测值yi,k更新集合成员：(16) Ci,ka = xi,kCf + KCyi,k?HxCi,kf。卡尔曼增益还用于更新投影状态：(17)xC?ka =xCf?k + KCyk?HxC?kf。

5. 细网格状态的更新：将粗网格状态的更新结果使用投影算子Φr投影回全序空间，以更新细网格解：(18) xfka = xf + ΦrxC?ka ?xCf?k。

MGEnKF算法的定性表示如图7所示。

**下载：** 下载高分辨率图像（372KB） / 下载全尺寸图像

图7. MGEnKF外循环的示意图。

**3.3. 具有降阶控制变量的多保真度集合卡尔曼滤波器**
Popov等人[22] 最近提出了一种专门的数据辅助框架，称为多保真度集合卡尔曼滤波器，该框架使用方差缩减策略进行蒙特卡洛方法。该方法结合了全阶物理模型（FOM）和降阶替代模型（ROM）层次结构，以提高数据同化的计算效率。这种方法的一个关键特点是它基于线性控制变量的理论，这是一种在蒙特卡洛方法中常用的方差缩减技术[35]。

考虑一个高保真度模型，其状态由χ表示，以及一个低保真度模型，其状态由γ表示。如果γ与χ有强相关性并且其期望值μγ已知，则γ被称为χ的控制变量。控制变量的主要目的是估计高保真度状态的统计特性，包括从低保真度状态获得的信息，后者在计算上成本较低。为此，控制变量算法使用一个总变量ξ，定义为χ和γ的集合均值μχ和μγ的线性组合。在这些条件下，ξ定义为：(19) ξ = χ ? λγ ? μγ。随机变量χ称为主变量，γ是控制变量，而ξ是总变量[22]。在方程(19)中，ξ的方差比χ小，而其均值等于主变量的均值μξ = μχ。这里选择λ∈R是为了最小化总变量的方差。在大多数实际情况下，μγ是未知的，因此引入了第三个变量η（称为辅助变量）。选择这个变量是为了进行集合计算，使得μγ = μη。使用这种近似，方程(19)变为：(20) ξ = χ ? λγ ? η。

现在将方程(20)转换为一个多级系统，该系统使用细网格-粗网格模型进行模拟，类似于MGEnKF中使用的模型。这两种方法的算法非常相似，但必须定义三个子集合。第一个子集合对应于主变量χ，对其运行Nχ个成员的模拟。这些模拟在细网格上进行。对于控制变量γ，生成Nγ个集合成员的模拟，这些模拟在粗网格上进行。最后，在粗网格上对辅助变量η运行Nη个集合成员的模拟。可以看出，在细网格上进行了Nχ次模拟，在粗网格上进行了Nγ+Nη次模拟。因此，集合成员的总数Ne等于：(21) Ne = 2Nχ + Nη。这意味着可以通过在细网格上进行有限的Nχ次模拟以及在粗网格上进行大量Nγ+Nη次模拟来获得状态估计。现在讨论MFEnKF算法的步骤（见图8中的定性表示）：

1. 预测步骤：根据相应的非线性模型M将三个集合的成员进行时间推进：(22) xi,kχf = Mχxi,k?1 + aχ, xi,kγf = Mγxi,k?1 + aγ, xi,kηf = Mηxi,k?1 + aη。

2. 投影 - 从χ到γ：由于Mχ和Mγ之间的差异，使用主变量和控制变量成员得到的物理场将解耦。使用投影算子Φ?将主变量集合的解投影到控制变量集合：(23) xi,kγf = Φ?xi,kχf。

3. 总变量均值：在该当前的集合框架中，使用三个预测集合变量μ?xχf、μ?xγf、μ?xηf的集合平均值来计算总变量均值估计：(24) μ?kxξf = μ?kxχf ? λμ?kxγf ? μ?kxηf。

4. 总变量观测算子：我们对主变量使用全空间观测算子H，对于通过算子Φ投影的控制变量和辅助变量也适用相同的算子。然而，总变量xξf不是直接观测的。为了克服这个限制，建议使用以下公式的算子Hˉ：(25) Hxi,kξfˉ = Hxi,kχf?1 + 1/2HΦxi,kγf + 1/2HΦxi,kηf。

5. 协方差矩阵的计算：使用异常矩阵来构建假设主变量和辅助变量为独立变量的总变量的近似协方差矩阵：(26) Σxξf,Hxξfˉ = Σxχf,Hxχf + 1/4ΦΣxγf,HΦxγf + ΦΣxηf,HΦxηf?1/2Σxχf,HΦxγf + ΦΣxγf,HxχfΣHxξfˉ,Hxξfˉ = ΣHxχf,Hxχf + 1/4ΣHΦxγf,HΦxγf + ΣHΦxηf,HΦxηf?1/2Σxχf,HΦxγf+ΣHΦxγf,Hxχf。最后，使用独立随机变量构建观测误差协方差矩阵，使得：θχ ～ N0,?χ 和 θγ ～ N0,?η。根据Popov等人[22]，这导致：(28) Σθχ,θχ = Σθγ,θγ = Σθη,θη = 2Σθξ,θξ。使用相应的随机变量获得扰动观测值yχ,yη,yγ。

6. 卡尔曼增益计算：然后计算卡尔曼增益来估计同化后的总变量：(29) K? = Σxξb,HxξbˉΣHxξb,Hxξb + Σθξ,θξ?1xξa = xξb + K?Hxξbˉ?yξ。

7. 变量集合更新：MFEnKF按如下方式同化集合：(30) xi,kχa = xi,kχf + K?yi,kχ?Hxi,kχf, xi,kγa = xi,kγf + ΦK?yi,kγ?Hxi,kγf, xi,kηa = xi,kηf + ΦK?yi,kη?Hxi,kηf。控制变量算法使总变量xξ具有与主变量xχ相同的均值，具有较小的协方差，因此对总变量结果的信心更高。

**下载：** 下载高分辨率图像（437KB） / 下载全尺寸图像

图8. MFEnKF的示意图。

**3.4. 卷积神经网络**
第3.2节已经说明，MGEnKF的内循环可以用来优化一个子网格校正项C，该校正项由一组系数ψ控制，并在粗网格上执行的模拟中集成。Moldovan等人[34]研究了系数ψ的数据辅助优化，表明可以提高在粗网格上进行模型运行的准确性。然而，对于复杂流动，这些函数形式通常难以推导出来。在本研究中，C是通过训练完整的数据驱动模型得到的，因此无需提供函数形式。所研究的数据驱动策略依赖于机器学习（ML），更具体地说，依赖于一种众所周知的神经网络家族。最近在文献中提出了数学框架，以严格整合数据辅助和机器学习方法，目的是模拟动态稀疏模型[36]、推断未解决的尺度[37]或改进湍流建模[15]。

本文介绍了所使用的ML工具。卷积神经网络（CNN）是一类主要应用于计算机视觉领域的机器学习方法[38]。最初由Lecun等人[39]概念化并成功评估，CNN在过去25年中快速发展，应用范围从各种面部和对象识别[40]、[41]到医学成像[42]、[43]，以及地震的地球物理研究[44]。流体力学研究中也应用了CNN，包括实验[45]和数值研究[46]、[47]、[48]。

在处理CNN时，有几种不同的架构可供选择。Zhao等人[24]对主要架构及相关技术进行了全面描述。Ronneberger等人[49]开发了一种全卷积架构，因其草图时的视觉表示而常被称为U-net。这种架构的特点是连续的空间维度缩减（下采样）和数据扩展（上采样）操作。下采样过程（或编码器）在不同尺度提取信息，同时减少空间维度并增加特征深度。上采样过程（或解码器）逐步恢复空间维度，同时增强流动特征的空间定位。CNN的主要结构如图9所示。

如前所述，U-net由一个收缩过程和一个扩展过程组成。一方面，收缩过程包括卷积层，这些层通过修正线性单元（ReLU）函数激活，然后通过最大池化操作进行下采样。在本研究中，这个过程执行了两次。另一方面，扩展过程通过上采样层进行，然后是拼接层，之后应用新的卷积层，这些层也由ReLU函数激活。由于在编码过程中进行了两次收缩操作，解码过程中必须进行相同次数的上采样操作。最后一层是一个通道深度为一的卷积层，目的是获得输出尺寸与输入相同的输出。黑色从左到右的箭头指明了拼接执行的U-net位置。根据图9和表4给出的通道尺寸，这种CNN U-net架构具有592513个可训练参数。最后，这种CNN架构是使用TensorFlow Python库[50]和Keras后端[51]构建的。

**下载：** 下载高分辨率图像（378KB） / 下载全尺寸图像

图9. CNN架构及其对应的通道尺寸。

**表4. CNN U-net结构描述**

| 层 | 输出形状 | 参数数量 |
|---------------|------------|-----------|
| Conv2D | (150,24,1) | 0 |
| Conv2D | (150,24,32) | 118 |
| ReLU | (150,24,32) | 0 |
| Conv2D | (150,24,32) | 368 |
| MaxPool2D | (50,8,32) | 0 |
| Conv2D | (50,8,32) | 368 |
| ReLU | (50,8,32) | 0 |
| Conv2D | (50,8,32) | 368 |
| MaxPool2D | (25,4,32) | 0 |
| Conv2D | (25,4,64) | 737 |
| ReLU | (25,4,64) | 0 |
| UpSampling | (50,8,64) | 0 |
| Concatenate | (50,8,96) | 0 |
| Conv2D | (50,8,32) | 110 |
| ReLU | (50,8,32) | 0 |
| Conv2D | (50,8,32) | 368 |
| UpSampling | (25,4,64) | 0 |
| Concatenate | (50,8,96) | 0 |
| Conv2D | (150,24,32) | 737 |
| ReLU | (150,24,32) | 0 |
| Conv2D | (150,24,32) | 368 |
| Conv2D | (150,24,1) | 115 |

**4. 多级EnKF策略的评估**
第3节中介绍的EnKF模型用于研究NACA 0012轮廓周围的流动。将进行四次数据辅助（DA）运行，并比较结果在准确性和所需计算资源方面的差异。

**4.1. 用于EnKF的设置**
现在介绍用于数据辅助研究的关键元素：
- 观测数据是通过使用参考网格运行的数值模拟的瞬时速度场获得的。该模拟被称为基准真值，选择攻角α=25°进行运行。使用三个传感器采样顺流向速度ux、法向速度uy和压力p，如图10所示。传感器的坐标分别为n1→=[0.032,0.064]、n2→=[0.6,0.12]和n3→=[1.39,0.039]（单位为c）。第一个传感器位于NACA 0012前缘附近。在这个位置附有边界层。第二个传感器位于轮廓的更远位置，根据α的值，它可以观察流动循环或不循环现象。第三个传感器位于尾流区域。观测向量包含六个标量量：(31) y = ux,n1, ux,n2, ux,n3, uy,n1, uy,n2, uy,n3。

通过用加性高斯噪声Ny,4扰动测量值来获得观测的不确定性，使得每个传感器的噪声相同，并且观测值的归一化方差为0.012。数据每ΔDA=13.6tA采样一次，相当于参考模拟的800个时间步长。

- 模型是OpenFOAM原生的非稳态可压缩求解器rhoPimpleFoam，可以使用细网格或粗网格运行。状态向量x是一个大小为2×n的向量，其中n是离散域的单元数，因此使用细网格或粗网格时n的值是不同的。状态向量包含两个速度分量（ux和uy），如下所示：(32) x = ux,0,…ux,i,…ux,n,…uy,0,…uy,i,…uy,n。

现在介绍数据辅助算法。集合方法依赖于生成多个成员，这可能会迅速增加计算成本。此外，通过求解器的停止和重启进行大量EnKF分析阶段并不是一个长期可行的解决方案，从计算效率和人工干预的角度来看。本研究使用自主研发的CONES框架（Coupling OpenFOAM with Numerical EnvironmentS）来执行数据辅助任务，该框架将OpenFOAM的CFD求解器与EnKF数值代码结合，以进行在线数据辅助分析[52], [53]。CONES通过MPI通信并行化数据辅助操作。该库根据用户选择的模拟数量和域划分自动为OpenFOAM分配CPU核心。额外的多个CPU核心被预留用于执行数据同化算法。这种架构通过消除人工操作和减少每次分析阶段之间手动模拟重启所需的机器排队时间，从而提高了整体时间效率。此外，通过自动化数据同化过程，它消除了与离线数据同化分析相关的ROM内存开销，将这些工作转移到了RAM上。库的详细信息见图11，源代码可在https://gitlab.ensam.eu/pe431/cones-dev下载。下载：下载高分辨率图像（104KB）下载：下载全尺寸图像

图10. NACA 0012剖面周围的传感器位置。

将使用不同版本的EnKF进行数据同化运行，以优化进气边界条件的特性。更具体地说，攻角α的值在参考模拟中固定为α=25°，但现在假设它是未知的，并将使用数据同化技术来推断它。这里需要提醒的是，这种优化的目标是减少模型预测与观测之间的差异，同时考虑两种信息来源中规定的不确定性。选择两个连续分析步骤之间的周期ΔDA=13.6tA，以便将更新的进气特性传递给传感器[14]。在这段时间内，流动平流将更新后的α值传递到大约13.6c的长度范围内，这超过了所研究物理域长度的50%。为了获得满意的α系数优化，将修改传播到模拟域中，并消散与进气条件变化相关的不稳定效应是重要的。特别是在最初的分析步骤中，一些集合成员运行的α值可能与目标值相差甚远。最后，本工作中进行的数据同化分析没有使用定位和膨胀程序[9]。对这个测试案例的初步分析表明，这些技术的敏感性非常低，因此为了简化起见，没有使用这些技术。

下载：下载高分辨率图像（519KB）下载：下载全尺寸图像

图11. 在线顺序数据同化的CONES示意图。

4.2. 使用单一网格进行数据同化模型的敏感性分析

前两次数据同化运行使用传统的EnKF进行，所有集合成员都使用相同的CFD求解器和数值网格。它们分别被称为EnKF-CG（模型使用粗网格进行模拟）和EnKF-FG（这里使用细网格）。这两种数据同化应用都依赖于Ne=39个集合成员。这一选择与文献中使用CFD模型的先前数据同化分析一致[54]，可以为状态协方差提供一个合理的近似值，同时也是评估数据同化技术计算成本的参考基准。集合成员初始化为随机值α∈[17°,19°]。集合成员同时在39个Intel(R) Xeon(R) Gold 5218R CPU核心上运行，每个集合成员分配一个CPU核心。此外，EnKF算法本身也在一个相同的CPU核心上执行。图12显示了EnKF-CG和EnKF-FG运行中α的优化值的演变。引入了一个收敛率参数?α，其数学定义为：(33) ?α=‖αk?1?αk‖/‖αk‖。这里αk?1表示前一次EnKF分析的更新后的攻角，αk是当前EnKF迭代优化的攻角。移动平均值定义为?10，平均窗口包括十个连续分析步骤的值。EnKF分析的收敛准则设定为?α<10?4。

首先，两种数据同化程序都收敛到了一个合适的α值。更具体地说，EnKF-FG运行在大约100个分析阶段后获得了α=25.45°的值，而EnKF-CG运行在大约240个分析阶段后优化的攻角收敛到α=25.33°。尽管两次数据同化运行获得的α值并不完全等于目标值α=25°，但它们与数据同化程序选定的传感器测量不确定性是一致的。使用膨胀技术或减少观测的不确定性可以获得更高的精度，但为了简洁起见，这里没有研究这一点。因此，这些结果表明，使用粗网格需要更多的数据同化迭代才能达到收敛。这一点并不奇怪，因为EnKF-FG运行使用的细网格更接近用于产生观测的参考网格的分辨率水平。模型与观测之间的流动快照相似度越高，收敛速度越快。然而，这一有利特性与由于使用了更细的网格而显著增加的计算成本相关。表5提供了执行本节中介绍的EnKF分析所需的计算成本的详细估计。

下载：下载高分辨率图像（484KB）下载：下载全尺寸图像

图12. 上图：应用于进气边界条件的集合平均攻角α的演变。下图：收敛准则?α的演变。

虽然EnKF-CG运行的平均成本为0.42 CPU小时来完成一个完整的EnKF预测和分析周期，但需要240次迭代才能实现最优参数的正确收敛，从而导致总计算成本接近100 CPU小时。另一方面，EnKF-FG运行平均需要10.42 CPU小时的计算负载来完成一个预测和分析周期，而需要86个周期才能达到收敛，总计算时间为900 CPU小时。可以得出结论，使用细网格进行EnKF运行所需资源大约是使用粗网格的九倍。

表5. 执行数据同化运行所需的计算成本

数据同化运行次数
CPU小时数
EnKF-CG：240 0.42
EnKF-FG：86 10.42
MGEnKF：33 2.73
MFEnKF：34 1.05
ConvMGEnKF：10 2.73
MFEnKF：11 2.33
MFEnKF：12 3.13

现在研究了数据同化运行预测的流动物理特性。图13显示了瞬时归一化速度幅度‖uˉ‖/a∞场的等值线，并与时间t=25时的参考模拟进行了比较。即使两种EnKF运行都得到了令人满意的优化攻角α值，但EnKF-CG预测的瞬时速度场与参考模拟相比在尾流区域存在明显差异。观察图14中呈现的时间平均速度幅度剖面也可以得出类似的结论。这种分辨率不足是导致体积流量量预测出现误差的原因，如表6所示。所呈现的结果包括升力系数CL、阻力系数CD、Strouhal数St以及用于测量再循环区域大小的参数Δxr/c、Δyr/c和xsep/c的时间平均值和标准差。与第2.3节中使用的α=25°的CFD运行结果相比，EnKF-CG运行的体积流量量结果与使用粗网格的模拟结果相同。这是预期之中的，因为数据同化通过的数据更新每13.6tA进行一次，因此更新后的解会周期性地向下游传输而不影响统计数据。考虑到EnKF-FG的情况，可以在图13和图14中看到流动拓扑被很好地表示出来，并且与参考模拟在质量上非常相似。对于这一情况，表6再次显示体积流量量几乎与第2.3节中使用细网格的α=25°的CFD模拟提供的结果相同。然而，在这两种情况下，与参考模拟的结果之间的差异明显较小。尽管流动快照之间的相似性不强，但数据同化算法仍然足够准确，能够识别出一个合适的攻角。然而，用于EnKF-FG和EnKF-CG运行的数据同化实现无法改善底层模型的预测。可以说，如果能够在更短的时间窗口内获得观测数据，那么来自粗网格和细网格上的数据同化过程的状态估计可以提高流动的统计量[55]、[56]。然而，实际系统中的采样频率有限，其特征时间显著大于数值模拟中通常的时间步长[57]。另一种可能性是包括子网格校正项以提高模型的准确性。这项任务将在第5节中明确讨论。

下载：下载高分辨率图像（226KB）下载：下载全尺寸图像

图13. 在t=25时计算的瞬时归一化速度幅度‖uˉ‖/a∞的等值线，(a) EnKF-CG运行（α=25.33°）；(b) EnKF-FG运行（α=25.45°）；(c) 参考模拟（α=25°）。

总之，通过使用不同网格进行的两次EnKF敏感性分析表明，该算法能够系统地优化影响测试案例的自由参数。然而，当考虑优化的收敛速率时，数据同化过程对模型网格的选择非常敏感。使用细网格虽然提高了性能，但却增加了计算开销。另一方面，尽管收敛性有所下降，但在粗网格上进行的分析仍然能够得到可接受的结果，同时保持了较低的计算努力。在下一节中，将在多网格和多保真度分析的框架下，结合使用细网格和粗网格的集合实现，采用之前介绍的MGEnKF和MFEnKF技术。

表6. 时间平均升力系数CLˉ、阻力系数CDˉ、Strouhal数St以及再循环区域Δxr/c、Δyr/c、xsep/c的特性。数据提供了用于数据同化运行的结果，以及用于获得观测的参考模拟的结果。

空单元格
α[°]
CLˉ：
CL′：
CDˉ：
CD′：
St：
Δxr/c：
Δyr/c：
xsep/c：

EnKF-CG：
25.33：0.39
0.01
0.35
0.00
4：
10.79
0.08
6：
0.42

EnKF-FG：
25.45：
0.37
0.00
8：
0.30
0.00
3：
1.07
0.98
0.07
8：
0.2

MGEnKF：
25.37：
0.37
0.00
8：
0.31
0.00
3：
1.08
0.98
0.07
9：
0.2

ConvMGEnKF：
25.09：
0.37
0.00
7：
0.30
0.00
3：
1.09
0.97
0.07
7：
0.2

ConvMFEnKF：
25.24：
0.37
0.00
8：
0.31
0.00
3：
1.09
0.98
0.07
8：
0.2

ConvMFEnKF：
24.97：
0.37
0.00
7：
0.30
0.00
3：
1.09
0.98
0.07：
8：
0.2

地面真实值：
25：
0.37
0.00
8：
0.30
0.00
4：
1.11
1.04
0.08
0.16

下载：下载高分辨率图像（593KB）下载：下载全尺寸图像

图14. 在t=25时计算的瞬时归一化速度幅度‖uˉ‖/a∞的等值线，(a) EnKF-CG运行（α=25.33°）；(b) EnKF-FG运行（α=25.45°）；(c) 参考模拟（α=25°）。

4.3. 使用多级技术进行数据同化应用

现在使用MGEnKF和MFEnKF技术进行了两次额外的数据同化运行。为了与EnKF-FG和EnKF-CG的运行进行适当的比较，使用了与前一项研究相同的机器，共进行了39次集合成员的模拟。具体来说，对于MGEnKF，在粗网格上模拟了38个集合成员，在细网格上模拟了一个成员；对于MFEnKF，使用细网格运行了两个主要成员，在粗网格上运行了两个控制成员以及35个辅助成员。数据同化运行的其他参数，包括连续分析阶段之间的时间窗口和入口条件的参数α的初始化，都遵循了EnKF-FG和EnKF-CG运行使用的相同标准。此外，方程(13)中的C项设置为零，即不执行内循环。

图15显示了参数α的数据同化优化的演变，其中系数?α用于衡量参数推断的收敛速率（见方程(33)）。报告了MGEnKF和MFEnKF运行的结果。与之前分析的两次数据同化运行类似，获得了正确的攻角值，MGEnKF为α=25.37°，MFEnKF为α=25.09°。这两种方法分别需要332（MGEnKF）和345（MFEnKF）个数据同化分析阶段才能达到收敛。

现在将多级方法所需的计算成本与EnKF-FG和EnKF-CG运行的结果进行比较。表5提供了总结。可以得出的第一个信息是，MGEnKF和MFEnKF需要更多的分析阶段才能达到参数α所需的收敛水平。这一观察可能与需要结合来自不同数值模型的信息有关，这需要额外的迭代来进行平滑和平衡。然而，在DA程序的预测和分析周期的成本方面，MGEnKF的平均计算时间为1.05 CPU小时，而MFEnKF非常接近，为1.09 CPU小时。这些要求大约是EnKF-CG运行计算时间的两到三倍，但仍然比EnKF-CG运行提供的估计值小十倍。如果考虑执行参数α优化的总CPU成本，那么EnKF-CG是四种DA运行中最好的。然而，需要提醒的是，这是唯一一个完全依赖于粗网格使用的应用，并且对整体流量量的全局预测并不令人满意。接下来，将研究依赖于粗网格和细网格的MGEnKF和MFEnKF运行的流量特征。下载：下载高分辨率图像（507KB）下载：下载全尺寸图像图15. 上图：应用于入口边界条件的平均攻角α的演变。下图：收敛标准?α的演变。图16显示了参考模拟与MGEnKF和MFEnKF运行之间的归一化瞬时速度场‖u‖/a∞的比较。与EnKF-FG的情况类似，瞬时流场与参考解极其相似。通过观察图17中显示的时间平均归一化速度大小，可以得出相同的结论。表6中报告的对整体流量量的分析也确认，MGEnKF和MFEnKF运行的统计量与使用细网格的α=25°的CFD实现是等效的。因此，现在提供了更多关于DA方法预测准确性的结果，可以再次分析计算成本的问题。如果需要准确预测流量的整体量，那么多级方法显然优于仅基于一个网格的两种策略。MGEnKF和MFEnKF的运行比EnKF-CG的实现更精确，而且它们的准确性与EnKF-FG相同，但它们的总计算成本仅为其两到三倍。总之，先进的多级方法在高效结合高分辨率和低分辨率来源的信息方面比标准EnKF具有优势。然而，一个潜在的危险是，细网格和粗网格得到的流量表示之间的差异可能是不可调和的，导致滤波器发散。对于MGEnKF和MFEnKF来说，例如在预测步骤中，如果在粗网格上模拟的成员的预测与在细网格上运行的成员的预测表现出非常不同的演变，由于子网格误差的累积效应，这种情况可能发生。这一考虑表明，在保持较低计算成本的同时，获得足够的低保真度模型的准确性至关重要。在下一节中，将研究机器学习工具是否可以提高在粗网格上模拟的集合成员的低保真度解的准确性，从而提高算法的性能。下载：下载高分辨率图像（228KB）下载：下载全尺寸图像图16. 在t=25时，（a）在细网格上优化的MGEnKF角度α=25.37°；（b）在细网格上优化的MFEnKF角度α=25.09°；（c）参考模拟（α=25°）的瞬时归一化速度大小‖u‖/a∞。下载：下载高分辨率图像（591KB）下载：下载全尺寸图像图17. 在t=25时计算的时间平均归一化速度大小‖uˉ‖/a∞的等值线，对于（a）MGEnKF运行（α=25.37°）；（b）MFEnKF运行（α=25.09°）；（c）参考模拟（α=25°）。5. 基于CNN增强的多级策略在本节中，使用机器学习工具来减少基于粗网格和细网格的模型对流量状态预测的差异。Moldovan等人在MGEnKF框架内首次尝试了这项任务[34]。在内部循环中，为粗网格上运行的模型方程（13）中的C项提出了一个函数形式，并使用来自细网格模拟的数据流对其自由系数进行了优化。虽然该程序是成功的，但主要缺点是无法为复杂流动情况数学上推导出合适的子网格修正项C的函数形式。因此，随着获得此类项的有效描述所需的自由参数数量的增加，计算成本可能迅速变得难以承受。为此，这里使用深度学习工具代替函数形式来表示项C。这些数据密集型算法由于通过EnKF生成了完整的流量快照而更容易训练。这里使用的DA和ML方法的互补特性通过以下方案结合：•首先进行一次DA运行，如第4.3节中介绍的那样。数据在多个时间点存储，尤其是在最初的DA分析期间，其中不同集合成员的参数α的值被分散。•使用这些数据来训练一个卷积神经网络（CNN）模型，该模型模仿项C。然后将该黑盒工具集成到在粗网格上执行的数值模拟中。因此，在时间步骤k时，用于项C的CNN工具必须提供一个模型修正：(34)CkCNNxkC=ΦMxkF?MrxkCin，以便获得ΦxkF=xkC，即在使用细网格和粗网格时最小化流量预测的差异（见方程（13）–（18））。这里需要提醒的是，Φ是从细网格到粗网格的投影算子。通过ML工具进行的优化只是近似的，因为它应用于CFD求解器在粗网格上预测的最终状态。此操作不包括在求解器的递归步骤中，因此不一定符合方程的守恒性。然而，由于与DA算法一起进行的状态更新，已经验证状态修正很小，与不守恒相关的效应在几次迭代后迅速消散，类似于初始条件效应。选择直接进行状态修正是因为实现的简单性，这与CNN算法处理图像（状态）的操作特性有关。然而，可以设想更复杂的架构，其中ML强制项可以包含在动态方程中[15]，这些将在结论中讨论。通过ML工具进行的状态修正是通过为每个流量物理量生成几个不同的CNN模型来获得的。这种架构在CONES库中得到了形式化和实现，遵循图18中呈现的方案。这种方法利用CNN的特征提取能力来解决由子网格误差引起的差异，从而提高了使用粗网格的计算成本较低的模拟的预测能力。下载：下载高分辨率图像（505KB）下载：下载全尺寸图像图18. CONES的示意图，用于在线顺序数据同化，并集成预训练的CNN工具。5.1. CNN训练现在描述CNN模型的训练过程。在机器学习中，数据质量对模型性能至关重要。次优数据可能导致训练效率显著下降，进而影响结果CNN的预测能力。用于训练的数据集包括32000对输入和输出场实例。具体来说，输入数据对应于从粗网格模拟中获得的流量场，而根据方程（34），输出数据是从细网格投影的解决方案与输入解决方案之间的计算差异。图19给出了视觉表示。上述32000对输入/输出是从MFEnKF运行中生成的。该模拟采用了特别配置的架构，包括10个主要成员、10个控制成员和10个辅助成员。这种设计选择非常重要，因为它影响了生成数据的多样性和代表性，直接影响训练模型的泛化能力。这些数据生成的初始攻角随机设置在17到19度之间，旨在优化到接近25度的值，如先前研究中所述。从DA周期中获取数据的过程严格地在四个不同的预测/分析阶段进行，每个离散时间步骤都进行了测量。这种全面的采样策略共产生了4×800×10=32000个快照，确保了足够大且多样化的数据集，以进行稳健的模型训练。需要注意的是，每对输入数据完全来自控制成员数据，为模型的预测任务提供了一致的基线。相比之下，输出CCNN是从主要成员和控制成员中计算得出的量。还可以看到，在四次DA分析期间，训练程序可以访问不同攻角α的数据，从而丰富了提供给ML算法的信息。下载：下载高分辨率图像（349KB）下载：下载全尺寸图像图19. CNN架构及相应的输入和输出数据。所展示的案例研究使用了适合身体的网格，这意味着从OpenFOAM中提取的流量场不是适合CNN的格式。CNN和神经网络通常执行矩阵运算，而非结构化网格本身不提供这些运算。为了解决这个问题，使用了网格投影策略。将CFD求解器中的流量场投影到一个2D结构化网格上。这个结构化网格与计算域的维度对齐。具体来说，原本由2600个单元格组成的粗网格网格被映射到一个150×24的2D结构化网格上，从而产生了3600个结构化单元格。这种网格投影是使用Python Scipy库中的griddata函数执行的。训练了三种不同的模型，分为两个功能家族：1. 速度模型（x分量）：该模型用于预测项CuxCNN。2. 速度模型（y分量）：该模型用于预测项CuyCNN。3. 压力模型：该模型用于预测项CpCNN。这三个模型按照状态变量顺序CkCNN的顺序排列：(35)CkCNN=CuxCNN,CuyCNN,CpCNN，并用于同步预测的流量状态，使得方程（13）和方程（22）现在具有以下形式：(36)xCkf=Mxk?1f,θk?1+CkCNNxCk?1f (37)xkγ=Mγxk?1f,θk?1+CkCNNxk?1χ方程（36）、（37）使用的模型是使用Adam优化器[58]编译的，而方程（35）中的压力模型使用RMSProp优化器[59]。RMSProp和Adam都是深度学习工作中广泛使用的自适应学习率优化器。它们的实用性在于它们能够自动调整训练期间各个网络参数（权重和偏置）的学习率，与全局固定的学习率相反。这种自适应特性显著提高了它们在复杂损失景观中的效率。为了减轻过拟合并增强泛化能力，速度模型结合了 dropout层[60]，dropout率为10%，并战略性地放置在激活层之后。之所以选择直接进行状态修正，是因为实现的简单性，这与CNN算法处理图像（状态）的操作特性相关。然而，可以设想更复杂的架构，其中ML强制项可以包含在动态方程中[15]，这些将在结论中讨论。通过ML工具进行的状态修正是通过生成几个不同的CNN模型来获得的，每个流量物理量一个模型。这种架构在CONES库中得到了形式化和实现，遵循图18中呈现的方案。这种方法利用CNN的特征提取能力来解决由子网格误差引起的差异，从而提高了使用粗网格的计算成本较低的模拟的预测能力。下载：下载高分辨率图像（505KB）下载：下载全尺寸图像图18. CONES的示意图，用于在线顺序数据同化，并集成预训练的CNN工具。5.1. CNN训练现在描述CNN模型的训练过程。在机器学习中，数据质量对模型性能至关重要。次优数据可能导致训练效率显著下降，从而影响最终CNN的预测能力。用于训练的数据集包括32000对输入和输出场的实例。具体来说，输入数据对应于从粗网格模拟中获得的流量场，而输出数据根据方程（34），是从细网格投影的解决方案与输入解决方案之间的计算差异。图19给出了视觉表示。上述32000对输入/输出是从MFEnKF运行中生成的。该模拟采用了由10个主要成员、10个控制成员和10个辅助成员组成的特别配置的架构。这种设计选择非常重要，因为它影响了生成数据的多样性和代表性，直接影响训练模型的泛化能力。这些数据生成的初始攻角随机设置在17到19度之间，旨在优化到接近25度的值，如先前研究中所描述的。从DA周期中获取数据的过程在四个不同的预测/分析阶段严格执行，每个离散时间步骤都进行了测量。这种全面的采样策略产生了总共4×800×10=32000个快照，确保了足够大且多样的数据集，以进行稳健的模型训练。需要注意的是，每对的输入专门来自控制成员数据，为模型的预测任务提供了一致的基线。相比之下，输出CCNN是从主要成员和控制成员中计算得出的量。还可以看到，在四次DA分析期间，训练程序可以访问不同攻角α的数据，从而丰富了提供给ML算法的信息。下载：下载高分辨率图像（349KB）下载：下载全尺寸图像图19. CNN架构及其相应的输出数据。所展示的案例研究使用了适合身体的网格，这意味着从OpenFOAM中提取的流量场不是适合CNN的格式。CNN和神经网络通常执行矩阵运算，而非结构化网格本身不提供这些运算。为了解决这个问题，使用了网格投影策略。将CFD求解器中的流量场投影到一个2D结构化网格上。这个结构化网格与计算域的维度对齐。具体来说，最初由2600个单元格组成的粗网格网格被映射到一个150×24的2D结构化网格上，从而产生了3600个结构化单元格。这种网格投影是使用Python Scipy库中的griddata函数执行的。训练了三种不同的模型，分为两个功能家族：1. 速度模型（x分量）：该模型用于预测项CuxCNN。2. 速度模型（y分量）：该模型用于预测项CuyCNN。3. 压力模型：该模型用于预测项CpCNN。这三个模型按照状态变量顺序CkCNN的顺序排列：(35)CkCNN=CuxCNN,CuyCNN,CpCNN，并用于同步预测的流量状态，使得方程（13）和方程（22）现在具有以下形式：(36)xCkf=Mxk?1f,θk?1+CkCNNxCk?1f (37)xkγ=Mγxk?1f,θk?1+CkCNNxk?1χ方程（36）、（37）使用的模型是使用Adam优化器[58]编译的，而方程（35）中的压力模型使用RMSProp优化器[59]。RMSProp和Adam都是在深度学习工作中广泛使用的自适应学习率优化器。它们的实用性在于它们能够在训练期间自主调整各个网络参数（权重和偏置）的学习率，与全局固定的学习率相反。这种自适应特性显著提高了它们在处理复杂损失景观时的效率。为了减轻过拟合并增强泛化能力，速度模型结合了dropout层[60]，dropout率为10%，并战略性地放置在激活层之后。模型泛化进一步使用k折交叉验证方法[61]进行评估，k=4折。数据集被划分为75%用于训练，25%用于测试。此外，从主要数据集中单独预留了一组50个输入/输出对的随机批次用于训练后评估。所有模型都使用均方误差（MSE）作为损失函数进行优化，并在400个周期内进行训练。训练计算过程使用了Nvidia Rtx A6000 GPU，并使用Tensorflow和Keras进行Python并行处理。CNN训练的具体计算成本在表7中列出。首先在CPU上执行的预处理步骤对于准备CNN的训练数据至关重要。此阶段包括将OpenFOAM场从非结构化网格插值到结构化网格等操作。这一步需要0.3 CPU小时。随后，由于神经网络操作的高度并行性，机器学习工具的训练阶段主要由GPU承担。速度模型的训练需要0.42 GPU小时，同时压力模型的训练需要0.63 GPU小时。表7. CNN训练过程的计算成本。模型预处理（CPU小时）训练（GPU小时）速度模型0.30.42压力模型0.30.635.2. 带有CNN模型的EnKF运行现在描述了之前介绍的与DA框架一起引入的ML工具。特别是，将进行两次新的DA运行。它们展示了与MGEnKF和MFEnKF相同的运行设置，但除此之外，还包含了CNN工具，以提高在粗网格上运行的集合成员预测的准确性。这些新的运行分别被称为ConvMGENKF和ConvMFEnKF，其程序的示意图如图20所示。现在提供了关于CNN工具使用的更多细节。CONES中使用的并行计算架构是多程序多数据（MPMD）[62]。MPMD结构非常适合需要不同算法（CFD求解器、数据同化代码和CNN模型）处理计算任务的数值问题，这些算法需要相互作用并交换数据。通过使用MPI通信例程，将预训练的CNN模型直接集成到CONES中。CNN模型被设置为每53个时间步长δt对低保真模型提供的流场应用状态修正CkCNN。连续应用CNN工具之间的时间窗口被选择为等于参考仿真计算的斯特劳哈尔数Stref=1.11的平均时间。考虑到数据同化分析每800个时间步长进行一次，这意味着在连续应用EnKF之间大约进行了15次状态修正。敏感性分析表明，状态更新的次数减少会逐渐降低解决方案的准确性，最终导致使用粗网格的模型的原始预测。另一方面，更频繁的状态更新会提高成本，但并不会显著提高准确性。CNN工具的应用仅限于特定的空间子域，其大小为x/c∈[?2.5,10]，y/c∈[?0.24,0.24]，如图21所示。这种选择是为了考虑在状态预测中观察到较大差异的区域，该区域对应于机翼区域。这样做是为了降低计算成本并避免与边界条件相互作用可能引起的虚假效应。

下载：下载高分辨率图像（786KB）
下载：下载全尺寸图像
图20. 使用CNN工具提高在粗网格上运行的模型预测能力的DA运行（a）ConvMGEnKF和（b）ConvMFEnKF的示意图。

现在将展示并讨论这两种DA运行的结果。
下载：下载高分辨率图像（209KB）
下载：下载全尺寸图像
图21. 在粗网格上进行的CFD运行中应用CNN工具的区域。该区域用绿色标出以便于查看。

图22显示，增强后的DA运行ConvMGEnKF和ConvMFEnKF能够获得与MGEnKF和MFEnKF非常相似的参数α的精确估计值。然而，这两种增强后的DA运行在所需分析阶段方面显著加快了优化速度。如图22（a）所示，根据之前提出的收敛标准?α<10?4，ConvMGEnKF运行在102次DA分析后获得了优化后的攻角α=25.24°，而MGEnKF运行则需要300次迭代。然而，这种算法收敛速度的加速与运行CNN工具所需的计算资源增加有关。表5中报告了总结。可以看出，单个预测-分析DA步骤的成本对于ConvMGEnKF运行来说是2.73 CPUh。这个成本高于MGEnKF运行所需的成本（1.05 CPUh），但低于EnKF-FG运行所需的成本（10.42 CPUh）。尽管与MGEnKF运行相比每个周期的成本增加了，但由于收敛速度的大幅加快，总体计算成本显著降低。考虑到优化过程需要102步才能收敛，因此该运行的总计算成本约为278.46 CPUh。因此，增强了CNN的DA运行显示出与仅使用细网格的EnKF运行相似的收敛速率，但其总体计算成本明显更低。对于这个应用，MGEnKF的总体计算成本减少了近20%。在这种情况下，增加的每个周期成本与减少的总迭代次数之间的权衡显然是有利的。通过比较图22（b）中的ConvMFEnKF和MFEnKF运行，可以得出类似的结论。在这种情况下，CNN工具的集成使得在112次DA迭代内获得了优化后的α=24.7°的值。这种优化过程的快速收敛突显了机器学习模型在多保真度框架中的有效性。与ConvMGEnKF相比，这种情况下的每个周期的计算开销略有增加。在这种情况下，总体计算成本相对减少了约7%。

下载：下载高分辨率图像（502KB）
下载：下载全尺寸图像
图22. 顶部：进气边界条件攻角α的演变。底部：收敛标准?α的演变。（a）MGEnKF与ConvMGEnKF的比较；（b）MFEnKF与ConvMFEnKF的比较。

最后，表6中展示了总体流量量的分析结果。结果确认，ConvMGEnKF和ConvMFEnKF的运行对于每个调查的量提供了与它们对应的MGEnKF和MFEnKF相当的预测，但计算成本降低，优化效果更高。通过对图23中显示的时间平均值速度大小的预测流场进行分析，也可以得出类似的结论。CNN增强型DA运行的结果与参考仿真数据显示出非常相似的等值线。因此，当前结果强调了使用CNN工具增强CFD模拟和表示子网格效应的潜力。

下载：下载高分辨率图像（663KB）
下载：下载全尺寸图像
图23. 在t=25时计算的（a）ConvMGEnKF运行（α=25.24°）；（b）ConvMFEnKF运行（α=24.97°）；（c）参考仿真（α=25°）的时间平均值速度大小‖uˉ‖/a∞的等值线。

6. 结论

本研究提出了一种改进多层次和多保真度集合数据同化技术的新方法。为此，使用了来自机器学习领域的卷积神经网络（CNN）工具来学习校正项，以补偿由于使用粗网格而产生的子网格误差。这是通过结合不同保真度级别的仿真数据来完成的，这些数据的预测被提供给ML工具进行训练。所得工具旨在减少高保真度和低保真度实现之间的差异，同时只需要适度的计算资源增加。这些技术已经在Re=1000和Ma=0.5的NACA 0012机翼级联的仿真中得到了验证。这个测试案例展示了由机翼入射角α控制的非稳态特性。数值仿真使用了OpenFOAM库中的rhoPimpleFoam求解器进行。数据同化过程的目标是优化被认为是未知的参数α。使用来自机翼周围三个传感器的观测数据流，从α=25°的参考仿真中对该参数进行了校准。进行了多次DA运行，使用了在细网格和粗网格上的模型运行，以及多层次方法。对于后者，还进行了包括CNN状态修正的调查研究。观察到后者在计算资源需求上总体减少了20%，同时提供了准确的结果。此外，所有DA运行都展示了与使用最细网格的EnKF程序相比的流量量和流场时间平均特性的预测准确性。

未来将对此主题进行进一步研究，以提高DA + ML工具的准确性和多功能性。首先，本工作中使用的CNN工具是通过初步的DA调查进行训练的，然后作为黑盒模型应用于后续运行。其中一个目标是在DA过程中进行在线训练，而不是在此之前。虽然这在算法收敛性方面更具挑战性，但计算成本也较低。未来的另一个研究点涉及CNN工具的特性以及它们如何集成到数值CFD求解器中。有人评论说，目前的校正是在求解器预测之上叠加的，因此，如果校正很重要或在每个时间步骤都应用，可能会对动态方程产生不稳定性。可以设想训练CNN工具作为方程中的强迫项，以提供更大的算法稳定性。然而，这里的一个挑战与非线性动态方程的递归性质有关，这将在每个时间步骤中需要多次评估CNN模型。最后一个挑战是DA + ML架构的准确性。观察到该方案的准确性与在最细网格上运行的仿真相同。讨论认为，模型的缺乏优化以及连续DA分析之间的时间窗口是造成这一结果的原因。未来的研究可以考虑训练第二个CNN工具，以改善使用细网格的模型的预测。在这种情况下，目标是由EnKF产生的状态估计。

所有这些关于未来活动的讨论点突显了潜力，同时也强调了探索数据同化和机器学习结合方法的必要性。目前正在进行的研究是关于通过压缩机的三维流动。数据同化方法将数值仿真与实验结果相结合，以训练神经网络，目的是通过有效改进湍流建模来增强数值解。这种新的分析还将克服当前工作的一个主要限制，即DA和ML方法的耦合目前仅限于一个相对简单的测试案例。这些研究对于预见工业应用（如运输工程、城市流动或能量收集）的应用是必要的。

CRediT作者贡献声明：
Tom Moussie：写作 - 审查与编辑、原始草稿编写、可视化、软件、方法论、调查。
Paolo Errante：写作 - 审查与编辑、验证、监督。
Marcello Meldi：写作 - 审查与编辑、验证、监督、项目管理、方法论、概念化。