基于神经网络与整体凸限制的浅水方程亚网格参数化方法

1. 项目概述：当机器学习遇上流体力学

在计算流体力学和气候建模领域，我们常常面临一个根本性的矛盾：物理世界中的流动现象，从微小的涡旋到全球尺度的洋流，跨越了巨大的时空尺度。然而，计算机的算力是有限的，我们不可能在模拟中解析每一个分子、每一个厘米尺度的湍流。这就好比你想画一幅城市地图，不可能把每一块地砖的纹路都画出来，但你又必须让地图使用者感受到不同街区道路的宽窄、建筑的疏密所带来的交通影响。这个“地砖纹路”对“街区交通”的影响，在流体模拟中，就是亚网格尺度过程对可解析尺度流动的影响。

传统的做法是，科学家们根据物理直觉和经验，构建参数化公式来近似这种影响。比如，在大气模型中用一个简单的公式来估算云的形成对热量输送的影响。但问题在于，很多小尺度过程极其复杂、高度非线性，用简单的公式去“猜”，往往不准，还会把误差累积起来，导致长期气候预测出现偏差。这就像用一把刻度粗糙的尺子去测量精密零件，短期看个大概还行，长期累积的误差足以让整个机器运行失常。

近年来，机器学习，特别是神经网络，为解决这个老难题提供了新思路。神经网络像个拥有超强学习能力的“黑匣子”，给它看足够多的高分辨率模拟结果（相当于“标准答案”），它就能自己摸索出那些难以用公式描述的、小尺度与大尺度之间的复杂关系。我们的工作，正是将这种思路应用于浅水方程的长期湍流模拟中。浅水方程是描述大气和海洋大尺度运动的基础模型之一，可以看作是一个简化版的“地球流体运动模拟器”。

但直接套用神经网络会带来新问题：它可能为了拟合数据而“放飞自我”，输出一些物理上不可能的结果，比如出现负的水深（现实中水怎么可能有负的深度？），或者在激波（类似水跃、海啸波前）附近产生剧烈的数值振荡，让模拟结果失真甚至崩溃。

因此，我们项目的核心创新在于“结合”。我们不仅用前馈神经网络去学习亚网格通量，更关键的是，将神经网络的输出与一个名为整体凸限制的传统数值方法“焊接”在一起。MCL就像一个严谨的“物理规则审查官”，它会检查神经网络给出的修正值，确保最终的数值解始终落在物理允许的集合内（例如水深非负）。这样，我们既利用了机器学习捕捉复杂模式的能力，又用坚实的数值方法保障了结果的物理合理性和数值稳定性。最终，我们得到的是一个局部四点的参数化方案——每个粗网格点的计算，只依赖它自己和左右邻居的信息。这种局部性带来了巨大的优势：计算高效，易于并行，并且可以针对不同区域（比如近岸湍流剧烈区和开阔海平静区）进行独立的训练和优化，而不必重新训练整个全球模型。

简单来说，我们做的是给粗网格的流体模拟装上一个由“AI学习+物理规则校验”共同驱动的“智能插件”。这个插件能让粗网格模拟表现出接近细网格的复杂流动特征，特别是改善不同尺度间的能量传递，从而在计算成本可控的前提下，大幅提升长期模拟的可靠性。无论你是从事气候模拟、海洋预报，还是对高精度计算流体力学感兴趣的研究者或工程师，这套方法都提供了一条融合数据驱动与物理约束的新路径。

2. 核心思路与方案设计

2.1 问题根源：从“完全解析”到“不得不粗化”

我们面对的起点是浅水方程的精细离散化。想象一个长度为L的计算域，我们用Nf=1024个非常细的网格单元将其剖分，网格尺寸Δx = L/1024。在这个精细网格上，我们使用局部Lax-Friedrichs格式进行空间离散，这是一种经典且稳健的数值格式，其核心思想是在计算网格界面通量时，引入一个与当地波速相关的数值粘性项，以抑制解在间断（如激波）附近的非物理振荡。时间推进则采用二阶精度的Heun方法（一种改进的欧拉法）。在这个设置下，模拟可以认为是“完全解析”的，它捕捉了我们关心的所有物理尺度，我们将它视作“真理”或“参考解”，记作DNS。

然而，1024个网格点的计算成本对于单次短期模拟或许可接受，但对于需要运行数十年、数百年的气候模式，或者需要进行大量参数敏感性研究的场景，这个成本就变得难以承受。我们必须进行粗化。

2.2 粗化策略：空间平均与通量闭合难题

粗化的直观做法是“合并单元格”。我们每8个连续的细网格单元（n=8）合并成一个粗网格单元，这样粗网格单元数Nc = Nf / n = 128，网格尺寸ΔX = 8Δx。粗网格上的变量，如平均水深HI和平均流量QI，就是其包含的8个细网格变量的算术平均。

现在，关键问题来了：当我们对精细网格上的方程进行空间平均，试图推导出粗网格变量满足的方程时，会发现方程中出现了一些“多出来的项”。这些项来源于细网格变量在粗网格界面处的亚网格尺度通量。具体来说，在粗网格界面处，真实的通量依赖于界面两侧细网格上的精确解，而我们只有粗网格的平均值。这就导致粗网格方程不封闭——方程中出现了未知的亚网格通量项。

传统的参数化方法（如增加一个简单的涡粘性项）就是试图用一个基于粗网格变量的经验公式来近似这个未知项。但正如前所述，对于复杂的非线性湍流，这种经验公式往往力不从心。

2.3 我们的方案：神经网络作为“通量学习器”

我们的核心思路是：让神经网络来学习这个从粗网格状态到亚网格通量的映射关系。

1. 输入与输出设计：我们不是让神经网络去猜整个复杂的通量，而是聚焦于非线性通量部分。我们将一个粗网格界面两侧各两个粗网格单元的状态（共4个点，每个点包含水深h和流量q两个变量，总计8个输入）作为神经网络的输入。神经网络的输出是这个界面处的非线性亚网格通量的两个分量（对应质量和动量）。这就是所谓的“四点模板”。

   • 为什么是四点，而不是两点？初期实验表明，仅使用界面左右两个粗网格点（两点模板）的信息，神经网络无法准确学习通量，尤其是在解梯度较大的区域。四点模板（左二、右二）提供了关于解局部曲率的更多信息，让网络能更好地区分平缓区域和可能存在激波的陡峭区域。这有严格的数学支撑，在随机多尺度系统的模态约简理论中，四点模板是保证近似精度所必需的。

2. 网络架构选择：我们采用标准的前馈神经网络，结构为输入层（8神经元）→ 3个隐藏层（每层128神经元）→ 输出层（2神经元）。激活函数选用GELU，而非更常见的ReLU。

   • 为什么是GELU？ReLU函数在负半轴直接输出0，这可能导致“神经元死亡”问题，影响梯度流动。GELU函数是ReLU的平滑版本，它根据输入值的大小，以概率方式对输入进行门控（缩放），对于负值不是简单截断，而是进行平滑抑制。这种特性使其在深度网络中通常能获得更好的训练效果和泛化能力，尤其适合我们这种需要精确拟合复杂函数映射的任务。

3. 损失函数与训练技巧：我们使用Focal Loss，而非标准的均方误差。

   • 为什么用Focal Loss？在我们的数据中，激波（不连续区域）是稀有事件，但恰恰是模拟的关键和难点。标准MSE损失会被大量平滑区域的数据所主导，导致网络对激波附近的通量预测不准。Focal Loss通过引入一个聚焦参数γ，降低那些预测误差小的样本（平滑区）在损失函数中的权重，迫使网络更专注于学习那些难以预测的样本（激波附近）。这相当于在训练中给“难题”加了更高的分数。

2.4 安全卫士：整体凸限制

神经网络学到的映射关系，本质是一个复杂的非线性函数。在绝大多数情况下它工作良好，但在极端情形下（如非常陡的梯度），它可能输出一个“过猛”的修正值，导致下一步计算出的水深变为负数，或者产生数值振荡。

整体凸限制就是用来防止这种情况的“安全阀”。它的核心思想基于一个深刻的数学原理：对于双曲守恒律，其真理解具有一些不变性（如水深非负）。MCL算法通过构造一个“限制器”，对神经网络输出的通量修正值进行“裁剪”。

具体操作如下：

1. 首先，我们计算一个低阶的通量（这里就是LLF格式的非线性部分）和一个由神经网络给出的“高阶”通量修正。
2. MCL算法会根据当前网格单元及其邻居的状态，计算出一组物理上允许的局部边界（例如，水深的上下界、速度的上下界）。
3. 然后，它检查神经网络给出的通量修正值是否会导致下一时间步的预测解超出这些边界。
4. 如果会超出，MCL就“温和地”将修正值拉回到边界允许的最大或最小值内；如果不会，则完全采用神经网络的输出。

这个过程确保了无论神经网络输出什么，最终用于推进计算的通量都是“物理可接受的”。更重要的是，MCL的施加是局部的，且与时间步长满足CFL条件时，可以数学上严格证明数值解保持在容许集内。这就将数据驱动的灵活性与物理约束的严谨性完美结合了起来。

2.5 方案流程图解

整个方案的执行流程可以概括为以下步骤：

1. 数据生成：在精细网格（Nf=1024）上运行DNS模拟，生成长时间序列数据。
2. 数据预处理：将细网格数据空间平均到粗网格（Nc=128），并计算出每个粗网格界面处“真实的”亚网格非线性通量作为训练目标。
3. 训练神经网络：以四点粗网格状态为输入，对应的亚网格通量为输出，使用Focal Loss训练前馈神经网络。
4. 构建粗网格模型：在粗网格上，将LLF通量（包含基础的数值粘性）、神经网络预测的非线性通量修正（经MCL限制）组合起来，形成封闭的粗网格方程。
5. 预测与验证：用训练好的模型（无需重新训练）在不同分辨率的粗网格（Nc=64, 256, 512）上进行模拟，并与DNS结果和纯LLF粗网格模拟结果对比，评估其精度和泛化能力。

3. 关键实现细节与实操要点

3.1 训练数据生成：质量决定上限

神经网络的性能天花板由训练数据决定。我们的数据生成策略遵循了物理模拟和数据科学结合的最佳实践：

• 高保真源：使用足够精细的网格（Nf=1024）和稳健的数值格式（LLF+Heun）进行DNS模拟，确保源数据本身尽可能准确，数值耗散极小。我们将此结果视为“地面真值”。
• 丰富的初始场：为了避免模型只学会某种特定流态，我们随机生成了100条不同的模拟轨迹。初始水深场由平均水深加上两个随机相位正弦波的扰动构成，初始流速也在一定范围内随机。这保证了训练数据能覆盖解空间中的一个有代表性的区域。
• 长期统计平衡：每条轨迹模拟总时长达T=400（无量纲时间），并舍弃初始瞬态过程的数据，主要采集系统达到统计稳态后的数据。这对于学习湍流的长期统计特性至关重要。
• 时间采样策略：模拟时间步长是∆t=0.01，但采样间隔是∆t_sample=0.2。这意味着我们不是每步都采样，而是跳着采。这样做有两个好处：一是减少相邻样本间的相关性，让数据更“独立”；二是降低数据存储和训练时的内存开销。注意，我们只采集Heun方法完成完整步进后的数据，不采集中间预测步的数据，以保证数据的一致性。

实操心得：生成训练数据是耗时最长的步骤，但绝不能偷工减料。初始条件的多样性、模拟时长的充足性，直接决定了学到的参数化模型的泛化能力。在实际项目中，需要根据具体问题调整随机参数的范围，确保覆盖目标应用场景可能出现的所有流态。

3.2 平滑性指示器与数据过滤：聚焦关键区域

在湍流模拟中，大部分区域是相对平滑的，只有少数位置存在陡峭梯度或激波。如果直接用所有数据训练，网络会主要学习平滑区域简单的映射关系，而对激波区域的拟合效果差。

为此，我们引入了平滑性指示器，其定义如下：β_i = (13/12) * (u_{i-1} - 2u_i + u_{i+1})^2 + (1/4) * (u_i - u_{i+1})^2这个公式计算了围绕第i个点的局部解的“粗糙度”。第一项近似二阶导数，捕捉解的曲率；第二项近似一阶导数的平方。β值越大，说明该区域变化越剧烈，越可能存在间断。

我们计算整个数据集中所有样本β值的分位数（例如0.6和0.8分位数），然后只保留β值落在这个区间的数据用于训练。这相当于一个数据过滤器，它：

1. 剔除了大量过于平滑、信息量低的样本，加快了训练速度。
2. 保留了所有激波区域和部分中等变化区域的样本，确保了模型对关键难点的学习。
3. 避免了使用β值最高的极端样本，这些样本可能对应罕见的强激波，数据量少且可能包含数值误差，强行学习容易导致过拟合。

这种方法本质上是困难样本挖掘的一种形式，与Focal Loss的思想一脉相承，从数据选择层面就强化了对关键区域的学习。

3.3 神经网络训练配置：稳定与高效

• 优化器：使用Adam优化器，学习率设为0.001。Adam结合了动量和自适应学习率的优点，在大多数深度学习任务中都是默认的稳健选择。
• 正则化与早停：设置早停耐心为100轮，最大训练轮数为2000轮。这意味着如果连续100个epoch在验证集上的损失没有下降，训练就自动终止。这是防止过拟合最简单有效的方法。我们并没有使用Dropout或L2正则化，因为对于当前规模的网络和数据集，早停通常已足够。
• 硬件与框架：训练在GPU集群上进行，使用PyTorch框架，并借助skorch库（一个将PyTorch模型包装成scikit-learn风格接口的库）来管理训练流程，这方便了数据集划分、验证等操作。
• 输入标准化：一个至关重要的预处理步骤是对输入数据进行标准化。我们将训练集中每个输入特征（8个分量）减去其均值，除以其标准差，将其变换为均值为0、标准差为1的分布。必须注意：在后续使用模型进行预测时，对新的粗网格状态进行相同的变换（使用训练集计算得到的均值和标准差）。这能加速训练收敛，并提高数值稳定性。

3.4 MCL限制器的具体实施步骤

将MCL整合到神经网络粗网格模型中的过程需要细致操作：

1. 计算低阶通量与神经网络通量：

   • 对于粗网格界面I+1/2，计算低阶LLF通量的非线性部分：¯G_{I+1/2} = (f(U_I) + f(U_{I+1})) / 2。
   • 将相邻四点状态输入神经网络，得到原���的高阶通量修正预测：G_{I+1/2}^{NN}。
   • 计算原始修正量：ΔG_{I+1/2} = G_{I+1/2}^{NN} - ¯G_{I+1/2}。这个ΔG就是需要被限制的对象。

2. 计算局部极值：

   • 计算界面两侧的LLF平均状态（bar states）：¯u_{I+1/2}和¯u_{I-1/2}。
   • 从这两个平均状态中，提取水深¯h和速度¯v。
   • 定义当前网格单元I的局部极值：
     • h_min^I = min(¯h_{I-1/2}, ¯h_{I+1/2})
     • h_max^I = max(¯h_{I-1/2}, ¯h_{I+1/2})
     • v_min^I = min(¯v_{I-1/2}, ¯v_{I+1/2})
     • v_max^I = max(¯v_{I-1/2}, ¯v_{I+1/2})
   • 同理定义网格单元I+1的局部极值。

3. 应用限制器：

   • 对水深通量修正的限制：计算允许的修正边界。
     • 上界：UB_h = λ_{I+1/2} * min(¯h_{I+1/2} - h_min^I, h_max^{I+1} - ¯h_{I+1/2})
     • 下界：LB_h = λ_{I+1/2} * max(¯h_{I+1/2} - h_max^I, h_min^{I+1} - ¯h_{I+1/2})
     • 最终被限制的水深通量修正：ΔG_{I+1/2}^{h, *} = min(max(ΔG_{I+1/2}^h, LB_h), UB_h)（实际上需根据正负号分别处理，见正文公式，此处为简化说明）。
   • 对流量通量修正的限制：过程类似，但需基于已限制的水深修正ΔG^{h,*}计算出的新平均状态¯u^{*}，再对流量修正的差值部分进行限制，确保动量通量也满足局部极值原理。

4. 组合最终通量：将限制后的修正通量ΔG^{*}加回低阶通量，并与粘性通量部分G^{visc}组合，得到用于粗网格方程推进的最终界面通量。

注意事项：MCL算法的实现需要仔细处理边界条件（周期边界下索引的循环）。同时，计算局部波速λ_{I+1/2}时，使用的是粗网格上的平均速度和水深。整个限制过程在每个时间步、每个网格界面都需要执行，但由于完全是局部操作，并行效率极高。

4. 数值结果分析与模型评估

我们通过一系列数值实验，从多个维度评估了所提出的神经网络结合MCL参数化方案的有效性。

4.1 能量谱分析：捕捉跨尺度能量传递

能量谱是分析湍流模拟质量的金标准。它描述了动能（或势能）在不同空间尺度（对应波数k）上的分布。一个理想的参数化方案，应该能在粗网格上复现精细网格模拟的能量谱，特别是在代表惯性子区的中间波数段。

我们对比了四种模拟结果：

1. DNS：精细网格（Nf=1024）参考解。
2. LLF-Nc：直接在粗网格（Nc=64,128,256,512）上用LLF格式模拟。这是传统的“无参数化”粗网格方案。
3. NN-Nc：使用我们训练的神经网络参数化（无MCL限制）的粗网格模型。
4. NN-MCL-Nc：使用神经网络参数化并施加MCL限制的粗网格模型。

核心发现：

• 显著优于传统粗网格方案：在所有测试的粗网格分辨率下，NN和NN-MCL模型的能量谱都比同分辨率的LLF模型更接近DNS谱。LLF模型由于数值耗散过大，在高波数（小尺度）部分能量衰减过快，无法维持正确的能量级串过程。
• 分辨率越高，匹配越好：在训练所用的分辨率（Nc=128）上，NN-128模型能较好地复现DNS谱直到波数k≈13。当使用更细的粗网格（Nc=256）时，NN-256模型复现的惯性子区范围扩展到k≈37。这说明我们的参数化模型具有良好的分辨率外推能力。即使网络是在Nc=128的数据上训练的，将其直接应用于Nc=256的网格，依然能有效工作。
• MCL不影响统计特性：对比NN-128和NN-MCL-128的能量谱，可以发现两者几乎完全重合。这是一个非常重要的结论：MCL限制器在保证了数值解物理可接受性的同时，并没有扭曲其长期的统计平衡特性。它主要作用于局部瞬时的激波结构，而对全局的统计能量分布影响微乎其微。

下表总结了不同模型在复现DNS能量谱惯性子区范围方面的能力：

4.2 瞬时解结构：激波处理与MCL的价值

能量谱衡量统计性能，而瞬时解则检验模型对具体物理结构的捕捉能力，尤其是在激波附近。

• 无MCL的潜在风险：在个别测试案例中，当初始条件包含非常强的间断时，纯NN模型（未加MCL）在激波附近偶尔会产生微小的数值过冲或下冲，虽然未必导致计算崩溃，但破坏了解的单调性，这在物理上是不可接受的。
• MCL的稳定化作用：引入MCL后，这些非物理的振荡被完全消除。激波被清晰地捕捉，且解在激波前后保持单调变化，水深严格非负。MCL像一把“软剪刀”，只在神经网络输出可能引发问题时进行干预，在平滑区域则完全信任网络的预测。
• 对光滑区无影响：在解的光滑区域，NN-MCL的解与NN的解几乎没有区别。这证明了MCL的限制是“最小侵入式”的，它不损害神经网络在大部分区域已经学到的精度。

4.3 泛化能力测试：超越训练集

一个强大的参数化模型必须能够处理训练数据未曾覆盖的场景。我们设计了严格的泛化测试：

1. 更强的随机强迫：将强迫项的振幅A从训练时的0.1增加到0.15甚至0.2。这相当于给流体系统输入了更多的能量。
2. 加入地形与摩擦：在方程中引入底部地形变化和曼宁摩擦项，这是实际海洋和河流模拟中必须考虑的物理过程，而训练数据是在平坦无摩擦底床上生成的。

测试结果表明，我们的NN-MCL模型在这些“陌生”的动力学体制下，依然能产生稳定、物理合理的模拟结果。虽然能量谱的细节与DNS会有偏差，但其整体形态、激波捕捉能力都显著优于同分辨率的纯LLF模拟。这证明了四点模板的局部参数化具有捕捉底层物理机制的能力，而不仅仅是记忆训练数据中的模式。

4.4 计算效率分析

虽然训练神经网络需要离线的成本（生成DNS数据、训练网络），但一旦训练完成，在线推理的成本极低。

• 在线成本：在粗网格模拟的每个时间步，每个界面只需调用一次小型神经网络的前向传播（8输入→2输出），以及执行一次本地的MCL限制计算。与一次LLF通量计算相比，增加的计算开销很小。相比于将网格加密8倍（从Nc=128到Nf=1024）所带来的计算量增长，神经网络的额外开销几乎可以忽略不计。
• 离线成本可控：数据生成和训练确实需要计算资源，但这属于“一次投入，长期受益”。并且，由于我们的参数化是局部的，如果需要针对特定区域（如海峡、复杂海岸线）优化，可以只在该区域生成高分辨率数据并重新训练局部网络，无需重做全球模拟，这大大提升了方法的灵活性。

5. 常见问题、挑战与应对策略

在实际实现和应用这套方法时，会遇到一些典型问题。以下是我在复现和研究过程中总结的经验和排查思路。

5.1 训练不收敛或效果差

• 问题表现：训练损失震荡不降，或验证损失远高于训练损失（过拟合），或模型在测试集上预测通量误差巨大。
• 排查清单：
  1. 数据质量：首先检查DNS数据本身是否可靠。确保精细网格模拟是收敛的，时间步长满足CFL条件，数值耗散足够小。可以对比不同分辨率DNS的结果，验证其一致性。
  2. 数据预处理：是否进行了正确的输入标准化？务必使用训练集的均值和标准差来标准化验证集和测试集，这是一个常见的错误源。
  3. 输入输出定义：确认四点模板的索引顺序是否正确，输入输出的物理量（水深、流量）单位是否一致。一个简单的检查方法是：选取一个已知的平滑解析解，手动计算粗网格平均和界面通量，看与网络输出是否量级相符。
  4. 损失函数：Focal Loss中的γ参数可能需调整。如果激波数据非常稀少，可以尝试增大γ（如3或4），进一步聚焦困难样本。同时监控训练过程中“困难样本”（损失大的样本）的比例变化。
  5. 平滑性过滤：检查β分位数的选择是否合理。画出β值的分布直方图，确保过滤后的数据既包含了足够的激波样本，又没有包含过多噪声极大的异常点。可以尝试不同的分位数（如[0.5, 0.9]）进行对比实验。
  6. 网络架构：虽然128x3的隐藏层对当前问题是够用的，但如果问题更复杂，可以尝试增加深度或宽度。同时，可以尝试将平滑性指示器β作为一个额外的输入特征（共9维输入），显式地告知网络当前位置的“粗糙度”，有时能提升性能。

5.2 MCL限制器过于“激进”或“保守”

• 问题表现：“激进”指MCL几乎总是触发，导致神经网络的作用被严重削弱，结果退化成接近纯LLF格式。“保守”指MCL几乎从不触发，但在强激波下依然出现非物理解。
• 排查与调整：
  1. 检查局部波速λ：λ的计算是否正确？它应该是基于粗网格状态U_I和U_{I+1}计算的最大特征速度。如果λ计算偏小，会导致允许的修正边界UB/LB变窄，MCL更容易被触发。
  2. 验证CFL条件：确保粗网格模拟使用的时间步长∆t满足CFL条件：∆t ≤ CFL * ∆X / max(λ)，其中CFL数通常取0.5左右。过大的时间步长会破坏MCL所依赖的凸组合性质，导致限制失效。
  3. 分析神经网络输出：在激波区域，输出ΔG的数值是否量级异常大？如果是，可能是训练数据中激波样本不足或网络未学好。可以针对性增加激波丰富的训练数据，或调整Focal Loss参数。
  4. 可视化限制过程：在调试阶段，可以输出每个界面处ΔG的原值、MCL计算出的边界值以及限制后的值。观察在哪些位置、多大程度上发生了限制。这有助于理解MCL的行为。

5.3 模型泛化到更高维或更复杂方程

• 挑战：本文工作聚焦于一维浅水方程。如何扩展到二维SWE或更复杂的方程（如三维Navier-Stokes方程）？
• 思路与策略：
  1. 模板扩展：在二维中，四点模板需扩展为九点模板（中心点及其八个邻居），输入维度会显著增加（例如，每个点有h, hu, hv三个变量，9*3=27维）。网络规模和训练数据需求会增长。
  2. 卷积神经网络：对于二维或三维问题，CNN是更自然的选择。它可以利用图像的局部空间相关性，参数共享也能减少参数量。输入可以是粗网格场的一个局部图像块。
  3. 物理约束的融入：对于更复杂的方程，MCL的限制条件需要重新推导。核心是找到该方程组的物理容许集（如密度、压力为正，速度有界等），并设计相应的限制算法。这可能涉及更复杂的特征分解和局部极值传播。
  4. 分步策略：不要试图一步到位。可以从二维无旋转SWE开始，再加入科氏力，然后考虑非静压近似等。每一步都需仔细验证能量守恒、位涡守恒等关键物理性质。

5.4 与其它机器学习参数化方法的对比

• 对比传统回归方法：我们尝试过用多项式回归、随机森林等方法来拟合同样的四点通量映射关系，但效果均不如神经网络。神经网络在捕捉高度非线性关系上的能力是这些方法无法比拟的。
• 对比全局参数化：有些方法使用一个巨大的神经网络，以整个计算域的粗网格场作为输入，输出所有界面的通量。这种“全局耦合”参数化虽然理论上容量更大，但计算成本高，难以并行，且泛化到不同网格分辨率或域大小时非常困难。我们的局部四点参数化在这些方面具有明显优势。
• 对比基于方程发现的方法：有些工作尝试用符号回归等方法直接从数据中发现参数化的解析形式。这类方法可解释性强，但通常只能发现相对简单的公式，对于湍流中复杂的非局部、非线性效应，目前还难以与深度神经网络竞争。

最后，我个人在实践中的体会是，成功的关键在于“平衡”。平衡数据驱动与物理约束，平衡模型复杂度与泛化能力，平衡离线训练成本与在线推理效率。我们的“神经网络+MCL”框架提供了一个很好的平衡点。它不是一个用AI完全取代物理的黑箱，而是一个让AI在严格物理规则辅助下，去解决那些传统方法难以处理的“灰色地带”问题的有力工具。对于从事高分辨率模拟受限领域的研究者和工程师来说，这种混合方法无疑是值得深入探索和适配的方向。