New method for constructing physical constraint oil reservoir agent model by coupling numerical simulator

Abstract

The invention discloses a novel method for constructing a physical constraint oil reservoir agent model by coupling a numerical simulator, and relates to the field of intelligent oil and gas field exploitation, and the method comprises the following steps: establishing an oil reservoir model, randomly disturbing injection-production parameters, and generating an injection-production parameter data set; the injection-production parameters are input into an oil reservoir model for simulation, state parameters are obtained, and the injection-production parameters and the state parameters form label data; building an R-U-Net neural network as an agent model; coupling a numerical simulator, and calculating mixing loss; and updating network parameters of the proxy model by using mixed loss back propagation to complete proxy model training. According to the method, the quality residual generated by discrete difference and the physical loss calculated by the Jacobian matrix are fused as constraint conditions, and the agent model training process is optimized; the precise modeling capability of a numerical simulator and the efficient prediction advantage of an agent model are fused, the problem that a traditional simulator consumes long calculation time in the linear solving process is effectively solved, and the limitation that the agent model is difficult to embed into a complex physical equation is made up.

Description

Technical Field

[0001] The present invention relates to the technical field of intelligent oil and gas field exploitation, and particularly to a new method for constructing a physically constrained reservoir surrogate model by coupling a numerical simulator. Background Art

[0002] In the construction of traditional reservoir simulation surrogate models, there are problems such as high data acquisition costs, limited computing resources, and the lack of physical meaning in the output of surrogate models. Existing solutions aim to reduce the dependence on training data and enhance the physical consistency of surrogate models. One is the Physics-Informed Neural Networks (PINN), which can construct a surrogate model that conforms to physical laws without a large amount of training data by directly embedding physical constraints (such as partial differential equations) into the network structure. However, when dealing with complex practical problems, PINN involves the numerical solution of partial differential equations in the loss function, resulting in a complex and costly calculation process and great implementation difficulty; another method is the Theory-Guided Neural Networks (TGNN), which embeds physical rules to constrain the training process of the surrogate model and combines physical loss with label loss, thereby improving the physical consistency and reliability of the surrogate model while ensuring the data fitting accuracy of the surrogate model.

[0003] Currently, there is a common problem in constructing a reservoir simulation surrogate model with physical meaning: for each complex reservoir flow problem, when constructing a physically constrained reservoir surrogate model, it is necessary to write complex control equations and constrain the training of the surrogate model through a difficult solution process. This method is theoretically feasible, but in practical applications, especially for complex reservoir models, it is generally difficult to promote due to the high complexity of physical principles and large computational difficulty. Current research mostly focuses on writing control equations for simple problems and fails to fully utilize the advantages of numerical simulators, resulting in insufficient applicability to complex reservoir systems.

[0004] In summary, existing technical solutions generally have problems such as high computational complexity, great implementation difficulty, and insufficient adaptability to complex reservoir models when constructing a physically constrained reservoir simulation surrogate model. These defects limit the wide application of physically constrained reservoir surrogate models in actual oil and gas field development. Therefore, a more efficient and practical solution is urgently needed. Summary of the Invention

[0005] To solve the problems of high computational complexity, difficulty in embedding physical equations, and insufficient adaptability to complex reservoir models in the existing methods for constructing reservoir simulation surrogate models, the present invention discloses a new method for constructing a physically constrained reservoir surrogate model by coupling a numerical simulator. This method incorporates the mass residuals and Jacobian matrices generated during the discrete difference process into the surrogate model training to calculate the physical loss as a constraint condition, thereby optimizing the training process of the surrogate model.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A new method for constructing a physically constrained reservoir surrogate model by coupling a numerical simulator, comprising the following steps:

[0008] S1. Establish a reservoir model: Use a numerical simulator to construct a three-dimensional oil-water two-phase seepage model, and perform spatial and temporal implicit discretization and solution of the governing equations by the finite volume method;

[0009] S2. Randomly perturb the injection-production parameters to generate an injection-production parameter dataset;

[0010] S3. Input the injection-production parameter dataset generated in step S2 into the reservoir model for simulation to obtain the corresponding reservoir state parameters. The injection-production parameters and state parameters form a training dataset, i.e., labeled data;

[0011] S4. Build an R-U-Net neural network as a surrogate model, including an input layer, an encoding module, a temporal module, a decoding module, and an output layer. The encoding module extracts spatial features through convolutional operations for downsampling, the temporal module uses long short-term memory networks to capture temporal dependencies, the decoding module restores the spatial dimension through upsampling, the input layer receives the permeability field and injection-production parameter field, and the output layer couples to predict the pressure field and saturation field;

[0012] S5. Couple the numerical simulator to calculate the hybrid loss;

[0013] S6. Use the hybrid loss for backpropagation to update the network parameters of the surrogate model and complete the training of the surrogate model.

[0014] Optionally, in step S1, the governing equation is:

[0015] ;

[0016] where x, y, and z represent three-dimensional spatial coordinates, the subscript l represents the fluid, represents the phase density, represents the phase viscosity, k represents the absolute permeability, represents the relative permeability, p represents the pressure, represents the phase saturation, represents the porosity, Indicates the flow rate of the fluid in the injection well or production well in phase l, calculated using the Peaceman equation: , where represents the thickness of the grid, represents the effective radius of the well, represents the absolute permeability at the well location, represents the actual radius of the well, represents the bottom-hole pressure;

[0017] After discrete difference, the control equation is simplified to:

[0018] ;

[0019] In the formula, g represents the mass residual, and A, F, and Q represent the cumulative phase, flowing phase, and source-sink phase respectively; and represent the reservoir states at times n + 1 and n respectively, including the distributions of saturation and pressure; represents the injection-production parameters at time n + 1. The numerical simulator can efficiently calculate the mass residual g based on the static reservoir model parameters, dynamic reservoir state parameters, and injection-production parameters of the reservoir.

[0020] Optionally, in step S2, for different problem perturbation variable parameters, taking the surrogate model for production optimization as an example, the injection-production parameters of the time series are randomly perturbed to obtain different datasets of injection-production parameters.

[0021] Optionally, in step S3, the injection-production parameters generated by random perturbation in step S2 are input into the reservoir model, and the reservoir state parameters (distributions of the pressure field and water saturation field) under different injection-production parameters are obtained through simulation. A set of injection-production parameters corresponds to a set of state parameters, jointly forming the training dataset, that is, the labeled data. The input features are selected as the permeability distribution and the injection-production parameters set by the sparse matrix. During the setting of the injection-production parameters, in the form of a sparse matrix, a positive well control value is set for the grid containing the injection well, and a negative well control value is set for the grid of the production well. The prediction result is the state distribution of pressure and water saturation at different injection-production parameters and different time steps. For the state parameters, water saturation and pressure are selected. Since the oil saturation can be obtained from the water saturation considering the two-phase problem, the saturation mentioned later is the water saturation unless otherwise specified.

[0022] Optionally, in step S4, according to the training data set generated in step S3, a corresponding R-U-Net neural network is designed and built as a surrogate model. The surrogate model uses a convolutional neural network (CNN) module to extract and process the spatial feature information of the reservoir state; captures the temporal dependence relationship of injection-production parameters and state parameters through a long short-term memory network (LSTM) module; combines the CNN and LSTM modules to achieve cyclic and continuous prediction of the reservoir state under different injection-production parameter conditions; the surrogate model is constructed by using a strategy of coupled solution of the pressure field and the water saturation field, and reduces the computational complexity by sharing features, laying a model foundation for subsequent coupled numerical simulation training.

[0023] In the input layer, the key permeability field and injection-production parameter field in the prediction of the time-varying well-controlled reservoir state are selected as the main input features. In the processing of injection-production parameters, a sparse matrix form is used. A positive injection parameter is assigned to the position of the production well, and a negative injection-production parameter is assigned to the position of the production well.

[0024] The encoding module uses convolutional operations to perform downsampling on the input features, compresses the spatial dimension through convolutional kernels with strides, thereby extracting higher-level feature representations and reducing the computational amount. Let C represent the number of input feature channels, H and D represent the height and width respectively, then there is an input feature , and the convolutional kernel is , where and represent the spatial dimensions of the convolutional kernel, represents the number of output channels, and the output feature map at this time is , and are the height and width after convolution. The process of the convolutional operation is:

[0025] ;

[0026] In the formula, i, j, and are the indexes of the output feature map in the height, width, and channels respectively, represents the value of the output feature at the channel, i, j position, is the bias term, s represents the stride, represents the (m, n)th element of the convolutional kernel in the c input channel and output channel. At this time, the spatial dimensions of the output feature O are respectively expressed as and ; when s = 2 is set, the spatial dimension of the output feature will be reduced by nearly half, achieving downsampling.

[0027] In the temporal module, the encoded hidden feature F is along the time dimension Input one by one to the input end of each time step of the time series module. This design can ensure that the proxy model captures the dynamic changes of the input sequence over time, thus better understanding and processing time series data. In the time series module, a long short-term memory network (LSTM) is adopted. Through its unique gating mechanism, LSTM can alleviate the problem of gradient disappearance in traditional RNNs. For the input sequence , where represents the input feature at time step n, the hidden state is , and the cell state is . The calculation process of the long short-term memory network LSTM is as follows:

[0028] ;

[0029] In the formula, represents the overall calculation function of the LSTM cell; represents all the gates and candidate states, represents the weight parameter of the forget gate, represents the weight parameter of the input gate, represents the weight parameter of the candidate cell, represents the weight parameter of the output gate, represents the biases corresponding to the forget gate, input gate, candidate cell, and output gate. The output and the cell state at the current time step are jointly determined by the input and the state and at the previous time step. The calculation process is as follows: , , where represents the cell state update function, represents the hidden state generation function.

[0030] In the decoding module, convolution operations are used to separately decode the features of each time step output by the time series module. The convolution process is similar to that of the encoding module. The difference is that s = 1 is set during convolution, and the upsampling process is implemented using interpolation.

[0031] In the prediction result of the output layer, a coupled solution method of the pressure field and water saturation field is adopted to calculate, and the distributions of pressure and water saturation are output simultaneously. While reducing computational resources, the prediction accuracy of the proxy model is increased.

[0032] Optionally, in step S5, the steps of calculating the hybrid loss by coupling the numerical simulator include:

[0033] S51. Agent model initialization and data input, specifically: load the R-U-Net neural network built in step S4 as the agent model, input the training data set in step S3 into the agent model, and set the hyperparameters required for training, including the learning rate, batch size, and number of iterations.

[0034] S52. Calculate the physical loss, specifically: first, predict the state parameters through the agent model and input them into the numerical simulator; under the current state parameters, use the numerical simulator to solve the control equations of the discrete grid and calculate the corresponding residuals, thereby obtaining the physical loss; at the same time, retain the Jacobian matrix in the numerical simulator and use it for the gradient calculation of the physical loss to improve the stability and computational efficiency of the optimization;

[0035] The calculation formula for the physical loss is:

[0036] ;

[0037] In the formula, the gradient of the physical loss with respect to the prediction result is:

[0038] ;

[0039] In the formula, represents the total number of training samples, represents the simulation time step of a single sample, represents the total number of grids included in the reservoir model, represents the reservoir state predicted by the agent model, including the field map distribution of pressure p and water saturation of, and respectively represent the reservoir states predicted by the agent model at the nth and n+1th time steps, represents the gradient information of the physical loss with respect to the prediction result, represents the transpose of the gradient of the loss with respect to the residual in the chain differentiation process, represents the derivative of the residual with respect to the state parameters in the chain differentiation process, that is, the Jacobian matrix J of the control equation.

[0040] S53. Calculate the label loss according to the label data. At this time, the loss gradient information is automatically generated during the calculation of the label loss. The calculation formula is:

[0041] ;

[0042] The gradient of the label loss with respect to the prediction result is:

[0043] ;

[0044] In the formula, Represents the gradient information of the label loss with respect to the prediction result, which is used to update the network parameters during the training process of the surrogate model. Represents the process of automatically differentiating to generate gradients during the training process of the surrogate model.

[0045] S54. Calculate the hybrid loss, and the calculation formula is as follows:

[0046] ;

[0047] In the formula, L represents the hybrid loss.

[0048] Optionally, in step S6, the surrogate model completes the training process under the constraint of the hybrid loss. The permeability field and the injection-production parameter field are used as input features and input into the surrogate model. The reservoir state under different injection-production parameters is predicted by the surrogate model , calculate the label loss according to the label data and automatically generate the gradient information required for backpropagation . At the same time, input the predicted state parameters into the numerical simulator, calculate the mass residual g of the control equation, and represent the gradient of the physical loss with the Jacobian matrix of the control equation . Finally, use the hybrid loss for backpropagation, that is, the formula for updating the network parameters of the surrogate model is:

[0049] ;

[0050] ;

[0051] In the formula, represents the adjustable parameters of the surrogate model, represents the parameter update process, represents the learning rate, represents the gradient of the prediction result with respect to the network parameters, which is calculated by automatic differentiation;

[0052] Thus, the process of constructing a reservoir surrogate model with physical constraints by coupling a numerical simulator is completed.

[0053] The beneficial effects of the present invention are as follows. Aiming at the problems of high computational complexity, difficulty in embedding physical equations, and insufficient adaptability to complex reservoir models in the existing methods for constructing physical-constrained reservoir surrogate models, the present invention aims to provide a new method for constructing a physical-constrained reservoir surrogate model by coupling a numerical simulator. This method incorporates the mass residuals and Jacobian matrices generated during the discrete difference process of the numerical simulator into the surrogate model training to calculate the physical loss as a constraint condition, thereby optimizing the training process of the surrogate model. The present invention combines the accurate modeling ability of the numerical simulator with the efficient prediction advantage of the surrogate model, effectively overcoming the problem of long computational time in the linear solution process of traditional simulators, and at the same time making up for the limitation that it is difficult for the surrogate model to embed complex physical equations, providing an innovative solution for the intelligent numerical simulation of complex reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 It is a flowchart of a new method for constructing a physical-constrained reservoir surrogate model by coupling a numerical simulator according to the present invention;

[0055] Figure 2 It is a distribution diagram of permeability and well positions shown in an embodiment of the present invention;

[0056] Figure 3 It is a schematic diagram of a convolutional temporal neural network for coupled prediction of pressure and water saturation according to the present invention;

[0057] Figure 4A It is a trend diagram of the prediction error of water saturation varying with the number of samples without physical constraints shown in an embodiment of the present invention;

[0058] Figure 4B It is a trend diagram of the prediction error of pressure varying with the number of samples without physical constraints shown in an embodiment of the present invention;

[0059] Figure 5A It is a trend diagram of the prediction error of water saturation varying with the number of samples with physical constraints shown in an embodiment of the present invention;

[0060] Figure 5B It is a trend diagram of the prediction error of pressure varying with the number of samples with physical constraints shown in an embodiment of the present invention;

[0061] Figure 6A It is the prediction of water saturation at each time step by a surrogate model trained with physical constraints added in a coupled numerical simulator based on 50 sets of training data shown in an embodiment of the present invention;

[0062] Figure 6B It is the prediction of pressure at each time step by a surrogate model trained with physical constraints added in a coupled numerical simulator based on 50 sets of training data shown in an embodiment of the present invention;

[0063] Figure 7 The error distribution of the surrogate model trained based on 50 groups of training data with and without physical constraints at the last moment for water saturation and pressure shown in an embodiment of the present invention;

[0064] Figure 8A The comparison chart of the error distribution of the surrogate model for water saturation on 50 groups of random test samples under the conditions of adding physical constraints and not adding physical constraints shown in an embodiment of the present invention;

[0065] Figure 8B The comparison chart of the error distribution of the surrogate model for pressure on 50 groups of random test samples under the conditions of adding physical constraints and not adding physical constraints shown in an embodiment of the present invention. Specific embodiments

[0066] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0067] A new method for constructing a physical constraint reservoir surrogate model by coupling a numerical simulator, as Figure 1 shown, includes the following steps:

[0068] S1. Establish a reservoir model: Use a numerical simulator to construct a three-dimensional oil-water two-phase seepage model. The grid scale of the reservoir model is 40×40×5, the size of each grid is 25 m×25 m×10 m, it includes 1 injection well and 3 production wells, all wells in the reservoir are controlled by bottom-hole pressure for production, the initial oil saturation and water saturation of the reservoir are 0.8 and 0.2 respectively, the porosity is set to 0.15, the initial reservoir pressure is 20 Mpa, the distribution of permeability and well positions is as Figure 2 shown. The control equation is discretized and solved implicitly in space and time by the finite volume method. The control equation is:

[0069] ;

[0070] where x, y, and z represent three-dimensional spatial coordinates, and the subscript l represents the fluid. represents the phase density, represents the phase viscosity, k represents the absolute permeability, represents the relative permeability, p represents the pressure, represents the phase saturation, represents the porosity, represents the flow rate of the fluid in phase l of the injection well or production well, which is calculated using the Peaceman equation: , where represents the thickness of the grid, represents the effective radius of the well, represents the absolute permeability at the well location, represents the actual radius of the well, represents the bottom-hole pressure;

[0071] S2. Randomly perturb the injection and production parameters to generate a dataset of injection and production parameters. The reservoir model is simulated for a total of 1080 days, divided into 6 time steps, with each time step controlling 180 days. During the perturbation, the bottom-hole pressure of the injection well is 25 - 30 Mpa, and the pressure of the production well is 13 - 18 Mpa.

[0072] S3. Input the dataset of injection and production parameters generated by random perturbation in step S2 into the reservoir model, and obtain the reservoir state parameters (distributions of the pressure field and water saturation field) under different injection and production parameters through simulation. A set of injection and production parameters corresponds to a set of state parameters, which together form the training dataset. The input features are selected as the permeability distribution and the injection and production parameters set using a sparse matrix. During the setting of the injection and production parameters, in the form of a sparse matrix, a positive well control value is set for the grid containing the injection well, and a negative well control value is set for the grid of the production well. The prediction results are the distribution states of the pressure and water saturation at different injection and production parameters and different time steps.

[0073] S4. According to the training data format generated in step S3, build an R - U - Net neural network as a surrogate model, as Figure 3As shown in the figure, it includes an input layer, an encoding module, a temporal module, a decoding module, and an output layer. The encoding module extracts spatial features through downsampling by convolution operations. The temporal module uses a long short-term memory network to capture temporal dependencies. The decoding module restores the spatial dimension through upsampling. The input layer receives the permeability field and injection-production parameter field, and the output layer couples and predicts the pressure field and water saturation field. The R-U-Net neural network uses the convolutional neural network (CNN) module to extract and process the spatial feature information of the reservoir state. The long short-term memory network (LSTM) module is used to capture the temporal dependencies of injection-production parameters and state parameters. By combining the CNN and LSTM modules, cyclic and continuous prediction of the reservoir state under different injection-production parameter conditions is achieved. The surrogate model is constructed using the strategy of coupling and solving the pressure field and water saturation field, and the computational complexity is reduced by sharing features, laying a model foundation for subsequent coupled numerical simulation training.

[0074] In the input layer, the key permeability field and injection-production parameter field in the prediction of the time-varying well-controlled reservoir state are selected as input features. In the processing of injection-production parameters, the sparse matrix form is used. Positive injection parameters are assigned to the positions of production wells, and negative injection-production parameters are assigned to the positions of production wells.

[0075] The encoding module uses convolution operations to perform downsampling on the input features, compressing the spatial dimension through convolutional kernels with strides to extract higher-level feature representations and reduce the computational amount. Let C represent the number of input feature channels, H and D represent the height and width respectively, then the input feature is , and the convolutional kernel is , where and represent the spatial dimensions of the convolutional kernel, represents the number of output channels. At this time, the output feature map is , and are the height and width after convolution. The process of the convolution operation is:

[0076] ;

[0077] In the formula, i, j, and are the indices of the output feature map in height, width, and channel respectively. represents the value of the output feature at the position of channel , i, j. is the bias term, s represents the stride, represents the (m, n)th element of the convolutional kernel in the c input channel and output channel. At this time, the spatial dimensions of the output feature O are respectively represented as and ; when s = 2 is set, the spatial dimensions of the output feature will be reduced by nearly half, achieving downsampling processing.

[0078] In the timing module, the encoded hidden feature F is input one by one into the input end of each time step of the timing module along the time dimension This design can ensure that the surrogate model captures the dynamic changes of the input sequence in time, so as to better understand and process the timing data. In the timing module, a long short-term memory network (LSTM) is adopted. Through its unique gating mechanism, the LSTM can alleviate the problem of gradient disappearance in traditional RNNs. For the input sequence , where represents the input feature at time step n, the hidden state is , and the cell state is , the calculation process of the long short-term memory network (LSTM) is as follows:

[0079] ;

[0080] In the formula, represents the overall calculation function of the LSTM cell; represents all the gates and candidate states, represents the weight parameter of the forget gate, represents the weight parameter of the input gate, represents the weight parameter of the candidate cell, represents the weight parameter of the output gate, represents the bias corresponding to each part, and the output and the cell state at the current time step are jointly determined by the input and the state and at the previous time step, and the calculation process is as follows: , , where represents the cell state update function, represents the hidden state generation function.

[0081] In the decoding module, convolution operations are used to separately decode the features of each time step output by the timing module. The convolution process is similar to that of the encoding module. The difference is that s = 1 is set during convolution, and the upsampling process is implemented using interpolation.

[0082] In the prediction result of the output layer, a coupled solution method of the pressure field and the water saturation field is adopted to calculate, and the distributions of pressure and water saturation are output simultaneously, reducing the computing resources while increasing the prediction accuracy of the surrogate model.

[0083] S5. Calculate the hybrid loss by coupling the numerical simulator; including:

[0084] S51. Agent model initialization and data input, specifically: Load the R-U-Net neural network built in step S4 as the agent model, input the training dataset in step S3 into the network, and set the hyperparameters required for training the agent model, including the learning rate, batch size, and number of iterations.

[0085] S52. Calculate the physical loss, specifically: First, predict the state parameters through the agent model and input them into the numerical simulator; under the current state parameters, use the numerical simulator to solve the control equations of the discrete grid and calculate the corresponding residuals, thereby obtaining the physical loss; at the same time, retain the Jacobian matrix in the numerical simulator and use it for the gradient calculation of the physical loss to improve the stability and computational efficiency of the optimization.

[0086] Taking three-dimensional oil-water two-phase seepage as an example, the control equations were given in step S1. After discrete difference, it was simplified to: , where g represents the mass residual, and A, F, and Q represent the cumulative phase, flowing phase, and source-sink phase respectively. and represent the reservoir states at times n+1 and n respectively, including the distributions of saturation and pressure; represents the injection-production parameters at time n+1. The numerical simulator efficiently calculates the mass residual g based on the static parameters of the reservoir and the injection-production parameters, and then calculates to constrain the training of the agent model. Thus, during the training of the agent model, the physical loss is:

[0087] ;

[0088] In the formula, the gradient of the physical loss with respect to the prediction result is:

[0089] ;

[0090] In the formula, represents the total number of training samples, represents the simulation time step of a single sample, represents the total number of grids included in the reservoir model, represents the reservoir state predicted by the agent model, including the field map distributions of pressure p and water saturation , and represent the reservoir states predicted by the agent model at the nth and (n+1)th time steps respectively, represents the gradient information of the physical loss with respect to the prediction result, represents the transpose of the gradient of the loss with respect to the residual during the chain differentiation process, represents the derivative of the residual with respect to the state parameters during the chain differentiation process, that is, the Jacobian matrix J of the control equation.

[0091] The numerical simulator calculates the residual g of the control equation and the Jacobian matrix J, and then can efficiently calculate the physical loss, which is used to ensure that the reservoir proxy model constructed by coupling the numerical simulator with physical constraints has good applicability.

[0092] S53. Calculate the label loss. At this time, the loss gradient information is automatically generated during the calculation of the label loss, and the calculation formula is:

[0093] ;

[0094] The gradient of the label loss with respect to the prediction result is:

[0095] ;

[0096] In the formula, represents the gradient information of the label loss with respect to the prediction result, which is used to update the network parameters during the training of the proxy model. represents the process of automatically differentiating to generate the gradient during the training of the proxy model.

[0097] S54. Calculate the hybrid loss. The calculation formula is as follows:

[0098] ;

[0099] In the formula, L represents the hybrid loss.

[0100] Use the numerical simulator to build a reservoir model, which provides a model basis for the training of the proxy model. Generate label data by randomly perturbing the injection and production parameters, input the label data into the proxy model, and the proxy model gives the predicted reservoir state parameters. Calculate the label loss according to the label data , and at this time, the gradient of the label loss will be automatically generated. In addition, input the predicted reservoir state parameters into the numerical simulator, and calculate the mass residual g of the control equation and the Jacobian matrix J under the current state. Thus, the calculation of the hybrid loss is completed.

[0101] S6. Use the gradient of the hybrid loss for backpropagation to update the network parameters of the proxy model, and complete the training of the proxy model. The proxy model completes the training process under the constraint of the hybrid loss. The formula for using the hybrid loss to constrain the update of the network parameters of the proxy model is:

[0102] ;

[0103] ;

[0104] In the formula, represents the adjustable parameters of the proxy model. represents the parameter update process. denotes the learning rate, denotes the gradient of the prediction result with respect to the network parameters, which is obtained by automatic differentiation;

[0105] Thus, the process of constructing a physically constrained reservoir surrogate model by the coupled numerical simulator is completed.

[0106] Figure 4A 、 Figure 4B 、 Figure 5A 、 Figure 5B respectively show the variation of the prediction errors of the unconstrained surrogate model and the constrained surrogate model with respect to water saturation and pressure with the amount of training data. It can be clearly seen from the figure that the prediction error of the surrogate model gradually decreases with the increase in the number of training samples. Among them, after introducing physical constraints, only 50 groups of training samples are required to complete the training of the surrogate model. At this time, the prediction accuracy of the surrogate model is close to the accuracy level when using 200 groups of training data without constraints. It can be obtained that the method of constructing a physically constrained reservoir surrogate model based on the coupled numerical simulator can reduce the required number of training samples by nearly 75%.

[0107] Figure 6A 、 Figure 6B show the prediction results of the water saturation and pressure of a set of random test samples by the surrogate model trained with 50 groups of samples in this embodiment. It can be seen from the figure that this method can efficiently construct a high-precision surrogate model, where the prediction error of the surrogate model for water saturation is mostly less than 0.01, and the prediction error for pressure is mostly less than 0.25 MPa. This shows that the method of constructing a physically constrained reservoir surrogate model based on the coupled numerical simulator has good applicability in three-dimensional oil-water two-phase seepage problems.

[0108] Figure 7 show the distribution of the pressure and water saturation errors at the last time step when using the same 50 groups of training samples. It can be seen from the figure that the proposed method of introducing physical constraints to construct the surrogate model can effectively improve the prediction accuracy of the surrogate model. Figure 8A 、 Figure 8B Further compare the prediction error distributions of the surrogate models constructed by the two methods with and without physical constraints on 50 groups of random test samples. The results show that by using the proposed physical constraint method, on the premise of maintaining the same training cost, the average prediction error of pressure is reduced from 1.24% to 0.95%, a reduction of about 25%; while the average prediction error of water saturation is reduced from 2.98% to 0.87%, a reduction of about 70%.

[0109] Certainly, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the scope of the essence of the present invention shall also fall within the protection scope of the present invention.

Claims

1. A new method for constructing a physically constrained reservoir proxy model by coupling a numerical simulator, characterized in that: The steps include: S1. Establishing a reservoir model: Using a numerical simulator to construct a three-dimensional oil-water two-phase flow model, the control equations are implicitly discretized in space and time and solved using the finite volume method; S2, randomly perturb the injection and production parameters to generate an injection and production parameter data set; S3, inputting the injection-production parameter data set generated in step S2 into the reservoir model for simulation to obtain corresponding reservoir state parameters, where the injection-production parameters and the state parameters together constitute label data; S4. Build RU-Net neural network as the proxy model, including input layer, encoding module, timing module, decoding module and output layer. The encoding module extracts spatial features by downsampling through convolution operation. The timing module uses long short-term memory network to capture temporal dependencies. The decoding module restores the spatial dimension by upsampling. The input layer receives the permeability field and the injection and production parameter field. The output layer couples and predicts the pressure field and saturation field. S5, coupled numerical simulator to calculate mixing losses; S6. Use mixed loss back propagation to update the network parameters of the proxy model to complete the training of the proxy model.

2. A novel method for constructing a physically constrained reservoir proxy model by coupling a numerical simulator as claimed in claim 1, characterized in that: In step S1, the control equation is: ； In the formula, x, y and z represent the three-dimensional space coordinates, the subscript l represents the fluid, represents the phase density, represents phase viscosity, k represents absolute permeability, represents relative permeability, p represents pressure, represents the phase saturation, represents the porosity, represents the flow rate of the l-phase fluid in the injection well or production well, calculated using the Peaceman equation: ,in, represents the thickness of the grid, represents the effective radius of the well, represents the absolute permeability at the well location, represents the actual radius of the well, Indicates the bottom hole pressure; The control equation is simplified to: ； In the formula, g represents the mass residual, A, F and Q represent the accumulation phase, mobile phase and source and sink phase respectively; and Represent the reservoir status at time n+1 and time n, including the distribution of saturation and pressure; Indicates the injection and production parameters at time n+1.

3. A novel method for constructing a physically constrained reservoir proxy model by coupling a numerical simulator as claimed in claim 1, characterized in that: In step S3, the injection and production parameters in the label data are represented in the form of a sparse matrix, the well control value of the grid corresponding to the injection well is positive, and the well control value of the grid corresponding to the production well is negative.

4. A novel method for constructing a physically constrained reservoir proxy model by coupling a numerical simulator as claimed in claim 1, characterized in that: In step S4, the encoding module uses convolution operation to downsample the input features, and compresses the spatial dimension through convolution kernel with stride. C represents the number of input feature channels, H and D represent the height and width respectively, then the input feature , the convolution kernel is ,in, and represents the spatial size of the convolution kernel, Indicates the number of output channels, the output feature map at this time , and is the height and width after convolution, and the convolution operation process is: ； In the formula, i, j and are the indices of the output feature map in height, width and channel, respectively. Indicates that the output feature is The value at the channel, i, j position, is the bias term, s represents the stride, Indicates that the convolution kernel is in c input channels and The (m,n)th element on the output channel, the spatial size of the output feature O is expressed as and ; In the timing module, for the input sequence ,in, represents the input features of time step n, and the hidden state is , the cell state is , the calculation process of the long short-term memory network LSTM is: ； In the formula, Represents the overall calculation function of the LSTM unit; represents all gates and candidate states, represents the weight parameter of the forget gate, represents the weight parameter of the input gate, represents the weight parameter of the candidate cells, represents the weight parameter of the output gate, Represents the bias corresponding to the forget gate, input gate, candidate cell and output gate, and the output of the current time step and cell status By input and the state at the previous time step and Jointly decided, the calculation process is: , ,in, represents the cell state update function, represents the hidden state generating function.

5. A novel method for constructing a physically constrained reservoir proxy model by coupling a numerical simulator as claimed in claim 1, characterized in that: In step S5, the step of calculating the mixing loss by coupling the numerical simulator includes: S51, proxy model initialization and data input, specifically: load the RU-Net neural network built in step S4 as the proxy model, input the label data in step S3 into the proxy model, and set the hyperparameters required for proxy model training, including learning rate, batch size, and number of iterations; S52, calculating the physical loss, specifically: first predicting the state parameters through the proxy model and inputting them into the numerical simulator; under the current state parameters, solving the control equation of the discrete grid using the numerical simulator, and calculating the corresponding residual, so as to obtain the physical loss; at the same time, retaining the Jacobian matrix of the control equation in the numerical simulator, and using it for the gradient calculation of the physical loss; The calculation formula for physical loss is: ； In the formula, the gradient of physical loss to the prediction result is: ； In the formula, represents the total number of training samples, represents the simulation time step of a single sample, represents the total number of grids contained in the reservoir model, Represents the reservoir state predicted by the proxy model, including pressure p and water saturation The field distribution of and Represent the reservoir state predicted by the proxy model at the nth and n+1th time steps, respectively. Represents the gradient information of physical loss for the prediction result, represents the gradient transpose of the loss with respect to the residual in the chain derivation process, It represents the derivative of the residual with respect to the state parameter in the chain derivation process, that is, the Jacobian matrix J of the control equation; S53. Calculate label loss based on label data. The calculation formula is: ； The gradient of the label loss with respect to the prediction result is: ； In the formula, Represents the gradient information of the label loss for the prediction result, which is used to update the network parameters during the training process of the proxy model. Represents the process of automatically differentiating and generating gradients during the training of the surrogate model; S54. Calculate the mixed loss. The calculation formula is as follows: ； Where L represents the mixing loss.

6. A novel method for constructing a physically constrained reservoir proxy model by coupling a numerical simulator as claimed in claim 1, characterized in that: In step S6, the formula for updating the network parameters of the proxy model with the mixed loss constraint is: ； ； In the formula, represents the adjustable parameters of the surrogate model, represents the parameter update process, represents the learning rate, Represents the gradient of the prediction result with respect to the network parameters, which is calculated by automatic differentiation; Therefore, the physical loss and label loss jointly constrain the training process of the proxy model.

Images