An intelligent dynamic analysis and prediction method for oil reservoirs based on physical guiding neural operators

Abstract

The application discloses a kind of based on physical guiding neural operator oil reservoir intelligent dynamic analysis prediction method, belong to oil and gas field development numerical simulation and artificial intelligence field, first based on the data set containing the space-time dynamics of multiphase flow of heterogeneous reservoir model is constructed, and the dimension is reduced by layer slicing method.Processing.Using neural network to extract the space-time evolution characteristics of reservoir fluid, and a double-track multi-head prediction structure is designed, including pressure prediction track and saturation prediction track, wherein the space-time evolution kernel representing non-steady seepage diffusion law is embedded in the pressure prediction track.The design directly uses the physical law as the guiding benchmark for network forward propagation, avoiding the computational disaster caused by solving complex multiphase partial differential equations, and finally optimizing the model by combining multiple constraints such as spatial gradient, boundary and initial conditions.The application realizes the guidance of physical prior to data-driven neural operator, greatly improves the analysis and prediction accuracy, generalization ability and computational efficiency of the model.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation and artificial intelligence in oil and gas field development, specifically involving a method for intelligent dynamic analysis and prediction of oil reservoirs based on physical guided neural operators. Background Technology

[0002] Accurate prediction of reservoir dynamic parameters is the core foundation for oil and gas field development planning and production optimization. Traditional reservoir numerical simulation methods achieve prediction by solving seepage control equations, but numerical simulation requires discrete grids and iterative solutions, resulting in high computational complexity and long processing time, making it difficult to meet the needs of rapid prediction for complex heterogeneous reservoirs. In recent years, neural operators based on deep learning have been introduced into the reservoir field. These methods learn high-dimensional spatial mapping relationships through data-driven learning, significantly improving prediction speed. However, the purely data-driven approach lacks physical constraints, and the prediction results are prone to violating the basic laws of seepage mechanics, leading to physical distortions. It also has weak generalization ability for complex heterogeneous reservoirs, and pressure prediction does not incorporate the physical nature of seepage field evolution over time, making it prone to global offset and local oscillation problems. On the other hand, although the spatiotemporal evolution kernel based on diffusion equations can rigorously characterize the spatiotemporal transmission law of pressure waves, its direct application faces the difficulty of finding explicit source terms, and the analytical kernel requires a completely homogeneous medium, making it difficult to adapt to real complex heterogeneous reservoir scenarios. Therefore, how to design a novel network architecture that uses the analytical physical transmission laws as the underlying benchmark to guide neural operators, and how to build a complete physical constraint system to ensure the physical consistency of prediction results while overcoming strong heterogeneity obstacles, has become a key technical problem that urgently needs to be solved in the field of reservoir dynamic prediction. Summary of the Invention

[0003] To address the issues of long computation time in traditional reservoir numerical simulation and the tendency of pure data-driven deep learning prediction models to violate seepage physics and have weak generalization ability, this invention proposes an intelligent dynamic analysis and prediction method for reservoirs based on physical guided neural operators. The method is rationally designed, overcomes the shortcomings of existing technologies, and has good results.

[0004] A method for intelligent dynamic analysis and prediction of oil reservoirs based on physical-guided neural operators includes the following steps: Step 1: Generate a large number of heterogeneous reservoir numerical models based on sequential Gaussian, and obtain time-series dynamic data of pressure field and three-phase saturation field through production simulation to form a reservoir dynamic prediction dataset. Step 2: Divide the 3D reservoir grid in the dataset into multiple 2D grid slices along the depth Z-axis, normalize and stitch the features of each slice, and construct the multiphase flow feature tensor corresponding to each slice; Step 3: For the multiphase flow feature tensor of each two-dimensional mesh slice, construct a Fourier neural operator model with embedded frequency domain analytical prior. The model includes a multi-layer spatiotemporal feature extraction module and a dual-track multi-head prediction structure. First, extract the spatiotemporal evolution features of the slice through the multi-layer spatiotemporal feature extraction module to obtain the latent feature tensor. Step 4: Input the latent feature tensor into the dual-track multi-head prediction structure, and output the pressure field and multiphase saturation field with seepage law in parallel by decoupling the pressure prediction track and the saturation prediction track. Step 5: Construct a joint loss function for data constraints, gradient constraints, initial condition constraints, and boundary condition constraints, and train and optimize the model. Step 6: Utilize the trained Fourier neural operator models to perform deduction, and execute the three-dimensional dynamic field reconstruction and output in all spatiotemporal dimensions.

[0005] Furthermore, step 1 includes the following sub-steps: Step 1.1: Generate a batch random heterogeneous permeability field and a heterogeneous porosity field using the sequential Gaussian simulation method; Step 1.2: Establish a three-dimensional three-phase reservoir numerical benchmark model, configure the grid size and reservoir physical property parameters, set the injection and production well pattern and well control conditions, and import the generated heterogeneous porosity field into the benchmark model; Step 1.3: Import each generated heterogeneous permeability field into the reservoir numerical benchmark model one by one, replacing the initial homogeneous permeability field to form multiple heterogeneous reservoir numerical models; perform production simulation on each model, extract the time-series dynamic data of pressure field and oil, gas and water three-phase saturation field, and save them together with the corresponding permeability field to form a reservoir dynamic prediction dataset.

[0006] Furthermore, step 2 includes the following sub-steps: Step 2.1: Obtain the static geological features and dynamic temporal features of the 3D reservoir grid from the dataset. The static geological features include absolute permeability features in the X, Y, and Z directions, while the dynamic temporal features include historical features of pressure, gas saturation, oil saturation, and water saturation. The grid with a size of [size missing] is then [details missing]. The three-dimensional reservoir grid is divided into layers along the depth Z-axis. Each size is Two-dimensional mesh slices; Step 2.2: For each 2D mesh slice, extract the previous... The dynamic time series characteristics of each historical time step are obtained, and Gaussian normalization is used to perform dimensionless processing on static and dynamic data respectively. Step 2.3: Concatenate the normalized data along the channel dimension, including the absolute penetration features of the three channels, and each channel contains... The historical characteristics of pressure, gas saturation, oil saturation, and water saturation in each channel are used to generate a high-dimensional input feature tensor of the physical field as the feature tensor of the multiphase flow.

[0007] Furthermore, step 3 includes the following sub-steps: Step 3.1: Concatenate the multiphase flow feature tensor obtained in Step 2 with a 3D spatial and temporal relative coordinate grid tensor in the channel dimension; then, through a fully connected mapping layer, map the number of channels of the concatenated tensor to the preset hidden layer feature dimension W; apply zero-padding to the tensor in the spatiotemporal dimension to obtain the initial hidden feature tensor. ; Step 3.2: Initialize the hidden feature tensor The input consists of a stacked multi-layer spatiotemporal feature extraction module. Within each module, a two-branch parallel computation is employed, with the first branch being a three-dimensional frequency domain spectral convolution branch. The specific operations are as follows: Step 3.2.1: Process the input feature tensor of the current layer. By applying a three-dimensional fast real Fourier transform to the spatial domain (X, Y) and the time domain (T), the frequency domain feature tensor is obtained. : ; in, Represents the three-dimensional fast real Fourier transform; They are respectively Frequency wavenumber of dimension for Discrete angular frequencies corresponding to the time dimension; Step 3.2.2: Perform mode truncation in the frequency domain, retaining only the low-frequency mode components in three dimensions; Step 3.2.3: Perform cross-channel tensor multiplication on the channel dimension of the truncated frequency domain feature tensor and the complex weight tensor based on the Einstein summation convention: ; in, For the frequency domain complex intermediate feature tensor, For the truncated feature tensor, This represents a cross-channel tensor multiplication operation based on Einstein's summation convention. It is a learnable complex weight matrix; Step 3.2.4, for The truncated high-frequency components are padded, and then an inverse three-dimensional real Fourier transform is applied to map them back to the physical spatiotemporal domain, yielding the spectral convolution feature tensor. ; Step 3.3: The second branch is a spatial domain local convolution branch, using a 3D convolution kernel with a kernel size of 1×1×1, a stride of 1, and padding of 0, to process the input feature tensor. By directly performing a linear mapping in the spatial domain, we obtain a tensor that captures local residual features. ; Step 3.4: Convert the spectral convolution feature tensor obtained in Step 3.2.4 into... The local residual tensor obtained in step 3.3 The data is added element-wise, then the data distribution is stabilized through a three-dimensional batch normalization layer, and nonlinear processing is performed using the GELU activation function to obtain the output feature tensor of the current module. ; Step 3.5: Output the feature tensor As input to the next layer feature extraction module, steps 3.2 to 3.4 are repeated L times to prune the final output tensor, removing the boundary zero-padding components added in step 3.1, thus obtaining the final high-dimensional spatiotemporal evolution latent feature tensor. .

[0008] Furthermore, step 4 includes the following sub-steps: Step 4.1: The dual-track multi-head prediction structure will use the latent feature tensor The data are input into the pressure prediction orbital and the saturation prediction orbital, respectively; the pressure prediction orbital adopts a dual-branch parallel structure with a data residual branch and a physical prior branch. Step 4.2: In the data residual branch, the latent feature tensor The input consists of a multilayer perceptron composed of two fully connected layers and the GELU activation function, and the output is a purely data-driven local pressure residual field tensor. ; Step 4.3: In the physical prior branch, perform analytical injection of the frequency domain spatiotemporal evolution kernel; the specific process is as follows: Step 4.3.1: Use a single-layer fully connected network to perform channel compression on the latent feature tensor, and output a single-channel virtual percolation source term. This is used to equivalently characterize formation heterogeneous disturbances and multiphase fluid interference; by introducing the operator equivalent substitution principle, an equivalent diffusion equation is constructed: ; in, As the reference pressure field; The overall pressure conductivity coefficient of an ideal medium; For the Laplace operator; Step 4.3.2: First, apply zero-causal padding of length equal to the prediction time step to the end of the time axis of the virtual source term tensor. Then, perform a three-dimensional fast real Fourier transform to obtain the frequency domain source term tensor. : ; Step 4.3.3: Construct a discrete two-dimensional spatial wavenumber grid and a one-dimensional time angular frequency grid, and define a spatiotemporal evolution kernel carrying the dual physical properties of time delay and spatial radial diffusion: ; in, As the core of spatiotemporal evolution; The imaginary unit; A small constant that is set; Step 4.3.4: Perform a complex element-wise product operation between the frequency domain source term tensor and the frequency domain spatiotemporal evolution kernel to solve for the frequency domain reference pressure field: ; in, The frequency domain reference pressure field; This represents the element-wise multiplication operation of a tensor with complex numbers. Step 4.3.5: For the frequency domain reference pressure field... An inverse three-dimensional real Fourier transform is applied to map back to the physical spacetime domain, and then the filling portion introduced in step 4.3.2 is trimmed along the time axis to obtain the reference pressure field. ; Step 4.4: Output the pressure residual field from the data residual branch. The reference pressure field output by the physical prior branch To achieve fusion, learnable dynamic equilibrium weight parameters are introduced. Output the final predicted pressure field: ; in, The pressure residual field is output by the data residual branch; The reference pressure field output by the physical prior branch; This represents the final pressure field after the merger. These are learnable, dynamically balanced weights. Step 4.5: In the saturation prediction trajectory, a multi-head parallel prediction structure is adopted, using a two-layer fully connected network to process the latent feature tensor. Mapping and decoding are performed; the first layer expands the number of channels to 128 and activates them using GELU, the second layer reduces the number of channels to 3, and outputs the gas saturation prediction field for future time steps in parallel. Oil saturation prediction field and water saturation prediction field ; Step 4.6: Concatenate the pressure field tensor predicted in Step 4.4 with the three-phase saturation field tensor predicted in Step 4.5 in the channel dimension to obtain the final output multiphase flow spatiotemporal dynamic field prediction tensor with 4 channels.

[0009] Furthermore, step 5 includes the following sub-steps: Step 5.1: For any physical field variable predicted by the network and its corresponding true label variable, calculate the composite relative norm error, which includes the fundamental physical quantity and the spatiotemporal gradient. : ; in, Represents any real physical field; This represents the corresponding prediction field; Represents the L2 norm of a tensor; Indicates in A comprehensive difference gradient operator in dimensionality; These are the gradient regularization weight coefficients; To prevent constants with a denominator of zero; The loss was calculated separately for the pressure field, gas-bearing field, oil-bearing field, and water-saturation field, and then weighted and summed to obtain the global data principal loss. : ; in, These represent the actual pressure, gas content, oil content, and water saturation fields, respectively. For the corresponding prediction field; This is a weighting coefficient specific to gas saturation. Step 5.2: Calculate the physical boundary condition loss. : ; in, The normal partial derivative operator is perpendicular to the boundary. This represents the total number of grid nodes involved in the boundary condition calculation. To predict physical fields At the outer boundary of the reservoir The values ​​at each grid node; Step 5.3: Calculate the initial condition loss. : ; in, The actual physical field at the initial time step; The physical field predicted by the network at the initial time step; Step 5.4: Construct the final total loss function : ; in, and These are the weighting coefficients for the boundary conditions and the initial conditions, respectively.

[0010] Furthermore, step 6 includes the following sub-steps: Step 6.1: After slicing the 3D data in the test set into layers, input them in parallel into the corresponding pre-trained 2D models of each layer; the models automatically perform spatiotemporal feature extraction and physical dual-track fusion, directly inferring and generating the future of each layer. Prediction results of pressure field and three-phase saturation field at each time step; Step 6.2: Perform denormalization on the dimensionless prediction tensor output by the network to restore it to its true physical dimensions; then, in the spatial dimension, The two-dimensional slices are spliced ​​together along the depth Z-axis to reassemble into a complete shape. The three-dimensional multiphase flow dynamic field, The number of samples in the test set.

[0011] The beneficial technical effects of this invention are as follows: This invention offers significant advantages over existing reservoir dynamic prediction technologies in multiple dimensions. By integrating physical priors of the seepage field into the neural operator architecture and employing multi-dimensional loss functions, it ensures that the prediction results conform to the fundamental laws of reservoir seepage mechanics, avoiding the problem of physical distortion easily encountered in purely data-driven models. The pressure branch adopts a dual-track fusion design of the data residual field and the physically guided reference field generated in the frequency domain by the spatiotemporal evolution kernel. This effectively solves the global offset and local oscillation defects that are prone to occur in traditional data-driven pressure prediction models. Combined with gradient-based L2 loss, it further enhances the spatiotemporal continuity of the prediction results, significantly improving the prediction accuracy of reservoir pressure and multiphase saturation dynamic parameters. Compared with traditional reservoir numerical simulation methods, this invention significantly improves prediction efficiency. At the same time, through the deep integration of physical priors and the joint optimization of multi-dimensional loss functions, this invention significantly reduces the model's dependence on the distribution of training data, demonstrating good adaptability to heterogeneous reservoirs with different permeability distributions and different development stages. It significantly improves the model's prediction accuracy and generalization ability, ultimately realizing end-to-end intelligent analysis and prediction from geological models to multiphase flow dynamic evolution. Attached Figure Description

[0012] Figure 1 This is a flowchart of the intelligent dynamic analysis and prediction method for oil reservoirs based on physical guided neural operators, as described in this invention.

[0013] Figure 2 This is a flowchart of the neural network model structure of the present invention.

[0014] Figure 3 This is a permeability field map of a random heterogeneous oil reservoir used in the dataset of this embodiment of the invention; Wherein, (a) is random permeability field 1; (b) is random permeability field 2; (c) is random permeability field 3; and (d) is random permeability field 4.

[0015] Figure 4 This is a schematic diagram of the injection-production well network location in an embodiment of the present invention.

[0016] Figure 5 These are the actual pressure field map, predicted pressure field map, and absolute error map for 1500 days in this embodiment of the invention.

[0017] Figure 6 These are the actual oil saturation field map, predicted oil saturation map, and absolute error map for 1500 days in this embodiment of the invention.

[0018] Figure 7 The images show the actual gas saturation field map, predicted gas saturation map, and absolute error map for 1500 days in this embodiment of the invention.

[0019] Figure 8 The images show the actual water saturation field map, predicted water saturation map, and absolute error map for 1500 days in this embodiment of the invention. Detailed Implementation

[0020] The specific embodiments of the present invention will be further described below with reference to specific examples: A method for intelligent dynamic analysis and prediction of oil reservoirs based on physical-guided neural operators, such as Figure 1 As shown, it includes the following steps: Step 1: Generate a large number of heterogeneous reservoir numerical models based on sequential Gaussian, and obtain time-series dynamic data of pressure field and three-phase saturation field through production simulation to form a reservoir dynamic prediction dataset. Includes the following sub-steps: Step 1.1: Generate a batch of random heterogeneous permeability field and a heterogeneous porosity field using the sequential Gaussian simulation method. Fix the global random seed to ensure that the generation process is repeatable. Complete the configuration of basic parameters such as planar grid size, number of vertical layers, and total number of samples. Step 1.1.1: For the batch random permeability field, set the planar grid size, assign an independent random sub-seed to each sample to ensure the difference in heterogeneous characteristics between samples, independently and randomly configure the core simulation parameters such as the principal range, principal azimuth, and logarithmic space variance of the exponential covariance function for each sample, and perform sequential Gaussian simulation in a loop to generate a logarithmic space Gaussian random field. After overall offset, upper and lower limit truncation, and exponential transformation, the physical range constraint is completed. Then, the permeability field is smoothed by three-dimensional Gaussian filtering and scaled to a reasonable range for reservoir engineering by min-max normalization. Finally, a batch of permeability fields with independent random heterogeneous characteristics is obtained and saved as a .mat format file. Step 1.1.2: For the heterogeneous porosity field, set the planar grid size, configure the simulation parameters of the corresponding exponential covariance function, perform sequential Gaussian simulation to obtain the Gaussian random field, and perform absolute value processing, three-dimensional Gaussian filtering smoothing and maximum-maximum normalization scaling to a reasonable porosity range to obtain a single-sample heterogeneous porosity field and save it as a .mat file.

[0021] Step 1.2: Establish a three-dimensional three-phase reservoir numerical benchmark model using a numerical simulator, configure the grid size and reservoir physical property parameters including the initial homogeneous permeability field, initial homogeneous porosity field, initial pressure, initial saturation, rock compressibility coefficient, density, volume factor, and relative permeability curve, set the injection-production well pattern and well control conditions, disregard capillary force, and import the generated heterogeneous porosity field into the benchmark model; In this embodiment, the total number of grids is 36,864, the reservoir model size is 96×96×4, the grid size is 10×10×5, and 5 injection wells and 4 production wells are set up. The specific well locations are as follows: Figure 4 As shown, the injection wells and production wells were subjected to well control conditions of constant injection rate and constant bottomhole flowing pressure, respectively, and a numerical simulator was used to simulate development for 1500 days; the permeability field map of the stochastic heterogeneous reservoir used is shown in the figure. Figure 3 As shown in (a) to (d); Step 1.3: Import each generated heterogeneous permeability field into the reservoir numerical benchmark model one by one, replacing the initial homogeneous permeability field to form multiple heterogeneous reservoir numerical models; perform production simulation on each model, extract the time-series dynamic data of pressure field and oil, gas and water three-phase saturation field, and save them together with the corresponding permeability field as a .mat format file to form a reservoir dynamic prediction dataset.

[0022] Step 2: To eliminate cross-layer gradient interference caused by extreme vertical heterogeneity, a quasi-three-dimensional decoupling strategy is adopted. The three-dimensional (oil, gas, water) reservoir data in the dataset is divided into multiple two-dimensional grid slices according to the depth Z-axis. Each slice is normalized and its features are stitched together to construct the multiphase flow feature tensor corresponding to each slice. Includes the following sub-steps: Step 2.1: Obtain the static geological features and dynamic temporal features of the 3D reservoir grid from the dataset. The static geological tensor includes absolute permeability features in the X, Y, and Z directions, while the dynamic temporal tensor includes historical features of pressure, gas saturation, oil saturation, and water saturation. The original 3D reservoir grid size is [size missing]. It is divided into layers along the depth Z-axis. A two-dimensional grid slice (size is) ); Step 2.2: For each 2D mesh slice, extract the previous... The dynamic time series characteristics of each historical time step are obtained, and the static and dynamic data are dimensionless by using Unit Gaussian Normalizer. Step 2.3: Concatenate the normalized data along the channel dimension, including the absolute penetration features of the three channels, and each channel contains... The historical characteristics of pressure, gas saturation, oil saturation, and water saturation in each channel are used to generate a high-dimensional input feature tensor of the physical field as the feature tensor of the multiphase flow.

[0023] Step 3: For the multiphase flow feature tensor of each two-dimensional mesh slice, construct a Fourier neural operator model with embedded frequency domain analytical priors, such as... Figure 2 As shown, the model includes a multi-layer spatiotemporal feature extraction module and a dual-track multi-head prediction structure. First, the spatiotemporal evolution features of the slice are extracted through the multi-layer spatiotemporal feature extraction module to obtain the latent feature tensor. The multi-layer spatiotemporal feature extraction module achieves long-distance spatiotemporal dependency modeling and local residual feature fusion through parallel three-dimensional frequency domain spectral convolution and spatial domain local convolution. Includes the following sub-steps: Step 3.1: The multiphase flow feature tensor (including permeability, pressure, and three-phase saturation) obtained in Step 2 is concatenated with a three-dimensional spatial and temporal relative coordinate grid tensor in the channel dimension. This includes a horizontal spatial coordinate tensor, a vertical spatial coordinate tensor, and a future time step coordinate tensor. The spatial coordinate tensor corresponds to the horizontal position of the grid point within the two-dimensional slice, the vertical spatial coordinate tensor corresponds to the vertical position of the grid point within the two-dimensional slice, and the future time coordinate tensor is used to represent the temporal position of the future prediction time step, thus giving the model an explicit spatiotemporal position awareness. This constructs the final high-dimensional input feature tensor of the operator network (the total number of channels includes the sum of static geological features, historical dynamic features, and spatiotemporal coordinate features). Subsequently, a fully connected mapping layer is used to map the number of channels of the concatenated tensor to the preset hidden layer feature dimension W. To avoid the periodic boundary effect generated by the Fourier transform, zero padding is applied to the tensor in the spatiotemporal dimension to obtain the initial hidden feature tensor. ; Step 3.2: Initialize the hidden feature tensor The input consists of a stacked multi-layer spatiotemporal feature extraction module. Within each module, a two-branch parallel computation is employed, with the first branch being a three-dimensional frequency domain spectral convolution branch. The specific operations are as follows: Step 3.2.1: Process the input feature tensor of the current layer. Applying a three-dimensional fast real Fourier transform (3D rFFT) to the spatial domain (X, Y) and time domain (T) maps the physical space-time domain to the frequency domain, obtaining the frequency domain feature tensor. : ; in, Represents the three-dimensional fast real Fourier transform; They are respectively Frequency wavenumber of dimension for Discrete angular frequencies corresponding to the time dimension; Step 3.2.2: In order to filter out high-frequency noise and reduce computational complexity, mode truncation is performed in the frequency domain, retaining only the first M1, M2, and M3 low-frequency mode components in three dimensions. Step 3.2.3: Define the learnable complex weight tensor The truncated frequency domain feature tensor and the complex weight tensor are multiplied across channels using a cross-channel tensor multiplication based on Einstein's summation convention, thereby achieving cross-channel frequency domain feature mixing and filtering. ; in, For the frequency domain complex intermediate feature tensor, For the truncated feature tensor, This represents a cross-channel tensor multiplication operation based on Einstein's summation convention; Step 3.2.4, for The truncated high-frequency components are padded, and then an inverse three-dimensional real Fourier transform (3D irFFT) is applied to map them back to the physical spatiotemporal domain, obtaining the spectral convolution feature tensor. ; Step 3.3: The second branch is a spatial domain local convolution branch, using a 3D convolution kernel (Conv3d) with a kernel size of 1×1×1, a stride of 1, and padding of 0, to process the input feature tensor. By directly performing a linear mapping in the spatial domain, we obtain a tensor that captures local residual features. ; Step 3.4: Convert the spectral convolution feature tensor obtained in Step 3.2.4 into... The local residual tensor obtained in step 3.3 The data is added element by element, then the data distribution is stabilized by a three-dimensional batch normalization layer (BatchNorm3d), and nonlinear processing is performed using the GELU activation function to obtain the output feature tensor of the current module. ; Step 3.5: Output the feature tensor As input to the next feature extraction module, steps 3.2 to 3.4 are repeated L times (steps 3.2 to 3.4 correspond to one Fourier layer, and there are 5 Fourier layers in this embodiment). The final output tensor is pruned to remove the boundary zero-padding components added in step 3.1, resulting in the final high-dimensional spatiotemporal evolution latent feature tensor. .

[0024] Step 4: Convert the latent feature tensor after spatiotemporal depth encoding The dual-track multi-head prediction structure is input, and the pressure prediction track and saturation prediction track are decoupled to output a pressure field and a multiphase saturation field with seepage law in parallel. Includes the following sub-steps: Step 4.1: The dual-track multi-head prediction structure will use the latent feature tensor The data are input into the pressure prediction orbital and the saturation prediction orbital, respectively; the pressure prediction orbital adopts a branch-parallel structure of data residual branch and physical prior branch. Step 4.2: In the data residual branch, the latent feature tensor The input consists of a multilayer perceptron composed of two fully connected layers and the GELU activation function; the first layer expands the number of channels from W to 128 to enhance the nonlinear representation, while the second layer compresses the number of channels to 1. This branch outputs a purely data-driven local pressure residual field tensor. Its physical meaning is the strong nonlinear pressure disturbance caused by the absolute pressure change and the local phase permeation change of the fitting system. Step 4.3: In the physical prior branch, perform analytical injection of the frequency domain spatiotemporal evolution kernel; the specific process is as follows: Step 4.3.1: Define and extract the implicitly represented virtual seepage source term q. virtual The flow control equations of real oil reservoirs are highly nonlinear. If they are directly used as loss function constraints, the permeability K and multiphase mobility will change rapidly with high frequency. This can lead to severe curses of difference differentiation and gradient conflicts. To solve this problem, we introduce the operator equivalent substitution principle, moving all nonlinear terms to the right-hand side of the equation to construct an equivalent diffusion equation: ; in, As the reference pressure field; The overall pressure conductivity coefficient of an ideal medium; For the Laplace operator; In this embodiment, a single-layer fully connected network is used to perform channel compression on the latent feature tensor, outputting a single-channel virtual percolation source term. The physical function of this term is to implicitly absorb and equivalently transform complex formation heterogeneous disturbances and multiphase fluid disturbances using neural networks, converting them all into local source-sink disturbances for ideal isotropic media, thereby solving nonlinear disturbances for subsequent application of linear physical operators. Step 4.3.2: If directly solving the above diffusion equation in the spatiotemporal domain requires expensive three-dimensional spatiotemporal convolution integration, to improve computational efficiency, the partial differential operation is transformed into frequency domain algebra operation using Fourier transform. Simultaneously, to prevent the periodicity of the discrete Fourier transform from causing future time-instance pressures to violate physical laws, a zero-causal padding of length equal to the prediction time step is first applied to the end of the time axis of the virtual source term tensor. Then, a three-dimensional fast real Fourier transform (3DrFFT) is performed to obtain the frequency domain source term tensor. : ; in, Represents the three-dimensional fast real Fourier transform; They are respectively Discrete spatial wavenumber of direction; The discrete angular frequency corresponding to the time dimension; Step 4.3.3: Construct a discrete two-dimensional spatial wavenumber grid and a one-dimensional time angular frequency grid, and define a spatiotemporal evolution kernel carrying the dual physical properties of time delay and spatial radial diffusion. This kernel function is the equivalent analytical operator of the above partial differential equation in the frequency domain: ; in, As the core of spatiotemporal evolution; The imaginary unit; For the set small constant (such as This kernel is used to prevent division-by-zero singularity. This kernel function will serve as a fixed physical prior operator, directly participating in the forward propagation computation of the subsequent network. Step 4.3.4: Obtain the frequency domain virtual source term tensor With frequency domain spatiotemporal evolution kernel Subsequently, according to the convolution theorem of Fourier transform, the global convolution operation in the spatiotemporal domain is strictly equivalent to element-wise algebraic multiplication in the frequency domain. Therefore, this invention performs a complex element-wise multiplication operation between the frequency domain source term tensor and the frequency domain spatiotemporal evolution kernel to solve for the frequency domain reference pressure field: ; in, The frequency domain reference pressure field; This represents the element-wise multiplication operation of a tensor with complex numbers. Step 4.3.5: For the frequency domain reference pressure field... Applying an inverse three-dimensional real Fourier transform to map back to the physical spacetime domain, and then trimming away the filling portion introduced in step 4.3.2 along the time axis, a reference pressure field with strict causality, time lag effect, and spatial diffusion physical properties is obtained. ; Step 4.4: Output the local pressure residual field tensor from the data residual branch. The reference pressure field output by the physical prior branch To achieve fusion, learnable dynamic equilibrium weight parameters are introduced. The output is a final predicted pressure field that combines the laws of physical spatiotemporal evolution with local accuracy. ; in, The pressure residual field is output by the data residual branch; The reference pressure field output by the physical prior branch; This represents the final pressure field after the merger. These are learnable, dynamically balanced weights. Step 4.5: In the saturation prediction branch, a multi-head parallel prediction structure is adopted for the latent feature tensor. To perform mapping decoding, first convert the latent feature tensor... The input is processed by a multilayer perceptron (MLP). The MLP employs a two-layer fully connected network structure. The first layer expands the number of channels to 128 and performs a nonlinear transformation using GELU activation. The second layer reduces the number of channels to 3. Based on the output of the MLP, a three-phase saturation data residual field is constructed. Finally, the gas saturation prediction field for future time steps is output in parallel. Oil saturation prediction field and water saturation prediction field ; Step 4.6: Concatenate the pressure field tensor predicted in Step 4.4 with the three-phase saturation field tensor predicted in Step 4.5 in the channel dimension to obtain the final output multiphase flow spatiotemporal dynamic field prediction tensor with 4 channels.

[0025] Step 5: During network training, construct a joint loss function that includes data constraints, gradient constraints, initial condition constraints, and boundary condition constraints; Includes the following sub-steps: Step 5.1: For any physical field variable predicted by the network... With the corresponding real label variable Calculate the composite relative norm error, which includes fundamental physical quantities and spatiotemporal gradients. : ; in, Represents any real physical field (pressure or saturation); This represents the corresponding prediction field; Represents the L2 norm of a tensor; Indicates in A comprehensive difference gradient operator in dimensionality; These are the gradient regularization weight coefficients; To prevent tiny constants with a denominator of zero; The loss was calculated separately for the pressure field, gas-bearing field, oil-bearing field, and water-saturation field, and then weighted and summed to obtain the global data principal loss. : ; in, These represent the actual pressure, gas content, oil content, and water saturation fields, respectively. For the corresponding prediction field; It is a dedicated weighting coefficient for gas saturation, used to balance the dimensional deviations caused by the highly compressible gas phase; Step 5.2: Calculate the physical boundary condition loss. By employing the Neumann closed boundary condition, i.e., the fluid flow rate at the outer boundary of the reservoir is zero, the spatial first derivative of the prediction field at the four boundaries is extracted to make its mean square error approach zero. ; in, The normal partial derivative operator is perpendicular to the boundary. This represents the total number of grid nodes involved in the boundary condition calculation. To predict physical fields At the outer boundary of the reservoir The values ​​at each grid node; Step 5.3: Calculate the initial condition loss. Extract the initial time step data from the prediction result sequence and calculate its mean square error with the true initial physical field: ; in, The actual physical field (pressure, saturation) at the initial time step; The physical fields (pressure, saturation) predicted by the network at the initial time step. Step 5.4: Construct the final total loss function The above losses are then combined using a linear weighted average: ; in, and These are the weight coefficients for the boundary conditions and the initial conditions, respectively. This loss function utilizes data constraints, gradient constraints, boundary condition constraints, and initial condition constraints to achieve stable and efficient parameter updates.

[0026] Step 6: Utilize the trained Fourier neural operator models to perform deduction, and execute the three-dimensional dynamic field reconstruction and output in all spatiotemporal dimensions; Includes the following sub-steps: Step 6.1: After slicing the 3D data in the test set into layers, input them in parallel into the corresponding pre-trained 2D models of each layer; the models automatically perform spatiotemporal feature extraction and physical dual-track fusion, directly inferring and generating the future of each layer. Prediction results of pressure field and three-phase saturation field at each time step; Step 6.2: Perform denormalization on the dimensionless prediction tensor output by the network to restore it to its true physical dimensions. Then, in the spatial dimension, [the following steps are performed]. The two-dimensional slices are spliced ​​together along the depth Z-axis to reassemble into a complete shape. The three-dimensional multiphase flow dynamic field, The number of samples in the test set.

[0027] To demonstrate the feasibility and superiority of this invention, a 96×96×4 three-dimensional three-phase heterogeneous reservoir dataset was established for verification. A numerical simulator was used to simulate development for 1500 days. A deep learning dataset was then built based on the dynamic production data and static geological data after development. The initial dataset contained 500 independent samples, randomly divided into 425 training samples and 75 test samples. Each time step was 30 days, with a total of 50 time steps per sample (1500 days). The input time steps were selected... =10, predict future time steps =40. During the model training phase, the initial learning rate was set to 0.001, the optimizer was Adam, and step learning rate decay (StepLR, decaying by 0.85 every 20 epochs) was configured. The maximum number of training iterations was set to 500.

[0028] Based on the generated dataset, the model was trained for 500 rounds. After training, the best model was saved. Cases were randomly selected from the test set for testing. The predicted results for pressure, oil saturation, gas saturation, and water saturation fields on day 1500 were compared with the actual results, and the absolute errors were analyzed as follows: Figures 5-8 As shown, by Figure 5 It can be seen that the overall accuracy of the three-dimensional full-field pressure prediction reaches 99.5%, and the accuracy of the near-wellbore local area reaches 96%, indicating a good overall prediction effect; Figure 6 It can be seen that the overall prediction accuracy of the three-dimensional full-field oil saturation reaches 99.6%, and the local accuracy in the near-wellbore area reaches 95%, accurately depicting the complex characteristics of the oil-water front and the near-wellbore area; Figure 7 It can be seen that the overall prediction accuracy of gas saturation in the three-dimensional field reaches 99.6%, and the local accuracy in the near-wellbore area reaches 94%, accurately capturing the drastic changes in gas saturation caused by rapid diffusion; from Figure 8It can be seen that the overall prediction accuracy of the three-dimensional full-field water saturation reaches 99.7%, and the local accuracy near the well reaches 98%. The maximum error near the well remains at an extremely low level, accurately predicting local flow characteristics. Experimental results show that the model training time is extremely short and no gradient collapse occurs. In long-term (1500 days) prediction, the three-dimensional full-field pressure and three-phase saturation distribution predicted by the method of this invention are in high agreement with the high-fidelity calculation results of commercial numerical simulators, and the relative error of the entire field is controlled at an extremely low level. While avoiding the solution of complex three-phase physical equations, this invention reduces the entire prediction time from several hours in traditional simulators to seconds, greatly improving the analysis and prediction efficiency of multiphase flow reservoirs.

[0029] 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 method for intelligent dynamic analysis and prediction of oil reservoirs based on physical-guided neural operators, characterized in that, Includes the following steps: Step 1: Generate a large number of heterogeneous reservoir numerical models based on sequential Gaussian, and obtain time-series dynamic data of pressure field and three-phase saturation field through production simulation to form a reservoir dynamic prediction dataset. Step 2: Divide the 3D reservoir grid in the dataset into multiple 2D grid slices along the depth Z-axis, normalize and stitch the features of each slice, and construct the multiphase flow feature tensor corresponding to each slice; Step 3: For the multiphase flow feature tensor of each two-dimensional mesh slice, construct a Fourier neural operator model with embedded frequency domain analytical prior. The model includes a multi-layer spatiotemporal feature extraction module and a dual-track multi-head prediction structure. First, extract the spatiotemporal evolution features of the slice through the multi-layer spatiotemporal feature extraction module to obtain the latent feature tensor. Step 4: Input the latent feature tensor into the dual-track multi-head prediction structure, and output the pressure field and multiphase saturation field with seepage law in parallel by decoupling the pressure prediction track and the saturation prediction track. Step 5: Construct a joint loss function for data constraints, gradient constraints, initial condition constraints, and boundary condition constraints, and train and optimize the model. Step 6: Utilize the trained Fourier neural operator models to perform deduction, and execute the three-dimensional dynamic field reconstruction and output in all spatiotemporal dimensions; Step 4 includes the following sub-steps: Step 4.1: The dual-track multi-head prediction structure will use the latent feature tensor The data are input into the pressure prediction orbital and the saturation prediction orbital, respectively; the pressure prediction orbital adopts a dual-branch parallel structure with a data residual branch and a physical prior branch. Step 4.2: In the data residual branch, the latent feature tensor The input consists of a multilayer perceptron composed of two fully connected layers and the GELU activation function, and the output is a purely data-driven local pressure residual field tensor. ; Step 4.3: In the physical prior branch, perform analytical injection of the frequency domain spatiotemporal evolution kernel; the specific process is as follows: Step 4.3.1: Use a single-layer fully connected network to perform channel compression on the latent feature tensor, and output a single-channel virtual percolation source term. This is used to equivalently characterize formation heterogeneous disturbances and multiphase fluid interference; by introducing the operator equivalent substitution principle, an equivalent diffusion equation is constructed: ； in, As the reference pressure field; The overall pressure conductivity coefficient of an ideal medium; For the Laplace operator; Step 4.3.2: First, apply zero-causal padding of length equal to the prediction time step to the end of the time axis of the virtual source term tensor. Then, perform a three-dimensional fast real Fourier transform to obtain the frequency domain source term tensor. : ； Step 4.3.3: Construct a discrete two-dimensional spatial wavenumber grid and a one-dimensional time angular frequency grid, and define a spatiotemporal evolution kernel carrying the dual physical properties of time delay and spatial radial diffusion: ； in, As the core of spatiotemporal evolution; The imaginary unit; A small constant that is set; Step 4.3.4: Perform a complex element-wise product operation between the frequency domain source term tensor and the frequency domain spatiotemporal evolution kernel to solve for the frequency domain reference pressure field: ； in, The frequency domain reference pressure field; This represents the element-wise multiplication operation of a tensor with complex numbers. Step 4.3.5: For the frequency domain reference pressure field... An inverse three-dimensional real Fourier transform is applied to map back to the physical spacetime domain, and then the filling portion introduced in step 4.3.2 is trimmed along the time axis to obtain the reference pressure field. ; Step 4.4: Output the pressure residual field from the data residual branch. The reference pressure field output by the physical prior branch To achieve fusion, learnable dynamic equilibrium weight parameters are introduced. Output the final predicted pressure field: ； in, The pressure residual field is output by the data residual branch; The reference pressure field output by the physical prior branch; This represents the final pressure field after the merger. These are learnable, dynamically balanced weights. Step 4.5: In the saturation prediction trajectory, a multi-head parallel prediction structure is adopted, using a two-layer fully connected network to process the latent feature tensor. Mapping and decoding are performed; the first layer expands the number of channels to 128 and activates them using GELU, the second layer reduces the number of channels to 3, and outputs the gas saturation prediction field for future time steps in parallel. Oil saturation prediction field and water saturation prediction field ; Step 4.6: Concatenate the pressure field tensor predicted in Step 4.4 with the three-phase saturation field tensor predicted in Step 4.5 in the channel dimension to obtain the final output multiphase flow spatiotemporal dynamic field prediction tensor with 4 channels.

2. The method for intelligent dynamic analysis and prediction of oil reservoirs based on physical guided neural operators according to claim 1, characterized in that, Step 1 includes the following sub-steps: Step 1.1: Generate a batch random heterogeneous permeability field and a heterogeneous porosity field using the sequential Gaussian simulation method; Step 1.2: Establish a three-dimensional three-phase reservoir numerical benchmark model, configure the grid size and reservoir physical property parameters, set the injection and production well pattern and well control conditions, and import the generated heterogeneous porosity field into the benchmark model; Step 1.3: Import each generated heterogeneous permeability field into the reservoir numerical benchmark model one by one, replacing the initial homogeneous permeability field to form multiple heterogeneous reservoir numerical models; perform production simulation on each model, extract the time-series dynamic data of pressure field and oil, gas and water three-phase saturation field, and save them together with the corresponding permeability field to form a reservoir dynamic prediction dataset.

3. The method for intelligent dynamic analysis and prediction of oil reservoirs based on physical guided neural operators according to claim 1, characterized in that, Step 2 includes the following sub-steps: Step 2.1: Obtain static geological features and dynamic temporal features from the dataset on the 3D reservoir grid. Static geological features include absolute permeability features in the X, Y, and Z directions, while dynamic temporal features include historical features of pressure, gas saturation, oil saturation, and water saturation. The data is then divided into three parts with a size of [missing information - likely a specific value]. The three-dimensional reservoir grid is divided into layers along the depth Z-axis. Each size is Two-dimensional mesh slices; Step 2.2: For each 2D mesh slice, extract the previous... The dynamic time series characteristics of each historical time step are obtained, and Gaussian normalization is used to perform dimensionless processing on static and dynamic data respectively. Step 2.3: Concatenate the normalized data along the channel dimension, including the absolute penetration features of the three channels, and each channel contains... The historical characteristics of pressure, gas saturation, oil saturation, and water saturation in each channel are used to generate a high-dimensional input feature tensor of the physical field as the feature tensor of the multiphase flow.

4. The method for intelligent dynamic analysis and prediction of oil reservoirs based on physical guided neural operators according to claim 1, characterized in that, Step 3 includes the following sub-steps: Step 3.1: Concatenate the multiphase flow feature tensor obtained in Step 2 with a three-dimensional spatial and temporal relative coordinate grid tensor in the channel dimension; then, through a fully connected mapping layer, map the number of channels of the concatenated tensor to the preset hidden layer feature dimension W. Applying zero-padding to the tensor along the spatiotemporal dimension yields the initial latent feature tensor. ; Step 3.2: Initialize the hidden feature tensor The input consists of a stacked multi-layer spatiotemporal feature extraction module. Within each module, a two-branch parallel computation is employed, with the first branch being a three-dimensional frequency domain spectral convolution branch. The specific operations are as follows: Step 3.2.1: Process the input feature tensor of the current layer. By applying a three-dimensional fast real Fourier transform to the spatial domain (X, Y) and the time domain (T), the frequency domain feature tensor is obtained. : ； in, Represents the three-dimensional fast real Fourier transform; They are respectively Frequency wavenumber of dimension for Discrete angular frequencies corresponding to the time dimension; Step 3.2.2: Perform mode truncation in the frequency domain, retaining only the low-frequency mode components in three dimensions; Step 3.2.3: Perform cross-channel tensor multiplication on the channel dimension of the truncated frequency domain feature tensor and the complex weight tensor based on the Einstein summation convention: ； in, For the frequency domain complex intermediate feature tensor, For the truncated feature tensor, This represents a cross-channel tensor multiplication operation based on Einstein's summation convention. It is a learnable complex weight matrix; Step 3.2.4, for The truncated high-frequency components are padded, and then an inverse three-dimensional real Fourier transform is applied to map them back to the physical spatiotemporal domain, yielding the spectral convolution feature tensor. ; Step 3.3: The second branch is a spatial domain local convolution branch, using a 3D convolution kernel with a kernel size of 1×1×1, a stride of 1, and padding of 0, to process the input feature tensor. By directly performing a linear mapping in the spatial domain, we obtain a tensor that captures local residual features. ; Step 3.4: Convert the spectral convolution feature tensor obtained in Step 3.2.4 into... The local residual tensor obtained in step 3.3 The data is added element-wise, then the data distribution is stabilized through a three-dimensional batch normalization layer, and nonlinear processing is performed using the GELU activation function to obtain the output feature tensor of the current module. ; Step 3.5: Output the feature tensor As input to the next layer feature extraction module, steps 3.2 to 3.4 are repeated L times to prune the final output tensor, removing the boundary zero-padding components added in step 3.1, thus obtaining the final high-dimensional spatiotemporal evolution latent feature tensor. .

5. The method for intelligent dynamic analysis and prediction of oil reservoirs based on physical guided neural operators according to claim 1, characterized in that, Step 5 includes the following sub-steps: Step 5.1: For any physical field variable predicted by the network and its corresponding true label variable, calculate the composite relative norm error, which includes the fundamental physical quantity and the spatiotemporal gradient. : ； in, Represents any real physical field; This represents the corresponding prediction field; Represents the L2 norm of a tensor; Indicates in A comprehensive difference gradient operator in dimensionality; These are the gradient regularization weight coefficients; To prevent constants with a denominator of zero; The loss was calculated separately for the pressure field, gas-bearing field, oil-bearing field, and water-saturation field, and then weighted and summed to obtain the global data principal loss. : ； in, These represent the actual pressure, gas content, oil content, and water saturation fields, respectively. For the corresponding prediction field; This is a weighting coefficient specific to gas saturation. Step 5.2: Calculate the physical boundary condition loss. : ； in, The normal partial derivative operator is perpendicular to the boundary. This represents the total number of grid nodes involved in the boundary condition calculation. To predict physical fields At the outer boundary of the reservoir The values ​​at each grid node; Step 5.3: Calculate the initial condition loss. : ； in, The actual physical field at the initial time step; The physical field predicted by the network at the initial time step; Step 5.4: Construct the final total loss function : ； in, and These are the weighting coefficients for the boundary conditions and the initial conditions, respectively.

6. The method for intelligent dynamic analysis and prediction of oil reservoirs based on physical guided neural operators according to claim 1, characterized in that, Step 6 includes the following sub-steps: Step 6.1: After slicing the 3D data in the test set into layers, input them in parallel into the corresponding pre-trained 2D models of each layer; the models automatically perform spatiotemporal feature extraction and physical dual-track fusion, directly inferring and generating the future of each layer. Prediction results of pressure field and three-phase saturation field at each time step; Step 6.2: Perform an inverse normalization operation on the dimensionless prediction tensor output by the network to restore it to the true physical dimensions; Subsequently, in the spatial dimension, The two-dimensional slices are spliced ​​together along the depth Z-axis to reassemble into a complete shape. The three-dimensional multiphase flow dynamic field, The number of samples in the test set.

Images