Automatic History Matching Method for Large-Scale Reservoirs Based on Stochastic Low-Rank Momentum Acceleration

Abstract

The present invention discloses a large-scale reservoir automatic history matching method based on random low-rank momentum acceleration, belonging to the technical field of reservoir numerical simulation and optimization, and comprising the following steps: initializing reservoir data and partitioning, decomposing the reservoir parameter field into multiple levels of local regions by adopting a macro-block hierarchical strategy, and performing low-rank covariance compression on each macro-block in combination with a Gaussian random projection matrix; designing a time adaptive sampling mechanism, dynamically selecting key time windows based on time sensitivity; calculating the hierarchical random low-rank gain; designing an asynchronous momentum acceleration and adaptive learning rate update strategy; globally synchronously aggregating the results of each macro-block and imposing geophysical constraints; terminating the iteration when the convergence condition is satisfied, and outputting a parameter field that meets the convergence condition. The present invention solves the technical bottlenecks of high computational complexity and large memory consumption existing in traditional methods when dealing with ultra-large-scale reservoir models.

Description

Technical Field

[0001] The present invention belongs to the technical field of reservoir numerical simulation and optimization, and particularly relates to a large-scale reservoir automatic history matching method based on stochastic low-rank momentum acceleration. Background Art

[0002] As the core link of reservoir numerical simulation, automatic history matching aims to achieve the best match between the numerical simulation results and the actual production dynamic data of the oilfield (including pressure, water cut, oil production, etc.) by adjusting geological model parameters (such as permeability, porosity, etc.), so as to improve the prediction accuracy of the reservoir model. This process essentially belongs to the solution of high-dimensional non-linear inverse problems, which is of great value for optimizing the development plan, reducing decision-making risks, and increasing oil recovery.

[0003] In terms of solution methods, gradient optimization methods (such as the Levenberg-Marquardt algorithm, adjoint method) have a fast convergence rate, but they need to calculate high-dimensional Jacobian matrices, with high computational complexity and being prone to falling into local optimal solutions; stochastic search methods (such as Markov chain Monte Carlo, particle swarm optimization) can ensure global convergence, but they require thousands of numerical simulations and are only applicable to extremely small-scale models. Ensemble Kalman Filter (EnKF) updates parameters through ensemble statistics and is naturally suitable for parallel computing, but it needs to perform multiple restart calculations. If the model update is unreasonable, non-physical situations will occur, leading to restart failures. Currently, the Ensemble Smoother with Multiple Data Assimilation (ES-MDA) method is generally adopted as the mainstream technology for efficient history matching. ES-MDA enhances the stability of traditional ensemble methods through multiple data assimilation steps and covariance correction mechanisms, and can better balance computational efficiency and inversion accuracy.

[0004] However, when facing large-scale complex reservoirs with a huge grid scale, a dense number of wells, and a long-term accumulation of production history, the ES-MDA method exposes a series of bottleneck problems. First, the ultra-large grid parameters lead to an exponential growth in the storage and calculation requirements of the covariance matrix, far exceeding the bearing capacity of existing hardware resources. Second, the interaction between massive observation data and high-dimensional parameter spaces makes it difficult to construct and update the gain matrix, and traditional matrix operation methods cannot meet the real-time requirements. In addition, the sharp increase in communication overhead under the parallel computing framework and the cumulative time consumption of data assimilation steps further limit the scalability of the algorithm. These problems significantly reduce the efficiency of traditional methods when dealing with grid models of tens of millions, and it is difficult to meet the urgent need of the oilfield development for real-time history matching. Summary of the Invention

[0005] To address the problems of high computational complexity, large memory consumption, and slow convergence speed in traditional ensemble smoothing multi-data assimilation methods when dealing with ultra-large-scale reservoir models, the present invention proposes a large-scale reservoir automatic history matching method based on stochastic low-rank momentum acceleration. This method decomposes the global problem into local sub-problems through macro-block layering, reducing the dimension of a single covariance matrix; uses stochastic projection to map high-dimensional parameters to a low-dimensional space, significantly reducing the computational amount; dynamically adjusts the projection direction according to parameter sensitivity, preferentially updating high-impact regions; and realizes the efficient utilization of computing resources through a hybrid architecture of task parallelism and data parallelism.

[0006] The technical solution of the present invention is as follows:

[0007] A large-scale reservoir automatic history matching method based on stochastic low-rank momentum acceleration, comprising the following steps:

[0008] Step 1, initialize reservoir data and partitioning, decompose the reservoir parameter field into multiple levels of local regions using the macro-block layering strategy, and perform low-rank covariance compression on each macro-block in combination with a Gaussian random projection matrix;

[0009] Step 2, design a time-adaptive sampling mechanism, dynamically select key time windows based on time sensitivity to optimize the allocation of computing resources;

[0010] Step 3, calculate the hierarchical stochastic low-rank gain;

[0011] Step 4, design an asynchronous momentum acceleration and adaptive learning rate update strategy;

[0012] Step 5, globally synchronously aggregate the results of each macro-block and impose geophysical constraints;

[0013] Step 6, terminate the iteration when the convergence condition is met, and output the parameter field that meets the convergence condition.

[0014] Further, the specific process of step 1 is as follows:

[0015] Step 1.1, load the reservoir parameter field, observation data, and hyperparameters to complete global initialization, and determine the dimension of the parameter field and the scale of the observation data; the hyperparameters include the number of iterations, learning rate, and momentum decay factor;

[0016] Step 1.2, decompose the reservoir parameter field into multiple levels of local regions using the macro-block layering strategy; the specific process of the macro-block layering strategy is: obtain the sedimentary facies map division result of the reservoir, and divide the reservoir parameter field into connected regions according to the sedimentary facies map division result, and one region is one macro-block; obtain the permeability field within each macro-block and calculate the permeability gradient, and perform clustering on each macro-block using the K-means algorithm according to the permeability gradient, and divide each macro-block into Sub - blocks, where one sub - block is a medium - sized block; obtain all grid cells within each medium - sized block, evenly group them according to the grid index, and each group cells, and one cell is a micro - block;

[0017] Step 1.3: Use the Johnson - Lindenstrauss lemma to generate a Gaussian random projection matrix for each macro - block, and combine the Gaussian random projection matrix to perform low - rank covariance compression on each macro - block, mapping the high - dimensional parameters to a low - dimensional space; each element in the matrix satisfies independent and identical distribution: ; where, is the Gaussian random projection matrix of the th macro - block; , are the row index and column index of the Gaussian random projection matrix respectively, identifying the row position and column position of the element in the matrix; is the normal distribution; is the dimension of the low - dimensional space;

[0018] Step 1.4: Initialize the momentum and second - moment estimate of each macro - block to 0, that is , ; where, is the initial momentum of the th macro - block; is the initial second - moment estimate of the th macro - block.

[0019] Furthermore, the specific process of Step 2 is as follows:

[0020] Step 2.1: Calculate the parameter - response covariance matrix at each time point, and the formula is:

[0021] ;

[0022] ;

[0023] ;

[0024] where, is the parameter - response covariance matrix at time point , , is the total number of time points; is the total number of samples; is the index of the sample; is the parameter vector of the th sample; is the mean vector of all sample parameters; is the model response function at time point ; Denote the mean vector of the model response at time point ; is the transpose symbol;

[0025] Step 2.2: Calculate the sensitivity at each time point through the Frobenius norm. The formula is:

[0026] ;

[0027] where, is the sensitivity at time point ; is the Frobenius norm; is 's row index, identifying the element position in the row direction of the matrix; is the total number of macroblocks; is 's column index, identifying the element position in the column direction of the matrix; is 's total dimension in the column direction;

[0028] Step 2.3: Generate the sampling probability of the time window according to the sensitivity calculated in Step 2.2, and randomly select some time periods without replacement from all time points based on the probability distribution;

[0029] The formula for calculating the sampling probability of the time window is:

[0030] ;

[0031] where, is the sampling probability at time point ; is the index of the time point, used to traverse all time points, and the value range is ; is the sensitivity at time point ;

[0032] Then, sample time points without replacement from time points according to the probability of time point . Denote the set of time points obtained by the -th sampling as , where is the index of the sampling times, used to record the time points selected each time.

[0033] Furthermore, the specific process of Step 3 is as follows:

[0034] Step 3.1: Use the Gaussian random projection matrix generated in Step 1.3 to perform low-rank projection on the parameters of each macroblock, compressing the parameters of each macroblock into a low-dimensional space. The formula is:

[0035] ;

[0036] where, is the low-dimensional parameter vector of the th macroblock under the th sample; is the original high-dimensional parameter vector of the th macroblock under the th sample;

[0037] Step 3.2: Calculate the low-dimensional parameter-response covariance and low-dimensional response covariance in the low-dimensional space, and then solve the low-dimensional Kalman gain matrix;

[0038] The formula for calculating the low-dimensional parameter-response covariance is:

[0039] ;

[0040] ;

[0041] ;

[0042] where, is the low-dimensional parameter-response covariance matrix of the th macroblock; is the mean value of the low-dimensional parameters of the th macroblock; is the model response function of the th macroblock under ; is the response mean vector of the th macroblock under all samples;

[0043] The formula for calculating the low-dimensional response covariance is:

[0044] ;

[0045] where, is the low-dimensional response covariance matrix of the th macroblock;

[0046] The formula for calculating the low-dimensional Kalman gain matrix is:

[0047] ;

[0048] where, is the low-dimensional Kalman gain matrix of the th macroblock; is the dilation factor at the th sampling time point; is the observation noise covariance matrix of the th macroblock;

[0049] Step 3.3, reconstruct the low-dimensional Kalman gain matrix to the original parameter space through the inverse mapping of the random projection matrix to obtain the reconstructed Kalman gain matrix. The formula is:

[0050] ;

[0051] where is the reconstructed Kalman gain matrix of the th macroblock.

[0052] Furthermore, the specific process of step 4 is as follows:

[0053] Step 4.1, calculate the momentum through the weighted average of the historical gradient directions. The specific process is as follows:

[0054] For each macroblock, first generate the perturbed observation data:

[0055] ;

[0056] where is the perturbed observation data of the th sampling time point and the th macroblock; is the observation data under ; is the th sampling time point and the th macroblock's random perturbation term, which is used to generate the perturbed observation data, and its distribution is ;

[0057] Then calculate the current parameter residual:

[0058] ;

[0059] where is the parameter residual of the th macroblock; is the th macroblock's parameter vector at the th iteration;

[0060] Finally, update the momentum:

[0061] ;

[0062] where , are respectively the The momentum of the macroblock at the -th and -th iterations; is the momentum decay factor;

[0063] Step 4.2: Dynamically adjust the learning rate based on the second-order moment estimation; the specific process is as follows:

[0064] For each macroblock, first update the second-order moment estimation:

[0065] ;

[0066] where , are the second-order moment estimations of the -th macroblock at the -th and -th iterations respectively; is the second-order moment estimation decay coefficient;

[0067] Then calculate the adaptive learning rate:

[0068] ;

[0069] where is the adaptive learning rate of the -th macroblock at the -th iteration; is the learning rate scaling factor; is the numerical stability constant;

[0070] Step 4.3: Each macroblock updates the parameters asynchronously according to the local reconstruction Kalman gain matrix and the observation residual;

[0071] For each macroblock, update the parameters in parallel:

[0072] ;

[0073] where , are the parameters of the -th macroblock at the -th and -th iterations respectively; and are the reconstruction Kalman gain matrix and the observation residual of the -th macroblock at the -th iteration, , is the dilation factor at the -th iteration.

[0074] Furthermore, the specific process of Step 5 is as follows:

[0075] Step 5.1: Aggregate the updated parameter fields of each macroblock into global parameters. The formula is:

[0076] ;

[0077] where is the global parameter at the th iteration;

[0078] Step 5.2: Impose geophysical constraints on the parameters, including non-negativity constraints and facies-controlled range limits;

[0079] The non-negativity constraint is: ; is the parameter update; is to find the maximum value;

[0080] The facies-controlled range limit is: Constrain the permeability in the channel area within a set range.

[0081] Furthermore, the specific process of Step 6 is as follows:

[0082] Step 6.1: The iteration termination judgment condition is: The norm of the observation residual is less than a preset threshold, that is , or the maximum number of iterations is reached; where is the observation data, which is the actual measured reservoir production dynamic data; is the model response function; is the preset threshold;

[0083] Step 6.2: Save and output the parameter field that meets the convergence condition, and verify its matching degree with the historical production data through a reservoir numerical simulator.

[0084] The beneficial technical effects brought by the present invention: The present invention proposes a parameter update framework, which combines a divide-and-conquer strategy based on macroblock stratification, stochastic low-rank covariance compression, momentum acceleration, and a dynamic sampling mechanism. It is applicable to the efficient inversion of geological parameters such as permeability fields and porosity fields, and solves the technical bottlenecks of high computational complexity and large memory consumption existing in traditional methods when dealing with ultra-large-scale reservoir models. Brief Description of the Drawings

[0085] Figure 1 is the flowchart of the large-scale reservoir automatic history matching method based on stochastic low-rank momentum acceleration of the present invention.

[0086] Figure 2 is the schematic diagram of the reservoir model adopted in the embodiment of the present invention.

[0087] Figure 3 is the fitting diagram of 100 groups of sample oil well and water well production dynamic data in the embodiment of the present invention.

[0088] Figure 4 This is the fitting graph of the production performance data of 40 groups of sample oil and water wells in the embodiments of the present invention.

[0089] Figure 5 This is the variation graph of the objective function of 100 groups of samples in the embodiments of the present invention.

[0090] Figure 6 This is the variation graph of the objective function of 40 groups of samples in the embodiments of the present invention. Detailed implementation manners

[0091] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners:

[0092] As Figure 1 shown, a large-scale reservoir automatic history matching method based on stochastic low-rank momentum acceleration includes the following steps:

[0093] Step 1. Initialize reservoir data and partitioning, decompose the reservoir parameter field into multiple levels of local regions by using a macroblock hierarchical strategy, and combine a Gaussian random projection matrix to perform low-rank covariance compression on each macroblock, significantly reducing the storage and computational complexity of the covariance matrix; the specific process is as follows:

[0094] Step 1.1. Load the reservoir parameter field, observation data, and hyperparameters to complete global initialization, and determine the dimension of the parameter field and the scale of the observation data; the hyperparameters include the number of iterations, learning rate, and momentum decay factor;

[0095] Step 1.2. Decompose the reservoir parameter field into multiple levels of local regions by using a macroblock hierarchical strategy; the specific process of the macroblock hierarchical strategy is as follows: First, divide the reservoir parameter field into multiple macroblocks based on geological features (such as sedimentary facies, fault distribution), and each macroblock covers a continuous geological unit; then, continue to divide within the macroblock into multiple mesoblocks; finally, continue to divide within the mesoblock into multiple microblocks;

[0096] A macroblock (MacroBlock) is divided according to sedimentary facies and reflects the large-scale geological structure. A mesoblock (MesoBlock) is divided within the macroblock according to the permeability gradient and reflects the medium-scale heterogeneity. A microblock (MicroBlock) is a uniform grid unit within the mesoblock that retains local details. Macroblocks, mesoblocks, and microblocks are the three-level structure of the reservoir.

[0097] The process of macroblock division is as follows: Obtain the division result of the sedimentary facies map of the reservoir, and divide the parameter field into connected regions according to the division result of the sedimentary facies map, and one region is one macroblock; define the th macroblock after division as , , , is the total number of macroblocks.

[0098] The process of medium block partitioning is as follows: Obtain the permeability field within each macroblock and calculate the permeability gradient. Then, use the K-means algorithm for clustering based on the permeability gradient to divide each macroblock into sub-blocks, where one sub-block is a medium block. The input of the K-means algorithm is the permeability gradient, and the output is the class label to which each permeability gradient belongs, thus realizing the classification of the regions within the macroblock. The specific steps include randomly initializing the centroids, assigning grid points to the nearest centroid, and updating the centroid positions until convergence. After the K-means algorithm clustering is completed, reorder the cluster labels according to the high and low permeability gradients into high permeability zones, medium permeability zones, and low permeability zones, and finally output the medium block partitioning result of each macroblock. Each medium block contains continuous grid regions, and the grid points in it have similar permeability change trends. Define the permeability field of the th macroblock as , and the calculated permeability gradient of the th macroblock is , is the gradient calculation; divide the th macroblock into sub-blocks, and define the th medium block after partitioning as , . In specific implementation, the high permeability zone is medium block 1, the medium permeability zone is medium block 2, and the low permeability zone is medium block 3.

[0099] The process of microblock partitioning is as follows: Obtain all grid cells within each medium block and directly group them evenly according to the grid index, with cells in each group, and one cell is a microblock; divide the th medium block into cells, and define the th medium block after partitioning as , .

[0100] The multi-level block division (macro-block, medium-block, micro-block) of the reservoir parameter field is mainly a hierarchical management strategy designed to meet different calculation requirements. In the present invention, the macro-block is the smallest unit directly operated by the algorithm and is used to decompose the global high-dimensional parameter field into local sub-problems. The calculation tasks of different macro-blocks can be independently assigned to different computing nodes to achieve distributed processing. The micro-block is the smallest division unit of the parameter field, corresponding to the grid of the reservoir, and the parameter field update actually acts at the micro-block level. The medium-block is an intermediate level between the macro-block and the micro-block. When the calculation tasks within a macro-block are unbalanced (for example, some macro-block parameters are highly sensitive and computationally intensive), the macro-block can be further split into medium-blocks to optimize task allocation. In addition, in an environment with limited memory, the medium-block can be used as the smallest unit for data loading. In this case, the medium-block is equivalent to the macro-block. For the sake of simplicity in description, in the subsequent steps, the smallest unit directly operated will be referred to as the macro-block.

[0101] Step 1.3: Use the Johnson–Lindenstrauss lemma (a high-dimensional data dimensionality reduction theory, abbreviated as the JL lemma) to generate a Gaussian random projection matrix for each macro-block, and combine the Gaussian random projection matrix to perform low-rank covariance compression on each macro-block, mapping the high-dimensional parameters to a low-dimensional space; each element in the matrix satisfies independent and identical distribution: ; where is the Gaussian random projection matrix of the -th macro-block; and are the row index and column index of the Gaussian random projection matrix respectively, identifying the row position and column position of the element in the matrix; is the normal distribution, describing the probability distribution characteristics of the random variable; is the dimension of the low-dimensional space, that is, the size of the low-dimensional space to which the high-dimensional parameters are mapped through random projection;

[0102] Step 1.4: Initialize the momentum and second-order moment of each macro-block simultaneously for subsequent adaptive parameter updates. Initialize the momentum and second-order moment estimates of each macro-block to 0, that is , ; where is the initial momentum of the -th macro-block; is the initial second-order moment estimate of the -th macro-block.

[0103] Step 2: Design a time-adaptive sampling mechanism to dynamically select key time windows based on time sensitivity and optimize the allocation of computing resources; the specific process is as follows:

[0104] Step 2.1: Calculate the parameter-response covariance matrix at each time point, and the formula is:

[0105] ;

[0106] ;

[0107] ;

[0108] wherein, is the parameter-response covariance matrix at time point , , is the total number of time points; is the total number of samples; is the index of the sample, and the value range is ; is the parameter vector of the th sample, describing the parameter configuration of the model; is the mean vector of all sample parameters; is the model response function at time point , which maps the parameters to the corresponding response values; represents the mean vector of the model responses at time point ; is the transpose symbol;

[0109] Step 2.2. Calculate the sensitivity at each time point through the Frobenius norm. The formula is:

[0110] ;

[0111] wherein, is the sensitivity at time point ; is the Frobenius norm; is 's row index, identifying the element position in the row direction of the matrix; total number of macro blocks; is 's column index, identifying the element position in the column direction of the matrix; is 's total number of dimensions in the column direction;

[0112] Step 2.3. Generate the sampling probability of the time window according to the sensitivity calculated in Step 2.2, and randomly select some time periods without replacement from all time points based on the probability distribution, giving priority to focusing on the observation data at highly sensitive moments.

[0113] The formula for calculating the sampling probability of the time window is:

[0114] ;

[0115] wherein, is the sampling probability at time point ; is the index of the time point, used to traverse all time points, and the value range is ; is the sensitivity at time point ;

[0116] Then, sample without replacement time points from time points according to the probability of time point . Denote the set of time points obtained in the th sampling as . Among them, is the index of the sampling times, used to record the time points selected in each sampling.

[0117] Step 3: Calculate the hierarchical random low-rank gain; the specific process is as follows:

[0118] Step 3.1: Use the Gaussian random projection matrix generated in Step 1.3 to perform low-rank projection on the parameters of each macroblock, and compress the parameters of each macroblock into a low-dimensional space. The formula is:

[0119] ;

[0120] where is the low-dimensional parameter vector of the th macroblock under the th sample, which is the representation after compressing the original high-dimensional parameters into a low-dimensional space through the Gaussian random projection matrix; is the original high-dimensional parameter vector of the th macroblock under the th sample, used to describe the geological parameter configuration of this macroblock;

[0121] Step 3.2: Calculate the low-dimensional parameter-response covariance and low-dimensional response covariance in the low-dimensional space, and then solve the low-dimensional Kalman gain matrix.

[0122] The calculation formula for the low-dimensional parameter-response covariance is:

[0123] ;

[0124] ;

[0125] ;

[0126] where is the low-dimensional parameter-response covariance matrix of the th macroblock; is the The mean value of the low-dimensional parameters of a macroblock; At the model response function of the th macroblock maps the low-dimensional parameters to the corresponding response values; is the response mean vector of the

[0127] The calculation formula for the low-dimensional response covariance is:

[0128] ;

[0129] where, is the low-dimensional response covariance matrix of the th macroblock;

[0130] The calculation formula for the low-dimensional Kalman gain matrix is:

[0131] ;

[0132] where, is the low-dimensional Kalman gain matrix of the th macroblock; is the dilation factor at the th subsampling time point; is the observation noise covariance matrix of the th macroblock;

[0133] Step 3.3, reconstruct the low-dimensional Kalman gain matrix to the original parameter space through the inverse mapping of the random projection matrix to obtain the reconstructed Kalman gain matrix, and the formula is:

[0134] ;

[0135] where, is the reconstructed Kalman gain matrix of the th macroblock.

[0136] Step 4, design an asynchronous momentum acceleration and adaptive learning rate update strategy, calculate the momentum and second-order moment estimation through the weighted average of the historical gradient directions, and improve the parameter update stability and convergence speed; the specific process is:

[0137] Step 4.1, calculate the momentum through the weighted average of the historical gradient directions to reduce the oscillation of parameter updates; the specific process is:

[0138] For each macroblock, first generate perturbed observation data:

[0139] ;

[0140] where, is the perturbed observation data at the -th sampling time point and for the -th macroblock; is the observation data under ; is the random perturbation term at the -th sampling time point and for the -th macroblock, used to generate the perturbed observation data, and its distribution is ;

[0141] Then calculate the current parameter residual:

[0142] ;

[0143] where is the parameter residual for the -th macroblock; is the parameter vector for the -th macroblock at the -th iteration;

[0144] Finally, update the momentum:

[0145] ;

[0146] where , are the momentum at the -th macroblock at the -th and -th iterations respectively; is the momentum decay factor, with a value range in [0, 1), used to control the weight of the historical momentum in the current momentum calculation;

[0147] Step 4.2, Dynamically adjust the learning rate based on the second moment estimation, and use a larger step size for the highly sensitive region; the specific process is as follows:

[0148] For each macroblock, first update the second moment estimation:

[0149] ;

[0150] where , are the second moment estimations at the -th macroblock at the -th and -th iterations respectively; is the second moment estimation decay coefficient, with a value range in [0, 1), used to control the weight of the historical second moment estimation in the current estimation;

[0151] Then calculate the adaptive learning rate:

[0152] ;

[0153] Among them, is the adaptive learning rate at the th iteration of the th macroblock; is the learning rate scaling factor, which overall controls the magnitude of the learning rate; is the numerical stability constant, which is used to avoid the case of division by zero and ensure the stability of the calculation;

[0154] Step 4.3: Each macroblock updates the parameters asynchronously according to the local reconstruction Kalman gain matrix and the observation residual, without waiting for global synchronization.

[0155] For each macroblock, update the parameters in parallel:

[0156] ;

[0157] Among them, , are respectively the parameters at the th iteration of the th and th iterations of the th macroblock; and are respectively the reconstruction Kalman gain matrix and the observation residual at the th iteration of the , is the dilation factor at the th iteration.

[0158] Step 5: Globally synchronize and aggregate the results of each macroblock and impose geophysical constraints to ensure the rationality and spatial continuity of the inversion parameters; the specific process is as follows:

[0159] Step 5.1: Aggregate the updated parameter fields of each macroblock into global parameters, and eliminate the boundary differences in the overlapping areas through weighted averaging to ensure spatial continuity. The formula is:

[0160] ;

[0161] Among them, is the global parameter at the th iteration;

[0162] Step 5.2: Introduce geological rules such as non-negativity constraints (such as non-negative permeability) and phase control range limitations to correct possible non-physical solutions. Impose geophysical constraints on the parameters, including non-negativity constraints and phase control range limitations, which are respectively:

[0163] The non-negativity constraint is: (such as non-negative permeability). For parameter update; For finding the maximum value;

[0164] The phased range is limited to: restricting the permeability of the channel area within a set range.

[0165] Step 6: When the convergence condition is met, the iteration terminates, and the parameter field that meets the convergence condition is output; the specific process is as follows:

[0166] Step 6.1: Monitor the decline rate of the objective function (observation residual norm). When the preset threshold or the maximum number of iterations is reached, the convergence condition is met, and the calculation terminates.

[0167] The iteration termination judgment condition is: the observation residual norm is less than the preset threshold, that is , or the maximum number of iterations is reached. Among them, is the observed data, which is the actual measured reservoir production dynamic data; is the model response function, which maps the parameters to the corresponding response values; is the preset threshold;

[0168] Step 6.2: Save and output the parameter field that meets the convergence condition, and verify its matching degree with the historical production data through a reservoir numerical simulator.

[0169] In order to prove the feasibility and superiority of the present invention, the following specific embodiments are given.

[0170] The embodiment of the present invention adopts a classic channel-type reservoir model (EGG model), and the model structure is as Figure 2 shown. The model grid system is 60×60×1, the total number of grids is 3600, and the heterogeneous channel sand body distribution is simulated. The reservoir contains 4 production wells (numbered production well 1 to production well 4 in sequence) and 8 injection wells (numbered injection well 1 to injection well 8 in sequence). The production period is 750 days, divided into 25 time points (30 days per step). The target inversion parameter is the permeability field, and an initial prior model set (100 realizations) is generated through geological modeling. The observed data includes the oil production rate (WOPR) and water production rate (WWPR) of each production well, with a total of 200 observation points (2 observation values for each production well at each time point, 4 wells × 25 time points × 2 indicators = 200). The channel area (high permeability zone) is distributed in a curved strip shape, with a permeability range of 500 - 3500 mD, and the non-channel area (low permeability zone) has a permeability of 1 - 100 mD.

[0171] Initialize the macro-block layering according to the process of Step 1: First, load the parameter field, observation data, and hyperparameters, set the number of iterations to 4, the learning rate to 1, and the momentum decay factor to 0.9. Then, divide the reservoir into 12 macro-blocks according to the formation sedimentary facies characteristics (channel, levee, floodplain), with each macro-block covering a continuous geological unit. Perform K-means algorithm clustering within the macro-block according to the permeability gradient, further subdivide the macro-block into medium-blocks, and finally form a grid parameter field at the micro-block level. Finally, generate a Gaussian random projection matrix for each macro-block to compress the original parameter space to a low-dimensional space.

[0172] Calculate the time sensitivity and sampling according to the process of Step 2: First, calculate the parameter-response covariance matrix at each time point, and evaluate the sensitivity at each time point through the Frobenius norm. Then, generate a probability distribution based on the sensitivity, and sample 5 highly sensitive time windows (each group contains 5 consecutive time points) without replacement from 25 time points, covering the stage with drastic changes in oil / water production rate.

[0173] Calculate the hierarchical low-rank gain according to the process of Step 3: First, perform low-rank projection on the parameters of each macro-block; then, calculate the low-dimensional response covariance and low-dimensional Kalman gain matrix in the low-dimensional space; finally, restore the low-dimensional Kalman gain matrix to the original space through inverse mapping.

[0174] Perform asynchronous momentum acceleration update according to the process of Step 4: For each macro-block, first, update the momentum term by weighted average of the historical gradient direction; then, update the second moment, and update the parameters in combination with the adaptive learning rate, where the adaptive learning rate is dynamically adjusted based on the second moment estimation, and the step size in the highly sensitive area is increased by 30%.

[0175] Perform global synchronization and constraints according to the process of Step 5: First, perform parameter field aggregation, that is, merge the update results of each macro-block, eliminate the boundary differences through weighted average, and ensure the spatial continuity of the permeability field. Then, perform physical constraints, that is, impose non-negativity constraints (for example, the permeability in the channel area ≥ 0) and phase control range (10 mD ≤ the permeability in the channel area ≤ 500 mD) to correct non-physical solutions.

[0176] Perform iteration termination and result verification according to the process of Step 6: First, perform convergence judgment, and reach the preset threshold (1e-4) after 4 iterations by calculating the decline rate of the monitoring objective function (observation residual norm). Then, perform simulation verification through a numerical simulator, and the matching error between the inversion model and the historical production data (oil production, water production) is less than 5%.

[0177] Select 100 groups of samples and 40 groups of production dynamic data of oil wells and water wells respectively for production data fitting experiments, and the fitting results are Figure 3 and Figure 4as shown Figure 3 and Figure 4 In Figure 4 , the initial prior model simulation results in a fitted production dynamic data curve (dark gray) obtained after 4 iterations of the production dynamic data curve (light gray), which tightly wraps the actual observed production dynamic data curve (black). For cases with fewer samples Figure 4 compared to those with more samples Figure 3 , the fitted production data after iteration matches the actual observed data well, with a smaller fitting error.

[0178] The convergence of the objective function for 100 - sample and 40 - sample cases is shown respectively in Figure 5 and Figure 6 as shown Figure 5 and Figure 6 show that the historical fitted objective function value continuously decreases with the increase of the number of iterations and finally converges. For cases with fewer samples Figure 6 compared to those with more samples Figure 5 , the objective function is easier to converge and the value is smaller.

[0179] The above results verify the efficiency and stability of the method of the present invention in large - scale reservoir history matching. By macro - block layering and low - rank covariance compression, the computational complexity is significantly reduced, and time - adaptive sampling and momentum acceleration further improve the convergence speed, providing a feasible solution for real - time fitting of ten - million - grid models.

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

Claims

1. A large-scale reservoir automatic history matching method based on random low-rank momentum acceleration, characterized in that: The steps include: Step 1: Initialize reservoir data and blocks, use macroblock hierarchical strategy to decompose reservoir parameter field into multi-level local regions, and perform low-rank covariance compression on each macroblock in combination with Gaussian random projection matrix; Step 2: Design a time-adaptive sampling mechanism to dynamically select key time windows based on time sensitivity to optimize computing resource allocation; Step 3, calculate the hierarchical random low-rank gain; Step 4: Design asynchronous momentum acceleration and adaptive learning rate update strategy; Step 5: Globally and synchronously aggregate the results of each macroblock and apply geophysical constraints; Step 6: When the convergence condition is met, the iteration is terminated, and the parameter field that meets the convergence condition is output.

2. The large-scale reservoir automatic history matching method based on random low-rank momentum acceleration according to claim 1 is characterized in that: The specific process of step 1 is as follows: Step 1.1, load the reservoir parameter field, observation data and hyperparameters to complete global initialization, determine the dimension of the parameter field and the scale of observation data; hyperparameters include the number of iterations, learning rate, and momentum decay factor; Step 1.2: Use the macroblock stratification strategy to decompose the reservoir parameter field into multiple local regions. The specific process of the macroblock stratification strategy is as follows: obtain the sedimentary phase map division results of the reservoir, and divide the reservoir parameter field into connected areas, one area is a macroblock; obtain the permeability field in each macroblock and calculate the permeability gradient, cluster using the K-means algorithm according to the permeability gradient, and divide each macroblock into sub-blocks, one sub-block is a middle block; get all grid cells in each middle block, group them evenly by grid index, and Units, one unit is a microblock; Step 1.3: Use the Johnson-Lindenstrauss lemma to generate a Gaussian random projection matrix for each macroblock, and combine the Gaussian random projection matrix to perform low-rank covariance compression on each macroblock to map high-dimensional parameters to low-dimensional space; each element in the matrix satisfies independent and identical distribution: ;in, For the Gaussian random projection matrix of macroblocks; , They are Gaussian random projection matrices The row index and column index of the matrix identify the row and column positions of the elements in the matrix; is a normal distribution; is the dimension of the low-dimensional space; Step 1.4: Initialize the momentum and second-order moment estimates of each macroblock to 0, that is, , ;in, For the The initial momentum of the macroblock; For the The initial second-order moment estimate of the macroblock.

3. The large-scale reservoir automatic history matching method based on random low-rank momentum acceleration according to claim 2 is characterized in that: The specific process of step 2 is: Step 2.1, calculate the parameter-response covariance matrix at each time point, the formula is: ； ； ； in, For time point The parameter-response covariance matrix, , is the total number of time points; is the total number of samples; is the index of the sample; For the The parameter vector of samples; is the mean vector of all sample parameters; For at time point The model response function of Indicates at a point in time The mean vector of model responses of ; is the transpose symbol; Step 2.2: Calculate the sensitivity at each time point using the Frobenius norm. The formula is: ； in, For time point sensitivity; is the Frobenius norm; for The row index of , identifies the element position in the row direction of the matrix; is the total number of macroblocks; for The column index of , which identifies the element position in the column direction of the matrix; for The total number of dimensions in the column direction; Step 2.3, generate the sampling probability of the time window according to the sensitivity calculated in step 2.2, and extract some time periods from all time points without replacement based on the probability distribution; The sampling probability calculation formula of the time window is: ； in, For time point The sampling probability of is the index of the time point, used to traverse all time points, and the value range is ; For time point sensitivity; Then, by time point Probability from time points without replacement Time point, remember The time point set obtained by sampling is ,in It is the index of the sampling times, which is used to record the time point selected for each sampling.

4. The large-scale reservoir automatic history matching method based on random low-rank momentum acceleration according to claim 3 is characterized in that: The specific process of step 3 is as follows: Step 3.1: Use the Gaussian random projection matrix generated in step 1.3 to perform low-rank projection on the parameters of each macroblock, and compress the parameters of each macroblock into a low-dimensional space. The formula is: ； in, For the The macroblock is Low-dimensional parameter vector for samples; For the The macroblock is The original high-dimensional parameter vector under samples; Step 3.2, calculating the low-dimensional parameter-response covariance and the low-dimensional response covariance in the low-dimensional space, and then solving the low-dimensional Kalman gain matrix; The formula for calculating the low-dimensional parameter-response covariance is: ； ； ； in, For the The low-dimensional parameters of the macroblocks - the response covariance matrix; For the The mean of the low-dimensional parameters of macroblocks; For Next, The model response function of each macroblock; For the The mean response vector of a macroblock under all samples; The formula for calculating the low-dimensional response covariance is: ； in, For the The low-dimensional response covariance matrix of macroblocks; The calculation formula of the low-dimensional Kalman gain matrix is: ； in, For the The low-dimensional Kalman gain matrix of macroblocks; For the dilation factor at the subsampling time point; For the The observed noise covariance matrix of the macroblocks; Step 3.3, reconstruct the low-dimensional Kalman gain matrix to the original parameter space through the inverse mapping of the random projection matrix to obtain the reconstructed Kalman gain matrix, the formula is: ； in, For the The reconstructed Kalman gain matrix of the macroblock.

5. The large-scale reservoir automatic history matching method based on random low-rank momentum acceleration according to claim 4 is characterized in that: The specific process of step 4 is as follows: Step 4.1: Calculate momentum by taking the weighted average of the historical gradient directions. The specific process is: For each macroblock, first generate the perturbed observation data: ； in, For the Sampling time point, The perturbed observation data of macroblocks; For The following observation data; For the Sampling time point, The random perturbation term of macroblocks is used to generate perturbation observation data, and its distribution is ; Then calculate the current parameter residual: ； in, For the The parameter residual of each macroblock; For the The macroblock is The parameter vector at the iteration; Last updated momentum: ； in, , Respectively Macroblock sequence Momentum at iteration ; is the momentum decay factor; Step 4.2: Dynamically adjust the learning rate based on the second-order moment estimation; the specific process is: For each macroblock, first update the second-order moment estimate: ； in, , Respectively Macroblock sequence The second-order moment estimate at the iteration; Estimate the decay coefficient for the second-order moment; Then calculate the adaptive learning rate: ； in, For the Macroblock Adaptive learning rate at iterations; is the learning rate scaling factor; is a numerical stability constant; Step 4.3, each macroblock updates parameters asynchronously based on the local reconstructed Kalman gain matrix and the observed residual; For each macroblock, update the parameters in parallel: ； in, , Respectively Macroblock sequence The parameters at the iteration; and Respectively Macroblock The reconstructed Kalman gain matrix and observation residual at the iteration, , For the The dilation factor for the iteration.

6. The large-scale reservoir automatic history matching method based on random low-rank momentum acceleration according to claim 5 is characterized in that: The specific process of step 5 is as follows: Step 5.1: Aggregate the updated parameter fields of each macroblock into global parameters. The formula is: ； in, For the Global parameters at iterations; Step 5.2, applying geophysical constraints to the parameters, including non-negativity constraints and phase control range restrictions; The non-negativity constraints are: ; is parameter update; To find the maximum value; The phase control range is limited to: constraining the permeability of the river channel area to be within the set range.

7. The large-scale reservoir automatic history matching method based on random low-rank momentum acceleration according to claim 6 is characterized in that: The specific process of step 6 is as follows: Step 6.1, the iteration termination judgment condition is: the observed residual norm is less than the preset threshold, that is , or the maximum number of iterations is reached; where, Observation data are the reservoir production dynamic data actually measured; is the model response function; is the preset threshold; Step 6.2: Save and output the parameter field that meets the convergence conditions, and verify its matching degree with the historical production data through the reservoir numerical simulator.

Images