Micro-tomography method based on physical information graph neural network and neural level set

Abstract

The invention belongs to the technical field of crossing of oil and gas reservoir geophysics and artificial intelligence, and particularly relates to a micro-tomography method based on a physical information graph neural network and a neural level set. The method comprises the following steps: modeling a reservoir into a graph structure; deploying a physical information neural network at each node to solve a pressure and saturation field; constructing a neural level set for the potential micro-fault edge to dynamically represent the geometry and closure of the potential micro-fault edge; calculating a fault transmission coefficient through an edge network, carrying out message transmission, and updating a node state; constructing a graph neural network; and finally, through optimizing a comprehensive loss function, synchronously inverting pressure-saturation field and micro fault parameters. According to the method, the physical information graph neural network and the neural level set are combined for the first time, dynamic production data are utilized, the resolution limitation of a traditional seismic method is broken through, real-time, high-precision and physically-credible inversion imaging of meter-scale micro-faults is achieved, and core technical support is provided for intelligent regulation and control of oil and gas fields.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of oil and gas reservoir geophysics and artificial intelligence, specifically a micro-fault imaging method based on physical information graph neural networks and neural level sets. Background Technology

[0002] In oil and gas field development, especially in mature oilfields, accurate characterization of subsurface fault systems is crucial for simulating fluid flow, optimizing injection and production schemes, and improving ultimate recovery. Microfaults, with scales of 1-10 meters, are far below the wavelength limit of conventional 3D seismic exploration. Although difficult to identify by traditional seismic exploration, they largely control injected water and water intrusion breakthrough, as well as inefficient injection-production cycles, severely impacting development outcomes. Existing technical solutions suffer from the following fundamental defects:

[0003] 1. Seismic interpretation technology: Based on three-dimensional seismic data volume, faults are identified through attributes such as coherence and curvature. However, its resolution is limited by the seismic wavelength (usually >30 meters), and it cannot effectively identify micro-faults at the meter level.

[0004] 2. Well test interpretation techniques: Pressure transient analysis, tracer testing, and pulse well tests require well shut-in operations, which seriously affect normal production;

[0005] 3. Numerical inversion techniques: such as grid-based adjoint method inversion, which adjusts model parameters by fitting production data. This method is extremely computationally expensive, with a single inversion taking several days or even weeks, making it difficult to meet real-time decision-making needs, and it is highly dependent on the initial model.

[0006] In recent years, Physical Information Neural Networks (PINNs) have achieved meshless solutions to forward problems by embedding the residuals of partial differential equations into a loss function; Graph Neural Networks (GNNs) have demonstrated powerful modeling capabilities in non-Euclidean domains. However, current techniques only combine PINNs and GNNs for fast forward approximation, failing to address the inverse problem of retrieving hidden geological structures from real-time production data. Therefore, there is an urgent need to develop a novel method that can integrate high-frequency IoT data, strictly satisfy physical conservation laws, and achieve real-time imaging of meter-level micro-faults. Summary of the Invention

[0007] To overcome the shortcomings of existing technologies, the present invention aims to provide a fault inversion method based on neural level sets and Graph-PINN (Physical Information Graph Neural Network). This method aims to achieve automated, real-time, and high-precision imaging of concealed meter-level faults in reservoirs, and simultaneously invert their closure characteristics. It provides a high-precision, low-latency, and loop-closed-loop micro-fault imaging method based on Graph-PINN, offering direct evidence for intelligent control and optimization of oilfields.

[0008] The technical solution adopted by this invention to achieve the above objectives is: a micro-tomographic imaging method based on physical information graph neural networks and neural level sets, comprising the following steps:

[0009] Step S1: Construct a reservoir model and obtain dynamic and static monitoring and production datasets including well flow rate and pressure; model the reservoir as a graph structure, where fault blocks are nodes, known faults and potential micro-faults generated through sampling are edges, and encode reservoir physical properties for each node;

[0010] Step S2: Construct a physical information neural network inside each node, input spatiotemporal coordinates to solve and output the pressure field and saturation field of the internal space of the node;

[0011] Step S3: Encode the geological and physical attributes of each node to obtain a feature vector representing the state of that node;

[0012] Step S4: Construct a neural level set function for each edge representing a potential micro-twist, and perform activation judgment by sampling in its vicinity to dynamically learn and characterize the geometry, closure and existence of the micro-twist corresponding to the edge;

[0013] Step S5: Based on the node feature vectors and neural level sets, design edge network message passing rules, calculate the tomographic transmission coefficient and interface throughput, update node states, and ensure mass conservation at the interface.

[0014] Step S6: Construct a graph neural network, update the physical state of all nodes in the graph based on the interface throughput message passing; and decode the updated node feature vectors to output the final node pressure and saturation prediction values.

[0015] Step S7: Model training and optimization: Construct a comprehensive loss function and use an optimizer to iteratively update all learnable parameters, and simultaneously invert to obtain the reservoir pressure-saturation field and microfault parameters.

[0016] Step S1 specifically includes:

[0017] The edges include: fixed edges and candidate edges;

[0018] Step S1-1: Treat the interpreted faults between adjacent blocks as fixed edges, and set the fault transmission coefficient β. e Based on the original geological model, potential microfaults were generated in the reservoir as candidate edges using the Poisson disk sampling method.

[0019] Step S1-2: For candidate edges, their initial fault transmission coefficients Set to 0;

[0020] The nodes are fault blocks that divide the reservoir; node characteristics include: pressure P, saturation S, and permeability. Porosity Rock compressibility PVT data and relative permeability data;

[0021] Convert the corner coordinates to the center point coordinates of the unit block r = (x, y, z). Number each block as ID b(r)∈{1, ..., N}.

[0022] Step S2 specifically includes:

[0023] Step S2-1: The physical information neural network is a fully connected network. Its input is the spatiotemporal coordinates ξ=(x, y, z,t), and its output is the pressure P and the saturation S, expressed as:

[0024]

[0025] in, The neural network representing node v, For network parameters;

[0026] Step S2-2: The training constraint of the physical information neural network is to minimize the residual of the network output that satisfies the reservoir seepage control equation. For oil-water two-phase systems, the PDE residual loss within the nodes is achieved by minimizing the residual, i.e.:

[0027] ;

[0028] ;

[0029] in, , Residuals for oil and water phases, respectively. Here, k is the absolute permeability, and k is the vector differential operator. , These are the relative permeabilities of oil and water, respectively. , These are the oil and water phase pressures, respectively. , These represent the oil and water phase saturation, respectively. , These are the volume coefficients of the oil and water phases, respectively. , The densities of the oil and water phases are respectively. , The viscosity of the oil phase and the viscosity of the water phase are respectively. , These are the source and sink of oil and water, respectively. Porosity This is the acceleration due to gravity.

[0030] Step S4 specifically includes:

[0031] Step S4-1: Construct a scalar level set function for each candidate edge e, implemented by a fully connected neural network. ,in These are learnable parameters;

[0032] Step S4-2: During training, sampling points are randomly set in the space near each candidate edge, and the level set function value of each point is calculated; and a non-differentiable hard decision is performed during forward propagation, that is: activation is determined based on the maximum and minimum values ​​of the function, and if the maximum value is... and minimum value Then activate the candidate edge so that it corresponds to a real micro-fault, and its corresponding fault transmission coefficient β e It will participate in subsequent optimizations; otherwise, the edge is considered non-existent, and its fault transmission coefficient β... e It was forcibly set to 1;

[0033] Step S4-3: Optimize the level set parameters during backpropagation In this case, a pass-through estimator is used to handle the non-differentiability of activation decisions, and an identity mapping is used to replace the non-existent gradient. , that is to say This allows the gradient to propagate back to the level set network parameters. .

[0034] Step S4 further includes introducing a regularization term into the total loss function, specifically:

[0035] Introducing a geometric smoothing regularization term: It is used to punish drastic changes in the level set function and encourage smooth fault interfaces;

[0036] Introducing a sparse regularization term: Encourage the sparsity of fault sets;

[0037] in, , All of these are hyperparameter control factors.

[0038] Step S5 includes the following steps:

[0039] Step S5-1: For the connection node and edge Its fault transmission coefficient Through a fully connected network Calculation by combining node latent vectors and activation states:

[0040] ;

[0041] in, Given the preset transmission coefficients of the fixed side, , The latent vectors of nodes are obtained by the node feature encoder. For edge networks, the Softplus activation function is defined as Softplus(x) = log(1 + exp(x));

[0042] Step S5-2: For the edge From node To the node News Calculated based on fault transmission coefficient, pressure gradient, and interface normal vector:

[0043] ;

[0044] Step S5-3: Updated state of node i Based on its current state With all neighboring nodes The sum of the messages is obtained by aggregation, and satisfies ,Right now:

[0045] ;

[0046] in, Let be the unit normal vector of the interface, ∇P be the pressure gradient of node j, and N(i) be the set of neighbors of node i. Let be the feature vector of node i in the l-th layer.

[0047] In step S6, constructing the graph neural network specifically involves:

[0048] Step S6-1: Using the node feature latent vector and edge features obtained from the node feature encoder as input, calculate the interface flux based on the fault transmission coefficient, update the physical state of each node according to the message passing rules, and construct a graph neural network containing at least 2 hidden layers.

[0049] Step S6-2: Use the updated node state of the graph neural network as the input of the decoder network, and the output of the decoder network is the predicted pressure P and saturation S of each node.

[0050] In step S7, constructing the total loss function specifically involves:

[0051] Step S7-1: In PINN, the residual loss within the node To enable neural network prediction Minimize the residuals that satisfy the physical governing equations, i.e.:

[0052]

[0053] in, , These are the residuals for the oil and water phases, respectively.

[0054] Step S7-2: The data loss in the PINN network inside the node is as follows:

[0055]

[0056] in, Pressure values ​​predicted by the PINN (Physical Information Neural Network) model; Real, observed stress values ​​were used as training labels; Saturation values ​​predicted by the PINN (Physical Information Neural Network) model; The actual, observed saturation values ​​are used as training labels;

[0057] The node data loss function output by the graph neural network for:

[0058]

[0059] in, To predict node stress using a graph neural network; To predict node saturation using a graph neural network; The characteristic value of the fault block node pressure; The saturation characteristic value of the block node;

[0060] The total data loss of the graph neural network nodes is:

[0061]

[0062] Step S7-3: To ensure that pressure and flow are continuous at the interface, the interface connection condition loss must meet the following conditions:

[0063]

[0064] in, and The pressure on both sides of the interface; and The fluid velocity of adjacent cells; The interface normal vector;

[0065] Step S7-4: Fit the pressure and flow rate observation data at the well point, and obtain the fitting loss of the observation data as follows:

[0066]

[0067] in, These are predicted pressure and observed pressure, respectively. These are predicted flow and observed flow, respectively.

[0068] Step S7-5: Loss through the internal physical equations of the node Graph Neural Network Node Loss Interface connection condition loss Observational data loss Based on the calculated regularization loss, construct the total loss function. for:

[0069]

[0070] in, and These are the geometric smoothing regularization term and the sparse regularization term, respectively.

[0071] In step S7, synchronous inversion is performed by optimizing the total loss function, specifically as follows:

[0072] Parameters of the neural network for the internal physical information of nodes Edge network Parameters, edge transfer coefficients and level set function parameters All parameters are updated using the Adam optimizer.

[0073] In the later stages of optimization, the L-BFGS algorithm is used for convergence to obtain more accurate inversion results.

[0074] An imaging system based on a physical information graph neural network and neural level set for microtomography includes:

[0075] Data acquisition and preprocessing module: used to build reservoir models and acquire and process dynamic and static monitoring and production datasets, including well flow and pressure.

[0076] Graph data structure construction module: used to convert the reservoir model into a graph data structure, map fault blocks in the reservoir as nodes, map known fault interfaces and potential micro-fault interfaces as edges connecting nodes, and encode reservoir physical properties for nodes;

[0077] Physical information neural network module: includes physical information neural network units deployed at each node, used to input spatiotemporal coordinates and solve for the pressure field and saturation field inside the node;

[0078] Node feature encoding module: used to encode the geological and physical attributes of each node and generate node feature vectors;

[0079] The neural level set construction module is used to construct and train neural level set functions for each potential micro-fault edge to dynamically characterize the geometry, closure, and existence of the micro-fault.

[0080] Edge network and message passing module: It is used to calculate the tomographic transmission coefficient through the edge network based on the node feature vector and the neural level set, and to perform message passing and update the node state according to the transmission coefficient, node pressure gradient and interface geometry information.

[0081] Graph Neural Network and Decoding Module: Used to build graph neural networks, update the physical state of all nodes based on message passing, decode the updated node features, and output the predicted values ​​of node pressure and saturation.

[0082] Inversion optimization module: used to construct a comprehensive loss function and use an optimizer to iteratively update all learnable parameters to achieve synchronous inversion of pressure-saturation field and microfault parameters.

[0083] The present invention has the following beneficial effects and advantages:

[0084] 1. Ultra-high resolution and recognition capability: By integrating high-frequency monitoring data with physical information graph neural networks and combining neural level sets for fine characterization of fault interfaces, this invention can effectively detect and characterize 1-10 meter level micro-faults that cannot be distinguished by traditional 3D seismic data, significantly improving the recognition accuracy of hidden geological structures.

[0085] 2. Real-time and efficient, meeting closed-loop control requirements: After model training converges, the inference time on a single GPU can be reduced to less than 1 second, supporting data assimilation and model update closed loops on hourly or even shorter time scales. Compared with traditional grid-based numerical inversion methods (such as the finite volume method, which takes about 35 seconds per iteration), the computational efficiency is improved by tens of times, truly meeting the timeliness requirements of real-time optimization and intelligent control of oilfield production.

[0086] 3. Strong physical consistency and interpretability: This invention directly embeds the reservoir seepage control equation as a physical constraint into the neural network training loss function, ensuring that the inversion process strictly follows the fluid flow law and the inversion results have clear physical meaning, overcoming the "black box" drawbacks of pure data-driven models that lack physical constraints and whose results are difficult to interpret.

[0087] 4. Good scalability and parallel efficiency: The graph-based modeling method of this invention is naturally adapted to complex fault network topologies and supports large-scale models (such as 128 fault blocks and 1000 candidate edges). In a multi-GPU parallel computing environment (such as 16 GPUs), it can still maintain a parallel efficiency of up to 85%, which has good scalability.

[0088] 5. Dynamic adaptive representation capability: This invention introduces a learnable neural level set, which does not require prior assumptions such as the location, shape and number of faults. It can adaptively learn and evolve the geometry and closure (transmission coefficient) of potential micro-faults based on dynamic data, breaking through the limitations of traditional parametric inversion methods that rely on fixed geological models.

[0089] 6. Significant cost and operational advantages: This invention only requires the use of high-frequency monitoring data such as wellbore flow and pressure during normal production processes, without the need for well shut-in operations (such as well testing) or the deployment of expensive high-density seismic observations. The monitoring cost is reduced by more than an order of magnitude compared to four-dimensional seismic acquisition, while the spatial resolution is improved by an order of magnitude, achieving "cost reduction and efficiency improvement".

[0090] 7. Framework versatility and broad application prospects: The method framework of this invention is open source and modular. By replacing the corresponding control equations (such as heat conduction equations, multiphase flow equations, etc.), it can be seamlessly extended to other underground resource and energy fields such as carbon dioxide geological storage, geothermal energy extraction, and underground hydrogen storage operation monitoring, providing general technical support for high-precision characterization and safe control of various underground spaces. Attached Figure Description

[0091] Figure 1 This is a flowchart of the microtomographic imaging method of the present invention. Detailed Implementation

[0092] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments. Those skilled in the art should understand that the following description is intended to clearly demonstrate the implementation process of the present invention, and not to limit the scope of protection of the present invention. Reasonable modifications or substitutions to the details of this embodiment without departing from the core principles of the present invention should be considered to fall within the scope of protection of the present invention.

[0093] This embodiment uses the geological structure of layers 20 to 30 of the industry-standard SPE10 model as a background to construct a conceptual reservoir model that includes known faults and potential micro-faults. The model deploys four injection wells and five production wells, employing a constant-fluid production regime to simulate 500 days of production dynamics, and generates corresponding bottomhole pressure and production monitoring data as the dynamic observation data source required for the inversion method of this invention. The following will step-by-step description of the specific steps for achieving micro-fault imaging based on a Graph-PINN neural network and neural level sets.

[0094] like Figure 1 The diagram shown is a flowchart of the microtomographic imaging method of the present invention. The specific steps of this embodiment are as follows:

[0095] 1. Reservoir Model Construction

[0096] In this embodiment, layers 20 to 30 (60×220×11 units) of the SPE10 model were selected as the basic geological layers. Three faults of 2-15 meters perpendicular to the strike of the rock strata were set. The reservoir fluid and rock physical parameters followed the original model. Four injection wells and five production wells were operated using a constant fluid production scheme, simulating production for 500 days.

[0097] 2. Construction of Graph Data Structure

[0098] Edge: Represents the interface. (1) Fixed edge: The interpreted fault between adjacent blocks is used as a fixed edge, and its fault transmission coefficient β e Determined by the original geological model. (2) Candidate edges: Candidate edges are generated by the Poisson disk sampling method with a minimum spacing of 5 meters within each pair of blocks. For candidate edge e, the initial transfer coefficient is... and instantiate a learnable level set. .

[0099] Nodes: Represent fault blocks. Fixed edges and candidate edges divide the entire reservoir into several small fault blocks, each of which, b, serves as a graph node. Node characteristics include pressure P, saturation S, and permeability tensor. Porosity Rock compressibility PVT data and relative permeability data.

[0100] Mesh transformation: Convert corner coordinates to cell center coordinates r = (x, y, z). Number each cell block ID b(r)∈{1, ..., N}.

[0101] 3. Node Physical Information Neural Network

[0102] Each node The model incorporates a pressure P and a saturation S, and is an eight-layer fully connected network. The input is the spatiotemporal coordinates ξ=(x, y, z, t), and the output is the pressure P and saturation S.

[0103]

[0104] The governing equations adopt the seepage equations of the black oil model. Taking the oil-water two-phase flow as an example, the specific formulas are as follows:

[0105]

[0106]

[0107] in, Here, k is the absolute permeability, and k is the vector differential operator. , These are the relative permeabilities of oil and water, respectively. , These are the oil and water phase pressures, respectively. , These represent the oil and water phase saturation, respectively. , These are the volume coefficients of the oil and water phases, respectively. , The densities of the oil and water phases are respectively. , The viscosity of the oil phase and the viscosity of the water phase are respectively. , These are the source and sink of oil and water, respectively. Porosity This is the acceleration due to gravity.

[0108] In PINN, the residual loss within a node is what the neural network predicts. Once the residuals satisfying the physical governing equations are minimized, the residual calculation formula for the oil-water phase is as follows:

[0109]

[0110]

[0111] in, , These are the residuals for the oil and water phases, respectively.

[0112] 4. Node Feature Encoder

[0113] The original features of each node are transformed by a fully connected network MLP, which maps them from the initial high-dimensional original feature space to a latent vector h. This learned vector can more effectively capture the high-level semantic information of the node, thus providing high-quality input features for subsequent tasks.

[0114] 5. Construction of neural level sets

[0115] 5.1 Network Construction

[0116] Each candidate edge corresponds to a scalar function. The function is constructed using a fully connected MLP neural network. During training, 5000 sampling points are randomly placed near each edge. If the maximum value is greater than 0 and the minimum value is less than 0, the zero surface will cut the interface and activate the edge, representing a real micro-fault. The transmission coefficient of this edge... Participate in message passing and optimization; otherwise It will be forcibly set to 1, which is equivalent to the edge not existing physically.

[0117] The conditions for activating an edge are as follows:

[0118]

[0119] 5.2 STE (Straight-through Estimator)

[0120] The non-differentiability problem of the activation mechanism is solved by using a direct-through estimator (STE). During the propagation of the current term, a hard check for non-differentiability is performed.

[0121]

[0122] The hard-determined output A (0 or 1), where 0 represents no fault and 1 represents a fault, is used to determine the transmission coefficient of the edge. It then participates in all subsequent physical calculations and the evaluation of the loss function. At this point, the system's behavior strictly conforms to the physical definition.

[0123] During backpropagation, the loss function φ is calculated with respect to the level set parameters ω. e When calculating the gradient, the chain rule gets stuck at activation decision A because... Does not exist:

[0124]

[0125] Here, the STE identity approximation is used to replace the real, non-existent ones. :

[0126]

[0127] Assume activation state A is paired with The changes are extremely sensitive, and the rate of change is always 1, so the gradient passes through directly and without decay.

[0128] 5.3 Regularization Terms

[0129] To obtain physically reasonable and stable solutions, and to ensure geometric reasonableness and sparsity, two regularization terms are introduced: a geometric smoothing regularization term and a sparsity regularization term.

[0130] (1) Geometric Smoothing Regularization: Penalizes drastic changes in the level set function and encourages smooth geometric interfaces. This prevents neural networks from fitting non-physical, oscillating fault shapes. Similar to the "curvature flow" smoothing effect in traditional level set methods, the formula is as follows:

[0131]

[0132] (2) Sparsity Regularization: This encourages fault sets to be sparse, meaning only a few candidate edges are ultimately activated. Geological priors indicate that faults are not ubiquitous but rather locally developed structures. This term will increase the fault transport coefficient of many edges. The value approaches zero, thus working in conjunction with the activation mechanism to achieve "pruning while pruning," as shown in the following formula:

[0133]

[0134] 6. Side-network message passing

[0135] For the edge Its interface area is A e Unit normal vector n e Build a fully connected network N β The network will splice the latent feature vector h of the encoder node. i and h j Mapped to nonnegative fault transmission coefficients:

[0136]

[0137] in, The activation function is Softplus(x) = log(1 + exp(x)); for the interpreted fault (fixed edge), its It is a fixed value pre-set based on prior information, which is not updated during training to ensure that known geological structures are respected; for potential micro-faults (candidate edges). Initialized to 0, its final transmission coefficient Relying on the first part of the formula So Learned from data.

[0138] Message passing Propagation along the edge, node updates are as follows:

[0139]

[0140] in, Let be the set of all neighboring nodes; let be the feature vector (hidden state) of node i at layer l; let be the set of all neighboring nodes of node i; let be the message passed from node i to node j; let be the updated features of node i at layer l+1; and let ∇P be the pressure gradient of node j. The physical meaning represents the interface flux.

[0141] For each edge e = (i, j), we have This ensures the conservation of mass at the interface.

[0142] 7. Construction of Graph Neural Networks

[0143] Using the node feature latent vectors and edge features obtained from the node feature encoder as input, the interface flux is calculated based on the fault transmission coefficient, i.e., following the... The message passing rules are used to update the physical state of each node and construct a graph neural network with two hidden layers.

[0144] The updated node state of the graph neural network is used as input, and the node pressure P and node saturation S are used as output to construct a decoder network (MLP).

[0145] The loss function is calculated as follows:

[0146]

[0147] in, Graph neural networks predict node stress; To predict node saturation using a graph neural network; The characteristic value of the fault block node pressure; The saturation characteristic value of the block node;

[0148] 8. Model Training and Optimization

[0149] 8.1 Design the loss function

[0150] The loss function is determined by the loss of the internal physical equations of the node. Graph Neural Network Node Loss Interface connection condition loss Observational data loss Regularization loss and The composition and calculation formula are as follows:

[0151]

[0152] (1) PDE residual loss

[0153] In PINN, the residual loss within a node It's about letting the neural network predict. Minimize the residuals that satisfy the physical governing equations:

[0154]

[0155]

[0156]

[0157] (2) Data loss

[0158] The data loss within the PINN network inside the node is as follows:

[0159]

[0160] in, Pressure values ​​predicted by the PINN (Physical Information Neural Network) model; Real, observed stress values ​​were used as training labels; Saturation values ​​predicted by the PINN (Physical Information Neural Network) model; The actual, observed saturation values ​​are used as training labels;

[0161] The loss function for a graph neural network is calculated as follows:

[0162]

[0163] in, To predict node stress using a graph neural network; To predict node saturation using a graph neural network; The characteristic value of the fault block node pressure; The saturation characteristic value of the block node;

[0164] The total data loss of the graph neural network nodes is:

[0165]

[0166] (3) Interface connection condition loss

[0167] To ensure that pressure and flow are continuous on the interface, the following conditions must be met:

[0168]

[0169] in, and The pressure on both sides of the interface; and The fluid velocity of adjacent cells; The interface normal vector;

[0170] (4) Loss of fitting the observed data

[0171] Fitting pressure and flow observation data at the well point:

[0172]

[0173] in, These are predicted pressure and observed pressure, respectively. These are predicted flow and observed flow, respectively.

[0174] (5) Regularization loss

[0175] To ensure geometric rationality and sparsity, geometric smoothing regularization and sparsity regularization terms are introduced:

[0176]

[0177]

[0178] λ1, λ2 — Regulatory factors, hyperparameters.

[0179] 8.2 Optimization Algorithm

[0180] Parameters of the neural network for the internal physical information of nodes Edge network Parameters, edge transfer coefficients and level set function parameters All algorithms use the Adam optimizer for parameter updates and the L-BFGS algorithm to achieve fine convergence.

[0181] Through the above steps, the entire process from data preparation, graph structure construction, physical information constraints, dynamic representation of neural level sets, graph neural network message passing and updating, to final loss function optimization and parameter inversion is completed. Experimental results show that this method successfully identified and finely characterized the preset meter-level micro-faults on the test model. The geometric morphology and closure parameter inversion results are in good agreement with the preset model, and the pressure and saturation field inversion accuracy is significantly higher than that of traditional methods. This invention, by organically integrating the spatial topology modeling capability of graph neural networks, the physical law constraint capability of physical information neural networks, and the dynamic interface representation capability of neural level sets, achieves for the first time high-resolution, real-time, and physically reliable inversion imaging of hidden micro-faults using only production dynamic data. The computational efficiency meets the needs of real-time optimization in the field, providing a practical and feasible new generation of technical means for intelligent development and control of oil and gas fields.

[0182] In summary, this invention overcomes the limitations of traditional methods that rely on static geological data to interpret faults, and achieves meter-level high-precision fault inversion based on dynamic data. It significantly improves the accuracy and efficiency of fault identification, location, and closure estimation, providing core technical support for intelligent regulation and development optimization of oil and gas fields. It has broad application prospects in oil and gas reservoir management, geothermal exploitation, and geological carbon sequestration.

[0183] Those skilled in the art will understand that the above description is merely a preferred embodiment of the present invention, and the features described in the various embodiments and / or claims of this disclosure can be combined or combined in various ways, even if such combinations or combinations are not explicitly described in this disclosure. This 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.

[0184] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention. Clearly, those skilled in the art can make various alterations and modifications to the invention without departing from its spirit and scope. Thus, if these modifications and modifications of the invention fall within the scope of the claims and their equivalents, the invention is also intended to include these modifications and modifications.

Claims

1. A microtomographic imaging method based on physical information graph neural networks and neural level sets, characterized in that, Includes the following steps: Step S1: Construct a reservoir model and obtain dynamic and static monitoring and production datasets, including well flow rate and pressure; The reservoir is modeled as a graph structure, where fault blocks are nodes, known faults and potential micro-faults generated by sampling are edges, and reservoir physical properties are encoded for each node. Step S2: Construct a physical information neural network inside each node, input spatiotemporal coordinates to solve and output the pressure field and saturation field of the internal space of the node; Step S3: Encode the geological and physical attributes of each node to obtain a feature vector representing the state of that node; Step S4: Construct a neural level set function for each edge representing a potential micro-twist, and perform activation judgment by sampling in its vicinity to dynamically learn and characterize the geometry, closure and existence of the micro-twist corresponding to the edge; Step S5: Based on the node feature vectors and neural level sets, design edge network message passing rules, calculate the tomographic transmission coefficient and interface throughput, update node states, and ensure mass conservation at the interface. Step S6: Construct a graph neural network, update the physical state of all nodes in the graph based on the interface throughput message passing; and decode the updated node feature vectors to output the final node pressure and saturation prediction values. Step S7: Model training and optimization: Construct a comprehensive loss function and use an optimizer to iteratively update all learnable parameters, and simultaneously invert to obtain the reservoir pressure-saturation field and microfault parameters.

2. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, Step S1 specifically includes: The edges include: fixed edges and candidate edges; Step S1-1: Treat the interpreted faults between adjacent blocks as fixed edges, and set the fault transmission coefficient β. e Based on the original geological model, potential microfaults were generated in the reservoir as candidate edges using the Poisson disk sampling method. Step S1-2: For candidate edges, their initial fault transmission coefficients Set to 0; The nodes are fault blocks that divide the reservoir; node characteristics include: pressure P, saturation S, and permeability. Porosity Rock compressibility PVT data and relative permeability data; Convert the corner coordinates to the center point coordinates of the unit block r = (x, y, z). Number each block as ID b(r)∈{1, ..., N}.

3. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, Step S2 specifically includes: Step S2-1: The physical information neural network is a fully connected network. Its input is the spatiotemporal coordinates ξ=(x, y, z, t), and its output is the pressure P and the saturation S, expressed as: in, The neural network representing node v, For network parameters; Step S2-2: The training constraint of the physical information neural network is to minimize the residual of the network output that satisfies the reservoir seepage control equation. For oil-water two-phase systems, the PDE residual loss within the nodes is achieved by minimizing the residual, i.e.: ； ； in, , Residuals for oil and water phases, respectively. Here, k is the absolute permeability, and k is the vector differential operator. , These are the relative permeabilities of oil and water, respectively. , These are the oil and water phase pressures, respectively. , These represent the oil and water phase saturation, respectively. , These are the volume coefficients of the oil and water phases, respectively. , The densities of the oil and water phases are respectively. , The viscosity of the oil phase and the viscosity of the water phase are respectively. , These are the source and sink of oil and water, respectively. Porosity This is the acceleration due to gravity.

4. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, Step S4 specifically includes: Step S4-1: Construct a scalar level set function for each candidate edge e, implemented by a fully connected neural network. ,in These are learnable parameters; Step S4-2: During training, sampling points are randomly set in the space near each candidate edge, and the level set function value of each point is calculated; and a non-differentiable hard decision is performed during forward propagation, that is: activation is determined based on the maximum and minimum values ​​of the function, and if the maximum value is... and minimum value Then activate the candidate edge so that it corresponds to a real micro-fault, and its corresponding fault transmission coefficient β e It will participate in subsequent optimizations; otherwise, the edge is considered non-existent, and its fault transmission coefficient β... e It was forcibly set to 1; Step S4-3: Optimize the level set parameters during backpropagation In this case, a pass-through estimator is used to handle the non-differentiability of activation decisions, and an identity mapping is used to replace the non-existent gradient. , that is to say This allows the gradient to propagate back to the level set network parameters. .

5. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 4, characterized in that, Step S4 further includes introducing a regularization term into the total loss function, specifically: Introducing a geometric smoothing regularization term: It is used to punish drastic changes in the level set function and encourage smooth fault interfaces; Introducing a sparse regularization term: Encourage the sparsity of fault sets; in, , All of these are hyperparameter control factors.

6. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, Step S5 includes the following steps: Step S5-1: For the connection node and edge Its fault transmission coefficient Through a fully connected network Calculation by combining node latent vectors and activation states: ； in, Given the preset transmission coefficients of the fixed side, , The latent vectors of nodes are obtained by the node feature encoder. For edge networks, the Softplus activation function is defined as Softplus(x) = log(1 + exp(x)); Step S5-2: For the edge From node To the node News Calculated based on fault transmission coefficient, pressure gradient, and interface normal vector: ； Step S5-3: Updated state of node i Based on its current state With all neighboring nodes The sum of the messages is obtained by aggregation, and satisfies ,Right now: ； in, Let be the unit normal vector of the interface, ∇P be the pressure gradient of node j, and N(i) be the set of neighbors of node i. Let be the feature vector of node i in the l-th layer.

7. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, In step S6, constructing the graph neural network specifically involves: Step S6-1: Using the node feature latent vector and edge features obtained from the node feature encoder as input, calculate the interface flux based on the fault transmission coefficient, update the physical state of each node according to the message passing rules, and construct a graph neural network containing at least 2 hidden layers. Step S6-2: Use the updated node state of the graph neural network as the input of the decoder network, and the output of the decoder network is the predicted pressure P and saturation S of each node.

8. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, In step S7, constructing the total loss function specifically involves: Step S7-1: In PINN, the residual loss within the node To enable neural network prediction Minimize the residuals that satisfy the physical governing equations, i.e.: in, , These are the residuals for the oil and water phases, respectively. Step S7-2: The data loss in the PINN network inside the node is as follows: in, Pressure values ​​predicted by the PINN (Physical Information Neural Network) model; Real, observed stress values ​​were used as training labels; Saturation values ​​predicted by the PINN (Physical Information Neural Network) model; The actual, observed saturation values ​​are used as training labels; The node data loss function output by the graph neural network for: in, To predict node stress using a graph neural network; To predict node saturation using a graph neural network; The characteristic value of the pressure at the fault node; The saturation characteristic value of the block node; The total data loss of the graph neural network nodes is: Step S7-3: To ensure that pressure and flow are continuous at the interface, the interface connection condition loss must meet the following conditions: in, and The pressure on both sides of the interface; and The fluid velocity of adjacent cells; The interface normal vector; Step S7-4: Fit the pressure and flow rate observation data at the well point, and obtain the fitting loss of the observation data as follows: in, These are predicted pressure and observed pressure, respectively. These are predicted flow and observed flow, respectively. Step S7-5: Loss through the internal physical equations of the node Graph Neural Network Node Loss Interface connection condition loss Observational data loss Based on the calculated regularization loss, construct the total loss function. for: in, and These are the geometric smoothing regularization term and the sparse regularization term, respectively.

9. The microtomographic imaging method based on physical information graph neural networks and neural level sets according to claim 1, characterized in that, In step S7, synchronous inversion is performed by optimizing the total loss function, specifically as follows: Parameters of the neural network for the internal physical information of nodes edge network Parameters, edge transfer coefficients and level set function parameters All parameters are updated using the Adam optimizer. In the later stages of optimization, the L-BFGS algorithm is used for convergence to obtain more accurate inversion results.

10. The imaging system of the microtomographic imaging method based on physical information graph neural networks and neural level sets according to any one of claims 1-9, characterized in that, include: Data acquisition and preprocessing module: used to build reservoir models and acquire and process dynamic and static monitoring and production datasets, including well flow and pressure. Graph data structure construction module: used to convert the reservoir model into a graph data structure, map fault blocks in the reservoir as nodes, map known fault interfaces and potential micro-fault interfaces as edges connecting nodes, and encode reservoir physical properties for nodes; Physical information neural network module: includes physical information neural network units deployed at each node, used to input spatiotemporal coordinates and solve for the pressure field and saturation field inside the node; Node feature encoding module: used to encode the geological and physical attributes of each node and generate node feature vectors; The neural level set construction module is used to construct and train neural level set functions for each potential micro-fault edge to dynamically characterize the geometry, closure, and existence of the micro-fault. Edge network and message passing module: It is used to calculate the tomographic transmission coefficient through the edge network based on the node feature vector and the neural level set, and to perform message passing and update the node state according to the transmission coefficient, node pressure gradient and interface geometry information. Graph Neural Network and Decoding Module: Used to build graph neural networks, update the physical state of all nodes based on message passing, decode the updated node features, and output the predicted values ​​of node pressure and saturation. Inversion optimization module: used to construct a comprehensive loss function and use an optimizer to iteratively update all learnable parameters to achieve synchronous inversion of pressure-saturation field and microfault parameters.

Images