A rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints

Abstract

This application relates to the technical field of reservoir numerical simulation, and discloses a rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints. The sedimentary facies boundaries are identified according to the distribution characteristics of sedimentary facies, and a connectivity network model considering sedimentary facies constraints is established based on the sedimentary facies boundaries, well positions, and sedimentary facies attributes. The connectivity network model is divided into one-dimensional grids and converted into the grid connection format of a general simulator, and the general simulator is used for solving to obtain pressure and saturation solutions. An automatic history matching mathematical model is established based on sedimentary facies and reservoir constraint conditions. Taking the economic net present value as an index and the injection volume, production liquid volume, and polymer injection volume of each well as optimization variables, a reservoir injection-production optimization mathematical model is established by using the differential evolution algorithm and constraining the total polymer injection volume and reservoir pressure. This application comprehensively considers the characteristics of the reservoir after waterflooding to polymer flooding to optimize the reservoir, so as to achieve rapid decision-making for oilfields and guide oilfield construction.

Description

Technical Field

[0001] This application relates to the technical field of reservoir numerical simulation, and particularly to a rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints. Background Art

[0002] Oilfields that have undergone years of waterflooding development generally face problems such as high water cut, high decline rate, and low recovery rate, and the water injection efficiency of waterflooding development has decreased significantly. To improve the oilfield development benefit, various enhanced oil recovery (EOR) methods have been developed currently. Among them, polymer flooding technology has become one of the important EOR means because it can effectively improve the migration characteristics of reservoir fluids and significantly increase the recovery rate. The core principle of polymer injection is to improve the viscosity of the injected water and reduce the seepage capacity of the water phase during the displacement process, thereby improving the swept volume and enhancing the oil displacement efficiency. Compared with traditional waterflooding, polymer flooding can effectively reduce water channeling, lower the water cut, and significantly increase the crude oil recovery rate. However, due to the high cost of polymers, how to optimize the injection-production plan and reduce the development cost while ensuring the increase in recovery rate has become an urgent problem to be solved in the current oilfield development process. On the other hand, with the development of intelligent oilfield technology, production optimization has become an important direction of oilfield management. The core idea of production optimization is to construct an optimization objective function and reasonably adjust control variables such as well control, so that the oilfield can maximize the crude oil production and economic benefits while reducing the development cost. In the process of converting from waterflooding to polymer flooding, how to achieve efficient numerical simulation and intelligent optimization decision-making to more scientifically guide the polymer injection strategy is the key to improving the reservoir development effect. Summary of the Invention

[0003] The purpose of this application is to provide a rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints. By building a connected network model with sedimentary facies constraints, numerical simulation solution, automatic history matching, and intelligent production optimization, an efficient optimization framework for converting from waterflooding to polymer flooding is constructed. This method can quickly predict the polymer flooding effect, optimize the injection-production strategy, and reduce the development cost, providing scientific support for the efficient development of oilfields.

[0004] Generally speaking, first, according to the sedimentary facies distribution of the reservoir, an image recognition method combining edge detection and level set function is used to identify the boundaries between sedimentary facies, generating a binary map. Then, based on the actual geological information of each well point, an initial well pattern subdivision map is generated according to the angle and distance discrimination criteria. Wells are connected by a one-dimensional connected unit, and grids are divided within the one-dimensional connected unit. On this basis, combined with the binary map, parameters such as the permeability of each phase are assigned within each grid to form an inter-well connected network model. The constructed model is mapped into a 2D simulation grid and coupled with a commercial numerical simulator (such as ECLIPSE), combining simple modeling with a fast calculation simulator. The historical fitting algorithm ES-MDA is used to update parameters such as the model grid size and permeability, and constraint conditions are added to fit the production data on the premise that the model conforms to geological characteristics. The fitted model is regarded as a real reservoir model. Considering the characteristics of polymer flooding after water flooding, the particle swarm algorithm is used, with the economic net present value as the objective function, to obtain a development strategy that can maximize the development benefits of the oilfield while constraining the parameters of each well and block, and a production plan is formulated to guide the subsequent production process of the oilfield.

[0005] To achieve the above objectives, the following technical solutions are adopted:

[0006] This application provides a rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary facies constraints. The method includes:

[0007] Identify sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and based on the sedimentary facies boundaries, well positions, and attributes of each sedimentary facies, establish a connected network model considering sedimentary facies constraints;

[0008] Perform one-dimensional grid division on the connected network model, convert it into the grid connection format of a general simulator, and use the general simulator for solution to obtain pressure and saturation solutions;

[0009] Based on sedimentary facies and reservoir constraint conditions, adopt the integrated smoothing multi-data assimilation algorithm to establish an automatic history matching mathematical model to realize the update and inversion of key reservoir parameters; among them, the key reservoir parameters include the permeability field and relative permeability curve;

[0010] Taking the economic net present value as an index, with the injection volume, production fluid volume, and polymer injection volume of each well as optimization variables, use the differential evolution algorithm to establish an oil reservoir injection-production optimization mathematical model by constraining the total polymer injection volume and reservoir pressure.

[0011] Furthermore, identifying sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and based on the sedimentary facies boundaries, well positions, and attributes of each sedimentary facies, establishing a connected network model considering sedimentary facies constraints includes:

[0012] According to the sedimentary facies distribution of the block, use the edge detection method to extract the facies boundaries in the sedimentary facies distribution map;

[0013] Take the facies boundaries in the sedimentary facies distribution map as the initialization boundaries of the level set function, and perform regional evolution through the level set method to obtain a binary image;

[0014] According to the actual geological information of each well point, generate an initial connected network model according to the angle and distance discrimination criteria, and connect the wells with a one-dimensional connected unit.

[0015] Furthermore, according to the sedimentary facies distribution of the block, using the edge detection method to extract the facies boundaries in the sedimentary facies distribution map includes:

[0016] Convert the obtained sedimentary facies map into a grayscale image, and use the Gaussian function to perform weighted averaging according to the gray values of the pixel points to be filtered and their neighborhoods to filter out the superimposed high-frequency noise in the image. The Gaussian function is as follows:

[0017] (1)

[0018] In the formula, is the standard deviation controlling the Gaussian filter; f ( x , y ) is the Gaussian function; e is the natural constant;

[0019] Use the Sobel operator to calculate the magnitude and direction of the image gradient through the following formulas (2), (3) and (4):

[0020] (2)

[0021] (3)

[0022] (4)

[0023] In the formula, are the gradients in the x and y directions of the image respectively, is the original image, is the Sobel filter, calculate the gradients in the horizontal and vertical directions respectively, and the total gradient G of each pixel point is calculated by the formula:

[0024] (5)

[0025] The expression of the gradient direction is:

[0026] (6)

[0027] For each pixel, check whether its gradient magnitude is a local maximum along the gradient direction, accurately locate the edges, connect the scattered edge points into a complete boundary, obtain a binary image, and use the binary image as the phase boundary in the sedimentary facies distribution map.

[0028] Further, use the phase boundary in the sedimentary facies distribution map as the initial boundary of the level set function, and perform regional evolution through the level set method to obtain a binary map, including:

[0029] Select the signed Euclidean distance from a point on the plane to the contour curve as the level set function :

[0030] (7)

[0031] where represents a point on the image to the level set contour curve distance;

[0032] Express the level set function that changes with time as , and the evolution of the surface is described by the following level set evolution equation:

[0033] (8)

[0034] where is the velocity function, which determines the evolution direction and velocity of the surface, represents the gradient magnitude of the level set function, t represents the actual time;

[0035] Before evolving the level set, convert the level set function into an SDF function through the following formula:

[0036] (9)

[0037] (10)

[0038] where is the re-initialization virtual time variable, independent of the actual time t, is the initial level set function, is the sign function, used to keep the positive and negative signs unchanged, is a constant, used to avoid the denominator being zero;

[0039] Perform level set evolution based on the level set evolution equation to obtain a binary map, which is used to divide the regions corresponding to multiple sedimentary facies.

[0040] Furthermore, according to the actual geological information of each well point, an initial connected network model is generated according to the angle and distance discrimination criteria, and the wells are connected by a one-dimensional connected unit, including:

[0041] According to the location of the well points in the reservoir, each well is connected to form an initial well-to-well fully connected network;

[0042] When the maximum angle of the triangle formed by the connection is greater than the set angle, it is determined that the two farther wells are indirectly connected through the middle well, and the longest connection is deleted;

[0043] When the distance between two wells is greater than the set distance, it is determined that there will be no direct impact effect between the two wells, and the corresponding connection is deleted.

[0044] Furthermore, the connected network model is divided into one-dimensional grids and converted into a grid connection format of a general simulator, and solved using a general simulator to obtain pressure and saturation solutions, including:

[0045] Dividing a one-dimensional grid in each connected unit in the connected network model, combining the binary graph with the connected network model, and assigning values ​​to each grid parameter in the one-dimensional grid;

[0046] The one-dimensional connection unit of the connected network model is converted into the grid connection format of the general simulator, and the general simulator is used to solve the problem to obtain the pressure and saturation solutions.

[0047] Furthermore, a one-dimensional grid is divided in each connected unit in the connected network model, and the binary graph is combined with the connected network model to assign values ​​to each grid parameter in the one-dimensional grid, including:

[0048] Based on the initial well pattern segmentation diagram, each connected unit is subdivided into a series of one-dimensional grids, each of which is characterized by grid parameters; wherein the grid parameters include permeability, porosity and / or water saturation;

[0049] According to the binary map, the one-dimensional grids corresponding to the different sedimentary phases are assigned values ​​respectively, so that the connected network model is combined with the sedimentary phase to form a connected network model based on the sedimentary phase.

[0050] Furthermore, the one-dimensional connection unit of the connected network model is converted into the grid connection format of the general simulator, and the general simulator is used to solve the pressure and saturation solutions, including:

[0051] Map the connected network model to a 2D Cartesian coordinate system. Each row represents a connected unit, and the total number of rows is the total number of connected units. Each column represents the grid corresponding to the positions of different connection units, and the total number of columns is equal to the number of grids divided by each connected unit. Represent the well point grids corresponding to non-adjacent connections with dead grids;

[0052] Bring the connected network model mapped to the 2D Cartesian coordinate system into a general simulator to quickly obtain pressure and saturation solutions.

[0053] Furthermore, based on sedimentary facies and reservoir constraint conditions, an integrated smooth multi-data assimilation algorithm is adopted to establish an automatic history matching mathematical model to realize the updated inversion of key reservoir parameters, including:

[0054] Construct a history matching objective function and determine constraint conditions in combination with sedimentary facies and reservoirs;

[0055] Based on the integrated smooth multi-data assimilation algorithm, establish an automatic history matching mathematical model to realize the two-step parameter inversion process.

[0056] Furthermore, for oil-water two-phase flow, the oil-water relative permeability curve is described by the following formula:

[0057] (11)

[0058] In the formula, and respectively represent the endpoints of the oil-water relative permeability curve, and respectively represent the exponents of the oil-water relative permeability curve, used to describe the curvature of the curve, and respectively represent the residual oil saturation and the irreducible water saturation, and respectively represent the oil saturation and the water saturation, represents the relative permeability of the oil phase, represents the relative permeability of the water phase.

[0059] Furthermore, construct a history matching objective function and determine constraint conditions in combination with sedimentary facies and reservoirs, including:

[0060] Construct a history matching objective function, expressed as:

[0061] (12)

[0062] In the formula, are adjustable parameters of the model, is the objective function; is the simulation result of the connected network model; represents the actual observed data; represents the covariance matrix of the observation error; min represents taking the minimum value; T is the symbol for matrix transpose operation;

[0063] The adjustable parameters of the model are expressed as:

[0064] (13)

[0065] In the formula, is the grid volume, is the grid permeability, is the depth of the oil-water contact surface;

[0066] The grid properties and reservoir conditions in each phase are constrained by the following formula:

[0067] (14)

[0068] In the formula, is the grid permeability of the sedimentary facies; is the grid permeability of the sedimentary facies; and are the minimum and maximum permeabilities of the sedimentary facies respectively; and are the minimum and maximum permeabilities of the sedimentary facies respectively; and are the grid volumes of the sedimentary facies and sedimentary facies respectively; and are the total reservoir volumes of the sedimentary facies and sedimentary facies respectively; and respectively represent the grid numbers of the sedimentary facies and and represent the lowest and highest depths of the oil-water contact surface respectively.

[0069] Furthermore, based on the integrated smooth multi-data assimilation algorithm, an automatic history matching mathematical model is established to realize the two-step parameter inversion process, including:

[0070] For a model m, there are parameters of dimension, observation data and A set, in the nth assimilation step, the update of the parameter in ensemble member j is as follows:

[0071] (15)

[0072] where is the covariance matrix between the parameter and the predicted data vector in the nth assimilation step, is the autocovariance matrix of the predicted data in the nth assimilation step, is the covariance matrix of the measurement error, is the inflation coefficient, represents the perturbed observed data, is the predicted data, is the parameter to be updated in ensemble member j in the (n + 1)th assimilation step, is the parameter to be updated in ensemble member j in the nth assimilation step;

[0073] and The specific calculation formulas of are as follows:

[0074] (16)

[0075] (17)

[0076] where , , is the mean of all ensemble parameters to be updated in the nth assimilation step, is the mean of all predicted data in the nth assimilation step.

[0077] Furthermore, the two-step parameter inversion process is a process of fitting the overall data of the oilfield block and then fitting the single-well production data; wherein, the overall data of the oilfield block includes block liquid production, block oil production, block pressure, and / or block water injection volume, and the single-well production data includes single-well oil production and / or single-well water cut.

[0078] Furthermore, taking the economic net present value as an index, with the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, using the differential evolution algorithm, by constraining the total polymer injection volume and reservoir pressure, an oil reservoir injection-production optimization mathematical model is established, including:

[0079] Taking the economic net present value as the objective function, with the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, an optimal control mathematical model is established;

[0080] Taking the total polymer injection volume and reservoir pressure as constraint conditions, the optimization variables are controlled;

[0081] The particle swarm optimization algorithm is used to solve the production optimization problem and find the optimal oilfield variables that can maximize the economic net present value.

[0082] Furthermore, the objective function of the optimization control mathematical model is expressed as:

[0083] (18)

[0084] In the formula, max represents taking the maximum value, is the objective function value; are the optimization variables, including the water injection volume, polymer injection volume, and liquid production volume; is the total number of control steps; are the number of production wells, injection wells, and polymer injection wells respectively; represent the oil price, sewage treatment cost, water injection price, and polymer injection price respectively; and are the average oil production rate and water production rate of the jth production well within the nth simulation time step respectively; is the average water injection rate of the th injection well within the nth simulation time step; The average polymer injection rate of the kth polymer injection well within the nth simulation time step; b is the annual discount rate; is the unit time step length within the nth simulation time step; is the cumulative time up to the nth simulation time step;

[0085] Furthermore, with the total polymer injection volume and reservoir pressure as the constraint conditions, the process of constraining the optimization variables is expressed as:

[0086] (19)

[0087] In the formula, indicates that the total injection-production ratio of the oilfield needs to be between 0.9 and 1.1; are the lower and upper limits of the optimization variable u respectively; and are the polymer injection volume and the total polymer injection volume respectively; and are the lower and upper limits of the injection-production well control pressure respectively.

[0088] Furthermore, the particle swarm optimization algorithm is used to solve the production optimization problem and find the optimal oilfield optimization variables that can maximize the economic net present value, including:

[0089] Randomly initialize the positions and velocities of the particles, ensure within the specified range, and calculate the objective function value corresponding to the position of each particle ;

[0090] Set the initial position of each particle to its individual historical optimal position , and find the global optimal position in the population ;

[0091] In each iteration process, update the velocity of each particle and position , calculate the objective function value of the new position of the particle , update the individual optimal position :

[0092] If , then set , update the global optimal position : If , then set ;

[0093] When the maximum number of iterations is reached or the objective function value meets the preset threshold, output the global optimal position and its corresponding objective function value ;

[0094] Update the velocity of the particle through the following two core formulas and position :

[0095] (20)

[0096] (21)

[0097] S.t.

[0098] (22)

[0099] In the formula, is the inertia weight; are the learning factors of the individual and the population respectively; are random numbers between 0 and 1 respectively; is the historical optimal position of the th particle; is the global optimal position of all particles in the population; is the velocity boundary of the particle movement; is the position boundary of the particle to prevent the particle from crossing the boundary.

[0100] The beneficial effects of this application are:

[0101] In this application, an image recognition method is used to combine the sedimentary facies with the connectivity network model for modeling. While simplifying the model, it can comprehensively consider its geological characteristics, improving the accuracy of the simplified model. At the same time, this application also establishes a reservoir automatic history matching model considering sedimentary facies constraints and a mathematical model for reservoir production optimization combining water flooding and polymer flooding, which helps the oilfield to better formulate production plans, improve the oil production efficiency of the oilfield, and provide effective technical support for oilfield decision-making. Description of the Drawings

[0102] Figure 1 It is a flowchart of a rapid simulation and optimization method for a polymer flooding reservoir considering sedimentary facies constraints provided by an embodiment of this application;

[0103] Figure 2 It is a flowchart for constructing a connectivity network model provided by an embodiment of this application;

[0104] Figure 3 It is a flowchart for solving the connectivity network model provided by an embodiment of this application;

[0105] Figure 4 It is a flowchart for inverse inversion of updating key reservoir parameters provided by an embodiment of this application;

[0106] Figure 5 It is a flowchart for constructing a mathematical model for reservoir injection-production optimization provided by an embodiment of this application;

[0107] Figure 6 It is a schematic diagram of mapping from a connectivity unit to a 2D Cartesian coordinate system provided by an embodiment of this application;

[0108] Figure 7 It is a schematic diagram of sedimentary facies provided by an embodiment of this application;

[0109] Figure 8 It is a schematic diagram of modeling using a connectivity network model provided by an embodiment of this application; where, (a), schematic diagram of the connectivity network model; (b), schematic diagram after mapping the connectivity network model to a 2D Cartesian coordinate system;

[0110] Figure 9 It is a distribution diagram of the prior model during the history matching process provided by an embodiment of this application;

[0111] Figure 10 It is a distribution diagram of the posterior model during the history matching process provided by an embodiment of this application;

[0112] Figure 11 It is an optimal optimization variable diagram provided by an embodiment of this application;

[0113] Figure 12 It is a production optimization result diagram provided by an embodiment of this application. Detailed Implementation Modes

[0114] The following describes the implementation manners of the present application through specific specific examples. Those skilled in the art can easily understand other advantages and effects of the present application from the content disclosed in this specification. The present application can also be implemented or applied through other different specific implementation manners. Various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present application. It should be noted that, without conflict, the following embodiments and the features in the embodiments can be combined with each other.

[0115] The specific implementation manners of the present application will be further described in detail below in conjunction with the accompanying drawings and embodiments.

[0116] Please refer to Figure 1 , which is a flowchart of a fast simulation optimization method for a polymer flooding reservoir considering sedimentary facies constraints provided by an embodiment of the present application. An embodiment of the present application provides a fast simulation optimization method for a polymer flooding reservoir considering sedimentary facies constraints, which can be implemented through the following steps S100 to S400.

[0117] S100. Identify the sedimentary facies boundary according to the distribution characteristics of the sedimentary facies, and establish a connectivity network model considering sedimentary facies constraints based on the sedimentary facies boundary, well positions, and various sedimentary facies attributes.

[0118] In some embodiments, as Figure 2 shown, step S100 is implemented through the following steps S101 to S103.

[0119] S101. Extract the facies boundary in the sedimentary facies distribution map by using an edge detection method according to the sedimentary facies distribution of the block.

[0120] In this embodiment, the edge detection method is specifically as follows: First, convert the obtained sedimentary facies map into a grayscale image, and then use a Gaussian filter to denoise the image, that is, perform weighted averaging according to the gray values of the pixel points to be filtered and their neighborhoods. Therefore, the high-frequency noise superimposed in the image can be effectively filtered out. The Gaussian function is as follows:

[0121] (1)

[0122] In the formula, is the standard deviation controlling the Gaussian filter; f ( x , y ) is the Gaussian function; e is the natural constant.

[0123] Then use the Sobel operator to calculate the magnitude and direction of the image gradient:

[0124] (2)

[0125] (3)

[0126] (4)

[0127] Wherein, are the gradients in the x and y directions of the image respectively, is the original image, is the Sobel filter, which calculates the gradients in the horizontal and vertical directions respectively. The total gradient of each pixel point G is calculated by the formula:

[0128] (5)

[0129] Gradient direction The expression of:

[0130] (6)

[0131] For each pixel, check whether the gradient magnitude along the gradient direction is a local maximum, accurately locate the edge, and connect the scattered edge points into a complete boundary. The result is a binary image, which can be used as the initial boundary of the level set function in step S102.

[0132] S102. Use the phase boundary in the sedimentary facies distribution map as the initial boundary of the level set function, and perform regional evolution through the level set method to obtain a binary map.

[0133] In this embodiment, the core idea of the level set method is to represent the position and shape of a surface with a high-dimensional scalar function (i.e., the level set function) without explicitly describing the surface. Usually, the signed Euclidean distance from a point on the plane to the contour curve is selected as the level set function:

[0134] (7)

[0135] Wherein, represents the point on the image to the level set contour curve the distance of.

[0136] The level set function that changes with time is expressed as , and the evolution of the surface is usually described by the core equation of the level set method (the level set evolution equation):

[0137] (8)

[0138] Wherein, is the velocity function, which determines the evolution direction and velocity of the surface, represents the gradient magnitude of the level set function, t represents the actual time.

[0139] Before evolving the level set, the level set function is usually converted into an SDF function. Because during the level set evolution, after a period of time, the level set function may deviate from the distance function (i.e., ), resulting in numerical instability or abnormal surface evolution. The SDF reconstruction can restore its distance function characteristics:

[0140] (9)

[0141] (10)

[0142] In the formula, is the virtual time variable for re-initialization, which has nothing to do with the actual time t, is the initial level set function, is the sign function, which is used to keep the positive and negative signs unchanged, is a constant, which is used to avoid the denominator being zero.

[0143] Assume that there are two sedimentary facies in the reservoir model, which are respectively. After evolution, the level set can divide the entire sedimentary facies map into two regions. Therefore, the obtained clear binary map can clearly divide the regions corresponding to the two sedimentary facies.

[0144] S103. According to the actual geological information of each well point, generate an initial connected network model according to the angle and distance discrimination criteria, and the wells are connected by a one-dimensional connected unit.

[0145] In this embodiment, the specific process of generating the initial well pattern dissection map according to the angle and distance discrimination criteria is as follows: First, according to the positions of the well points in the reservoir, all pairs of wells are connected to form an initial fully connected network between wells. However, not every pair of wells in the reservoir has a connection relationship. Therefore, the introduction of the angle and distance limit conditions can screen the connections in the fully connected network. That is, when the largest angle in the triangle formed by the connections is greater than the established angle limit, it is considered that the two relatively far wells are indirectly connected through the middle well, and the longest connection is deleted; when the distance between two wells is greater than the established distance limit, it is considered that there will be no sweep effect between the two wells, and the connection is deleted. The maximum angle limit is generally taken as 120° - 160°, and the maximum distance limit is generally taken as 2 - 3 times the minimum well spacing or the average well spacing. Through these two constraint conditions, the initial well pattern dissection map can be generated.

[0146] The rules for generating the initial well network are independent of the well type and are determined only by the geodetic coordinates and relative positions of the well points. Therefore, connections between wells of the same type are allowed in the network, such as connections between injection wells and connections between production wells. In addition, once the initial connected network is established, the structure of the network will no longer change during subsequent simulations or optimization processes. If new wells need to be drilled during the simulation, they should be included in the network generation process. These wells can remain closed during the simulation until they are actually drilled.

[0147] S200, dividing the connected network model into a one-dimensional grid and converting it into a grid connection format of a general simulator, and solving it using a general simulator to obtain pressure and saturation solutions.

[0148] In some embodiments, Figure 3 As shown, step S200 is implemented by the following steps S201 to S202.

[0149] S201, dividing a one-dimensional grid in each connected unit in the connected network model, combining the binary graph with the connected network model, and assigning values ​​to each grid parameter in the one-dimensional grid.

[0150] In this embodiment, based on the initial well pattern segmentation map, each connected unit is subdivided into a series of one-dimensional grid units, each grid unit is characterized by parameters such as permeability, porosity, and water saturation, and combined with the sedimentary facies binary map obtained in step 1.2, The grids in the corresponding areas are assigned values, so that the connected network model is combined with the sedimentary phase to form a connected network model based on the sedimentary phase.

[0151] S202, converting the one-dimensional connection unit of the connected network model into the grid connection format of the general simulator, and solving with the general simulator to obtain pressure and saturation solutions.

[0152] The currently established model cannot be directly brought into a commercial numerical simulator for calculation and must be processed into a more standard format before it can be recognized. In this embodiment, all grids are mapped to a 2D Cartesian coordinate system, such as Figure 6 As shown in the figure, each row represents a connected unit, the total number of rows is the total number of connected units, each column represents the grid corresponding to different connected units, and the total number of columns is equal to the number of grids divided by each connected unit. Here, 25 grids (including well point grids, i.e. the first or last grid in each row) are taken to better characterize the sedimentary phase distribution characteristics. Because the same well exists in more than one connected unit, but the grid of the well point is fixed, it is necessary to use non-adjacent connections to connect a well grid to multiple connected units. In order to maintain the neatness of the grid matrix, we use dead grids to represent the well point grids corresponding to non-adjacent connections.

[0153] The model after being mapped to the 2D Cartesian coordinate system can be directly imported into ECLIPSE. With its powerful computing function, the pressure and saturation solutions can be obtained quickly. At the same time, since the reconstructed model is only characterized by the grids within the connected units, which is significantly lower than the parameter quantity of the full-scale reservoir model, it saves a great deal of time for the hundreds and thousands of calls to the reservoir numerical simulator during the subsequent history matching and production optimization processes.

[0154] S300. Based on sedimentary facies and reservoir constraint conditions, adopt the integrated smooth multi-data assimilation algorithm to establish an automatic history matching mathematical model and realize the updated inversion of key reservoir parameters; among them, the key reservoir parameters include main parameters such as the permeability field and relative permeability curve.

[0155] In some embodiments, as Figure 4 shown, step S300 is implemented by the following steps S301 to S303.

[0156] S301. For oil-water two-phase flow, use the Corey formula to describe the oil-water relative permeability curve:

[0157] (11)

[0158] In the formula, and respectively represent the endpoints of the oil-water relative permeability curve, and respectively represent the exponents of the oil-water relative permeability curve, which are used to describe the curvature of the curve, and respectively represent the residual oil saturation and irreducible water saturation, and respectively represent the oil saturation and water saturation, represents the relative permeability of the oil phase, represents the relative permeability of the water phase.

[0159] S302. Construct a history matching objective function and determine the constraint conditions in combination with sedimentary facies and the reservoir.

[0160] First, construct a history matching objective function. Our goal is to adjust the model parameters in step 3.1 so that the simulation results of the connected network model can match the actual production dynamics of the oilfield. Therefore, our objective function is:

[0161] (12)

[0162] In the formula, is the adjustable parameter of the model, is the objective function; For the simulation results of the connected network model; Represents the actual observed data; Represents the covariance matrix of the observation error; min represents taking the minimum value; T is the matrix transpose operation symbol.

[0163] Usually, we do not uniformly adjust all variables. This will increase the degrees of freedom of the model, and the uncertainty of the model will also increase in the case of difficult fitting. Therefore, the parameters adjusted during the fitting process mainly include the volume of the grid, grid permeability, relative permeability curve, and the depth of the oil-water contact surface. The parameters Can be expressed as follows:

[0164] (13)

[0165] In the formula, Is the grid volume, Is the grid permeability, Is the depth of the oil-water contact surface.

[0166] Similarly, assume there are two sedimentary facies in the reservoir , To reduce the randomness in the modeling process and ensure that the model is consistent with the actual geological understanding, we constrain the grid properties and reservoir conditions within each facies:

[0167] (14)

[0168] In the formula, Is the Grid permeability of the sedimentary facies ; Is the Grid permeability of the sedimentary facies ; And Are respectively the minimum and maximum permeabilities of the sedimentary facies ; And Are respectively the minimum and maximum permeabilities of the sedimentary facies ; And Are respectively the Grid volumes of the sedimentary facies and Sedimentary facies ; And Are respectively the total reservoir volumes of the sedimentary facies And Sedimentary facies And Respectively represent the Grid numbers of the sedimentary facies and Sedimentary facies And respectively represent the lowest and highest depths of the oil-water contact surface.

[0169] S303. Based on the integrated smoothing multi-data assimilation algorithm, establish an automatic history matching mathematical model to realize the two-step parameter inversion process.

[0170] The essence of the integrated smoothing multi-data assimilation algorithm (ES-MDA) method can be described as follows: For a model m, there are dimensional parameters, n observation data and M sets. In the nth assimilation step, the update of the parameter in the set member j is:

[0171] (15)

[0172] where, Cn is the covariance matrix between the parameter and the predicted data vector in the nth assimilation step, Pn is the autocovariance matrix of the predicted data in the nth assimilation step, R is the covariance matrix of the measurement error, α is the inflation coefficient, δyn represents the perturbed observation data, yn is the predicted data, θn+1,j is the parameter to be updated in the set member j in the (n + 1)th assimilation step, θn,j is the parameter to be updated in the set member j in the nth assimilation step. and The specific calculation formulas of are as follows:

[0173] (16)

[0174] (17)

[0175] where, , , θn is the mean value of all the parameters to be updated in the sets in the nth assimilation step, yn is the mean value of all the predicted data in the nth assimilation step.

[0176] Among them, the two-step parameter inversion process in step 3.3 means that in the fitting process, first, the overall data of the oilfield block is fitted, such as the liquid production volume of the block, the oil production volume of the block, the block pressure, the water injection volume of the block, etc. Secondly, the production data of a single well is fitted, such as the oil production volume of a single well, the water cut of a single well, etc. The history matching process is designed as a hierarchical iterative process, first making the model fit the overall level of the block, and then performing fine fitting.

[0177] S400. Using the economic net present value as an index, with the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, and using the differential evolution algorithm, by constraining the total polymer injection volume and reservoir pressure, an optimization mathematical model for reservoir injection and production is established.

[0178] In some embodiments, as Figure 5 shown, step S400 is implemented by the following steps S401 to S403.

[0179] S401. Taking the economic net present value as the objective function and the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, an optimization control mathematical model is established.

[0180] The objective function of the optimization control mathematical model is:

[0181] (18)

[0182] In the formula, max represents taking the maximum value, is the objective function value; are the optimization variables, including the water injection volume, polymer injection volume, and liquid production volume; is the total control step number; are the number of production wells, injection wells, and polymer injection wells, respectively; represent the oil price, sewage treatment cost, water injection price, and polymer injection price, respectively; and are the average oil production rate and water production rate of the jth production well in the nth simulation time step, respectively; is the average water injection rate of the th injection well in the nth simulation time step; The average polymer injection rate of the kth polymer injection well in the nth simulation time step; b is the annual discount rate; is the unit time step length of the nth simulation time step; is the cumulative time up to the nth simulation time step;

[0183] The physical meaning of the above formula is the profit from oil production minus the costs of water production, water injection, and polymer injection, that is, the economic net present value.

[0184] S402. Taking the total polymer injection volume and reservoir pressure as constraint conditions, control the optimization variables.

[0185] The constraint conditions for the optimization process are as follows:

[0186] (19)

[0187] In the formula, indicates that the total injection-production ratio of the oilfield needs to be between 0.9 and 1.1; are the lower and upper limits of the optimization variable u, respectively, and this value is determined according to the actual situation of the oilfield; and are the polymer injection volume and the total polymer injection volume, respectively; and are the lower limit and the upper limit of the injection-production well control pressure respectively.

[0188] S403. Use the particle swarm optimization algorithm to solve the production optimization problem and find the optimal optimization variables of the oilfield that can maximize the economic net present value.

[0189] The basic principle of the particle swarm optimization algorithm is as follows:

[0190] The particle swarm optimization algorithm is an optimization algorithm based on swarm intelligence proposed by Kennedy and Eberhart. It imitates the social behavior of flocks of birds foraging or schools of fish moving, and realizes the search for the global optimal solution through information sharing among individuals. Each particle in the algorithm represents a possible solution in the search space, and the particles dynamically adjust their positions through their own and the group's experience, thus gradually approaching the optimal solution.

[0191] The particle swarm optimization algorithm updates the velocity of the particle through the following two core formulas and position ;

[0192] (20)

[0193] (21)

[0194] S.t.

[0195] (22)

[0196] In the formula, is the inertia weight; are the learning factors of the individual and the group, respectively; are random numbers between 0 and 1, respectively; is the historical optimal position of the th particle; is the global optimal position of all particles in the group; is the velocity boundary of the particle movement; is the position boundary of the particle to prevent the particle from crossing the boundary.

[0197] The physical meanings of the three terms on the right side of the particle update velocity formula are as follows: retaining the original motion inertia of the particle helps to explore new areas; encouraging the particle to approach its own historical optimal position; and pushing the particle to approach the global optimal position of the group.

[0198] First, randomly initialize the position and velocity , ensure that within the specified range, the objective function values corresponding to the positions of each particle are calculated , set the initial position of each particle as its individual historical optimal position , and find the global optimal position in the population . In each iteration process, update the velocity of each particle and position , calculate the objective function value of the new position of the particle , update the individual optimal position : If , then set , update the global optimal position : If , then set . When the maximum number of iterations is reached or the objective function value meets the preset threshold, the algorithm terminates, and the global optimal position and its corresponding objective function value are output , which is the optimization result that makes the economic net present value optimal.

[0199] To prove the feasibility and superiority of this application, the following embodiments are given.

[0200] First, use a full-scale reservoir numerical simulator to construct a channel reservoir model with a grid scale of 40*40*1, a model size of 400m*400m*5m. The reservoir sedimentary facies map is as shown in Figure 7 . The reservoir is divided into three parts by a high-permeability sedimentary facies belt from the upper left to the lower right. The white channel is the high-permeability sedimentary facies belt with a porosity of 0.2 and a permeability of 100 mD. The remaining black background part is the low-permeability sedimentary facies belt with a porosity of 0.15 and a permeability of 10 mD. The average depth of the reservoir is 4000m, the depth of the oil-water interface is 4200m, and the viscosities of oil and water are 20 mPa*s and 1 mPa*s respectively. In this model, we did not fit the relative permeability of oil and water, but modeled it as known conditions. There are a total of four wells in the reservoir. One injection well is located in the lower left, and three production wells are located in the other three corners respectively. The reservoir simulator was used to simulate for 6000 days. The injection well injected water in the first 3000 days, and polymer with a concentration of 1000 mg / L was injected starting from the 3001st day. The production history of the first 5000 days was used as the observed data for history matching, and the production history of the last 1000 days was used as the prediction data to verify the prediction effect of the model.

[0201] Similarly, using the method described in the present invention, based on the sedimentary facies map shown in Figure 7 , after obtaining the binary sedimentary facies map by combining edge detection and level set function, establish as shown in Figure 8The connected network model shown, its connected network schematic diagram and the model after mapping to the 2D Cartesian coordinate system are respectively as Figure 8 shown in (a) and (b) in the figure. The ES-MDA algorithm is used for history matching. The population size is set to 150, and the number of iterations is 4. Taking production well 1 as an example to show the model effect, its prior distribution is as Figure 9 shown. Among them, the dotted line divides the history matching period and the prediction period. The black line represents the observed result, and the gray line represents the prior model simulation result; its posterior distribution is as Figure 10 shown. In the figure, the black line represents the observed result, and the gray line represents the posterior model simulation result; it can be seen that after history matching, both the history matching result and the prediction result of the model are good, verifying the ability of the model in water flooding and polymer flooding simulation.

[0202] Finally, the particle swarm optimization algorithm is used to optimize the polymer injection volume and liquid production volume. Taking 100 days as an optimization time step, optimizing for ten steps, a total of 1000 days are optimized. The optimal optimization variables are as Figure 11 shown. The comparison of the cumulative oil production after optimization is as Figure 12 shown. It can be seen that after optimization, it is obvious that a higher cumulative oil production can be obtained than before optimization. Therefore, the reliability of the method proposed in this application in the process of considering sedimentary facies modeling and optimization is verified.

[0203] The above embodiments are only used to illustrate the present application, rather than limiting the present application. Those of ordinary skill in the relevant technical fields can also make various changes and modifications without departing from the spirit and scope of the present application. Therefore, all equivalent technical solutions also belong to the scope of the present application. The patent protection scope of the present application shall be defined by the claims.

Claims

1. A rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints, characterized in that The method includes: Identifying sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and establishing a connectivity network model considering sedimentary facies constraints based on the sedimentary facies boundaries, well locations, and sedimentary facies attributes; Performing one-dimensional grid division on the connectivity network model, converting it into the grid connection format of a general simulator, and using the general simulator for solution to obtain pressure and saturation solutions; Based on sedimentary facies and reservoir constraints, using an integrated smooth multi-data assimilation algorithm to establish an automatic history matching mathematical model and realize the updated inversion of key reservoir parameters; where the key reservoir parameters include grid size, grid permeability, and relative permeability curves; Taking the economic net present value as an index, with the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, using the differential evolution algorithm, and establishing a reservoir injection-production optimization mathematical model by constraining the total polymer injection volume and reservoir pressure; Identifying sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and establishing a connectivity network model considering sedimentary facies constraints based on the sedimentary facies boundaries, well locations, and sedimentary facies attributes, including: Extracting the facies boundaries in the sedimentary facies distribution map using an edge detection method according to the sedimentary facies distribution of the block; Taking the facies boundaries in the sedimentary facies distribution map as the initial boundaries of the level set function, and performing regional evolution through the level set method to obtain a binary image, which is used to divide the regions corresponding to multiple sedimentary facies; Generating an initial connectivity network model according to the actual geological information of each well point based on the angle and distance discrimination criteria, and connecting wells with a one-dimensional connectivity unit; Extracting the facies boundaries in the sedimentary facies distribution map using an edge detection method according to the sedimentary facies distribution of the block, including: Converting the obtained sedimentary facies map into a grayscale image, and using a Gaussian function to perform weighted averaging according to the gray values of the pixel points to be filtered and their neighborhoods to filter out the superimposed high-frequency noise in the image. The Gaussian function is as follows: (1) In the formula, is the standard deviation for controlling the Gaussian filter; f ( x , y ) is the Gaussian function; e is the natural constant; Using the Sobel operator to calculate the magnitude and direction of the image gradient through the following formulas (2), (3), and (4): (2) (3) (4) In the formula, are the gradients in the x and y directions of the image respectively, is the original image, is the Sobel filter, which calculates the gradients in the horizontal and vertical directions respectively. The total gradient G of each pixel is calculated by the formula: (5) Gradient direction The expression of: (6) For each pixel, checking whether its gradient amplitude is a local maximum along the gradient direction, accurately locating the edges, connecting the scattered edge points into a complete boundary to obtain a binary image, and taking the binary image as the facies boundaries in the sedimentary facies distribution map; Performing one-dimensional grid division on the connectivity network model, converting it into the grid connection format of a general simulator, and using the general simulator for solution to obtain pressure and saturation solutions, including: Dividing one-dimensional grids in each connectivity unit of the connectivity network model, combining the binary image with the connectivity network model, and assigning values to each grid parameter in the one-dimensional grids; converting the one-dimensional connection units of the connectivity network model into the grid connection format of a general simulator, and using the general simulator for solution to obtain pressure and saturation solutions; Based on the initial well pattern dissection diagram, each connected unit is subdivided into a series of one-dimensional grids, and each one-dimensional grid is characterized by grid parameters; the specific parameters include permeability, porosity, and / or water saturation; according to the binary map, the one-dimensional grids in the areas corresponding to different sedimentary facies are assigned values respectively, so as to combine the connected network model with the sedimentary facies for modeling, and a connected network model based on sedimentary facies is formed; the connected network model is mapped into a 2D Cartesian coordinate system, each row represents a connected unit, the total number of rows is the total number of connected units, each column represents the grids at the corresponding positions of different connection units, and the total number of columns is equal to the number of grids divided by each connected unit, and the well point grids corresponding to non-adjacent connections are represented by dead grids; the connected network model after being mapped into the 2D Cartesian coordinate system is brought into a general simulator to quickly obtain pressure and saturation solutions.

2. The rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints as described in claim 1, wherein Taking the phase boundary in the sedimentary facies distribution map as the initialization boundary of the level set function, regional evolution is carried out by the level set method to obtain a binary map, including: The signed Euclidean distance from a point on the plane to the contour curve is selected as the level set function : (7) In the formula, represents a point on the image to the level set contour curve distance; The level set function that changes with time is expressed as , and the evolution of the surface is described by the level set evolution equation shown below: (8) wherein, is the velocity function, which determines the evolution direction and velocity of the surface, represents the gradient magnitude of the level set function, t represents the actual time; Before evolving the level set, the level set function is transformed into an SDF function by the following formula: (9) (10) In the formula, is the reinitialized virtual time variable, which has nothing to do with the actual time t. is the initial level set function. is the sign function, which is used to keep the positive and negative signs unchanged. is a constant, which is used to avoid a zero denominator. Based on the level set evolution equation, level set evolution is carried out to obtain a binary map.

3. The rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints as described in claim 1, wherein Based on the sedimentary facies and reservoir constraint conditions, an integrated smooth multi-data assimilation algorithm is adopted to establish an automatic history matching mathematical model to realize the updated inversion of key reservoir parameters, including: Construct a history matching objective function and determine the constraint conditions in combination with the sedimentary facies and reservoir; Based on the integrated smooth multi-data assimilation algorithm, establish an automatic history matching mathematical model to realize the two-step parameter inversion process.

4. The rapid simulation optimization method for polymer flooding reservoirs considering sedimentary facies constraints according to claim 3, characterized in that Construct a history matching objective function and determine the constraint conditions in combination with the sedimentary facies and reservoir, including: The constructed history matching objective function is expressed as: (11) In the formula, is the adjustable parameter of the model, is the objective function; is the simulation result of the connected network model; represents the actual observed data; represents the covariance matrix of the observation error; min represents taking the minimum value; T is the symbol for matrix transpose operation; The adjustable parameters of the model are represented as: (12) In the formula, is the grid volume, is the grid permeability, is the depth of the oil-water contact surface; and respectively represent the endpoints of the relative permeability curves of oil and water, and respectively represent the exponents of the relative permeability curves of oil and water, which are used to describe the curvature of the curves; The grid attributes and reservoir conditions in each phase are constrained by the following formula: (13) Wherein, is the grid permeability of the sedimentary facies; is the grid permeability of the sedimentary facies; and are respectively the minimum and maximum permeabilities of the sedimentary facies; and are respectively the minimum and maximum permeabilities of the sedimentary facies; and are respectively the grid volumes of the sedimentary facies and sedimentary facies; and are respectively the total reservoir volumes of the sedimentary facies and sedimentary facies; and respectively represent the numbers of grids of the sedimentary facies and and respectively represent the minimum and maximum depths of the oil-water contact surface.

5. The polymer flooding reservoir rapid simulation optimization method considering sedimentary facies constraints according to claim 4, characterized in that, The two-step parameter inversion process is a process of fitting the overall data of the oilfield block first and then fitting the single-well production data; among them, the overall data of the oilfield block includes block liquid production, block oil production, block pressure, and / or block water injection volume, and the single-well production data includes single-well oil production and / or single-well water cut.

6. The polymer flooding reservoir rapid simulation optimization method considering sedimentary facies constraints according to claim 1, characterized in that Taking the economic net present value as an index, the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, using the differential evolution algorithm, and by constraining the total polymer injection volume and reservoir pressure, establish an oil reservoir injection-production optimization mathematical model, including: Taking the economic net present value as the objective function and the water injection volume, liquid production volume, and polymer injection volume of each well as optimization variables, establish an optimal control mathematical model; Taking the total polymer injection volume and reservoir pressure as constraint conditions, the optimization variables are constrained; The particle swarm algorithm is used to solve the production optimization problem to find the optimal optimization variables of the oilfield that can maximize the economic net present value.

7. The polymer flooding reservoir rapid simulation optimization method considering sedimentary facies constraints according to claim 6, characterized in that The objective function of the optimal control mathematical model is expressed as: (14) Where max represents taking the maximum value, is the objective function value; are the optimization variables, including water injection volume, polymer injection volume, and liquid production volume; is the total number of control steps; are the number of production wells, injection wells, and polymer injection wells respectively; represent the oil price, sewage treatment cost, water injection price, and polymer injection price respectively; and are the average oil production rate and water production rate of the j - th production well in the n - th simulation time step respectively; is the average water injection rate of the - th injection well in the n - th simulation time step; is the average polymer injection rate of the k - th polymer injection well in the n - th simulation time step; b is the annual discount rate; is the unit time step length of the n - th simulation time step; is the cumulative time up to the n - th simulation time step.

Images