Hydrate reservoir fluid-solid output multi-field coupling prediction model construction method
By constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs and employing unstructured grids and deep learning methods, the problems of insufficient adaptability and accuracy in hydrate reservoir prediction were solved, and high-precision long-term prediction was achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- GUANGZHOU NANSHADI BINHAI RESEARCH INSTITUTE
- Filing Date
- 2026-01-30
- Publication Date
- 2026-05-08
AI Technical Summary
Existing deep learning methods are difficult to flexibly adapt to the irregular boundaries of hydrate reservoirs. The lack of physical constraints leads to insufficient prediction accuracy and consistency, and long-term predictions often produce non-monotonic output errors that violate physical facts.
A coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction, and sediment constitutive relations was constructed. The reservoir was discretized using an unstructured grid, and the evolution data was generated by an alternating iterative algorithm. A graph neural network and a long short-term memory network were constructed, and features that integrate physical properties and spatial structure were extracted through graph convolutional layers and gating mechanisms. A loss function of prediction error and physical constraint terms was introduced for training.
It significantly improves the physical consistency and generalization accuracy of long-term prediction of heterogeneous reservoirs, and achieves high-precision full-dimensional evolution prediction, while taking into account both mechanism fidelity and computational real-time performance.
Smart Images
Figure CN121996974A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of intelligent hydrate mining technology, and in particular to a method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs. Background Technology
[0002] Natural gas hydrates, as a potentially high-potential clean alternative energy source, involve the extraction of a complex interplay of multiple physical mechanisms, including phase transformation decomposition, multiphase flow, heat conduction, and the mechanical response of sedimentary skeletons. Accurate prediction of cumulative gas production and sand production throughout the entire extraction cycle is crucial for assessing production potential, developing scientific extraction plans, and preventing wellbore geological hazards. With the development of artificial intelligence, utilizing deep learning models to mine complex reservoir evolution data has become a mainstream technological trend, replacing traditional time-consuming numerical simulations and enabling efficient real-time production prediction.
[0003] However, existing deep learning methods are mostly limited by regular grid architectures, making it difficult to flexibly adapt to irregular reservoir boundaries and local encryption requirements. Furthermore, pure data-driven models lack physical conservation laws, often producing non-monotonic output errors that violate physical facts in long-term predictions, resulting in insufficient model generalization ability and physical consistency. Summary of the Invention
[0004] To overcome the above shortcomings, this invention provides a method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs, aiming to improve the problems of existing methods being difficult to adapt to unstructured grids and having insufficient prediction accuracy and consistency due to the lack of physical constraints.
[0005] This invention provides the following technical solution: a method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs, comprising:
[0006] S1. Establish a coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction and sediment constitutive relations, and use unstructured grids to discretize the reservoir;
[0007] S2. Solve the coupled mathematical model using an alternating iterative algorithm to generate evolution data of pore pressure, temperature, phase saturation, solid displacement, and corresponding cumulative gas production and cumulative sand production of the grid nodes of the unstructured grid at continuous time steps.
[0008] S3. Construct a topology graph based on the unstructured mesh, map the mesh nodes to graph nodes, map the geometric adjacency relationships between the mesh nodes to graph edges, and assign the pore pressure, temperature, saturation of each phase, and solid displacement as initial node features.
[0009] S4. Construct a graph neural network, and perform message passing and aggregation on the initial node features through graph convolutional layers to extract local embedded features that fuse physical properties and spatial structure.
[0010] S5. Construct a long short-term memory network to receive the time series of the local embedded features and generate a temporal latent vector representing the fluid-solid transport state through a gating mechanism;
[0011] S6. Using a loss function containing prediction error and physical constraint terms, jointly train the graph neural network and the long short-term memory network to establish a nonlinear mapping between the time-series hidden vector and the cumulative gas production and cumulative sand production, and generate a prediction model.
[0012] Preferably, in step S1, the step of establishing a coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction, and sediment constitutive relations includes:
[0013] Define the hydrate decomposition rate affected by saturation and contact area, and establish a mass source term expression for the solid-to-fluid phase transition.
[0014] By combining relative permeability and capillary pressure function, a set of partial differential equations for mass conservation of multiphase fluids containing gravity and pressure gradient is constructed.
[0015] An elastoplastic mechanical model is used to describe the stress response of the skeleton, and a fluid-solid output flux model is constructed using a joint criterion of plastic strain and flow velocity.
[0016] Preferably, in step S1, the step of discretizing the reservoir using an unstructured grid includes:
[0017] The reservoir geometry is discretized using a spatial partitioning algorithm, and local mesh refinement is performed in regions where the physical field gradient exceeds a preset threshold.
[0018] Based on geological logging data, an interpolation algorithm is used to generate a spatial attribute field, and continuous geological attributes are projected onto the grid nodes of the unstructured grid.
[0019] The coupled mathematical model is transformed into a system of algebraic equations using a numerical discretization method, and the flux conservation of the unstructured mesh interface is handled.
[0020] Preferably, in step S2, the step of solving the coupled mathematical model using an alternating iterative algorithm includes:
[0021] Set the time step parameters and construct a nonlinear solver that uses the solution from the previous time step as the initial value;
[0022] A segmented iterative strategy is adopted to solve the fluid and temperature fields first, and then update the solid mechanical and chemical field parameters.
[0023] Calculate the residual norm of each physical field between iterations, and output the state vector and advance the time step when the residual satisfies the preset tolerance.
[0024] Preferably, in step S3, the step of constructing a topology graph based on the unstructured mesh includes:
[0025] Using the grid nodes as vertices, the normalized multiphysics state data is combined as the initial node attribute vector;
[0026] Based on the topological connection relationships between the grid nodes, connected edges are established to form an adjacency matrix describing spatial connectivity;
[0027] The weight coefficients are calculated based on the spatial distance and transmission characteristics between adjacent graph nodes to characterize the interaction strength between grids.
[0028] Preferably, in step S4, the step of performing message passing and aggregation on the initial node features through a graph convolutional layer includes:
[0029] Perform a nonlinear transformation on the source node features and edge weights to generate a message vector containing local physical gradients;
[0030] The influence coefficients of neighboring nodes are calculated using an attention mechanism, and the message vectors on the adjacent edges are weighted and aggregated.
[0031] By fusing aggregated features with its own historical features, the network is updated to generate new node embedding vectors that contain spatial topology information.
[0032] Preferably, in step S5, the step of generating the temporal implicit vector representing the fluid-solid transport state through a gating mechanism includes:
[0033] A loop structure containing multiple gating units is configured to process the time series of the local embedded features;
[0034] The gating operator is used to selectively forget the historical hidden states and to perform weighted fusion of the current input features;
[0035] A nonlinear mapping is performed on the updated cell state to output the temporal implicit vector that reflects the evolution of fluid-solid transport.
[0036] Preferably, in step S6, the step of jointly training using a loss function containing prediction error and physical constraint terms includes:
[0037] By comparing the model's predicted values with the corresponding cumulative gas production and cumulative sand production in the evolution data, a basic loss term characterizing the data fitting accuracy is constructed.
[0038] Based on the physical monotonicity characteristics of the cumulative gas production and cumulative sand production, a penalty is imposed on the difference values in the prediction sequence that violate physical laws.
[0039] An integrated objective function is constructed by using adaptive weights to balance the basic loss term and the physical constraint term.
[0040] Preferably, in step S6, the step of establishing the nonlinear mapping between the time-series implicit vector and the cumulative gas production and cumulative sand production includes:
[0041] The temporal latent vector is mapped to a low-dimensional output space using a fully connected network;
[0042] The gradient of the synthetic objective function with respect to the trainable parameters of the model is calculated based on the automatic differentiation framework;
[0043] The weights of the graph neural network and the long short-term memory network are updated synchronously according to the gradient using an adaptive optimization strategy.
[0044] The present invention has the following beneficial effects:
[0045] 1. In this invention, by constructing an unstructured grid topology map to adapt to the complex geometric boundaries of the reservoir, and introducing physical monotonicity constraints on output during training, the problem of poor adaptability of traditional regular grid models and violation of physical conservation by pure data-driven methods is overcome, and the physical consistency and generalization accuracy of long-term prediction of heterogeneous reservoirs are significantly improved.
[0046] 2. In this invention, complete physical evolution data generated by numerical simulation is used to drive a deep learning model. Lightweight neural networks replace time-consuming numerical discrete iterative calculations. While solving the problem of scarce field samples, this invention achieves a leapfrog improvement in both mechanism fidelity and computational real-time performance.
[0047] 3. In this invention, the spatial topology extraction of coupled graph convolutional networks and the temporal deduction capability of long short-term memory networks simultaneously capture the nonlocal correlation of fluid conduction and the long time lag characteristics of phase transition production, solving the accuracy loss caused by the separation of spatiotemporal information and providing high-precision full-dimensional evolution prediction for hydrate mining. Attached Figure Description
[0048] Figure 1 This is a flowchart of the method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs proposed in this invention. Detailed Implementation
[0049] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0050] In embodiments of the present invention, the present invention provides a method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs, such as... Figure 1 As shown, it includes:
[0051] S1. Establish a coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction and sediment constitutive relations, and use unstructured grids to discretize the reservoir;
[0052] Further, in step S1, the step of establishing a coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction, and sediment constitutive relations includes: defining the hydrate decomposition rate affected by saturation and contact area; establishing a mass source term expression for the solid-to-fluid phase transition; constructing a multiphase fluid mass conservation partial differential equation system containing gravity and pressure gradient by combining relative permeability and capillary pressure function; describing the skeleton stress response using an elastoplastic mechanical model; and constructing a fluid-solid output flux model using a joint criterion of plastic strain and flow velocity.
[0053] Further, in step S1, the step of discretizing the reservoir using an unstructured grid includes: using a spatial partitioning algorithm to perform unstructured discretization of the reservoir geometry, and performing local grid refinement in areas where the physical field gradient exceeds a preset threshold; generating a spatial attribute field based on geological logging data using an interpolation algorithm, and projecting continuous geological attributes onto the grid nodes of the unstructured grid; using a numerical discretization method to transform the coupled mathematical model into a system of algebraic equations, and handling the flux conservation at the unstructured grid interface.
[0054] Specifically, a coupled mathematical model describing the physicochemical behavior of hydrate reservoirs is first constructed. This model's core encompasses hydrate decomposition kinetics, gas-water two-phase flow, heat transfer, and the mechanical response of the sediment skeleton. For the hydrate decomposition process, a kinetic-based phase transition model is used to define the decomposition rate, which is jointly controlled by the current pore pressure, temperature, and hydrate saturation. The mass source term for the solid-phase hydrate to fluid phase transition is calculated using a modified Kim-Bishnoi model, and the decomposition rate is... The calculation formula is:
[0055] ;
[0056] In the formula, This indicates the rate of gas production from the decomposition of hydrates per unit volume, expressed in kilograms per cubic meter per second. This is the kinetic constant for hydrate decomposition, which is related to activation energy and temperature; is the molar mass of the hydrate; is the specific surface area of the hydrate particles; The hydrate saturation of the current grid node; This is an empirical index of the effect of saturation on the reaction area; This represents the equilibrium fugacity at the current temperature. This represents the gas fugacity corresponding to the current gas phase pressure. This formula clarifies that when the reservoir pressure is lower than the equilibrium pressure, the hydrate decomposes and becomes a mass source term in the subsequent flow equation.
[0057] Regarding fluid transport, a system of partial differential equations for mass conservation in multiphase fluids, incorporating gravity and pressure gradients, is established. Darcy's law is used to describe the macroscopic flow of gas and water phases in porous media, and relative permeability curves and capillary pressure functions are introduced to characterize multiphase disturbances. For any fluid phase... Take gas phase or aqueous phase Its mass conservation governing equation is constructed as follows:
[0058] ;
[0059] In the formula, Reservoir porosity; The density of the fluid phase; This refers to the fluid phase saturation. For time; For Hamiltonian operators; The absolute permeability tensor of sediments; Relative permeability, saturation The function; For fluid viscosity; For phase pressure; This is the vector of gravitational acceleration; For source and sink terms, it includes the mass of fluid produced by the decomposition of hydrates. And the wellbore production. This set of equations provides a precise conservation description of fluid pressure propagation and mass transport within the reservoir.
[0060] To address reservoir mechanical response and sand production prediction, constitutive relations and sand production criteria are defined. An elastoplastic mechanical model is used to describe the stress-strain relationship of the sedimentary framework under effective stress variations, and the plastic shear strain of the framework is calculated. Based on this, a fluid-solid production flux model is constructed using a joint criterion of plastic strain and fluid velocity. Sand production is only considered to occur when the framework undergoes plastic yielding and the pore flow velocity is sufficient to carry particles. Sand production flux is then calculated. The calculation formula is set as follows:
[0061] ;
[0062] In the formula, For sand output quality throughput; The erosion coefficient is used to characterize the bonding strength of the skeleton. The density of sand particles; The cumulative equivalent plastic strain is obtained from the mechanical constitutive equations; This represents the true velocity modulus of the pore fluid. The critical flow velocity for particle initiation is given. This formula directly correlates the macroscopic sand output with the microscopic skeletal failure and hydrodynamic conditions, achieving a quantitative characterization of fluid-solid production.
[0063] After constructing the mathematical model, the unstructured grid discretization of the reservoir is performed. First, a spatial partitioning algorithm is used to discretize the reservoir geometry into a series of unstructured control volumes. During this process, the physical field gradients in the pre-model are calculated, primarily the pressure gradient. and hydrate saturation gradient Set a preset gradient threshold .exist In areas such as around production wells and at the hydrate decomposition front, the grid node density is automatically increased, and local grid refinement is performed to ensure that numerical calculations can capture drastic physical changes.
[0064] Subsequently, the heterogeneous physical properties are mapped. Geological logging data of the target reservoir is acquired, and the spatial distribution fields of porosity, permeability, and initial saturation are generated using Kriging interpolation or inverse distance weighting algorithms. These continuous geological properties are projected onto each grid node of the unstructured grid to form an initial parameter matrix. Finally, the coupled mathematical model is numerically discretized using the finite volume method (FVM), transforming the partial differential equations into a system of algebraic equations with unknowns at the grid nodes as variables. When dealing with the unstructured grid interface, the normal flux of the control volume interface is calculated using flux reconstruction technology, strictly ensuring the conservation of mass and energy at the irregular grid interface, ultimately forming a matrix system that can be solved iteratively by a computer.
[0065] This implementation step, by constructing a refined physical model that includes phase change, seepage, and mechanical damage, and combining it with highly adaptable unstructured mesh technology, can realistically reproduce the complex mining response of hydrate reservoirs under heterogeneous geological conditions, especially in the high-stress zone around the well and the high-flow zone at the decomposition front. It can simultaneously ensure physical conservation and computational accuracy, providing high-fidelity physical evolution data for subsequent neural network training.
[0066] S2. Solve the coupled mathematical model using an alternating iterative algorithm to generate evolution data of pore pressure, temperature, phase saturation, solid displacement, and corresponding cumulative gas production and cumulative sand production of the grid nodes of the unstructured grid at continuous time steps.
[0067] Furthermore, in step S2, the step of solving the coupled mathematical model using the alternating iterative algorithm includes: setting time step parameters, constructing a nonlinear solver that uses the solution of the previous time step as the initial value; adopting a segmented iterative strategy to solve the fluid and temperature fields first, and then updating the solid mechanics and chemical field parameters; calculating the residual norm of each physical field between iterations, and outputting the state vector and advancing the time step when the residual meets the preset tolerance.
[0068] Specifically, the time stepping strategy for the numerical simulation is first initialized, and the total simulation duration is set. With the initial time step The time partial differential equation established in step S1 is discretized using a fully implicit backward Euler difference scheme to ensure unconditional stability of the numerical computation. When constructing the nonlinear solver, the previous time step is used... The convergent solution vector is used as the current time step. initial guess value This improves the convergence speed of nonlinear iteration, and the converged solution vector includes the pressure, temperature, saturation and displacement of all grid nodes.
[0069] Entering the multiphysics decoupling loop stage, an operator partitioning iterative strategy is employed to handle the strongly coupled equations. In each nonlinear iteration step... In this process, the coupled subsystems of fluid flow and heat conduction are solved first. The Jacobian matrix is constructed. With residual vector The pressure is solved using the Newton-Raphson algorithm. ,temperature With saturation The increment. The iterative correction formula is expressed as:
[0070] ;
[0071] ;
[0072] In the formula, This is the state variable vector of the heat flow subsystem, which includes nodal pore pressure, gas-water saturation, and temperature. The tangent stiffness matrix is composed of the partial derivatives of the governing equations with respect to the state variables. This is the discrete residual vector of the heat flow equations; This is the correction amount for the current iteration step.
[0073] After solving the fluid-thermal field, the effective stress change is calculated using the updated pore pressure, followed by solving the solid mechanics and chemical decomposition subsystems. The updated fluid state is then substituted into the force equilibrium equations to solve for the displacement increments of the mesh nodes. And the kinetic equations for hydrate decomposition. During this process, the stiffness matrix is updated based on the sediment constitutive model. And calculate the nodal displacements:
[0074] ;
[0075] In the formula, is the tangent stiffness matrix of the sediment skeleton; The nodal displacement vector; This represents the increment of the external load. This is the strain-displacement matrix; For Biot coefficients; This is the integration region. Simultaneously, based on the updated temperature and pressure, the hydrate saturation is updated using the phase transition equation. and porosity The current instantaneous sand output rate is calculated based on the fluid-solid production flux model.
[0076] Finally, convergence is determined and the state is advanced. The relative residual norms of the fluid field, temperature field, and mechanical field in the current iteration step are calculated. The convergence criterion is defined as follows: the residual norm of all physical fields is less than a preset tolerance. :
[0077] ;
[0078] In the formula, Represents the infinite norm; This is a reference force or flux vector; The value is usually taken as If the convergence condition is not met, update the iteration step. Repeat the above decoupling loop; if the convergence condition is met, output the current time step. The status data of all grid nodes is collected. Based on the flow velocity and sand concentration at the boundary nodes, the cumulative gas production is calculated through time integration. With cumulative sand output :
[0079] ;
[0080] In the formula, Represents gaseous or solid phase sand; This is the set of nodes at the wellbore outflow boundary; For nodes The instantaneous mass flow rate at the given point. After confirming convergence, the current solution is used as the initial value for the next time step, and the system time is advanced. Until the total simulation time is reached.
[0081] This implementation step transforms a highly nonlinear, multi-field strongly coupled problem into a series of solvable linear subproblems through alternating iteration and operator segmentation techniques. This effectively solves the numerical non-convergence problem caused by hydrate phase transitions and large framework deformations, and can stably generate high spatiotemporal resolution reservoir evolution datasets, providing complete physical samples for the training of subsequent artificial intelligence models.
[0082] S3. Construct a topological graph based on the unstructured mesh, map mesh nodes to graph nodes, map the geometric adjacency relationships between mesh nodes to graph edges, and assign pore pressure, temperature, saturation of each phase and solid displacement as initial node features.
[0083] Furthermore, in step S3, the step of constructing a topological graph based on an unstructured grid includes: using grid nodes as vertices, combining normalized multiphysics state data as initial node attribute vectors; establishing connected edges based on the topological connection relationships between grid nodes to form an adjacency matrix describing spatial connectivity; and calculating weight coefficients based on the spatial distance and conduction characteristics between adjacent graph nodes to characterize the interaction strength between grids.
[0084] Specifically, the physical field data of the unstructured mesh generated in step S2 is first extracted at the current time step to construct the graph structure data for the graph neural network input. This process first defines the vertex set and feature vectors of the graph, directly mapping each mesh node in the unstructured mesh to a vertex in the graph structure. To eliminate the influence of different physical dimensions and their numerical magnitudes on the convergence of subsequent neural network training, the extracted pore pressure, temperature, hydrate saturation, gas phase saturation, water phase saturation, and solid displacement components are subjected to maximum and minimum normalization. The normalization calculation formula is:
[0085] ;
[0086] In the formula, Indicates the first The first grid node Normalized values of physical properties; These are the original physical values; and These represent the minimum and maximum values of the physical attribute across the entire dataset. After normalization, all processed physical attributes of each node are combined to construct the initial state attribute vector for that node. This vector serves as the initial node feature of the input layer of the graph neural network.
[0087] Subsequently, the edge set and adjacency matrix of the graph are constructed. The topological connections of the unstructured mesh are traversed. For any two mesh nodes in 3D space, if they share the same mesh edge in a finite element or finite volume mesh, it is determined that these two nodes have a connection in the graph structure, and an undirected edge is established between the graph vertices. An adjacency matrix is generated based on this connection. If node With nodes If connected, then matrix elements Set to 1 otherwise set to 0. This adjacency matrix fully preserves the geometric connectivity of the unstructured mesh in physical space.
[0088] Finally, edge attribute weights are calculated to characterize the interaction strength between grids. This interaction depends not only on geometric distance but also on reservoir properties such as permeability. The Euclidean distance between adjacent graph nodes and the geometric mean of node permeability are calculated to construct weighting coefficients that incorporate spatial distance attenuation and fluid conductivity. The weight calculation formula is set as follows:
[0089] ;
[0090] In the formula, and They are nodes With nodes Absolute penetration rate at the location; and Let be the spatial coordinate vectors of the two nodes; The Euclidean norm of a vector is used to represent spatial distance. To prevent tiny constants with a denominator of zero; This is a Gaussian kernel bandwidth parameter used to control the degree of distance attenuation. The larger this weighting coefficient is, the closer the two nodes are physically and the stronger their fluid conduction ability, which should result in higher efficiency for subsequent message transmission.
[0091] This implementation step transforms unstructured physical meshes into graph-structured data containing rich physical features and connection weights, solving the problem that traditional convolutional neural networks cannot directly process non-Euclidean spatial data. At the same time, by integrating distance and conductivity weights, the graph structure not only reflects geometric topology but also reflects the physical priority channels of underground fluid flow, laying a data foundation for subsequent graph neural networks to accurately extract fluid-structure interaction features.
[0092] S4. Construct a graph neural network, and use graph convolutional layers to perform message passing and aggregation on the initial node features to extract local embedded features that fuse physical properties and spatial structure.
[0093] Furthermore, in step S4, the steps of message passing and aggregation are performed on the initial node features through the graph convolutional layer, including: performing nonlinear transformation on the source node features and edge weights to generate a message vector containing local physical gradients; calculating the influence coefficients of neighboring nodes using an attention mechanism, and weighting and aggregating the message vectors on adjacent edges; fusing the aggregated features with its own historical features, and generating a new node embedding vector containing spatial topology information by updating the network.
[0094] Specifically, a graph neural network is constructed to process the graph data containing physical properties and topological structure generated in step S3. To effectively capture the non-uniformity and anisotropy of fluid flow and stress transfer in hydrate reservoirs, a graph attention network (GAT) is adopted as the core architecture, which has been proven to perform well in processing non-Euclidean space data. First, the initial node feature vector generated in step S3 is used as the input to the 0th layer of the graph neural network, denoted as... This vector specifically contains the first Normalized values of pore pressure, temperature, phase saturation, and solid displacement for each grid node.
[0095] The initial node features are processed through graph convolutional layers, involving message passing and aggregation. This begins with non-linear transformations of the features and message generation. In each layer of the GAT... In this process, a shared, learnable linear transformation matrix is utilized. Projecting the feature vectors of all nodes onto a higher-dimensional feature space extracts the implicit physical gradient information. At this point, for each edge in the graph, i.e., connecting nodes... and nodes The edges, combined with the edge attribute weights calculated in step S3. This process generates message vectors containing local physical gradients, which are represented as unweighted transformation features in GAT. Next, an attention mechanism is used to calculate the influence coefficients of neighboring nodes. A single-layer wavefront neural network is constructed as the attention mechanism function. compute nodes Its neighboring nodes Attention coefficient between neighboring nodes, which characterizes the attention coefficient between neighboring nodes. The physical state of the central node The degree of importance. The formula for calculating the attention coefficient is:
[0096] ;
[0097] In the formula, The weight matrix is a linear transformation matrix; and These are the feature vectors of the center node and its neighboring nodes, respectively. This represents a vector concatenation operation; This represents the weight vector for the attention mechanism; This is a non-linear activation function. To ensure the comparability of these coefficients across different nodes, the attention coefficients are normalized using the Softmax function, resulting in the final normalized attention weights. :
[0098] ;
[0099] In the formula, For nodes The set of first-order neighbor nodes. Then, weighted aggregation is performed, using the calculated... As weights, the transformation characteristics of all neighboring nodes are summed in a weighted manner to aggregate the fluid-structure interaction patterns within the local neighborhood.
[0100] Finally, the aggregated features are fused with the node's own historical features to update the network and generate new node embedding vectors. To prevent the gradient vanishing problem caused by an overly deep network and to preserve the original physical information, residual connections are introduced into the network. The aggregated neighborhood features are then integrated with the node's... The node embedding vector is updated by summing its own historical features and then passing the non-linear activation function ELU. The complete node update formula is:
[0101] ;
[0102] In the formula, This is a non-linear activation function. To enhance the model's expressive power, a multi-head attention mechanism is used in practice, i.e., independent computation. Group attention coefficients and aggregate them, ultimately... The features are concatenated or averaged. After being processed by stacking multiple GAT layers, the final output feature vector is the local embedded feature that integrates physical properties and spatial structure.
[0103] This implementation step adaptively learns the dependencies between grid nodes through a graph attention mechanism, which can automatically focus on key nodes in the reservoir where the physical field changes drastically or in the direction of the main flow channel. This allows for the efficient extraction of high-order spatial topological features reflecting fluid-solid transport laws from unstructured grid data, solving the technical problem that traditional convolution methods struggle to handle irregular grids and anisotropic flows.
[0104] S5. Construct a long short-term memory network to receive the time series of locally embedded features and generate a temporal latent vector representing the fluid-solid transport state through a gating mechanism.
[0105] Furthermore, in step S5, the step of generating a temporal latent vector representing the fluid-structure transport state through a gating mechanism includes: setting up a loop structure containing multiple gating units to process the time series of locally embedded features; selectively forgetting historical latent states using gating operators and weighting and fusing the current input features; performing nonlinear mapping on the updated cell state and outputting a temporal latent vector reflecting the evolution of fluid-structure transport.
[0106] Specifically, a standard Long Short-Term Memory (LSTM) network is constructed as the core architecture for time-series processing, used to receive and parse the local embedded feature time series output by the graph neural network in step S4. First, the graph convolutional layers are applied in consecutive simulation time steps... The sequence of locally embedded feature vectors that fuse physical properties and spatial structure output by the network is defined as the input sequence of the network. Each input vector here All of these contain the spatial topological state of the hydrate reservoir and the convergence information of the fluid-structure interaction physical fields at that moment. The LSTM network maintains a cell state. and a hidden state To capture the long-term and short-term dependencies of hydrate decomposition front propulsion and pressure wave diffusion.
[0107] The steps for generating temporal latent vectors representing fluid-structure interaction (FSI) states through a gating mechanism first involve selectively forgetting historical latent states using a forgetting gate operator. The forgetting gate then retrieves the previous time-series latent state representing the historical FSI trend. and the current moment input features that characterize the current physical state of space The sigmoid activation function outputs a forgetting coefficient vector with values between 0 and 1. This coefficient determines which old physical history information in the cell state should be retained or discarded. The formula for calculating the forgetting gate is:
[0108] ;
[0109] In the formula, Use the Sigmoid activation function; Here is the weight matrix for the forget gate; It is the bias vector; This represents the concatenation operation of two vectors.
[0110] The input gate operator is then used to weight and fuse the current input features, updating the cell state of the memory unit. The input gate consists of two parts: one part uses a Sigmoid layer to determine which new physical feature values will be used for updating; the other part uses a Tanh layer to create a new candidate cell state vector. By changing the old cell state Multiplying the result by the forgetting coefficient and adding the new candidate state after input gate weighting, we obtain the cell state updated at the current time step. This process enables the iteration of physical evolution information. The calculation formula is as follows:
[0111] ;
[0112] ;
[0113] ;
[0114] In the formula, The input gate activation vector; and These are the weight matrices for the input gate and the candidate state, respectively; and This is the corresponding bias vector; This represents the Hadamard product.
[0115] Finally, a nonlinear mapping is performed on the updated cell state, outputting a temporal implicit vector reflecting the evolution of fluid-structure transport. The output gate first calculates the output activation vector based on the current input and historical states. The activation vector is then used to process the current cell state after the Tanh function. Filtering is performed to generate the final hidden state at the current moment. This vector This is a high-dimensional temporal implicit vector characterizing the evolution of the entire fluid-solid production process of the hydrate reservoir from the initial moment to the current moment. The calculation formula is:
[0116] ;
[0117] ;
[0118] In the formula, The output gate activation vector; This is the output gate weight matrix; This is the bias vector.
[0119] This implementation step effectively solves the gradient vanishing or exploding problem that traditional recurrent neural networks are prone to when processing long-cycle hydrate mining simulation data through the unique gating mechanism of LSTM. It can accurately capture the time lag effect and cumulative effect in the hydrate decomposition and sand production process, and integrate discrete instantaneous spatial physical features into a continuous fluid-solid transport state representation with historical memory, providing feature inputs that include all spatiotemporal dimensions for subsequent accurate prediction of cumulative gas production and sand production.
[0120] S6. Using a loss function containing prediction error and physical constraint terms, jointly train a graph neural network and a long short-term memory network to establish a nonlinear mapping between the time-series latent vector and the cumulative gas production and cumulative sand production, and generate a prediction model.
[0121] Furthermore, in step S6, the step of jointly training using a loss function containing prediction error and physical constraint terms includes: comparing the model's predicted values with the corresponding cumulative gas production and cumulative sand production in the evolution data to construct a basic loss term characterizing the data fitting accuracy; penalizing the difference values in the prediction sequence that violate physical laws based on the physical monotonicity characteristics of cumulative gas production and cumulative sand production; and constructing a comprehensive objective function by using adaptive weights to balance the basic loss term and physical constraint term.
[0122] Furthermore, in step S6, the step of establishing the nonlinear mapping between the time-series latent vector and the cumulative gas production and cumulative sand production includes: mapping the time-series latent vector to a low-dimensional production space using a fully connected network; calculating the gradient of the comprehensive objective function relative to the trainable parameters of the model based on an automatic differentiation framework; and updating the weights of the graph neural network and the long short-term memory network synchronously according to the gradient using an adaptive optimization strategy.
[0123] Specifically, a decoding module is first constructed to establish a nonlinear mapping relationship. A fully connected neural network is used as the decoder to receive the temporal hidden vector output from step S5. This fully connected network typically includes an input layer, several hidden layers, and an output layer. The hidden layers use the ReLU activation function to increase nonlinear expressiveness, while the output layer uses a linear activation function to adapt to the regression task. The temporal hidden vector generated in step S5 is then processed... The input is fed into this fully connected network, which maps it to a low-dimensional output space and directly outputs the predicted value vector at the current time step. The predicted value vector contains two components: the predicted cumulative gas production. and predicted cumulative sand output The mapping calculation formula for a fully connected network is as follows:
[0124] ;
[0125] In the formula, For predicting output vectors; and These are the weight matrices for the output layer and the hidden layer, respectively. and This is the corresponding bias vector; This is the temporal hidden vector output by the LSTM.
[0126] Next, a joint loss function incorporating prediction error and physical constraint terms is constructed. First, the basic data fitting loss term is calculated and compared with the model's predicted values. Compared with the real evolution data generated by the simulation in step S2 This corresponds to the actual cumulative gas production and cumulative sand production. The mean squared error is used as a metric to calculate the Euclidean distance between the two values, characterizing the model's fitting accuracy at the data level. Data fitting loss term. The calculation is as follows:
[0127] ;
[0128] In the formula, This represents the total number of time steps. Let L2 norm be the square of a vector.
[0129] Subsequently, a physical constraint loss term is constructed. This is based on the physical monotonicity of cumulative output, meaning that both cumulative gas production and cumulative sand production must maintain a non-decreasing trend over time. A penalty is imposed on the difference values in the predicted sequence that violate the physical law. The difference between predicted values at adjacent time steps is calculated; if the difference is negative, it is included in the loss; otherwise, it is ignored. Physical constraint loss term. The calculation is as follows:
[0130] ;
[0131] In the formula, Indicates the first Time of the first The predicted cumulative amount of the product gas or sand; The function here only retains positive values, that is, it generates an error penalty when the accumulated amount at the previous time step is greater than that at the current time step.
[0132] Finally, an integrated objective function is constructed by using adaptive weights to balance the basic loss term and the physical constraint term. :
[0133] ;
[0134] In the formula, These are the weighting coefficients for the physical constraint terms, which can be adaptively adjusted based on the gradient magnitude during training or set as fixed hyperparameters.
[0135] Finally, joint training and parameter updates of the models are performed. The comprehensive objective function is calculated based on the automatic differentiation framework. Relative to the gradient tensors of all trainable parameters in the model, the trainable parameters include the graph convolution weights of the graph neural network, the gated weights of the long short-term memory network, and the mapping weights of the fully connected network. The model parameters are updated using the Adam adaptive optimization strategy, which combines first-order momentum and second-order momentum estimation to adaptively adjust the learning rate for different parameters. The parameter update formula is:
[0136] ;
[0137] In the formula, Represents any weight parameter in the model; The initial learning rate; The corrected first-moment estimate characterizes the gradient mean; The corrected second-moment estimate characterizes the uncentered variance of the gradient; To prevent smoothing terms with a denominator of zero, repeated iterative training is performed until the overall objective function converges.
[0138] This implementation step introduces physical monotonicity constraints, forcing the deep learning model to follow basic physical conservation laws while learning the data distribution. This effectively avoids the physical distortion that may occur in pure data-driven models in areas with scarce samples, thereby significantly improving the robustness and interpretability of the hydrate mining prediction model in long-term extrapolation.
[0139] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.
Claims
1. A method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs, characterized in that, Includes the following steps: S1. Establish a coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction and sediment constitutive relations, and use unstructured grids to discretize the reservoir; S2. Solve the coupled mathematical model using an alternating iterative algorithm to generate evolution data of pore pressure, temperature, phase saturation, solid displacement, and corresponding cumulative gas production and cumulative sand production of the grid nodes of the unstructured grid at continuous time steps. S3. Construct a topology graph based on the unstructured mesh, map the mesh nodes to graph nodes, map the geometric adjacency relationships between the mesh nodes to graph edges, and assign the pore pressure, temperature, saturation of each phase, and solid displacement as initial node features. S4. Construct a graph neural network, and perform message passing and aggregation on the initial node features through graph convolutional layers to extract local embedded features that fuse physical properties and spatial structure. S5. Construct a long short-term memory network to receive the time series of the local embedded features and generate a temporal latent vector representing the fluid-solid transport state through a gating mechanism; S6. Using a loss function containing prediction error and physical constraint terms, jointly train the graph neural network and the long short-term memory network to establish a nonlinear mapping between the time-series hidden vector and the cumulative gas production and cumulative sand production, and generate a prediction model.
2. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, Step S1, the step of establishing a coupled mathematical model covering hydrate decomposition, multiphase flow, heat conduction, and sediment constitutive relations, includes: Define the hydrate decomposition rate affected by saturation and contact area, and establish a mass source term expression for the solid-to-fluid phase transition. By combining relative permeability and capillary pressure function, a set of partial differential equations for mass conservation of multiphase fluids containing gravity and pressure gradient is constructed. An elastoplastic mechanical model is used to describe the stress response of the skeleton, and a fluid-solid output flux model is constructed using a joint criterion of plastic strain and flow velocity.
3. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, In step S1, the step of discretizing the reservoir using an unstructured grid includes: The reservoir geometry is discretized using a spatial partitioning algorithm, and local mesh refinement is performed in regions where the physical field gradient exceeds a preset threshold. Based on geological logging data, an interpolation algorithm is used to generate a spatial attribute field, and continuous geological attributes are projected onto the grid nodes of the unstructured grid. The coupled mathematical model is transformed into a system of algebraic equations using a numerical discretization method, and the flux conservation of the unstructured mesh interface is handled.
4. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, In step S2, the step of solving the coupled mathematical model using an alternating iterative algorithm includes: Set the time step parameters and construct a nonlinear solver that uses the solution from the previous time step as the initial value; A segmented iterative strategy is adopted to solve the fluid and temperature fields first, and then update the solid mechanical and chemical field parameters. Calculate the residual norm of each physical field between iterations, and output the state vector and advance the time step when the residual satisfies the preset tolerance.
5. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, Step S3, the step of constructing a topology graph based on the unstructured mesh, includes: Using the grid nodes as vertices, the normalized multiphysics state data is combined as the initial node attribute vector; Based on the topological connection relationships between the grid nodes, connected edges are established to form an adjacency matrix describing spatial connectivity; The weight coefficients are calculated based on the spatial distance and transmission characteristics between adjacent graph nodes to characterize the interaction strength between grids.
6. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, In step S4, the step of performing message passing and aggregation on the initial node features through graph convolutional layers includes: Perform a nonlinear transformation on the source node features and edge weights to generate a message vector containing local physical gradients; The influence coefficients of neighboring nodes are calculated using an attention mechanism, and the message vectors on the adjacent edges are weighted and aggregated. By fusing aggregated features with its own historical features, the network is updated to generate new node embedding vectors that contain spatial topology information.
7. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, In step S5, the step of generating the temporal implicit vector representing the fluid-structure transport state through a gating mechanism includes: A loop structure containing multiple gating units is configured to process the time series of the local embedded features; The gating operator is used to selectively forget the historical hidden states and to perform weighted fusion of the current input features; A nonlinear mapping is performed on the updated cell state to output the temporal implicit vector that reflects the evolution of fluid-solid transport.
8. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, Step S6, the step of jointly training using a loss function containing prediction error and physical constraint terms, includes: By comparing the model's predicted values with the corresponding cumulative gas production and cumulative sand production in the evolution data, a basic loss term characterizing the data fitting accuracy is constructed. Based on the physical monotonicity characteristics of the cumulative gas production and cumulative sand production, a penalty is imposed on the difference values in the prediction sequence that violate physical laws. An integrated objective function is constructed by using adaptive weights to balance the basic loss term and the physical constraint term.
9. The method for constructing a multi-field coupled prediction model for fluid-solid production in hydrate reservoirs according to claim 1, characterized in that, Step S6, the step of establishing the nonlinear mapping between the time-series implicit vector and the cumulative gas production and cumulative sand production, includes: The temporal latent vector is mapped to a low-dimensional output space using a fully connected network; The gradient of the synthetic objective function with respect to the trainable parameters of the model is calculated based on the automatic differentiation framework; The weights of the graph neural network and the long short-term memory network are updated synchronously according to the gradient using an adaptive optimization strategy.