利用尺度分离和随机闭合方法进行混沌动力学的数据驱动预测

《Data-Centric Engineering》：Leveraging scale separation and stochastic closure for data-driven prediction of chaotic dynamics

　　摘要：模拟湍流流动是一项计算上极具挑战性的任务，因为它需要解析微小尺度结构并在多个尺度上捕捉复杂的非线性相互作用。因此，大量的研究都集中在从数据驱动的角度分析湍流流动。然而，由于这些系统的复杂和混沌性质，传统模型往往变得不稳定。为了克服这些限制，我们提出了一种纯粹的随机方法，该方法分别处理大尺度相干结构的演变和高保真统计数据的闭合。为此，使用自回归模型来学习过滤后数据的动态。这种方法结合了变分自编码器（VAE）和Transformer架构。VAE的投影是概率性的，确保了模型的随机性与流动的统计特性之间的一致性。我们模型中随机采样轨迹的平均实现显示出相对的$L_1$和$L_2$距离分别为6%和10%。此外，我们的框架能够构建有意义的置信区间，实现了80%的预测区间覆盖概率，并且区间宽度最小。为了从过滤后的空间中恢复高保真速度场，采用了高斯过程（GP）回归。这种策略已经在表现出混沌行为的Kolmogorov流动的背景下进行了测试。我们将我们模型的性能与最先进的概率基线进行了比较，包括VAE和扩散模型。我们证明了基于高斯过程的闭合方法在该特定测试床中捕获一阶和二阶统计特性的性能优于这些基线，提供了稳健且自适应的置信区间。

影响声明：在流体动力学中模拟湍流是出了名的困难。数据驱动的方法可以降低计算成本，但往往无法捕捉到复杂的时空依赖性，导致即使在短期预测中轨迹也会迅速发散。我们提出学习流动的潜在统计结构而不是确定性轨迹。通过首先解析大尺度结构动态，我们降低了空间复杂性，并使用自回归降阶模型（ROM）实现了更稳定的长期预测。最后，我们引入了一种基于高斯过程的闭合策略来恢复小尺度结构。

1. 引言：在模拟湍流时，确定性工具箱遇到了两个显著的挑战。第一个挑战与动态系统的混沌和不可预测性有关，这种系统表现出看似没有重复模式的混乱轨迹，并且对初始条件非常敏感。正如Lorenz所说：“现在决定了未来，但近似的现在并不能近似决定未来。”第二个障碍源于问题的复杂性和多尺度性质。涉及的各个尺度对直接数值模拟（DNS）来说是一个难题，因为需要过度的分辨率约束来考虑所有尺度；对于机器学习模型来说，它们也难以表示这种内在的复杂性。从降阶建模的角度来看（ML策略是这一类别的一个子集），我们认为将任务分解为两个连续的子任务将克服上述两个障碍。“第一个任务”是从低通过滤数据中推导出动态，重点关注对总动能有显著贡献的大结构的演变。这通常是流体力学界采用的方法，依赖于大涡模拟（LES）（参见Deardorff，参考文献Deardorff1970），而不是直接数值模拟（DNS）。然而，与LES一样，闭合问题仍然存在。因此，在“第二个任务”中，一旦在过滤空间中重构了轨迹，我们就假设存在一个从该空间到全分辨率域的时间独立且随机的映射。第一个任务涉及构建一个自回归和预测性的降阶模型（ROM）。ROM框架已经在包括混沌系统在内的各种动态系统中展示了有希望的能力（Lassila等人，参考文献Lassila, Manzoni, Quarteroni, Rozza, Quarteroni和Rozza2014）。例如，适当的正交分解（POD）使得系统状态能够投影到由系统的时空相干结构（称为“模式”）所跨越的更低维空间。作为一种分析推导出的技术，POD揭示了可解释的模式，并被广泛采纳为物理和湍流建模中降维的强大基线（Schmid，参考文献Schmid和Durbin2021）。在非线性动力学领域，人们对Koopman理论越来越感兴趣，该理论通过无限维的Koopman算子支持即使对于强非线性系统也存在线性动力学（Brunton等人，参考文献Brunton, Budi?i?, Kaiser和Kutz2022）。许多工作和算法专注于找到这种算子的可行、有限维的近似，以捕捉本质的动力学。其中，动态模态分解（DMD）是最突出的例子之一（Schmid，参考文献Schmid2010）。另一种强大的方法涉及识别由线性化系统在固定点附近的特征空间所跨越的最慢不变流形，并通过系统的非线性将其扩展。这种方法称为光谱子流形（SSMs），揭示了允许模型简化和动态简化的低维非线性结构（Haller和Ponsioen，参考文献Haller和Ponsioen2016）。或者，这些低维结构也可以使用机器学习（ML）来识别，从而放宽了POD固有的线性或正交性约束。这些ROM方法通常可以分为两大类：侵入式的，它们通过例如Galerkin投影（Menier等人，参考文献Menier, Bucci, Yagoubi, Mathelin和Schoenauer2023）针对Navier–Stokes方程进行定制；和非侵入式的，它们完全是数据驱动的（Gupta等人，参考文献Gupta, Schmid, Sipp, Sayadi和Rigas2023；Simpson等人，参考文献Simpson, Vlachas, Garland, Dervilis和Chatzi2024；Menier等人，参考文献Menier, Kaltenbach, Yagoubi, Schoenauer和Koumoutsakos2025；?zalp等人，参考文献?zalp, Nóvoa和Magri2025b）。一些非侵入式方法甚至以与网格无关的方式运行（Serrano等人，参考文献Serrano, Migus, Yin, Mazari和Gallinari2023）。嵌入式存储架构的最新进展，特别是注意力机制和Transformer（Vaswani等人，参考文献Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser, Polosukhin, Guyon, Luxburg, Bengio, Wallach, Fergus, Vishwanathan和Garnett2017）在时间序列应用中的发展，极大地改变了这一领域（Shokar等人，参考文献Shokar, Kerswell和Haynes2024；Solera-Rico等人，参考文献Solera-Rico, Sanmiguel Vila, Gómez-López, Wang, Almashjary, Dawson和Vinuesa2024；Zighed等人，参考文献Zighed, Nicolas, Gallinari和Sayadi2025）。无论采用哪种策略，ROM通常依赖于某种形式的降维，如适当的正交分解（POD）、自编码器或其他压缩技术。这是基于一个假设，即系统的动态在相空间中的低维流形上演变（Haller和Ponsioen，参考文献Haller和Ponsioen2016；Buza等人，参考文献Buza, Jain和Haller2021）。值得注意的是，这一假设通常适用于混沌系统（Liu等人，参考文献Liu, Axas和Haller2024），臭名昭著的Lorenz蝴蝶就是这种低维底层结构的典型例子。然而，识别吸引子并不容易，许多尝试表明，由于固有的湍流，生成的模型往往与预测模式中使用的真实轨迹发散（Liu等人，参考文献Liu, Axas和Haller2024；?zalp等人，参考文献?zalp, N’ovoa和Magri2025a）。因此，在本文中，我们采用了一种随机方法来克服这一挑战。为此，我们使用了一种涉及编码-处理-解码架构的概率降维策略，在潜空间中使用Transformer并在VAE中生成随机性。这使我们能够学习过滤空间中系统的时间演变（即大尺度动态），并生成一组统计特性与过滤流动密切相似的随机轨迹。然后使用第二个概率模型来闭合较小涡旋的动态，生成在时间和空间上与真实分布（DNS数据）统计一致的新样本。这种“闭合”与超分辨率的核心思想相关，即从较粗的测量或模拟中重建高分辨率数据或场（Dong等人，参考文献Dong, Loy, He和Tang2016）。类似的超分辨率方法已成功应用于许多不同领域，包括图像处理（Kim等人，参考文献Kim, Lee和Lee2016；Wang等人，参考文献Wang, Yu, Wu, Gu, Liu, Dong, Qiao和Loy2018）、气候建模（Jiang等人，参考文献Jiang, Yang, Wang, Huang, Xue, Chakraborty, Chen和Qian2023；Soares等人，参考文献Soares, Johannsen, Lima, Lemos, Bento和Bushenkova2024）、流体动力学（Fukami等人，参考文献Fukami, Fukagata和Taira2023；Ward，参考文献Ward2023）以及医学成像（Huang等人，参考文献Huang, Shao和Frangi2017；Chen等人，参考文献Chen, Shi, Christodoulou, Xie, Zhou, Li, Frangi, Schnabel, Davatzikos, Alberola-López和Fichtinger2018；Zhang等人，参考文献Zhang, Hu, Philbrick, Conte, Sobek, Rouzrokh和Erickson2022），在这些领域中捕获详细的空间或时间模式对于准确分析和预测至关重要。对于这种小尺度闭合或“超分辨率”任务，许多研究利用了图像生成领域的突破性进展，如生成对抗网络（GAN）（De Souza等人，参考文献De Souza, Marques, Batagelo和Gois2023）和变分自编码器（VAE）（Kingma和Welling，参考文献Kingma和Welling2013）。在2020年代，去噪扩散模型作为重大突破出现（Ho等人，参考文献Ho, Jain和Abbeel2020；Song和Ermon，参考文献Song和Ermon2020）。它们为高分辨率图像合成提供了最先进的工具，尽管它们在基于物理的应用中的渗透仍然有限。值得注意的例子包括带有嵌入式Transformer的扩散模型。例如，Serrano等人（参考文献Serrano, Wang, Naour, Vittaut和Gallinari2024）在他们的AROMA框架中成功利用了这种方法来开发一个适应不同几何形状和尺度的强大且灵活的偏微分方程（PDE）仿真器。Mardani等人（参考文献Mardani, Brenowitz, Cohen, Pathak, Chen, Liu, Vahdat, Nabian, Ge, Subramaniam, Kashinath, Kautz和Pritchard2025）将确定性降阶模型（ROM）与基于扩散的闭合集成，以改进气候网格细化。扩散模型及其伴随的广泛研究在条件下的复杂概率分布学习中展示了令人印象深刻的能力；然而，它们在训练和推理方面仍然计算成本高昂（Salimans等人，参考文献Salimans, Mensink, Heek和Hoogeboom2024）。在这种情况下，我们提出了一种用于概率闭合的替代方法，利用高斯过程来实现小尺度动态的高效且富有表现力的表示。本文分为三个主要部分。首先，在第2节中，我们描述了Kolmogorov流动、从模拟中得到的数据以及过滤策略。接下来，在第3节中，我们介绍了降阶模型（ROM）及其在过滤轨迹上的性能。最后，在第4节中，我们介绍了高斯过程闭合并详细说明了其理论基础。我们还将其结果与其他概率基线进行了比较，这些基线的实现细节在附录A.3中。第5节介绍了结合了ROM和闭合的整个建模流程。第6节提供了结论性的评论。

2. 问题表述
2.1. Kolmogorov流动
考虑的测试案例是Kolmogorov流动，这是一种二维、不可压缩、非稳态的Navier–Stokes方程解，具有双重周期性边界条件，$\left[x,y\right]\in {\mathrm{\mathbb{R}}}^2=\left[0,2\pi \right]\times \left[0,2\pi \right]$。它由在$x$方向上施加的正弦驱动项驱动。控制方程为：
$$ {\displaystyle \begin{array}{c}\frac{\partial \mathbf{u}}{\partial t}=-\nabla p-\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}+\frac{1}{\mathit{\operatorname{Re}}}{\nabla}^2\mathbf{u}+\mathbf{f},\\ {}\nabla \cdot \mathbf{u}=0,\\ {}\mathbf{f}\left(x,y\right)=A\sin \left({n}_fy\right)\hskip0.1em {\mathbf{e}}_x\end{array}}
$$
其中$\mathbf{u}=\left(u,v\right)$是速度场，$p$是压力，$\nu$是运动粘度，$\mathit{\operatorname{Re}$是雷诺数，$\mathbf{f}$是外部驱动力。驱动函数$\mathbf{f}$在$x$方向引入了一个稳定的、空间周期性的体力，其波数为${n}_f$，固定振幅为$A$。这种流动在图1中展示，非常适合研究流体不稳定性以及向湍流的过渡（Kolmogorov等人，参考文献Kolmogorov, Levin, Hunt, Phillips和Williams1991），根据驱动波数${n}_f$和雷诺数$\mathit{\operatorname{Re}}$的不同，会显示出不同的状态（Platt等人，参考文献Platt, Sirovich和Fitzmaurice1991）。对于足够高的$\mathit{\operatorname{Re}}$，会表现出能量级联，其中 large coherent 结构衰变为较小的涡旋，这些涡旋贡献的能量较少，并最终被粘度耗散。这是本工作关注的领域。这里，我们考虑$\mathit{\operatorname{Re}}=34$和驱动波数${n}_f=4$，这一选择受到了Mo和Magri（参考文献Mo和Magri2025）以及Racca（参考文献Racca2023）最近研究的启发。这些参数导致了一个弱湍流，不同尺度以混沌的方式相互作用，激励我们通过统计视角进行研究。该模拟是在一个$64\times 64$的欧拉网格上进行的，使用的是公开可用的伪谱求解器kolSol，该求解器基于Canuto等人描述的傅里叶-伽辽金方法（参考文献Canuto, Quarteroni, Hussaini和Zang2007）。我们采用了$32$个傅里叶模式来进行分辨率处理。方程在傅里叶空间内通过四阶显式龙格-库塔（Runge-Kutta）积分方案求解，时间步长为$\delta t=0.01$，并且每$\delta t=0.12$就存储一次结果。速度场如图1a所示，是定义在${\mathrm{\mathbb{R}}}^{N_x\times {N}_y\times T}$上的矩阵，其中${N}_x={N}_y=64$，$T=1000$，这定义了总体的快照数量。所产生的动能信号$K$通过方程(2.2)计算，并在图1b中展示。(2.2) $$ {k}_t=\frac{1}{2{N}_x{N}_y}\sum \limits_{d=1}^{N_x{N}_y}\left({u}_{d,t}^2+{v}_{d,t}^2\right) $$ 图1. 速度场及其对应的动能信号。表1. 概率预测的评估指标。

能量级联是Kolmogorov和Richardson湍流理论的核心概念，它基于将流动的能量分解为不同大小的结构，或者等效地，分解为不同的频率，这也被称为波数。较大的结构对应于较低的波数，通常携带更多的能量。我们的方法依赖于在存在一个大尺度相干结构与它们分解成的较小涡流之间的二分法的基础上对这些尺度进行分离。这些结构的空间分布，或者更准确地说，它们在波数上的分布，是通过对方程进行二维（2D）傅里叶变换得到的，表示为：$$ U\left(x,y\right)\mapsto \hat{U}\left(m,n\right)={\mathcal{F}}_{2D}\left[U\right]\left(m,n\right) $$ 其中变换定义为：(2.3) $$ \mathcal{F}\left(m,n\right)={\int}_{-\infty}^{\infty }{\int}_{-\infty}^{\infty }U\left(x,y\right)\hskip0.1em {e}^{-j2\pi \left( mx+ ny\right)}\hskip0.1em dx\hskip0.1em dy. $$ 这里，$m$和$n$分别是$x$和$y$方向上的空间频率（或波数）。量$\mid \mathcal{F}\left(m,n\right)\mid$代表动能场$K$的幅度谱。它提供了在任何给定时间步下的波数分布。时间平均谱在对数尺度上显示在图2中。图2. 在对数尺度上显示的时间平均能量谱。频谱图表明，能量最高的结构集中在中心附近，对应于小波数。随着波数的增加，远离中心，幅度通常会减小，反映了较小尺度结构对动能的贡献逐渐减少。我们计算出的能量级联在图3中展示。图3. 保持90%能量的低通滤波器阈值。

2.2. 尺度分离策略 在这项工作中，我们提出将流动尺度分为两大类：大尺度结构和小尺度结构。这两类结构将采用两种不同的策略进行处理。我们的目标是确定级联中的一个能量阈值，这个阈值能够保留大约90%总能量的较大、“较简单”的结构。这个阈值在图3中进行了说明，对应于一个临界波数${k}_c$，然后使用这个阈值来构建一个低通滤波器$H$，以便将大尺度运动与小尺度波动分开。我们强调，选择这种滤波操作是有意为了满足尺度分离的需求。基于模态压缩的滤波器，例如投影到有限数量的主导POD模式上，并不会选择性地移除高频内容，因为这些频率可能仍然存在于保留的模式中。在这里考虑的低雷诺数（$\mathit{\operatorname{Re}}=34$）下，能量谱的惯性范围较窄。因此选择了90%的能量阈值，因为它有效地分割了这个惯性范围，保留了维持混沌动力学的模式，同时移除了向耗散尺度衰减的弱湍流波动。如果应用稍微激进一些的滤波器，比如85%，将会消除动态上重要的模式，包括那些接近驱动波数的模式，从而将动力学压缩到一个确定性的状态。对于在更高雷诺数下研究的系统，其中惯性范围明显更宽，最佳截止值的选择已经得到了广泛的研究。在这样的情况下，通常发现保留高达${k}_c\approx 0.2{\eta}^{-1}$的模式就足以确保随机稳定性（参考文献Inubushi, Saiki, Kobayashi和Goto2023；Matsumoto等人，参考文献Matsumoto, Inubushi和Goto2024）。(2.4) $$ H\left(m,n\right)=\left\{\begin{array}{ll}1,& \mathrm{if}\sqrt{m^2+{n}^2}\le {k}_c,\\ {}0,& \mathrm{otherwise}.\end{array}\right. $$ 将这个滤波器应用于傅里叶变换后的场，我们得到了过滤后的谱：(2.5) $$ {\hat{U}}_{\mathrm{low}}\left(m,n\right)=H\left(m,n\right)\cdot \hat{U} $$ 最后，通过在物理空间中对逆2D傅里叶变换，可以恢复过滤后的速度场${U}_{\mathrm{low}}\left(x,y\right)$：(2.6) $$ {U}_{\mathrm{low}}\left(x,y\right)={\int}_{-\infty}^{\infty }{\int}_{-\infty}^{\infty }{\hat{U}}_{\mathrm{low}}\left(m,n\right)\hskip0.1em {e}^{j2\pi \left( mx+ ny\right)}\hskip0.1em dm\hskip0.1em dn. $$ 这个操作隔离了对应于波数低于${k}_c$的大尺度结构，过滤掉了小尺度波动，如图4所示。虽然这种表述是在连续环境中提出的，但低通滤波器实际上是在离散傅里叶空间中使用离散傅里叶变换实现的。图4. 在波数空间的对数尺度上显示的低通滤波器和过滤后的能量谱。滤波器减少了流动的空间复杂性和细节特征，这些特征对于降阶模型（ROM）来说难以准确捕捉（Wang等人，参考文献Wang, Akhtar, Borggaard和Iliescu2012；Ahmed等人，参考文献Ahmed, Pawar, San, Rasheed, Iliescu和Noack2021）。滤波器对流动的影响在参考文献Chen, Shi, Christodoulou, Xie, Zhou, Li, Frangi, Schnabel, Davatzikos, Alberola-López和Fichtinger的图5中通过速度场和动能图（分别显示在图5a、b、c中）进行了说明。图5. 低通滤波器对速度场和动能的影响。过滤后的动能是从过滤后的速度场计算得出的，并不是直接对动能场应用低通滤波器得到的。这个过程产生的过滤速度场将用于训练下一节中描述的动态降阶模型（ROM）。

3. 第一项任务：预测过滤后的动力学 3.1. 模型架构 我们使用变分自编码器（VAE）来识别一个简化的潜在空间，然后通过变换器在这个空间上学习动力学，利用了嵌入式记忆和注意力机制的最新进展（Vaswani等人，参考文献Vaswani, Shazeer, Parmar, Uszkoreit, Jones, Gomez, Kaiser, Polosukhin, Guyon, Luxburg, Bengio, Wallach, Fergus, Vishwanathan和Garnett2017；Shokar等人，参考文献Shokar, Kerswell和Haynes2024）。Mori–Zwanzig形式主义和Takens定理（Mori，参考文献Mori1965；Zwanzig，参考文献Zwanzig1973；Takens，参考文献Takens, Rand和Young1981）强调了结合记忆项的重要性，以有效考虑未观测到的正交子空间对通过VAE投影学习到的可观测量集合的影响。此外，Gupta等人（参考文献Gupta, Schmid, Sipp, Sayadi和Rigas2023）展示了在有限的、不变的Koopman空间内使用纯马尔可夫表示来展开动力学的局限性。我们的架构如图6所示，与UP-dROM框架中使用的架构相同（Zighed等人，参考文献Zighed, Nicolas, Gallinari和Sayadi2025）。图6. 降阶模型架构，用于预测过滤后的轨迹。该模型在涉及二维不可压缩流动过圆柱体的确定性环境中展示了强大的性能和泛化能力。该ROM的理论基础和实现细节在GitHub上公开可用（Zighed，参考文献Zighed2025a）和UP-dROM中。特别有趣的是变分自编码器（VAE）的使用，因为它不仅能产生结构良好且密度高的潜在空间，从而增强泛化性能，还能学习潜在变量的分布。这使得生成多样化的同时物理上合理的轨迹成为可能。因此，这一概率框架自然支持随机采样和不确定性量化，正如UP-dROM中所展示的。因此，这种压缩策略比传统的POD等技术表现更好。它的优势源于其本质上的概率性表述，这允许我们量化不确定性，并利用其方差成分；同时通过放宽POD固有的线性约束，VAE学习了非线性的潜在流形，从而实现了复杂流动动力学的先进压缩性能（Zighed等人，参考文献Zighed, Nicolas, Gallinari和Sayadi2025）。VAE学习了一个映射$f:{\mathrm{\mathbb{R}}}^{64\times 64\times 2}\to {\mathrm{\mathbb{R}}}^{64}$，并且变换器通过两个注意力块共同训练来学习这个潜在流形上的时间动力学。损失函数、架构和超参数的细节在附录A.2（表A1）中提供，同时还提供了训练的算法伪代码。为了给模型提供足够的潜在数据分布的先验知识，我们使用八个具有相同统计特性但在初始条件略有不同的流动实现来训练它。每个实现中，只有前80%的1000个快照被用来构建训练和验证数据集，而剩余的20%用于评估模型再现长期统计行为的能力。训练-验证分割采用奇偶数方案进行，即一半选定的快照用于训练，另一半用于验证。因此，训练和验证集每个都包含400个快照，并覆盖相同的时间范围96个无量纲时间单位。我们用$\tau$表示这个时间范围。测试集在验证序列的基础上额外增加了100个快照，相当于$0.25\tau$，专门用于长期评估。一旦模型生成了超出$\tau$标记的预测，我们就将预测的流动统计与真实统计数据进行比较，以评估长期性能。为了保持表述的清晰简洁，我们展示了单个流动实例的结果。模型的任务是从其初始条件重建流动，并随后在扩展的预测范围内展开动力学，同时保持与它试图模仿的真实流动的统计一致性。我们强调，该模型也可以从更生成的角度进行评估，因为它完全有能力在提供之前未见过的初始条件时生成新的流动实例。然而，我们认为这种实验设置超出了本研究的范围。

3.2. 模型预测 在训练之后，模型可以根据过滤后的初始条件预测过滤后的轨迹。对于每个投影，变分自编码器（VAE）在潜在坐标上建模一个高斯分布。这个分布的熵取决于训练过程中学习的方差，反映了投影的不确定性。在每个时间步骤，都会抽取一个潜在向量以引入不确定性。然后这种不确定性通过模型的自回归结构传播。因此，每个预测都包含了之前步骤积累的不确定性以及当前随机投影引入的不确定性。然后我们解码这些潜在样本轨迹，将它们投影回物理的高维空间。这些因此构成了一个与流动的潜在概率分布一致的合理轨迹族。我们将这个族称为一个集合，其平均值代表最可能的流动演变，其标准差代表不确定性。如图7所示，随着动力学展开的进行，不确定性呈现出增长趋势。该图显示了平均轨迹以及由集合的标准差$\sigma$定义的置信区间。我们还展示了真实的验证轨迹，以便更好地评估生成集合的两个统计矩。图7. 在给定时间（$0.25\tau$）下，比较了两个流动场的预测轨迹和验证轨迹。该图还展示了在$2\tau$期间由ROM在不同位置生成的样本轨迹集合。实线红线代表ROM的平均预测，虚线黑色线条表示验证轨迹，阴影区域表示$\pm \sigma$和$\pm 3\sigma$的置信区间。随着时间的推移，由于模型的自回归性质，不确定性趋于累积，导致预测的方差增加。图7中报告的不确定性结合了随机效应和认知效应的贡献。认知不确定性源于模型外推，即当对远离训练分布的输入状态进行推理时。相比之下，偶然不确定性是由预测任务的内在难度以及混沌动力学的不可预测性引起的。因此，即使对于来自训练分布的输入，它也代表了可实现不确定性的下限。在预测范围$\tau$内，由于奇偶数的训练-验证分割策略，验证分布实际上与训练分布相同。尽管如此，仍然观察到一个非零的方差，这表明这种不确定性主要源于任务和数据的偶然性。随着时间的推移，预测轨迹可能会逐渐偏离训练分布。在这种情况下，可能会产生额外的认知不确定性，反映了与分布外状态相关的模型不确定性增加。尽管存在这种方差，我们观察到预测均值与验证轨迹之间有很强的一致性。此外，在不确定性显著增加之前，置信区间会适应演变的动态和转变。这种适应有效地平衡了区间宽度与真实值的覆盖范围。这种关系通过预测区间覆盖概率（PICP）进行了正式量化，如表1所示。PICP衡量了落在预测置信区间内的真实值的比例，由方程(3.1)给出：$$ \mathrm{PICP}=\frac{1}{N}\sum \limits_{i=1}^N\mathbf{1}\left({y}_i\in \left[{\mu}_i-{\sigma}_i,{\mu}_i+{\sigma}_i\right]\right), $$其中${\mu}_i$和${\sigma}_i$分别是预测均值和标准差。表1中报告的经验置信区间接近于理想高斯分布所期望的值，即$\pm \sigma$内为$68\%$，$\pm 3\sigma$内为$99\%$，这表明模型找到了遵循高斯行为的误差分布的均值和方差。这些区间代表了集成分布与其捕捉真实值能力之间的实际权衡。这种权衡进一步通过连续排名概率分数（CRPS）进行了量化，也在表1中报告。CRPS衡量概率预测的准确性，由方程(3.2)给出：$$ \mathrm{CRPS}\left(\mu, \sigma; x\right)=\sigma \left[\frac{1}{\sqrt{\pi }}-2\phi \left(\frac{x-\mu }{\sigma}\right)-\frac{x-\mu }{\sigma}\left(2\Phi \left(\frac{x-\mu }{\sigma}\right)-1\right)\right], $$其中$\phi$和$\Phi$分别表示标准正态概率密度函数（PDF）和累积分布函数（CDF）。较低的CRPS值表示更准确的概率预测，惩罚了预测不确定性的偏差和过度或不足的信心。这里报告的低CRPS值表明区间保守性和预测准确性之间有良好的平衡。表1还包括相对的${\mathrm{\ell}}_1$和${\mathrm{\ell}}_2$误差，这两个指标都突出了集成均值与真实轨迹之间的一致性。这些指标在方程(3.3)中定义：$$ \mathrm{Rel}\hbox{-} \mathrm{L}1=\frac{\sum \mid {y}_{\mathrm{true}}-{y}_{\mathrm{mean}}\mid }{\sum \mid {y}_{\mathrm{true}}\mid}\times 100\%,\hskip1em \mathrm{Rel}\hbox{-} \mathrm{L}2=\frac{\sqrt{\sum {\left({y}_{\mathrm{true}}-{y}_{\mathrm{mean}}\right)}^2}}{\sqrt{\sum {y}_{\mathrm{true}}^2}}\times 100\%. $$由于系统的混沌性质和模型中固有的随机性，直接在个别快照和真实值之间进行点对点比较并不总是有意义的。此外，我们的目标是评估超出真实流时间窗口的预测范围。因此，我们将通过统计一致性来评估模型的长期有效性。具体来说，我们分析生成轨迹在统计上是否与真实动态吻合，使用概率密度函数。我们计算生成的去噪速度场的概率密度函数（PDFs），并将其与真实去噪PDFs进行比较。这种频率主义方法确保了即使轨迹在实现上有所偏差，它仍然来自相同的底层统计分布，并且长期的动态展开与该分布保持一致。图8展示了这些PDFs的比较，显示出可接受的统计一致性。图8. 字段U和V的概率密度函数（PDFs）。最后，对生成数据和真实去噪轨迹计算了Wasserstein距离，分别为$U$的$0.0081$和$V$的$0.0269$。这表明生成的和参考分布之间存在强烈的统计一致性。这证实了这样一个结论：对于字段$U$和$V$，即使某个特定轨迹在点对点上偏离了真实轨迹，或者源自更长的预测展开，它仍然遵循生成真实字段的底层统计分布。4. 第二个任务：高斯过程回归—概率小尺度闭合 这第二部分描述了本工作的第二个任务：通过采用高斯过程回归（GPR）框架来提高作用在去噪流场上的降阶模型的保真度。高斯过程的优势在于可以使用矩阵运算进行精确推理。还值得注意的是，控制高斯过程形式的超参数是稀疏的，可以直接从数据中估计出来（Williams和Rasmussen，参考文献Williams, Rasmussen, Touretzky, Mozer和Hasselmo1995）。这导致了一个轻量级模型，能够快速训练以及准实时闭合和推理。目标是从低维去噪表示建立到全分辨率流场的概率映射。我们想强调的是，这种方法可以连续应用于几乎任何低保真数据，无论是由机器学习模型生成的还是像URANS这样的数值模拟等。第3节中介绍的方法与这里呈现的方法论无关，仅作为低保真数据的必要底层动力模型。4.1. 高斯过程 高斯过程（GP）定义了一个函数上的分布，表示为$$ f(x)\sim \mathcal{GP}\left(m(x),\mathcal{K}\left(x,{x}^{\prime}\right)\right), $$其中$m(x)$是均值函数（通常中心位于0），$\mathcal{K}\left(x,{x}^{\prime}\right)$是一个正定连续协方差（核）函数。给定训练数据$\mathcal{D}={\left\{\left({x}_i,{y}_i\right)\right\}}_{i=1}^N$，其中${y}_i=f\left({x}_i\right)+{\epsilon}_i$，且${\epsilon}_i\sim \mathcal{N}\left(0,{\sigma}_{\epsilon}^2\right)$，我们假设${\overline{y}}_i=0$。真实的底层函数$f$是未知的，${\sigma}_{\epsilon}^2$是表示观测噪声方差的超参数。设$X={\left[{x}_1,\dots, {x}_N\right]}^{\top}\in {\mathrm{\mathbb{R}}}^{N\times d}$为输入矩阵，$Y={\left[{y}_1,\dots, {y}_N\right]}^{\top}\in {\mathrm{\mathbb{R}}}^{N\times 1}$为相应的输出向量。核矩阵$\mathcal{K}_{XX}\in {\mathrm{\mathbb{R}}}^{N\times N}$编码了训练输入之间的成对相似性，给定(4.1)$$ \mathcal{K}\left({x}_i,{x}_j|\tau \right)={\sigma}^2\exp \left(-\frac{1}{2{l}^2}\parallel {x}_i-{x}_j{\parallel}^2\right), $$其中$\tau =\left\{{\sigma}^2,l\right\}$是模型的超参数；输出方差${\sigma}^2$和长度尺度$l$。长度尺度$l$控制函数的平滑度，输出方差${\sigma}^2$决定了函数值的典型尺度。这里使用了平方指数（RBF）核，因为它的平滑性很好，尽管对于表现出更粗糙或周期性行为的问题，Matérn或周期性核可能更合适。设$\mathcal{D}}^{\ast }={\left\{\left({x}_i^{\ast },{y}_i^{\ast}\right)\right\}}_{i=1}^M$表示一组新的输入，其中${X}^{\ast }={\left[{x}_1^{\ast },\dots, {x}_M^{\ast}\right]}^{\top}\in {\mathrm{\mathbb{R}}}^{M\times d}$。因此，高斯过程中的关键思想是训练输出$Y$和测试输入处的隐函数值${f}^{\ast }=f\left({X}^{\ast}\right)$作为多变量高斯联合分布，如(4.2)$$ \left[\begin{array}{c}Y\\ {}{f}^{\ast}\end{array}\right]\sim \mathcal{N}\left(0,\left[\begin{array}{cc}{\hat{\mathcal{K}}}_{XX}& {\mathcal{K}}_{XX^{\ast }}\\ {}{\mathcal{K}}_{X^{\ast }X}& {\mathcal{K}}_{X^{\ast }X}^{\ast }}\end{array}\right]\right), $$其中${\hat{\mathcal{K}}}_{XX}={\mathcal{K}}_{XX}+{\sigma}_{\epsilon}^2{I}_N$包括观测噪声。这意味着所有数据点都是从同一个高斯分布中抽取的部分观测值。给定联合高斯分布的部分观测值，可以计算未观测变量的条件分布，从而得到预测分布。新输入${X}^{\ast }$处的后验预测分布由方程(4.3)给出：$$ {f}^{\ast}\mid {X}^{\ast },\mathcal{D}\sim \mathcal{N}\left({\mu}_f,{\Sigma}_f\right), $$其中：$$ {\mu}_f={\mathcal{K}}_{X^{\ast }X}{\hat{\mathcal{K}}}_{XX}^{-1}Y,\hskip1em {\Sigma}_f={\mathcal{K}}_{X^{\ast }{X}^{\ast }}-{\mathcal{K}}_{X^{\ast }X}{\hat{\mathcal{K}}_{XX}^{-1}{\mathcal{K}}_{XX^{\ast }}. $$图9展示了一个2D示例，其中$N=M=1$。图9. 单个训练和推理点的高斯过程预测示例（$N=M=1$）。通过最大化给定输入$X$的训练输出$Y$的对数边缘似然来优化超参数$\tau$和${\sigma}_{\epsilon }$，使用方程(4.4)和(4.5)：$$ \log p\left(Y|X,{\sigma}_{\epsilon}^2,\tau \right)=\log \left[\int p\left(Y|\tilde{f},{\sigma}_{\epsilon}^2\right)p\Big(\tilde{f}|X,\tau \Big)d\tilde{f}\right] $$ $$ \log p\left(Y|X\right)=\log \left[\int \mathcal{N}\left(Y|\tilde{f},{\sigma}_{\epsilon}^2I\right).\mathcal{N}\Big(\tilde{f}|0,{\mathcal{K}}_{xx}\Big)d\tilde{f}\right] $$ $$ \log p\left(Y|X\right)=\log \mathcal{N}\left(Y|0,{\hat{\mathcal{K}}}_{XX}\right) $$其中(4.7)$$ {\displaystyle \begin{array}{c}{\hat{\mathcal{K}}}_{XX}={\mathcal{K}}_{XX}+{\sigma}_{\epsilon}^2I\\ {}\log p\left(Y|X\right)=-\frac{1}{2}{Y}^{\top }{\hat{\mathcal{K}}}_{XX}^{-1}Y-\frac{1}{2}\log \left|{\hat{\mathcal{K}}}_{XX}\right|-\frac{N}{2}\log \left(2\pi \right)} $$这个训练过程在算法1中总结。算法1 高斯过程超参数训练 1：输入：训练数据集$\mathcal{D}={\left\{\left({x}_i,{y}_i\right)\right\}}_{i=1}^N$ 2：初始化超参数$\tau$，${\sigma}_{\epsilon}^2$ 3：当对数似然未最大化时进行：$$ \log p\left(Y|X,\tau, {\sigma}_{\epsilon}^2\right) $$ 4：计算协方差矩阵$\mathcal{K}_{XX}$ 5：计算对数似然：$$ \log p\left(Y|X,\tau, {\sigma}_{\epsilon}^2\right) $$ 6：更新超参数$\tau$，${\sigma}_{\epsilon}^2$ 7：结束while在推理时，新点${x}^{\ast }$的预测分布如算法2所示计算。算法2 高斯过程推理 1：输入：训练数据$\mathcal{D}={\left\{\left({x}_i,{y}_i\right)\right\}}_{i=1}^N$，新的样本输入${x}^{\ast }$，以及训练好的超参数$\tau$ 2：计算核子矩阵：$$ \mathcal{K}=\left[\begin{array}{cc}{\hat{\mathcal{K}}}_{XX}& {\mathcal{K}}_{XX^{\ast }}\\ {}{\mathcal{K}}_{X^{\ast }X}& {\mathcal{K}}_{X^{\ast }X}^{\ast }}\right] $$ 3：计算后验均值和方差：$$ {\mu}_f={\mathcal{K}}_{X^{\ast }X}{\hat{\mathcal{K}}}_{XX}^{-1}Y,\hskip1em {\Sigma}_f={\mathcal{K}}_{X^{\ast }{X}^{\ast }}-{\mathcal{K}}_{X^{\ast }X}{\hat{\mathcal{K}}_{XX}^{-1}{\mathcal{K}}_{XX^{\ast }} $$ 4：预测输出：抽样${y}^{\ast}\sim \mathcal{N}\left({\mu}_f,{\Sigma}_f\right)$ 4.2. 应用于闭合的GPR 相比在物理空间中寻找映射（这将导致太多的自由度），我们利用数据固有的低维结构，并使用适当的正交分解（POD）在一个降维空间中提取映射。POD是一种数据驱动的线性降维技术，从系统中提取出主要的相干结构。我们将POD应用于快照矩阵$\varphi \in {\mathrm{\mathbb{R}}}^{T\times D}$，其中$\varphi$的每一列代表给定时间下的D维系统状态。使用奇异值分解（SVD），我们写成：(4.8) $$ \varphi =U\Sigma {V}^T $$其中$U\in {\mathrm{\mathbb{R}}}^{D\times r}$包含空间模式，$\Sigma \in {\mathrm{\mathbb{R}}}^{r\times r}$是对角矩阵，包含奇异值，$V\in {\mathrm{\mathbb{R}}}^{T\times r}$包含潜在坐标，且$r\le \min \left(T,D\right)$。我们定义$\Psi =U\Sigma$，将空间结构与其能量权重结合在一起。在秩$r=3$处进行截断，因为这些包含了过滤和全尺度数据集中98%的总能量。通过考虑更多的模式可以获得更好的结果，尽管训练成本更高，但不会修改实现。我们假设有一个不变的降维基$U$，以确保在长期动态展开过程中模态和统计的一致性。最后，我们使用符号${}^{\wedge }$来表示过滤数据，而全尺度数据则不加任何符号。图10展示了高斯过程回归（GPR）执行的映射。图10. 低通过滤POD空间和全尺度POD空间之间的映射。因此，系统在时间t的状态可以通过组合三个识别的POD模式来近似，对于过滤数据有(4.9) $$ \hat{\varphi}\approx {\hat{\alpha}\hat{\psi}}_1+{\hat{\beta}\hat{\psi}}_2+{\hat{\gamma}\hat{\psi}}_3, $$对于全状态数据有(4.10) $$ \varphi \approx {\alpha \psi}_1+{\beta \psi}_2+{\gamma \psi}_3, $$因此，我们的目标是学习一个时间独立的映射，即：(4.11) $$ f:\left(\hat{\alpha},\hat{\beta},\hat{\gamma}\right)\mapsto \left(\alpha, \beta, \gamma \right). $$而不是学习一个确定性函数，我们对$f$施加了一个先验：(4.12) $$ f\sim \mathcal{GP}\left(0,\mathcal{K}\right). $$因此，给定训练数据$\mathcal{D}$，我们推断出一个函数的后验(4.13) $$ p\left(f|\mathcal{D},\hat{\alpha},\hat{\beta},\hat{\gamma}\right). $$高斯过程回归（GPR）相比传统回归方法的主要优势在于其本质上是概率性的。通过从预测分布${y}^{\ast}\sim {f}^{\ast }$中抽样，GPR产生了一组输出，其统计属性与流动的底层分布一致。相关的方差自然量化了模型置信度，不仅提供了点对点预测，还提供了围绕它们的可信区间。这种考虑不确定性的建模方法与最近的风流建模发展一致，后者越来越关注学习流的统计表示而不是确定性轨迹。它也与本工作的目标一致，即在统计框架内进行重构。高斯过程在DNS轨迹的$\tau$部分进行训练和评估，采用85/15的序列分割方法。与使用奇偶快照划分的ROM不同，GP使用连续的、按时间顺序划分的方法。这样，训练和评估数据集分别包含680个和120个快照。鉴于可训练的超参数数量非常有限，使用如此有限的训练集是足够的。更重要的是，传统的GP随着训练样本数量的增加而表现不佳：训练需要计算核相似性矩阵$\mathbf{K}\in {\mathrm{\mathbb{R}}}^{m\times m}$的逆矩阵（或Cholesky分解），其中$m$是样本数量，导致计算成本为$\mathcal{O}\left({m}^3\right)$。因此，随着$m$的增加，计算变得非常困难。此外，在高维环境中，过大的训练集可能会由于成对相似性的崩溃而降低核的质量。由于这些原因，通常会在精心挑选的代表性子集上训练GP，并通过重新训练、下采样或特征选择（如自动相关性确定（ARD）来解决可扩展性问题。图11展示了我们实验设置中的高斯过程回归。图11a显示了训练和验证集中三个 reduced manifold 变量的高斯过程回归。这三个系数定义了一个在由前三个POD模式张成的三维空间中演化的轨迹，如图11b所示。图11显示了3个POD系数上的高斯过程回归（GPR）训练和验证。在反向投影到高维输入空间后，使用多种指标来评估模型性能，这些指标量化了预测的准确性和验证集上不确定性估计的质量。我们使用了相对L1误差、相对L2误差以及第3.2节中介绍的PICP来计算置信区间$\pm \sigma$和$\pm 3\sigma$，以及CRPS。它们记录在表2中。表2显示了验证集上的评估指标。高斯过程的表现表明，它在训练过程中获得了强有力的先验知识，使其能够做出具有可靠置信区间的有意义预测。验证集上的第一时刻得分略高于10％，这表明其对数据的忠实度相当高。预测区间覆盖概率（PICP）达到了82％，使用的是一个标准差构造的区间。较低的连续排名概率得分（CRPS）表明，这种覆盖率是在保持最小区间宽度的情况下实现的。训练后，模型确定了三个最优超参数集${\sigma}^2$、$l$和${\sigma}_{\epsilon}^2$，这允许构建图12所示的相似性核。因此，可以对新到达的数据进行推理。图12显示了训练样本核$\mathcal{K}_{XX}$。相似性核的结构表明，即使快照在时间上相隔较远，它们之间也高度相关。这种强相关性表明它具有很好的泛化到未见数据的潜力。如前一节所述，高斯过程被训练用来将过滤后的数据映射到全尺度的直接数值模拟（DNS）数据。最初在真实过滤轨迹上进行训练的高斯过程，在全模型格式中应用时，可以从由动态降阶模型（ROM）生成的过滤轨迹中预测全尺度数据。这建立了一个端到端的数据管道，其中输入是过滤后的初始条件，输出是生成的全尺度轨迹，如图13所示。图13显示了数据管道。当我们将高斯过程（GP）模型条件化为新的、未观察到的推断数据${X}^{\ast}$（由降阶模型（ROM）生成）时，我们得到了后验分布：$$ {f}^{\ast}\sim \mathcal{GP}\left({\mu}_f,{\mathcal{K}}_{XX^{\ast }}\right)，$$其中${\mu}_f$表示后验均值，${\mathcal{K}}_{XX^{\ast }}$代表更新的协方差矩阵。这个矩阵捕获了训练数据$X$和推断点${X}^{\ast}$之间的相似性。从这个后验中，我们可以抽取样本${f}^{\ast}$，这些样本代表预测的全尺度模式系数。GPR将这些系数视为来自联合分布的独立样本；因此这一步完全不依赖于时间。生成的输出成功重建了过滤掉的频率，如图14所示的平均频谱图所示。图14显示了波数空间的平均能量频谱图。降阶模型（ROM）和高斯过程（GP）的结合能够生成新的轨迹，尽管这些轨迹可能与真实轨迹有所不同，但在统计上与底层动力学保持一致。在任何给定的时间快照中，GP都能够重建缺失或未解析的流动结构，如图15所示的小尺度涡旋。图15显示了T = 120时刻的能量。显然，与图15a,c所示的流动进行一对一的图像比较看起来是不同的。超出某个预测范围后，像素对像素的比较变得不相关，因为我们的框架并不是简单地学习和复制特定轨迹，而是学习流动的底层统计特性，因此不应该期望重建完全相同的实现。5.1. 对长期预测的评估为了评估所提出框架在长期预测范围内的稳定性，我们评估了生成流动的统计特性是否随时间保持不变。具体来说，我们比较了预测的全尺度流动的概率密度函数（PDF）与从测试数据集获得的PDF。训练和验证数据集的时间范围为96个无量纲时间单位，我们用$\tau$表示。测试数据集覆盖了更长的时间范围$1.25\tau$，而ROM+GP框架用于生成延伸到$4\tau$的轨迹。图16显示了在二维网格上计算的动力能量$k$，并通过对$x$方向平均将其简化为一维信号。预测的能量PDF在长度为$\tau$的连续时间窗口上进行评估，并与从测试数据集计算的参考PDF进行比较。图16还报告了在应用高斯过程闭合之前获得的动能密度。图16显示了沿$x$方向平均的动能$k=\frac{1}{2}\left({u}^2+{v}^2\right)$的时空演化。颜色刻度表示$\left\langle k\left(y,t\right)\right\rangle_x$作为墙法向坐标$y$和无量纲时间的函数。每个$\tau$间隔的数据分布与测试集的PDF进行了比较。我们观察到，即使在预测的时间范围是训练/评估窗口的四倍时，闭合也固定了过滤后的展开能量统计特性，并在预测模式中保持不变，没有显著漂移。这一点得到了图17的支持，为了量化这一观察，我们在图17a中计算了预测能量密度与测试集分布之间的Wasserstein距离，使用的是$0.5\tau$的窗口。图17显示了长期展开的统计一致性，通过能量分布的Wasserstein距离来量化，以及混沌动力学的保持。图17a的结果表明，流动的统计特性在整个预测窗口内保持稳定。尽管在比训练/评估窗口长得多的时间范围内演化，模型仍然保持了空间和时间复杂性，同时保持了动力学的混沌特性。图17b表明，在耗散信号的延迟嵌入空间中表示的混沌吸引子上，动力学得到维持，其中$\Delta t$为2个快照。耗散计算为$\varepsilon (t)=\frac{1}{\operatorname{Re}}\left\langle \parallel \nabla \mathbf{u}{\parallel}^2\right\rangle$，在整个域上取平均。从动力学角度来看，长期预测也与数值模拟的Kolmogorov流动的预期行为一致，这一点得到了图18的支持。图18a比较了直接数值模拟（DNS）和生成流的长期展开得到的湍流动能（TKE）频谱图。在惯性范围内，两种谱图之间有很好的一致性。图18a还显示，生成的流动尽管表现出略微更多的耗散行为。然而，这种差异是适度的。考虑到为了实现闭合而进行的降维，这解释了这种差异的起源。闭合是在将动力学投影到前三个POD模式之后应用的。尽管这些模式捕获了几乎所有的系统能量，但由于保留的模式数量较少，以及投影的线性特性，截断了与高波数结构相关的模式。图18a还显示，这种截断精确再现了我们模型重建后观察到的波长损失。这意味着在重建的流动中观察到的额外耗散不能归因于过滤策略、ROM或闭合机制本身，而是因为选择了用于闭合的降维基础。图18b显示了使用三个POD模式进行闭合时从ROM+GP程序获得的二维耗散流形。该流形是使用耗散信号的延迟嵌入构建的，其中$\Delta t=4$。尽管全局结构和拓扑与湍流吸引子的结构一致，但耗散的差异变得更加明显。图18c和18d显示，将POD模式的数量从3增加到18可以校正湍流动能谱图和相关的耗散流形。这种改进是以更高的计算成本和第一时刻准确性的降低为代价的（${L}_1$和${L}_2$误差分别为$73.6\%$和$75.5\%$）。然而，第二时刻的准确性仍然令人满意，${PICP}_{\pm 3\sigma }=87.7\%$。尽管第一时刻误差不是主要目标，但仅保留三个POD模式是在重建能力和统计及动力学一致性之间的一个合适折中。5.2. 与替代机器学习架构的比较我们将基于GPR的模型性能与最先进的概率机器学习预测器进行了比较，即变分自编码器（VAE）和扩散模型。选择这些统计分布模拟器的原因是它们在图像生成和超分辨率任务中的出色能力。在这种背景下，多尺度闭合任务是使用这些ML架构来执行的。附录A.3提供了每种架构的实现细节，以及选择它们作为基准的进一步理由。表3总结了VAE和扩散模型在验证集上的结果。两者都是在与GP相同的流动分辨率上训练的，并且具有相同的训练和验证集。所有性能测量都是在配备了NVIDIA RTX 4000 Ada Generation GPU的工作站上获得的，该GPU用于训练VAE和扩散模型。高斯过程回归是在第13代Intel i9-13950HX CPU（2.20 GHz）上训练的。表3比较了我们的高斯过程回归、变分自编码器（VAE）和扩散模型之间的评估指标。高斯过程在验证集上的第一和第二时刻指标上都优于这两个基准。然而，尽管扩散模型在图像生成方面表现出最先进的性能，并且不依赖于关于数据结构的任何先验假设，但其训练和推理成本显著较高。具体来说，生成单个图像需要许多去噪步骤（在我们的实现中为1000步），这意味着模型必须进行一千次推理才能产生一个独特样本。这导致了显著的计算开销。我们假设，使用更大的模型和更彻底的搜索最优超参数，扩散模型的表现可能会比这里展示的更好，尤其是在有更多数据的情况下。然而，这样的搜索将产生高昂的计算成本。VAE需要对每个样本进行一次推理，虽然可以实现更快的预测，但潜在的计算成本也会很高，尤其是对于深度编码器和解码器而言。相比之下，GP可以一次性生成所有样本，直接在降维流形空间中，将不同的快照视为单一高斯分布的部分观测。值得注意的是，在训练过程中只优化了GP的三个超参数。由于我们预测的是三维的模态状态，我们为每个输出维度训练一个GP。这样总共就有三个GP和九个优化的超参数，从而得到了一个更加高效和可解释的模型。

6. 结论

本研究介绍了一个用于长期生成和预测代表湍流状态的混沌轨迹的预测建模框架。所提出的方法采用了一种完全随机的公式，利用尺度分离来再现原始系统的统计特性。在这个框架中，首先通过数据过滤识别出连贯的大尺度动态，然后通过随机表示纳入小尺度波动的影响。这种简化使得降阶模型（ROM）能够专注于时间依赖动态的演化，而不是精细尺度特征的详细重建。一旦过滤后的轨迹被推断出来，一个次级的、时间独立的模型会将这些结果概率性地映射回全尺度空间。值得注意的是，这两个模型是独立的，并且被设计成模块化的：ROM可以模拟大尺度动态，同时与替代的闭合方案耦合，而小尺度学习策略可以作为替代的低保真动态模型的闭合机制。这种模块化架构旨在为表现出尺度分离的应用提供灵活性和适应性。在这个框架中使用的降阶模型（ROM）基于我们最近的工作UPdROM（Zighed等人，参考文献Zighed, Nicolas, Gallinari and Sayadi2025）。尽管它最初是在确定性的流体经过钝物体时进行训练和验证的，但模型的随机性质使其能够生成轨迹集合，其统计特性与Kolmogorov流分布非常吻合。这里将这种分布视为基准验证案例。模型使用变分自编码器（VAE）学习到低维流形的概率映射。这个潜在空间中的动态通过变换器架构来捕捉。此外，通过在每个时间步从VAE中随机采样，模型生成了轨迹集合。集合均值忠实于真实数据，而方差显示了强大的覆盖特性，这一点通过预测区间覆盖概率（PICP）得到了确认。集合扩散与覆盖之间的权衡通过接近零的连续等级概率分数（CRPS）来量化。为了闭合空间结构，我们使用高斯过程（GP）将过滤空间映射到各自POD子空间中的全尺度空间。仅使用三个POD模式，我们就可以忠实地重建多尺度数据，这一点通过强大的测试集性能得到了证明。这包括准确的一阶统计特性（具有较低的${\mathrm{\ell}_1$和${\mathrm{\ell}_2$误差）以及稳健的二阶统计特性（如PICP和CRP所证明的）。基于这种尺度分离的随机架构，我们展示了一个预测模型，它在扩展的预测范围内保持统计一致性和预测准确性。所提出的随机预测框架的效率为实时训练、微调和推理打开了新的可能性，特别是在流控制应用的背景下。我们将这些结果与其他概率映射基线进行了基准测试，具体来说是VAE和去噪扩散模型，使用了它们的标准实现。GP在所有考虑的指标上都优于这两种方法。鉴于扩散模型通常被认为是图像生成和超分辨率等任务的黄金标准，这与我们的问题密切相关，这一点尤其令人惊讶。然而，标准的去噪扩散算法通过递归去噪白噪声图像来推断单个快照，通常每个快照需要一千次去噪步骤。这意味着，为了生成一个干净的样本，模型必须被推断一千次。相比之下，VAE每个快照只需要一次推理，使其更具成本效益。此外，GP通过单次推理生成整个动态过程，因为所有快照都被视为来自同一多元高斯分布的部分观测。此外，这个分布仅由三个优化的参数参数化，确保了快速训练和显著更快的推理。我们承认，所提出程序的强大性能部分归因于测试案例的简单性，该案例由低雷诺数下的二维不可压缩流动组成。因此，应将本研究视为湍流流动建模的概念验证。我们相信，所提出的方法可以扩展到更复杂的配置中，尽管计算成本会增加，并且可能需要不同的闭合策略，例如基于分数的模型（例如扩散模型），而无需对框架进行任何根本性的概念重构。未来的工作将集中在三维的、现实世界的湍流流动上，以进一步证明我们解耦的概率ROM闭合框架的鲁棒性和可扩展性。

数据可用性声明

第3节中使用的降阶模型（ROM）在GitHub上公开可用，地址为https://github.com/IsmaelZig-SU/p-DynSys_EncoderTransformerDecoder。包括源代码和本研究中使用的实验数据集的完整建模流程在另一个GitHub仓库中也可获取，地址为https://github.com/IsmaelZig-SU/Scale-Separation-and-stochastic-closure-for-Chaotic-Dynamics。

作者贡献

概念化-平等：I.Z., T.S.；数据 Kuratorial-平等：I.Z.；形式分析-平等：I.Z.；调查-平等：I.Z., T.S.；方法论-平等：I.Z., T.S.；资源-平等：N.T., T.S.；监督-平等：N.T., P.G., T.S.；验证-平等：I.Z., N.T., P.G., T.S.；可视化-平等：I.Z.；写作 – 原稿-平等：I.Z.；写作 – 审查 & 编辑-平等：N.T., P.G., T.S.

资金声明

该项目由欧盟的Horizon Europe研究与创新计划Cofund SOUND.AI共同资助，根据Marie Sklodowska-Curie资助协议编号101081674。

利益冲突

作者声明没有利益冲突。

伦理标准

该研究符合所有伦理准则，包括遵守研究国家的法律要求。

附录A.1. 频率截止选择

我们选择的截止频率，在波数空间的总动能的90%，是由我们的流动配置所决定的。本研究考虑了Racca（参考文献Racca2023）在低雷诺数（$\mathit{\operatorname{Re}}=34$）下的Kolmogorov流动配置，即接近湍流的起始阶段。在这种状态下，粘性效应与惯性效应相当，导致强烈的耗散和能量谱上一个非常狭窄的惯性范围。能量在强迫波数处注入，并转移到更高波数，然后迅速耗散。因此，动态由少数相互作用的模式控制。选择这种状态是为了在一个大尺度连贯结构与弱混沌波动共存的设置中测试尺度分离假设。图A1显示，在这种低雷诺数状态下，耗散流形$\left(E,\varepsilon \right)$对频谱截止非常敏感。当低通滤波器保留了90%的总动能时，得到的吸引子保持了一个高维的混沌结构，这在定性上与未过滤的动态一致。相比之下，如果将阈值降低到80%，耗散流形会从混沌云状结构塌缩为几乎一维的结构。这表明过滤过程不仅去除了小尺度波动，还去除了与维持混沌行为相关的动态关键模式，使系统进入确定性状态。在更高雷诺数下，惯性范围更宽，能量和耗散尺度更好地分离，基于能量的截止选择（例如80%对90%）更符合经典的尺度分离论点，这些论点通常建议保留直到${k}_c\approx 0.2\hskip0.1em {\eta}^{-1}$（其中$\eta$是Kolmogorov长度尺度）的模式足以确保随机稳定性（Matsumoto等人，参考文献Matsumoto, Inubushi and Goto2024）。图A1显示了在保留100%总动能（左）、90%（中）和80%（右）的频谱截止下，耗散流形（$k,\epsilon$）的演化。注意从混沌云到一维极限环的转变，表明了不稳定动态的抑制。

A.2. 过滤后的动态—ROM架构

总损失是VAE和变换器损失组分的加权和，包括两个模型的联合前向传递。联合训练的动机是将学习到的投影与潜在空间动态联系起来。在这种设置中，VAE不仅为了重建精度进行优化，还为了产生具有动态信息性的潜在表示，适合自回归建模。这促进了潜在空间的平滑、可预测的时间演化，有利于长期预测。从概念上讲，这相当于学习一组有限的可观测量，这些可观测量近似于Koopman不变子空间，符合Takens嵌入定理的精神，同时保持有限维度和数据驱动。损失是在长度为$q$的回望窗口和长度为$h$的预测范围内计算的。它结合了重建、潜在空间正则化和动态预测误差，定义为：
$$ {\displaystyle \begin{array}{l}\mathrm{\mathcal{L}}=\frac{\lambda }{m}\sum \limits_{k=0}^m\left(\parallel {\phi}_{n-q:n}-\mathcal{D}\mathrm{\mathcal{E}}\left({\phi}_{n-q:n}\right)\parallel +\beta \hskip0.1em \mathrm{KLD}\right)\\ {}\hskip3em +\frac{1}{m}\sum \limits_{k=0}^m\parallel {z}_{n+1:n+h}-{\hat{z}}_{n+1:n+h}\parallel \\ {}\hskip3em +\frac{1}{m}\sum \limits_{k=0}^m\parallel {\phi}_{n+1:n+h}-\mathcal{D}\left({\hat{z}}_{n+1:n+h}\right)\parallel .\end{array} $$
这里，$\mathrm{\mathcal{E}}$和$\mathcal{D}$分别表示编码器和解码器，$\phi$代表过滤后的状态，$z$代表相应的潜在变量，$\hat{z}$代表变换器模型产生的潜在预测。假设潜在变量遵循高斯分布$$ \hat{z}\sim \mathcal{N}\left(\mu, \Sigma \right)，$$其中$\Sigma$是对角协方差矩阵，对角线元素为${\sigma}^2$，编码器将输入状态映射到潜在分布参数：$$ \mathrm{\mathcal{E}}:\phi \mapsto \left(\mu, \sigma \right)。$$系数$\lambda$和$\beta$是通过参数扫描选择的正则化参数。损失的第一项强制编码器和解码器之间的一致性，其行为类似于传统的降阶建模中的非线性投影。Kullback–Leibler散度项促进了潜在空间的平滑和良好结构，定义为：
$$ \mathrm{KLD}=-\frac{1}{2}\sum \limits_{i=1}^Z\left({\sigma}_i^2+{\mu}_i^2-1-\log \left({\sigma}_i^2\right)\right)，$$其中$Z$表示潜在空间维度。第二项惩罚变换器在预测范围内产生的潜在变量的预测误差，第三项衡量解码预测潜在变量后的重建误差。训练程序在算法A1中概述。

算法A1 VAE–Transformer ROM框架的训练程序

要求：
- 过滤后的快照${\left\{{\phi}_t\right\}}_{t=1}^T$，
- 回望窗口$q$，
- 预测范围$h$，
- 批量大小$m$，
- 正则化权重$\lambda, \beta$。

确保：
- 训练好的编码器$\mathrm{\mathcal{E}$，
- 解码器$\mathcal{D}$，
- 和变换器模型$\mathcal{T}$。

1: 对于每个训练周期执行以下操作：
2: 对于每个$m$个快照的批次执行以下操作：
3: 提取批次过滤后的状态${\phi}_{n-q:n}$和相应的未来潜在状态${\phi}_{n+1:n+h}$；
4: 对过滤后的快照进行编码：${\z}_{n-q:n}=\mathrm{\mathcal{E}}\left({\phi}_{n-q:n}\right)$；
5: 对未来的过滤后快照进行编码：${\z}_{n+1:n+h}=\mathrm{\mathcal{E}}\left({\phi}_{n+1:n+h}\right)$以评估变换器；
6: 解码以重构过滤后的快照：${\hat{\phi}}_{n-q:n}=\mathcal{D}\left({z}_{n-q:n}\right)$；
7: 为变换器初始化自回归输入：${\z}_{\mathrm{input}}={z}_{n-q:n}$；
8: 初始化预测的潜在状态列表：${\hat{z}}_{n+1:n+h}=\left[\hskip0.1em \right]$；
9: 对于$t=1$到$h$执行以下操作：
10: 预测下一个潜在状态：${\hat{z}}_{n+t}=\mathcal{T}\left({z}_{\mathrm{input}}\right)$；
11: 将预测添加到列表中：${\hat{z}}_{n+1:n+h}.\mathrm{append}\left({\hat{z}}_{n+t}\right)$；
12: 更新下一步的输入窗口：${\z}_{\mathrm{input}}=\left[{\hat{z}}_{n+1:n+h}\left[1:\right],{\hat{z}}_{n+t}\right]$；
13: 结束；
14: 计算总损失$\mathrm{\mathcal{L}}$；
15: 反向传播梯度并更新$\mathrm{\mathcal{E}$、$\mathcal{D}$和$\mathcal{T}$的参数；
16: 结束；
17: 结束；
18: 返回训练好的模型$\mathrm{\mathcal{E}$、$\mathcal{D}$和变换器：

表A1. 模型和数据配置参数

A.3.全面关闭——替代架构
第5.2节中作为基准模型使用的扩散模型受到了原始去噪扩散概率模型（DDPM）框架的启发（Ho等人，参考文献Ho, Jain和Abbeel2020）。该模型经过训练，能够迭代地去除受高斯噪声污染的数据。去噪网络采用了U-Net风格的自动编码器，使用最大池化层和大小为3的卷积核，逐步将输入图像（尺寸为64×64）和两个输入通道（u,v）下采样到尺寸为4×4的潜在表示，其中包含256个特征通道。正式来说，前向去噪过程由噪声调度{\left\{{\beta}_t\right\}}_{t=1}^T定义，T=1000表示去噪步骤数，通过逐步向干净数据{x}_0中添加噪声来生成噪声样本{x}_t，具体公式为：
q(x_t|x_0)=N(x_t;√{\overline{\alpha}_t}x_0,(1-{\overline{\alpha}_t)\mathbf{I})，
其中{\overline{\alpha}_t={\prod}_{s=1}^t(1-{\beta}_s)。该模型通过最小化均方误差损失来预测每一步添加的噪声ε：
\mathrm{\mathcal{L}}={\unicode{x1D53C}}_{x_0,\epsilon, t}[\parallel \epsilon -{\epsilon}_{\theta}(x_t,t){\parallel}^2]，
其中{\epsilon}_{\theta}是网络在时间步t预测的噪声。推理过程是通过从条件为过滤后数据的高斯噪声开始的逆过程采样来执行的，通过迭代去噪重建干净样本。训练过程中使用了Adam优化器，共进行了2000个周期。从概念上讲，扩散模型学习了逐渐去除噪声的逆条件高斯转移的参数化表示。这些高斯转移是直接在原始的高维输入空间中定义的。

第5.2节中作为另一基准模型使用的变分自动编码器（VAE）的出发点是假设数据分布位于一个低维流形上，这在动态系统应用中很常见。与扩散模型不同，VAE还旨在通过学习参数(\mu_{\theta},{\Sigma}_{\theta})来近似给定观测数据{x}_{\mathrm{data}}的潜在变量的后验分布，这些参数对应于一个维度远低于原始输入空间的潜在空间中的高斯分布（Kingma和Welling，参考文献Kingma和Welling2013）。这种低维表示有助于更高效的推理，并能够捕捉数据的基本结构。然而，它依赖于这样一个低维投影存在的强假设。本研究中使用的VAE是一个全连接神经网络，它逐渐将输入数据从64×64的维度压缩到4维的潜在空间。该模型通过优化变分下界来进行训练，表达式为：
\mathrm{\mathcal{L}}=\mathrm{KLD}(q_{\theta}(z|X)\parallel p(z))-\unicode{x1D53C}_{z\sim q_{\theta}}[\log p_{\theta}(X|z)]，
其中q_{\theta}(z|X)是近似后验，p(z)是先验，p_{\theta}(X|z)是似然。实际上，我们在Kullback–Leibler散度项上引入了一个权重因子，以优先考虑重构的准确性而不是严格遵循潜在空间中的高斯性，从而在保真度和正则化之间取得平衡。本例中使用的权重因子为0.01。选择0.01这个值是通过实验得出的，旨在平衡重构准确性和潜在空间正则化。研究发现，较小的KL权重足以在不降低重构质量的情况下促进潜在分布的平滑性，这对于准确的下游动态预测和关闭映射至关重要。

注释：
本文因透明实践获得了开放数据和开放材料的认证。详情请参阅数据可用性声明。