A method for predicting reservoir production dynamics by fusing discrete gradient information

Abstract

This invention discloses a method for dynamic prediction of reservoir production that integrates discrete gradient information, belonging to the interdisciplinary field of reservoir development and artificial intelligence. The steps are as follows: First, a numerical simulation dataset for oil-water two-phase flow in heterogeneous reservoirs is constructed based on numerical simulation methods. Second, a dual-branch network structure is designed to extract the spatial and physical features of the input field data in parallel, including a main feature encoding branch and a difference operator branch. Third, a backbone network is designed for deep nonlinear modeling. Fourth, an efficient pressure and saturation field prediction neural network model is constructed based on the dual-branch network structure and the backbone network. During the model training phase, spatial observation points from some time steps are used to participate in the calculation of data loss terms, while the residuals of the physical control equations are introduced as physical loss terms at all time steps and throughout the entire space. Finally, a well-trained efficient pressure and saturation field prediction neural network model is obtained, achieving high-precision prediction of pressure and saturation fields across the entire time series.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of reservoir development and artificial intelligence, and specifically relates to a method for dynamic prediction of reservoir production by integrating discrete gradient information. Background Technology

[0002] Traditional reservoir production dynamics prediction typically relies on numerical simulation methods. However, when solving multiphase flow problems in highly heterogeneous reservoirs, these methods often suffer from high computational costs and poor convergence, particularly limiting their applicability in rapid multi-scenario extrapolation and real-time optimization. While purely data-driven methods have improved computational efficiency to some extent in recent years, their predictions often lack physical consistency, making it difficult to guarantee the reliability of extrapolated predictions, and they usually require a large amount of labeled data for training. Therefore, how to integrate physical mechanisms with the advantages of data-driven methods to enhance the model's generalization ability and physical interpretability while ensuring efficiency is a pressing technical problem to be solved in this field. Summary of the Invention

[0003] To address the aforementioned problems in existing technologies, this invention proposes a reservoir production dynamic prediction method that integrates discrete gradient information. The method is rationally designed, overcomes the shortcomings of existing technologies, and has good results.

[0004] To achieve the above-mentioned objectives, the following technical solution is adopted:

[0005] A method for predicting reservoir production dynamics by fusing discrete gradient information includes the following steps:

[0006] Step 1: Construct a numerical simulation dataset for oil-water two-phase flow in heterogeneous reservoirs based on numerical simulation methods;

[0007] Step 2: Design a dual-branch network structure to extract the spatial and physical features of the input field data in parallel, and achieve effective feature fusion. The dual-branch network structure includes a main feature encoding branch and a difference operator branch.

[0008] Step 3: Design a backbone network for deep nonlinear modeling to capture the evolution of field data. The backbone network includes a Patch Embedding module, a backbone Transformer network, and a decoder.

[0009] Step 4: Construct an efficient neural network model for predicting stress and saturation fields based on a dual-branch network structure and a backbone network. During the model training phase, spatial observation points from some time steps are used to participate in the calculation of data loss. At the same time, the residuals of the physical control equations are introduced as physical loss terms in all time steps and the entire space.

[0010] Step 5: Obtain the trained high-efficiency neural network model for predicting pressure and saturation fields, and achieve high-precision prediction of pressure and saturation fields throughout the entire time series.

[0011] Further, step 1 includes the following sub-steps:

[0012] Step 1.1: Establish a two-dimensional oil-water two-phase reservoir model using a numerical simulator, and set different reservoir parameter attribute data, including: grid size, reservoir thickness, initial reservoir pressure, initial water saturation, rock compressibility coefficient, oil and water phase density, oil and water phase volume coefficient, oil and water phase viscosity, and oil-water relative permeability curves, and set capillary force to be ignored.

[0013] Step 1.2: Using the sequential Gaussian simulation method, combined with the Gaussian kernel function and related length parameters, covariance modeling is performed on spatial points to generate a logarithmic permeability field with spatial heterogeneity. The logarithmic permeability field is then transformed by exponential transformation to obtain a physically reasonable positive permeability field.

[0014] Subsequently, based on empirical formulas, a linear transformation is performed on the permeability field to generate a porosity field associated with the permeability distribution, expressed as:

[0015] (1)

[0016] in, Porosity; For penetration rate; , These are empirical parameters;

[0017] Convert the permeability field and porosity field data according to the numerical simulator specifications and save them as inc format files;

[0018] Step 1.3: After importing the permeability field and porosity field data into the numerical simulator, set the location and well control conditions of the injection well and production well, and then run the numerical simulator to obtain the complete data sequence of the pressure field and saturation field evolution over time. Organize the results into a tensor structure [C, T, H, W] and save it as a mat format data, where C represents the number of data channels, T represents the number of time steps, and H and W represent the height and width of the spatial dimension, respectively.

[0019] Furthermore, step 2 includes the following sub-steps:

[0020] Step 2.1: Read the .mat and .inc data files, and concatenate the initial pressure and saturation field data, permeability field data, and porosity field data along the channel dimension to form an input tensor with the structure [4, 100, 100]. The MinMaxScaler normalization method is used to normalize the data from different channels.

[0021] Step 2.2: The spatial and physical features of the input field data are extracted in parallel using a dual-branch network structure, which includes a main feature encoding branch and a difference operator branch.

[0022] The main feature encoding branch extracts spatial features of pressure, saturation, and geological parameter fields through multi-layer convolution and pooling operations, and reduces the dimensionality of the features to a low-dimensional space to enhance the representation of the main features. Then, the transposed convolution method is used to restore the features to the original input spatial size, and the output spatial feature tensor is output. ;

[0023] The difference operator branch uses weighted, fixed convolution kernels to perform first- and second-order difference operations on the input data, extracting physics-related features including physical gradients and Laplace operators, effectively encoding the physical mechanisms, and outputting a physical constraint feature tensor with the same size as the original input. ;

[0024] Step 2.3: Convert the input tensor The spatial feature tensor output by the main feature encoding branch Physical constraint feature tensor output by the difference operator branch The components are concatenated along the channel dimension to obtain the fused feature tensor. ; then on Applying a 1×1 convolution operation, combined with BatchNorm and ReLU, achieves deep fusion and channel dimensionality reduction of dual-branch features, expressed as:

[0025] (3)

[0026] in, To fuse the feature tensors after dimensionality reduction, These are the weights of the convolutional layer; The bias vector of the convolutional layer;

[0027] Introducing a simplified channel attention mechanism, firstly... Global average pooling is performed along the spatial dimension to obtain a global description of each channel. ,Right now:

[0028] (4)

[0029] in, For the first A global description of each channel; They are respectively Height and width; To fuse the feature map The spatial location of each channel The value;

[0030] Subsequently, the channel weights are learned through two fully connected layers and a sigmoid activation function, expressed as:

[0031] (5)

[0032] in, For the first Adaptive attention weights for each channel, This represents the Sigmoid activation function. , These are the weights of the two fully connected layers; , These are the bias vectors for the two fully connected layers, respectively.

[0033] Finally, the attention weights This is applied to each channel of the fused feature tensor to achieve adaptive feature enhancement, expressed as:

[0034] (6)

[0035] in, For the final output feature map, the first... The spatial location of each channel The value of .

[0036] Furthermore, the specific process of the main feature encoding branch is as follows:

[0037] The input tensor obtained in step 2.1 Input the main feature encoding branch to prepare for spatial feature extraction;

[0038] right A first 2D convolution operation is applied, with a kernel size of 3×3, a stride of 1, and padding of 1, followed by non-linear processing using the ReLU activation function. Next, a 2×2 max pooling operation is performed with a stride of 2 to compress the spatial dimension, reducing the data feature space size to half its original size. A second 2D convolution operation is then applied to the pooled feature tensor, with a kernel size of 3×3, a stride of 1, and padding of 1, again using the ReLU activation function for non-linear processing, resulting in a compressed feature tensor. At this point, the data feature space size is H / 2 × W / 2.

[0039] right A transposed convolution operation is applied to upsample the feature tensor. The transposed convolution kernel size is 2×2, the stride is 2, and the padding is 0. The number of output channels is consistent with the number of channels after the first convolution. The transposed convolution operation restores the feature tensor spatial size from H / 2×W / 2 to the original size H×W. A 3×3 convolution operation is applied to the transposed convolution result with a stride of 1 and padding of 1 to obtain the feature data. The Sigmoid activation function is used to process the feature data, and the output is a spatial feature tensor with the same spatial size as the original input. .

[0040] Furthermore, the specific process of the difference operator branching is as follows:

[0041] The input tensor obtained in step 2.1 The input difference operator branch is used to extract physical constraint features; the boundaries of the input tensor are copied and filled according to the convolution kernel size to expand the spatial size of the input tensor.

[0042] Discrete approximation of the first derivative of the padded input tensor using a first-order difference convolution kernel. That is, the spatial first-order gradient feature tensor of the pressure field with the same size as the original input. The kernel size is 3×3, the stride is 1, and the padding is 0.

[0043] for The reservoir model uses a fourth-order precision 9-point difference convolution kernel to calculate the second derivative. That is, the spatial second derivative characteristic tensor of the pressure field with the same size as the original input. The convolution kernel size is 5×5, the stride is 1, and the padding is 0. The expression is:

[0044] (7)

[0045] (8)

[0046] in, and These represent the grid sizes in the x and y directions, respectively. coordinates Grid pressure; and These represent the grid sizes in the x and y directions, respectively.

[0047] for The reservoir model employs a 25-point finite difference convolution kernel to convolve the filled pressure field, and then uses a more refined fourth-order precision 25-point hybrid grid difference calculation to approximate the spatial second derivative. That is, a higher-precision pressure field spatial second derivative characteristic tensor The convolution kernel size is 5×5, and the stride is 1; the specific formula is as follows:

[0048] (9)

[0049] in, These are the weighting coefficients. , ;

[0050] Based on Darcy's seepage equation The coupling relationship, for the output results of the above physical operators, will or With liquid phase flow rate and input tensors The penetration rate in or Perform calculations to obtain the characteristic tensor with strong physical constraints. ;wherein liquid phase mobility The relative permeability tensor is calculated from relative permeability, volume index, and viscosity. It is obtained by interpolation based on the relative permeability curve, with the volume index and viscosity set as constants. The specific formula is as follows:

[0051] (10)

[0052] in, , These are the relative permeabilities of the aqueous and oil phases, respectively. , These are the volume factors for the aqueous phase and the oil phase, respectively. , These are the viscosities of the aqueous phase and the oil phase, respectively.

[0053] Will , or ,as well as The data is stitched and fused along the channel dimension to generate a physical constraint feature tensor with the same size as the input space. .

[0054] Furthermore, step 3 includes the following sub-steps:

[0055] Step 3.1: [The text appears to be incomplete and contains several grammatical errors. A more accurate translation would require The Patch Embedding module takes a 5-kernel, 5-stride convolution operation as input and divides the 2D spatial features into a sequence of N tokens, each with dimension d. A learnable positional encoding is then added to each token, resulting in a position-encoded sequence. ;

[0056] Step 3.2: Sequence The backbone Transformer network is fed into a stacked multi-layered encoding module. Each encoding module includes a multi-head self-attention mechanism and a feedforward neural network. Specifically, the query, key, and value are first obtained through linear transformation and divided into h heads. At each head, global dynamic weighted aggregation is used to achieve long-distance dependent feature interaction, as shown in the formula:

[0057] (11)

[0058] in, These are the query, key, and value, used for matching and weighting in the self-attention mechanism. This is the transpose of the key matrix; For the number of bullish self-attention heads; Feature dimensions for each attention head; Indicates attention output;

[0059] Each encoding module will The input sequence is obtained by residual connection and LayerNorm normalization. Then the sequence Input to feedforward neural network and sequence The residual connections are then processed by LayerNorm to obtain the output sequence of a single encoding module. This process is repeated for multiple encoding modules, with the output of each module serving as the input of the next module, ultimately resulting in a highly integrated token sequence.

[0060] Step 3.3: The highly integrated token sequence is restored to a low-resolution spatial feature tensor through linear transformation and reshape. The original field size is restored through deconvolution operation with a kernel size of 5 and a stride of 5. After 3×3 convolution and ReLU refinement, the refined feature tensor is output and connected to the pressure branch and saturation branch respectively. The pressure branch uses 1×1 convolution to output the pressure field prediction of the next time step. The saturation branch uses 1×1 convolution followed by a Sigmoid activation function to output the saturation field prediction of the next time step with the value limited to the interval [0,1].

[0061] Furthermore, step 4, training the model, includes the following sub-steps:

[0062] Step 4.1: Starting from the initial field state, the model sequentially uses the prediction result at each time step as the input for the next time step, continuously rolling and recursively, and finally outputs a complete multi-step physical field dynamic prediction sequence;

[0063] Step 4.2: During recursive prediction, for certain moments within a specified time step range, the real physical field data of the corresponding moment is used to replace the model's prediction result for that moment as the network input for the next time step. For other moments, the model's autoregressive output is used. By gradually reducing the probability of using real physical field data for replacement, the model relies more on its own prediction in the later stages of training.

[0064] Step 4.3: Data Loss Items It is computed only at a specified training time step and a partial set of spatial sampling points. The local spatial sampling points, expanded around the well center at the specified training time step, are selected, and the expression is:

[0065] (12)

[0066] in, They are respectively Time of the first Predicted and actual values ​​for each spatial sampling point The number of spatial sampling points at that time. To calculate the number of time steps for data loss;

[0067] Physical governing equations Constraints are applied across all time steps and the entire space grid to ensure the overall dynamic physical rationality and conservation. The physical governing equations are derived from the oil-water two-phase flow equations. The mathematical equations for the oil-water two-phase flow are as follows:

[0068] (13)

[0069] (14)

[0070] in, For absolute penetration rate, and These are the relative permeabilities of the oil phase and the water phase, respectively. and The viscosity of the oil phase and the viscosity of the water phase are respectively. and These are the formation volume factors for the oil and water phases, respectively. and The pressures of the oil phase and the water phase are respectively. and These are the fluid weights of the oil phase and the water phase, respectively. For spatial height gradient, Porosity and These represent the saturation levels of the oil and water phases, respectively. and These are the formation volume factors for the oil and water phases, respectively. and These represent the volumetric flow rates of the oil and water phases under standard ground conditions, respectively. For point source functions, Indicates the first The coordinates of each well point;

[0071] For the oil-water two-phase flow equation, where permeability is isotropic and gravity and capillary forces are neglected, the volume integral formula for the aqueous phase flow equation is as follows:

[0072] (15)

[0073] in, For the unit volume size, The relative permeability of the aqueous phase;

[0074] Discretization formula (15), after rearranging terms, yields the physical information residual function. As shown below:

[0075] (16)

[0076] in, For the first The volume of each grid, For time step, Use the The grid, the first Each time step , and calculate, For the reservoir model The first grid and its adjacent first grid The water phase conductivity between grids is the product of the flow coefficient and the geometric factor, where the geometric factor is the ratio of the grid contact surface to the center point spacing, and its expression is:

[0077] (17)

[0078] in, For the first The first grid and its adjacent first grid The contact area of ​​each grid, For the first The first grid and its adjacent first grid The center point spacing of each grid Use the The first grid and its adjacent first grid In each grid, the upstream grid properties are arranged according to the pressure difference direction. , , and calculate;

[0079] The physical loss term calculates the mean square error of the physical information residuals across all time steps and the entire spatial grid, and its expression is:

[0080] (18)

[0081] in, for Time of the first Predicted values ​​for each grid point The total number of grid points in a single time step. To predict the number of data points in a dataset;

[0082] Step 4.4: The total loss function is designed as a weighted combination of data item loss and physical equation residual loss, expressed as:

[0083] (19)

[0084] in, Total loss; The weights are for the data loss items.

[0085] The beneficial technical effects of this invention are as follows:

[0086] This invention proposes a novel neural network model that integrates discrete gradient information from multiphase flow in oil reservoirs. Without relying on traditional numerical simulations, it achieves rapid and high-precision prediction of pressure and saturation fields in heterogeneous reservoirs, significantly improving the efficiency and accuracy of reservoir dynamic evolution simulation. This method extracts spatial distribution features and physical gradient features in parallel through a dual-branch network, embedding physical constraints from the flow equations into a deep learning process, effectively enhancing the model's physical consistency and extrapolation stability. Simultaneously, it introduces a channel attention mechanism and a feature fusion layer to enhance the expressive power of multi-source data features. The model utilizes convolution blocks and positional encoding to achieve sequential representation of spatial information, combining a Transformer backbone network to capture long-range spatiotemporal dependencies in the reservoir field, and reconstructs the output pressure and saturation fields through deconvolution and convolution refinement. Compared to traditional numerical simulation methods, this invention significantly improves computational efficiency, greatly reduces data dependency, and possesses both physical interpretability and generalization ability, providing a reliable and efficient new technical approach for dynamic prediction of multiphase flow in complex oil reservoirs. Attached Figure Description

[0087] Figure 1 This is a flowchart of a reservoir production dynamic prediction method that integrates discrete gradient information in this invention.

[0088] Figure 2 This is a flowchart of the neural network model calculation in this invention.

[0089] Figure 3 This is a diagram of the differential convolution kernel weight parameters used in this invention;

[0090] Among them, (a) is the weight parameter diagram of the first-order differential convolution kernel; (b) is the weight parameter diagram of the 9-point differential convolution kernel; and (c) is the weight parameter diagram of the 25-point differential convolution kernel.

[0091] Figure 4 This is a reservoir permeability and porosity diagram used in the embodiments of the present invention;

[0092] (a) is a permeability diagram; (b) is a porosity diagram.

[0093] Figure 5 These are the actual pressure field map, predicted pressure field map, and relative error map for some time steps in this embodiment of the invention;

[0094] Among them, (a), (b) and (c) are the actual pressure field map, the predicted pressure field map and the relative error map when T=50, respectively; (d), (e) and (f) are the actual pressure field map, the predicted pressure field map and the relative error map when T=50, respectively; (g), (h) and (i) are the actual pressure field map, the predicted pressure field map and the relative error map when T=100, respectively.

[0095] Figure 6 These are the actual water saturation field map, predicted water saturation map, and relative error map for some time steps in this embodiment of the invention;

[0096] Among them, (a), (b) and (c) are the actual water saturation field map, the predicted water saturation map and the relative error map at T=50, respectively; (d), (e) and (f) are the actual water saturation field map, the predicted water saturation map and the relative error map at T=50, respectively; (g), (h) and (i) are the actual water saturation field map, the predicted water saturation map and the relative error map at T=100, respectively. Detailed Implementation

[0097] The specific embodiments of the present invention will be further described below with reference to specific examples:

[0098] A method for predicting reservoir production dynamics by fusing discrete gradient information, such as Figure 1 As shown, it includes the following steps:

[0099] Step 1: Construct a numerical simulation dataset for oil-water two-phase flow in heterogeneous reservoirs based on numerical simulation methods. The database contains reservoir parameter attribute data (permeability field, porosity field), well network operating system data (constant bottomhole flowing pressure of production wells, constant water injection volume of injection wells), pressure field variation data over time, and saturation field variation data over time.

[0100] Step 1 includes the following sub-steps:

[0101] Step 1.1: Establish a two-dimensional oil-water two-phase reservoir model using a numerical simulator, and set different reservoir parameter attribute data, including: grid size, reservoir thickness, initial reservoir pressure, initial water saturation, rock compressibility coefficient, oil and water phase density, oil and water phase volume coefficient, oil and water phase viscosity, and oil-water relative permeability curves, and set capillary force to be ignored.

[0102] Step 1.2: Using the sequential Gaussian simulation method, combined with the Gaussian kernel function and related length parameters, covariance modeling is performed on spatial points to generate a logarithmic permeability field with spatial heterogeneity. The logarithmic permeability field is then transformed by exponential transformation to obtain a physically reasonable positive permeability field.

[0103] Subsequently, based on empirical formulas, a linear transformation is performed on the permeability field to generate a porosity field associated with the permeability distribution, expressed as:

[0104] (1)

[0105] in, Porosity; For penetration rate; , These are empirical parameters, determined by specific data, core tests, or historical analysis, and their values ​​vary depending on the rock type and sedimentary environment in different regions.

[0106] The permeability and porosity data were converted and saved as inc format files according to the numerical simulator specifications for subsequent numerical simulation and model training.

[0107] Step 1.3: After importing the permeability field and porosity field data into the numerical simulator, set the location and well control conditions of the injection well and production well (e.g., constant injection rate for injection wells and constant bottom-hole flowing pressure for production wells). Then run the numerical simulator to obtain the complete data sequence of the pressure field and saturation field evolution over time. Organize the results into a tensor structure [C, T, H, W] and save it as a mat format data, where C represents the number of data channels (reservoir pressure, water saturation), T represents the number of time steps, and H and W represent the height and width of the spatial dimension, respectively.

[0108] Step 2: Design a dual-branch network structure to extract the spatial and physical features of the input field data in parallel, and achieve effective feature fusion. The dual-branch network structure includes a main feature encoding branch and a difference operator branch.

[0109] Step 2 includes the following sub-steps:

[0110] Step 2.1: Read the .mat and .inc data files, and concatenate the initial pressure and saturation field data, permeability field data, and porosity field data along the channel dimension to form an input tensor with the structure [4, 100, 100]. Where 4 represents the number of channels (reservoir pressure field, water saturation field, permeability field, and porosity field), and 100 represents the spatial size; the MinMaxScaler normalization method is used to normalize the data of different channels;

[0111] The expression for the MinMaxScaler normalization method is:

[0112] (2)

[0113] in, These are the parameters for each channel dimension; for The normalized value; for The minimum value; for The maximum value;

[0114] Step 2.2: The spatial and physical features of the input field data are extracted in parallel using a dual-branch network structure, which includes a main feature encoding branch and a difference operator branch.

[0115] The main feature encoding branch extracts spatial features of pressure, saturation, and geological parameter fields through multi-layer convolution and pooling operations, and reduces the dimensionality of the features to a low-dimensional space to enhance the representation of the main features. Then, the transposed convolution method is used to restore the features to the original input spatial size, and the output spatial feature tensor is output. ;

[0116] Specifically, the input tensor obtained in step 2.1 is... Input the main feature encoding branch to prepare for spatial feature extraction;

[0117] right A first 2D convolution operation is applied, with a kernel size of 3×3, a stride of 1, and padding of 1, followed by non-linear processing using the ReLU activation function. Next, a 2×2 max pooling operation is performed with a stride of 2 to compress the spatial dimension, reducing the data feature space size to half its original size. A second 2D convolution operation is then applied to the pooled feature tensor, with a kernel size of 3×3, a stride of 1, and padding of 1, again using the ReLU activation function for non-linear processing, resulting in a compressed feature tensor. At this point, the data feature space size is H / 2 × W / 2.

[0118] right A transposed convolution operation is applied to upsample the feature tensor. The transposed convolution kernel size is 2×2, the stride is 2, and the padding is 0. The number of output channels is consistent with the number of channels after the first convolution. The transposed convolution operation restores the feature tensor spatial size from H / 2×W / 2 to the original size H×W. A 3×3 convolution operation is applied to the transposed convolution result with a stride of 1 and padding of 1 to obtain the feature data. The Sigmoid activation function is used to process the feature data, and the output is a spatial feature tensor with the same spatial size as the original input. .

[0119] The difference operator branch uses weighted, fixed convolution kernels to perform first- and second-order difference operations on the input data, extracting physics-related features including physical gradients and Laplace operators, effectively encoding the physical mechanisms, and outputting a physical constraint feature tensor with the same size as the original input. ;

[0120] Specifically, the input tensor obtained in step 2.1 is... The input difference operator branch is used to extract physical constraint features; the boundaries of the input tensor are copied and filled according to the convolution kernel size to expand the spatial size of the input tensor, so as to ensure that the boundaries are effective when the physical operator is applied.

[0121] Discrete approximation of the first derivative of the padded input tensor using a first-order difference convolution kernel. That is, the spatial first-order gradient feature tensor of the pressure field with the same size as the original input. The weight parameters of the first-order difference convolution kernel are shown in the figure. Figure 3 As shown in (a), the kernel size is 3×3, the stride is 1, and the padding is 0;

[0122] for The reservoir model uses a fourth-order precision 9-point difference convolution kernel to calculate the second derivative. That is, the spatial second derivative characteristic tensor of the pressure field with the same size as the original input. The weight parameters of the 9-point difference convolution kernel are shown in the figure. Figure 3 As shown in (b), the convolution kernel size is 5×5, the stride is 1, and the padding is 0. The expression is:

[0123] (3)

[0124] (4)

[0125] in, and These represent the grid sizes in the x and y directions, respectively. coordinates Grid pressure; and These represent the grid sizes in the x and y directions, respectively.

[0126] for The reservoir model employs a 25-point finite difference convolution kernel to convolve the filled pressure field, and calculates the second derivative of the difference approximation space using a finer fourth-order precision 25-point hybrid mesh. That is, a higher-precision pressure field spatial second derivative characteristic tensor The weight parameters of the 25-point difference convolution kernel are shown in the figure below. Figure 3 As shown in (c), a higher-precision pressure field spatial second-order derivative characteristic tensor is obtained. The convolution kernel size is 5×5, and the stride is 1; the specific formula is as follows:

[0127] (5)

[0128] in, These are the weighting coefficients. , ;

[0129] Based on Darcy's seepage equation The coupling relationship, for the output results of the above physical operators, will or With liquid phase flow rate and input tensors The penetration rate in or Perform calculations to obtain the characteristic tensor with strong physical constraints. The fluid parameters include liquid phase mobility. The relative permeability tensor is calculated from the relative permeability, volume factor, and viscosity. It is obtained by interpolation based on the relative permeability curve, with the volume factor and viscosity set as constants. The specific formula is as follows:

[0130] (6)

[0131] in, , These are the relative permeabilities of the aqueous and oil phases, respectively. , These are the volume factors for the aqueous phase and the oil phase, respectively. , These are the viscosities of the aqueous phase and the oil phase, respectively.

[0132] Will , or ,as well as The data is stitched and fused along the channel dimension to generate a physical constraint feature tensor with the same size as the input space. .

[0133] Step 2.3: Convert the input tensor The spatial feature tensor output by the main feature encoding branch Physical constraint feature tensor output by the difference operator branch The components are concatenated along the channel dimension to obtain the fused feature tensor. ; then on Applying a 1×1 convolution operation, combined with BatchNorm and ReLU, achieves deep fusion and channel dimensionality reduction of dual-branch features, expressed as:

[0134] (7)

[0135] in, To fuse the feature tensors after dimensionality reduction, These are the weights of the convolutional layer; The bias vector of the convolutional layer;

[0136] Introducing a simplified channel attention mechanism, firstly... Global average pooling is performed along the spatial dimension to obtain a global description of each channel. ,Right now:

[0137] (8)

[0138] in, For the first A global description of each channel; They are respectively Height and width; To fuse the feature map The spatial location of each channel The value;

[0139] Subsequently, the channel weights are learned through two fully connected layers and a sigmoid activation function, expressed as:

[0140] (9)

[0141] in, For the first Adaptive attention weights for each channel, This represents the Sigmoid activation function. , These are the weights of the two fully connected layers; , These are the bias vectors for the two fully connected layers, respectively.

[0142] Finally, the attention weights This is applied to each channel of the fused feature tensor to achieve adaptive feature enhancement, expressed as:

[0143] (10)

[0144] in, For the final output feature map, the first... The spatial location of each channel The value of .

[0145] Through the above process, spatial feature tensors and physical constraint feature tensors are fused, and the representation of key information is adaptively adjusted through a channel attention mechanism, providing higher-quality feature inputs with greater discriminativeness and physical consistency for the subsequent backbone network. .

[0146] Step 3: Design a backbone network for deep nonlinear modeling to capture the evolution of field data. The backbone network includes a Patch Embedding module, a backbone Transformer network, and a decoder. The high-dimensional features processed by fusion and channel attention mechanisms are input for deep representation of spatial structure and physical constraints, dynamic encoding, and field resolution output. Spatial features are segmented and dimensionality reduced through convolutional Patch Embedding. After adding position encoding, the token sequence is fed into the stacked backbone Transformer network for long-distance dependency modeling and nonlinear dynamic feature learning. Finally, after decoding and outputting the head, the pressure field and saturation field prediction results with the same spatial resolution as the input are obtained.

[0147] Step 3 includes the following sub-steps:

[0148] Step 3.1: [The text appears to be incomplete and contains several grammatical errors. A more accurate translation would require The Patch Embedding module is input, and through convolution operations with a kernel size of 5 and a stride of 5, the 2D spatial features are divided into a sequence of tokens of length N, where each token has a dimension of d, thus mapping from a spatial tensor to a sequence vector. To improve the model's spatial awareness, a learnable positional encoding is added to each token, resulting in a position-encoded sequence. This enables the sequence to have positional distinguishability, forming a token input that enhances spatial information;

[0149] Step 3.2: Sequence The backbone Transformer network is fed into a stacked multi-layered encoding module. Each encoding module includes a multi-head self-attention mechanism and a feedforward neural network. Specifically, the query, key, and value are first obtained through linear transformation and divided into h heads. At each head, global dynamic weighted aggregation is used to achieve long-distance dependent feature interaction, as shown in the formula:

[0150] (11)

[0151] in, These are the query, key, and value, used for matching and weighting in the self-attention mechanism. This is the transpose of the key matrix; For the number of bullish self-attention heads; Feature dimensions for each attention head; Indicates attention output;

[0152] Each encoding module will and input sequence The sequence is obtained through residual connection and LayerNorm normalization. Then the sequence Input to feedforward neural network and sequence The residual connections are then processed by LayerNorm to obtain the output sequence of a single encoding module. This process is repeated for multiple encoding modules, with the output of each module serving as the input of the next module, ultimately resulting in a highly integrated token sequence.

[0153] Step 3.3: The highly integrated token sequence is restored to a low-resolution spatial feature tensor through linear transformation and reshape. The original field size is restored through deconvolution operation with a kernel size of 5 and a stride of 5. After 3×3 convolution and ReLU refinement, the refined feature tensor is output and connected to the pressure branch and saturation branch respectively. The pressure branch uses 1×1 convolution to output the pressure field prediction of the next time step. The saturation branch uses 1×1 convolution followed by a Sigmoid activation function to output the saturation field prediction of the next time step with the value limited to the interval [0,1].

[0154] Step 4: Construct an efficient neural network model for predicting stress and saturation fields based on a dual-branch network structure and a backbone network, such as... Figure 2 As shown, during network training, the initial pressure field, saturation field, and geological parameters are used as model inputs. After neural network inference, the predicted pressure and saturation at the next time step are output. The results are then recycled as new inputs, and this process is repeated to obtain long-term pressure and saturation field predictions. During the model training phase, spatial observation points from some time steps are used to participate in the calculation of data loss. At the same time, the residuals of the physical control equations are introduced as physical loss terms in all time steps and the entire space to ensure physical consistency constraints across the entire spatiotemporal region.

[0155] Step 4, training the model, includes the following sub-steps:

[0156] Step 4.1: Starting from the initial field state, the model sequentially uses the prediction result at each time step as the input for the next time step, continuously rolling and recursively, and finally outputs a complete multi-step physical field dynamic prediction sequence;

[0157] Step 4.2: During recursive prediction, for certain moments within a specified time step range, the real physical field data of the corresponding moment is used to replace the model's prediction result for that moment as the network input for the next time step. For other moments, the model's autoregressive output is used. By gradually reducing the probability of using real physical field data for replacement, the model relies more on its own prediction in the later stages of training, thereby enhancing the stability and generalization ability under long-term recursion.

[0158] Step 4.3: Data Loss Items It is computed only at a specified training time step and a partial set of spatial sampling points. The local spatial sampling points, expanded around the well center at the specified training time step, are selected, and the expression is:

[0159] (12)

[0160] in, They are respectively Time of the first Predicted and actual values ​​for each spatial sampling point The number of spatial sampling points at that time. To calculate the number of time steps for data loss;

[0161] Physical governing equations Constraints are applied across all time steps and the entire space grid to ensure the overall dynamic physical rationality and conservation. The physical governing equations are derived from the oil-water two-phase flow equations. The mathematical equations for the oil-water two-phase flow are as follows:

[0162] (13)

[0163] (14)

[0164] in, For absolute penetration rate, and These are the relative permeabilities of the oil phase and the water phase, respectively. and The viscosity of the oil phase and the viscosity of the water phase are respectively. and These are the formation volume factors for the oil and water phases, respectively. and The pressures of the oil phase and the water phase are respectively. and These are the fluid weights of the oil phase and the water phase, respectively. For spatial height gradient, Porosity and These represent the saturation levels of the oil and water phases, respectively. and These are the formation volume factors for the oil and water phases, respectively. and These represent the volumetric flow rates of the oil and water phases under standard ground conditions, respectively. For point source functions, Indicates the first The coordinates of each well point;

[0165] For the oil-water two-phase flow equation, where permeability is isotropic and gravity and capillary forces are neglected, the volume integral formula for the aqueous phase flow equation is as follows:

[0166] (15)

[0167] in, For the unit volume size, The relative permeability of the aqueous phase;

[0168] Discretization formula (15), after rearranging terms, yields the physical information residual function. As shown below:

[0169] (16)

[0170] in, For the first The volume of each grid, For time step, Use the The grid, the first Each time step , and calculate, For the reservoir model The first grid and its adjacent first grid The water phase conductivity between grids is the product of the flow coefficient and the geometric factor, where the geometric factor is the ratio of the grid contact surface to the center point spacing, and its expression is:

[0171] (17)

[0172] in, For the first The first grid and its adjacent first grid The contact area of ​​each grid, For the first The first grid and its adjacent first grid The center point spacing of each grid Use the The first grid and its adjacent first grid In each grid, the upstream grid properties are arranged according to the pressure difference direction. , , and calculate;

[0173] The physical loss term calculates the mean square error of the physical information residuals across all time steps and the entire spatial grid, and its expression is:

[0174] (18)

[0175] in, for Time of the first Predicted values ​​for each grid point The total number of grid points in a single time step. To predict the number of data points in a dataset;

[0176] Step 4.4: The total loss function is designed as a weighted combination of data item loss and physical equation residual loss. The weight parameters can be dynamically adjusted to balance the model's fitting accuracy and physical reliability. The expression is:

[0177] (19)

[0178] in, Total loss; The weights are for the data loss items.

[0179] Step 5: Iteratively train the network model established in the above steps until the maximum number of iterations is reached to obtain a well-trained, efficient neural network model for predicting pressure and saturation fields. Perform multi-time-step rolling predictions on the input initial pressure field, saturation field, and geological parameters to achieve high-precision prediction and post-processing of pressure and saturation fields throughout the entire time series.

[0180] Step 5 includes the following sub-steps:

[0181] Step 5.1: Use the data from the first 10% of time steps in the dataset from Step 1 to calculate the data loss term. Use the data from the first 20% of time steps to replace the prediction results with the probability set in Step 4.2 as the input for the next time step. Iterate and train the neural network model until the maximum number of iterations is reached. The maximum number of iterations is set to 2000, the initial learning rate is set to 1e-3, the optimizer is the Adam optimizer, and an adaptive learning rate adjustment method is used with a learning rate decrease rate of 0.95.

[0182] Step 5.2: Use the trained new neural network model that integrates the discrete gradient information of multiphase seepage in the reservoir to predict the pressure field and water saturation field of the complete time step. Compare the model prediction results with the real data of numerical simulation and calculate the relative error distribution of the whole time series.

[0183] The generated full-time step pressure field, saturation field, and relative error distribution are organized into a standardized output structure. A comparison chart of the actual and predicted values ​​is plotted in chronological order and the output is saved to intuitively reflect the accuracy of reservoir pressure and saturation prediction.

[0184] To demonstrate the feasibility and superiority of this invention, a 100×100 two-dimensional two-phase heterogeneous reservoir was constructed for verification. The grid size was 10m×10m×5m, the initial reservoir pressure was 150 bar, and the initial water saturation was 0.4. Permeability and porosity were generated using sequential Gaussian simulation, and the permeability and saturation fields are as follows. Figure 4 As shown in (a) and (b), the time step is 10 days, and data loss is calculated using real data from the previous 100 days. A central injection well is set with a fixed injection volume of 50 m³ / s. 3 / d, the bottomhole flowing pressure of the four surrounding production wells is set at 100 bar. Pressure and saturation prediction results are as follows: Figure 5 (a)~(i) and Figure 6 As shown in (a) to (i), the prediction errors at each time step remain at a low level: the average relative error for pressure on day 50 is 0.15%, and the average relative error for water saturation is 0.12%; the average relative errors for pressure and saturation on day 100 are 0.15% and 0.12%, respectively; and by day 500, the average relative errors for both are 0.12% and 0.14%, respectively. The error is only slightly larger at the injection-production well grid, while the error in the well perimeter area is significantly smaller, indicating that the model has good predictive consistency in most areas. Experimental results show that the method of this invention can effectively adapt to the spatial variation of parameters in heterogeneous reservoirs and the complex multiphase flow field coupling problem. It possesses strong generalization stability and achieves efficient data-driven training, with high training efficiency and low data requirements, greatly improving the solution efficiency and prediction accuracy of reservoir seepage problems.

[0185] 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 of reservoir performance prediction integrating discrete gradient information, characterized in that, The method comprises the following steps: Step 1: building a numerical simulation data set of oil-water two-phase flow in a heterogeneous reservoir based on a numerical simulation method; Step 2: designing a double-branch network structure to extract spatial and physical features of input field data in parallel, realizing effective fusion of features, wherein the double-branch network structure comprises a main feature coding branch and a difference operator branch; Step 3: designing a backbone network to perform deep nonlinear modeling and capture the evolution law of field data, wherein the backbone network comprises a Patch Embedding module, a backbone Transformer network and a decoder; Step 4: building an efficient pressure and saturation field prediction neural network model based on the double-branch network structure and the backbone network; in the model training stage, part of the time step spatial region observation points are used to participate in data item loss calculation, and the residual of the physical control equation is introduced as a physical loss item at all time steps and in the whole space; Step 5: obtaining the trained efficient pressure and saturation field prediction neural network model to realize high-precision prediction of the full-time-series pressure field and saturation field.

2. The method of claim 1, wherein, The step 1 comprises the following sub-steps: Step 1.1: establishing a two-dimensional oil-water two-phase reservoir model by using a numerical simulator, and setting different reservoir parameter attribute data, including: grid size, reservoir thickness, initial reservoir pressure, initial water saturation, rock compressibility, oil phase and water phase density, oil phase and water phase volume coefficient, oil phase and water phase viscosity, oil-water relative permeability curve, and setting to not consider capillary force; Step 1.2: using a sequential Gaussian simulation method, combining a Gaussian kernel function and a correlation length parameter to model the covariance of spatial points, generating a log permeability field with spatial heterogeneity, and obtaining a physically reasonable positive permeability field by exponential transformation of the log permeability field; Then, a linear transformation is performed on the permeability field according to an empirical correlation formula to generate a porosity field associated with the permeability distribution, and the expression is: ；(1) wherein, is the porosity; is the permeability; , is an empirical parameter; The permeability field and porosity field data are converted and saved as inc format files according to the numerical simulator specification; Step 1.3: after importing the permeability field and porosity field data into the numerical simulator, setting the positions and well control conditions of the injection wells and production wells, then running the numerical simulator to obtain complete data sequences of the pressure field and saturation field evolving with time, arranging the results into a tensor structure [C, T, H, W], and saving them as mat format data, wherein C represents the number of data channels, T represents the number of time steps, and H and W represent the height and width of the spatial dimensions, respectively.

3. The method of claim 2, wherein the discrete gradient information is fused. The step 2 comprises the following sub-steps: Step 2.1: Read mat data files and inc data files, splice initial pressure and saturation field data, permeability field and porosity field data in the channel dimension to form an input tensor with a structure of [4, 100, 100] ; normalize the data of different channels using the MinMaxScaler normalization method; Step 2.2: using a double-branch network structure to extract spatial and physical features of input field data in parallel, wherein the double-branch network structure comprises a main feature coding branch and a difference operator branch; The main feature encoding branch extracts spatial features of the pressure, saturation and geological parameter field through multi-layer convolution and pooling operation, and reduces the features to a low-dimensional space to strengthen the expression of the main features, and then restores the features to the original input space size through the transpose convolution method, and outputs the spatial feature tensor ; The differential operator branch uses a weight fixed convolution kernel to realize first-order and second-order differential operations on input data, extract physical related features including physical gradients and Laplace operators, effectively encode physical mechanisms, and output a physical constraint feature tensor with the same size as the original input ; Step 2.3: Convert the input tensor The spatial feature tensor output by the main feature encoding branch Physical constraint feature tensor output by the difference operator branch The components are concatenated along the channel dimension to obtain the fused feature tensor. ; then on Applying a 1×1 convolution operation, combined with BatchNorm and ReLU, achieves deep fusion and channel dimensionality reduction of dual-branch features, expressed as: ；(3) wherein, is the fused reduced dimensionality feature tensor, is the weight of the convolution layer; is the bias vector of the convolution layer; The simple channel attention mechanism is introduced, and first, the Global average pooling is performed in the spatial dimension to obtain the global description of each channel That is: ；(4) wherein, is a global description of the th channel; are respectively the height and the width of the th feature map; is the value of the th channel of the fused feature map at the spatial position . Then, the channel weight is learned through two fully connected layers and a Sigmoid activation function, and the expression is: ；(5) wherein, is the adaptive attention weight for the th channel, denotes a Sigmoid activation function, , are the weights of the two fully connected layers, respectively; , are the bias vectors of the two fully connected layers, respectively. Finally, the attention weight acts on each channel of the fusion feature tensor to achieve adaptive feature enhancement, and the expression is: ；(6) wherein, is the value of the final output feature map at the spatial location .​ 4. The method of claim 3, wherein the discrete gradient information is fused. The specific process of the main feature coding branch is as follows: The input tensor obtained in step 2.1 is input into the main feature encoding branch, and spatial feature extraction is prepared. The input tensor obtained in step 2.1 is input into the main feature encoding branch, and spatial feature extraction is prepared. right A first 2D convolution operation is applied, with a kernel size of 3×3, a stride of 1, and padding of 1, followed by non-linear processing using the ReLU activation function. Next, a 2×2 max pooling operation is performed with a stride of 2 to compress the spatial dimension, reducing the data feature space size to half its original size. A second 2D convolution operation is then applied to the pooled feature tensor, with a kernel size of 3×3, a stride of 1, and padding of 1, again using the ReLU activation function for non-linear processing, resulting in a compressed feature tensor. At this point, the data feature space size is H / 2 × W / 2. right A transposed convolution operation is applied to upsample the feature tensor. The transposed convolution kernel size is 2×2, the stride is 2, and the padding is 0. The number of output channels is consistent with the number of channels after the first convolution. The transposed convolution operation restores the feature tensor spatial size from H / 2×W / 2 to the original size H×W. A 3×3 convolution operation is applied to the transposed convolution result with a stride of 1 and padding of 1 to obtain the feature data. The Sigmoid activation function is used to process the feature data, and the output is a spatial feature tensor with the same spatial size as the original input. .

5. The method of predicting reservoir performance by integrating discrete gradient information according to claim 4, wherein, The specific process of the difference operator branch is as follows: The input tensor obtained in step 2.1 is copied and padded around the four sides according to the size of the convolution kernel, expanding the spatial size of the input tensor. An input difference operator branch is used for the extraction of physical constraint features; the input tensor is copied and padded around the four sides according to the size of the convolution kernel, expanding the spatial size of the input tensor. A first-order difference convolution kernel is used to calculate a discrete approximation of the first derivative on the padded input tensor i.e. the spatial first-order gradient feature tensor of the pressure field consistent with the original input size with a kernel size of 3x3, a stride of 1 and padding of 0; For the reservoir model of , the second-order derivative is calculated using a fourth-order precision 9-point difference convolution kernel , that is, the spatial second-order derivative feature tensor of the pressure field consistent with the original input size , the convolution kernel size is 5x5, the step is 1, and the padding is 0, and the expression is: ；(7) ； (8) wherein, and are the x and y grid sizes, respectively, is the grid pressure at coordinates ; and are the x and y grid sizes, respectively. For the oil reservoir model of , the filled pressure field is convolved using a 25-point finite difference convolution kernel, and the approximate spatial second-order derivative , i.e. the higher-precision pressure field spatial second-order derivative characteristic tensor , is calculated by a more precise fourth-order 25-point hybrid grid difference, with a convolution kernel size of 5x5 and a step size of 1; the specific formula is as follows: ；(9) wherein is a weight coefficient, , ; Based on Darcy's seepage equation The coupling relationship, for the output results of the above physical operators, will or With liquid phase flow rate and input tensors The penetration rate in or Perform calculations to obtain the characteristic tensor with strong physical constraints. ;wherein liquid phase mobility The relative permeability tensor is calculated from relative permeability, volume index, and viscosity. It is obtained by interpolation based on the relative permeability curve, with the volume index and viscosity set as constants. The specific formula is as follows: ；(10) wherein, , Krwand Korare the relative permeability of the water and oil phases, respectively; , Kvwand Korare the volume factor of the water and oil phases, respectively; , Kvwand Korare the viscosity of the water and oil phases, respectively; Will , or , and splicing fusion in the channel dimension, generating a physical constraint feature tensor consistent with the input space size .

6. The method of predicting reservoir performance by integrating discrete gradient information according to claim 5, wherein, In the step 3, the following sub-steps are included: Step 3.1: [The text appears to be incomplete and contains several grammatical errors. A more accurate translation would require the full context.] The Patch Embedding module takes a 5-kernel, 5-stride convolution operation as input and divides the 2D spatial features into a sequence of N tokens, each with dimension d. A learnable positional encoding is then added to each token, resulting in a position-encoded sequence. ; Step 3.2: The sequence The main trunk Transformer network of the input stacked multi-layer encoding module, each encoding module includes a multi-head self-attention mechanism and a feedforward neural network; specifically: first, linear transformation is performed to obtain queries, keys and values respectively and divided into h heads, and long-distance dependent feature interaction is realized by globally dynamically weighted aggregation on each head, and the formula is: ；(11) wherein, query, key, value, respectively, for matching and weighting of the self-attention mechanism; is the transpose of the key matrix; is the number of multi-head self-attention heads; is the feature dimension of each attention head; denotes the attention output; Each encoding module will and the input sequence through the residual connection and LayerNorm normalization sequence , and the sequence input into the feedforward neural network and sequence residual connection, and finally through LayerNorm to get the single encoding module output sequence, repeat multiple encoding modules, the output of each module as the input of the next module, and finally get the highly integrated Token sequence; Step 3.3: The highly integrated Token sequence is recovered to a low-resolution spatial feature tensor by linear transformation and reshape, restored to the original field size by a deconvolution operation with a convolution kernel size of 5 and a step size of 5, refined by 3x3 convolution and ReLU, and the output refined feature tensor is connected to the pressure branch and the saturation branch. The pressure branch uses a 1x1 convolution to output the next time pressure field prediction, and the saturation branch uses a 1x1 convolution followed by a Sigmoid activation function to output the next time saturation field prediction with a value limited to the interval [0, 1].

7. The method of predicting reservoir performance by integrating discrete gradient information according to claim 6, wherein, In step 4, the training of the model includes the following sub-steps: Step 4.1: The model starts from the initial field state and feeds back the prediction result of each time step as the input of the next time step, and continuously rolls and recursively predicts, finally outputs the complete multi-step physical field dynamic prediction sequence; Step 4.2: During the recursive prediction, for some time steps within the specified time step range, the real physical field data of the corresponding time step is used to replace the model's prediction result of that time step as the network input of the next time step, and the rest of the time steps use the model's self-predicted output results. By gradually reducing the probability of using real physical field data for replacement, the model relies more on its own prediction in the later training; Step 4.3: Data loss term Only calculated on the specified training time steps and partial spatial sampling points set, the local spatial sampling points are selected on the specified training time steps and expanded around the well center, expressed as: ；(12) wherein, are respectively the prediction value and the true value of the th spatial sampling point at the is the number of spatial sampling points at the is the number of time steps for calculating data loss; Physical control equation terms The constraints at all time steps and the full spatial grid ensure the overall dynamics to be physically reasonable and conservative. The physical control equation terms are derived from the oil-water two-phase seepage equation. First, the oil-water two-phase mathematical equation is as follows: ；(13) ；(14) wherein, is the absolute permeability, and are the relative permeability of the oil and water phases, respectively, and are the viscosities of the oil and water phases, respectively, and are the formation volume factors of the oil and water phases, respectively, and are the pressures of the oil and water phases, respectively, and are the fluid fractions of the oil and water phases, respectively, is the spatial height gradient, is the porosity, and are the saturations of the oil and water phases, respectively, and are the formation volume factors of the oil and water phases, respectively, and are the volumetric flow rates of the oil and water phases at the surface standard conditions, respectively, is the point source function, denotes the well point coordinate position of the th well; For the oil-water two-phase flow equation, the permeability is isotropic, and gravity and capillary force are not considered. After volume integration of the water phase flow equation, the formula is as follows: ；(15) wherein, is the unit volume size, is the water phase relative permeability; Discretizing equation (15), we obtain the physical information residual function as follows: ；(16) wherein, V is the volume of the grid, V is the volume of the grid, is the time step, using the grid, the grid, the grid, , and , is the water phase transmissibility between the grid, the grid, the grid, adjacent to the grid, and is the product of the flow coefficient and the geometric factor, which is the ratio of the grid face to the center point spacing, and is expressed as: ；(17) wherein, the contact area of the first grid with its adjacent second grid, the contact area of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, the center point distance of the first grid with its adjacent second grid, The physical loss term calculates the mean square error of the physical information residual error of all time steps and full spatial grid, and the expression is: ；(18) wherein, is the predicted value of the grid point at time step , is the total number of grid points for a single time step, is the number of data points in the prediction data set; Step 4.4: The total loss function is designed as a weighted combination of data item loss and physical equation residual loss, and the expression is: ；(19) wherein, is the total loss; is the data loss term weight.

Images