Reservoir remaining oil prediction method based on physical information neural network
By constructing an equidistant three-dimensional convolutional neural network model and combining discrete Ritz energy terms with single-point implicit algebraic polishing, the spatial aliasing and mass drift problems in the prediction of remaining oil in strongly heterogeneous reservoirs were solved, achieving high-precision prediction of bottom hole flowing pressure and providing a basis for the exploration and development of heterogeneous reservoirs.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- QINGDAO UNIV OF TECH
- Filing Date
- 2026-05-25
- Publication Date
- 2026-07-31
AI Technical Summary
When applied to predicting remaining oil in highly heterogeneous reservoirs, existing neural network models are prone to problems such as spatial aliasing, mass drift, and inaccurate bottomhole pressure response, making it difficult to meet the requirements of mass conservation, shock wave stability, and accurate local response.
A three-dimensional convolutional neural network model with equidistant voids was constructed. The deep features of the reservoir model were extracted by combining the symmetric expansion scheduling strategy. The discrete Ritz energy term and single-point implicit algebraic polishing were integrated. An explicit hyperbolic divergence anti-leakage protection layer was introduced. The training model was optimized to predict the dimensionless pressure drop and local Darcy flux, so as to achieve accurate prediction of the distribution of remaining oil in the whole area.
It effectively mitigates the optimization gradient explosion caused by multi-scale heterogeneity and well-control mutation, improves the geometric characteristics of the reservoir model optimization space, avoids spatial resampling errors, enhances the accuracy of bottom hole flowing pressure prediction, and achieves accurate prediction of the remaining oil distribution in heterogeneous reservoirs.
Smart Images
Figure CN122333419B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of oil and gas reservoir development technology, and specifically to a method for predicting remaining oil in oil reservoirs based on physical information neural networks. Background Technology
[0002] Currently, multiphase flow in porous media is a core physical process in petroleum engineering applications such as enhanced oil recovery. Accurately characterizing the evolution of the saturation field across the entire region and tracking the fluid front position are the theoretical foundation for precise prediction of remaining oil distribution. In numerical simulations of industrial-scale reservoirs, the implicit pressure-explicit saturation scheme decouples the pressure and saturation fields in time, maintaining a sharp displacement front during convection-dominated transport, thus clearly defining the boundaries of the enriched areas of remaining oil. However, when dealing with strongly heterogeneous real-world reservoirs, its computational efficiency is severely limited by the iterative solution of the ill-conditioned elliptic pressure equation at each time step.
[0003] In recent years, methods such as Physical Information Neural Networks (PINNs) have provided new avenues for constructing continuous, differentiable surrogate solvers. These methods can significantly accelerate forward modeling processes and support gradient-based automatic history fitting and closed-loop parameter inversion. However, when directly applied to predict the remaining oil distribution in industrial-scale hyperbolic-elliptic reservoir systems, they encounter severe physical inconsistencies, specifically facing three main technical bottlenecks: First, spatial aliasing caused by resampling contaminates the pressure gradient and disrupts the stability of upwind transport, thereby completely destroying the hyperbolic saturation shock wave structure controlling the remaining oil boundary; Second, multi-scale heterogeneity and discontinuous well control scheduling generate severe scale imbalances during optimization, leading to severe residual fluctuations and gradient explosions in the network during optimization training, making it difficult to accurately characterize the local real physical response near the well point, thus affecting the dynamic simulation of remaining oil around a single well; Third, pressure approximation errors may accumulate into non-physical mass drift during explicit saturation updates, resulting in severe distortion in the prediction of well water cut and the overall remaining oil distribution.
[0004] Therefore, there is an urgent need for a reservoir residual oil prediction method based on physical information neural networks, which can meet the requirements of mass conservation, shock wave stability and accurate local response while maintaining end-to-end differentiability, and predict the distribution of residual oil in the reservoir with high accuracy. Summary of the Invention
[0005] This invention aims to solve the above-mentioned problems and proposes a reservoir residual oil prediction method based on physical information neural network. It solves the problems that existing neural network models are prone to spatial aliasing, mass drift and inaccurate bottom hole pressure response when applied to the prediction of residual oil in strongly heterogeneous reservoirs. It realizes accurate prediction of residual oil distribution in heterogeneous reservoirs and provides a basis for the exploration and development of heterogeneous reservoirs.
[0006] The present invention adopts the following technical solution: A method for predicting remaining oil in a reservoir based on a physical information neural network includes the following steps: Step 1: Construct a reservoir model based on the basic physical parameters and initial state of the reservoir to be predicted, calculate the initial residual field of the reservoir model, normalize the initial residual field and the initial permeability field, and then splice them with the initial saturation field and the initial total mobility field of the reservoir model to construct a dimensionless input dataset. Step 2: Construct an equidistant three-dimensional convolutional neural network model with equidistant holes. Combine the symmetric dilation scheduling strategy to extract deep features of the reservoir model while maintaining the topological consistency of the physical grid space of the reservoir model. Use the equidistant three-dimensional convolutional neural network model with equidistant holes to predict the dimensionless pressure drop. Step 3: Construct a discrete Ritz energy term by integrating the conductivity between adjacent grids of the reservoir model. Combine the global total residual and local well grid residual of the reservoir model to determine the loss function of the equidistant hollow 3D convolutional neural network model. Perform physical constraints and optimization training on the equidistant hollow 3D convolutional neural network model. Step 4: Perform single-point implicit algebraic polishing on the grid containing the well network in the reservoir model, reconstruct the Dirac singularity at the bottom of the well, and then perform local correction on the pressure field of the reservoir model to obtain the real physical pressure field of the reservoir model. Step 5: Calculate the local Darcy flux using the real physical pressure field of the modified reservoir model, and introduce an explicit hyperbolic divergence leak-proof protection layer in the explicit saturation update for transport update, obtain the whole-area water saturation of the reservoir model, and obtain the prediction result of the whole-area remaining oil distribution of the reservoir to be predicted.
[0007] Preferably, step 1 includes the following steps: Step 1.1: Obtain the basic physical parameters and initial state of the reservoir to be predicted and construct the reservoir model. The basic physical parameters include the reservoir's permeability field, porosity field, total mobility field and relative permeability curve. The initial state of the reservoir includes the reservoir's initial water saturation field and initial pressure field. Step 1.2: Calculate the original algebraic residuals based on the initial pressure field of the reservoir to obtain the initial residual field of the reservoir model; The formula for calculating the primitive algebraic residual is as follows: ; In the formula, The primitive algebraic residual; For source terms; It is a discrete Laplace matrix; Total flowability; This represents the initial pressure field; The initial residual field and initial permeability field of the reservoir model are normalized to obtain: , ; In the formula, The residuals after normalization; The maximum value of the residual. ,in, It is a function for maximizing the value; It is a function for maximizing the value; For rule domains; This is the error term; This represents the normalized penetration rate. This represents the initial permeability field; This represents the maximum penetration rate. ; Step 1.3 involves stitching together the initial residual field, initial saturation field, initial total mobility field, and normalized residual field of the reservoir model to form multiple channels with 4 channels and dimensions [missing information]. The input tensor is used to construct the input dataset. , ,in, The initial water saturation level, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grid cells in each direction.
[0008] Preferably, in step 2, the equidistant dilated 3D convolutional neural network model includes a data preprocessing module, a core network module, and an output processing module. The data preprocessing module normalizes the input physical field parameters and then concatenates the data to form an input tensor. It includes an input layer, a normalization layer, and a concatenation layer. In the data preprocessing module, the input layer is used to input the physical field, the normalization layer is used to normalize the input physical field, and the concatenation layer is used to concatenate the normalized physical field to form the input tensor. The core network module extracts deep features from the input tensor and includes an input convolutional layer, multiple dilated residual blocks, and an output convolutional layer. The input convolutional layer is used to input the input tensor obtained from the data preprocessing module. The dilated residual block includes an input layer, a dilated convolutional layer, an activation function layer, a regular convolutional layer, and an output layer, with residual connections between the input and output layers within the dilated residual block. The output processing module predicts the dimensionless pressure drop based on the extracted deep features.
[0009] Preferably, in step 2, the input physical field is input to the data preprocessing module. After normalization and data concatenation, an input tensor is constructed to form the input dataset. This input tensor is then input to the core network module. The input convolutional layer in the core network module expands the number of channels in the input tensor from 4 to 32. Edge padding is then performed on the input tensor using a copying mode, and the tensor is input to the core network module. Using continuous, non-pooling and non-upsampling dilated residual blocks in the core network module, deep features are extracted from the input tensor according to a symmetric dilation scheduling strategy based on the dilation rate. These features are then input to the output processing module. The output processing module performs dimensionality reduction on the 32-channel deep features, outputting a predicted value of the dimensionless pressure drop. The dimensionless pressure drop prediction value The number of channels is 1, and the dimension is ,in, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grid cells in each direction.
[0010] Preferably, in step 3, in order to avoid spurious convergence caused by the loss of pure discrete residuals, a discrete Ritz energy term is constructed by introducing the conductivity between adjacent grids of the reservoir model, which is used to force the output of the reservoir model network to satisfy the physical space smoothness. The expression for the discrete Ritz energy term is: ; In the formula, For discrete Ritz energy terms; , All are grid numbers; Adjacent grids in the reservoir model With grid The conductivity between them; Mesh in reservoir model The predicted pressure correction value after dimensional restoration; Mesh in reservoir model The predicted pressure correction value after dimensional restoration; This represents the algebraic residual corresponding to the initial pressure field; This represents the total number of grid cells in the reservoir model; The total residual and local residual of the reservoir model are calculated and combined with the discrete Ritz energy term to determine the loss function of the equidistant three-dimensional convolutional neural network model. The loss function of the equidistant three-dimensional convolutional neural network model is then used to perform backpropagation and optimization training on the equidistant three-dimensional convolutional neural network model. The loss function of the equidistant, dilated 3D convolutional neural network model is a normalized hybrid Ritz-residual loss function used to maintain dimensional balance, and the calculation formula is as follows: ; in, ; ; In the formula, The normalized mixed Ritz-residual loss function; This represents the maximum residual value. The average conductivity rate for the entire region; It is a diagonal matrix of discrete operator matrices. ,in, It is a diagonal matrix. It is a discrete Laplace matrix; The penalty weighting coefficient is used to control the global algebraic consistency residual; The penalty weighting coefficient for the local response residuals of the well; It is a 2-norm; The total residual for the entire region; For source terms; The current forecast is for total pressure; As the reference pressure field; This is the pressure correction value for the reservoir grid.
[0011] Preferably, in step 4, during the optimization training of the equidistant dilated 3D convolutional neural network model, when the model is in the inference stage where it has stopped backpropagation, the output tensor of the model is extracted and its dimensions are restored to obtain the pressure field. Combined with relaxation vectors Set the reservoir model's regular domain mesh, set the standard relaxation factor, and include source terms in the reservoir model. A strong relaxation factor is set at the well node mesh to enforce strict algebraic coverage; Based on the single-point implicit algebraic polishing formula, the pressure field is analyzed using discrete linear operators. Local exact inversion and algebraic correction are performed to reconstruct the Dirac singularity point at the bottom grid of the reservoir model. This is used to eliminate the pressure response hysteresis when the fluid mobility changes suddenly inside the reservoir model, and to obtain the real physical pressure field of the reservoir model. The single-point implicit algebraic polishing formula is as follows: ; In the formula, The physical pressure field after polishing correction; Let it be a relaxation vector; This indicates element-wise multiplication. It is a diagonal matrix; For source terms; It is a discrete Laplace matrix.
[0012] Preferably, in step 5, the residual elliptical approximation error in the equidistant three-dimensional convolutional neural network model is obtained, and the divergence defect at each grid of the reservoir model is determined as follows: ; In the formula, For grid number; For grid The divergence defect; For grid The source item; It is a discrete Laplace matrix; The current forecast is for total pressure; The local Darcy flux was calculated using the real physical pressure field of the modified reservoir model, and an explicit hyperbolic divergence leak-proof protection layer operator was embedded in the explicit finite volume substep cycle to update the water saturation at all grids of the reservoir model. The update formula for the water saturation of the reservoir model grid is as follows: ; In the formula, for Time Grid Water saturation; for Time Grid Water saturation; For grid porosity; For grid Volume; The step size for explicit time substeps; For grid Net windward flux; For grid The collection of items; It is a fractional flow function; For limiting functions; This represents the largest source term amplitude in the entire region.
[0013] An autoregressive training method was used to train an equidistant, three-dimensional convolutional neural network model. Based on the water saturation at each grid of the reservoir model, the remaining oil at all grids of the reservoir model was obtained, and the prediction results of the distribution of remaining oil in the entire reservoir at each time step were obtained.
[0014] The present invention has the following beneficial effects: This invention proposes a reservoir residual oil prediction method based on a physical information neural network. By constructing an equidistant, three-dimensional convolutional neural network model, it effectively alleviates the optimization gradient explosion caused by multi-scale heterogeneity and well-control mutations, improves the geometric characteristics of the reservoir model's optimization space, and the equidistant cavity topology used in the equidistant, three-dimensional convolutional neural network model maintains its geometric consistency with the physical discrete grid, expands the global receptive field, and avoids aliasing errors caused by spatial resampling. By introducing local single-point implicit algebraic polishing, it compensates for the severe local gradients in the neural network model's prediction process. This significantly improves the prediction accuracy of bottom hole flowing pressure by addressing the shortcomings of traditional methods. Furthermore, by introducing an embedded explicit hyperbolic divergence anti-leakage protection layer in the explicit saturation transmission step, the artificial mass accumulation path of small pressure residuals over time is directly cut off, achieving macroscopic zero mass drift. This enhances the accurate prediction of the dynamic evolution of residual oil in long-term water-drive operations in highly heterogeneous reservoirs. It also solves the problems of spatial aliasing, mass drift, and inaccurate bottom hole pressure response that easily occur when existing neural network models are applied to predict residual oil in highly heterogeneous reservoirs. This enables accurate prediction of the distribution of residual oil in heterogeneous reservoirs, providing a basis for the exploration and development of heterogeneous reservoirs. Attached Figure Description
[0015] Figure 1 This is a schematic diagram of a reservoir residual oil prediction method based on a physical information neural network according to the present invention.
[0016] Figure 2 This is a permeability parameter distribution map of the intermediate layer in a three-dimensional reservoir according to the present invention. In the map, white triangles represent water injection wells, and black circles represent production wells.
[0017] Figure 3 This is a schematic diagram of the isometric, dilated three-dimensional convolutional neural network model of the present invention.
[0018] Figure 4 This is a schematic diagram of the structure of the void residual block in the equidistant void 3D convolutional neural network model of the present invention.
[0019] Figure 5The figure shows the prediction of the remaining oil distribution using the method of the present invention. In the figure, (a) is the remaining oil distribution of the heterogeneous reservoir predicted by the method of the present invention on day 350, (b) is the remaining oil distribution of the heterogeneous reservoir predicted by the method of the present invention on day 660, and (c) is the remaining oil distribution of the heterogeneous reservoir predicted by the method of the present invention on day 1000.
[0020] Figure 6 The figures show the predicted distribution of remaining oil in the heterogeneous reservoir using the finite element volumetric method. In the figure, (a) shows the distribution of remaining oil in the heterogeneous reservoir on day 350, predicted using the finite element volumetric method; (b) shows the distribution of remaining oil in the heterogeneous reservoir on day 660, predicted using the finite element volumetric method; and (c) shows the distribution of remaining oil in the heterogeneous reservoir on day 1000, predicted using the finite element volumetric method.
[0021] Figure 7 The figure shows the absolute error of the method of the present invention. In the figure, (a) is the absolute error of the method of the present invention for predicting the distribution of remaining oil in a heterogeneous reservoir on day 350, (b) is the absolute error of the method of the present invention for predicting the distribution of remaining oil in a heterogeneous reservoir on day 660, and (c) is the absolute error of the method of the present invention for predicting the distribution of remaining oil in a heterogeneous reservoir on day 1000.
[0022] Figure 8 This is a schematic diagram of the water cut variation curve of the production well throughout the entire cycle using the method of the present invention. In the diagram, the red dashed line represents the limit of the production regime change, the blue solid line represents the water cut variation curve of the production well throughout the entire cycle obtained by the finite element volume method, and the orange dashed line represents the water cut variation curve of the production well throughout the entire cycle obtained by the method of the present invention.
[0023] Figure 9 This is a graph showing the change in bottom hole flowing pressure curve throughout the entire cycle during the prediction process using the method of this invention. In the graph, the red dashed line represents the boundary of the production system change, the gray solid line represents the relative bottom hole flowing pressure change curve throughout the entire cycle obtained using the finite element volume method, the black dashed line represents the relative bottom hole flowing pressure change curve throughout the entire cycle obtained using the method of this invention, the blue solid line represents the relative bottom hole flowing pressure change curve throughout the entire cycle of the injection well, and the black solid line represents the relative bottom hole flowing pressure change curve throughout the entire cycle of the production well. Detailed Implementation
[0024] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings.
[0025] This invention proposes a reservoir residual oil prediction method based on physical information neural networks, such as... Figure 1 As shown, the method for accurately predicting the remaining oil distribution in three-dimensional heterogeneous reservoirs includes the following sub-steps: Step 1: Construct a reservoir model based on the fundamental physical parameters and initial state of the reservoir to be predicted. Calculate the initial residual field of the reservoir model. After normalizing the initial residual field and the initial permeability field, concatenate them with the initial saturation field and the initial total mobility field of the reservoir model to construct a dimensionless input dataset. This specifically includes the following steps: Step 1.1: Obtain the basic physical parameters and initial state of the reservoir to be predicted. The basic physical parameters include the reservoir's permeability field, porosity field, total mobility field, and relative permeability curve. The initial state of the reservoir includes the initial water saturation field and initial pressure field. Construct a reservoir model based on the basic physical parameters and initial state of the reservoir to be predicted.
[0026] In this embodiment, a reservoir model is constructed in three-dimensional space based on the basic physical parameters and initial state of the reservoir to be predicted. The reservoir model includes four water injection wells and one production well. The water injection wells are located near the four corners of the physical domain of the reservoir model, and the production well is located at the center of the physical domain. The porosity of the reservoir model is set to 0.2, and the net-to-gross ratio is set to 1. , , The total lengths in the directions are 500m, 500m, and 20m respectively. The reservoir model is then meshed to determine its position within the model. Total number of grids in each direction For 31, in Total number of grids in each direction For 31, in Total number of grids in each direction The number is 3, with a total of 2883 grids.
[0027] In this embodiment, the reservoir model is located at the middle position along... Permeability distribution in the directional plane is as follows Figure 2 As shown in the figure. In this embodiment, the production wells operate on a fixed production rate, with an initial production rate of 400 cubic meters per day, increasing to 500 cubic meters per day on day 300. The water injection wells operate on a fixed injection rate, with an initial injection rate of 100 cubic meters per day. On day 300, the injection rate of the two wells on the right side is changed to 50 cubic meters per day, and the injection rate of the two wells on the left side is changed to 200 cubic meters per day. By day 600, the injection rate of the two wells on the right side is changed to 200 cubic meters per day, and the injection rate of the two wells on the left side is changed to 50 cubic meters per day. The reservoir model considers two phases of oil and water, with both rock and fluid being incompressible. The production time is set to 1000 days, divided into 200 time steps for remaining oil prediction.
[0028] Step 1.2: Calculate the original algebraic residuals based on the initial pressure field of the reservoir to obtain the initial residual field of the reservoir model.
[0029] Specifically, the formula for calculating the original algebraic residual is as follows: ; In the formula, The primitive algebraic residual; For source terms; It is a discrete Laplace matrix; Total flowability; This represents the initial pressure field.
[0030] Directly inputting unnormalized residuals leads to severe ill-conditioned Hessian matrix and triggers gradient explosion during neural backpropagation. Normalizing the initial residual field and initial permeability field of the reservoir model yields: , ; In the formula, The residuals after normalization; The maximum value of the residual. ,in, It is a function for maximizing the value; It is a function for maximizing the value; This is the rule domain, i.e., the sample space; This is the error term; This represents the normalized penetration rate. This represents the initial permeability field; This represents the maximum penetration rate. .
[0031] Step 1.3 involves stitching together the initial residual field, initial saturation field, initial total mobility field, and normalized residual field of the reservoir model to form multiple channels with 4 channels and dimensions [missing information]. The input tensor is used to construct the input dataset. , ,in, The initial water saturation level, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grid cells in each direction.
[0032] Step 2: Construct an equidistant, three-dimensional convolutional neural network model with equidistant holes. Combine this with a symmetric expansion scheduling strategy to extract deep features of the reservoir model while maintaining the topological consistency of the physical grid space. Then, use the equidistant, three-dimensional convolutional neural network model with equidistant holes to predict the dimensionless pressure drop.
[0033] In this embodiment, an equidistant, dilated 3D convolutional neural network model is constructed, such as... Figure 3 As shown, the equidistant dilated 3D convolutional neural network model includes a data preprocessing module, a core network module, and an output processing module. The data preprocessing module normalizes the input physical field parameters and then concatenates the data to form an input tensor. It includes an input layer, a normalization layer, and a concatenation layer. In the data preprocessing module, the input layer is used to input the physical field, the normalization layer is used to normalize the input physical field, and the concatenation layer is used to concatenate the normalized physical field to form the input tensor. The core network module extracts deep features from the input tensor and includes an input convolutional layer, multiple dilated residual blocks, and an output convolutional layer, such as... Figure 4 As shown, the input convolutional layer is used to input the input tensor obtained by the input data preprocessing module, the dilated residual block includes an input layer, a dilated convolutional layer, an activation function layer, a normal convolutional layer and an output layer, and there is a residual connection between the input layer and the output layer in the dilated residual block; the output processing module is used to predict the dimensionless pressure drop of the output based on the extracted deep features.
[0034] The input physical field is fed into the data preprocessing module, which normalizes the input physical field and concatenates the data to form an input tensor. After constructing the input dataset, the input tensor in the input dataset is fed into the core network module. The input convolutional layer in the core network module expands the number of channels of the input tensor from 4 to 32. Then, the edge padding of the input tensor is performed using a copy mode to ensure that the size of the output feature map remains unchanged.
[0035] The output feature map of the data preprocessing module is then input into the core network module. Using the continuous dilated residual blocks in the core network module that do not contain any pooling or upsampling operations, deep features are extracted from the input tensor according to the dilation rate following a symmetric dilation scheduling strategy. During the deep feature extraction process, the dilated convolutional layers in the dilated residual blocks are used to fill holes of the appropriate size according to the current dilation rate to expand the receptive field, and the holes are added to the residuals of the input tensor before being input into the output processing module.
[0036] The output processing module performs dimensionality reduction on the 32-channel deep features, outputting the predicted value of the dimensionless pressure drop. The dimensionless pressure drop prediction value The number of channels is 1, and the dimension is ,in, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grid cells in each direction.
[0037] Step 3: Construct a Discrete Ritz Energy Term by fusing the conductivity between adjacent grids in the reservoir model. Combine the global total residual and local well grid residual of the reservoir model to determine the loss function of the equidistant hollow 3D convolutional neural network model, which is a normalized hybrid Ritz-residual loss function used to maintain dimensional balance. Then, perform physical constraints and optimization training on the equidistant hollow 3D convolutional neural network model.
[0038] In this embodiment, to avoid pseudo-convergence caused by pure discrete residual loss, a discrete Ritz energy term is constructed by introducing the conductivity between adjacent grids of the reservoir model, which is used to force the output of the reservoir model network to satisfy physical space smoothness.
[0039] The expression for the discrete Ritz energy term is: ; In the formula, For discrete Ritz energy terms; , All are grid numbers; Adjacent grids in the reservoir model With grid The conductivity between them; Mesh in reservoir model The predicted pressure correction value after dimensional restoration; Mesh in reservoir model The predicted pressure correction value after dimensional restoration; This represents the algebraic residual corresponding to the initial pressure field; This represents the total number of grids in the reservoir model.
[0040] The total and local residuals of the reservoir model were calculated, and combined with the discrete Ritz energy term, the loss function of the equidistant, three-dimensional convolutional neural network model was determined as follows: ; in, ; ; In the formula, The normalized mixed Ritz-residual loss function; This represents the maximum residual value. The average conductivity rate for the entire region; It is a diagonal matrix of discrete operator matrices. ,in, It is a diagonal matrix. It is a discrete Laplace matrix; The penalty weighting coefficient is used to control the global algebraic consistency residual; The penalty weighting coefficient for the local response residuals of the well; It is a 2-norm; The total residual for the entire region; For source terms; The current forecast is for total pressure; As the reference pressure field; This is the pressure correction value for the reservoir grid.
[0041] We use the loss function of the equidistant dilated 3D convolutional neural network model to perform backpropagation and optimization training on the equidistant dilated 3D convolutional neural network model.
[0042] Step 4: Perform single-point implicit algebraic polishing on the grid containing the well network in the reservoir model, reconstruct the Dirac singularity at the bottom of the well, and then locally correct the pressure field of the reservoir model to obtain the real physical pressure field of the reservoir model.
[0043] In this embodiment, during the optimization training of the equidistant dilated 3D convolutional neural network model, when the model is in the inference stage where it has exited backpropagation, the output tensor of the model is extracted and its dimensions are restored to obtain the pressure field. Combined with position-dependent relaxation vectors Set the reservoir model's regular domain mesh, set the standard relaxation factor, and include source terms in the reservoir model. A strong relaxation factor with a value close to 1 is set at the well node grid to enforce strict algebraic coverage.
[0044] Based on the single-point implicit algebraic polishing formula, the pressure field is analyzed using discrete linear operators. Local precise inversion and algebraic correction are performed to reconstruct the Dirac singularity point at the bottom grid in the reservoir model. This is used to eliminate the pressure response hysteresis when fluid mobility changes occur inside the reservoir model, and to obtain the true physical pressure field of the reservoir model.
[0045] Specifically, the single-point implicit algebraic polishing formula is: ; In the formula, The physical pressure field after polishing correction; Let it be a relaxation vector; This indicates element-wise multiplication. It is a diagonal matrix; For source terms; It is a discrete Laplace matrix.
[0046] Step 5: Calculate the local Darcy flux using the real physical pressure field of the modified reservoir model, and introduce an explicit hyperbolic divergence leak-proof protection layer in the explicit saturation update for transport update, obtain the whole-area water saturation of the reservoir model, and obtain the prediction result of the whole-area remaining oil distribution of the reservoir to be predicted.
[0047] In this embodiment, the residual elliptical approximation error in the equidistant hollow 3D convolutional neural network model is obtained, and the divergence defect at each grid of the reservoir model is determined as follows: ; In the formula, For grid number; For grid The divergence defect; For grid The source item; It is a discrete Laplace matrix; This represents the current forecast of total pressure.
[0048] The local Darcy flux is calculated using the real physical pressure field of the modified reservoir model, and an explicit hyperbolic divergence leak-proof protection layer operator is embedded in the explicit finite volume substep cycle to update the water saturation at all grids of the reservoir model.
[0049] Specifically, the update formula for the water saturation of the reservoir model grid is as follows: ; In the formula, for Time Grid Water saturation; for Time Grid Water saturation; For grid porosity; For grid Volume; The step size for explicit time substeps; For grid Net windward flux; For grid The collection of items; It is a fractional flow function; For limiting functions; This is the largest source term amplitude value in the entire region.
[0050] In this embodiment, the explicit hyperbolic divergence leak-proof protection layer operator is used to explicitly neutralize the artificial spurious phase volume caused by local pressure defects through algebraic correction, thereby cutting off the accumulation path of elliptic error in hyperbolic transport.
[0051] Finally, an autoregressive training was performed on the equidistant, three-dimensional convolutional neural network model. Based on the water saturation at each grid point of the reservoir model, the remaining oil content at each grid point was determined using the saturation conservation formula. The formula for calculating the remaining oil content is as follows: ,in, Remaining oil content, The water saturation level is used to obtain the remaining oil at all grid points in the reservoir model, thus obtaining the predicted distribution of the remaining oil in the entire reservoir at each time step.
[0052] To verify the prediction effect of the reservoir residual oil prediction method based on physical information neural network proposed in this invention, the method of this invention is compared with the finite volume method in the prior art. The residual oil distribution of the heterogeneous reservoir in this embodiment is predicted using both the method of this invention and the finite volume method. The prediction results of the residual oil distribution of the heterogeneous reservoir in this embodiment at day 350, day 660, and day 1000 are obtained. Days 350, 660, and 1000 correspond to the intermediate propagation period before breakthrough, after the first injection-production switch, and the later state after two switches take effect, respectively. Figure 5 and Figure 6 As shown; further, the absolute error in the prediction process of the method of the present invention is obtained, such as Figure 7 As shown, comparing the remaining oil distribution prediction results of the two methods at different times reveals that at day 350, the results of the two methods are highly consistent, with the high-saturation zone advancing towards the central production well in a leaf-like pattern, and the minimal error is concentrated only at the steep gradient leading edge. At day 660, although the first asymmetric switch changed the pressure topology, the method of this invention can still accurately capture the displacement branches guided by the flow channel and the leading edge reforming phenomenon. Although the absolute error has increased, it is still strictly limited to the narrow high-gradient zone and no global diffusion has occurred. At day 1000, the late-stage sweep morphology and leading edge position predicted by this method remain highly consistent with the reference solution obtained by the finite element volume method, and the error is always limited to the leading edge neighborhood, without producing any spatial fragmentation or non-physical artifacts.
[0053] Furthermore, the historical evolution of water cut in production wells was obtained using both the method of this invention and the finite volume method, as follows: Figure 8 As shown, after comparison, it was found that the prediction results of the method of the present invention and the prediction results of the finite volume method were highly consistent during the complete simulation period of 1000 days. That is, the equidistant three-dimensional convolutional neural network model used in the method of the present invention can not only accurately predict the timing of the initial water injection breakthrough, but also keenly and without delay capture the water cut fluctuation characteristics caused by the subsequent two injection and production scheme switching. This verifies that the equidistant three-dimensional convolutional neural network model constructed by the method of the present invention can simultaneously take into account the stability of numerical calculation and the strong constraint of physical conservation laws when dealing with long-term, large-step explicit integral evolution.
[0054] Furthermore, by comparing the relative bottomhole pressure change history of production wells and injection wells within the reservoir model during the remaining oil prediction process using the method of this invention, as follows: Figure 9 As shown, the comparison revealed significant reconstructing of the well response near the production regime switching points on days 300 and 600. Furthermore, due to the heterogeneity of permeability distribution, the responses of each injection well exhibited strong differences. Although bottomhole pressure is sensitive to local singular gradients around the well, the method of this invention closely matches the baseline solution of the finite volume method for both injection and production wells. Moreover, this method not only captures smooth long-term trends but also accurately captures sudden pressure redistribution caused by production regime switching. This indicates that even under rapidly changing reservoir conditions, its pressure reconstruction results maintain high accuracy.
[0055] In summary, the method of this invention can effectively alleviate the optimization gradient explosion caused by reservoir multi-scale heterogeneity and well-control abrupt changes, improve the prediction accuracy of the dynamic evolution of residual oil in long-term waterflooding in strongly heterogeneous reservoirs, and provide a basis for the exploration and development of heterogeneous reservoirs.
[0056] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.
Claims
1. A physical information neural network-based remaining oil prediction method for an oil reservoir, characterized by, Includes the following steps: Step 1: Construct a reservoir model based on the basic physical parameters and initial state of the reservoir to be predicted, calculate the initial residual field of the reservoir model, normalize the initial residual field and the initial permeability field, and then splice them with the initial saturation field and the initial total mobility field of the reservoir model to construct a dimensionless input dataset. Step 2: Construct an equidistant three-dimensional convolutional neural network model with equidistant holes. Combine the symmetric dilation scheduling strategy to extract deep features of the reservoir model while maintaining the topological consistency of the physical grid space of the reservoir model. Use the equidistant three-dimensional convolutional neural network model with equidistant holes to predict the dimensionless pressure drop. Step 3: Construct a discrete Ritz energy term by integrating the conductivity between adjacent grids of the reservoir model. Combine the global total residual and local well grid residual of the reservoir model to determine the loss function of the equidistant hollow 3D convolutional neural network model. Perform physical constraints and optimization training on the equidistant hollow 3D convolutional neural network model. Step 4: Perform single-point implicit algebraic polishing on the grid containing the well network in the reservoir model, reconstruct the Dirac singularity at the bottom of the well, and then perform local correction on the pressure field of the reservoir model to obtain the real physical pressure field of the reservoir model. Step 5: Calculate the local Darcy flux using the real physical pressure field of the modified reservoir model, and introduce an explicit hyperbolic divergence leak-proof protection layer in the explicit saturation update for transport update, obtain the whole area water saturation of the reservoir model, and obtain the whole area remaining oil distribution prediction results of the reservoir to be predicted. In step 3, to avoid pseudo-convergence caused by pure discrete residual loss, the conductivity between adjacent grids of the reservoir model is introduced to construct a discrete Ritz energy term, which is used to force the output of the reservoir model network to satisfy physical space smoothness. The expression for the discrete Ritz energy term is: ； In the formula, For discrete Ritz energy terms; , All are grid numbers; Adjacent grids in the reservoir model With grid The conductivity between them; Mesh in reservoir model The predicted pressure correction value after dimensional restoration; Mesh in reservoir model The predicted pressure correction value after dimensional restoration; This represents the algebraic residual corresponding to the initial pressure field; This represents the total number of grid cells in the reservoir model; The total residual and local residual of the reservoir model are calculated and combined with the discrete Ritz energy term to determine the loss function of the equidistant three-dimensional convolutional neural network model. The loss function of the equidistant three-dimensional convolutional neural network model is then used to perform backpropagation and optimization training on the equidistant three-dimensional convolutional neural network model. The loss function of the equidistant, dilated 3D convolutional neural network model is a normalized hybrid Ritz-residual loss function used to maintain dimensional balance, and the calculation formula is as follows: ； in, ； ； In the formula, The normalized mixed Ritz-residual loss function; This represents the maximum residual value. The average conductivity rate for the entire region; It is a diagonal matrix of discrete operator matrices. ,in, It is a diagonal matrix. It is a discrete Laplace matrix; The penalty weighting coefficient is used to control the global algebraic consistency residual; The penalty weighting coefficient for the local response residuals of the well; It is a 2-norm; The total residual for the entire region; For source terms; The current forecast is for total pressure; As the reference pressure field; This is the pressure correction value for the reservoir grid. 2.The physical information neural network-based reservoir remaining oil prediction method according to claim 1, characterized in that, Step 1 includes the following steps: Step 1.1: Obtain the basic physical parameters and initial state of the reservoir to be predicted and construct the reservoir model. The basic physical parameters include the reservoir's permeability field, porosity field, total mobility field and relative permeability curve. The initial state of the reservoir includes the reservoir's initial water saturation field and initial pressure field. Step 1.2: Calculate the original algebraic residuals based on the initial pressure field of the reservoir to obtain the initial residual field of the reservoir model; The formula for calculating the primitive algebraic residual is as follows: ； In the formula, The primitive algebraic residual; For source terms; It is a discrete Laplace matrix; Total flowability; This represents the initial pressure field; The initial residual field and initial permeability field of the reservoir model are normalized to obtain: ， ； In the formula, The residuals after normalization; The maximum value of the residual. ,in, It is a function for maximizing the value; It is a function for maximizing the value; For rule domains; This is the error term; This represents the normalized penetration rate. This represents the initial permeability field; This represents the maximum penetration rate. ; Step 1.3 involves stitching together the initial residual field, initial saturation field, initial total mobility field, and normalized residual field of the reservoir model to form multiple channels with 4 channels and dimensions [missing information]. The input tensor is used to construct the input dataset. , ,in, Initial water saturation For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grid cells in each direction. 3.The physical information neural network based reservoir remaining oil prediction method of claim 1, wherein, In step 2, the equidistant dilated 3D convolutional neural network model includes a data preprocessing module, a core network module, and an output processing module. The data preprocessing module normalizes the input physical field parameters and then concatenates the data to form an input tensor. It includes an input layer, a normalization layer, and a concatenation layer. The input layer is used to input the physical field, the normalization layer normalizes the input physical field, and the concatenation layer concatenates the normalized physical field to form the input tensor. The core network module extracts deep features from the input tensor and includes an input convolutional layer, multiple dilated residual blocks, and an output convolutional layer. The input convolutional layer is used to input the input tensor obtained from the data preprocessing module. The dilated residual block includes an input layer, a dilated convolutional layer, an activation function layer, a regular convolutional layer, and an output layer, with residual connections between the input and output layers within the dilated residual block. The output processing module predicts the dimensionless pressure drop based on the extracted deep features.
4. The reservoir residual oil prediction method based on physical information neural network according to claim 3, characterized in that, In step 2, the input physical field is input to the data preprocessing module. The data preprocessing module normalizes the input physical field and concatenates the data to form an input tensor. After constructing the input dataset, the input tensor from the input dataset is input to the core network module. The input convolutional layer in the core network module expands the number of channels of the input tensor from 4 to 32. Then, edge padding is performed on the input tensor using a copying mode, and it is input to the core network module. Using the continuous, non-pooling and non-upsampling dilated residual blocks in the core network module, deep features are extracted from the input tensor according to a symmetric dilation scheduling strategy based on the dilation rate. These features are then input to the output processing module. The output processing module performs dimensionality reduction on the 32-channel deep features and outputs the predicted value of the dimensionless pressure drop. The dimensionless pressure drop prediction value The number of channels is 1, and the dimension is ,in, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grids in the direction, For reservoir model in Total number of grid cells in each direction. 5.The physical information neural network based reservoir remaining oil prediction method of claim 1, wherein, In step 4, during the optimization training of the equidistant dilated 3D convolutional neural network model, when the model is in the inference stage where it has stopped backpropagation, the output tensor of the model is extracted and its dimensions are restored to obtain the pressure field. Combined with relaxation vectors Set the reservoir model's regular domain mesh, set the standard relaxation factor, and include source terms in the reservoir model. A strong relaxation factor is set at the well node mesh to enforce strict algebraic coverage; Based on the single-point implicit algebraic polishing formula, the pressure field is analyzed using discrete linear operators. Local exact inversion and algebraic correction are performed to reconstruct the Dirac singularity point at the bottom grid of the reservoir model. This is used to eliminate the pressure response hysteresis when the fluid mobility changes suddenly inside the reservoir model, and to obtain the real physical pressure field of the reservoir model. The single-point implicit algebraic polishing formula is as follows: ； In the formula, The physical pressure field after polishing correction; Let it be a relaxation vector; This indicates element-wise multiplication. It is a diagonal matrix; For source terms; It is a discrete Laplace matrix. 6.The physical information neural network based reservoir remaining oil prediction method of claim 1, wherein, In step 5, the residual elliptical approximation error in the equidistant three-dimensional convolutional neural network model is obtained, and the divergence defect at each grid of the reservoir model is determined as follows: ； In the formula, For grid number; For grid The divergence defect; For grid The source item; It is a discrete Laplace matrix; The current forecast is for total pressure; The local Darcy flux was calculated using the real physical pressure field of the modified reservoir model, and an explicit hyperbolic divergence leak-proof protection layer operator was embedded in the explicit finite volume substep cycle to update the water saturation at all grids of the reservoir model. The update formula for the water saturation of the reservoir model grid is as follows: ； In the formula, for Time Grid Water saturation; for Time Grid Water saturation; For grid porosity; For grid Volume; The step size for explicit time substeps; For grid Net windward flux; For grid The collection of items; It is a fractional flow function; For limiting functions; This represents the largest source term amplitude in the entire region. An autoregressive training method was used to train an equidistant, three-dimensional convolutional neural network model. Based on the water saturation at each grid of the reservoir model, the remaining oil at all grids of the reservoir model was obtained, and the prediction results of the distribution of remaining oil in the entire reservoir at each time step were obtained.