A fault leakage deduction method based on graph network and tree search

Abstract

This invention relates to the field of carbon geological storage technology and discloses a fault leakage inference method based on graph networks and tree search. The method includes: acquiring discrete fault network geological data; constructing a two-layer heterogeneous graph containing hydraulically connected edges and mechanically dispersed edges and performing initial masking; extracting features using a multi-head graph attention network to infer the prior probability of fluid movement and stress correction factors; constructing a Monte Carlo search tree with the injection node as the root node and introducing differentiated logs to record parameter changes; performing fluid frontal expansion in conjunction with the prior probability; when local pressure induces fault rupture, removing the masking and initiating asynchronous stress expansion to update the hydraulic and mechanical states; and extracting the node sequence with the highest access frequency as the dominant leakage channel after stopping the search. This invention avoids the computational burden of full-grid solution and significantly improves inference efficiency and accuracy while ensuring the physical fidelity of multi-field coupling.

Description

Technical Field

[0001] This invention relates to the field of carbon geological storage technology, specifically a fault leakage inference method based on graph networks and tree search. Background Technology

[0002] In carbon sequestration projects, injecting carbon dioxide into deep formations is a crucial means of reducing greenhouse gas emissions. Subsurface reservoirs typically possess complex, discrete fault networks, and the injected fluids alter formation pore pressure and local formation stress. If the fluid migrates along faults and breaches safety boundaries, it poses a risk to shallow groundwater and the surface environment. Accurately predicting the fluid migration path within complex fault networks and assessing leakage risks are fundamental to ensuring the safety of carbon sequestration projects.

[0003] For fault leakage evolution, multiphysics coupled numerical simulation methods with full mesh partitioning are commonly used. This method, based on the finite element method or finite volume method, obtains fluid migration trajectories by discretizing the study area into a three-dimensional spatial mesh and jointly solving the seepage and rock mechanics equations. However, when dealing with complex geological networks containing numerous fault segments, this method requires iterative solving of the global mesh stiffness matrix. The massive computational scale results in significant memory consumption and excessively long individual calculation times, making it difficult to meet the needs of rapid risk simulation under multiple pressure conditions in engineering projects.

[0004] To improve computational efficiency, some conventional techniques abstract fault networks into a basic graph data structure and use standard tree search or path planning algorithms to find fluid connectivity channels. These simplified methods mostly rely solely on initial permeability and spatial geometric intersections to infer fluid direction, failing to reflect the true fluid-structure interaction mechanisms within the formation medium. In actual physical evolution, the increase in pore pressure caused by fluid influx can induce fault slippage and rupture. The mechanical energy released by fault rupture can be transmitted across space through impermeable rock media to distant, unconnected faults, triggering cascading mechanical failures. Simultaneously, the compaction and abrasion effects of fault slip significantly reduce local permeability, forming physical blockages that impede further fluid migration. Conventional graph search models lack the ability to quantitatively track the asynchronous cross-space stress propagation, errors in local heterogeneous stress transmission, and dynamic secondary blocking effects. Due to the lack of underlying multiphysics interaction mechanisms, these simplified models easily overlook mechanically induced hidden leakage channels, leading to deviations between the inferred paths and actual geological responses, and failing to guarantee the accuracy of leakage risk assessment. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a fault leakage inference method based on graph networks and tree search. This method solves the problems of long computation time and high memory consumption in traditional full-mesh coupled numerical simulation methods, as well as the low accuracy of fault leakage path inference caused by the lack of fluid-structure interaction mechanism and inability to quantify asynchronous stress wave and dynamic secondary blocking effect in conventional graph search simplified models.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a fault leakage inference method based on graph networks and tree search, comprising:

[0007] Obtain discrete fault network geological data for the study area. The discrete fault network geological data includes the three-dimensional spatial geometric coordinates and initial physical parameters of the fault segments. The fault segments are mapped as nodes, hydraulic connectivity edges are constructed based on spatial intersection relationships, and mechanical remote diffusion edges are constructed based on static Coulomb stress transfer coefficient and preset wave sweep threshold. A two-layer heterogeneous graph is generated based on the hydraulic connectivity edges and mechanical remote diffusion edges, and a masking process is applied to the mechanical remote diffusion edges. The two-layer heterogeneous graph after applying a mask is input into a multi-head graph attention network, and the prior probability distribution of fluid movement between nodes and the dynamic sweep coefficient are output. Construct a Monte Carlo search tree by setting the node corresponding to the injection well location as the root node, and perform fluid front expansion along the hydraulically connected edge; When a node satisfies the Mohr-Coulomb rupture criterion, fluid front extension is suspended, the masking of the corresponding node's mechanical remote dispersion edge is removed, stress wave asynchronous front extension is initiated, and the physical parameters of the far-end connected node are updated using the dynamic sweep efficiency coefficient. After completing a single branch search, calculate the reward value of the search branch and perform backpropagation. Based on the differentiated log of the node modification operation, undo the changes in node pressure and penetration rate layer by layer to restore the initial physical state of the two-layer heterogeneous graph. When the maximum number of iterations is reached, the Monte Carlo search tree is traversed to extract the sequence of nodes with the highest access frequency, and the node sequence is output as the dominant leakage channel.

[0008] Furthermore, a two-layer heterogeneous graph is generated and a masking process is applied to the mechanically divergent edges, including: Assign each node a physical feature vector containing initial permeability, porosity, initial normal stress, initial shear stress, cohesion, friction coefficient, and initial formation pore pressure; The static Coulomb stress transfer coefficient between nodes is calculated using an elastic analytical model. When the absolute value of the static Coulomb stress transfer coefficient is greater than the preset swirl threshold, a directed edge is established between the nodes as a mechanical remote diffusion edge. The hydraulically connected edges and the mechanically dispersed edges are integrated to generate a two-layer heterogeneous graph, and the mechanically dispersed edges are assigned visibility flags with initial values ​​of false to achieve masking.

[0009] Furthermore, the prior probability distribution of the output fluid's movement between nodes and the dynamic sweep efficiency include: By utilizing a multi-head graph attention network to perform aggregation operations on node neighborhood features within the graph attention aggregation layer, a hidden layer representation that integrates local connected topology and heterogeneous attributes is extracted. The aggregated hidden layer representation is mapped to the prior probability distribution of the action of the fluid front moving to connected neighbor nodes using the classification prediction head; The regression head is used to output the non-homogeneous stress correction factor for the node pairs connected by the mechanical remote dispersion edge, and the non-homogeneous stress correction factor is multiplied by the pre-calculated static Coulomb stress transfer coefficient to generate the dynamic sweep coefficient.

[0010] Furthermore, the Monte Carlo search tree is constructed, including: The tree nodes of the Monte Carlo search tree contain a composite state view, which includes a fluid leading edge pointer, a local graph attribute view, a differential log, and a global flag snapshot. The global flag snapshot is a one-dimensional bitmap array used to record the mechanical steady-state flag state of each node; The differentiated log is a tuple containing the target node identifier, the changed attribute category, the original value, and the updated value. It is used to record the incremental modification operations of the underlying physical parameters caused by the current search step.

[0011] Furthermore, fluid front extension is performed along the hydraulically connected edge, including: By combining the prior probability distribution of actions and using the upper confidence bound algorithm of the predictor, the upper confidence bound of each node in the action candidate set is calculated, and the node with the largest upper confidence bound is selected as the expansion target node for fluid leading edge expansion. Calculate and update the pore pressure of the extended target node, and calculate the effective normal stress in combination with the initial normal stress of the target node; The effective normal stress is substituted into the Mohr-Coulomb fracture criterion for judgment. If the judgment is valid, the mechanical steady-state flag value of the expanded target node is modified to the locked state and the fluid leading edge expansion is suspended.

[0012] Furthermore, initiating asynchronous propagation of the stress wave and its leading edge includes: Change the visibility flag of the corresponding node's mechanical remote edge from false to true, dynamically remove the mask of the mechanical remote edge, and construct a local mechanical topology view. Calculate the static stress drop of the expanded target node that has fractured, and push the expanded target node as a disturbance source node into the asynchronous expansion queue; The induced Coulomb stress change at the affected connected node is calculated by multiplying the dynamic sweep efficiency by the static stress drop. Update the equivalent shear stress of the affected connected nodes and call the Mohr-Coulomb rupture criterion for judgment. If the rupture condition is met, push the affected connected nodes to the tail of the asynchronous expansion queue. Repeat this process until the asynchronous expansion queue is empty, then terminate the asynchronous expansion process and resume fluid leading edge expansion.

[0013] Furthermore, the process of initiating the asynchronous propagation of stress waves and leading edges also includes performing secondary blocking determination and cascading cutoff control: The permeability multiplier method was used to calculate the updated permeability values ​​of the wave-affected and connected nodes that experienced slip rupture. If the updated permeability value is lower than the preset critical closure permeability threshold, it is determined that secondary blockage has occurred, and the hydraulic connection edge connected to the affected connected node is set to a temporary shielded state. Before the affected and connected nodes are pushed into the asynchronous expansion queue, the propagation depth and the induced Coulomb stress change are checked. If the propagation depth reaches the preset cascade depth threshold, or the absolute value of the induced Coulomb stress change is less than the preset minimum stress cutoff threshold, the operation of queuing the affected and connected nodes is forcibly discarded, and the mechanical transmission chain is cut off.

[0014] Furthermore, the reward value of the search branch is calculated and backpropagation is performed, including: Extract the multiphysics state parameters corresponding to the leaf nodes, and calculate the comprehensive reward value of the current leaf node by combining the total amount of injected fluid, the number of fault nodes that have slipped and ruptured, and the Euclidean distance between the fluid front and the preset safety boundary through normalization and weighted summation. The algorithm backtracks from bottom to top to the root node of the Monte Carlo search tree, updates the node access frequency and historical access frequency of the action branch of the tree node on the backtracking path, and uses the incremental mean update formula to iteratively calculate the expected value of the action.

[0015] Furthermore, by reversing the changes in node pressure and penetration rates layer by layer according to the differentiated logs, the initial physical state of the two-layer heterogeneous graph is restored, including: Based on the depth of the Monte Carlo search tree, log popping operations are performed layer by layer from the leaf node to the root node, and within a single node, operation records in the differentiated logs are extracted according to the last-in-first-out reverse order rule. The physical parameters corresponding to the corresponding node in the operation record are overwritten from the updated values ​​back to the original values; Change the state of the node that was set to 1 in the current branch in the global flag snapshot back to 0, and undo the Coulomb slip and cascade breakage events that occurred in the current branch.

[0016] Furthermore, after outputting the node sequence as the dominant leakage channel, the process also includes an indicator calculation step: Extract the geological and geometric parameters and hydraulic properties of adjacent nodes in the dominant leakage channel, and calculate the hydraulic conductivity between adjacent nodes; The flow distribution weight of each downstream node is determined by the ratio of the hydraulic conductivity of the corresponding branch to the total conductivity. The actual volumetric flow rate allocated to each downstream node is calculated by combining the initial total flow rate of the injection well with the flow rate allocation weight. The time increment of fluid breakthrough at each downstream node is calculated based on the actual volumetric flow rate and the effective pore volume of the corresponding node, and the time increment along the dominant leakage channel is accumulated step by step to generate a leakage time series curve.

[0017] This invention provides a method for fault leakage prediction based on graph networks and tree search. It has the following advantages: 1. This invention constructs a two-layer heterogeneous graph containing hydraulically connected edges and mechanically dispersed edges, and introduces a local graph attribute view and a differentiated log data structure. During the multi-branch exploration phase of tree search, the system does not copy the full basic geological map data, but records incremental modifications of physical parameters such as pressure and permeability through differentiated logs, and performs reverse rollback during search backtracking. The above technical solution effectively isolates the multi-field coupled physical states between independent search branches, reduces the memory allocation overhead during multi-branch parallel computing, avoids the computational burden of traditional full-grid matrix solving, and improves the overall inference efficiency of the system.

[0018] 2. This invention utilizes a multi-head graph attention network to extract the local connectivity topology and heterogeneity attributes of discrete fracture networks, outputting the action prior probability distribution and the heterogeneity stress correction factor. The action prior probability is introduced into the action selection formula of the tree search, providing heuristic guidance for the Monte Carlo search and reducing invalid search branches. Simultaneously, the heterogeneity stress correction factor is used to correct the static Coulomb stress transfer coefficient calculated by the elastic analytical model, compensating for the bias of conventional theoretical models when dealing with heterogeneous rock masses, and improving the hit rate and accuracy of dominant leakage path prediction.

[0019] 3. This invention designs a cross-field deduction mechanism for the asynchronous propagation of fluid front and stress wave front, and, in conjunction with the dynamic mask removal operation of mechanical remote dispersion edge, restores the true fluid-structure interaction evolution sequence. When local pressure rise induces fault rupture, the system suspends fluid propagation, initiates stress-dominated asynchronous wave calculation, and simultaneously performs hydraulic property updates and secondary blockage determination. Through the above asynchronous wave and blockage determination, the system can accurately map the local physical sealing effect caused by fault cascade rupture and fault slip compaction in the graph search architecture, ensuring the physical fidelity of leakage risk assessment in complex geological storage scenarios. Attached Figure Description

[0020] Figure 1 This is a diagram illustrating the deductive system architecture of the present invention; Figure 2 This is a flowchart of the derivation method of the present invention; Figure 3 This is a comparison chart of the computational efficiency of the present invention; Figure 4 This is a comparison chart of the prediction accuracy of the present invention. Detailed Implementation

[0021] 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.

[0022] See attached document Figure 1 , Figure 1 This is a system architecture diagram for inference according to an embodiment of the present invention. The present invention provides a fault leakage inference system based on graph networks and tree search, deployed in a computer device equipped with a processor and memory. The inference system includes a data acquisition module, a graph initialization module, a network inference module, a tree search module, and a result analysis module.

[0023] See attached document Figure 2 , Figure 2 This is a flowchart of a deduction method according to an embodiment of the present invention. The present invention provides a tomography leakage deduction method based on graph networks and tree search, comprising the following steps: S10, the data acquisition module acquires the discrete fault network geological data of the study area. The discrete fault network geological data includes the three-dimensional spatial geometric coordinates of the fault segments, the initial physical parameters, and the preset safety boundary coordinates for risk assessment. S20, the graph initialization module receives the geological data of the discrete fracture network, maps the fault segments to nodes in the graph data structure, and constructs hydraulically connected edges based on the spatial intersection relationship of the fault segments; S30, the graph initialization module uses the elastic analytical model to calculate the static Coulomb stress transfer coefficient between nodes, constructs the mechanical remote divergence edge based on the sweep threshold, generates a two-layer heterogeneous graph based on the hydraulic connectivity edge and the mechanical remote divergence edge, and applies a masking process to the mechanical remote divergence edge. S40, the network inference module reads the masked double-layer heterogeneous graph, inputs it into the multi-head graph attention network for feature aggregation, and outputs the prior probability distribution of the fluid's movement between nodes. S50, the network inference module uses the regression head of the multi-head graph attention network to output the inhomogeneous stress correction factor, and multiplies the inhomogeneous stress correction factor with the static Coulomb stress transfer coefficient to generate the dynamic sweep coefficient. S60, the tree search module receives the action prior probability distribution and the two-layer heterogeneous graph, sets the node corresponding to the injection well position as the root node, and constructs a Monte Carlo search tree; S70, the tree search module performs fluid front extension along the hydraulically connected edge, calculates the pore pressure increment and effective normal stress of the newly visited node, and suspends fluid front extension when the effective normal stress satisfies the Mohr-Coulomb rupture criterion; S80, the tree search module removes the masking of the mechanical remote dispersion edge of the corresponding node, starts the asynchronous expansion of the stress wave front edge, calculates the effective stress disturbance of the far-end connected node using the dynamic wave coefficient, and updates the permeability value of the far-end connected node. S90, after the tree search module completes a single branch search, it calculates the reward value of the search branch and performs backpropagation. During the backpropagation process, it cancels the changes in node pressure and penetration rate layer by layer according to the differentiated log, and restores the initial physical state of the two-layer heterogeneous graph. S100, the result analysis module stops searching when the set maximum number of iterations is reached, traverses the Monte Carlo search tree to extract the node sequence with the highest access frequency, outputs the node sequence as the dominant leakage channel and calculates the flow allocation weight and propagation time distribution.

[0024] In a specific implementation, the graph initialization module receives discrete fracture network geological data transmitted by the data acquisition module, specifically including the following steps: The graph initialization module performs geometric discretization mapping on fault segments in a discrete fracture network. A discrete fracture network is typically composed of fault planes distributed in three-dimensional space. The graph initialization module maps each fault segment in the discrete fracture network to a node in a graph data structure. The system sets the generated set of fault segment nodes as follows: ,in This represents the total number of fault segments. Represents the set of fault segment nodes The first in Each node. For the discretization and coordinate extraction of the fault plane geometry, as well as the calculation and storage of the corresponding fault grid physical volume, those skilled in the art can use conventional geological modeling grid partitioning algorithms to complete the task. The related operations are well-known technologies in this field and will not be described in detail here.

[0025] The graph initialization module is a set of fault segment nodes. The graph initialization module assigns physical feature vectors to nodes. These physical feature vectors characterize the initial geomechanical state of the fault segment in multi-field coupled simulations. Considering the normalization requirements of the input data tensor dimension for subsequent multi-head graph attention networks, the graph initialization module extracts the geometric center point location attribute of the fault segment, using the parameters at that center point location to represent the overall physical state of the corresponding node. The physical feature vector is configured as a seven-dimensional array structure. Specifically, this vector consists of the initial permeability, porosity, initial normal stress, initial shear stress acting on the fault plane, cohesion of the rock at the fault plane, friction coefficient, and initial formation pore pressure corresponding to the fault segment. At the physical principle level, permeability and porosity determine the basic conduction conditions for fluid transport between nodes, while the anisotropic stress parameters and friction coefficient jointly determine the mechanical stability of the fault plane. To avoid interference from differences in physical dimensions in subsequent feature aggregation calculations, the graph initialization module standardizes the initial physical parameter values ​​input from the data acquisition module and then writes the processed values ​​into the corresponding dimensions of the aforementioned physical feature vector.

[0026] The graph initialization module instantiates mechanical steady-state flags in the node data structure. In rock strata physics, after stress release and slippage occur in a fault segment, the local area enters an unloading steady state. It is generally believed that the fault segment is unlikely to undergo secondary fracturing within the short duration of this unloading steady state. Based on this principle, the graph initialization module initializes each node... Add a separate data field to store the mechanical steady-state flag. This flag is used to restrict any node from repeatedly acting as a stress spillover source during algorithm iteration, preventing the computational logic from falling into a deadlock state of infinite cascading triggers.

[0027] The graph initialization module initializes the mechanical steady-state flag. The data type is set to integer variable. Its value range is preferably limited to: A value of 0 indicates that the corresponding node is in an initial bearing state without slippage, while a value of 1 indicates that the corresponding node has triggered Mohr-Coulomb fracture and is in a stress-locked state. The initialization module generates the fault segment node set. At that time, traverse all nodes. and set all mechanical steady-state flags The initial value is set to 0. Through the above structure definition and parameter filling, the graph initialization module completes the memory allocation and initialization of heterogeneous graph nodes, including physical feature parameters and state control words, providing structured data input for subsequent feature extraction of the graph network.

[0028] In this embodiment, after completing the memory allocation for the nodes, the graph initialization module continues to read the structural information from the geological data to establish topological relationships between the nodes. To accurately characterize the physical differences between fluid transport and stress transmission in the fault network, the graph initialization module constructs a two-layer heterogeneous edge structure incorporating both hydraulic and mechanical connectivity characteristics. The specific implementation steps are as follows: The graph initialization module constructs hydraulically connected edges based on the spatial geometric intersections of fault segments in the discrete fault network. In actual geological environments, fluid transport within the fault network depends on the physical connectivity between fault planes. Specifically, the graph initialization module traverses the set of fault segment nodes. All nodes in the set, when any two nodes and When the corresponding fault segments intersect or coincide at their endpoints in three-dimensional space, the graph initialization module initializes the nodes. and Establish undirected edges between them, thereby generating a set of hydraulically connected edges. The aforementioned undirected edges indicate that, in the presence of a pressure gradient difference, fluid can conduct bidirectionally between these intersecting fault segments. The determination of three-dimensional geometric intersections can be accomplished by those skilled in the art using conventional bounding box collision detection algorithms; the related operations are well-known techniques in this field and will not be elaborated upon here.

[0029] Considering the complex multi-field coupling characteristics of the formation medium, the local stress disturbance generated by fault slip is not limited to intersecting faults but can also be transmitted to distant, unconnected fault segments through the surrounding impermeable rock mass. To characterize this cross-spatial mechanical coupling effect, the graph initialization module introduces a mechanism for constructing mechanical remote dispersion edges. The graph initialization module uses an elastic half-space analytical model (such as the commonly used Okada analytical model) to calculate the static Coulomb stress transfer coefficient between nodes. This coefficient characterizes the change in Coulomb fracture stress on the fault plane of the target node when the fault where the source node is located slips by a unit area. In specific implementation, the graph initialization module targets the fault segment node set... Any two non-overlapping nodes and Geometric positional parameters such as their three-dimensional spatial distance, relative azimuth, fault strike, and dip angle are extracted and substituted into the elastic half-space analytical model as input data to obtain the static Coulomb stress transfer coefficient. .

[0030] The graph initialization module sets a sweep threshold based on the static Coulomb stress transfer coefficient to generate mechanically divergent edges. Since fully connected graphs often result in excessively high tensor dimensions during subsequent graph network aggregation calculations and contain numerous physically meaningless weak perturbations, the graph initialization module configures a sweep threshold within the system. The sweep threshold is typically determined based on the minimum measurable stress accuracy obtained from geological exploration within the study area; for example, its value can be set between 0.01 MPa and 0.05 MPa. When the static Coulomb stress transfer coefficient... The absolute value is greater than the spurious threshold. At that time, the graph initialization module is at the node With nodes Establish a link between them point to The directed edges that satisfy this criterion. All directed edges that satisfy this criterion together constitute the set of mechanical remote edges. The introduction of directed edges reflects the unidirectional propagation characteristic of stress release sources transferring mechanical energy to specific target areas in the surrounding region. This helps to eliminate weak connections that transfer energy, thereby reducing the computational load of the system.

[0031] To form a complete computational topology, the graph initialization module integrates the aforementioned topological relationships to generate a two-layer heterogeneous graph and configures an initial masking mechanism. In this embodiment, the graph initialization module sets the hydraulically connected edges... With mechanical remote edge set Integrate them into the same data framework to generate a two-level heterogeneous graph. In the initial operational phase, before conventional formation pressure triggers fault mechanical slippage, the mechanical stress transmission channels are inactive. To achieve physical isolation at the algorithm's underlying layer, the graph initialization module initializes the set of mechanically dispersed edges. A global masking process is applied. This process is performed by storing the mechanical remote edge set in computer memory. An additional Boolean visibility flag is allocated to achieve this. The graph initialization module sets the initial value of this visibility flag to false. When the subsequent network inference module reads the two-layer heterogeneous graph for data feature extraction, due to the limitation of this visibility flag, the network inference module can only identify and read the set of hydraulically connected edges at this time. The topology connections within the system. This masking process is designed to enable the simulation system to dynamically switch between different topology views based on actual physical triggering conditions.

[0032] In this embodiment, the network inference module undertakes the core computational task of joint feature inference. Considering the accuracy limitations of analytical solutions in elasticity when dealing with heterogeneous rock masses, the network inference module introduces a deep learning architecture to achieve joint extraction and nonlinear mapping of multi-physics features. Its specific inference process includes the following steps: The network inference module constructs a multi-head graph attention network architecture and performs data preprocessing. This multi-head graph attention network consists of an input layer, two to three parallel graph attention aggregation layers, and dual-branch output head connections. The network inference module uses a masked, two-layer heterogeneous graph as input data. Because the set of mechanically dispersed edges is masked by visibility flags, the data entering the network essentially degenerates into a locally connected view containing only hydraulically connected edges. To eliminate dimensional differences between different geophysical quantities (such as permeability and formation pressure), the network inference module performs max-min normalization on each dimension of the physical feature vector at the input layer. After preprocessing, the input tensor is mapped to a standard distribution seven-dimensional data format to receive basic physical parameters such as initial permeability, porosity, and stress of each fault segment node.

[0033] The network inference module performs aggregation operations on the neighborhood features of nodes within the graph attention aggregation layer. The permeability and conduction capacity of fluid in a fault depends not only on the physical properties of the nodes themselves but also on the potential energy gradient of surrounding connected nodes. To capture this spatial dependency, the network inference module uses a multi-attention mechanism to assign different attention weights to each neighboring node on a hydraulically connected edge. From a physics perspective, these attention weights dynamically reflect the tendency of fluid to flow towards a specific neighboring node under the influence of local porosity and permeability differences. The network inference module calculates the feature interaction matrix between nodes using multiple independent attention heads and maps the concatenated feature vectors using a nonlinear activation function, thereby extracting a hidden layer representation that integrates local connected topology and heterogeneous properties. For the basic weight allocation and gradient update process of the multi-head attention mechanism, those skilled in the art can refer to the conventional forward propagation implementation in standard deep learning frameworks; its underlying computational logic is well-known in the field and will not be elaborated upon here.

[0034] After multi-layer feature aggregation, the network inference module generates inference results through a dual-branch output head. One branch is a classification prediction head, which the network inference module uses to output the action prior probability distribution in a multi-head graph attention network. This classification prediction head typically consists of fully connected layers and a normalized exponential function (Softmax function). Its function is to map the aggregated feature vectors to the probability distribution values ​​of the fluid front moving towards each connected neighbor node in the next time step. In practical business scenarios, this action prior probability distribution intuitively reflects the physical trend of which specific fault path the fluid tends to preferentially propagate along under the current local pressure gradient and geostress configuration. This probability distribution will serve as an important heuristic guide for the subsequent tree search module's exploration direction.

[0035] While outputting the fluid conduction probability, the network inference module simultaneously outputs a heterogeneous stress correction factor using the regression head in the multi-head graph attention network. Conventional elastic half-space analytical models are based on the ideal assumption of homogeneous and isotropic media, making it difficult to reflect stress transmission deviations caused by abrupt changes in lithofacies within complex fault zones. To correct this theoretical error, a fully connected layer with a linear or sigmoid activation function is configured inside the regression head. For node pairs connected by mechanically divergent edges, a dimensionless scalar characterizing the degree of heterogeneity is output—the heterogeneous stress correction factor. The value of this correction factor can be limited to the range of [0.5, 1.5] to represent a moderate fine-tuning of the theoretical stress analytical value. Subsequently, the network inference module will apply the non-homogeneous stress correction factor... Compared with the pre-calculated static Coulomb stress transfer coefficient Multiply to generate the dynamic sweep efficiency. The dynamic sweep efficiency coefficient integrates the baseline physical constraints of the macroscopic analytical model with the local heterogeneous perturbation features extracted from geological data by the graph neural network, which to some extent helps to improve the reliability of the quantitative calculation of mechanical sweep efficiency effects.

[0036] To ensure the multi-head graph attention network possesses strong inference capabilities, the network inference module requires offline model training before online inference. The module acquires historical microseismic monitoring records of the study area or uses high-precision calculation results generated by finite element multi-field coupled numerical simulation software as training samples. The sample data includes node feature matrices under various injection pressure conditions. The corresponding label data consists of two parts: first, the probability distribution of fluid-dominant leakage paths calculated based on the actual fluid flux ratio between nodes in the numerical simulation; second, the ratio of the actual effective stress increment caused by fault slip extracted from the actual model to the value calculated by the static analytical model, serving as the regression label for the stress correction factor. During training, the network inference module uses a cross-entropy loss function to optimize the classification prediction head to fit the prior probability of the action, and a mean squared error loss function to optimize the regression head to reduce the deviation between the predicted dynamic sweep efficiency and the actual stress disturbance extreme value. The network inference module linearly adds the two loss functions with set weight coefficients to form a joint loss function, and iteratively updates the network's internal parameters using a backpropagation algorithm until the joint loss function converges to a set tolerance range, completing the model construction and parameter solidification.

[0037] In this embodiment, after the network inference module completes the inference of prior probabilities and dynamic sweep efficiency, the tree search module executes the subsequent tree search deduction process. Fluid transport within a fault network causes dynamic changes in the local pressure and stress environment. To adapt to this dynamic evolution requiring multi-field coupling, the tree search module improves the underlying data structure based on the conventional tree search algorithm, specifically including the following steps: The tree search module receives the action prior probability distribution and the two-layer heterogeneous graph, sets the node corresponding to the injection well location as the root node, and begins constructing the Monte Carlo search tree. During construction, the tree search module configures a composite state view within each node of the Monte Carlo search tree. Specifically, any tree node in the Monte Carlo search tree... In addition to storing the regular node access frequency Action value and the prior probability of the action output by the network inference module In addition, it is configured to include four data components: fluid leading edge pointer, local graph attribute view, differential log, and global flag snapshot. When initializing the root node, the tree search module points the fluid leading edge pointer to the physical node number where it is injected and initializes an empty differential log list.

[0038] To achieve dynamic tracking of spatial position and mechanical state, the tree search module is used at tree nodes. The code defines a fluid leading edge pointer and a global snapshot flag. The fluid leading edge pointer records the specific physical node number in the current simulation where the fluid propagates to the bilayer heterogeneous graph. This pointer physically reflects the current spatial exploration depth of this search branch. Meanwhile, stress release caused by fault slip is irreversible. The tree search module is located at tree nodes. A global flag snapshot is instantiated to record the mechanical steady-state flag status of all fault segment nodes under the current search path. To optimize memory usage, the tree search module configures the global flag snapshot as a one-dimensional bitmap array. The array index corresponds to the node number in the fault segment node set, and the bitmap values ​​0 and 1 correspond to the non-slip and locked states, respectively. The global flag snapshot ensures that subsequent stress waves and leading edges can accurately identify valid target nodes that have not been locked during asynchronous expansion.

[0039] During fault leakage simulation, multiple exploration branches will generate their own independent pressure and stress distribution states. Completely copying the basic geological map data with each branch expansion would significantly increase computer memory consumption and reduce search efficiency. To address the data redundancy issue caused by multi-branch exploration, the tree search module employs a local graph attribute view mechanism. (Tree nodes...) Instead of actually storing the complete two-layer heterogeneous graph structure data, it holds a static memory reference to the underlying basic two-layer heterogeneous graph. During a single branch search, the inference system's reading operations of basic geological environment parameters are all mapped to the globally unique underlying graph data structure through this local graph attribute view.

[0040] To complement the local graph attribute view and achieve physical state isolation between search branches, the tree search module uses tree nodes. The system defines a differential log data structure. The differential log records incremental modifications to the underlying physical parameters caused by the current search step. The tree search module configures the differential log data format as a tuple containing four dimensions: target node identifier, changed attribute category, original value, and updated value. For example, when fluid expansion causes an increase in the pore pressure of a physical node, the differential log will generate an operation record containing the node number, pore pressure identifier, pressure value before the change, and pressure value after the change. In terms of data reading logic, when the tree search module queries the physical attributes of a node through the local graph attribute view, the system prioritizes searching the differential log on the current branch path. If a change record for that node exists in the log, the system reads and returns the updated value; otherwise, the system reads the initial parameters from the underlying graph data structure through static memory references. Through this incremental recording method, the tree search module can track the dynamic evolution trajectory of complex multi-field parameters with less memory consumption when multiple search branches are running in parallel, providing data support for subsequent state restoration.

[0041] In this embodiment, after constructing a composite state view of the Monte Carlo search tree, the tree search module selects an exploration direction from the currently feasible fluid paths based on the predictor's confidence bound algorithm, and then performs fluid front extension calculations. To simulate the gradual transport of fluid in complex porous media, this extension process specifically includes the following steps: The tree search module selects actions based on prior guidance from the multi-head graph attention network. It reads the physical node pointed to by the fluid leading edge pointer in the current tree node and retrieves all adjacent unvisited nodes of that physical node under the hydraulic connectivity edge set through the local graph attribute view, generating an action candidate set. In actual fault connectivity networks, the fluid leading edge may advance to a blind end with no subsequent path. The tree search module pre-determines whether the action candidate set is empty; if the candidate set is empty, it indicates that the search branch has reached the local physical boundary, and the tree search module suspends the current fluid expansion process and directly enters the subsequent reward calculation stage. If the candidate set is not empty, for each candidate node in the action candidate set, the tree search module calculates its corresponding confidence upper bound value based on the action prior probability output by the network inference module. To accurately quantify the balance between exploration and exploitation, the tree search module uses the predictor upper confidence bound algorithm to calculate the confidence upper bound value of each node in the action candidate set. The specific calculation formula is configured as follows: ; In the formula, Indicates the current state of the tree node. Select action branch below The upper confidence limit; The expected value of the current action in this branch; A constant used to control the degree of exploration; The prior probability of the corresponding action is output by the network inference module; The current tree node Total access frequency; This represents the historical access frequency of this specific action branch. The constant controlling the exploration level can be configured between 1.0 and 2.0 depending on the fault mesh size. The tree search module traverses the action candidate set and selects the node with the largest confidence upper limit as the target node for fluid front extension.

[0042] The tree search module performs spatial transfer of the fluid front and calculates pressure field disturbances. After determining the target node, the tree search module updates the fluid front pointer of the current tree node to the target node number, indicating that the fluid migration front has reached a new position. The injection of fluid disrupts the initial hydrostatic equilibrium of the target node. Based on preset injection and flow rate parameters and the porosity and permeability of the fault segment where the target node is located, the tree search module calculates the corresponding pore pressure increment. Based on this increment, the tree search module calculates the updated pore pressure of the target node. To maintain the purity of the bottom map data, the tree search module writes the pressure change record of the node into a differential log, which sequentially includes the target node number, pore pressure attribute identifier, initial formation pore pressure value, and updated pore pressure value. For the specific method of calculating the pore pressure increment based on the continuity equation or empirical formulas of seepage mechanics, those skilled in the art can use conventional reservoir engineering algorithms, and the related operations are well-known techniques in the field, and will not be elaborated here.

[0043] The tree search module performs a mechanical state check on the target node. In the principle of multi-field coupling in rock formations, an increase in formation fluid pressure can offset part of the normal stress acting on the fault plane, thus reducing the shear strength of the fault. The tree search module extracts the initial normal stress of the target node and calculates the effective normal stress using the updated pore pressure. This effective normal stress is numerically equal to the initial normal stress minus the updated pore pressure. After completing the stress conversion, the tree search module substitutes this effective normal stress into the Mohr-Coulomb fracture criterion for judgment. The tree search module uses the Mohr-Coulomb fracture criterion to construct a specific mechanical discrimination formula, which is expressed as: ; In the formula, The initial shear stress borne by the fracture surface of the target node; The cohesive force of the rock at the fault plane; The friction coefficient of the fracture surface; The initial normal stress acting on the fault plane; The pore pressure is renewed after fluid filling; This characterizes the effective normal stress after the fluid pressure is offset. Logically, if the above condition is not met, it indicates that the current fluid injection only changes the pressure field and does not induce mechanical instability of the rock mass structure. The tree search module continues to select actions in the next round along the hydraulically connected edge. If the above condition is met, it indicates that the increase in pore pressure causes the effective normal stress to drop to a critical level, and the fault triggers slip fracture at this location.

[0044] The tree search module responds to the rupture event and suspends the current fluid front extension process. Because fault slip radiates elastic strain energy outward, the transfer rate of this mechanical energy is typically greater than the infiltration conduction rate of fluid in the micropores of rock. To align with the real-world multiphysics evolution timeline, the tree search module pauses the advancement of the fluid front pointer after detecting that the rupture condition has been met. The tree search module modifies the mechanical steady-state flag of the target node from 0 to 1 and synchronously updates this state change to the global flag snapshot and the differential log. This suspension operation establishes the timing node for the transition from the fluid field to the stress field, providing the triggering condition for subsequent algorithms to activate the mechanical remote divergence edge and initiate the asynchronous stress extension mechanism.

[0045] In this embodiment, after detecting fault slip induced by a rise in local pore pressure, the tree search module switches the simulation system from a fluid-dominated unidirectional expansion mode to a stress-dominated asynchronous sweep mode in order to accurately reconstruct the cascading failure effect caused by stress release. The tree search module performs the following specific steps for this physical process: S420, the tree search module dynamically removes the masking of the mechanical remote divergence edges, activating the local mechanical topology view. After the fluid leading edge node triggers Mohr-Coulomb rupture in the preceding steps, this node becomes a stress disturbance source. The tree search module accesses the underlying two-layer heterogeneous graph data structure and changes the visibility flag of all mechanical remote divergence edges originating from this stress disturbance source node from false to true. Through this underlying memory flag overwrite operation, the tree search module constructs a local mechanical connected subgraph for this rupture node, enabling the deductive logic to proceed with subsequent exploration along the mechanical transmission path.

[0046] The tree search module calculates the stress drop of the disturbance source and initializes an asynchronous extended queue. In rock mechanics, the mechanical energy released by fault slip depends primarily on the stress drop at the moment of rupture. The tree search module extracts the initial shear stress of the fault plane at the disturbance source node and calculates the residual shear stress after rupture by multiplying the current effective normal stress of the node by a set dynamic friction coefficient. The value of this dynamic friction coefficient is usually smaller than the static friction coefficient, and its specific value can be calibrated using rock mechanics experimental data. Based on this, the tree search module uses the difference between the initial shear stress and the residual shear stress as the static stress drop. To handle possible multi-node cascading ruptures, the tree search module instantiates an asynchronous extended queue of a first-in-first-out data structure in computer memory and pushes the current disturbance source node number into this queue as the initial sweep point.

[0047] The tree search module calculates the induced Coulomb stress change at the target node using the dynamic sweep efficiency coefficient. The module pops a node from the head of the asynchronous expansion queue as the current source node and traverses all its adjacent target nodes along the activated mechanical divergence edges. For each target node, the tree search module, combined with the dynamic sweep efficiency coefficient pre-output by the network inference module, calculates the induced Coulomb stress change caused by the slippage of the current source node on the target node. In practice, the increment of Coulomb fracture stress borne by the target node is numerically equal to the product of the dynamic sweep efficiency coefficient from the source node to the target node inferred by the network inference module and the static stress drop released by the source node upon fracture. This calculation process applies the inference results of the neural network to the quantitative assessment of physical propagation.

[0048] The tree search module updates the mechanical state of the target node and checks for cascade failure. For a target node that receives a stress disturbance, the tree search module reads a snapshot of the global flag bit in its current tree node. If the mechanical steady-state flag bit value of the target node in the snapshot is 1, it indicates that the node has already slipped in the current search branch and is in a stress unloading locked state. To avoid the algorithm getting stuck in a deadlock loop triggered by two nodes, the tree search module will skip this node. If the flag bit value is 0, the tree search module adds its initial shear stress to the calculated induced Coulomb stress change, updates its equivalent shear stress, and writes this change record to the differentiation log. After completing the numerical update, the tree search module calls the Mohr-Coulomb failure criterion again to determine the mechanical instability of the target node. If the determination result is that the failure condition is met, the tree search module modifies its mechanical steady-state flag bit to 1, synchronously updates the global flag bit snapshot and the differentiation log, and pushes the target node to the tail of the asynchronous expansion queue as a secondary stress disturbance source for the next level.

[0049] The tree search module terminates the asynchronous expansion process and resumes fluid front expansion. The tree search module repeatedly executes the queue popping and checking logic until the asynchronous expansion queue is empty. An empty queue physically indicates that the strain energy released by the local fault slip has been absorbed or attenuated by the surrounding rock mass, preventing further rupture events, and a new transient equilibrium has been established in the stratigraphic mechanical field. At this point, the tree search module ends the mechanically driven asynchronous ripple flow and unsuspends the previously paused fluid front expansion process. Based on the updated local graph attribute view and differentiated logs, the system continues the Monte Carlo tree search action selection for the next time step along the hydraulic connectivity edge. Through this cross-field asynchronous interaction mechanism, the system, to a certain extent, helps to map the spatiotemporal evolution sequence of complex multiphysics fields without increasing the time slice density.

[0050] In this embodiment, the tree search module, during the asynchronous expansion of stress waves and leading edges, not only needs to track changes in the mechanical state but also needs to handle the nonlinear influence of fault slip on local fluid conduction capacity. To prevent computational bottlenecks caused by loop non-convergence that may occur under complex network topologies and to accurately reflect the physical blocking effect caused by fault slippage, the tree search module introduces a secondary blocking determination and anti-cascading truncation control mechanism. The specific implementation steps are as follows: The tree search module updates the hydraulic properties of target nodes that have experienced slip fractures. In the intersection of rock mechanics and seepage mechanics, shear slip at fault planes can grind and generate fault gouge or alter the original micro-protrusion support structure, thereby causing changes in local permeability and porosity. When the tree search module determines that the fracture condition of a node is met in the asynchronous expansion queue and modifies its mechanical steady-state flag to 1, the system synchronously triggers the correction of the hydraulic properties of that node. The tree search module uses the permeability multiplier method to extract the current initial permeability of the target node and multiplies it by a preset slip permeability multiplier to calculate the updated permeability value. In practical applications, this slip permeability multiplier is usually configured according to the lithological characteristics of the target reservoir. For example, for faults rich in clay minerals, the compaction effect caused by slip will significantly reduce conductivity, and the multiplier can be between 0.01 and 0.1. After the calculation is completed, the tree search module writes the permeability change record to the differential log so that the latest seepage parameters can be called when the fluid front expands subsequently. For the specific numerical conversion of permeability change based on the stress-seepage coupling constitutive equation, those skilled in the art can use conventional reservoir numerical simulation theoretical formulas to complete it. The relevant operations are well-known techniques in this field and will not be elaborated here.

[0051] The tree search module performs secondary blockage determination based on updated hydraulic properties. Under specific geological conditions, the compaction effect caused by fault slippage can significantly reduce local permeability, forming a physical barrier that hinders further fluid migration. To map this phenomenon into the graph search model, the tree search module configures a critical closure permeability threshold within the system. This threshold can be set to 1×10-1. -16 m 2 (Approximately 0.1 millidarcy) is used to characterize the boundary condition where the rock mass loses its effective permeability. The tree search module compares the updated permeability value with this critical closure permeability threshold. If the updated permeability is lower than this threshold, it indicates that secondary blockage has occurred in the fault segment. The tree search module, through the local graph attribute view, configures all hydraulically connected edges connected to this node to a temporary shielded state. This operation cuts off the possible path of fluid through this node at the algorithm's underlying level, thereby preventing the subsequent upper confidence bound algorithm of the predictor from including it in the effective action candidate set during action selection.

[0052] The tree search module implements anti-cascading truncation control in the asynchronous expansion queue. Because underground fault networks may have highly dense topological connections, local stress release could trigger a large-scale chain reaction of ruptures, causing the algorithm to fall into a resource-exhausting logical operation loop. To ensure the computational robustness and convergence of the simulation system, the tree search module sets a cascading depth threshold and a minimum stress truncation threshold. The cascading depth threshold can be set to an integer between 5 and 10 to limit the continuous triggering levels of secondary stress disturbance sources; the minimum stress truncation threshold is set to a specific value between 0.01 MPa and 0.05 MPa based on the required computational accuracy. Before pushing a new node into the asynchronous expansion queue, the tree search module verifies the current cascading propagation depth and the induced Coulomb stress change. In the specific algorithm implementation, the tree search module appends a depth label to each node upon enqueueing; this label value is equal to the depth label of its triggering source node plus one, thus achieving accurate tracking of the propagation depth. If the propagation depth reaches the cascade depth threshold, or if the absolute value of the induced Coulomb stress change calculated is less than the minimum stress cutoff threshold, the tree search module will forcibly discard the enqueue operation of the node and directly cut off the mechanical transmission chain of the branch.

[0053] The tree search module solidifies the overall state of the current node and advances the global search rounds. After the aforementioned secondary blocking judgment and truncation control, the differentiated logs and global flag snapshots within the current tree node have effectively recorded the evolution results of the fluid field and mechanical field interaction triggered by a single fluid injection. The tree search module solidifies these updated composite state views in the current search branch. Through this mechanism, the system utilizes local asynchronous loops to handle the transient changes of multi-field coupling during unidirectional Monte Carlo time-step advancement, enabling the inference logic to objectively reflect the complex dynamic response characteristics of the geological medium while maintaining the algorithm's deadlock prevention security.

[0054] In this embodiment, after the tree search module expands to the boundary at the fluid leading edge or reaches the preset maximum search depth, it needs to evaluate the value of the current exploration results and transmit the information back along the search tree. To ensure that the next exploration branch can operate based on the accurate underlying physical environment, the system needs to roll back the underlying parameters using the previously recorded differential logs. This process specifically includes the following steps: The tree search module calculates the comprehensive reward value of the current leaf node. When a search branch triggers a termination condition, the tree search module extracts the multiphysics state parameters corresponding to that leaf node and evaluates the merits of the current path using a preset reward value function. This reward value function comprehensively considers the total amount of injected fluid, the number of fault nodes experiencing slip rupture, and the Euclidean distance between the fluid front and the preset safety boundary. In the specific evaluation logic, the tree search module configures the total amount of injected fluid as a positive reward term to incentivize the algorithm to find injection paths with greater capacity; simultaneously, it configures the number of fault nodes experiencing slip rupture and the degree to which the fluid front approaches the safety boundary as negative penalty terms to constrain the safety boundary of the deduction process. After normalizing the above three indicators, the tree search module assigns corresponding weight coefficients and performs a weighted summation. In conventional geological storage safety assessment scenarios, the weight of safety indicators is usually higher than that of capacity indicators; for example, the weight of the rupture penalty term can be set to 0.5 to 0.7, and the weight of the injection volume reward term can be set to 0.1 to 0.3. The specific design of the multi-objective reward function in the reinforcement learning framework can be flexibly adjusted by those skilled in the art according to the risk tolerance of specific business scenarios. Its underlying logic is a well-known technology in the field and will not be elaborated here.

[0055] The tree search module performs backpropagation calculations of action value and access frequency along the search path. After obtaining the reward value of the leaf nodes, the tree search module backtracks from the bottom up to the root node of the search tree. At each tree node traversed on the backtracking path, the tree search module increments the access frequency of its internal nodes and the historical access frequency of the action branch by one. Based on the frequency update, the tree search module updates the expected action value using the incremental mean method. To ensure the convergence of the value assessment, the tree search module iteratively calculates the expected action value using the incremental mean update formula, specifically configured as follows: ; In the formula, The updated expected value of the action; The expected value of the original actions recorded at the current node; This is the original access frequency; This is the total reward value returned by backpropagation from the leaf node; This refers to the updated access frequency. This smooth update mechanism helps to gradually approach the true long-term value of a node through multiple iterations.

[0056] The tree search module utilizes differentiated logs to perform reverse rollback of physical states. Monte Carlo tree search requires exploring different action branches under the same base state. Early fluid injection and asynchronous stress spillover alter the pressure, stress, and permeability parameters of the underlying two-layer heterogeneous graph. To prevent data contamination, the tree search module reads the differentiated log list within the current node while tracing back to the previous node. To ensure the accuracy of multi-level parameter recovery, the tree search module performs log popping operations layer by layer from leaf nodes to the root node according to the depth of the search tree, and within a single node, follows a last-in-first-out (LIFO) reverse order rule to extract operation records from the differentiated logs one by one. For each operation record, the tree search module parses the target node identifier and change attribute category, and overwrites the corresponding physical parameters back to the initial values ​​in the record with the updated values. For example, if a record shows that the node pore pressure changed from the initial 15 MPa to 18 MPa, the system resets the pore pressure of that node in the underlying view to 15 MPa.

[0057] The tree search module resets the global flag snapshot to achieve effective isolation of physical states. In addition to rolling back continuous physical quantities, the discrete states of fault mechanics steady-state flags also need to be synchronously restored. After rolling back attribute parameters, the tree search module modifies the state of nodes in the global flag snapshot that were set to 1 in the current branch back to 0, indicating the cancellation of Coulomb slip and cascade rupture events occurring in that branch. Through the above log-based incremental rollback and state clearing mechanism, the tree search module restores the underlying data structure to its initial environment before the current decision time step began. This isolation mechanism largely avoids the need for large-scale copying of the entire geological map data in memory for each search branch, thereby reducing computer memory allocation overhead while ensuring the non-interference of multi-field coupled evolution processes between independent search branches. At this point, the tree search module completes a single search iteration, and the system uses this as a basis to launch the next round of action candidate and confidence boundary exploration.

[0058] In this embodiment, after numerous Monte Carlo search iterations, the tree search module constructs a search tree containing multiple potential fluid transport branches in computer memory. To transform the abstract node data of the search tree into evaluation results with physical guidance, the result analysis module in the deduction system initiates the extraction process of leakage path deduction results. The specific implementation steps are as follows: The tree search module determines the termination conditions of the Monte Carlo search iteration process and locks the tree structure. When the number of search iterations reaches a preset maximum threshold, or the computation time exceeds the set time window, the tree search module stops the expansion and backpropagation process of the Monte Carlo tree and locks the topology and node attributes of the current search tree. In actual engineering-level carbon dioxide geological storage simulations, this maximum threshold is usually set between 1,000 and 10,000 iterations depending on the size of the fault grid; the time window is configured between 10 and 60 minutes depending on the computing power of the underlying computing server. After the above thorough exploration, each node in this search tree records its accumulated access frequency and expected action value in multiple rounds of exploration, providing a statistical basis for the subsequent extraction of advantageous paths.

[0059] The results analysis module greedily extracts dominant path nodes based on the converged search tree. At the root node of the search tree, the module traverses all adjacent child nodes and selects the child node with the highest access frequency as the optimal action for the current time step by comparing the historical access frequencies stored within the node. The specific mathematical filtering logic for this extraction process is configured as follows: under the current tree node, find the target node that maximizes the historical access frequency of the action branch. In the actual comparison process, if two or more child nodes have the same maximum access frequency, to avoid the algorithm getting stuck in a selection deadlock or random blindness, the results analysis module introduces the expected action value as a second-level judgment index and selects the node with the highest expected action value. Compared with the traditional maximum action value strategy, the system's priority use of the maximum access frequency strategy helps to filter out noise interference from single, accidental high returns during evaluation and improves the statistical robustness of the selected path. After determining the dominant action node, the results analysis module uses this node as a new starting point and continues to execute the same frequency comparison logic downwards along its child node set until it reaches a leaf node or touches a preset safety boundary. This continuous greedy selection process forms a sequence of tree nodes from the root node to the terminal node. For the greedy traversal extraction algorithm of the node sequence in the tree structure, those skilled in the art can use a conventional depth-first traversal combined with conditional selection. The related operations are well-known techniques in this field and will not be elaborated upon here.

[0060] The results analysis module performs a physical mapping of the abstract tree node sequence to three-dimensional geological space. After acquiring the tree node sequence, the module sequentially reads the fluid leading edge pointer data within each node. This fluid leading edge pointer directly corresponds to the physical node number in the two-layer heterogeneous map in the underlying data structure. By extracting and concatenating these numbers, the module transforms the abstract tree node sequence into a physical node trajectory composed of a series of discrete spatial coordinate points. Physically, this trajectory typically represents a spatial channel in a complex fault network where fluid transport resistance is relatively low and the probability of leakage is high—that is, the system's dominant leakage path.

[0061] The results analysis module extracts multiphysics dynamic evolution data along the dominant leakage path. Simply obtaining the spatial trajectory is often insufficient for a comprehensive assessment of the leakage risk of geological sealing projects; the system also needs to reconstruct the fluid-rock mechanical interaction state along this path. The results analysis module analyzes the differential logs and global flag snapshots within each corresponding tree node along the extracted physical node trajectory. By reading the pore pressure and permeability change records in the differential logs and combining them with the mechanical steady-state flag values ​​of each fault segment node in the global flag snapshot, the system reconstructs the stress distribution evolution and cascade fracture response sequence of the underlying fault medium at the corresponding time step as the fluid advances along the dominant leakage path. These accompanying multiphysics state data, freed from the complex search tree background, are extracted separately and stored in a structured manner, thus providing fundamental data support for subsequent quantitative calculations of physical indicators and risk warning.

[0062] In this embodiment, after extracting the node sequence and physical field evolution data of the dominant leakage path, the results analysis module needs to further quantify the specific timescale and flow rate of the fluid along this path. To assess the potential risk of the fluid breaching a specific geological boundary, the results analysis module performs the following index calculation steps: The results analysis module calculates the hydraulic conductivity between adjacent nodes along the dominant leakage path. To reflect the fluid transport capacity of the fault channel to a greater extent, the module extracts the geological geometric parameters and hydraulic properties of adjacent nodes along the dominant leakage path. These geological geometric parameters include the center-to-center distance between nodes and the equivalent conduction cross-sectional area; the hydraulic properties utilize the permeability updated during the simulation. In the physics of seepage, hydraulic conductivity characterizes the volumetric flow rate of fluid passing through the fault segment under a unit pressure difference. Combining the fluid's dynamic viscosity, the results analysis module calculates the hydraulic conductivity between adjacent nodes using seepage mechanics principles. Based on seepage mechanics theory, the results analysis module constructs the following formula for calculating the hydraulic conductivity between adjacent nodes: ; In the formula, For nodes To node Hydraulic conductivity; This is the harmonic average of the permeability at both nodes; It is the equivalent flow-guiding cross-sectional area of ​​the fluid passing through the cross-section of the two nodes; The dynamic viscosity of the injected fluid; The distance between the centers of the two nodes is denoted as . For the aforementioned fluid dynamic viscosity parameters, those skilled in the art can typically obtain them by referring to tables or using empirical formulas based on the equation of state of the injected fluid (such as carbon dioxide) under local geological temperature and pressure conditions.

[0063] The results analysis module determines the flow distribution weights of nodes along the dominant path. In actual underground porous media, fluid transport typically exhibits a multi-branching divergent pattern. To determine the fluid share flowing along the dominant path, the results analysis module traverses each node on the dominant path, reading all adjacent connected edges of the node through a local hydraulic connectivity view. The results analysis module then uses the above logic to calculate the sum of conductivity from the node to all connected child nodes. Considering that fault media may experience local closure due to stress compaction, if the calculated sum of conductivity approaches zero, subsequent division operations will fall into a dead zone. Therefore, the system configures a very small constant at the bottom layer (e.g., set to 1×10). -15 The lower limit of the denominator is used as the cutoff value. After verification, the flow allocation weight of the branch corresponding to the dominant path is defined as the ratio of the hydraulic conductivity of that dominant branch to the total conductivity. This weight allocation mechanism based on the inverse ratio of conductivity to resistance, to some extent, conforms to the physical law that fluids tend to move along high-permeability, low-resistance channels.

[0064] The results analysis module calculates the node passage time of the fluid front on the dominant path. After obtaining the flow allocation weights, the results analysis module, combined with the set injection and initial total flow, calculates the actual volumetric flow rate allocated to the current target node, which is equal to the product of the total flow rate and the flow allocation weight of that node. In this step, if the allocated actual volumetric flow rate is less than the set flow cutoff threshold, it indicates that the fluid has essentially stopped advancing in the current branch, and the system will assign a maximum time constant to that node to characterize the blockage state, thereby avoiding calculation overflow. When the flow rate is within the normal range, the results analysis module extracts the physical volume of the fault grid represented by the target node and the current porosity, and multiplies the two to obtain the effective pore volume of that node. Based on the principle of mass balance, the time increment required for the fluid front to fill the node and break through to the next node is obtained by the following formula: ; In the formula, For fluid to break through the target node Required time increment; For the target node The physical volume of the fault grid it represents; This represents the current porosity of the node. To combine the total traffic with the assigned weights and actually distribute it to the nodes Volumetric flow rate.

[0065] The results analysis module accumulates the time increments and generates a leakage time-series curve. Along the extracted dominant leakage path, starting from the injection source node, the module progressively accumulates the time increments of each subsequent node. The accumulated absolute time represents the estimated time taken for the fluid front to reach a specific spatial location on the path since injection began. The results analysis module correlates these time data with the spatial coordinate sequences of the nodes, constructing a spatiotemporal evolution model of the fluid front's propagation. This model provides a time window reference for fluid reaching shallow aquifers or the surface for engineering applications. The specific derivation process for the conventional volumetric method or based on Darcy's law to derive fluid transport velocity can be completed using conventional reservoir engineering seepage theory formulas by those skilled in the art; the relevant calculation details are well-known in the field and will not be elaborated here.

[0066] Specific application example: Risk simulation of carbon dioxide geological storage leakage (I) Project Background: An oil and gas field company is conducting a demonstration project for carbon dioxide geological storage in a deep saline aquifer in a basin. The study area contains a complex discrete fault network (containing approximately 500 cutting fault segments). The goal is to simulate whether carbon dioxide fluid will leak along the faults to the shallow freshwater layer (i.e., the preset safety boundary) under the condition of continuous injection of 1 million tons of carbon dioxide.

[0067] (II) Implementation of the simulation process: Graph Initialization: The simulation system maps 500 fault segments to 500 graph nodes. Due to the geometric intersections of the faults in three-dimensional space, the system generates 1200 hydraulically connected edges. Simultaneously, using an elastic half-space analytical model, the Coulomb stress between nodes is calculated, and nodes with stress transmission greater than 0.02 MPa are connected to directed edges, generating 3000 mechanically dispersed edges. In the initial stage, these mechanically dispersed edges are globally masked and remain invisible.

[0068] Network inference: The multi-head graph attention network inside the module reads the local hydraulic connectivity edge view and basic physical parameters such as permeability and initial geostress of each fault node, infers the highest prior probability of fluid migration to the northeast fault (e.g., reaching 85%), and uses the regression head to output a heterogeneous stress correction factor (e.g., an output value of 1.2, to correct the stress concentration error caused by local mudstone interlayers).

[0069] Tree search and multi-field coupling: Fluid front extension: The simulation system begins constructing a Monte Carlo search tree with the node corresponding to the injection well location as the root node. As carbon dioxide fluid continues to be injected, the fluid front advances to node 12, causing the pore pressure at that target node to rise from 20 MPa to 28 MPa.

[0070] Mechanical instability and asynchronous propagation: A significant increase in local pressure leads to a decrease in the effective normal stress at node 12, triggering the Mohr-Coulomb rupture criterion. At this point, the system immediately suspends the fluid leading edge propagation process, dynamically removing the masking effect of the mechanically dispersed edges around node 12. Node 12 ruptures and releases stress, and this stress disturbance affects the distant, unconnected node 35, causing a surge in its shear stress and subsequently leading to cascading rupture.

[0071] Secondary obstruction determination: After the slippage of node 12, due to the grinding and compaction effect of the fault gouge, its local permeability dropped sharply to 0.05 millidarcy, successfully triggering the secondary obstruction determination threshold set by the system. The system cut off the subsequent path of fluid through node 12 at the bottom layer. The algorithm then used differentiated logs to perform state reverse rollback, undo parameter changes, and redirect the branch to explore other feasible paths.

[0072] Path Output: After 5000 Monte Carlo search iterations, the system outputs the most frequently accessed dominant leakage path from the converged search tree (its physical trajectory is mapped as: node 1 to node 8 to node 22 to node 89). Combined with flow allocation weight calculations, it is shown that due to the compaction and sealing effect caused by local fault slippage, the time for carbon dioxide fluid to eventually break through to the shallow freshwater layer is estimated to be 3.5 years, rather than the 0.8 years predicted by traditional static numerical models.

[0073] (III) Experimental verification, effect comparison and illustration: To verify the effectiveness of the technical solution of this invention, researchers conducted a comparative experiment using real three-dimensional geological data of the aforementioned basin and historical microseismic monitoring data. This experiment included three technical solutions for comparison testing, the specific meanings of which are respectively attached in the appendix. Figure 3 and attached Figure 4 The horizontal axis labels (bottom text) are divided into three categories from left to right: Traditional numerical simulation: refers to the basic comparative method of using traditional full-mesh generation and multi-field coupled finite element and finite volume joint simulation software for full-time integral solution; Standard tree search: refers to a comparative method that only uses the basic Monte Carlo search algorithm to find hydraulic connectivity paths, without introducing prior guidance from graph neural networks and without considering the asynchronous ripple effect of stress; The method of this invention refers to the deduction method proposed in this patent, which is jointly driven by graph network prior guidance and dynamic asynchronous expansion of two-layer heterogeneous graphs.

[0074] Combined with appendix Figure 3 Comparison of computational efficiency and conclusions are shown in the appendix. Figure 3The vertical axis on the left side of the chart, labeled "Computation Time (hours)," represents the actual physical time consumed by the algorithm from input geological data to output of the final result. This is further analyzed in conjunction with experimental data. The dark gray bars in the figure and the top value (124.5) represent the computation time of the "traditional numerical simulation" method, which is 124.5 hours. The highest bar indicates that this method is extremely inefficient due to the need for full-grid matrix solving (experiments recorded a peak memory usage of up to 128 gigabytes).

[0075] The light gray bar in the figure and its top value (8.2) represent the computation time of 8.2 hours for the "standard tree search" method. This bar is the lowest, with a peak memory usage of only 16 gigabytes. Although the computation speed is fast, the derivation process is overly simplified due to the complete lack of mechanical coupling calculations.

[0076] The medium-gray bars in the figure and their top value (11.5) represent the computation time of the "method of this invention" as 11.5 hours (with a peak memory usage of only 22 gigabytes). Conclusion (efficiency): The height of the medium-gray bars corresponding to the method of this invention is significantly lower than that of the dark gray bars in traditional simulations. This indicates that by introducing a two-layer graph network mask and a local differential logging mechanism, this invention successfully avoids the cumbersome process of solving the full-mesh stiffness matrix, significantly reducing the computation time by approximately 90% and greatly improving the running efficiency of engineering simulations.

[0077] Combined with appendix Figure 4 The comparison of prediction accuracy and high fidelity conclusions are shown in the appendix. Figure 4 The vertical axis on the left side of the chart, labeled "Path Hit Rate (%)", represents the degree of consistency between the predicted paths and actual microseismic monitoring records. A higher value indicates more accurate predictions. Further analysis is needed, combined with experimental data. The dark gray bars in the figure and the text at the top (94%) represent a path hit rate of 94% for the "traditional numerical simulation" method (88% capture rate for microseismic rupture events). Because its computational model is the most complex and comprehensive, it serves as a high-fidelity benchmark.

[0078] The light gray bars in the figure and the text at the top (58%) represent the path hit rate of the "standard tree search" method, which is only 58% (the microseismic rupture event capture rate is 0%). This extremely low bar indicates that conventional search methods that do not incorporate graph networks and mechanical coupling mechanisms are unable to accurately predict complex fault leakage patterns and are prone to deviating from reality.

[0079] The medium-gray bars in the figure and the text at the top (92%) represent a path hit rate of 92% for the "method of this invention" (while the microseismic rupture event capture rate reaches 85%). Conclusion (high fidelity): The height of the medium-gray bars corresponding to the method of this invention is close to the height of the dark gray baseline bars. This indicates that, while significantly reducing computation time, this invention, through prior guidance from a multi-head graph attention network and the introduction of mechanical remote dispersion edge calculations to compute cascade rupture and secondary sealing effects, achieves a main leakage path hit rate and physical evolution fidelity highly consistent with the full numerical simulation. This invention can still highly reproduce the real physical state of multi-field coupling under complex strata, ensuring high accuracy and reliability of the deduction results.

Claims

1. A method for fault leakage inference based on graph networks and tree search, characterized in that, include: Obtain discrete fault network geological data for the study area, which includes the three-dimensional spatial geometric coordinates of fault segments, initial physical parameters, and preset safety boundary coordinates for risk assessment. The fault segment is mapped as a node, a hydraulic connectivity edge is constructed based on the spatial intersection relationship, and a mechanical remote dispersion edge is constructed based on the static Coulomb stress transfer coefficient and the preset wave sweep threshold. A two-layer heterogeneous graph is generated based on the hydraulic connectivity edge and the mechanical remote dispersion edge, and a masking process is applied to the mechanical remote dispersion edge. The two-layer heterogeneous graph is input into a multi-head graph attention network, and the action prior probability distribution and dynamic sweep coefficient are output. A Monte Carlo search tree is constructed by setting the node corresponding to the injection well location as the root node, and fluid front expansion is performed along the hydraulically connected edge. When a node satisfies the Mohr-Coulomb rupture criterion, the fluid leading edge extension is suspended, the mask shielding of the corresponding node's mechanical remote dispersion edge is removed, the stress wave asynchronous leading edge extension is initiated, and the physical parameters of the far-end connected node are updated using the dynamic sweep efficiency coefficient. After completing a single branch search, the reward value of the search branch is calculated and backpropagation is performed. The changes in node pressure and penetration rate are revoked layer by layer according to the differential log, and the initial physical state of the two-layer heterogeneous graph is restored. When the maximum number of iterations is reached, the Monte Carlo search tree is traversed to extract the sequence of nodes with the highest access frequency, and the sequence of nodes is output as the dominant leakage channel.

2. The fault leakage inference method based on graph networks and tree search according to claim 1, characterized in that, The process of generating a two-layer heterogeneous graph and applying a mask to the mechanically divergent edge includes: Assign each node a physical feature vector containing initial permeability, porosity, initial normal stress, initial shear stress, cohesion, friction coefficient, and initial formation pore pressure; The static Coulomb stress transfer coefficient between nodes is calculated using an elastic analytical model. When the absolute value of the static Coulomb stress transfer coefficient is greater than the preset swirl threshold, a directed edge is established between the nodes as the mechanical remote diffusion edge. The hydraulically connected edge and the mechanically dispersed edge are integrated to generate the two-layer heterogeneous graph, and a visibility flag bit with an initial value of false is assigned to the mechanically dispersed edge to achieve the masking process.

3. The fault leakage inference method based on graph networks and tree search according to claim 1, characterized in that, The prior probability distribution and dynamic sweep coefficient of the output fluid movement between nodes include: The multi-head graph attention network is used to perform aggregation operations on node neighborhood features within the graph attention aggregation layer to extract hidden layer representations that integrate local connected topology and heterogeneous attributes. The aggregated hidden layer representation is mapped to the prior probability distribution of the action of moving from the fluid front to connected neighbor nodes using a classification prediction head; The regression head outputs a non-homogeneous stress correction factor for the node pairs connected by the mechanical remote dispersion edge, and the non-homogeneous stress correction factor is multiplied by the pre-calculated static Coulomb stress transfer coefficient to generate the dynamic sweep coefficient.

4. The fault leakage inference method based on graph networks and tree search according to claim 1, characterized in that, The construction of the Monte Carlo search tree includes: The tree nodes of the Monte Carlo search tree contain a composite state view, which includes a fluid leading edge pointer, a local graph attribute view, a differential log, and a global flag snapshot. The global flag snapshot is a one-dimensional bitmap array used to record the mechanical steady-state flag state of each node; The differentiated log is a tuple containing the target node identifier, the changed attribute category, the original value, and the updated value, used to record the incremental modification operations of the underlying physical parameters caused by the current search step.

5. The fault leakage inference method based on graph networks and tree search according to claim 4, characterized in that, The fluid front extension along the hydraulically connected edge includes: The confidence upper bound value of each node in the action candidate set is calculated by combining the prior probability distribution of the action and using the upper confidence bound algorithm of the predictor. The node with the largest confidence upper bound value is selected as the expansion target node for fluid leading edge expansion. Calculate and update the pore pressure of the extended target node, and calculate the effective normal stress in combination with the initial normal stress of the extended target node; The effective normal stress is substituted into the Mohr-Coulomb fracture criterion for determination. If the determination is successful, the mechanical steady-state flag value of the extended target node is modified to a locked state and the fluid leading edge extension is suspended.

6. The fault leakage inference method based on graph networks and tree search according to claim 5, characterized in that, The initiation stress wave and leading edge asynchronous propagation include: The visibility flag of the mechanical remote divergence edge corresponding to the expanded target node that has broken is changed from false to true, the mask of the mechanical remote divergence edge is dynamically removed, and a local mechanical topology view is constructed. Calculate the static stress drop of the expanded target node that has fractured, and push the expanded target node as a disturbance source node into the asynchronous expansion queue; The induced Coulomb stress change of the affected connected node is calculated by multiplying the dynamic sweep efficiency by the static stress drop. The equivalent shear stress of the affected and connected nodes is updated and the Mohr-Coulomb rupture criterion is called for determination. If the rupture condition is met, the affected and connected nodes are pushed to the tail of the asynchronous expansion queue. This process is repeated until the asynchronous expansion queue is empty, and then the asynchronous expansion process is terminated and the fluid leading edge expansion is resumed.

7. The fault leakage inference method based on graph networks and tree search according to claim 6, characterized in that, The process of initiating the asynchronous propagation of the stress wave and its leading edge also includes performing secondary blocking determination and anti-cascading cutoff control: The updated permeability values ​​of the affected and connected nodes that experienced slip rupture were calculated using the permeability multiplier method. If the updated permeability value is lower than the preset critical closed permeability threshold, it is determined that a secondary blockage has occurred, and the hydraulic connection edge connected to the affected connected node is set to a temporary shielded state. Before the affected and connected nodes are pushed into the asynchronous expansion queue, the propagation depth and the induced Coulomb stress change are checked. If the propagation depth reaches a preset cascade depth threshold, or the absolute value of the induced Coulomb stress change is less than a preset minimum stress cutoff threshold, the operation of enqueuing the affected and connected nodes is forcibly discarded, and the mechanical transmission chain is cut off.

8. The fault leakage inference method based on graph networks and tree search according to claim 4, characterized in that, The calculation of the reward value of the search branch and the execution of backpropagation include: Extract the multiphysics state parameters corresponding to the leaf nodes, and calculate the comprehensive reward value of the current leaf node by combining the total amount of injected fluid, the number of fault nodes that have slipped and ruptured, and the Euclidean distance between the fluid front and the preset safety boundary through normalization and weighted summation. The algorithm backtracks from bottom to top, layer by layer, to the root node of the Monte Carlo search tree, updates the node access frequency and historical access frequency of the action branch of the tree node on the backtracking path, and uses the incremental mean update formula to iteratively calculate the expected value of the action.

9. The fault leakage inference method based on graph networks and tree search according to claim 8, characterized in that, The step of reversing the changes in node pressure and permeability values ​​layer by layer according to the differentiated logs to restore the initial physical state of the two-layer heterogeneous graph includes: Based on the level depth of the Monte Carlo search tree, log popping operations are performed layer by layer from the leaf node to the root node, and operation records in the differentiated logs are extracted within a single node following the last-in-first-out reverse order rule. The physical parameters corresponding to the corresponding node in the operation record are overwritten from the updated values ​​back to the original values; The state of the node that was set to 1 in the current branch in the global flag snapshot is changed back to 0, and the Coulomb slip and cascade breakage events that occurred in the current branch are undone.

10. The fault leakage inference method based on graph networks and tree search according to claim 1, characterized in that, After outputting the node sequence as the dominant leakage channel, the method further includes an index calculation step: Extract the geological geometric parameters and hydraulic properties of adjacent nodes in the dominant leakage channel, and calculate the hydraulic conductivity between the adjacent nodes; The flow distribution weight is determined by the ratio of the hydraulic conductivity of the dominant branch to the sum of the conductivity of all connected sub-nodes. The actual volumetric flow rate allocated to the target node is calculated by combining the initial total flow rate of the injection well with the flow rate allocation weight; The time increment for fluid to break through the target node is calculated based on the actual volumetric flow rate and the effective pore volume of the corresponding node, and the time increment along the dominant leakage channel is accumulated step by step to generate a leakage time series curve.

Images