面向复杂地质构造中CO2迁移的快速GNN替代模型

# 面向复杂地质构造中CO₂运移的快速GNN代理模型

来源：https://arxiv.org/html/2606.17180

11institutetext:Rodrigo S\. Luna22institutetext:系统工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；22email:luna@cos\.ufrj\.br33institutetext:Thiago H\. N\. Coelho44institutetext:系统工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；44email:tcoelho@cos\.ufrj\.br55institutetext:Luiz S\. L\. Neto66institutetext:土木工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；66email:luizlealn@nacad\.ufrj\.br77institutetext:Roberto M\. Velho88institutetext:系统工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；88email:roberto\.velho@nacad\.ufrj\.br99institutetext:Adriano M\. A\. Cortes1010institutetext:系统工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；1010email:adriano@nacad\.ufrj\.br1111institutetext:Renato N\. Elias1212institutetext:土木工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；1212email:rnelias@nacad\.ufrj\.br1313institutetext:Alexandre G\. Evsukoff1414institutetext:土木工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；1414email:alexandre\.evsukoff@coc\.ufrj\.br1515institutetext:Fernando A\. Rochinha1616institutetext:机械工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；1616email:faro@mecanica\.coppe\.ufrj\.br1717institutetext:Mauricio Araya\-Polo1818institutetext:道达尔能源公司1818email:mauricio\.araya@totalenergies\.com1919institutetext:Herve Gross2020institutetext:道达尔能源公司2020email:herve\.gross@totalenergies\.com2121institutetext:Alvaro L\. G\. A\. Coutinho2222institutetext:土木工程与高性能计算中心，NACAD \- COPPE，里约热内卢联邦大学；2222email:alvaro@nacad\.ufrj\.brThiago H\. N\. CoelhoLuiz S\. L\. NetoRoberto M\. VelhoAdriano M\. A\. CortesRenato N\. EliasAlexandre G\. EvsukoffFernando A\. RochinhaMauricio Araya\-PoloHerve GrossAlvaro L\. G\. A\. Coutinho

## 1 引言

准确高效地预测复杂地质构造中CO₂羽流运移，对于安全部署碳捕集与封存（CCS）技术至关重要。虽然传统数值储层模拟器在物理上可靠，但在长期模拟和不确定性量化方面会带来高昂的计算成本。近期的科学机器学习方法，包括傅里叶神经算子、DeepONet和图神经网络（GNN）[Juet al. (2024)](https://arxiv.org/html/2606.17180#bib.bib5)，在多查询场景（如参数探索和情景分析）中表现出高精度和显著的加速效果，并越来越多地被应用于数字孪生中。GNN特别适用于地下应用，因为它们天然地适用于地质储层典型的非结构化网格。然而，大多数基于图的模拟器（包括MeshGraphNet [Pfaffet al. (2021)](https://arxiv.org/html/2606.17180#bib.bib8)）依赖于各向同性消息传递，即均匀地向相邻单元传播信息，这对于受网格几何、渗透率对比、重力效应和地质非均质性主导的强方向性输运的储层流动问题可能具有局限性。在本工作中，我们开发了一种先进的端到端GNN管道用于CO₂羽流运移预测，该管道引入了一种各向异性消息传递机制，该机制由几何和物理边属性显式调制，同时结合了潜在时间动态公式和自回归训练策略，以确保稳定的长期预测。所提出的方法在SPE11A基准上进行了评估，这是一个受控CO₂封存实验启发的二维实验室规模测试案例 [Landa-Marbán and Sandve (2025)](https://arxiv.org/html/2606.17180#bib.bib11); [Rasmussenet al.(2021)](https://arxiv.org/html/2606.17180#bib.bib14); [Nordbottenet al.(2024)](https://arxiv.org/html/2606.17180#bib.bib15)，其控制方程为由达西流、组分质量守恒和相平衡约束耦合的全非线性偏微分方程组。由于尖锐的气水界面以及伴随广泛指进的快速对流混合，该基准带来了显著挑战，使得气体饱和度和相组成的准确预测尤为困难。本文其余部分组织如下。第[2节](https://arxiv.org/html/2606.17180#Ch0.S2)简要描述了GNN引擎，该引擎基于[Juet al.(2024)](https://arxiv.org/html/2606.17180#bib.bib5)中提出的引擎构建。在第[3节](https://arxiv.org/html/2606.17180#Ch0.S3)中，我们描述了实验设置，包括考虑随机地质变化的SPE11A数据生成、数据准备、模型与训练设置、评估协议以及预测结果。本章最后总结了我们的主要结论。

## 2 基于图的预测模型

我们的预测模型基于MeshGraphNet [Pfaffet al. (2021)](https://arxiv.org/html/2606.17180#bib.bib8)，这是一个用于基于网格的物理系统的编码器-处理器-解码器图神经模拟器。在其原始公式中，MeshGraphNet采用各向同性消息传递，并使用单步监督进行训练，然后以自回归方式部署，这可能导致长期储层流动模拟中的误差累积 [Juet al.(2024)](https://arxiv.org/html/2606.17180#bib.bib5)。

储层流动由于网格几何、渗透率对比和地质非均质性而表现出强方向性输运。为了更好地捕捉这些效应，我们引入了受几何各向异性公式启发的各向异性消息传递 [Thürlemann and Riniker (2023)](https://arxiv.org/html/2606.17180#bib.bib4)。我们不采用张量多极表示，而是采用一种轻量级的定向边加权机制，其中相互作用权重由几何特征（如相对单元位置和边属性）计算得出，使模型能够优先考虑物理相关的输运方向。

此外，为了捕获超出单步转换的时间依赖性并提高长期预测的稳定性，我们结合了一个基于GraphConv–LSTM的循环潜在动态模块以及一个残差状态更新公式。这些扩展定义了所提出的*AnisoMeshGraph-LSTM*架构。

### 2.1 SPE11A基准的图表示

我们将SPE11A基准重新表述为图G=(V,E)，其中每个计算单元表示为一个节点，具有非零传导率的单元之间通过边连接。单元i和j之间每个连接的传导率根据网格几何和岩石渗透率计算，当两个单元被断层分隔时为零 [Juet al.(2024)](https://arxiv.org/html/2606.17180#bib.bib5)。对于给定的网格和渗透率，图结构在所有时间步保持不变。与节点i在时间t相关的属性称为节点特征。与边(i,j)相关的属性称为边特征，并且与时间无关。

节点特征编码每个单元处系统的物理状态，包括时间相关变量和外生输入，而边属性捕获几何和物理相互作用属性。虽然传导率用于定义图拓扑，但并不直接作为边特征包含在内；相反，单元中心之间的位移向量和欧几里得距离作为几何代理，隐式编码了连接结构。

### 2.2 编码器

令X^t ∈ R^(n×F_v) 表示系统在时间t的物理状态，定义在所有图节点上。每一行x_i^t 包含与节点i相关的时间相关状态变量（例如饱和度或密度），以及时间不变参数和控制输入。模拟器提供轨迹 {X^t}_{t=1}^N，该轨迹既作为模型的输入，也作为监督训练中的真实值。

边属性收集在矩阵E_attr ∈ R^(m×F_e) 中，其中每一行e_ij 编码节点i和j之间界面的几何和物理属性。节点和边特征使用多层感知器（MLP）嵌入到潜在空间中：

z_i^{t,0} = φ_v^0(x_i^t) ∈ R^{d_v}, \qquad z_ij^{t,0} = φ_e^0(e_ij) ∈ R^{d_e}, \tag{1}

其中 d_v 和 d_e 分别表示节点和边的潜在维度。φ_v^0 和 φ_e^0 都是带有残差连接的MLP。

### 2.3 处理器：各向异性消息传递

处理器应用L层消息传递。在每一层 l=1,...,L，边和节点的嵌入依次更新。为了捕捉由储层连通性和地质非均质性引发的方向性交互，我们采用具有可学习的、依赖于边的交互权重的各向异性消息传递。

##### 边更新

z_ij^{t,l} = φ_e^l( [z_ij^{t,l-1}, z_i^{t,l-1}, z_j^{t,l-1}] ), \tag{2}

其中 φ_e^l 是一个带有残差连接和ReLU激活函数的MLP。

##### 各向异性权重

对于每条有向边 (i,j)，计算一个标量兼容性分数：

s_ij^{t,l} = φ_α^l( [z_ij^{t,l}, z_i^{t,l}, z_j^{t,l}] ), \tag{3}

并在节点i的邻域上进行归一化：

α_ij^{t,l} = \frac{\exp(s_ij^{t,l})}{\sum_{k∈N(i)} \exp(s_ik^{t,l})}. \tag{4}

##### 节点更新

z_i^{t,l} = φ_v^l( [z_i^{t,l-1}, \sum_{j∈N(i)} α_ij^{t,l} z_ij^{t,l}] ), \tag{5}

其中 φ_v^l 是一个带有残差连接和ReLU激活函数的MLP。N(i) 表示节点v_i的邻域，[·] 表示沿特征维度的拼接。

尽管加权机制在结构上与边条件注意力相关，但其不同之处在于，式(3)中的兼容性分数以边嵌入 z_ij 为条件，而 z_ij 在第0层仅由几何特征（位移向量和欧几里得距离）初始化，且不包含任何语义节点特征信息。这种几何初始化使所学权重偏向于物理上合理的输运方向，这与标准注意力机制（其分数来源于无约束的节点表示）不同。

### 2.4 用于潜在时间动态的GraphConv–LSTM

为了减缓多个展开步上的误差累积，在每层消息传递块之后集成一个LSTM单元，使模型能够捕获并传播潜在节点嵌入中的时间依赖性。具体来说，我们遵循Seo等人[Seoet al.(2018)](https://arxiv.org/html/2606.17180#bib.bib9)提出的图卷积LSTM公式，其中ConvLSTM中的标准卷积被谱图卷积算子 *G 取代。

令 Z^{t,L} ∈ R^(n×d_v) 为从消息传递模块最后一个步骤L和时间t获得的节点嵌入矩阵，并令 H^{t-1}、C^{t-1} 为前一步LSTM的隐藏状态和细胞状态。GraphConv-LSTM单元的更新定义为：

i^t = σ(W_xi *G Z^{t,L} + W_hi *G H^{t-1} + b_i), \tag{6}
f^t = σ(W_xf *G Z^{t,L} + W_hf *G H^{t-1} + b_f),
C^t = f^t ⊙ C^{t-1} + i^t ⊙ tanh(W_xc *G Z^{t,L} + W_hc *G H^{t-1} + b_c),
o^t = σ(W_xo *G Z^{t,L} + W_ho *G H^{t-1} + W_c ⊙ C^t + b_o),
H^t = o^t ⊙ tanh(C^t).

这里，i^t、f^t 和 o^t 分别表示时间t的输入门、遗忘门和输出门；⊙ 是逐元素（Hadamard）乘积。每个门对其仿射变换应用sigmoid激活函数σ(·)，该变换由可学习的权重矩阵 W_* 和偏置向量 b_* 参数化。图卷积算子 *G 使用切比雪夫多项式实现 [Defferrardet al.(2016)](https://arxiv.org/html/2606.17180#bib.bib10)。

在L层消息传递之后，节点嵌入收集在 Z^t ∈ R^(n×d_v) 中。这些嵌入由GraphConv–LSTM处理，该LSTM通过图卷积捕获时间依赖性，同时保持空间耦合。该模块输出隐藏状态 H^t ∈ R^(n×d_v)。

### 2.5 解码器与时间积分

消息传递之后，最终的节点嵌入由GraphConv-LSTM处理，该LSTM直接在潜在空间中捕获时间依赖性，并产生用于预测的隐藏状态。与[Juet al.(2024)](https://arxiv.org/html/2606.17180#bib.bib5)直接预测下一状态不同，我们采用残差公式，其中解码器预测物理状态的增量而非绝对值。动态状态使用显式欧拉更新推进，如[Eliasofet al.(2024)](https://arxiv.org/html/2606.17180#bib.bib6)所提出：

X_dyn^{t+1} = X_dyn^t + Δt φ_dec(H^t), \qquad H^t = GraphConvLSTM(Z^t). \tag{7}

### 2.6 训练目标

模型使用自回归多步监督进行训练，预测视界长度为T。给定模拟器生成的轨迹 {Ŷ^t}_{t=1}^T，损失定义为所有节点和展开步上预测状态与真实状态之间的均方误差：

L(θ) = \frac{1}{T - T_train} \sum_{t=1}^{T - T_train} \frac{1}{T_train} \sum_{j=1}^{T_train} \frac{1}{n} \| Ŷ^{t+j} - X^{t+j} \|_2^2.