Method and system for predicting gas drive effect of compact low-permeability reservoir fracturing well pattern

Abstract

The invention provides a compact low-permeability reservoir fracturing well pattern gas drive effect prediction method and system, and relates to the technical field of gas drive effect prediction. The method comprises the following steps: firstly, determining basic geometric parameters of a rhombic inverted nine-point well pattern of a target block and fracturing design parameters to be optimized; and constructing a well pattern embedded into the main fracture and a fracture parameterized model, and reflecting a coupling relationship between the fracture and matrix flow. A training sample is generated through experimental design, rapid flow simulation is conducted through a parameterized model, and space collaboration quantitative indexes such as a gas shortest breakthrough time index, an effective wave and coefficient correction factor and a displacement pressure gradient strengthening factor are calculated. Performing extended production dynamic simulation on part of samples to obtain key indexes such as accumulative oil production and accumulative gas-oil ratio, and training a machine learning agent model; and with maximization of the cumulative oil production and the gas shortest breakthrough time index and control of the cumulative gas-oil ratio as optimization objectives, a multi-objective optimization algorithm is adopted to obtain a Pareto optimal solution set, and a final recommendation scheme is screened out.

Description

Technical Field

[0001] This invention relates to the field of gas drive effect prediction technology, specifically to a method and system for predicting the gas drive effect of fracturing well networks in tight, low-permeability oil reservoirs. Background Technology

[0002] Tight, low-permeability reservoirs, as an important replacement area for global oil and gas resources, rely heavily on the combination of hydraulic fracturing and enhanced oil recovery (EOR) technologies for their effective development. Among these, the development approach that combines fracturing to create complex fracture networks with gas injection displacement has become a key technological direction for improving single-well production and overall recovery in these reservoirs. Particularly in blocks using regular well patterns such as the diamond-shaped inverted nine-point well network for overall development, the synergy between fracturing design and injection-production well network layout directly determines the gas sweep efficiency, the stability of the displacement front, and ultimately, the economic benefits. In this specific application scenario, accurately predicting the gas drive development effect of different fracturing schemes (such as fracture length, azimuth, and conductivity) within a given well network before construction, and optimizing the design accordingly, is a core engineering challenge for oilfield decision-makers.

[0003] Currently, the industry mainly relies on two methods to predict the gas drive effect of fracturing well networks: one is the detailed simulation of the entire oil and gas reservoir based on numerical simulation software, and the other is the empirical analogy method based on field statistics. Although the former has a clear mechanism, the modeling process is complex and computationally time-consuming, making it difficult to use for rapid iteration and optimization design of a large number of parameter combinations. Moreover, its accuracy is heavily dependent on geological and fracture parameters that are difficult to obtain accurately. The latter relies too much on historical data of specific blocks, lacks universality, cannot guide the optimization design of new areas, and cannot quantify the dynamic interaction between the fracture system and the spatial layout of the well network.

[0004] Of particular concern is the lack of quantitative evaluation and optimization methods for the crucial factor of spatial synergy between fracturing networks and injection-production well networks in existing technologies. In a diamond-shaped inverted nine-point well network, improper design of the length and orientation of fracturing fractures in injection and production wells can easily lead to short-circuiting of injected gas along high-conductivity fractures, severely reducing macroscopic sweep efficiency and oil displacement effects. Traditional methods often address this problem in a fragmented manner: first, the well network is deployed, then individual well fracturing is designed independently, and finally, the effect is verified through numerical simulation. This sequential design process cannot proactively optimize fracture parameters to match the well network in the early stages, nor can it quickly assess systemic issues such as how the relative orientation of fractures between injection and side wells affects the effectiveness of corner wells. As a result, many schemes only reveal severe gas channeling and unsatisfactory results after investment and implementation.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a method and apparatus for predicting the gas drive effect of fracturing well networks in tight, low-permeability oil reservoirs, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs, comprising the following steps: Step 1: Determine the basic geometric parameters of the diamond-shaped inverted nine-point well pattern in the target block, including the well spacing and row spacing of injection and production wells; set the fracturing design parameters to be optimized, including the main fracture half-length of each well, the angle between the main fracture azimuth and the well row direction, and the conductivity of the main fracture. Step 2: Simplify the main fracture of each well in the well network into a linear high-conductivity channel, and embed the main fracture into the matrix grid model that characterizes the matrix properties to form a well network and fracture parameterization model that reflects the flow coupling relationship between the fracture and the matrix. Step 3: Set the range of values ​​for fracturing design parameters, and generate N sets of training samples within this range using experimental design methods; for each set of training samples, perform rapid flow simulation using a parameterized model of well network and fracture, and calculate the spatial synergy quantification index of the corresponding sample based on the simulation results; the spatial synergy quantification index includes the gas shortest breakthrough time index, the effective sweep efficiency correction factor, and the displacement pressure gradient enhancement factor. Step 4: For each training sample, a numerical simulator is used to perform extended production dynamic simulation to calculate the key indicators of the gas drive development effect of each training sample within the complete development cycle, including cumulative oil production and cumulative gas-oil ratio; using basic geometric parameters and fracturing design parameters to be optimized as input features, and spatial synergy quantification indicators and key indicators as output labels, a machine learning proxy model is trained. Step 5: With the optimization objectives of maximizing cumulative oil production and the shortest gas breakthrough time index, and with the cumulative gas-oil ratio controlled below the target threshold as a constraint, a multi-objective optimization algorithm is used to obtain the Pareto optimal solution set, and then the final recommended scheme is selected by combining key indicators.

[0008] Furthermore, the method for determining the value range of the fracturing design parameters to be optimized is as follows: For the range of the main fracture half-length, based on meeting the minimum engineering requirements for forming an effective modified zone during fracturing construction, its upper limit is set according to the well spacing and row spacing of the injection and production wells in the aforementioned diamond-shaped inverted nine-point well network, and its upper limit satisfies the constraint conditions. ,in, This indicates the upper limit of the half-length of the main crack; and This is the proportionality coefficient, with a value range of [0, 0.5]. The distance between injection and production wells, The spacing between injection and production wells; this constraint ensures that the main fractures do not physically intersect within the well network unit; The angle between the azimuth of the main fracture and the direction of the well network is called the fracture azimuth angle. Its value range is set based on the well network layout and the principle of minimizing gas breakthrough risk. ,in, The angle between the normal direction of the line connecting the gas injection well and the corner well and the well row direction. Allowable design deviation angle; The conductivity range of the main fracture is determined dimensionlessly based on the matrix permeability of the target reservoir, and its relationship is expressed as follows: ,in, This represents the dimensionless conductivity of the crack. This indicates the conductivity of the main fracture. The matrix permeability of the target reservoir is represented; the conductivity of the main fracture should be within the range of [1, 100]. This indicates the half-length corresponding to the main crack; Based on the function and spatial relationship of different well locations in the diamond-shaped inverted nine-point well network, the fracturing design parameters are differentiated into groups. Specifically, the gas injection well located at the center of the well network is divided into a gas injection well group, and its fracturing design parameter group is set, denoted as [group name missing]. ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the gas injection well; the four corner production wells located at the four corners of the well network are classified as corner production well groups, and their fracturing design parameter groups are set, denoted as... ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the corner well; the four side wells located at the midpoint of the four edges of the well network are divided into a side well production well group, and their fracturing design parameter group is set, denoted as . ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the side well, respectively. The fracturing design parameters to be optimized include the fracturing design parameter groups for the aforementioned gas injection well groups, corner well production well groups, and side well production well groups.

[0009] Furthermore, a minimum symmetric element of the rhombic inverse nine-point well network is selected as the physical modeling domain to construct the matrix mesh model; Within the physical modeling domain, a background Cartesian mesh is initialized, and a local coordinate system is established with the location of the gas injection well as the origin. Based on the main fracture design parameters of each well, the mesh is non-uniformly refined. The specific rule is as follows: for any fracture, along the well row direction and its perpendicular direction, the mesh is refined at the midpoint of the line connecting the injection and production wells and in the expected extension area of ​​the fracture. The mesh size satisfies the following constraints: in, The minimum side length of the grid on the x-axis. This represents the minimum side length of the grid along the y-axis. These are the main fracture half-lengths of the gas injection well, corner well, and side well, respectively. is the first encryption coefficient, with a value range of [10, 20]; In regions far from wells and fractures, the mesh size increases toward the model boundary in a geometric progression, with the growth rate ranging from (1,2]. Symmetric boundary conditions are applied to the boundaries of the physical modeling domain, specifically: the boundary parallel to the wellbore direction is set as a flow-free boundary, and the boundary perpendicular to the wellbore direction is set as a constant pressure gradient boundary or a connectivity boundary, in order to simulate the flow state of the symmetric unit in an infinitely large periodic well network. Based on the concept of an embedded discrete fracture model, the main fracture of each well in the well network is simplified into a linear high-conductivity channel. The embedding method and flow calculation rules of the linear high-conductivity channel in the matrix grid model are as follows: Each main crack is discretized into a series of interconnected crack units based on its length and orientation angle. Calculate the geometric intersection relationship between each crack element and the matrix mesh it passes through, and establish non-adjacent connections between the crack element and the matrix mesh it passes through based on the principle of embedded discrete crack model; The formula for calculating the conductivity of the non-adjacent connections is as follows: in, Indicates the conductivity of non-adjacent connections; The effective permeability of the matrix mesh in the direction normal to the fracture surface; This represents the area of ​​intersection between the crack element and the matrix mesh it traverses; This represents the vertical distance from the center of the matrix grid to the crack cell it passes through.

[0010] Furthermore, obtaining the training samples specifically includes: A complete set of fracturing design parameters to be optimized is defined as a 9-dimensional vector; the sampling range of the parameters in each dimension is set according to the constraints. In a 9-dimensional parameter space, generate a Latin hypercube design matrix containing N sample points, where N is not less than 10 times the parameter dimension, i.e., N≥90; N represents the number of training samples. Each row of the Latin hypercube design matrix is ​​combined with the basic geometric parameters to form a training sample; by traversing each row of the Latin hypercube design matrix, N training samples are obtained. Performing the aforementioned rapid flow simulation specifically refers to running a simplified, non-implicit flow calculation process based on a well pattern and fracture parameterization model, and satisfying at least one of the following conditions: Steady-state flow simulation is used, which solves for the pressure field distribution under a constant injection-production pressure difference, without simulating the change of saturation over time; this simulation is used to calculate the conductivity field, streamlines, and pressure gradient. Streamline simulation is used to calculate the streamline distribution based on the pressure field obtained from steady-state simulation, and one-dimensional tracer or leading-edge propulsion equations are solved along the streamlines to assess the breakthrough time and the affected area. The simulation timeframe is limited, focusing only on the early gas emergence stage, and does not include long-term simulations of the entire life cycle.

[0011] Furthermore, the calculation process for the shortest breakthrough time exponent of the gas is as follows: The set of effective flow paths connecting the injection wells and each production well is identified using the shortest path search algorithm in graph theory; each path consists of several matrix grid blocks or fracture elements, or is composed of matrix grids and fracture elements connected in series. For each path, calculate the quasi-steady-state flow time of gas breakthrough along that path: in, Indicates the quasi-steady-state flow time; For matrix porosity, The gas saturation is the average value taken during gas frontal displacement. The value is 0 for the pore volume of the matrix mesh to which the i-th path belongs; if the starting point is a crack element, its value is 0. This represents the conductivity between two connected computational units corresponding to the i-th path; when the connection is between a crack unit and a matrix mesh, this conductivity value is the conductivity of a non-adjacent connection; when the connection is between two adjacent matrix meshes, this conductivity value is the interface conductivity between them. The pressure difference between the two ends of the i-th path is calculated through a single steady-state flow under the initial injection-production pressure difference condition; i is the path index, and N represents the total number of paths; the path refers to the path through which gas flows from the center of one computational cell to the center of the next directly connected computational cell, and the computational cell refers to the matrix grid or fracture cell; For each production well, the minimum value of all its corresponding paths is denoted as the shortest gas breakthrough time index for that production well. The calculation steps for the effective sweep efficiency correction factor specifically include: In the parameterized model of well network and fracture, a tracer is continuously injected from the injection wells, while all production wells produce at a constant flow rate or constant pressure. The simulation stops when the tracer concentration in the produced fluid of any production well first reaches the preset threshold. This moment is recorded as the gas breakthrough time. At the moment of gas exposure, the total volume of all matrix grids in the statistical model with tracer mole fractions greater than zero is denoted as the effective swept volume. The theoretical swept volume of a symmetrical element in a rhombic inverse nine-point well network under ideal piston displacement is calculated using the following formula: in, Where H is the theoretical swept volume and H is the effective reservoir thickness; in, For effective sweep efficiency correction factor, For effective sweep volume; The specific steps for calculating the displacement pressure gradient enhancement factor include: A steady-state flow simulation was performed using a parameterized model of well grid and fracture to obtain the pressure values ​​of each matrix grid and each fracture element. For each matrix mesh, identify all crack elements that are not adjacent to it; calculate the vertical distance between these crack elements and the matrix mesh one by one, and regard the crack element corresponding to the minimum vertical distance as the associated crack element of the matrix mesh; The formula for calculating the effective displacement pressure gradient of the matrix mesh is: in, This represents the effective displacement pressure gradient of the matrix grid; The pressure value associated with the crack element; This represents the pressure value of the matrix mesh; This represents the vertical distance between the matrix mesh and the associated crack element, with the gradient direction perpendicular to the crack wall. All matrix meshes with a vertical distance less than a preset distance threshold are selected to form the crack's direct influence zone; the distance threshold range is [value missing]. ; The arithmetic mean of the effective displacement pressure gradients of all matrix meshes in the crack-affected zone is the average effective displacement pressure gradient. When there are no fractures, the average pressure gradient between the injection well and the farthest corner well is taken as the reference pressure gradient, and the calculation formula is as follows: in, The reference pressure gradient; This refers to the bottom pressure of the gas injection well. The bottom-hole flowing pressure of the corner well furthest from the injection well; This is the straight-line distance between the bottom of the gas injection well and the farthest angle well. The displacement pressure gradient enhancement factor is the ratio of the average effective displacement pressure gradient to the reference pressure gradient.

[0012] Furthermore, the preliminary screening of N training samples based on spatial synergy quantification indicators specifically includes: The screening criteria include that the gas shortest breakthrough time index is greater than the preset critical time, and the effective sweep efficiency correction factor is greater than the preset critical sweep efficiency. From the training samples that meet the screening criteria, Q groups of training samples are randomly selected according to a preset ratio as a subset of high-fidelity simulation samples for extended production dynamic simulation. For each training sample in the high-fidelity simulation sample subset, based on its corresponding fracturing design parameters and basic geometric parameters, the following operations are performed to construct a full-well-field numerical simulation model: Using a complete well group of a rhombic inverse nine-point well network as the simulation region, a three-dimensional structured mesh was constructed using a local mesh refinement method. Refinement was applied around the expected path of the main fracture to ensure that the mesh dimensions of the fracture-penetrating mesh in the fracture extension direction met the following requirements: ,in, For the extended dimension length; This is the second encryption factor; Using the crack parameters of the training samples, the main crack is characterized in the refined mesh by the equivalent conductivity assignment method. Specifically, the permeability of each matrix mesh through which the crack passes is modified to the ratio of the conductivity of the main crack to the mechanical width of the crack. A multi-component fluid model including methane, ethane, propane, and C7+ heavy fractions was adopted, and the Peng-Robinson equation of state was used to describe the phase changes of oil and gas during the gas injection process; initial reservoir pressure, temperature, and fluid saturation were set. Set the gas injection rate or bottom hole pressure of the gas injection well, set the production rate or bottom hole pressure of the production well, and define the total simulation time as the complete development cycle. Run the constructed full-well-field numerical simulation model and record the daily oil production and daily gas production of all production wells at each simulation time step; after the simulation is completed, calculate the cumulative oil production and cumulative gas-oil ratio of the training sample by time integration; For the remaining training samples that were not selected into the simulated sample subset, their key metrics were temporarily filled using the preliminary predictions of the built machine learning surrogate model or interpolation of similar samples for the initial training of the surrogate model.

[0013] Furthermore, the construction and iterative verification process of the machine learning agent model specifically includes: Using basic geometric parameters and fracturing design parameters to be optimized as input features, the calculated spatial synergy quantitative index and key index are used as output labels, and the machine learning proxy model is initially trained based on N sets of training sample data. Define a prediction uncertainty threshold; use the machine learning proxy model after initial training to predict regions in the training sample space that have not undergone high-fidelity simulation and estimate their prediction uncertainty; select several new sample points whose prediction uncertainty is higher than the prediction uncertainty threshold, calculate the spatial coherence quantification index of these new sample points, and add them to the training dataset to obtain an enhanced training dataset. Using the enhanced training dataset, the active learning process of the machine learning agent model is retrained until the prediction uncertainty of new sample points is lower than the prediction uncertainty threshold, or the total number of high-fidelity simulations reaches the preset upper limit. At this point, the machine learning agent model is considered to have been trained.

[0014] Furthermore, the selection of the final recommended solution specifically includes: The optimization variables are defined as the fracturing design parameters of the target block; a multi-objective optimization problem is constructed, which includes two maximization objective functions, namely maximizing the cumulative oil production and maximizing the shortest gas breakthrough time exponent, and the constraint that the cumulative gas-oil ratio is not greater than the target gas-oil ratio. Using the trained machine learning agent model as the fitness evaluation function, a non-dominated sorting genetic algorithm with an elitist strategy is used to iteratively solve the problem within the range of values ​​of the optimization variables to obtain a Pareto optimal solution set that satisfies the constraints. Each solution corresponds to a set of fracturing design parameters and their predicted key indicators and spatial synergy quantitative indicators. Set the effectiveness filtering conditions: the effective sweep coefficient correction factor is not less than the preset correction factor threshold, and the displacement pressure gradient enhancement factor is not less than the preset enhancement factor threshold. From the Pareto optimal solution set, remove all individuals that do not meet the above validity filtering conditions, and each individual corresponds to a solution; the set of the remaining individuals after removal is defined as the engineering feasible solution set; Clustering algorithms are used to cluster the fracture parameter patterns of the feasible solutions in the engineering project set and group them into several typical scheme categories. Within each scheme category, the scheme with the largest cumulative oil production that satisfies the constraint that the cumulative gas-oil ratio is not greater than the target gas-oil ratio is selected as the representative scheme of the category. If no solution satisfies this condition, the scheme with the smallest cumulative gas-oil ratio is selected as the representative scheme of the category. For the representative schemes of each final category, calculate their unit gas recovery efficiency; that is, the ratio of the cumulative oil recovery to the cumulative gas injection volume corresponding to the cumulative oil recovery is the unit gas recovery efficiency; output the representative scheme of the category corresponding to the maximum unit gas recovery efficiency as the final recommended scheme, and output its corresponding cumulative oil recovery, cumulative gas-oil ratio and unit gas recovery efficiency together to quantitatively characterize the gas drive development effect of the scheme.

[0015] This invention also provides a system for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs. This system is used to implement the aforementioned method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs, and includes: The basic parameters and constraint definition module is used to determine the basic geometric parameters of the diamond-shaped inverted nine-point well network of the target block, including the well spacing and row spacing of injection and production wells; and to set the range of fracturing design parameters to be optimized, including the half length of the main fracture of each well, the angle between the azimuth of the main fracture and the well row direction, and the conductivity of the main fracture. The well network and fracture model construction module is used to simplify the main fracture of each well in the well network into a linear high-conductivity channel and embed the main fracture into the matrix grid model that characterizes the matrix properties to form a well network and fracture parameterized model that reflects the flow coupling relationship between the fracture and the matrix. The spatial synergy assessment module is used to set the range of values ​​for fracturing design parameters and generate N sets of training samples through experimental design methods. For each set of training samples, rapid flow simulation is performed using a parameterized model of well network and fracture, and the spatial synergy quantification index of the corresponding sample is calculated based on the simulation results. The spatial synergy quantification index includes the gas shortest breakthrough time index, the effective sweep efficiency correction factor, and the displacement pressure gradient enhancement factor. The surrogate model training module is used to perform extended production dynamic simulation using a numerical simulator for each training sample, and calculate the key indicators of gas drive development effect during the complete development cycle of each training sample, including cumulative oil production and cumulative gas-oil ratio; using basic geometric parameters and fracturing design parameters to be optimized as input features, and spatial synergy quantification indicators and key indicators as output labels, to train the machine learning surrogate model. The multi-objective optimization and decision-making module is used to maximize the cumulative oil production and the shortest gas breakthrough time index as optimization objectives, and control the cumulative gas-oil ratio below the target threshold as a constraint. It uses a multi-objective optimization algorithm to obtain the Pareto optimal solution set, and then combines key indicators to select the final recommended solution.

[0016] The technical effects and advantages provided by the present invention in the above technical solution are as follows: This invention addresses the specific scenario of gas drive development following fracturing using a rhomboid inverted nine-point well pattern in tight, low-permeability oil reservoirs, achieving significant technological advancements and application benefits. By constructing a parameterized well pattern and fracture system that integrates an embedded discrete fracture model, and innovatively defining a set of spatial synergy quantification indicators, it enables rapid and quantitative evaluation of the matching relationship between the fracture system and the well pattern geometry. This directly solves the core deficiency in the prior art where the synergy between the well pattern and fracture is difficult to quantify and assess, allowing engineers to proactively identify and mitigate the risk of gas short-circuiting during the design phase, thus optimizing sweep efficiency.

[0017] Based on the aforementioned synergistic indicators and efficient proxy models, this invention constructs an intelligent optimization and decision-making system for fracturing parameters. It can rapidly complete traditional numerical simulations and massive scheme comparison and optimization, outputting a Pareto optimal solution set and a final recommended scheme that balances high production, stable production, and a low gas-oil ratio. This fundamentally changes the traditional sequential, trial-and-error design process, realizing integrated synergistic optimization of fracturing design and well network development. This significantly reduces the cost of numerous ineffective simulations while substantially improving the success rate and economic efficiency of gas drive schemes, providing a reliable technical decision-making tool for the efficient development of tight, low-permeability reservoirs. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a schematic diagram comparing the cumulative oil recovery under different fracturing modes of the present invention; Figure 3 This is a schematic diagram comparing the shortest gas breakthrough time index under different fracturing modes of the present invention; Figure 4 This is a schematic diagram comparing the effective sweep efficiency correction factors under different fracturing modes of the present invention; Figure 5 This is a schematic diagram comparing the cumulative gas-oil ratio under different fracturing modes of the present invention; Figure 6 This is a schematic diagram of the system structure of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0020] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0021] Example: Please see Figures 1 to 5 The present invention provides a technical solution: A method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs, comprising the following steps: Step 1: Determine the basic geometric parameters of the diamond-shaped inverted nine-point well pattern in the target block, including the well spacing and row spacing of injection and production wells; set the fracturing design parameters to be optimized, including the main fracture half-length of each well, the angle between the main fracture azimuth and the well row direction, and the conductivity of the main fracture.

[0022] In this embodiment, the method for determining the value range of the fracturing design parameters to be optimized is as follows: The lower limit of the range for the main fracture half-length must meet the minimum engineering requirements for forming an effective stimulation zone during fracturing. This is usually determined comprehensively based on the rock mechanical properties of the target reservoir, the fracturing fluid system, and the scale of the operation. For example, it can be determined through minimum net present value (NPV) analysis or by comparing successful fracturing cases of similar reservoirs. One specific method is to simulate different small-scale fracture half-lengths using fracturing simulation software, using the minimum economically feasible half-length that can establish an effective seepage channel from the wellbore to the reservoir as the lower limit. The upper limit is set based on the well spacing and row spacing of the injection and production wells in the aforementioned diamond-shaped inverted nine-point well network. Due to the geometric constraints of the well network, it is necessary to prevent physical intersection of fractures within the well network unit to avoid inter-well interference and gas channeling. The upper limit must meet the constraint conditions. ,in, This indicates the upper limit of the half-length of the main crack; and This is the proportionality coefficient, with a value range of [0, 0.5]. The well spacing between injection and production wells refers to the straight-line distance between the gas injection well located in the center of the diamond-shaped inverted nine-point well network and any production well located at one of the corners. The well spacing refers to the vertical distance between two adjacent rows of production wells. When the target block contains multiple geological units, parameters such as well spacing and row spacing can be normalized to fall within the [0,1] interval to eliminate the influence of dimensions and facilitate subsequent analysis.

[0023] proportionality coefficient and The method for determining the value is as follows: The first method is to determine the coefficient through geostatistics: analyze the fracture monitoring data (such as microseismic monitoring) of the target block and similar blocks, statistically analyze the ratio distribution of the actual fracture extension length to the well spacing / row spacing, and take its 95th percentile as the upper limit of the coefficient to ensure safety in most cases.

[0024] The second method is determined through numerical simulation experiments: establishing different... and A conceptual model is used to simulate the gas breakthrough time under various values. A critical value that significantly shortens the gas breakthrough time is selected. The value is then multiplied by a safety factor less than 1 (such as 0.8) to obtain the final coefficient. This multiplication by the safety factor is to further tighten the upper limit, ensuring that in actual design, the crack half-length does not approach the critical length that leads to gas channeling, thereby reducing the risk of gas channeling.

[0025] The third method is to determine the method through engineering experience: when detailed data is lacking, industry-standard practices can be followed, typically taking [a certain value]. For blocks with strong reservoir heterogeneity and high uncertainty in fracture propagation, a more conservative value (such as 0.3) should be taken; for blocks with good homogeneity and stable geostress field, a higher value (such as 0.4) can be taken.

[0026] The range of values ​​[0, 0.5] ensures that the fracture half-length does not exceed half the well spacing or row spacing. Geometrically, this is a sufficient condition to avoid two symmetrically arranged fractures directly connecting between wells; values ​​exceeding 0.5 would lead to a sharp increase in the risk of fractures intersecting between wells.

[0027] The angle between the azimuth of the main fracture and the direction of the well network is called the fracture azimuth angle. Its value range is set based on the well network layout and the principle of minimizing gas breakthrough risk. ,in, The angle between the normal direction of the line connecting the gas injection well and the corner well and the well row direction is calculated as follows: Establish a coordinate system: take the wellbore direction as the reference direction (for example, define it as the positive x-axis direction, corresponding to 0°); Determine the direction vector of the line connecting the centers of injection and production wells: for example, the direction from the central gas injection well to a corner well; Calculate the normal direction of the line connecting them. For a direction angle of... The straight line whose normal direction is or (Usually, the smaller angle between the wellbore direction and the wellbore direction is chosen); Calculate the angle between the normal direction and the wellbore direction (0°), which is... .

[0028] For example, assuming the wellbore direction is east-west (0°), and the line connecting the central injection well to a corner well is 45° east of north, then one of its normal directions is 45° west of north (i.e., 135°). The angle between 135° and 0° is 135°, but the other normal direction is 45° east of south (i.e., -45° or 315°), and its angle with 0° is 45°. Usually, the direction with an absolute value less than 90° is used, therefore... .

[0029] The allowable design deviation angle is determined as follows: The first method is based on the uncertainty analysis of the geostress direction: This is based on the standard deviation of geostress measurement data (such as wellbore collapse, DITF testing). ,set up ,in Typically, 2 to 3 is chosen to represent a 95% to 99% confidence interval, where n is selected based on empirical values.

[0030] The second method is based on fracturing construction error assessment: considering the accuracy of current directional fracturing technology, for example, if the tool face angle control accuracy is ±5°, then it can be taken as... .

[0031] In the absence of specific data, values ​​can be taken based on engineering practice or experience. As a reasonable default design deviation, set This is to accommodate uncertainties in geological understanding and engineering implementation. If the scope is too small, the optimized solution may not be feasible in practice; if the scope is too large, the optimization search space will be too large, reducing efficiency. Generally, a deviation of 10° can better balance the ideal design and realistic constraints.

[0032] The conductivity range of the main fracture is determined dimensionlessly based on the matrix permeability of the target reservoir, and the relationship is expressed as follows: in, This represents the dimensionless conductivity of the crack. This indicates the conductivity of the main fracture. This represents the matrix permeability of the target reservoir; the matrix permeability of the target reservoir is usually derived from core analysis or well logging interpretation data. Preprocessing operations include: depth correction and realignment to ensure accurate matching of permeability data with the corresponding reservoir depth; outlier removal to eliminate abnormally high or low values ​​caused by experimental errors or special minerals; and inter-interval averaging to calculate the geometric mean or flow unit weighted average of the permeability data for the target inter-interval to obtain a value representative of macroscopic flow capacity. The geometric mean better reflects the log-normal distribution characteristics of permeability.

[0033] Outliers can be identified using statistical methods. Specifically, for data that approximately follows a normal distribution, the arithmetic mean and standard deviation of all permeability data are calculated, and data points falling outside the range of (mean - 3 * standard deviation, mean + 3 * standard deviation) are considered outliers. For data that does not rely on strict distribution assumptions, a box plot method is used. Specifically, the upper quartile (Q3, 75th percentile), lower quartile (Q1, 25th percentile), and interquartile range (IQR = Q3 - Q1) of all permeability data are calculated; outlier boundaries are set: the lower bound is... The upper boundary is Data points below the lower bound or above the upper bound are identified as outliers. For permeability data, which typically follows a log-normal distribution, it is recommended to first take the natural logarithm and then apply the box plot method for better results. Similarly, methods based on geological knowledge can be used, combining well logging curves, lithological descriptions, and sedimentary facies analysis. For example, extremely high permeability values ​​in pure sandstone sections may be caused by fractures or pores and are not within the scope of matrix permeability; they should be removed as outliers. Likewise, extremely low values ​​measured in tight mudstone interlayers should also be treated differently. After removing one layer, the geometric mean of the remaining effective permeability data can be directly used to fill the gap.

[0034] This represents the half-length corresponding to the main fracture. When calculating the conductivity range of the main fracture, this value can be estimated by taking the median or upper limit of the range of half-length values ​​of the corresponding well group fracture. It is preferable to take the upper limit to ensure the completeness of the search space.

[0035] The range of dimensionless fracture conductivity is determined based on a comprehensive analysis of seepage mechanics principles, economic benefits, and engineering practice experience. Seepage mechanics studies indicate that when… At that time, the crack's flow capacity is insufficient, and its effect on improving fluid flow is limited, making it difficult to form an effective high-speed seepage channel; when At this point, the fracture conductivity already far exceeds the matrix's fluid supply capacity, and further increasing conductivity has a drastically reduced marginal contribution to production capacity. Further combining engineering economic analysis, a model is established to relate fracture conductivity to fracturing costs (proppant dosage, pumping procedure), and coupled with production capacity forecasting, a net present value analysis is performed to determine the economic optimum that maximizes net present value. It typically falls within a narrow range of 10 to 50. This range represents the option with the highest return on investment under specific oil price and cost conditions.

[0036] To ensure that this optimization method can fully explore various feasible designs, including the economically optimal solution, while avoiding searches in physically inefficient or extremely uneconomical regions, a reasonable and complete parameter space needs to be defined. Therefore, The optimized search range is set to [1, 100]. The lower limit of 1 ensures that the considered fractures have basic conductivity enhancement effects, excluding obviously ineffective designs; the upper limit of 100 covers the economically optimal range and is moderately extended to areas with very low marginal benefits, in case of situations requiring extremely high conductivity under certain special geological conditions or development strategies, thus ensuring the completeness of the search space.

[0037] This invention recognizes the functional differences between wells at different locations in a diamond-shaped inverted nine-point well network and accordingly sets differentiated fracturing design parameters for different groups. This is the foundation for achieving coordinated optimization of the well network and fracture space.

[0038] Based on the function and spatial relationship of different well locations in the diamond-shaped inverted nine-point well network, the fracturing design parameters are set in differentiated groups, specifically: The gas injection wells located at the center of the well network are divided into gas injection well groups, and their fracturing design parameter groups are set, denoted as . ,in, These represent the main fracture half-length, fracture azimuth, and conductivity of the gas injection well. As the source of displacement energy, the primary function of the fractures in a gas injection well group is to efficiently inject gas. A relatively long fracture half-length and high conductivity may be required to increase the injection contact area and reduce the injection pressure.

[0039] The four corner production wells located at the four corners of the well network are classified as corner production well groups, and their fracturing design parameter groups are set, denoted as . ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the corner well. The corner well production group is furthest from the injection well, representing the last point of gas breakthrough and crucial for controlling the overall swept volume. Its fracture design likely focuses more on the connectivity and directionality with the injection well fractures to form an effective displacement channel.

[0040] The four production wells located at the midpoint of the four edges of the well network are divided into a production well group, and their fracturing design parameter group is set, denoted as . ,in, These are the main fracture half-length, fracture azimuth, and conductivity of the side well. Side well production well groups are located between injection wells and corner wells, with a higher risk of gas breakthrough. Their fracture design may require a moderate half-length and conductivity, and may intentionally deviate their fracture direction from the direct line of connection with the injection well (adjustment). This is to delay the gas breakthrough while ensuring a certain level of production capacity.

[0041] The fracturing design parameters to be optimized include the fracturing design parameter sets for the aforementioned gas injection well group, corner production well group, and side production well group. This set comprises nine variables, denoted as... .

[0042] After setting the theoretical scope as described above, final adjustments are needed based on actual engineering constraints: Standardization adjustment: Adjust the calculated parameter range to conform to commonly used standard intervals or specifications in engineering. For example, round the crack half-length to a multiple of 5 meters, and round the conductivity to a multiple of 5 meters. Multiples of.

[0043] With the help of geological and fracturing engineering experts, and combined with their specific understanding of the block (such as fault distribution and lithological boundaries), the set parameter range was reviewed and fine-tuned.

[0044] Step 2: Simplify the main fracture of each well in the well network into a linear high-conductivity channel, and embed the main fracture into the matrix grid model that characterizes the matrix properties to form a well network and fracture parameterization model that reflects the flow coupling relationship between the fracture and the matrix.

[0045] In this embodiment, a minimum symmetric element of the rhombic inverse nine-point well network is selected as the physical modeling domain to construct the matrix mesh model. The method for determining this symmetric element is as follows: The rhomboid inverse nine-point well network exhibits high symmetry. Its smallest symmetric unit is typically a rectangular or rhomboid region containing one complete injection well, one-quarter of the corner wells, and one-half of the side wells. More specifically, it is the region enclosed by the lines connecting the midpoints of two adjacent side wells, the central injection well, and the midpoints of the two adjacent corner wells. Utilizing this symmetry, the flow within this unit can represent the flow characteristics of the entire infinitely large periodic well network. By applying appropriate boundary conditions to the boundaries of this unit, the interactions of the entire well network can be simulated.

[0046] For example, suppose the corner and side wells of a rhombus-shaped inverted nine-point well network are arranged on a rhombus with the major axis being the well spacing *d* and the minor axis being the row spacing *L*. Establish a local coordinate system with the central injection well as the origin (the x-axis is parallel to the well row direction). The smallest symmetric unit can then be a quadrilateral region enclosed by the following vertices: (0,0) (injection well), (d / 2,0) (side well position), (d / 4,L / 4) (a point inside the rhombus), and (0,L / 2) (midpoint of the line connecting the corner well and the origin). In practice, to simplify mesh generation, it is often approximated as a rectangular region with (0,0) and (d / 2,L / 2) as opposite vertices. This rectangular region contains complete flow information.

[0047] Within the physical modeling domain, a background Cartesian mesh is initialized, with the injection well location as the origin, the x-axis parallel to the wellbore direction, and the y-axis perpendicular to the wellbore direction, establishing a local coordinate system. The mesh boundary is aligned with the boundary of the physical modeling domain.

[0048] The selection of the basic grid size should take into account both computational efficiency and resolution requirements. The preferred basic grid size should not exceed 1 / 20 to 1 / 10 of the minimum well spacing or row spacing.

[0049] In order to accurately characterize the flow gradient near the crack, the background mesh needs to be locally non-uniformly refined according to the crack parameters (half-length, orientation) set in step 1. The refinement rules are the key to ensuring the accuracy of the model in this invention.

[0050] For each well and its designed fractures, the infill zone is the area of ​​the fracture's expected extension and the midpoint of the line connecting the injection and production wells (where flow changes may be drastic). Specifically, this is done based on the fracture half-length... and azimuth Determine the coordinates of the two endpoints of the crack, and using the crack segment as the center, expand outwards to both sides by a specific width (e.g., 10% of half the crack length or a fixed distance such as 10 meters) to form a rectangular mesh reinforcement zone. Within the reinforcement zone, the minimum size of the mesh must satisfy the following constraints: in, The minimum side length of the grid on the x-axis. This represents the minimum side length of the grid along the y-axis. These are the main fracture half-lengths of the gas injection well, corner well, and side well, respectively. The first encryption coefficient has a value range of [10, 20].

[0051] The value of determines the number of mesh segments into which the crack is discretized. Its selection should be based on numerical convergence analysis. A common practice is to establish a simple single-crack model and gradually increase . (i.e., denser mesh) to observe changes in key output parameters (such as bottom hole pressure and fracture flow). When If, after exceeding a certain value, the change in output parameters is less than a preset tolerance (e.g., 1%), then the numerical solution is considered to have converged. This is considered a suitable value. For most gas-driven simulation problems, to ensure the accuracy of pressure distribution within the fracture, each fracture is typically required to be characterized by at least 10-20 mesh elements. Therefore, Setting the value range to [10, 20] is reasonable and has been verified. If the crack half-length difference is large, taking the minimum value for conservative encryption can ensure that even the shortest crack has sufficient resolution.

[0052] For example, suppose the fracture half-lengths of the three wells are respectively , , ,but .Pick ,but: That is, within the encrypted area, the grid size should not exceed 3.33 meters.

[0053] In regions far from wells and fractures (unrefined zones), to control the total number of meshes, the mesh size increases exponentially towards the model boundary, with a growth rate ranging from (1,2). A preferred value is between 1.2 and 1.5. Smaller values ​​result in smoother mesh transitions and better computational stability, but also faster mesh growth; larger values ​​provide better mesh number control, but excessive gradients may lead to numerical oscillations. The optimal value is typically determined through calculation, ensuring no abrupt changes in the pressure field during the transition zone.

[0054] To accurately simulate the flow state of the smallest symmetric element in an infinitely large well network, symmetric boundary conditions need to be applied to the boundary of the physical modeling domain, specifically: Set the boundaries parallel to the wellbore direction (the upper and lower boundaries in this example) as flow-free boundaries. This is because, due to symmetry, these boundaries are the boundaries of streamlines, and no fluid crosses them.

[0055] Set the boundaries perpendicular to the wellbore direction (in this example, the left and right boundaries) as either constant pressure gradient boundaries or connectivity boundaries. Constant pressure gradient boundaries: suitable for simulating steady-state seepage. The gradient direction is parallel to the wellbore direction, and the gradient value can be estimated based on the overall injection-production pressure difference and well network scale. Periodic connectivity boundaries: more accurately simulate infinitely large well networks, treating the grid cells opposite the boundary as adjacent units for flow calculations. This is the more recommended method, but its implementation is slightly more complex.

[0056] Preprocessing of boundary condition data: If a constant pressure gradient boundary is used, a reasonable background pressure gradient value needs to be estimated as input based on geomechanical analysis or regional pressure data. This value can be obtained through regression analysis of adjacent well pressure test data.

[0057] Based on the concept of an embedded discrete fracture model, the main fracture of each well in the well network is simplified into a linear high-conductivity channel. The embedding method and flow calculation rules of the linear high-conductivity channel in the matrix grid model are as follows: Each main fracture is discretized into a series of interconnected fracture elements based on its length and orientation angle. The discretization principle is to match the length of the fracture elements to the size of the local matrix mesh. Typically, this requires partitioning along the fracture path during discretization to ensure that the length of each generated fracture element is no greater than twice the minimum side length of the smallest matrix mesh along its extension direction, thus guaranteeing the accuracy of conductivity calculations.

[0058] Calculate the geometric intersection relationship between each crack element and the matrix mesh it passes through, and establish non-adjacent connections between the crack element and the matrix mesh it passes through based on the principle of embedded discrete crack model; Geometric calculations are used to determine whether crack element segments intersect with matrix mesh polygons. For 3D models, It is usually the length of the intersecting line segment multiplied by the effective thickness of the reservoir.

[0059] The formula for calculating the conductivity of the non-adjacent connections is as follows: in, Indicates the conductivity of non-adjacent connections; This represents the area of ​​intersection between the crack element and the matrix mesh it traverses; This represents the vertical distance from the center of the matrix grid to the crack cell it passes through.

[0060] This represents the effective permeability of the matrix mesh along the normal direction of the fracture surface; important preprocessing and calculations are required here. Obtaining anisotropic data is crucial, as tight reservoirs often exhibit permeability anisotropy. This requires core experiments (in different directions) or imaging logging interpretation to obtain permeability data in the x (parallel to the wellbore) and y (perpendicular to the wellbore) directions. and .

[0061] Direction conversion: For azimuth angles of... The crack, whose normal direction is The effective permeability of the matrix along this normal direction is then... It can be estimated based on the total tensor permeability or using the following formula (if the principal permeability direction is consistent with the coordinate axis): in, It is the angle between the crack normal direction and the x-axis. If anisotropy data is lacking, a conservative value can be taken as... It is the geometric mean or the smaller of the permeability in the two main directions.

[0062] For example, suppose a matrix mesh has a center coordinate of (10,10) and a crack element segment is from (9,9) to (11,11) (i.e., azimuth angle 45°).

[0063] Calculate the equation of the crack element and determine the perpendicular distance from the mesh center to the line segment. Calculations show that the perpendicular distance from point (10,10) to line segment y=x is 0 because the point lies on the line. In practical applications, it is almost impossible for the mesh center to be exactly on the crack; therefore, we assume... .

[0064] Assuming the fracture element intersects the grid by a length of 2m and the reservoir thickness H = 10m, then .

[0065] Assuming the fracture normal direction is 135°, the permeability in this direction is calculated based on anisotropic permeability. .

[0066] Calculate conductivity: In practical numerical simulators, fluid viscosity and unit conversions are usually combined to convert it into a practical conductivity value.

[0067] Flow connections also need to be established between discretized fracture units and between fracture units and the wellbore. This depends on the fracture conductivity. The center distance between adjacent crack elements is calculated. The formula is: in, The conductivity between adjacent crack elements. It is fluid viscosity. It is the distance between the centers of adjacent crack units. As input parameters, their units need to be uniform (e.g.) Convert to The fluid viscosity u needs to be calculated based on the composition of the injected gas, formation temperature, and pressure.

[0068] Step 3: Set the range of values ​​for fracturing design parameters, and generate N sets of training samples within this range using experimental design methods; for each set of training samples, perform rapid flow simulation using a parameterized model of well network and fracture, and calculate the spatial synergy quantification index of the corresponding sample based on the simulation results; the spatial synergy quantification index includes the gas shortest breakthrough time index, the effective sweep efficiency correction factor, and the displacement pressure gradient enhancement factor. In this embodiment, obtaining the training samples specifically includes: A complete set of fracturing design parameters to be optimized is defined as a 9-dimensional vector, denoted as X. The sampling range of the parameters for each dimension is set according to the constraints, for example: etc.

[0069] In a 9-dimensional parameter space, a Latin hypercube (LHD) design matrix containing N sample points is generated, where N is not less than 10 times the parameter dimension (i.e., N ≥ 90). N represents the number of training samples. This design ensures that the value range of each parameter is uniformly divided into N intervals; each interval contains only one sample point; and the projections of all parameters in their respective dimensions are uniformly distributed and uncorrelated. This empirical rule for setting N is to ensure that, in the higher-dimensional space, the sample points have sufficient ability to capture the basic shape and complexity of the response surface. For systems with strong nonlinearity, this can be appropriately increased to 10-20 times (i.e., N = 135-180).

[0070] The generation process of the Latin hypercube design matrix includes: Interval partitioning: Divide the range of values ​​for each parameter into N non-overlapping intervals (sub-intervals).

[0071] Random sampling: Select a value independently and randomly within each interval of each parameter.

[0072] Random pairing: Randomly arrange and pair the values ​​of all parameters selected within each interval to form N 9-dimensional sample points. This step ensures that the distribution of each parameter is uniform across the entire range, and that the projected distributions between any two parameters are uncorrelated, thereby maximizing space filling.

[0073] Because the dimensions and numerical ranges of the parameters vary greatly (e.g., half-length is tens of meters, angle is tens of degrees, and flow capacity is hundreds), the parameters must be preprocessed by normalization before inputting samples into the model or for subsequent machine learning training. In this embodiment, max-min normalization is used to linearly transform each parameter to the [0,1] interval.

[0074] Each row of the Latin hypercube design matrix is ​​combined with the basic geometric parameters to form a training sample; by traversing each row of the Latin hypercube design matrix, N training samples are obtained.

[0075] Performing the aforementioned rapid flow simulation specifically refers to running a simplified, non-implicit flow calculation process based on a well pattern and fracture parameterization model, and satisfying at least one of the following conditions: Steady-state flow simulation is employed under constant injection-production pressure differential (e.g., the difference between the bottom-hole pressure of the injection well and the bottom-hole flowing pressure of the production well). The pressure field distribution is obtained by solving a single-phase (usually a pseudo-pressure function) equation, or, for computational convenience, by using an approximate gas-phase flow steady-state pressure equation. The change in saturation over time is not simulated. This steady-state pressure field is used to calculate the conductivity field, streamline distribution, and pressure gradient. This forms the basis for subsequent calculations of all spatial synergy quantification indicators. An initial injection-production pressure differential is required. This value can be determined through well test analysis of the target block, formation pressure testing, or engineering analogy. The bottom-hole pressure of each well must be input as a boundary condition before simulation. Steady-state flow simulation is preferred in this embodiment.

[0076] Streamline simulation is employed, using the pressure field obtained from steady-state simulation. The streamline distribution from the injection well to each production well is calculated using the particle tracking method. One-dimensional tracer convection equations or leading-edge propagation equations (such as the Buckley-Leverett equations) are solved along each streamline to rapidly assess the gas breakthrough time and the location of the ripple front. This is several orders of magnitude faster than a full three-dimensional two-phase simulation.

[0077] The simulation timeframe is limited, focusing only on the early gas-bearing stage and not performing long-term simulations of the entire lifecycle. A realistic two-phase (oil-gas) unsteady-state simulation is conducted, but only from the start of gas injection to the first production well achieving gas (or to a short, pre-set early timeframe, such as 30 days). Early dynamics are directly obtained, but the computational load is greater than the previous two methods.

[0078] In this embodiment, the gas shortest breakthrough time index aims to quantify the short-circuit risk of the gas-driven channel. A larger value indicates a later gas breakthrough, potentially leading to better development results. The specific calculation process is as follows: Constructing the flow network graph: Each matrix mesh and crack element in the parametric model is abstracted as a node in the graph. The connections between nodes (matrix-matrix adjacent connections, crack-crack connections, and crack-matrix connections established through NNC) are abstracted as edges, and the weight of the edge is related to the flow resistance of the connection.

[0079] Identifying the set of effective flow paths: Using a shortest path search algorithm in graph theory (such as Dijkstra's algorithm), with the node containing the injection well as the node containing each production well as the destination, search for the path with the minimum weight. Here, the weight is defined as a measure of flow time or flow resistance. More specifically, the weights of edges i and j can be defined as follows: or (Derived from Darcy's Law and mass equilibrium), where It is the conduction rate between the two nodes. Here, edges i and j are equivalent to the paths described below, and the paths are treated as edges.

[0080] For each path, calculate the quasi-steady-state flow time of gas breakthrough along that path: in, Indicates the quasi-steady-state flow time; The matrix porosity is the average value of well logging interpretation or core analysis of the target layer. The gas saturation is the average value taken during gas front displacement; this is a parameter that needs to be preset. For gas-driven miscible or near-miscible displacement, it can be taken as... For immiscible flooding, the gas saturation corresponding to the residual oil saturation can be used; this value needs to be determined through phase permeation experiments or analogy. The value is the pore volume of the matrix mesh to which the i-th path belongs. If the starting point is a fracture element, its value is 0 (because its storage capacity is negligible). This represents the conductivity between two connected computational units corresponding to the i-th path. When the connection is between a crack unit and a matrix mesh, this conductivity value is the conductivity of a non-adjacent connection. When the connection is between two adjacent matrix meshes, this conductivity value is the interface conductivity between them. The pressure difference between the two ends of the i-th path is calculated through a single steady-state flow under the initial injection-production pressure difference. i is the path index, and N represents the total number of paths; the path refers to the route taken by gas flowing from the center of one computational cell to the center of the next directly connected computational cell, where the computational cell refers to a matrix grid or fracture cell. This formula is the integral of the reciprocal of the Darcy velocity along the path, approximately estimating the time required for the gas front to advance along that path.

[0081] This quasi-steady-state flow time formula, by discretizing the flow path into a series of cascaded matrix and fracture cells, approximates the time required for the gas front to advance from the injection well to the production well along a specific channel, based on Darcy's law and the principle of mass balance. The formula integrates reservoir properties such as porosity, gas saturation, and pore volume, the flow capacity characterized by conductivity, and the flow potential energy driven by the steady-state pressure difference. Essentially, it is a time-quantified expression of the integral of the flow resistance along the path under the quasi-steady-state assumption. It provides a rapid and quantitative method for evaluating channel quality based on a steady-state pressure field, avoiding the time-consuming and lengthy full-dynamic two-phase simulation to predict breakthrough time. This allows for the efficient identification of high-risk flow paths prone to premature gas breakthrough through extensive sample screening and optimization iterations.

[0082] The larger this calculated value, the longer it takes for gas to break through along this path, indicating greater resistance or a longer distance in the flow channel. The spatial configuration of the well network and fractures helps to delay gas channeling, leading to more stable displacement and higher oil recovery. Conversely, a smaller value indicates the existence of low-resistance short-circuit channels, allowing gas to quickly break through to the production well. This usually means that fractures may directly connect injection and production wells, or form highly conductive finger-like protrusions, resulting in a sharp decline in gas drive efficiency, a rapid increase in the gas-oil ratio, and a significant deterioration in development performance.

[0083] For each production well, the shortest gas breakthrough time exponent is defined as the minimum of the quasi-steady-state flow time across all paths from the injection well to that well. Ultimately, the shortest gas breakthrough time exponent for the entire system can be the minimum of this value across all production wells, or the value of the corner well (the farthest and latest to reach gas), to characterize the worst-case scenario.

[0084] The effective sweep efficiency correction factor is used to quantify the effect of fracturing fractures on the original theoretical sweep volume at the gas breakthrough point; a larger value indicates higher sweep efficiency. The specific calculation steps are as follows: Based on the flow field established by steady-state flow simulation, an inert tracer (concentration of 1) is continuously injected from the injection wells, and all production wells produce at a constant flow rate or constant pressure. This is a transient convection-diffusion simulation, but the calculation speed is much faster than the full-component simulation. A constant flow rate production method is preferred.

[0085] Real-time monitoring of tracer concentration in the produced fluid of each production well. When the tracer concentration of any production well first reaches a preset threshold... Stop the simulation; this moment is recorded as the moment when the air is seen. The value of is determined as follows: It is usually set to a small value, such as 0.01 (1%), representing the initial signal of gas breakthrough. Its specific value can be adjusted according to the oilfield's engineering definition of gas breakthrough.

[0086] At the moment of gas exposure, the total volume of all matrix grids in the statistical model with tracer mole fractions greater than zero is denoted as the effective swept volume.

[0087] The theoretical swept volume of a symmetrical element in a rhombic inverse nine-point well network under ideal piston displacement is calculated using the following formula: in, H represents the theoretical swept volume, and H represents the effective reservoir thickness.

[0088] The calculation of this theoretical swept volume provides a standardized geometric reference volume, enabling the acquisition of the actual effective swept volume through rapid simulation. This allows for comparison, thereby quantifying the degree to which the fracture system corrects for sweep efficiency (i.e., the effective sweep efficiency correction factor). ).like This indicates that the fracture network expands the affected area and enhances the reservoir's mobilization capacity; if This means that cracks may lead to gas leakage or the formation of dead oil zones, which would reduce the theoretical sweep efficiency.

[0089] The independent variable of this formula is the basic geometric parameter: injection-production well spacing. Row spacing and reservoir effective thickness The dependent variable is the theoretical swept volume. Theoretical sweep volume With well spacing Row spacing and effective thickness The volume increases linearly with the increase of the injection-production well. This is because, under the ideal piston displacement assumption, the swept volume is determined by the geometry of the controlled drainage area of ​​the injection and production wells. The well spacing and drainage spacing directly define the projected area of ​​this symmetrical unit on the plane. ), and effective thickness This determines the vertical extension volume of the planar region. Therefore, the sparser the well network ( and The larger the reservoir (or the thicker the reservoir) The larger the unit (the larger the capacity), the greater the volume of crude oil it can theoretically access and displace. However, it is worth noting that excessively large units... and In practice, this may lead to insufficient displacement pressure gradient, thereby reducing the actual sweep efficiency. This is precisely the contradiction that this invention aims to overcome by optimizing fracturing design. Theoretical Sweep Volume The larger the value, the greater the resource base controlled by the well network unit under ideal conditions, providing greater potential for fracturing optimization and gas displacement; the smaller the value, the stronger the inherent geological or well network constraints of the unit, and the relatively limited absolute space for optimization and improvement.

[0090] The displacement pressure gradient enhancement factor quantifies the degree to which fractures enhance the displacement dynamics within the reservoir; a larger value indicates a stronger displacement effect of the fractures on the matrix. The specific calculation steps are as follows: A steady-state flow simulation was conducted using a parameterized model of well grid and fracture to obtain the pressure values ​​of each matrix grid and each fracture element.

[0091] For each matrix mesh, identify all crack elements that are not adjacent to it; calculate the perpendicular distance between each crack element and the matrix mesh, and find the minimum perpendicular distance. The corresponding crack element is considered its associated crack element. All elements satisfying this condition are filtered out. The matrix mesh forms the zone directly affected by the cracks.

[0092] The range of values ​​for the distance threshold is: ,in This is the maximum value of the upper limit of the half-length of all cracks in step 1. This range roughly corresponds to the region where the pressure disturbance generated by the cracks can significantly affect the matrix flow.

[0093] For each matrix mesh within the crack-affected zone, the formula for the effective displacement pressure gradient is: in, This represents the effective displacement pressure gradient of the matrix grid; The pressure value associated with the crack element; This represents the pressure value of the matrix mesh; This represents the vertical distance between the matrix mesh and the associated crack element, with the gradient direction perpendicular to the crack wall.

[0094] The arithmetic mean of the effective displacement pressure gradients of all matrix meshes in the crack-affected zone is the average effective displacement pressure gradient.

[0095] When there are no fractures, the average pressure gradient between the injection well and the farthest corner well is taken as the reference pressure gradient, and the calculation formula is as follows: in, The reference pressure gradient; This refers to the bottom pressure of the gas injection well. The bottom-hole flowing pressure of the corner well furthest from the injection well; This is the straight-line distance between the bottom of the gas injection well and the farthest angle well. The formula for calculating the enhancement factor is: in, As a displacement pressure gradient enhancement factor, It is the arithmetic mean of the pressure of all grid-related fracture elements in the area directly affected by the fracture. This displacement pressure gradient strengthening factor is a dimensionless value. A value greater than 1 indicates that the presence of hydraulic fracturing significantly enhances the local displacement dynamics within the matrix rock block, and the larger the value, the more obvious the strengthening effect.

[0096] In this embodiment, the preliminary screening of N sets of training samples based on the spatial synergy quantification index specifically includes: The screening criteria are set as follows: the shortest gas breakthrough time index is greater than the preset critical time and the effective sweep efficiency correction factor is greater than the preset critical sweep efficiency.

[0097] The methods for determining the critical time and critical sweep efficiency are as follows: Engineering analogy method: Refer to historical fracturing gas drive well group data of the target block or similar blocks, and use the index value corresponding to the case with poor effect (such as early gas channeling) as the lower limit of the critical value.

[0098] Quantile method: Calculate the distribution of the shortest breakthrough time index and effective sweep efficiency correction factor for N samples, and use the lower quartile (25th percentile) as the critical value. For example, if the 25th percentile of the shortest breakthrough time index for N=100 samples is 60 days, then set the critical time to 60 days. This automatically selects the top 75% of samples as candidates. For example, set the critical time to 90 days and the critical sweep efficiency to 0.8. Samples that meet both conditions enter the candidate pool.

[0099] From the candidate pool that meets the screening criteria, Q groups of training samples are randomly selected according to a preset ratio to form a high-fidelity simulated sample subset. The preset ratio is usually 20%-30% to ensure that the subset size Q is on the order of tens (e.g., selecting 20-30 from 100), thus achieving a balance between computational resources and data representativeness.

[0100] For each training sample in the high-fidelity simulation sample subset, based on its corresponding fracturing design parameters and basic geometric parameters, the following operations are performed to construct a full-well-field numerical simulation model: Using a complete well group of a rhombic inverse nine-point well network as the simulation area, a three-dimensional structured mesh is constructed using a local mesh refinement method. Refinement is applied around the expected path of the main fracture to ensure that the dimensions of the mesh penetrating the fracture in the fracture extension direction meet the following requirements: ,in, For the extended dimension length; This is the second encryption coefficient. The second encryption coefficient can be slightly smaller than the first encryption coefficient (e.g., 8-12), because this is the refined model used for the final evaluation.

[0101] Using the fracture parameters from the training samples, the main fracture is characterized in a refined mesh using the equivalent conductivity assignment method. The permeability of each matrix mesh traversed by the fracture is modified to the ratio of fracture conductivity to fracture mechanical width. The fracture mechanical width needs to be preset here, typically estimated based on fracturing design (proppant particle size, proppant concentration) or microseismic inversion; a typical value is 0.005-0.01 meters.

[0102] A multi-component fluid model including methane, ethane, propane, and C7+ heavy fractions, along with the Peng-Robinson equation of state, is used to describe the phase changes of oil and gas during gas injection. Preprocessing experimental PVT data (expansion experiments, flash evaporation experiments, etc.) is required to fit the parameters of the equation of state. PVT data refers to pressure-volume-temperature data, which is a series of experimental measurements describing the phase behavior and physical properties of reservoir fluids (crude oil, condensate, natural gas, and injected gas) under specific temperature and pressure conditions. For gas injection displacement processes (especially miscible or near-miscible displacement), accurate PVT data is fundamental for accurately simulating fluid phase changes, volume coefficients, viscosity, and miscibility. The specific preprocessing steps are as follows: Extended compositional analysis: This involves further fractionating the C7+ heavy fraction of crude oil into single-carbon number (SCC) components (e.g., C7, C8, C9, ...) up to C45+ or higher. This is typically based on gas chromatography (GC) analysis data or a standard distribution function (e.g., exponential distribution). If only the total amount, molecular weight, and density of C7+ are available, a component splitting model (e.g., Whitson's) is required. The distribution model divides it into multiple single-carbon fractions and estimates the mole fraction, molecular weight, and density of each single-carbon fraction.

[0103] Grouping and merging of pseudo-groups: To reduce computational load, dozens of single-carbon arrays are grouped into fewer (usually 3-5) pseudo-groups. Common merging methods include merging by carbon number range and merging by volatility. For example, C7-C12 are grouped into pseudo-group 1, C13-C21 into pseudo-group 2, and C22+ into pseudo-group 3. Alternatively, cluster analysis can be performed based on the critical properties (critical temperature Tc, critical pressure Pc) and eccentricity factor of each component, grouping components with similar properties into one class. Dedicated PVT simulation software (such as PVTi, WinProp, etc.) or scripts can be used to perform the merging calculations.

[0104] Using the pre-processed pseudo-component system and its initially estimated physical properties, based on the Peng-Robinson (PR) equation of state, these parameters were adjusted using a nonlinear regression algorithm to achieve the best match between the calculated EOS and the pre-processed experimental data.

[0105] Initial reservoir pressure, temperature, and fluid saturation are set. For gas injection wells, the maximum injection rate is typically set under the constraint of maximum bottom hole pressure, while for production wells, the maximum production rate is set under the constraint of minimum bottom hole flowing pressure. The total simulation time is typically the target development cycle (e.g., 10-15 years).

[0106] The constructed full-well-field numerical simulation model is run, and the cumulative oil production and cumulative gas production for the entire development cycle are obtained by integration. Then, the cumulative gas-oil ratio, that is, the ratio of cumulative gas production to cumulative oil production, is calculated.

[0107] For the remaining training samples (NQ samples) not selected for the high-fidelity subset, their key metrics (cumulative oil production and cumulative gas-oil ratio) can be temporarily filled using K-nearest neighbor (KNN) interpolation of the simulated samples during the initial training of the surrogate model. For example, for an unsimulated sample, its K nearest simulated samples (e.g., K=3) are found in the parameter space, and the distance-weighted average of the key metrics of these neighboring samples is used as its initial filling value. This is only a starting point for initializing the surrogate model, and its inaccuracies will be corrected in subsequent active learning iterations.

[0108] Step 4: For each training sample, a numerical simulator is used to perform extended production dynamic simulation to calculate the key indicators of the gas drive development effect of each training sample within the complete development cycle, including cumulative oil production and cumulative gas-oil ratio; using basic geometric parameters and fracturing design parameters to be optimized as input features, and spatial synergy quantification indicators and key indicators as output labels, a machine learning surrogate model is trained.

[0109] In this embodiment, the construction and iterative verification process of the machine learning agent model specifically includes: Using basic geometric parameters and fracturing design parameters to be optimized as input features, the calculated spatial synergy quantification index and key index are used as output labels, and the machine learning surrogate model is initially trained based on N sets of training sample data.

[0110] To ensure the stability and convergence speed of machine learning model training, preprocessing of input features and output labels is essential. This preprocessing can involve max-min normalization or Z-score standardization of input features. For the cumulative oil recovery in the output labels, which is typically large and positive, logarithmic transformation or Z-score standardization can be applied. Logarithmic transformation compresses the data range and makes its distribution closer to a normal distribution. The cumulative gas-oil ratio undergoes the same preprocessing as the cumulative oil recovery. All spatial synergy indicators can be uniformly standardized using Z-score. Initial values ​​of cumulative oil recovery and cumulative gas-oil ratio generated by interpolation are considered valid data during the initial training, but their source (i.e., whether it is interpolated or a true value) needs to be recorded. Outliers in the output labels (such as maxima or minima caused by simulation non-convergence) should be checked. These can be identified using the 3-standard-deviation principle or box plot method, and these outlier samples should be temporarily marked and removed from the training set. The cause will be analyzed before deciding whether to restore or remove them.

[0111] Proxy models need to handle regression problems with multiple inputs, multiple outputs, nonlinearity, and potentially complex interactions. Possible algorithms include: Gaussian process regression: particularly suitable for scenarios with small samples (N<200) and where prediction uncertainty needs to be provided. This is one of the most commonly used models in active learning because it can directly output the mean and variance (uncertainty) of the predicted values.

[0112] Gradient boosting decision trees, such as XGBoost and LightGBM, are suitable for medium sample sizes (N>200) and possess strong nonlinear fitting capabilities and high computational efficiency. However, their native output does not directly provide uncertainty estimates, which must be obtained indirectly through ensemble methods (such as quantile regression).

[0113] Artificial neural networks, such as multilayer perceptrons (MLPs), are suitable for large sample sizes and can fit extremely complex functional relationships. Uncertainty estimation can also be obtained through ensemble learning or Bayesian neural networks.

[0114] Support Vector Regression: Suitable for small to medium sample sizes, it has certain advantages for high-dimensional problems, but its scalability and uncertainty estimation are not as good as GPR.

[0115] In this embodiment, given that the initial sample size N is typically in the range of 100-200 and active learning requires reliable uncertainty estimation, the preferred algorithm is Gaussian process regression as a surrogate model.

[0116] The construction and training of the Gaussian process regression (GPR) model are as follows: Model definition: GPR hypothesis function value It follows a Gaussian process (GP) that is completely composed of the mean function. Sum of covariance functions (kernel functions) Decision. Usually set. .

[0117] Kernel function selection and hyperparameters: Commonly used kernel functions include radial basis function (RBF) kernels and Matérn-like kernels (such as Matérn5 / 2). RBF kernels assume that the function is infinitely smooth, while Matérn5 / 2 kernels can better handle functions with a certain degree of roughness and are more robust in engineering applications.

[0118] If using the RBF kernel function: in, The function value represents two sample points in the input control (i.e., the fracturing design parameter space). and The similarity or correlation between the two points. The output is a non-negative scalar value. When two points are exactly the same ( The function value is at its maximum when two points are far apart, indicating perfect correlation; when the two points are far apart, the function value approaches 0, indicating no correlation. This function forms the basic element of the GPR covariance matrix, determining the smoothness, volatility, and other characteristics of the prediction function throughout the input space. and These are two distinct input sample vectors, containing basic geometric parameters and fracturing design parameters to be optimized. i and j are indices of the sample vectors. It represents the signal variance, which controls the fluctuation range of the function value. is the length scale of the q-th input dimension, controlling how smoothly this dimension smooths the image of the function; this is a key parameter, and a larger value indicates that the feature is less sensitive. q represents the dimension index of the input feature, and D represents the total number of dimensions of the input feature. This represents the value of the q-th feature of the i-th sample vector. This represents the value of the q-th feature of the j-th sample vector. This is the noise variance, used to account for observational noise in the fitted data (such as errors in the numerical simulation itself). For Kroneko function.

[0119] Model training (learning hyperparameters): Using a training dataset (N samples), the hyperparameters of the kernel function are optimized by maximizing the marginal likelihood function. This is a nonlinear optimization problem, which is usually solved using the conjugate gradient method or the quasi-Newton method.

[0120] Marginal likelihood: in, Let X represent the conditional probability density function, given the input data matrix X and hyperparameters. Under the given conditions, the probability of observing the output data vector; the output data vector is usually an n×1 column vector containing the true value of a specific output label for all training samples. It is the covariance matrix calculated by the kernel function. Let N represent the n×n identity matrix. N represents the number of training samples.

[0121] The training process involves finding the hyperparameter that maximizes the likelihood value. .

[0122] The core of active learning is to leverage the uncertainty of the model to guide the selection of new samples. For the Gaussian Proportional (GPR) model, given a new input point, its predicted output follows a Gaussian distribution. The standard deviation of the prediction is a measure of the uncertainty of the model's prediction for that point. A threshold for prediction uncertainty is defined. For example, the relative threshold method can be used: calculate the standard deviation of the prediction error for all samples in the current training set. Set the threshold to a multiple of that value, for example... The rationale for this setting is that it is desirable for the uncertainty of the new sample points to be significantly higher than the average error level of the current model.

[0123] Using the initially trained machine learning agent model, predictions are made for regions within the training sample space that have not undergone high-fidelity simulation, and their prediction uncertainty is estimated. Several new sample points with prediction uncertainties higher than a prediction uncertainty threshold are selected, and the spatial coherence quantification index of these new sample points is calculated and added to the training dataset to obtain an enhanced training dataset. Using the enhanced training dataset, the active learning process of the machine learning agent model is retrained until the prediction uncertainty of all new sample points is lower than the prediction uncertainty threshold, or the total number of high-fidelity simulations reaches a preset upper limit. At this point, the machine learning agent model is considered to have completed training.

[0124] The active learning iterative process aims to obtain a high-precision surrogate model with the fewest high-fidelity simulations. Specifically, the initial state consists of a pre-trained surrogate model and a training set containing Q high-fidelity samples and (NQ) interpolated samples.

[0125] Iteration steps (kth iteration): Finding candidate points: Within the entire parameter space defined in step 1 (not just LHD sample points), use optimization algorithms (such as genetic algorithms or Bayesian optimization) to find new points that maximize the acquisition function. In selecting the acquisition function: the most common approach is expected improvement, aiming to find the expected value of the point most likely to improve the current optimal solution. Maximizing prediction uncertainty: directly select the point with the largest prediction variance. This method is simple and direct, suitable for the objective of reducing global model uncertainty in this invention. For example (using maximum uncertainty): densely sample within the parameter space (e.g., generate 10,000 candidate points using Sobol sequences), using the current model... Predict the standard deviation for each point. Select the top B points with the largest standard deviations as candidates. The batch size B is typically 1-5, and processing can be done serially or in parallel.

[0126] High-fidelity simulation and data update: For the selected new points and their corresponding basic geometric parameters, perform the extended production dynamic simulation in step 3 to calculate their five true output labels. Add this true sample (i.e., input and output) to the training set; obtain the updated dataset.

[0127] Model retraining: Using the updated dataset, retrain the agent model to obtain... .

[0128] Convergence Criteria: Check if the stopping conditions are met. Condition 1 (Accuracy Satisfaction): The prediction uncertainty of new points is below a preset threshold. Condition 2 (Resource Exhaustion): The total number of high-fidelity simulations (initial Q + number of iterations added) reaches a preset upper limit, such as 200. Condition 3 (Performance Saturation): On an independent validation set (some initial LHD samples can be reserved for non-training), the improvement in model prediction accuracy (such as a decrease in root mean square error RMSE) is less than a small amount (such as 1%) for several consecutive iterations (such as 3).

[0129] Loop: If the stopping condition is not met, return to find candidate points for the next iteration; if the condition is met, terminate the iteration and output the finally trained surrogate model.

[0130] After the active learning iterations are completed, the final trained surrogate model must undergo independent performance validation to ensure its reliability. Specifically, when generating the initial LHD samples in step 3, a certain percentage (e.g., 10%-20%) of the samples should be reserved as an independent test set, not participating in any training. Alternatively, a new set of LHD samples can be generated in the parameter space and subjected to high-fidelity simulation as the test set. The final trained surrogate model is then used to predict the test set, and the results are compared with the actual high-fidelity simulation results to calculate the coefficient of determination, root mean square error (RMSE), mean absolute percentage error (MAE), and maximum absolute error (MAE). A coefficient of determination closer to 1 is better, indicating the model's ability to explain data variation. RMSE represents an absolute error measure. MAE represents a relative error measure, particularly suitable for cumulative oil recovery. MAE is used to assess worst-case scenarios. When the coefficient of determination for key indicators (e.g., cumulative oil recovery) is greater than 0.9 and the MAE is less than 10%, the surrogate model is generally considered reliable and usable. If these criteria are not met, it may be necessary to adjust the active learning strategy (e.g., the acquisition function), increase the initial sample size N, or check the data quality.

[0131] Step 5: With the optimization objectives of maximizing cumulative oil production and the shortest gas breakthrough time index, and with the cumulative gas-oil ratio controlled below the target threshold as a constraint, a multi-objective optimization algorithm is used to obtain the Pareto optimal solution set, and then the final recommended scheme is selected by combining key indicators.

[0132] In this embodiment, selecting the final recommended solution specifically includes: The optimization variables are defined as the fracturing design parameters of the target block. A scheme is constructed that maximizes two objective functions, aiming to find the scheme that maximizes the benefits of gas-driven development. This typically manifests in two conflicting aspects: maximizing oil production and minimizing gas production (controlling gas channeling). Therefore, the two objective functions are to maximize cumulative oil production and maximize the shortest gas breakthrough time exponent. Optimization must satisfy key production technology constraints. The most important constraint is controlling the production gas-oil ratio to prevent ineffective gas circulation; therefore, the constraint condition is a multi-objective optimization problem where the cumulative gas-oil ratio does not exceed the target gas-oil ratio.

[0133] Determining the target gas-oil ratio requires a comprehensive consideration of technical, economic, and safety factors. The specific determination method is as follows: One approach is the technology limit method. This involves referencing empirical data from the target block or similar reservoirs used in gas-drive development to determine a technological upper limit. For example, statistical analysis might show that when the production gas-oil ratio exceeds 1500 m³ / m³ (standard conditions), the wellhead back pressure increases significantly, the lift efficiency drops sharply, and the production system becomes difficult to operate stably. Therefore, a target gas-oil ratio of 1500 can be set.

[0134] For example, the economic break-even point method involves a simple economic evaluation. This involves calculating the operating cost per unit of oil production and the cost of injected gas. When the value of the produced gas (usually with low reinjection or sales value) cannot cover the costs of its separation, compression, and reinjection, or when net revenue begins to decline, the corresponding cumulative gas-oil ratio represents the economic limit. The formula simplifies to: Net Revenue = Oil Price × Oil Production - Gas Processing Cost × Gas Production - Fixed Costs. By setting the derivative of net revenue with respect to gas production to zero, the critical gas production can be obtained, leading to the critical cumulative gas-oil ratio. For instance, if the calculated economic critical cumulative gas-oil ratio is 1200 m³ / m³, this can be set as the target gas-oil ratio.

[0135] Considerations may be based on minimum miscibility pressure (MMP). For miscible or near-miscible flooding, a higher gas injection rate is required to maintain formation pressure above the MMP, which may lead to a higher production gas-oil ratio. In this case, the target gas-oil ratio can be set higher, such as 2000 m³ / m³, but this needs to be verified in conjunction with the aforementioned economic considerations.

[0136] Using the trained machine learning agent model as the fitness evaluation function, a non-dominated sorting genetic algorithm with an elitist strategy is used to iteratively solve the problem within the range of values ​​of the optimization variables to obtain a Pareto optimal solution set that satisfies the constraints. Each solution corresponds to a set of fracturing design parameters and their predicted key indicators and spatial synergy quantification indicators.

[0137] In this embodiment, four fracturing design modes are defined: Mode 1 is a symmetrical design where all wells have the same fracture parameters; Mode 2 is a long fracture design for injection wells, where the fracture length of the injection well is significantly greater than that of the production well; Mode 3 is a corner well with strong conductivity, where the corner well has the highest conductivity, enhancing far-end drainage; and Mode 4 is a side well with deflection angle, where the fracture azimuth of the side well deviates from the normal direction, delaying gas channeling. The basic conditions are set as follows: well network parameters: injection-production well spacing d = 200m, row spacing L = 150m; matrix permeability 0.1mD; reservoir thickness 10m, porosity 0.12; injection pressure 35MPa; production well flowing pressure 15MPa; target cumulative gas-oil ratio 1200. This basic condition setting applies to all modes. A total of 30 sets of sample data were collected under the four modes mentioned above. Cumulative oil production, minimum gas breakthrough time index, effective sweep efficiency correction factor, and cumulative gas-oil ratio were used as comparison indicators to compare the effects of asymmetric fracture parameter configurations in different well groups. Specific data are shown in the table below: Table 1: Comparison of the effects of asymmetric configuration of fracture parameters in different well groups Based on Table 1 above, Figure 2-5 It can be seen that the four fracturing design modes exhibit distinctly different development effects in the diamond-shaped inverse nine-point well network gas drive system of tight, low-permeability reservoirs, reflecting the complex coupling relationship between fracture parameters and well network spatial synergy. From the perspective of cumulative oil production trends, Mode 2 (long fracture in the injection well) demonstrates the highest oil production capacity in most groups, with its oil production generally exceeding that of Mode 1 (symmetric design) by approximately 15%-25%. This is mainly due to the significant increase in the half-length of the fracture in the injection well, which expands the gas injection sweep range and enhances displacement energy. However, this high oil production advantage comes at the cost of sacrificing gas channeling risk control. Its shortest gas breakthrough time index is the shortest among the four modes in most cases, shortening by an average of approximately 30%-40% compared to Mode 1. This indicates that while long fracture injection wells improve injection capacity, they also establish a high-speed channel for gas to rapidly penetrate to the production well.

[0138] In terms of gas breakthrough time control, Mode 3 (strong conductivity in corner wells) and Mode 4 (side well deflection angle) showed significant advantages. Mode 4, in particular, effectively extended the seepage path of gas propagation from the injection well to the side well by deviating the fracture azimuth angle of the side well from the wellbore normal direction. Its breakthrough time was extended by an average of 35%-50% compared to Mode 1, and in some groups even exceeded 130 days, demonstrating excellent gas channeling delay capabilities. Mode 3, by enhancing the conductivity of the corner well, improved the drainage capacity of the distant production wells, maintaining a high oil production rate while keeping the breakthrough time at a high level, reflecting the design logic of strong distant end and stable near end. It is worth noting that although the effective sweep efficiency of Modes 3 and 4 was slightly lower than that of Mode 2, it was still generally higher than or close to that of Mode 1, indicating that these two asymmetric designs still have a positive effect on improving sweep efficiency, although their mechanism of action differs from simply increasing the fracture length of the injection well.

[0139] From the perspective of the key economic indicator of cumulative gas-oil ratio, the four modes exhibit a clear gradient characteristic. Mode 2, due to early gas channeling and rapid gas circulation, generally exceeds the target threshold (1200 m³ / m³) by approximately 10%-20%, with some groups even reaching 1450 m³ / m³. This means that more injected gas is required to produce one cubic meter of crude oil, significantly reducing economic efficiency. Mode 1's cumulative gas-oil ratio mostly fluctuates around the target value, showing robust performance but lacking optimization potential. Modes 3 and 4, especially Mode 4, with their excellent breakthrough time control capabilities, generally have cumulative gas-oil ratios that are approximately 10%-25% lower than the target value, reaching as low as 850 m³ / m³, demonstrating higher gas utilization efficiency and better economic benefits. This result intuitively shows that by differentiating well group fracture parameters, the gas-oil ratio characteristics of the gas-drive system can be substantially improved, reducing development costs.

[0140] A comparative analysis of the four indicators reveals that no single model can achieve optimal performance across all indicators simultaneously, reflecting the inherent multi-objective trade-offs in fracturing well network gas drive design. While Model 2 maximizes oil production, it comes at the cost of high-risk gas channeling and a high gas-oil ratio; Model 4 optimally controls gas channeling and the gas-oil ratio, but compromises on oil production; Model 3 achieves a better balance between oil production, breakthrough time, and the gas-oil ratio, demonstrating the most robust overall performance. These images and data collectively reveal that in optimizing fracturing well networks for gas drive in tight, low-permeability reservoirs, it is crucial to make targeted choices among various asymmetric design strategies based on the specific geological conditions of the block, gas resource accessibility, and economic objectives, rather than simply adopting a symmetric or single-enhancement strategy.

[0141] A non-dominated sorting genetic algorithm with an elitist strategy is used. Within the range set in step 1, an initial population of P individuals is randomly generated. P is typically set to 10-20 times the dimension of the optimization variable. For a 9-dimensional problem, P can be 100-200.

[0142] For the t-th generation population, perform the following steps until the preset maximum number of iterations is reached. The maximum number of iterations is one of the optimization termination conditions, and can be set to 100-300, depending on the convergence speed.

[0143] For each individual in the population, it is input into the trained machine learning agent model to obtain its predicted spatial synergy quantitative indicators and key indicators.

[0144] For individuals that violate the constraints, a penalty function value is calculated. If the difference between the predicted cumulative gas-oil ratio and the target gas-oil ratio is greater than 0, the penalty function value is the product of this difference and a preset penalty coefficient; otherwise, the penalty function value is 0. The penalty function value is deducted from the original fitness or used as an additional cost. The penalty coefficient is a large integer used to push individuals that violate the constraints out of the feasible region. The penalty coefficient should be large enough that the corrected fitness of any individual that violates the constraints is significantly worse than that of a feasible solution. For example, a penalty coefficient of [value missing]... .

[0145] Based on the objective function value of each individual after penalty function correction, the current generation t population is non-dominated and sorted, dividing it into multiple non-dominated frontier levels. Simultaneously, the crowding degree of individuals within the same frontier level is calculated to maintain the diversity of the solution set. Using a binary tournament selection mechanism, two individuals are randomly selected from the population, and their non-dominated level and crowding degree are compared. The individual with the lower non-dominated level is preferred; if the non-dominated levels are the same, the individual with the higher crowding degree is selected. This process is repeated until a parent individual is selected, and simulated binary crossover and polynomial mutation operations are applied to generate the offspring population. The parent and offspring populations are merged, and the non-dominated sorting and crowding degree calculation are re-performed on the merged population. The top P individuals are selected to form the new generation population. During crossover, the crossover probability controls the generation of new individuals, typically between 0.7 and 0.9. During mutation, the mutation probability is used to maintain population diversity, typically taking the reciprocal of n.

[0146] When the maximum number of iterations is reached, the Pareto optimal solution set is output for all individuals in the final population that belong to the first non-dominated frontier.

[0147] Although the optimization process takes into account CGR constraints, it may still be necessary to ensure that the scheme has basic displacement efficiency and driving force.

[0148] The effectiveness filtering conditions are set as follows: the effective sweep efficiency correction factor is not less than a preset correction factor threshold, and the displacement pressure gradient enhancement factor is not less than a preset enhancement factor threshold. All individuals that do not meet the above effectiveness filtering conditions are removed from the Pareto optimal solution set; each individual corresponds to one solution. The set of remaining individuals after removal is defined as the engineering feasible solution set. The effective sweep efficiency correction factor threshold can be set to 0.8 or 0.9. Values ​​below this indicate that the fracture system as a whole has reduced sweep efficiency, and the solution may not be attractive. The displacement pressure gradient enhancement factor threshold can be set to 1 or 1.2. Values ​​below 1 indicate that the fracture has not effectively enhanced the matrix displacement dynamics, and the solution may be ineffective.

[0149] The feasible solution set may still contain a large number of solutions with similar parameter patterns. To facilitate decision-making, cluster analysis is needed to group solutions with similar parameter design patterns into one category.

[0150] Clustering characteristics: The K-means clustering algorithm is widely used because of its simplicity and efficiency.

[0151] Determining the number of clusters: The elbow rule is used. Calculate the sum of squared errors of the clustering results for different numbers of clusters (e.g., 2 to 10), and plot the K-SSE curve. Select the K value corresponding to the inflection point (elbow) of the curve. Generally, for the optimization of a well group, 3-5 typical schemes are sufficient to cover different design strategies (e.g., "long-slot strong conductivity", "short-slot high density", "asymmetric slotting", etc.).

[0152] Clustering operation: K-means clustering is performed on the feasible solution set of the project, and each solution is assigned a class label.

[0153] Within each cluster category, a representative scheme is selected according to the following priority: First priority: Among the schemes that satisfy the requirement that the cumulative gas-oil ratio is not greater than the target gas-oil ratio, select the scheme with the largest cumulative oil production.

[0154] Second priority: If none of the schemes in this category meet the cumulative gas-oil ratio constraint (indicating that this type of design mode is prone to high gas-oil ratio), then select the scheme with the smallest predicted cumulative gas-oil ratio (i.e. the least relative harm).

[0155] For the selected representative schemes, a key economic indicator is further calculated: unit gas recovery efficiency, which is the ratio of cumulative oil recovery to the cumulative gas injection volume corresponding to the achievement of the cumulative oil recovery. The representative schemes corresponding to the category with the highest unit gas recovery efficiency are output as the final recommended schemes, along with their corresponding cumulative oil recovery, cumulative gas-oil ratio, and unit gas recovery efficiency, to quantitatively characterize the gas drive development effect of these schemes.

[0156] Please see Figure 6 The present invention also provides a system for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs. This system is used to implement the aforementioned method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs, and includes: The basic parameters and constraint definition module is used to determine the basic geometric parameters of the diamond-shaped inverted nine-point well network of the target block, including the well spacing and row spacing of injection and production wells; and to set the range of fracturing design parameters to be optimized, including the half length of the main fracture of each well, the angle between the azimuth of the main fracture and the well row direction, and the conductivity of the main fracture. The well network and fracture model construction module is used to simplify the main fracture of each well in the well network into a linear high-conductivity channel and embed the main fracture into the matrix grid model that characterizes the matrix properties to form a well network and fracture parameterized model that reflects the flow coupling relationship between the fracture and the matrix. The spatial synergy assessment module is used to set the range of values ​​for fracturing design parameters and generate N sets of training samples through experimental design methods. For each set of training samples, rapid flow simulation is performed using a parameterized model of well network and fracture, and the spatial synergy quantification index of the corresponding sample is calculated based on the simulation results. The spatial synergy quantification index includes the gas shortest breakthrough time index, the effective sweep efficiency correction factor, and the displacement pressure gradient enhancement factor. The surrogate model training module is used to perform extended production dynamic simulation using a numerical simulator for each training sample, and calculate the key indicators of gas drive development effect during the complete development cycle of each training sample, including cumulative oil production and cumulative gas-oil ratio; using basic geometric parameters and fracturing design parameters to be optimized as input features, and spatial synergy quantification indicators and key indicators as output labels, to train the machine learning surrogate model. The multi-objective optimization and decision-making module is used to maximize the cumulative oil production and the shortest gas breakthrough time index as optimization objectives, and control the cumulative gas-oil ratio below the target threshold as a constraint. It uses a multi-objective optimization algorithm to obtain the Pareto optimal solution set, and then combines key indicators to select the final recommended solution.

[0157] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0158] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0159] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0160] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs, characterized in that, The specific steps include: Step 1: Determine the basic geometric parameters of the diamond-shaped inverted nine-point well pattern in the target block, including the well spacing and row spacing of injection and production wells; set the fracturing design parameters to be optimized, including the main fracture half-length of each well, the angle between the main fracture azimuth and the well row direction, and the conductivity of the main fracture. Step 2: Simplify the main fracture of each well in the well network into a linear high-conductivity channel, and embed the main fracture into the matrix grid model that characterizes the matrix properties to form a well network and fracture parameterization model that reflects the flow coupling relationship between the fracture and the matrix. Step 3: Set the range of values ​​for fracturing design parameters, and generate N sets of training samples within this range using experimental design methods; for each set of training samples, perform rapid flow simulation using a parameterized model of well network and fracture, and calculate the spatial synergy quantification index of the corresponding sample based on the simulation results; the spatial synergy quantification index includes the gas shortest breakthrough time index, the effective sweep efficiency correction factor, and the displacement pressure gradient enhancement factor. Step 4: For each training sample, a numerical simulator is used to perform extended production dynamic simulation to calculate the key indicators of the gas drive development effect of each training sample within the complete development cycle, including cumulative oil production and cumulative gas-oil ratio; using basic geometric parameters and fracturing design parameters to be optimized as input features, and spatial synergy quantification indicators and key indicators as output labels, a machine learning proxy model is trained. Step 5: With the optimization objectives of maximizing cumulative oil production and the shortest gas breakthrough time index, and with the cumulative gas-oil ratio controlled below the target threshold as a constraint, a multi-objective optimization algorithm is used to obtain the Pareto optimal solution set, and then the final recommended solution is selected by combining key indicators. The method for determining the value range of the fracturing design parameters to be optimized is as follows: For the range of the main fracture half-length, based on meeting the minimum engineering requirements for forming an effective modified zone during fracturing construction, its upper limit is set according to the well spacing and row spacing of the injection and production wells in the aforementioned diamond-shaped inverted nine-point well network, and its upper limit satisfies the constraint conditions. ,in, This indicates the upper limit of the half-length of the main crack; and This is the proportionality coefficient, with a value range of [0, 0.5]. The distance between injection and production wells, The spacing between injection and production wells; this constraint ensures that the main fractures do not physically intersect within the well network unit; The angle between the azimuth of the main fracture and the direction of the well network is called the fracture azimuth angle. Its value range is set based on the well network layout and the principle of minimizing gas breakthrough risk. ,in, The angle between the normal direction of the line connecting the gas injection well and the corner well and the well row direction. Allowable design deviation angle; The conductivity range of the main fracture is determined dimensionlessly based on the matrix permeability of the target reservoir, and its relationship is expressed as follows: ,in, This represents the dimensionless conductivity of the crack. This indicates the conductivity of the main fracture. The matrix permeability of the target reservoir is represented by the dimensionless fracture conductivity falling within the range of [1, 100] as a constraint condition to determine the range of values ​​for the conductivity of the main fracture. This indicates the half-length corresponding to the main crack; Based on the function and spatial relationship of different well locations in the diamond-shaped inverted nine-point well network, the fracturing design parameters are differentiated into groups. Specifically, the gas injection wells located at the center of the well network are divided into gas injection well groups, and their fracturing design parameter groups are set, denoted as... ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the gas injection well; the four corner wells located at the four corners of the well network are divided into corner well production well groups, and their fracturing design parameter groups are set, denoted as... ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the corner well; the four side wells located at the midpoints of the four edges of the well network are divided into a side well production well group, and their fracturing design parameter group is set, denoted as . ,in, These are the main fracture half-length, fracture azimuth angle, and conductivity of the side well, respectively. The fracturing design parameters to be optimized include the fracturing design parameter groups of the above-mentioned gas injection well groups, corner well production well groups, and side well production well groups; A minimum symmetric element of a rhombic inverse nine-point well grid is selected as the physical modeling domain to construct the matrix grid model; Within the physical modeling domain, a background Cartesian mesh is initialized, and a local coordinate system is established with the location of the gas injection well as the origin. Based on the main fracture design parameters of each well, the mesh is non-uniformly refined. Specifically, for any fracture, along the wellbore direction and its perpendicular direction, the mesh is refined at the midpoint of the line connecting the injection and production wells and in the expected fracture extension area. The mesh size satisfies the following constraints: in, The minimum side length of the grid on the x-axis. This represents the minimum side length of the grid along the y-axis. These are the main fracture half-lengths of the gas injection well, corner well, and side well, respectively. is the first encryption coefficient, with a value range of [10, 20]; In regions far from wells and fractures, the mesh size increases toward the model boundary in a geometric progression, with the growth rate ranging from (1,2]. Symmetric boundary conditions are applied to the boundaries of the physical modeling domain, specifically: the boundary parallel to the wellbore direction is set as a flow-free boundary, and the boundary perpendicular to the wellbore direction is set as a constant pressure gradient boundary or a connectivity boundary, in order to simulate the flow state of the symmetric unit in an infinitely large periodic well network. Based on the concept of an embedded discrete fracture model, the main fracture of each well in the well network is simplified into a linear high-conductivity channel. The embedding method and flow calculation rules of the linear high-conductivity channel in the matrix grid model are as follows: Each main crack is discretized into a series of interconnected crack units based on its length and orientation angle. Calculate the geometric intersection relationship between each crack element and the matrix mesh it passes through, and establish non-adjacent connections between the crack element and the matrix mesh it passes through based on the principle of embedded discrete crack model; The conductivity calculation formula for the non-adjacent connections is as follows: in, Indicates the conductivity of non-adjacent connections; The effective permeability of the matrix mesh in the direction normal to the fracture surface; This represents the area of ​​intersection between the crack element and the matrix mesh it traverses; This represents the vertical distance from the center of the matrix grid to the crack cell it passes through.

2. The method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs according to claim 1, characterized in that, Obtaining the training samples specifically includes: A complete set of fracturing design parameters to be optimized is defined as a 9-dimensional vector; the sampling range of the parameters in each dimension is set according to the constraints. In a 9-dimensional parameter space, generate a Latin hypercube design matrix containing N sample points, where N is not less than 10 times the parameter dimension, i.e., N≥90; N represents the number of training samples. Each row of the Latin hypercube design matrix is ​​combined with the basic geometric parameters to form a training sample; by traversing each row of the Latin hypercube design matrix, N training samples are obtained. Performing the aforementioned rapid flow simulation specifically refers to running a simplified, non-implicit flow calculation process based on a well pattern and fracture parameterization model, and satisfying at least one of the following conditions: Steady-state flow simulation is used, which solves for the pressure field distribution under a constant injection-production pressure difference, without simulating the change of saturation over time; this simulation is used to calculate the conductivity field, streamlines, and pressure gradient. Streamline simulation is used to calculate the streamline distribution based on the pressure field obtained from steady-state simulation, and one-dimensional tracer or leading-edge propulsion equations are solved along the streamlines to assess the breakthrough time and the affected area. The simulation timeframe is limited, focusing only on the early gas emergence stage, and does not include long-term simulations of the entire life cycle.

3. The method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs according to claim 1, characterized in that, The calculation process for the shortest breakthrough time exponent of the gas is as follows: The set of effective flow paths connecting the injection wells and each production well is identified using the shortest path search algorithm in graph theory; each path consists of several matrix grid blocks or fracture elements, or is composed of matrix grids and fracture elements connected in series. For each path, calculate the quasi-steady-state flow time of gas breakthrough along that path: in, Indicates the quasi-steady-state flow time; For matrix porosity, The gas saturation is the average value taken during gas frontal displacement. The value is 0 for the pore volume of the matrix mesh to which the i-th path belongs; if the starting point is a crack element, its value is 0. This represents the conductivity between two connected computational units corresponding to the i-th path; when the connection is between a crack unit and a matrix mesh, this conductivity value is the conductivity of a non-adjacent connection; when the connection is between two adjacent matrix meshes, this conductivity value is the interface conductivity between them. The pressure difference between the two ends of the i-th path is calculated through a single steady-state flow under the initial injection-production pressure difference condition; i is the path index, and N represents the total number of paths; the path refers to the path through which gas flows from the center of one computational cell to the center of the next directly connected computational cell, and the computational cell refers to the matrix grid or fracture cell; For each production well, the minimum value of all its corresponding paths is denoted as the shortest gas breakthrough time index for that production well. The calculation steps for the effective sweep efficiency correction factor specifically include: In the parameterized model of well network and fracture, a tracer is continuously injected from the injection wells, while all production wells produce at a constant flow rate or constant pressure. The simulation stops when the tracer concentration in the produced fluid of any production well first reaches the preset threshold. This moment is recorded as the gas breakthrough time. At the moment of gas exposure, the total volume of all matrix grids in the statistical model with tracer mole fractions greater than zero is denoted as the effective swept volume. The theoretical swept volume of a symmetrical element in a rhombic inverse nine-point well network under ideal piston displacement is calculated using the following formula: in, Where H is the theoretical swept volume and H is the effective reservoir thickness; in, For effective sweep efficiency correction factor, For effective sweep volume; The specific steps for calculating the displacement pressure gradient enhancement factor include: A steady-state flow simulation was performed using a parameterized model of well grid and fracture to obtain the pressure values ​​of each matrix grid and each fracture element. For each matrix mesh, identify all crack elements that are not adjacent to it; calculate the vertical distance between these crack elements and the matrix mesh one by one, and regard the crack element corresponding to the minimum vertical distance as the associated crack element of the matrix mesh; The formula for calculating the effective displacement pressure gradient of the matrix mesh is: in, This represents the effective displacement pressure gradient of the matrix grid; The pressure value associated with the crack element; This represents the pressure value of the matrix mesh; This represents the vertical distance between the matrix mesh and the associated crack element, with the gradient direction perpendicular to the crack wall. All matrix meshes with a vertical distance less than a preset distance threshold are selected to form the crack's direct influence zone; the distance threshold range is [value missing]. ; The arithmetic mean of the effective displacement pressure gradients of all matrix meshes in the crack-affected zone is the average effective displacement pressure gradient. When there are no fractures, the average pressure gradient between the injection well and the farthest corner well is taken as the reference pressure gradient, and the calculation formula is as follows: in, The reference pressure gradient; This refers to the bottom pressure of the gas injection well. The bottom-hole flowing pressure of the corner well furthest from the injection well; This is the straight-line distance between the bottom of the gas injection well and the farthest angle well. The displacement pressure gradient enhancement factor is the ratio of the average effective displacement pressure gradient to the reference pressure gradient.

4. The method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs according to claim 1, characterized in that, The preliminary screening of N sets of training samples based on spatial synergy quantification indicators specifically includes: The screening criteria include that the gas shortest breakthrough time index is greater than the preset critical time, and the effective sweep efficiency correction factor is greater than the preset critical sweep efficiency. From the training samples that meet the screening criteria, Q groups of training samples are randomly selected according to a preset ratio as a subset of high-fidelity simulation samples for extended production dynamic simulation. For each training sample in the high-fidelity simulation sample subset, based on its corresponding fracturing design parameters and basic geometric parameters, the following operations are performed to construct a full-wellfield numerical simulation model: Using a complete well group of a rhombic inverse nine-point well network as the simulation region, a three-dimensional structured mesh was constructed using a local mesh refinement method. Refinement was applied around the expected path of the main fracture to ensure that the mesh dimensions of the fracture-penetrating mesh in the fracture extension direction met the following requirements: ,in, For the extended dimension length; This is the second encryption factor; Using the crack parameters of the training samples, the main crack is characterized in the refined mesh by the equivalent conductivity assignment method. Specifically, the permeability of each matrix mesh through which the crack passes is modified to the ratio of the conductivity of the main crack to the mechanical width of the crack. A multi-component fluid model including methane, ethane, propane, and C7+ heavy fractions was adopted, and the Peng-Robinson equation of state was used to describe the phase changes of oil and gas during the gas injection process; initial reservoir pressure, temperature, and fluid saturation were set. Set the gas injection rate or bottom hole flowing pressure of the gas injection well, set the production rate or bottom hole flowing pressure of the production well, and define the total simulation time as the complete development cycle; Run the constructed full-well-field numerical simulation model and record the daily oil production and daily gas production of all production wells at each simulation time step; after the simulation is completed, calculate the cumulative oil production and cumulative gas-oil ratio of the training sample by time integration; For the remaining training samples that were not selected into the simulated sample subset, their key metrics were temporarily filled using the preliminary predictions of the built machine learning surrogate model or interpolation of similar samples for the initial training of the surrogate model.

5. The method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs according to claim 4, characterized in that, The construction and iterative verification process of the machine learning agent model specifically includes: Using basic geometric parameters and fracturing design parameters to be optimized as input features, the calculated spatial synergy quantitative index and key index are used as output labels, and the machine learning proxy model is initially trained based on N sets of training sample data. Define a prediction uncertainty threshold; use the machine learning proxy model after initial training to predict regions in the training sample space that have not undergone high-fidelity simulation and estimate their prediction uncertainty; select several new sample points whose prediction uncertainty is higher than the prediction uncertainty threshold, calculate the spatial coherence quantification index of these new sample points, and add them to the training dataset to obtain an enhanced training dataset. Using the enhanced training dataset, the active learning process of the machine learning agent model is retrained until the prediction uncertainty of new sample points is lower than the prediction uncertainty threshold, or the total number of high-fidelity simulations reaches the preset upper limit. At this point, the machine learning agent model is considered to have been trained.

6. The method for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs according to claim 1, characterized in that, The specific steps to select the final recommended solution include: The optimization variables are defined as the fracturing design parameters of the target block; a multi-objective optimization problem is constructed, which includes two maximization objective functions, namely maximizing the cumulative oil production and maximizing the shortest gas breakthrough time exponent, and the constraint that the cumulative gas-oil ratio is not greater than the target gas-oil ratio. Using the trained machine learning agent model as the fitness evaluation function, a non-dominated sorting genetic algorithm with an elitist strategy is used to iteratively solve the problem within the range of values ​​of the optimization variables to obtain a Pareto optimal solution set that satisfies the constraints. Each solution corresponds to a set of fracturing design parameters and their predicted key indicators and spatial synergy quantitative indicators. Set the effectiveness filtering conditions: the effective sweep coefficient correction factor is not less than the preset correction factor threshold, and the displacement pressure gradient enhancement factor is not less than the preset enhancement factor threshold. From the Pareto optimal solution set, remove all individuals that do not meet the above validity filtering conditions, and each individual corresponds to a solution; the set of the remaining individuals after removal is defined as the engineering feasible solution set; Clustering algorithms are used to cluster the fracture parameter patterns of the feasible solutions in the engineering project set and group them into several typical scheme categories. Within each scheme category, the scheme with the largest cumulative oil production that satisfies the constraint that the cumulative gas-oil ratio is not greater than the target gas-oil ratio is selected as the representative scheme of the category. If no solution satisfies this condition, the scheme with the smallest cumulative gas-oil ratio is selected as the representative scheme of the category. For the representative schemes of each final category, calculate their unit gas recovery efficiency; that is, the ratio of the cumulative oil recovery to the cumulative gas injection volume corresponding to the cumulative oil recovery is the unit gas recovery efficiency; output the representative scheme of the category corresponding to the maximum unit gas recovery efficiency as the final recommended scheme, and output its corresponding cumulative oil recovery, cumulative gas-oil ratio and unit gas recovery efficiency together to quantitatively characterize the gas drive development effect of the scheme.

7. A system for predicting the gas drive effect of fracturing well patterns in tight, low-permeability oil reservoirs, characterized in that, The aforementioned system for predicting the gas drive effect of fractured well patterns in tight, low-permeability reservoirs is used to implement the method for predicting the gas drive effect of fractured well patterns in tight, low-permeability reservoirs as described in any one of claims 1-6, comprising: The basic parameters and constraint definition module is used to determine the basic geometric parameters of the diamond-shaped inverted nine-point well network of the target block, including the well spacing and row spacing of injection and production wells; and to set the range of fracturing design parameters to be optimized, including the half length of the main fracture of each well, the angle between the azimuth of the main fracture and the well row direction, and the conductivity of the main fracture. The well network and fracture model construction module is used to simplify the main fracture of each well in the well network into a linear high-conductivity channel and embed the main fracture into the matrix grid model that characterizes the matrix properties to form a well network and fracture parameterized model that reflects the flow coupling relationship between the fracture and the matrix. The spatial synergy assessment module is used to set the range of values ​​for fracturing design parameters and generate N sets of training samples through experimental design methods. For each set of training samples, rapid flow simulation is performed using a parameterized model of well network and fracture, and the spatial synergy quantification index of the corresponding sample is calculated based on the simulation results. The spatial synergy quantification index includes the gas shortest breakthrough time index, the effective sweep efficiency correction factor, and the displacement pressure gradient enhancement factor. The surrogate model training module is used to perform extended production dynamic simulation using a numerical simulator for each training sample, and calculate the key indicators of gas drive development effect during the complete development cycle of each training sample, including cumulative oil production and cumulative gas-oil ratio; using basic geometric parameters and fracturing design parameters to be optimized as input features, and spatial synergy quantification indicators and key indicators as output labels, to train the machine learning surrogate model. The multi-objective optimization and decision-making module is used to maximize the cumulative oil production and the shortest gas breakthrough time index as optimization objectives, and control the cumulative gas-oil ratio below the target threshold as a constraint. It uses a multi-objective optimization algorithm to obtain the Pareto optimal solution set, and then combines key indicators to select the final recommended solution.

Images