Rapid simulation and optimization method for polymer flooding reservoir considering sedimentary facies constraint

Abstract

The invention relates to the technical field of reservoir numerical simulation, and discloses a polymer flooding reservoir rapid simulation optimization method considering sedimentary facies constraints, which comprises the following steps: identifying sedimentary facies boundaries according to distribution characteristics of sedimentary facies, and establishing a connected network model considering the sedimentary facies constraints on the basis of the sedimentary facies boundaries, well positions and various sedimentary facies attributes; carrying out one-dimensional grid division on the connected network model, converting the connected network model into a grid connection format of a universal simulator, and solving the connected network model by adopting the universal simulator to obtain pressure and saturation solutions; establishing an automatic historical fitting mathematical model based on sedimentary facies and oil reservoir constraint conditions; the economic net present value is used as an index, the water injection rate, the liquid production rate and the polymer injection rate of each well are used as optimization variables, the differential evolution algorithm is used, the total polymer injection amount and the oil reservoir pressure are restrained, and an oil reservoir injection and production optimization mathematical model is established. According to the method, the oil reservoir is optimized by comprehensively considering the characteristics of the oil reservoir after water flooding is converted into polymer flooding, so that rapid decision-making of the oil field and guidance of oil field construction are realized.

Description

Technical Field

[0001] The present application relates to the technical field of reservoir numerical simulation, and in particular to a rapid simulation optimization method for polymer flooding reservoirs taking sedimentary phase constraints into consideration. Background Art

[0002] After years of water flooding, oil fields generally face problems such as high water content, high decline, and low recovery rate. The water flooding efficiency of water flooding development has dropped significantly. In order to improve the benefits of oil field development, various enhanced oil recovery (EOR) methods have been developed. Among them, polymer flooding technology has become one of the important EOR methods because it can effectively improve the migration characteristics of reservoir fluids and significantly improve the recovery rate. The core principle of polymer injection is to increase the viscosity of injected water and reduce the seepage capacity of the water phase during the displacement process, thereby improving the swept volume and improving the oil recovery efficiency. Compared with traditional water flooding, polymer flooding can effectively reduce water channeling, reduce water content, and significantly improve crude oil recovery. However, due to the high cost of polymers, how to optimize the injection and production plan and reduce the development cost while ensuring the improvement of recovery rate has become an urgent problem to be solved in the current oil field 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 maximize crude oil production and economic benefits while reducing development costs by constructing an optimization objective function and reasonably adjusting control variables such as well control. In the process of converting water flooding to polymer flooding, how to achieve efficient numerical simulation and intelligent optimization decision-making to more scientifically guide polymer injection strategies is the key to improving reservoir development results. Summary of the invention

[0003] The purpose of this application is to provide a rapid simulation optimization method for polymer flooding reservoirs considering sedimentary phase constraints, and to build an efficient water-to-polymer flooding optimization framework through connected network modeling, numerical simulation solution, automatic history matching and intelligent production optimization of sedimentary phase constraints. This method can quickly predict the effect of polymer flooding, optimize the injection-production strategy, and reduce development costs, providing scientific support for the efficient development of oil fields.

[0004] In general, firstly, according to the sedimentary phase distribution of the reservoir, the image recognition method combining edge detection and level set function is used to identify the boundaries between the sedimentary phases and generate a binary image. Then, according to the actual geological information of each well point, the initial well network segmentation map is generated according to the angle and distance judgment criteria. The wells are connected by a one-dimensional connection unit, and the grid is divided in the one-dimensional connection unit. On this basis, combined with the binary image, the parameters such as the permeability of each phase are assigned in each grid to form an inter-well connection network model. The constructed model is mapped to the 2D simulation grid and coupled with the commercial numerical simulator (such as ECLIPSE) for calculation, combining simple modeling with a fast-calculating simulator. The history fitting algorithm ES-MDA is used. The model parameters such as grid size and permeability are updated, and constraints are added to fit the production data under the premise that the model conforms to the geological characteristics; the fitted model is regarded as a real reservoir model, and the characteristics of water flooding followed by polymer flooding are considered. The particle swarm algorithm is used to constrain the parameters of each well and block, with the economic net present value as the objective function, while obtaining a development strategy that can maximize the benefits of oilfield development, and formulating a production plan to guide the subsequent production process of the oilfield.

[0005] In order to achieve the above purpose, the technical solutions adopted are as follows: The present application provides a rapid simulation optimization method for polymer flooding reservoirs considering sedimentary phase constraints, the method comprising: Identifying sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and establishing a connectivity network model taking sedimentary facies constraints into account based on the sedimentary facies boundaries, well locations and attributes of each sedimentary facies; Dividing the connected network model into a one-dimensional grid and converting it into a grid connection format of a general simulator, solving it using a general simulator, and obtaining pressure and saturation solutions; Based on sedimentary facies and reservoir constraints, an integrated smooth multi-data assimilation algorithm is used to establish an automatic history fitting mathematical model to achieve updated inversion of key reservoir parameters; wherein the key reservoir parameters include permeability field and relative permeability curve; Taking the economic net present value as the indicator, the water injection volume, liquid production volume and polymer injection volume of each well as the optimization variables, the differential evolution algorithm was used to establish a mathematical model for reservoir injection and production optimization by constraining the total polymer injection volume and reservoir pressure.

[0006] Furthermore, the sedimentary facies boundaries are identified according to the distribution characteristics of the sedimentary facies, and based on the sedimentary facies boundaries, well locations and attributes of each sedimentary facies, a connected network model considering sedimentary facies constraints is established, including: According to the sedimentary facies distribution of the block, the phase boundaries in the sedimentary facies distribution map are extracted using edge detection method; The phase boundary in the sedimentary phase distribution map is used as the initialization boundary of the level set function, and the regional evolution is performed by the level set method to obtain a binary map; According to the actual geological information of each well point, the initial connected network model is generated according to the angle and distance judgment criteria, and the wells are connected by a one-dimensional connected unit.

[0007] Furthermore, according to the sedimentary facies distribution of the block, the phase boundaries in the sedimentary facies distribution map are extracted using edge detection methods, including: The obtained sedimentary phase map is converted into a grayscale image, and a Gaussian function is used to perform weighted averaging according to the grayscale values ​​of the pixel to be filtered and its neighborhood to filter out the high-frequency noise superimposed in the image. The Gaussian function is as follows: (1) In the formula, is the standard deviation of the controlling Gaussian filter; f ( x , y ) is a Gaussian function; e is a natural constant; Using the Sobel operator, the magnitude and direction of the image gradient are calculated using the following formulas (2), (3) and (4): (2) (3) (4) In the formula, are the gradients of the image in the x and y directions, is the original image, It is a Sobel filter, which calculates the gradients in the horizontal and vertical directions respectively, and the total gradient of each pixel G The calculation formula is: (5) Gradient direction The expression is: (6) For each pixel, check whether its gradient amplitude is the local maximum along the gradient direction, accurately locate the edge, connect the scattered edge points into a complete boundary, and obtain a binary image, which is used as the phase boundary in the sedimentary phase distribution map.

[0008] Furthermore, the phase boundary in the sedimentary phase distribution map is used as the initialization boundary of the level set function, and the regional evolution is performed by the level set method to obtain a binary map, including: Select the signed Euclidean distance from the point on the plane to the contour curve as the level set function : (7) In the formula, Represents a point on an image To level set contour curve distance; The time-varying level set function is expressed as , the evolution of the surface is described by the level set evolution equation as follows: (8) In the formula, is the velocity function, which determines the evolution direction and speed of the surface. represents the gradient magnitude of the level set function, t Indicates actual time; Before evolving the level set, the level set function is converted into the SDF function by the following formula: (9) (10) In the formula, is a reinitialized virtual time variable, which has nothing to do with the actual time t. is the initial level set function, is a sign function, used to keep the positive and negative signs unchanged, is a constant used to avoid the denominator being 0; Level set evolution is performed based on the level set evolution equation to obtain a binary map, and the binary map is used to divide regions corresponding to multiple depositions.

[0009] 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: 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; 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; 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.

[0010] 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: 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; 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.

[0011] 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: 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; 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.

[0012] 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: The connected network model is mapped to a 2D Cartesian coordinate system, where 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. The well point grids corresponding to non-adjacent connections are represented by dead grids. The connected network model mapped to a 2D Cartesian coordinate system is then brought into a general simulator to quickly obtain pressure and saturation solutions.

[0013] Furthermore, based on the sedimentary facies and reservoir constraints, an integrated smooth multi-data assimilation algorithm is used to establish an automatic history matching mathematical model to achieve updated inversion of key reservoir parameters, including: Construct the history matching objective function and determine the constraints in combination with the sedimentary facies and reservoirs; Based on the integrated smooth multi-data assimilation algorithm, an automatic history fitting mathematical model is established to realize the two-step parameter inversion process.

[0014] Furthermore, for oil-water two-phase flow, the oil-water relative permeability curve is described by the following formula: (11) In the formula, and represent the endpoints of the oil-water relative permeability curve, and Respectively represent the index of the oil-water relative permeability curve, which is used to describe the curvature of the curve. and represent the residual oil saturation and irreducible water saturation, respectively. and represent oil saturation and water saturation respectively, represents the relative permeability of the oil phase, Represents the relative permeability of water phase.

[0015] Furthermore, the historical matching objective function is constructed, and the constraints are determined in combination with the sedimentary facies and reservoirs, including: Construct the history matching objective function, expressed as: (12) In the formula, are adjustable parameters of the model, is the objective function; Simulation results for the connected network model; represents the actual observed data; Represents the covariance matrix of the observation error; min means taking the minimum value; T is the matrix transposition operator; The model can adjust the parameters It is expressed as: (13) In the formula, is the grid volume, is the grid permeability, is the depth of oil-water contact surface; The grid attributes and reservoir conditions in each phase are constrained by the following formula: (14) In the formula, For Sedimentary Grid permeability; For Sedimentary Grid permeability; and They are respectively minimum and maximum permeabilities of sedimentary facies; and They are respectively minimum and maximum permeabilities of sedimentary facies; and They are respectively Sedimentary facies and Sedimentary Mesh volume; and They are respectively Sedimentary facies and Total reservoir volume of sedimentary facies; and Respectively indicate the Sedimentary facies and The number of grids for sedimentary facies; and Respectively represent the lowest and highest depths of the oil-water contact surface.

[0016] 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: For a model m, we have Dimension parameters, Observation data and There are sets. In the nth assimilation step, the update of the parameters in the set member j is: (15) In the formula, is the covariance matrix between the parameters 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 measurement errors, is the coefficient of expansion, represents the perturbed observation data, To predict data, is the parameter to be updated in the set member j in the n+1th assimilation step, is the parameter to be updated in the set member j in the nth assimilation step; and The specific calculation formula is as follows: (16) (17) In the formula, , , is the mean of all the parameters to be updated in the nth assimilation step, is the mean of all predicted data in the nth assimilation step.

[0017] 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 the block liquid production, the block oil production, the block pressure and / or the block water injection volume, and the single well production data includes the single well oil production and / or the single well water content.

[0018] Furthermore, taking the economic net present value as the indicator, the water injection volume, liquid production volume and polymer injection volume of each well as the optimization variables, the differential evolution algorithm is used to constrain the total polymer injection volume and reservoir pressure to establish a reservoir injection and 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 the optimization variables, an optimal control mathematical model is established; Taking the total amount of polymer injection and reservoir pressure as constraints, the optimization variables are controlled; The particle swarm algorithm is used to solve the production optimization problem and find the oilfield optimization variables that can maximize the economic net present value.

[0019] Furthermore, the objective function of the optimal control mathematical model is expressed as: (18) In the formula, max means taking the maximum value, is the objective function value; are optimization variables, including water injection rate, polymer injection rate, and liquid production rate; is the total number of control steps; They are the number of production wells, water injection wells, and polymer injection wells; represent 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 number of Average water injection rate of injection wells; The average polymer injection rate of the kth polymer injection well in the nth simulation time step; b is the annual discount rate; The unit time step is the nth simulation time step; is the cumulative time up to the n simulation time step; Furthermore, taking the total amount of polymer injection and reservoir pressure as constraints, the process of constraining the optimization variables is expressed as: (19) In the formula, This means that the total injection-production ratio of the oil field needs to be between 0.9 and 1.1; are the lower and upper bounds of the optimization variable u respectively; and are polymer injection amount and total polymer injection amount, respectively; and The control pressure of injection and production wells are lower and upper limits.

[0020] Furthermore, the particle swarm algorithm is used to solve the production optimization problem and find the optimal optimization variables of the oil field that can maximize the economic net present value, including: Randomly initialize the particle positions and speed , ensuring that the objective function value corresponding to the position of each particle is calculated within the specified range ; Set the initial position of each particle to its individual historical optimal position , and find the global optimal position in the group ; During each iteration, update the velocity of each particle and location , calculate the objective function value of the new position of the particle , update the individual optimal position : like , then suppose , update the global optimal position :like , then suppose ; When the maximum number of iterations is reached or the objective function value meets the preset threshold, the global optimal position is output. and its corresponding objective function value ; Update the particle's velocity using the following two core formulas and location : (20) (twenty one) St. (twenty two) In the formula, is the inertia weight; are the learning factors for individuals and groups, respectively; are random numbers between 0 and 1 respectively; For the The historical optimal position of a particle; is the global optimal position of all particles in the group; is the velocity boundary of particle motion; It is the position boundary of the particle to prevent the particle from crossing the boundary.

[0021] The beneficial effects of this application are: This application combines the sedimentary phase and the connected network model through image recognition methods to build a model, which can comprehensively consider its geological characteristics while simplifying the model, thereby improving the accuracy of the simplified model. At the same time, this application also establishes an automatic history fitting model for oil reservoirs that considers sedimentary phase constraints and a mathematical model for oil reservoir production optimization that combines water drive and polymer drive, which helps oil fields to better formulate production plans, improve oil field oil recovery efficiency, and provide effective technical support for oil field decision-making. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 A flowchart of a rapid simulation optimization method for polymer flooding reservoirs considering sedimentary phase constraints provided in an embodiment of the present application; Figure 2 A flow chart for constructing a connected network model provided in an embodiment of the present application; Figure 3 A flow chart for solving a connected network model provided in an embodiment of the present application; Figure 4 A flowchart for updating and inversion of key reservoir parameters provided in an embodiment of the present application; Figure 5 A flow chart for constructing a mathematical model for reservoir injection and production optimization provided in an embodiment of the present application; Figure 6 A schematic diagram of mapping a connected unit to a 2D Cartesian coordinate system provided in an embodiment of the present application; Figure 7 A schematic diagram of the deposition phase provided for the embodiment of the present application; Figure 8 A schematic diagram of modeling using a connected network model provided in an embodiment of the present application; wherein (a) is a schematic diagram of a connected network model; (b) is a schematic diagram of a connected network model mapped to a 2D Cartesian coordinate system; Fig. 9 A distribution diagram of a priori models in the history matching process provided by an embodiment of the present application; Fig.10 A distribution diagram of the posterior model in the history matching process provided in the embodiment of the present application; Fig.11 The optimal optimization variable diagram provided by the embodiment of the present application; Fig.12 This is a diagram of the production optimization results provided in an embodiment of the present application. DETAILED DESCRIPTION

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

[0024] The specific implementation of the present application is further described in detail below in conjunction with the drawings and examples.

[0025] See also Figure 1 , which is a flow chart of a rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints provided in an embodiment of the present application. An embodiment of the present application provides a rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints, which can be implemented through the following steps S100 to S400.

[0026] S100, identifying sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and establishing a connected network model that takes sedimentary facies constraints into consideration based on the sedimentary facies boundaries, well locations, and attributes of each sedimentary facies.

[0027] In some embodiments, Figure 2 As shown, step S100 is implemented by the following steps S101 to S103.

[0028] S101. According to the sedimentary phase distribution of the block, the phase boundary in the sedimentary phase distribution map is extracted by using an edge detection method.

[0029] In this embodiment, the edge detection method is as follows: first, the obtained sedimentary phase image is converted into a grayscale image, and then a Gaussian filter is used to denoise the image, that is, a weighted average is performed according to the grayscale values ​​of the pixel to be filtered and its neighborhood, so that the high-frequency noise superimposed in the image can be effectively filtered out. The Gaussian function is as follows: (1) In the formula, is the standard deviation of the controlling Gaussian filter; f ( x , y ) is a Gaussian function; e is a natural constant.

[0030] Then use the Sobel operator to calculate the size and direction of the image gradient: (2) (3) (4) In the formula, are the gradients of the image in the x and y directions, is the original image, It is a Sobel filter, which calculates the gradients in the horizontal and vertical directions respectively, and the total gradient of each pixel G The calculation formula is: (5) Gradient direction The expression is: (6) For each pixel, check whether its gradient amplitude is a local maximum along the gradient direction, 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.

[0031] S102, using the phase boundary in the sedimentary phase distribution map as the initialization boundary of the level set function, performing regional evolution through the level set method, and obtaining a binary map.

[0032] In this embodiment, the core idea of ​​the level set method is to use a high-dimensional scalar function (i.e., the level set function) to represent the position and shape of the surface 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: (7) In the formula, Represents a point on an image To level set contour curve distance.

[0033] The time-varying level set function is expressed as , the evolution of the surface is usually described by the core equation of the level set method (level set evolution equation): (8) In the formula, is the velocity function, which determines the evolution direction and speed of the surface. represents the gradient magnitude of the level set function, t Indicates the actual time.

[0034] Before evolving the level set, it is usually necessary to convert the level set function into an SDF function, because in the level set evolution, after a period of time, the level set function may deviate from the distance function (i.e. ), which leads to numerical instability or abnormal surface evolution. SDF reconstruction can restore its distance function characteristics: (9) (10) In the formula, is a reinitialized virtual time variable, which has nothing to do with the actual time t. is the initial level set function, is a sign function, used to keep the positive and negative signs unchanged, is a constant used to avoid the denominator being zero.

[0035] Assume that there are two sedimentary phases in the reservoir model, namely The evolved level set can divide the entire sedimentary phase map into two regions, so the obtained clear binary map can clearly divide the two corresponding sedimentary regions.

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

[0037] In this embodiment, the initial well network segmentation diagram is generated according to the angle and distance judgment criteria as follows: First, 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. However, not every well in the reservoir is connected. Therefore, the introduction of angle and distance constraints can screen the connections in the fully connected network. That is, when the maximum angle of the triangle formed by the connection is greater than the established angle limit, it is considered that the two farther wells are indirectly connected through the middle well, and the longest connection is deleted; when the distance between the two wells is greater than the established distance limit, it is considered that there will be no spillover effect between the two wells, and the connection is deleted. The maximum angle limit is generally 120°-160°, and the maximum distance limit is generally 2-3 times the minimum well spacing or the average well spacing. The initial well network segmentation diagram can be generated by these two constraints.

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

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

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

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

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

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

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

[0045] The model mapped to the 2D Cartesian coordinate system can be directly brought into ECLIPSE, and with its powerful computing function, pressure and saturation solutions can be quickly obtained. At the same time, since the reconstructed model is only represented by the grids in the connected units, the number of parameters is significantly lower than that of the full-scale reservoir model, which saves a lot of time for the subsequent history matching and production optimization process to call the reservoir numerical simulator hundreds or thousands of times.

[0046] S300, based on sedimentary facies and reservoir constraints, uses an integrated smooth multi-data assimilation algorithm to establish an automatic history fitting mathematical model to achieve updated inversion of key reservoir parameters; among them, key reservoir parameters include main parameters such as permeability field and relative permeability curve.

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

[0048] S301. For oil-water two-phase flow, the Corey formula is used to describe the oil-water relative permeability curve: (11) In the formula, and represent the endpoints of the oil-water relative permeability curve, and Respectively represent the index of the oil-water relative permeability curve, which is used to describe the curvature of the curve. and represent the residual oil saturation and irreducible water saturation, respectively. and represent oil saturation and water saturation respectively, represents the relative permeability of the oil phase, Represents the relative permeability of water phase.

[0049] S302, constructing a history matching objective function, and determining constraint conditions in combination with sedimentary facies and oil reservoirs.

[0050] First, we construct the 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 oil field. Therefore, our objective function is: (12) In the formula, are adjustable parameters of the model, is the objective function; Simulation results for the connected network model; represents the actual observed data; Represents the covariance matrix of the observation error; min means taking the minimum value; T is the matrix transpose operator symbol.

[0051] Usually we do not adjust all variables uniformly, which will increase the degree of freedom of the model and increase the uncertainty of the model when it is difficult to fit. Therefore, the parameters we adjust during the fitting process mainly include the volume of the grid, the grid permeability, the relative permeability curve and the depth of the oil-water contact surface. It can be expressed as follows: (13) In the formula, is the grid volume, is the grid permeability, is the oil-water contact depth.

[0052] Similarly, assuming there are two sedimentary phases in the reservoir In order to reduce the randomness of the modeling process and ensure that the model is consistent with the actual geological understanding, we constrain the grid attributes and reservoir conditions in each phase: (14) In the formula, For Sedimentary Grid permeability; For Sedimentary Grid permeability; and They are respectively minimum and maximum permeabilities of sedimentary facies; and They are respectively minimum and maximum permeabilities of sedimentary facies; and They are respectively Sedimentary facies and Sedimentary Mesh volume; and They are respectively Sedimentary facies and Total reservoir volume of sedimentary facies; and Respectively indicate the Sedimentary facies and The number of grids for sedimentary facies; and Respectively represent the lowest and highest depths of the oil-water contact surface.

[0053] S303. Based on the integrated smooth multi-data assimilation algorithm, an automatic history fitting mathematical model is established to realize the two-step parameter inversion process.

[0054] The essence of the ES-MDA method can be described as follows: for a model m, Dimension parameters, Observation data and There are sets. In the nth assimilation step, the update of the parameters in the set member j is: (15) In the formula, is the covariance matrix between the parameters 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 measurement errors, is the coefficient of expansion, represents the perturbed observation data, To predict data, is the parameter to be updated in the set member j in the n+1th assimilation step, is the parameter to be updated in the set member j in the nth assimilation step. and The specific calculation formula is as follows: (16) (17) In the formula, , , is the mean of all the parameters to be updated in the nth assimilation step, is the mean of all predicted data in the nth assimilation step.

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

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

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

[0058] 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 the optimization variables, an optimal control mathematical model is established.

[0059] The objective function of the optimal control mathematical model is: (18) In the formula, max means taking the maximum value, is the objective function value; are optimization variables, including water injection rate, polymer injection rate, and liquid production rate; is the total number of control steps; They are the number of production wells, water injection wells, and polymer injection wells; represent 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 number of Average water injection rate of injection wells; The average polymer injection rate of the kth polymer injection well in the nth simulation time step; b is the annual discount rate; The unit time step is the nth simulation time step; is the cumulative time up to the n simulation time step; The physical meaning of the above formula is the profit from oil production minus the cost of water production, water injection and polymer injection, that is, the economic net present value.

[0060] S402. Taking the total amount of polymer injection and reservoir pressure as constraints, the optimization variables are controlled.

[0061] The constraints of the optimization process are as follows: (19) In the formula, This means that the total injection-production ratio of the oil field needs to be between 0.9 and 1.1; are the lower and upper limits of the optimization variable u, respectively. This value is determined according to the actual situation of the oil field; and are polymer injection amount and total polymer injection amount, respectively; and The control pressure of injection and production wells are lower and upper limits.

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

[0063] The basic principle of particle swarm algorithm is as follows: Particle swarm optimization algorithm is an optimization algorithm based on swarm intelligence proposed by Kennedy and Eberhart. It imitates the social behavior of bird flocks foraging or fish schools moving, and realizes the search for the global optimal solution through information sharing between individuals. Each particle in the algorithm represents a possible solution in the search space, and the particle dynamically adjusts its position based on its own and the group's experience, thereby gradually approaching the optimal solution.

[0064] The particle swarm algorithm updates the particle speed through the following two core formulas: and location ; (20) (twenty one) St. (twenty two) In the formula, is the inertia weight; are the learning factors for individuals and groups, respectively; are random numbers between 0 and 1 respectively; For the The historical optimal position of a particle; is the global optimal position of all particles in the group; is the velocity boundary of particle motion; It is the position boundary of the particle to prevent the particle from crossing the boundary.

[0065] The physical meanings of the three terms on the right side of the particle update rate formula are: retaining the original motion inertia of the particles, which helps to explore new areas; encouraging particles to move closer to their historical optimal positions; and pushing particles to move closer to the global optimal position of the group.

[0066] First, randomly initialize the position of the particles and speed , ensuring that the objective function value corresponding to the position of each particle is calculated within the specified range , set the initial position of each particle to its individual historical optimal position , and find the global optimal position in the group In each iteration, the velocity of each particle is updated and location , calculate the objective function value of the new position of the particle , update the individual optimal position :like , then suppose , update the global optimal position :like , then suppose When the maximum number of iterations is reached or the objective function value meets the preset threshold, the algorithm terminates and outputs the global optimal position and its corresponding objective function value , That is to optimize the economic net present value.

[0067] In order to demonstrate the feasibility and superiority of the present application, the following examples are given.

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

[0069] Similarly, using the method of the present invention, based on Figure 7 The sedimentary phase map shown in the figure is obtained by combining edge detection with the level set function to obtain a sedimentary phase binary map. Figure 8 The connected network model shown in FIG. 1 and the connected network diagram and the model after being mapped to the 2D Cartesian coordinate system are shown in FIG. Figure 8 The ES-MDA algorithm is used for historical matching. The population size is set to 150 and the number of iterations is 4. The production well 1 is used as an example to demonstrate the model effect. Its prior distribution is as follows: Fig. 9 As shown in the figure, the dotted line divides the historical fitting period from the prediction period, the black line represents the observation results, and the gray line represents the simulation results of the prior model; its posterior distribution is shown in Fig.10 As shown in the figure, the black line represents the observation results, and the gray line represents the posterior model simulation results; it can be seen that after historical fitting, the historical fitting results and prediction results of the model are good, which verifies the ability of the model in water drive and polymer drive simulation.

[0070] Finally, the particle swarm optimization algorithm was used to optimize the injection rate and liquid production. 100 days was used as an optimization time step, and the optimization was carried out for ten steps, for a total of 1000 days. The optimal optimization variables are as follows: Fig.11 As shown in Figure 2, the optimized cumulative oil production is Fig.12 As shown, it can be seen that after optimization, significantly higher cumulative oil production can be obtained than before optimization, thus verifying the reliability of the method proposed in this application in considering sedimentary phase modeling and optimization process.

[0071] The above implementation modes are only used to illustrate the present application, and are not intended to limit the present application. Ordinary technicians in the relevant technical field may 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, and the scope of patent protection of the present application shall be limited by the claims.

Claims

1. A rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints, characterized in that: The method comprises: Identifying sedimentary facies boundaries according to the distribution characteristics of sedimentary facies, and establishing a connectivity network model taking sedimentary facies constraints into account based on the sedimentary facies boundaries, well locations and attributes of each sedimentary facies; Dividing the connected network model into a one-dimensional grid and converting it into a grid connection format of a general simulator, solving it using a general simulator, and obtaining pressure and saturation solutions; Based on sedimentary facies and reservoir constraints, an integrated smooth multi-data assimilation algorithm is used to establish an automatic history matching mathematical model to achieve updated inversion of key reservoir parameters; wherein the key reservoir parameters include grid size, grid permeability and relative permeability curve; Taking the economic net present value as the indicator, the water injection volume, liquid production volume and polymer injection volume of each well as the optimization variables, the differential evolution algorithm was used to establish a mathematical model for reservoir injection and production optimization by constraining the total polymer injection volume and reservoir pressure.

2. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 1, characterized in that: The sedimentary facies boundaries are identified according to the distribution characteristics of the sedimentary facies, and based on the sedimentary facies boundaries, well locations and attributes of each sedimentary facies, a connected network model considering sedimentary facies constraints is established, including: According to the sedimentary facies distribution of the block, the phase boundaries in the sedimentary facies distribution map are extracted using edge detection method; The phase boundary in the sedimentary phase distribution map is used as the initialization boundary of the level set function, and the regional evolution is performed by the level set method to obtain a binary map; According to the actual geological information of each well point, the initial connected network model is generated according to the angle and distance judgment criteria, and the wells are connected by a one-dimensional connected unit.

3. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 2, characterized in that: According to the sedimentary facies distribution of the block, the phase boundaries in the sedimentary facies distribution map are extracted using edge detection methods, including: The obtained sedimentary phase map is converted into a grayscale image, and a Gaussian function is used to perform weighted averaging according to the grayscale values ​​of the pixel to be filtered and its neighborhood to filter out the high-frequency noise superimposed in the image. The Gaussian function is as follows: (1) In the formula, is the standard deviation of the controlling Gaussian filter; f ( x , y ) is a Gaussian function; e is a natural constant; Using the Sobel operator, the magnitude and direction of the image gradient are calculated using the following formulas (2), (3) and (4): (2) (3) (4) In the formula, are the gradients of the image in the x and y directions, is the original image, It is a Sobel filter, which calculates the gradients in the horizontal and vertical directions respectively, and the total gradient of each pixel G The calculation formula is: (5) Gradient direction The expression is: (6) For each pixel, check whether its gradient amplitude is the local maximum along the gradient direction, accurately locate the edge, connect the scattered edge points into a complete boundary, and obtain a binary image, which is used as the phase boundary in the sedimentary phase distribution map.

4. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 2, characterized in that: The phase boundary in the sedimentary phase distribution map is used as the initialization boundary of the level set function, and the regional evolution is performed by the level set method to obtain a binary map, including: Select the signed Euclidean distance from the point on the plane to the contour curve as the level set function : (7) In the formula, Represents a point on an image To level set contour curve distance; The time-varying level set function is expressed as , the evolution of the surface is described by the level set evolution equation as follows: (8) In the formula, is the velocity function, which determines the evolution direction and speed of the surface. represents the gradient magnitude of the level set function, t Indicates actual time; Before evolving the level set, the level set function is converted into the SDF function by the following formula: (9) (10) In the formula, is a reinitialized virtual time variable, which has nothing to do with the actual time t. is the initial level set function, is a sign function, used to keep the positive and negative signs unchanged, is a constant used to avoid the denominator being 0; Level set evolution is performed based on the level set evolution equation to obtain a binary map, and the binary map is used to divide regions corresponding to multiple depositions.

5. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 1, characterized in that: The connected network model is divided into one-dimensional grids and converted into a grid connection format of a general simulator. The general simulator is used to solve the model to obtain pressure and saturation solutions, including: 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; the one-dimensional connected unit of the connected network model is converted into a grid connection format of a universal simulator, and the universal simulator is used to solve the problem to obtain pressure and saturation solutions; On the basis of the initial well network segmentation map, 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 corresponding to different deposits 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; the connected network model is mapped to 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 grid corresponding to the position of different connected 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 mapped to the 2D Cartesian coordinate system is brought into the universal simulator to quickly obtain pressure and saturation solutions.

6. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 1, characterized in that: Based on sedimentary facies and reservoir constraints, an integrated smooth multi-data assimilation algorithm is used to establish an automatic history matching mathematical model to achieve updated inversion of key reservoir parameters, including: Construct the history matching objective function and determine the constraints in combination with the sedimentary facies and reservoirs; Based on the integrated smooth multi-data assimilation algorithm, an automatic history fitting mathematical model is established to realize the two-step parameter inversion process.

7. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 6, characterized in that: Construct the history matching objective function and determine the constraints based on sedimentary facies and reservoirs, including: Construct the history matching objective function, expressed as: (11) In the formula, are adjustable parameters of the model, is the objective function; Simulation results for the connected network model; represents the actual observed data; Represents the covariance matrix of the observation error; min means taking the minimum value; T is the matrix transposition operator; The model can adjust the parameters It is expressed as: (12) In the formula, is the grid volume, is the grid permeability, is the depth of oil-water contact surface; The grid attributes and reservoir conditions in each phase are constrained by the following formula: (13) In the formula, For Sedimentary Grid permeability; For Sedimentary Grid permeability; and They are respectively minimum and maximum permeabilities of sedimentary facies; and They are respectively minimum and maximum permeabilities of sedimentary facies; and They are respectively Sedimentary facies and Sedimentary Mesh volume; and They are respectively Sedimentary facies and Total reservoir volume of sedimentary facies; and Respectively indicate the Sedimentary facies and The number of grids for sedimentary facies; and Respectively represent the lowest and highest depths of the oil-water contact surface.

8. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 7, characterized in that: 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 the block liquid production, the block oil production, the block pressure and / or the block water injection volume, and the single well production data includes the single well oil production and / or the single well water content.

9. The rapid simulation optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 1, characterized in that: Taking the economic net present value as the indicator, the water injection volume, liquid production volume and polymer injection volume of each well as the optimization variables, the differential evolution algorithm is used to constrain the total polymer injection volume and reservoir pressure to establish a reservoir injection and 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 the optimization variables, an optimal control mathematical model is established; The total amount of polymer injection and reservoir pressure are used as constraints to constrain the optimization variables; The particle swarm algorithm is used to solve the production optimization problem and find the optimal optimization variables of the oil field that can maximize the economic net present value.

10. The rapid simulation and optimization method for polymer flooding reservoirs considering sedimentary phase constraints according to claim 9, characterized in that: The objective function of the optimal control mathematical model is expressed as: (14) In the formula, max means taking the maximum value, is the objective function value; are optimization variables, including water injection rate, polymer injection rate, and liquid production rate; is the total number of control steps; They are the number of production wells, water injection wells, and polymer injection wells; represent 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 number of Average injection rate of injection wells; The average polymer injection rate of the kth polymer injection well in the nth simulation time step; b is the annual discount rate; The unit time step is the nth simulation time step; is the cumulative time up to the n simulation time step.

Images