Rapid simulation and layered injection-production real-time optimization method suitable for complex oil reservoir

Abstract

The invention discloses a rapid simulation and layered injection-production real-time optimization method suitable for a complex oil reservoir, and belongs to the field of unconventional oil reservoir numerical simulation, and the method comprises the following steps: 1, constructing a GMM-XGB hybrid model based on a Gaussian mixture model and an XGBoost algorithm; 2, constructing a layered connected network model based on a physical model; step 3, carrying out automatic history fitting by adopting an ES-MDA algorithm; and 4, establishing a layered injection and production real-time optimization objective function mathematical model, and performing layered injection and production real-time optimization based on an SHADE algorithm. According to the method, the recovery efficiency of the complex oil reservoir can be improved, the stable production period of an oil field is prolonged, and efficient utilization of resources and maximization of economic benefits are achieved.

Description

Technical Field

[0001] The invention belongs to the field of numerical simulation of unconventional oil reservoirs, and in particular relates to a fast simulation and stratified injection and production real-time optimization method suitable for complex oil reservoirs. Background Art

[0002] Reservoirs with fractures and faults have obvious geological complexity, which makes their effective development and simulation calculations face great challenges. Fractures will aggravate the lateral and vertical heterogeneity of the formation, develop water injection advantage channels in the plane, and increase the difference in interlayer utilization in the vertical direction. Therefore, it is a difficult problem to develop well group injection and production in reservoirs with fractures in a balanced and efficient manner. At the same time, the presence of faults will lead to the isolation or loss of local connectivity of the reservoir, which will affect the recovery efficiency of crude oil and the effect of water drive. In three-dimensional space, due to the discontinuity and irregularity of fractures and faults, traditional reservoir simulation methods require expensive costs to characterize fractures and faults, and then perform more accurate simulation calculations. Summary of the invention

[0003] In order to solve the above problems, the present invention proposes a fast simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs. The Gaussian mixture model and XGBoost algorithm are integrated to construct a GMM-XGB hybrid model to predict the injection and production volume of small layers, and the characterization of faults and fractures is embedded when constructing the flow network model, successfully constructing a set of layered connected network models applicable to complex oil reservoirs. Then, the ES-MDA algorithm considering formation constraints is used to automatically perform historical fitting inversion model well connectivity and the dynamic influence of fractures and faults. Based on this, an adaptive mechanism based on successful history and a SHADE algorithm with an archive mechanism are introduced, and the economic net present value is used as the optimization target. Considering multiple engineering constraints, real-time optimization of layered injection and production is achieved.

[0004] The technical solution of the present invention is as follows: A fast simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs comprises the following steps: Step 1: Construct a GMM-XGB mixture model based on the Gaussian mixture model and XGBoost algorithm; Step 2: Construct a hierarchical connected network model based on the physical model; Step 3: Use ES-MDA algorithm for automatic history matching; Step 4: Establish a mathematical model of the objective function for real-time optimization of stratified injection and production, and perform real-time optimization of stratified injection and production based on the SHADE algorithm.

[0005] Furthermore, in step 1, the GMM-XGB hybrid model uses a Gaussian mixture model to perform cluster analysis on the original sample set to obtain the probability distribution of the small layer injection and production volume, and integrates the probability distribution into the original sample set as the input of the XGBoost algorithm to predict and output the small layer injection and production volume; the specific working process of the GMM-XGB hybrid model is as follows: Step 1.1, collect the longitudinal splitting data of injection and production volume, and construct the original sample set; The vertical splitting data of injection and production volume include porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability extreme difference, porosity extreme difference, bottom hole flowing pressure, single well injection and production volume and water content; Step 1.2: Preprocess the original sample set. The specific process is as follows: First, data cleaning is performed. For missing values ​​of random missing types, the natural neighbor interpolation method in MATLAB software is used for interpolation according to the geodetic coordinates of the well points. For outliers, deletion or transformation processing is selected according to the actual situation. Transformation processing includes smoothing or scaling. For categorical data, the number of base classes of categorical data needs to be determined first. If there are few base classes, one-hot encoding is used for processing; if there are many base classes, label encoding or target encoding is used for processing. Then, the cleaned data was normalized using the Min-Max normalization method; Finally, the normalized sample set was subjected to correlation analysis; for linear relationships in the data, the Pearson correlation coefficient was used for correlation analysis; for sample sets with strong nonlinear relationships and no normal distribution, the Spearman rank correlation coefficient was used for correlation analysis; Step 1.3, using Gaussian mixture model to perform cluster analysis on the preprocessed original sample set to obtain the probability distribution of small layer injection and production volume; Step 1.4, the probability distribution of the small layer injection and production volume is fused with the preprocessed original sample set to obtain a fused sample set, and the fused sample set is divided into a training set and a validation set; Step 1.5: Input the fused sample set into the XGBoost algorithm and use the Bayesian optimization method to train until the requirements are met, and finally output the injection and production volume of all small layers in all wells and all time steps.

[0006] Furthermore, the specific process of step 1.3 is as follows: Step 1.3.1: Based on the preprocessed original sample set, a Gaussian mixture model is constructed; the Gaussian mixture model assumes that the injection and production volume data points come from multiple different Gaussian distributions, and the probability density function of the Gaussian distribution is: ; in, is the probability density function of Gaussian distribution; is the data point of injection and production volume; is the mean vector; is the covariance matrix; is the data variance; is the data dimension; is an exponential function with base e; Step 1.3.2: Use the expectation maximization algorithm to iteratively solve the established Gaussian mixture model until convergence. The specific process is as follows: First, initialize the parameters of each Gaussian distribution in the Gaussian mixture model; Then, the responsiveness of all injection and production volume data points is calculated. The responsiveness calculation formula is as follows: ; in, is responsiveness; is a binary indicator variable, indicating Injection and production volume data points Whether Gaussian components are generated; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; is the index variable of the Gaussian component, which is used to traverse all Gaussian components in the mixed Gaussian model for summation; is the total number of Gaussian components in the mixed Gaussian model; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; The parameters of each Gaussian distribution are updated according to the calculated responsiveness. The update formula is as follows: ; ; ; in, , , Respectively The first iteration The mean vector, covariance matrix, and prior probability of the Gaussian components; For the Responsiveness at iterations; is the total number of data points of injection and production volume; is the transpose symbol; Calculate the likelihood function of the model, the formula is as follows: ; in, is the likelihood function; Finally, the model parameters are continuously updated iteratively until the parameters of the Gaussian mixture model converge. At this time, the training ends and the final responsiveness is output. The responsiveness is the probability distribution of the injection and production volume of the small layer.

[0007] Furthermore, in step 1.5, Bayesian optimization is an automated hyperparameter tuning method, which establishes a proxy model and continuously adjusts the hyperparameters according to the historical performance of the model to find the optimal hyperparameter configuration; the specific process is: Step 1.5.1, initialize the Bayesian optimization algorithm framework, initialize the evaluation of the objective function in the framework and set the hyperparameter search space of the GMM-XGB hybrid model; the objective function uses the mean absolute error, and the formula is as follows: ; in, is the mean absolute error; For the The true value of samples; For the The predicted value of samples; is the total number of samples; Hyperparameters include learning rate, maximum tree depth, subsample ratio, and number of tree structures. The search space of learning rate is set to [0.01, 0.2], the search space of maximum tree depth is set to [5, 10], the search space of subsample ratio is set to [0.5, 1], and the search space of number of tree structures is set to [100, 500]. Step 1.5.2, iterative training, to find the optimal hyperparameters of the model in the search space; firstly, sample collection is performed, that is, a set of initial points are selected as the initial hyperparameters in the model hyperparameter combination space to be optimized, and then the Bayesian optimization algorithm evaluates the next set of hyperparameters to be selected based on the existing hyperparameters and objective function values; finally, the hyperparameter combination with the largest objective function is selected as the optimal hyperparameter combination; Step 1.5.3: Use the selected optimal hyperparameter combination to train the GMM-XGB hybrid model on the training set, calculate the mean absolute error between the predicted value and the true value of the GMM-XGB hybrid model on the validation set, and obtain the objective function value; Step 1.5.4: After obtaining the new hyperparameter-objective function pair, update the surrogate model. The surrogate model will continuously adjust its estimate of the distribution of the objective function in the hyperparameter space based on the new data. Step 1.5.5: Determine whether the pre-set termination condition is met; if the termination condition is not met, return to step 1.5.2 to collect samples again and continue iterating; when the termination condition is met, select a set of hyperparameters that optimize the objective function value from all evaluated hyperparameter combinations, and use this set of optimal hyperparameters to retrain the GMM-XGB hybrid model on the training set until the pre-set number of training times is reached; after the training is completed, the final GMM-XGB hybrid model for prediction is obtained.

[0008] Furthermore, the specific process of step 2 is as follows: Step 2.1, the oil reservoir is divided into several well points and connecting pipes, and an inter-well connecting network is constructed; in the inter-well connecting network, each well point represents a water injection well or oil production well on a small layer, and the well points are connected to each other through connecting pipes; the connecting pipes are characterized by two parameters: inter-well conductivity and connecting volume; The conductivity calculation formula is as follows: ; in, Well and well The conductivity between Well and well The seepage area between Well and well The average value of the permeability between Well and well The distance between The formula for calculating the connected volume is as follows: ; in, Well and well The connected volume between them; Well and well The average value of the effective thickness between Step 2.2, discretize the fracture into a group of uniformly arranged high permeability nodes; discretize the fault into two groups of uniformly arranged nodes, located on both sides of the fault, and characterize the flow capacity of the fluid near the fault by heterogeneous conductivity; embed the fracture and fault nodes into the well-to-well connection network to construct a dynamic flow channel; Step 2.3: Based on the hierarchical well connectivity relationship, multiple connectivity units are set in each layer and discretized into one-dimensional grids. Then, the grid connectivity relationship is mapped into a standardized model file to construct a physically driven hierarchical connectivity network model.

[0009] Furthermore, the specific process of step 2.3 is as follows: First, the reservoir is divided into multiple small layers according to its geological characteristics, and each small layer corresponds to an independent physical model. In each layer, multiple connected units are set based on geological data and well test data. The flow characteristics in each layer are refined by one-dimensional grid discretization, and each grid unit is equivalent to a part of the connected unit. The flow process of the reservoir is converted into a discretized numerical calculation problem; the basic idea of ​​grid discretization is to convert the continuous fluid flow process in each connected unit into a discrete system composed of multiple small units, and the permeability, volume and flow rate of each small unit are solved by numerical methods; the connection relationship between discrete units is mapped through a standardized model file to generate a hierarchical connected network model.

[0010] Furthermore, the specific process of step 3 is as follows: Step 3.1: Construct the history matching objective function as: ; in, is the historical fitting objective function value; is the model fitting parameter vector; Calculate the results for the model simulation; is the actual observation data; is the covariance matrix of the observation error; Step 3.2, select the conductivity between discrete grids, the volume of discrete grids, the heterogeneous conductivity and the number of characteristic points of fractures, the heterogeneous conductivity and the number of characteristic points of faults, and the phase permeability parameter as control variables to form a model fitting parameter vector; ; in, The initial moment grid and Grid The conductivity between The initial moment grid Volume; is the characteristic point of the crack at the initial moment and feature points exist Heterogeneous conductivity across the layer; For the initial moment The number of characteristic points of cracks on the layer; is the characteristic point of the fault at the initial moment and feature points exist Heterogeneous conductivity across the layer; For the initial moment The number of characteristic points of the fault on the layer; , , , are different phase permeability parameters; In the process of reservoir history matching, physical constraints are imposed while minimizing the objective function, including: ; in, For Grid and Grid conductivity; For Grid Volume; is the number of grids; is the total volume of the reservoir; Step 3.3: Use the ES-MDA algorithm, combined with historical production data and physical constraints, to achieve automatic inversion and real-time update of fracture and fault parameters. The specific process is as follows: Firstly, the initial heterogeneous conductivity distribution of fractures and faults is set in the reservoir numerical model, and the number of characteristic points is defined. Then, in the history fitting process, the actual production data of oil wells (daily oil production per well) is used as observation information to perform multiple data assimilation on the model. Through the iterative update mechanism of the ES-MDA algorithm, the matrix parameters as well as the parameters of fractures and faults are continuously adjusted.

[0011] Furthermore, the specific process of step 4 is as follows: Step 4.1: Select the economic net present value as the objective function of real-time optimization of injection and production; the following is the mathematical model of the optimization objective function: ; in, is the economic net present value; is the controlled variable, representing the regulation scheme of stratified water injection volume and stratified liquid production volume during the optimization process; for oil prices; For the Time step well Oil production; Cost of sewage treatment; For the Time step well The water production; is the water injection cost; For the Time step well The amount of water injected; is the total number of time steps; is the number of producing wells; is the number of water injection wells; Step 4.2: Set engineering constraints, including injection-production ratio constraint, bottom hole pressure constraint, maximum water cut constraint, and single well maximum liquid volume constraint. The specific constraints are as follows: ; In the formula, is the minimum injection-production ratio; is the maximum injection-production ratio; For the Time step well The amount of liquid produced or water injected; For the Time step well Minimum bottom hole flowing pressure; For the Time step well Maximum bottom hole flowing pressure; For the Time step well Bottom hole pressure; For the Time step well The moisture content; For the Time step well The ultimate moisture content; For the Time step well The amount of liquid produced or water injected; For the Time step well The maximum liquid production or maximum water injection volume; Step 4.3: Introduce the adaptive mechanism based on success history and the SHADE algorithm with archive mechanism to optimize reservoir injection and production.

[0012] Furthermore, in step 4.3, the specific process of the SHADE algorithm for reservoir injection and production optimization is as follows: Step 4.3.1, adjusting the injection and production flow rate to generate a new injection and production plan; Step 4.3.2, combine different injection and production schemes and calculate the economic net present value of each injection and production scheme combination; Step 4.3.3: Select the injection-production scheme combination that meets the engineering constraints and has the largest economic net present value as the individual to enter the next generation. The individual will perform steps 4.3.1 and 4.3.2 again according to the adaptive mechanism until the number of iterations is reached; after the iteration, the final injection-production scheme combination is the optimal injection-production scheme.

[0013] Beneficial technical effects brought by the present invention: The present invention adopts a Gaussian mixture model for cluster analysis, and based on the cluster analysis results, uses the XGBoost algorithm to train the data that integrates cluster features to obtain a GMM-XGB hybrid model. Compared with the conventional XGBoost method, the GMM-XGB hybrid model of the present invention can better discover the distribution pattern of split data, initialize the XGBoost model and thus reduce the overfitting problem. When establishing a flow network model, faults and cracks are embedded, which can characterize the shielding effect of faults and the conduction advantage of cracks, thereby improving the simulation accuracy of the flow path of fluids in complex oil reservoirs and optimizing the injection and production control effect. Compared with the traditional DE algorithm, the SHADE algorithm has more advantages in finding the global optimal solution. The model constructed by the present invention can improve the recovery rate of complex oil reservoirs, extend the stable production period of oil fields, and achieve efficient utilization of resources and maximize economic benefits. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 The present invention is a flow chart of the method for vertical production splitting and injection-production real-time optimization applicable to complex oil reservoirs.

[0015] Figure 2 The present invention is a flow chart for predicting the amount of injection and production fluid in a small layer.

[0016] Figure 3 Schematic diagram of a hierarchical connected network model in an embodiment of the present invention.

[0017] Figure 4 It is a schematic diagram of the historical matching results of the first iteration of sublayer 1 of production well 1 in an embodiment of the present invention.

[0018] Figure 5 It is a schematic diagram of the historical matching results of the second iteration of sublayer 1 of production well 1 in an embodiment of the present invention.

[0019] Figure 6 It is a schematic diagram of the historical matching results of the third iteration of sublayer 1 of production well 1 in an embodiment of the present invention.

[0020] Figure 7 This is a comparison diagram of the water content curves of the blocks after historical fitting in the embodiment of the present invention.

[0021] Figure 8 This is a comparison chart of the daily oil production curve of the small layer 1 after historical fitting in the embodiment of the present invention.

[0022] Fig. 9 This is the injection system control diagram of the water injection well 5 sublayer 1 in the embodiment of the present invention.

[0023] Fig.10 This is the injection system control diagram of the water injection well 5 sublayer 2 in the embodiment of the present invention.

[0024] Fig.11 This is a comparison chart of the cumulative oil production before and after optimization in the embodiment of the present invention.

[0025] Fig.12 This is a comparison chart of moisture content before and after optimization in the embodiment of the present invention. DETAILED DESCRIPTION

[0026] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments: like Figure 1 As shown, a fast simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs includes the following steps: Step 1: Construct a GMM-XGB hybrid model based on the Gaussian mixture model and the XGBoost algorithm; the GMM-XGB hybrid model uses the Gaussian mixture model to perform cluster analysis on the original sample set to obtain the probability distribution of the injection and production volume, and integrates the probability distribution into the original sample set as the input of the XGBoost algorithm to predict the output of the small layer injection and production volume. The specific working process of the GMM-XGB hybrid model is as follows: Step 1.1, collect the longitudinal splitting data of injection and production volume, and construct the original sample set; The vertical splitting data of injection and production volume include porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability extreme difference, porosity extreme difference, bottom hole flowing pressure, single well injection and production volume and water content and other factors affecting the liquid volume of small layers.

[0027] Step 1.2: Preprocess the original sample set, including data cleaning and correlation analysis; When constructing a sample set, it is necessary to ensure the integrity and consistency of the data. If the sample set has many data sources, there may be some noise and inconsistencies, so data cleaning must be performed to provide a reliable basis for subsequent modeling and analysis. For missing values, first determine the type of missing values, and select the appropriate missing value processing method based on the type of missing values. For outliers, first detect outliers and outliers, and choose to delete or transform them according to the actual situation. For categorical data in the sample set, first determine the base class of the categorical data to select the encoding method. Normalizing the cleaned data can adjust the data scale to better meet the requirements of the machine learning algorithm. In order to find variables that are strongly correlated with the small-layer injection and production fluid volume data, it is necessary to perform correlation analysis on the normalized sample set. For sample sets with strong nonlinear relationships and no normal distribution, the Spearman rank correlation coefficient is used for correlation analysis. The specific process is: Step 1.2.1, missing values ​​are mainly divided into three types: missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR). After analysis, the missing of reservoir data is mainly due to the missing of field testing or recording. This kind of missing often means that several parameters are missing at the same time, so it is mainly random missing. According to the geodetic coordinates of the well points, the natural neighbor interpolation ("natural") method in MATLAB software is used for interpolation. This method is based on Delaunay triangulation and achieves an effective balance between linear and cubic.

[0028] Step 1.2.2, outliers and outliers are important factors that affect data quality, so they need to be effectively detected and processed. First, use the box plot to detect outliers on the sample set and identify values ​​that deviate from the normal distribution. Specifically calculate the first quartile, median, and third quartile. After sorting the data from small to large, the data at the 25% position is the first quartile, the data at the 50% position is the median, and the data at the 75% position is the third quartile. The box range is the range of most of the data, and the calculation formula for the box is: (1); in, For the box; is the third quartile; is the first quartile; Outliers are usually data points outside the upper and lower edges. The specific calculation formulas for the upper and lower edges are: (2); (3); in, is the upper edge; for the lower edge; For detected outliers, you can choose to delete or transform them according to the actual situation. For example, for significant outliers, you may need to delete them directly; for some edge data, you can choose to smooth or scale the data to make it conform to the normal data distribution. The key to this step is to balance the rigor of outlier processing with the fidelity of the data.

[0029] Step 1.2.3: Some features in the sample set may be categorical data, such as geological type, injection and production type, etc. First, it is necessary to determine the number of base classes of categorical data. If there are fewer base classes, One-Hot Encoding can be used for processing; if there are more base classes, Label Encoding or Target Encoding is more efficient. Choose the appropriate encoding method according to different data types and task requirements. The appropriate encoding method helps improve the training effect and prediction ability of the model.

[0030] Step 1.2.4. Inconsistent data scales may affect the performance of machine learning algorithms, especially when using distance-based algorithms. To solve this problem, data normalization is a very important step. Normalization reduces the dimensional differences between different features by mapping the data to a uniform scale range (such as [0,1] or [-1,1]), allowing machine learning algorithms to converge better. In the normalization process, the commonly used method is Min-Max normalization. Min-Max normalization can accelerate model training and improve the stability and accuracy of the model. The specific formula is: (4); in, is the first in the normalized data set Column data; all columns after normalization constitute a data set , which is a matrix containing all features and labels; is the first in the data set before normalization Column data; For the The minimum value in the column data; For the The maximum value in the column data.

[0031] Step 1.2.5. In order to find variables that are strongly correlated with the small layer injection and production fluid volume data, it is necessary to perform correlation analysis on the normalized sample set. First, calculate the correlation between each feature to understand which variables have a significant relationship with the injection and production fluid volume. For linear relationships in the data, the Pearson correlation coefficient can be used for analysis. However, for some data with strong nonlinear relationships and no normal distribution, the Pearson correlation coefficient may not be able to effectively reflect the relationship between variables. Therefore, it is a better choice to use the Spearman rank correlation coefficient for analysis. The Spearman rank correlation coefficient can measure the monotonic relationship between variables, not just the linear relationship, thereby more comprehensively capturing the dependency between variables. The specific calculation formula is as follows: (5); in, is the Spearman rank correlation coefficient; For the The difference between the rank values ​​of the sample and the small layer injection and production volume data; is the total number of samples.

[0032] Step 1.3: Use Gaussian mixture model to perform cluster analysis on the preprocessed original sample set to obtain the probability distribution of small layer injection and production volume; the specific process is as follows: Step 1.3.1: Based on the preprocessed original sample set, construct a Gaussian mixture model (GMM). The Gaussian mixture model assumes that the injection and production volume data points come from multiple different Gaussian distributions.

[0033] The probability density function of the Gaussian distribution is: (6); in, is the probability density function of Gaussian distribution; is the data point of injection and production volume; is the mean vector; is the covariance matrix; is the data variance; is the data dimension; is an exponential function with base e; Step 1.3.2: Use the expectation-maximization (EM) algorithm to iteratively solve the established Gaussian mixture model until convergence, aiming to maximize the log-likelihood function of the observed data. The specific process of the EM algorithm is as follows: First, initialize the parameters of each Gaussian distribution in the Gaussian mixture model; Then, the responsiveness of all injection and production volume data points is calculated. The responsiveness calculation formula is as follows: (7); in, is responsiveness; is a binary indicator variable, indicating Injection and production volume data points Whether Gaussian components are generated; For the The prior probability of the Gaussian components is also called the mixing weight; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; is the index variable of the Gaussian component, which is used to traverse all Gaussian components in the mixed Gaussian model for summation; is the total number of Gaussian components in the mixed Gaussian model; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; The parameters of each Gaussian distribution are updated according to the calculated responsiveness. The update formula is as follows: (8); (9); (10); in, , , Respectively The first iteration The mean vector, covariance matrix, and prior probability of the Gaussian components; For the Responsiveness at iterations; is the total number of data points of injection and production volume; is the transpose symbol; The likelihood function of the model is calculated based on the calculated mixed weights. The likelihood function represents the fit of the model to the data. The likelihood function is used to determine whether the parameters of the model have converged. That is, if the likelihood function value between two iterations is less than 10 -5 , then the parameters of the model are determined to have converged. The likelihood function calculation formula is as follows: (11); in, is the likelihood function; Finally, the model parameters are continuously updated iteratively until the parameters of the Gaussian mixture model converge. At this time, the training ends and the final responsiveness is output. The responsiveness is the probability distribution of the injection and production volume of the small layer.

[0034] When using this method to deal with small sample problems such as splitting, GMM clustering has significant advantages over direct training with the XGBoost algorithm: from the perspective of data understanding and feature mining, the feature information of small sample data is limited, and it may be difficult to fully capture the data rules using the XGBoost algorithm directly. GMM can obtain the probability distribution and different clusters of split data, each cluster represents a subset of data with similar features, which helps to deeply understand the intrinsic characteristics of the data. Through clustering, hidden features that were originally difficult to detect can be mined, providing the XGBoost algorithm with richer and more discriminative feature information, so that the model can learn the key patterns in the data more accurately; in terms of improving the generalization ability of the model, the direct use of the XGBoost algorithm for split data is prone to overfitting, that is, the model overfits the training data and performs poorly when facing new data. After GMM clustering obtains different clusters and the probability distribution of split data, the data distribution of each cluster is relatively more concentrated and similar. Using the XGBoost algorithm for training in each cluster separately allows the model to learn and optimize more carefully according to the characteristics of different clusters. This clustering training method enables the model to better adapt to the diversity of data, reduce dependence on specific samples, and prevent overfitting, thereby enhancing the model's generalization ability on new data and making it more robust when facing new data.

[0035] Step 1.4: The probability distribution of the small layer injection and production volume output contains the potential distribution pattern of the liquid volume, and the constructed original sample set contains more features related to the liquid volume. Therefore, the probability distribution of the small layer injection and production volume output by the Gaussian mixture model in the previous step is fused with the preprocessed original sample set to obtain a fused sample set, and the fused sample set is divided into a training set and a validation set. The purpose of fusing the two is to combine the feature information of the probability distribution with the feature information in the actual sample to enhance the model's prediction ability for the distribution of injection and production volume. The fusion process mainly uses the probability distribution output by the Gaussian mixture model as an additional input feature to form a new training data set together with other feature data samples.

[0036] In order to reduce overfitting, the present invention adopts the 10-fold cross-validation method for division. Specifically, the 10-fold cross-validation divides the fusion sample set into 10 subsets, uses 9 subsets for training in each iteration, and the remaining 1 subset is used as a validation set. This process is repeated 10 times, and each subset will be used as a validation set once. Finally, the GMM-XGB hybrid model will be evaluated on all validation sets to ensure that the model performs evenly on different data subsets, thereby obtaining a more robust GMM-XGB hybrid model.

[0037] Step 1.5: Input the fused sample set into the extreme gradient boosting (XGBoost) algorithm for training. In order to improve the performance of the XGBoost algorithm, the Bayesian optimization method is used for training until the requirements are met, and finally the injection and production volume of all wells and all time steps in all small layers is output. By training different decision trees multiple times and combining them weightedly, the model performance is gradually improved. The specific process is as follows: Step 1.5.1. Bayesian optimization is an automated hyperparameter tuning method that establishes a proxy model and continuously adjusts the hyperparameters based on the historical performance of the model to find the optimal hyperparameter configuration. The key to Bayesian optimization is to build a proxy model (usually a Gaussian process, GP) to approximate the objective function. Initialize the evaluation of the objective function and the setting of the hyperparameter search space of the GMM-XGB hybrid model in the Bayesian optimization algorithm framework. The commonly used objective function is the mean absolute error (MAE), and the formula is as follows: (12); in, is the mean absolute error, For the The true value of the samples, For the The predicted value of the sample is the estimate of the true value by the GMM-XGB mixed model.

[0038] Commonly used hyperparameters in the GMM-XGB hybrid model include learning rate, maximum depth of tree structure, subsample ratio, number of tree structures, etc. When initializing the model, the search space of hyperparameters needs to be set: the search space of learning rate is [0.01, 0.2], the search space of maximum depth of tree structure is [5, 10], the search space of subsample ratio is [0.5, 1], and the search space of number of tree structures is [100, 500].

[0039] Step 1.5.2, after initialization, the GMM-XGB hybrid model is iteratively trained to find the optimal hyperparameters of the model in the search space. First, sample collection is performed, that is, a set of initial points are selected as the initial hyperparameters in the model hyperparameter combination space to be optimized (search space). Then the Bayesian optimization algorithm evaluates the next set of hyperparameters to be selected based on the existing hyperparameters and objective function values. The Bayesian optimization algorithm automatically balances utilization and exploration through the acquisition function (calculating expected improvement) and the Gaussian process (calculating prediction variance). This balance makes Bayesian optimization more efficient than grid search or random search, and is particularly suitable for hyperparameter tuning scenarios with high computational costs.

[0040] Step 1.5.3, then train and evaluate the GMM-XGB hybrid model. Use the selected optimal hyperparameter combination to train the GMM-XGB hybrid model. Train the GMM-XGB hybrid model on the training set, evaluate the performance of the GMM-XGB hybrid model on the validation set, calculate the mean absolute error between the predicted value and the true value of the GMM-XGB hybrid model on the validation set, and obtain the objective function value.

[0041] Step 1.5.4: After obtaining the new hyperparameter-objective function pair, update the surrogate model (usually a Gaussian process model). The surrogate model will continuously adjust the estimate of the distribution of the objective function in the hyperparameter space based on the new data, making the approximation of the objective function more accurate.

[0042] Step 1.5.5, check whether the preset termination condition is met. The common termination condition is set to reach the preset maximum number of iterations. If the termination condition is not met, return to step 1.5.2 to collect samples again and continue iteration; when the termination condition is met, select the set of hyperparameters that makes the objective function value optimal from all evaluated hyperparameter combinations. Use this set of optimal hyperparameters to retrain the GMM-XGB hybrid model on the training set until the preset number of training times is reached. After the training is completed, the final model for prediction or analysis is obtained. The performance of the model can then be evaluated on an independent test set to verify the generalization ability of the model.

[0043] The GMM-XGB hybrid model trained with Bayesian optimization will achieve good results on the given training set and validation set, and can accurately predict the distribution of injection and production fluid volume in small layers. Ultimately, the model will output the injection and production fluid volume of all wells at each time step, and refine the liquid volume distribution of each small layer. These output results will provide guidance for oilfield management and optimization, help reasonably allocate injection and production fluid volume, and improve the development efficiency of oilfields.

[0044] Step 2: Construct a hierarchical connected network model based on the physical model. The injection and production volume of small layers and other dynamic and static parameters of the reservoir are used as input, and the "feature point-connected network" method is proposed. The fractures and faults are abstracted as a series of feature points and embedded into the interconnected network between injection and production wells. The connected network that is more in line with the actual geological structure is reconstructed. The coupling of the matrix-fracture-fault multi-scale flow network is realized through hierarchical discretization technology. Then all the connections in the network are discretized into grids, mapped to a two-dimensional rectangular coordinate grid and input into the simulator for solution. The specific process is: Step 2.1: First, the reservoir can be equivalent to a well-to-well network formed by a series of well points and layered well-to-well connecting pipes. The well points are connected to each other through connecting pipes to form an inter-well connecting network. The complexity of this network can be flexibly adjusted by controlling the well spacing and well angle, so that the model can adapt to reservoirs of different scales and complexities. In this network, the inter-well pipe is regarded as a connecting unit, which is characterized by two parameters: inter-well conductivity and connecting volume. By calculating these parameters, the characteristics of fluid flow in the reservoir can be effectively captured.

[0045] In the process of reservoir development, studying the fluid flow characteristics of the reservoir and the connectivity between injection and production wells is a key link in optimizing oilfield production and improving oil recovery. In order to quickly and accurately simulate this process, an equivalent method can be used to regard the reservoir as a well-to-well interconnection network consisting of a series of well points and layered well-to-well interconnection pipelines. This network model simplifies and abstracts the complex physical processes of the reservoir, making it possible to analyze fluid flow and well-to-well interactions in a simpler and more intuitive way.

[0046] Specifically, in this network model, the reservoir is divided into several well points and connecting pipes. Each well point represents an injection well or oil production well on a small layer, and these well points are connected by a series of connecting pipes. Each inter-well pipe is equivalent to a flow channel, through which fluid can be transferred between the well points connected. The characteristics of each pipeline are characterized by two main parameters: inter-well conductivity and connected volume. Among them, inter-well conductivity reflects the ability of fluid to flow in the pipeline, and the connected volume represents the reserves between the well point and the surrounding well points. These two parameters are the core parameters for studying reservoir fluid dynamics and inter-well flow.

[0047] The conductivity calculation formula is as follows: (13); in, Well and well The conductivity between Well and well The seepage area between ; Well and well The average permeability between ; Well and well The distance between .

[0048] The formula for calculating the connected volume is as follows: (14); in, Well and well The connected volume between them; Well and well The average value of the effective thickness between.

[0049] The connecting pipe structure between well points is a model that is highly dependent on the specific conditions of the reservoir. The two geometric parameters, well spacing and well angle, directly determine the topological structure of the network. By adjusting these parameters, the complexity of the model can be flexibly controlled to adapt to reservoir conditions of different scales and different geological conditions. For example, in an oil reservoir, if the distance between well points is close and the angle is small, the connectivity between wells will be strong, and the exchange and conduction capacity of fluids will also be enhanced; while when the well spacing is far or the angle between wells is large, the flow and conduction effect of fluids are relatively weak. Therefore, controlling the well spacing and well angle enables the model to better simulate the actual well flow, thereby providing a scientific basis for production optimization.

[0050] In addition, this well-to-well network model can simplify complex reservoir conditions into a mathematical model, thereby providing support for subsequent real-time optimization of injection and production. By independently calculating each interconnected unit (inter-well pipeline), the flow characteristics of the fluid in the reservoir and the mutual influence between wells can be effectively captured.

[0051] Step 2.2, traditional flow network models are usually only roughly characterized by equivalent permeability parameters. The "feature point-connected network" method is proposed to abstract fractures and faults into discrete feature points and embed them into the well-to-well connected network to construct a dynamic flow channel that is closer to the actual geological structure. By giving feature points differentiated conductivity properties (such as high permeability of fractures and shielding effect of faults), the model can quantitatively characterize the high-speed seepage effect of fractures on injection and production fluids and the blocking effect of faults on flow paths.

[0052] In oil reservoirs, fractures and faults are two important geological features that significantly affect the flow path of oil and gas and the connectivity between wells, which in turn affects the recovery effect of the oil reservoir. In order to solve this problem, the influence of fractures and faults needs to be specially considered in reservoir simulation. The "feature point-connectivity network" method is proposed. By embedding feature points representing fractures and faults on the basis of the original well-to-well connectivity network model, the fluid flow of the oil reservoir can be described more accurately, and a basis for optimizing injection and production management can be provided.

[0053] Fractures are a type of high-speed flow channel in oil reservoirs, and their role in oil reservoirs is very significant. Fractures can greatly increase the regional permeability of oil reservoirs and become the main channel for crude oil flow. Due to the presence of fractures, the flow path of the fluid in the reservoir will change and the flow rate will increase significantly. In the original connected network, well points are connected through relatively low permeability channels, but in oil reservoirs with fractures, fractures are a high-speed channel and usually require special treatment.

[0054] To this end, on the basis of network segmentation, the fractures can be discretized into a set of uniformly arranged high permeability nodes. By introducing these fracture nodes on the basis of the original well points and network connecting pipes, the connectivity relationship between the well points and the fracture nodes is reconstructed. The high permeability characteristics of the fracture nodes mean that the fluid flows faster in these nodes. Therefore, the connectivity between these nodes is usually strong, the heterogeneous conductivity is large, and the flow of fluids is often more concentrated. This method can simulate the impact of fractures on reservoir permeability more quickly and accurately, and further improve the dynamic simulation and optimization capabilities of the reservoir.

[0055] Faults are another type of geological structure that affects the flow characteristics of oil reservoirs. They usually cause the reservoir to be separated or sealed. The presence of faults often results in restricted fluid flow between different areas of the reservoir, especially forming a certain barrier to the migration and accumulation of crude oil. In actual reservoirs, the nature of faults varies. Some faults are completely closed, resulting in the inability of oil and gas on both sides to flow to each other; while some faults are partially conductive and have a certain degree of permeability.

[0056] In the reservoir model considering the influence of faults, the fault is usually regarded as a closed structure. The present invention discretizes it into two groups of evenly arranged nodes, located on both sides of the fault, and characterizes the flow capacity of the fluid near the fault by heterogeneous conductivity. Through this discretization method, the closing effect of the fault on the connectivity between wells can be accurately reflected, and fault nodes can be formed in the network. When the well point and these fault nodes re-establish the connectivity relationship, the connectivity between the well point and the nodes on both sides of the fault will be affected by the fault. For closed faults, there is almost no fluid conduction between the well point and the fault, while for conductive faults, a certain flow connection may be maintained between the well point and the fault, but the permeability difference of the fault usually needs to be considered.

[0057] By discretizing faults into nodes and reconstructing the connectivity between well points and faults, the flow barriers or restrictions caused by faults in the reservoir can be simulated more realistically, which is of great significance for production optimization, injection and production scheduling, and pressure management of the reservoir.

[0058] Step 2.3, based on the hierarchical well-to-well connectivity relationship, multiple connected units are set in each layer and one-dimensional grid discretization is performed, and then the grid connection relationship is mapped into a standardized model file to construct a physically driven hierarchical connected network model. In the model construction, the coupling of the matrix-fracture-fault multi-scale flow network is realized through hierarchical discretization technology: for the flow of the matrix, one-dimensional grid discretization is used to describe the conventional seepage between wells; for the flow of fractures and faults, the conductivity of the fractures and the shielding characteristics of the faults are characterized by the heterogeneous conductivity parameters between the characteristic points. This coupling mechanism not only retains the computational efficiency of the flow network model, but also realizes the refined description of the complex flow path, solving the simulation deviation problem caused by the simplification of the flow path in the traditional method in fractured reservoirs. Finally, the pressure, saturation, injection and production rate and water content of each injection and production well are solved by the numerical simulator.

[0059] First, the reservoir is divided into multiple small layers according to the geological characteristics of the reservoir, and each small layer corresponds to an independent physical model. These small layers are usually divided according to the results of the stratum combination division given on site, and the well point connectivity between each layer is considered. In each layer, multiple connected units are set based on geological data and well test data. For the flow characteristics in each layer, they are refined by one-dimensional grid discretization. Each grid unit is equivalent to a part of the connected unit, which can more accurately simulate the fluid flow in the area.

[0060] The flow process of the reservoir is converted into a discretized numerical calculation problem. The basic idea of ​​grid discretization is to convert the continuous fluid flow process in each connected unit into a discrete system composed of multiple small units. The properties of each small unit such as permeability, volume and flow can be solved by numerical methods. The connection relationship between these discrete units can be mapped through a standardized model file to generate a grid model suitable for numerical solution. This network model is a hierarchical connected network model.

[0061] For the flow of reservoir matrix, the one-dimensional grid discretization method is used to discretize the well connectivity relationship of each layer into multiple one-dimensional grid units; for the flow in fractures, the fractures are abstracted as a series of characteristic points with high heterogeneous conductivity, and the conductivity of the fractures is characterized by heterogeneous conductivity parameters. The connectivity network between the fracture characteristic points can simulate the high-speed flow of fluid in the fractures and reflect the rapid conduction effect of the fractures on the injection and production fluids; for the flow near the faults, the faults are abstracted as characteristic points with shielding effects, and the characteristics of the fluid flow near the faults are characterized by heterogeneous conductivity parameters. The introduction of fault characteristic points can accurately characterize the blocking effect of the faults on the flow path.

[0062] Through the hierarchical discretization technology, the flow network of matrix, fractures and faults is coupled. The characteristic points of fractures and faults are embedded in the matrix network, and the fluid exchange between the matrix and fractures and faults is realized through conductivity. The fracture characteristic points serve as high-speed channels for matrix flow, and the fault characteristic points play a shielding role for the fluid near the fault, avoiding the poor simulation effect caused by the simplified characterization of fractures and faults in previous methods.

[0063] This grid discretization not only helps improve the accuracy of the simulation, but also facilitates parameter adjustment during the calculation process. The connection between each grid cell and the surrounding cells determines the path and speed of fluid flow, and changes in these connections will directly affect the production dynamics of the entire reservoir.

[0064] By mapping the discretized grid cells into standardized model files, a physically driven hierarchical connected network model can be constructed. This model not only processes the reservoir in layers, but also embeds the fractures and faults existing in the actual complex reservoir, further reducing the resource cost of numerical simulation.

[0065] Step 3, use the ES-MDA algorithm for automatic historical fitting. First, select the daily oil production of a single oil well as the fitting target. On the basis of selecting conductivity, grid volume and phase permeability parameters as fitting parameters, add the heterogeneous conductivity and characteristic points of fractures and faults so that it can adaptively correct the dynamic influence of fractures and faults. Construct a mathematical model of the automatic historical fitting objective function and add engineering constraints. Use the ES-MDA algorithm to iteratively solve the objective function so that the daily oil production simulated by the hierarchical connected network model is consistent with the actual observed daily oil production and conforms to the actual geological understanding. The specific process is: Step 3.1: In reservoir automatic history matching, the construction of the objective function is one of the core links, which determines the effect and efficiency of the fitting process. Traditional history matching methods often rely on manual adjustment of parameters and gradually optimize model parameters through trial and error. In order to improve the automation level and accuracy of fitting, the objective function of history matching is constructed based on the Bayesian Maximum A Posteriori Estimation (MAP), which can effectively integrate physical models and prior knowledge, thereby achieving more accurate history matching.

[0066] In the process of reservoir history matching, the daily oil production of oil wells is usually used as observation data to make the model output as close to the observation data as possible. The constructed history matching objective function is: (15); in, is the historical fitting objective function value, that is, the difference between the model simulation calculation result and the actual observation data; is the model fitting parameter vector; Calculate the results for the model simulation; is the actual observation data; is the covariance matrix of the observation errors.

[0067] By minimizing the historical fitting objective function value, the model calculation results are made consistent with the actual observed data.

[0068] Step 3.2: In the reservoir history fitting process, it is critical to select appropriate control variables. In the reservoir model, the conductivity between discrete grids, the volume of discrete grids, the heterogeneous conductivity and the number of characteristic points of fractures, the heterogeneous conductivity and the number of characteristic points of faults, and the phase permeability parameters are usually selected as control variables to form the model fitting parameter vector; (16); in, is the model fitting parameter vector, which includes the conductivity between all grids and the volume of the grid; The initial moment grid and Grid The conductivity, in units of ; The initial moment grid The volume in units of ; is the characteristic point of the crack at the initial moment and feature points exist The inhomogeneous conductivity of the layer is ; For the initial moment The number of characteristic points of cracks on the layer; is the characteristic point of the fault at the initial moment and feature points exist The inhomogeneous conductivity of the layer is ; For the initial moment The number of characteristic points of the fault on the layer; , , , are different phase permeability parameters; In the process of reservoir history fitting, it is necessary to impose some reasonable physical constraints while minimizing the objective function. This ensures that the model parameters obtained in the fitting process can match the observed data and follow the physical characteristics of the actual reservoir. Considering the reservoir characteristics constraining the range of model parameter changes, specific constraints may include: (17); in, For Grid and Grid conductivity; For Grid Volume; is the number of grids; is the total volume of the reservoir, in .

[0069] Ensure that the inter-well conductivity is always greater than 0, the grid volume is always greater than 0 and less than the sum of all grid volumes, and the sum of all grid volumes is equal to the total reservoir volume. , Always between 1 and 6, , Always between 0 and 1.

[0070] Step 3.3, in the process of complex reservoir development, fractures and faults have a crucial impact on fluid flow characteristics. However, since the spatial distribution of fractures and faults is highly uncertain, and their parameters are difficult to obtain directly through conventional logging or seismic data, it is a major challenge to accurately characterize and model them. In order to overcome the challenge of direct acquisition of fracture and fault parameters, the present invention adopts the ES-MDA (Ensemble Smoother with Multiple Data Assimilation) algorithm, combined with production history data and physical constraints, to achieve automatic inversion and real-time update of fracture and fault parameters.

[0071] The ES-MDA algorithm is a multiple data assimilation technology based on the ensemble method. It can continuously optimize reservoir parameters during the historical matching process, so that the prediction results of the model gradually approach the actual production situation. First, the initial heterogeneous conductivity distribution of fractures and faults is set in the reservoir numerical model, and the number of characteristic points is defined to describe the complex geometric characteristics of fractures and faults. Then, in the historical matching process, the actual production data of the oil well, the daily oil production of a single well, is used as observation information to perform multiple data assimilation on the model. Through the iterative update mechanism of the ES-MDA algorithm, the matrix parameters and the parameters of fractures and faults are continuously adjusted, so that the fitting degree of the simulation results with the actual production data is gradually improved.

[0072] Through the deep integration of data-driven and physical mechanisms, this method enables the reservoir model to adaptively adjust the dynamic impact of fractures and faults on fluid flow, thereby improving the model's ability to predict the injection and production response of complex reservoirs. Compared with traditional manual adjustment methods, the automatic inversion and real-time update mechanism of the ES-MDA algorithm not only improves the computational efficiency, but also can more accurately capture the changing characteristics of the internal structure of the reservoir, providing more reliable technical support for dynamic reservoir management and optimized development.

[0073] Step 4: Establish a mathematical model of the objective function for real-time optimization of stratified injection and production, and perform real-time optimization of stratified injection and production based on the SHADE algorithm. First, select the economic net present value as the real-time optimization target for injection and production, and the injection and production volume of a single well layer as the optimization variable. Considering the actual engineering constraints, establish a mathematical model for the objective function for real-time optimization of stratified injection and production, introduce an adaptive mechanism based on the success history and the SHADE algorithm with an archive mechanism, which significantly improves the solution efficiency and accuracy of the injection and production optimization model, and maximizes the economic benefits by iteratively maximizing the optimization target. The specific process is: Step 4.1, real-time optimization of injection and production is to maximize the recovery rate and economic benefits of the oil field by continuously monitoring the production status of the oil field and adjusting the water injection and oil production strategies. It can also make the oil field management more flexible and intelligent, and provide strong support for the long-term stable development of the oil field. This process involves dynamically adjusting the water injection volume, liquid production and the working status of the well, so as to optimize resource allocation and improve the overall economic benefits while ensuring the long-term stable production of the oil field. In the actual oil field production process, as the mining process advances, the geological, physical properties and fluid distribution of the reservoir will change, and traditional injection and production optimization methods are difficult to cope with these dynamic changes. Therefore, real-time optimization of injection and production has become an effective means to meet this challenge.

[0074] To achieve the above goals, it is reasonable to select the economic net present value (NPV) as the objective function of real-time optimization of injection and production. The economic net present value is obtained by discounting future cash flows to the current moment, taking into account the return on investment and production costs, and then obtaining the economic benefits of the oil field project. Using the economic net present value as the objective function can effectively balance the relationship between the improvement of recovery rate and economic benefits, ensuring that the economic return is maximized while pursuing long-term stable production. Through real-time optimization, the economic net present value can be used as a direct indicator to measure the optimization effect, providing a scientific decision-making basis for the sustainable development of the oil field. The following is the mathematical model of the optimization objective function established: (18); in, is the economic net present value, in yuan; is the controlled variable, representing the regulation scheme of stratified water injection volume and stratified liquid production volume during the optimization process; is the oil price, in yuan; For the Time step well The oil production in ; is the sewage treatment cost, in yuan; For the Time step well The water output, in units of ; is the water injection cost, in yuan; For the Time step well The water injection volume, in units of ; is the total number of time steps; is the number of producing wells; is the number of water injection wells; Step 4.2: During reservoir development, the complexity and variability of reservoirs make the formulation and optimization of development plans challenging. In order to ensure the rationality and feasibility of the development plan, it is necessary to reasonably constrain the relevant parameters according to the specific conditions of the reservoir in practical applications. These constraints can not only improve the accuracy of reservoir numerical simulation and optimization processes, but also avoid unrealistic results in actual operations. The following are some common reservoir development parameter constraints, which help improve production stability and economic benefits while ensuring the effect of reservoir development.

[0075] A reasonable injection-production ratio can ensure that water injection and oil production activities in the reservoir are coordinated with each other to avoid excessive water injection or excessive oil production. Too high a water injection ratio may lead to flooding of the reservoir and affect the long-term production potential of the reservoir; too low a water injection ratio may cause insufficient displacement of the reservoir and affect the recovery rate.

[0076] Bottom hole pressure directly affects the production capacity of the reservoir and the working status of the well. Reasonable control of bottom hole pressure can avoid wellbore collapse, downhole equipment damage and other problems, and can also ensure the stable exploitation of the reservoir. Too high bottom hole pressure may cause excessive wellbore pressure and affect the oil production effect; too low bottom hole pressure may cause the well's oil pumping capacity to decrease and fail to achieve the ideal production.

[0077] As production progresses, the proportion of water tends to gradually increase, which is also a common problem in the process of reservoir development. Too high a water content means that more water is produced, which makes oil-water separation difficult, reduces oil quality, and increases post-processing costs. Therefore, a maximum water content constraint is usually set in reservoir development to prevent the water content from exceeding a certain critical value, thereby ensuring the economic efficiency of the production process. This constraint is usually set based on the economic analysis of the oil field, equipment capacity, and production requirements.

[0078] The single well maximum liquid volume constraint is set based on the well's pumping capacity, the wellbore's fluid delivery capacity, and the wellhead facility's processing capacity. Too much liquid volume may exceed the wellbore's carrying capacity, leading to equipment failure or reduced production efficiency. At the same time, too little liquid volume may also lead to underdevelopment of the reservoir and failure to achieve optimal recovery results.

[0079] The main constraints are injection-production ratio constraint, bottom hole pressure constraint, maximum water cut constraint, and single well maximum liquid volume constraint. The specific constraints are as follows: (19); In the formula, is the minimum injection-production ratio; is the maximum injection-production ratio; For the The total water injection volume in the time step is ; For the Time step well The amount of liquid produced or water injected; For the The total liquid production in time steps, in units of ; For the Time step well The minimum bottom hole flowing pressure is in ; For the Time step well The maximum bottom hole flowing pressure, in ; For the Time step well Bottom hole pressure; For the Time step well The moisture content; For the Time step well The ultimate moisture content; For the Time step well The liquid production or water injection volume, in units of ; For the Time step well The maximum liquid production or maximum water injection volume, in units of .

[0080] Step 4.3, the traditional DE algorithm needs to manually set the scaling factor (F) and crossover probability (CR) during the optimization process, and its fixed parameters are difficult to adapt to the dynamic characteristics of complex injection-production optimization problems. In view of the limitations of the traditional differential evolution algorithm (DE) that it is easy to fall into local optimality and slow convergence when dealing with high-dimensional and nonlinear optimization problems, the SHADE algorithm (Success-History based Adaptive DifferentialEvolution) with an adaptive mechanism based on success history and an archive mechanism is introduced. At the same time, the engineering constraints are considered to solve the optimization objective function to obtain the maximum economic net present value, which significantly improves the solution efficiency and accuracy of the injection-production optimization model, helps to timely adjust the injection-production system according to different stages and real-time production data during oilfield development, solves problems such as inter-layer interference and well connectivity, and improves production efficiency.

[0081] The SHADE algorithm is one of the modern evolutionary algorithms. Its core idea is to adjust the search strategy by introducing an adaptive mechanism, thereby effectively improving the convergence speed and global search ability of the optimization algorithm. SHADE dynamically adjusts the mutation parameters by using historical successful experience (successful mutation operations) in each generation of the algorithm, which enables SHADE to show high robustness in complex optimization problems. The SHADE algorithm dynamically adjusts the values ​​of F and CR by introducing a historical successful parameter adaptive mechanism, so that it can automatically adapt to the problem characteristics according to the feedback information of the optimization process. This adaptive mechanism not only reduces the dependence on manual parameter adjustment, but also significantly improves the robustness and convergence efficiency of the algorithm in injection and production optimization. The SHADE algorithm also introduces a historical memory library (Memory Archive) to store the parameter information of historical successful individuals and generate new trial individuals based on this information. This mechanism enables the algorithm to make full use of historical optimization experience, avoid invalid search, and thus accelerate convergence. This feature is particularly important in injection and production optimization because it can quickly locate high potential solution areas and reduce the waste of computing resources.

[0082] The SHADE algorithm includes mutation operations, crossover operations, and selection operations. In reservoir injection and production optimization, mutation operations can generate new solutions by adjusting the injection and production flow rates, while crossover operations can combine different injection and production schemes to obtain more possible solutions. In this way, the SHADE algorithm can generate new solutions that meet the optimization requirements while ensuring the constraints.

[0083] The individuals of each generation will be evaluated for fitness and sorted according to the objective function value and the satisfaction of constraints. The fitness evaluation not only considers the objective function (maximizing the economic net present value), but also ensures that the engineering constraints (such as injection-production ratio, bottom hole pressure, maximum water content, etc.) are met. Those solutions that meet the engineering constraints and maximize the objective function will be given priority.

[0084] According to the fitness evaluation results, the selection operation will select individuals with higher fitness to enter the next generation. These individuals will undergo mutation and crossover operations according to the adaptive mechanism to generate new solutions. In reservoir injection and production optimization, the selection operation can select the best injection and production scheme by evaluating the objective function and constraints, thereby guiding the actual operation in oilfield development.

[0085] In order to demonstrate the feasibility and superiority of the present invention, the following specific examples are given.

[0086] In the embodiment of the present invention, a double-layer one-injection and four-production model containing fractures and faults is established using a finite volume method reservoir numerical simulator, with a grid scale of 25*25*2, a model size of 100m*100m*20m, a reservoir porosity of 0.3, an initial water saturation of the reservoir of 0.25, oil and water viscosities of 4mPa*s and 1mPa*s, respectively. There is 1 injection well and 4 production wells in the reservoir, the injection well number is 5, and the production well numbers are 1, 2, 3, and 4. The reservoir numerical simulator is used to simulate production for 4770 days, and the overall injection and production of the block is balanced.

[0087] Then as Figure 2 The method of the present invention is shown to predict the injection and production volume of a small layer. The specific process is as follows: first, the reservoir data is collected to construct a sample library for vertical splitting of the injection and production volume, the production and absorption profile in the data is used as the label data, and the porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability extreme difference, porosity extreme difference, bottom hole flow pressure, single well injection and production volume and water content are used as input features; then, a Gaussian mixture model is constructed to output the probability distribution of the injection and production volume of the small layer; the reservoir data is integrated with the probability distribution of the injection and production volume of the small layer to construct a GMM-XGB hybrid model, and the Bayesian optimization algorithm is used to train the model, and finally the predicted injection and production volume of the small layer is output. The embodiment of the present invention splits the injection and production volume of the whole well into two small layers based on the above process.

[0088] Then, the same reservoir parameters are used to establish a hierarchical connected network model, and the reservoir is equivalent to an inter-well connected network formed by a series of well points and hierarchical inter-well connecting pipes. The well points are connected to each other through connecting pipes to form an inter-well connected network. The complexity of this network can be flexibly adjusted by controlling the well spacing and well angle. In this network, the inter-well pipe is regarded as a connected unit, which is characterized by two parameters: inter-well conductivity and connected volume. Then, the fractures and faults are abstracted as discrete feature points and embedded in the inter-well connected network to reconstruct the connectivity relationship between the well points and the fracture and fault nodes. The specific processing method is as follows: Figure 3As shown in the figure, five well nodes are first used to represent five well points, the gray diamond structure is a fault, and the gray rectangular structure is a fracture. Then the fractures and faults are abstracted into three fracture nodes and six fault nodes respectively. Finally, the connectivity relationship between the five well nodes (numbered 1 to 5) and the fracture nodes and fault nodes is reconstructed to obtain the required hierarchical connectivity network model. Finally, the coupling of the matrix-fracture-fault multi-scale flow network is realized by hierarchical discretization technology: for the flow of the matrix, one-dimensional grid discretization is used to describe the conventional seepage between wells; for the flow of fractures and faults, the conductivity of the fracture and the shielding characteristics of the fault are characterized by the heterogeneous conductivity parameters between the feature points.

[0089] The model's single forward run time is 32 seconds, which is about 4 times faster than the finite volume reservoir numerical simulator, which takes 123.3 seconds. The ES-MDA algorithm is used for historical matching, with a population of 200 and three iterations. Considering the constraints, the output of the prior model calculation results of all small layers and the actual observation data comparison image during the historical matching period is shown. Here, only the prior model calculation results of the production well 1 small layer 1 and the actual observation data comparison image are shown. Figure 4-Figure 6 These are the historical fitting results for the 1st, 2nd and 3rd iterations respectively. It can be seen that during the historical fitting period, the daily oil production curve obtained by each set of fitting parameters (corresponding to the calculation results of the prior model) envelops the daily oil production curve obtained by the numerical simulator (corresponding to the actual observed data). Figure 4-Figure 6 In the figure, one population corresponds to one set of fitting parameters, so the figure contains the daily oil production curves obtained by 200 sets of fitting parameters corresponding to 200 populations, and some gray dotted lines overlap.

[0090] After historical matching, the output block moisture content fitting curve is as follows Figure 7 As shown, it can be seen that the true value is basically consistent with the fitting value output by the numerical simulator after history matching.

[0091] After historical matching, the fitting curve of production well 1 layer 1 is output as follows Figure 8 As shown, it can be seen that the daily oil production rate of production well 1 using the method of the present invention is basically consistent with the daily oil production rate output by the numerical simulator.

[0092] The model after historical matching of the present invention is used to optimize the layered injection and production by using the SHADE algorithm combined with engineering constraints. A total of 12 time steps are optimized, each time step is 30 days, and the total optimization time is 360 days. Fig. 9 and Fig.10 As shown, only the injection system control diagram of sublayer 1 and sublayer 2 of injection well 5 is shown. Fig. 9 and Fig.10It can be seen that, subject to the 20% upper and lower constraints of the single well liquid volume of the injection well, the optimized water injection rate is constrained to be between 0.8 and 1.2 times the water injection rate before optimization.

[0093] After adjusting the injection and production system for all layers of all wells, the overall optimization results are as follows: Fig.11 and Fig.12 As shown, Fig.11 and Fig.12 They are the optimization results of cumulative oil production and water content, respectively. It can be seen that the cumulative oil production after optimization is significantly improved, and the water content is slightly reduced, which verifies the reliability of the method of the present invention in complex reservoir simulation and optimization calculation.

[0094] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs, characterized in that: The steps include: Step 1: Construct a GMM-XGB mixture model based on the Gaussian mixture model and XGBoost algorithm; Step 2: Construct a hierarchical connected network model based on the physical model; Step 3: Use ES-MDA algorithm for automatic history matching; Step 4: Establish a mathematical model of the objective function for real-time optimization of stratified injection and production, and perform real-time optimization of stratified injection and production based on the SHADE algorithm.

2. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 1, characterized in that: In step 1, the GMM-XGB hybrid model uses a Gaussian mixture model to perform cluster analysis on the original sample set to obtain the probability distribution of the small layer injection and production volume, and integrates the probability distribution into the original sample set as the input of the XGBoost algorithm to predict the output of the small layer injection and production volume; the specific working process of the GMM-XGB hybrid model is as follows: Step 1.1, collect the longitudinal splitting data of injection and production volume, and construct the original sample set; The vertical splitting data of injection and production volume include porosity, permeability, oil layer thickness, initial water saturation, formation coefficient, number of perforations, perforation thickness, permeability extreme difference, porosity extreme difference, bottom hole flowing pressure, single well injection and production volume and water content; Step 1.2: Preprocess the original sample set. The specific process is as follows: First, data cleaning is performed. For missing values ​​of random missing types, the natural neighbor interpolation method in MATLAB software is used for interpolation according to the geodetic coordinates of the well points. For outliers, deletion or transformation processing is selected according to the actual situation. Transformation processing includes smoothing or scaling. For categorical data, the number of base classes of categorical data needs to be determined first. If there are few base classes, one-hot encoding is used for processing; if there are many base classes, label encoding or target encoding is used for processing. Then, the cleaned data was normalized using the Min-Max normalization method; Finally, the normalized sample set was subjected to correlation analysis; for linear relationships in the data, the Pearson correlation coefficient was used for correlation analysis; for sample sets with strong nonlinear relationships and no normal distribution, the Spearman rank correlation coefficient was used for correlation analysis; Step 1.3, using Gaussian mixture model to perform cluster analysis on the preprocessed original sample set to obtain the probability distribution of small layer injection and production volume; Step 1.4, the probability distribution of the small layer injection and production volume is fused with the preprocessed original sample set to obtain a fused sample set, and the fused sample set is divided into a training set and a validation set; Step 1.5: Input the fused sample set into the XGBoost algorithm and use the Bayesian optimization method to train until the requirements are met, and finally output the injection and production volume of all small layers in all wells and all time steps.

3. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 2, characterized in that: The specific process of step 1.3 is as follows: Step 1.3.1: Based on the preprocessed original sample set, a Gaussian mixture model is constructed; the Gaussian mixture model assumes that the injection and production volume data points come from multiple different Gaussian distributions, and the probability density function of the Gaussian distribution is: ； in, is the probability density function of Gaussian distribution; is the data point of injection and production volume; is the mean vector; is the covariance matrix; is the data variance; is the data dimension; is an exponential function with base e; Step 1.3.2: Use the expectation maximization algorithm to iteratively solve the established Gaussian mixture model until convergence. The specific process is as follows: First, initialize the parameters of each Gaussian distribution in the Gaussian mixture model; Then, the responsiveness of all injection and production volume data points is calculated. The responsiveness calculation formula is as follows: ； in, is responsiveness; is a binary indicator variable, indicating Injection and production volume data points Whether Gaussian components are generated; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; is the index variable of the Gaussian component, which is used to traverse all Gaussian components in the mixed Gaussian model for summation; is the total number of Gaussian components in the mixed Gaussian model; For the The prior probability of the Gaussian components; For the The mean vector of the Gaussian components; For the The covariance matrix of the Gaussian components; The parameters of each Gaussian distribution are updated according to the calculated responsiveness. The update formula is as follows: ； ； ； in, , , Respectively The iteration The mean vector, covariance matrix, and prior probability of the Gaussian components; For the Responsiveness at iterations; is the total number of data points of injection and production volume; is the transpose symbol; Calculate the likelihood function of the model, the formula is as follows: ； in, is the likelihood function; Finally, the model parameters are continuously updated iteratively until the parameters of the Gaussian mixture model converge. At this time, the training ends and the final responsiveness is output. The responsiveness is the probability distribution of the injection and production volume of the small layer.

4. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 3, characterized in that: In step 1.5, Bayesian optimization is an automated hyperparameter tuning method that establishes a proxy model and continuously adjusts hyperparameters based on the historical performance of the model to find the optimal hyperparameter configuration. The specific process is as follows: Step 1.5.1, initialize the Bayesian optimization algorithm framework, initialize the evaluation of the objective function in the framework and set the hyperparameter search space of the GMM-XGB hybrid model; the objective function uses the mean absolute error, and the formula is as follows: ； in, is the mean absolute error; For the The true value of samples; For the The predicted value of samples; is the total number of samples; Hyperparameters include learning rate, maximum tree depth, subsample ratio, and number of tree structures. The search space of learning rate is set to [0.01, 0.2], the search space of maximum tree depth is set to [5, 10], the search space of subsample ratio is set to [0.5, 1], and the search space of number of tree structures is set to [100, 500]. Step 1.5.2, iterative training, to find the optimal hyperparameters of the model in the search space; firstly, sample collection is performed, that is, a set of initial points are selected as the initial hyperparameters in the model hyperparameter combination space to be optimized, and then the Bayesian optimization algorithm evaluates the next set of hyperparameters to be selected based on the existing hyperparameters and objective function values; finally, the hyperparameter combination with the largest objective function is selected as the optimal hyperparameter combination; Step 1.5.3: Use the selected optimal hyperparameter combination to train the GMM-XGB hybrid model on the training set, calculate the mean absolute error between the predicted value and the true value of the GMM-XGB hybrid model on the validation set, and obtain the objective function value; Step 1.5.4: After obtaining the new hyperparameter-objective function pair, update the surrogate model. The surrogate model will continuously adjust its estimate of the distribution of the objective function in the hyperparameter space based on the new data. Step 1.5.5: Determine whether the pre-set termination condition is met; if the termination condition is not met, return to step 1.5.2 to collect samples again and continue iterating; when the termination condition is met, select a set of hyperparameters that optimize the objective function value from all evaluated hyperparameter combinations, and use this set of optimal hyperparameters to retrain the GMM-XGB hybrid model on the training set until the pre-set number of training times is reached; after the training is completed, the final GMM-XGB hybrid model for prediction is obtained.

5. The rapid simulation and real-time optimization method for layered injection and production applicable to complex oil reservoirs according to claim 4, characterized in that: The specific process of step 2 is: Step 2.1, the oil reservoir is divided into several well points and connecting pipes, and an inter-well connecting network is constructed; in the inter-well connecting network, each well point represents a water injection well or oil production well on a small layer, and the well points are connected to each other through connecting pipes; the connecting pipes are characterized by two parameters: inter-well conductivity and connecting volume; The conductivity calculation formula is as follows: ； in, Well and well The conductivity between Well and well The seepage area between Well and well The average value of the permeability between Well and well The distance between The formula for calculating the connected volume is as follows: ； in, Well and well The connected volume between them; Well and well The average value of the effective thickness between Step 2.2, discretize the fracture into a group of uniformly arranged high permeability nodes; discretize the fault into two groups of uniformly arranged nodes, located on both sides of the fault, and characterize the flow capacity of the fluid near the fault by heterogeneous conductivity; embed the fracture and fault nodes into the well-to-well connection network to construct a dynamic flow channel; Step 2.3: Based on the hierarchical well connectivity relationship, multiple connectivity units are set in each layer and discretized into one-dimensional grids. Then, the grid connectivity relationship is mapped into a standardized model file to construct a physically driven hierarchical connectivity network model.

6. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 5, characterized in that: The specific process of step 2.3 is as follows: First, the reservoir is divided into multiple small layers according to its geological characteristics, and each small layer corresponds to an independent physical model. In each layer, multiple connected units are set based on geological data and well test data. The flow characteristics in each layer are refined by one-dimensional grid discretization, and each grid unit is equivalent to a part of the connected unit. The flow process of the reservoir is converted into a discretized numerical calculation problem. The basic idea of ​​grid discretization is to convert the continuous fluid flow process in each connected unit into a discrete system composed of multiple small units. The permeability, volume and flow rate of each small unit are solved by numerical methods. The connection relationships between discrete units are mapped through standardized model files to generate a hierarchical connected network model.

7. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 6, characterized in that: The specific process of step 3 is as follows: Step 3.1: Construct the history matching objective function as: ； in, is the historical fitting objective function value; is the model fitting parameter vector; Calculate the results for the model simulation; is the actual observation data; is the covariance matrix of the observation error; Step 3.2, select the conductivity between discrete grids, the volume of discrete grids, the heterogeneous conductivity and the number of characteristic points of fractures, the heterogeneous conductivity and the number of characteristic points of faults, and the phase permeability parameter as control variables to form a model fitting parameter vector; ； in, The initial moment grid and Grid The conductivity between The initial moment grid Volume; is the characteristic point of the crack at the initial moment and feature points exist Heterogeneous conductivity across the layer; For the initial moment The number of characteristic points of cracks on the layer; is the characteristic point of the fault at the initial moment and feature points exist Heterogeneous conductivity across the layer; For the initial moment The number of characteristic points of the fault on the layer; , , , are different phase permeability parameters; In the process of reservoir history matching, physical constraints are imposed while minimizing the objective function, including: ； in, For Grid and Grid conductivity; For Grid Volume; is the number of grids; is the total volume of the reservoir; Step 3.3: Use the ES-MDA algorithm, combined with historical production data and physical constraints, to achieve automatic inversion and real-time update of fracture and fault parameters. The specific process is as follows: Firstly, the initial heterogeneous conductivity distribution of fractures and faults is set in the reservoir numerical model, and the number of characteristic points is defined. Then, in the history fitting process, the actual production data of oil wells (daily oil production per well) is used as observation information to perform multiple data assimilation on the model. Through the iterative update mechanism of the ES-MDA algorithm, the matrix parameters as well as the parameters of fractures and faults are continuously adjusted.

8. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 7, characterized in that: The specific process of step 4 is as follows: Step 4.1: Select the economic net present value as the objective function of real-time optimization of injection and production; the following is the mathematical model of the optimization objective function: ； in, is the economic net present value; is the controlled variable, representing the regulation scheme of stratified water injection volume and stratified liquid production volume during the optimization process; for oil prices; For the Time step well Oil production; Cost of sewage treatment; For the Time step well Water production; is the water injection cost; For the Time step well The amount of water injected; is the total number of time steps; is the number of producing wells; is the number of water injection wells; Step 4.2: Set engineering constraints, including injection-production ratio constraint, bottom hole pressure constraint, maximum water cut constraint, and single well maximum liquid volume constraint. The specific constraints are as follows: ； In the formula, is the minimum injection-production ratio; is the maximum injection-production ratio; For the Time step well The amount of liquid produced or water injected; For the Time step well Minimum bottom hole flowing pressure; For the Time step well Maximum bottom hole flowing pressure; For the Time step well Bottom hole pressure; For the Time step well The moisture content; For the Time step well The ultimate moisture content; For the Time step well The amount of liquid produced or water injected; For the Time step well The maximum liquid production or maximum water injection volume; Step 4.3: Introduce the adaptive mechanism based on success history and the SHADE algorithm with archive mechanism to optimize reservoir injection and production.

9. The method for rapid simulation and real-time optimization of stratified injection and production applicable to complex oil reservoirs according to claim 8, characterized in that: In step 4.3, the specific process of reservoir injection and production optimization by SHADE algorithm is as follows: Step 4.3.1, adjusting the injection and production flow rate to generate a new injection and production plan; Step 4.3.2, combine different injection and production schemes and calculate the economic net present value of each injection and production scheme combination; Step 4.3.3: Select the injection-production scheme combination that meets the engineering constraints and has the largest economic net present value as the individual to enter the next generation. The individual will perform steps 4.3.1 and 4.3.2 again according to the adaptive mechanism until the number of iterations is reached; after the iteration, the final injection-production scheme combination is the optimal injection-production scheme.

Images