Large-scale oil reservoir automatic history fitting method based on random low-rank momentum acceleration

Abstract

The invention discloses a large-scale oil reservoir automatic history fitting method based on random low-rank momentum acceleration, and belongs to the technical field of oil reservoir numerical simulation and optimizing.The large-scale oil reservoir automatic history fitting method comprises the following steps that oil reservoir data and blocks are initialized, an oil reservoir parameter field is decomposed into multi-level local areas through a macro block layering strategy, and the multi-level local areas are divided into multiple levels; performing low-rank covariance compression on each macro block in combination with a Gaussian random projection matrix; designing a time adaptive sampling mechanism, and dynamically selecting a key time window based on time sensitivity; calculating a layered random low-rank gain; designing an asynchronous momentum acceleration and adaptive learning rate updating strategy; globally and synchronously aggregating all macro block results and applying geological physical constraints; when the convergence condition is met, iteration is terminated, and the parameter field meeting the convergence condition is output. According to the method, the technical bottlenecks of high calculation complexity and large memory consumption when a traditional method is used for processing a super-large-scale oil reservoir model are solved.

Description

Technical Field

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

[0002] As the core link of reservoir numerical simulation, automatic history matching is to adjust geological model parameters (such as permeability, porosity, etc.) to achieve the best match between numerical simulation results and actual oilfield production dynamic data (including pressure, water content, oil production, etc.), thereby improving the prediction accuracy of reservoir models. This process is essentially a high-dimensional nonlinear inverse problem solution, which is of great value in optimizing development plans, reducing decision-making risks, and improving recovery.

[0003] At the level of solution methods, although gradient optimization methods (such as Levenberg-Marquardt algorithm and adjoint method) have fast convergence speed, they need to calculate high-dimensional Jacobian matrices, have high computational complexity, and are prone to fall into local optimal solutions; although random search methods (such as Markov chain Monte Carlo and particle swarm optimization) can guarantee global convergence, they require thousands of numerical simulations and are only applicable to very small-scale models. Ensemble Kalman filtering (EnKF) updates parameters through ensemble statistics and is naturally suitable for parallel computing, but it requires multiple restart calculations. If the model update is unreasonable, non-physical situations will occur, which will lead to restart failure. At present, the ensemble smoothing multiple data assimilation (ES-MDA) method is widely used as the mainstream technology for efficient history fitting. ES-MDA enhances the stability of traditional ensemble methods through multiple data assimilation steps and covariance correction mechanism, and can better balance computational efficiency and inversion accuracy.

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

[0005] In order to solve the problems of high computational complexity, large memory consumption and slow convergence speed existing in traditional ensemble smoothing multi-data assimilation methods when processing ultra-large-scale reservoir models, the present invention proposes a large-scale reservoir automatic history fitting method based on random low-rank momentum acceleration. This method decomposes the global problem into local sub-problems through macroblock stratification to reduce the dimension of a single covariance matrix; uses random projection to map high-dimensional parameters to low-dimensional space, significantly reducing the amount of calculation; dynamically adjusts the projection direction according to parameter sensitivity, and gives priority to updating high-impact areas; and achieves 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: A large-scale reservoir automatic history matching method based on random low-rank momentum acceleration includes the following steps: 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.

[0007] Furthermore, 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.

[0008] Furthermore, the specific process of step 2 is as follows: 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.

[0009] Furthermore, 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 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.

[0010] Furthermore, 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 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 the 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 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.

[0011] Furthermore, 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.

[0012] Furthermore, 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.

[0013] Beneficial technical effects brought by the present invention: The present invention proposes a parameter update framework that combines a macroblock-based divide-and-conquer strategy, random low-rank covariance compression, momentum acceleration, and a dynamic sampling mechanism. It is suitable for efficient inversion of geological parameters such as permeability fields and porosity fields, and solves the technical bottlenecks of high computational complexity and high memory consumption in traditional methods when processing ultra-large-scale reservoir models. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 The present invention is a flow chart of the large-scale reservoir automatic history matching method based on random low-rank momentum acceleration.

[0015] Figure 2 Schematic diagram of the reservoir model used in the embodiment of the present invention.

[0016] Figure 3 It is a fitting diagram of production dynamic data of 100 groups of sample oil and water wells in an embodiment of the present invention.

[0017] Figure 4 This is a fitting diagram of production dynamic data of 40 groups of sample oil and water wells in an embodiment of the present invention.

[0018] Figure 5 Graph showing changes in the objective function of 100 groups of samples in an embodiment of the present invention.

[0019] Figure 6 This is a graph showing changes in the objective function of 40 groups of samples in an embodiment of the present invention. DETAILED DESCRIPTION

[0020] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments: like Figure 1 As shown, a large-scale reservoir automatic history matching method based on random low-rank momentum acceleration includes the following steps: Step 1: Initialize reservoir data and block, use macroblock layering strategy to decompose reservoir parameter field into multi-level local areas, and combine Gaussian random projection matrix to perform low-rank covariance compression on each macroblock, which significantly reduces the storage and calculation complexity of covariance matrix; the specific process 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, using the macroblock stratification strategy to decompose the reservoir parameter field into multiple levels of local areas; the specific process of the macroblock stratification strategy is: first, based on geological characteristics (such as sedimentary phase, fault distribution), the reservoir parameter field is divided into multiple macroblocks, each macroblock covers a continuous geological unit; then, the macroblock is further divided into multiple medium blocks; finally, the medium block is further divided into multiple microblocks; MacroBlock refers to the division according to sedimentary phase, reflecting the large-scale geological structure. MesoBlock is divided according to the permeability gradient within the macroblock, reflecting the mesoscale heterogeneity. MicroBlock is a uniform grid unit within the mesoblock, retaining local details. Macroblock, mesoblock and microblock are the three-level structures of the reservoir.

[0021] The macroblock division process is as follows: obtain the sedimentary phase map division results of the reservoir, and divide the parameter field into connected regions, one region is a macroblock; define the first The macroblock is , , , is the total number of macroblocks.

[0022] The process of block division is as follows: obtain the permeability field in each macroblock and calculate the permeability gradient, cluster the permeability gradient using the K-means algorithm, and divide each macroblock into sub-blocks, one sub-block is a middle block; the input of the K-means algorithm is the permeability gradient, and the output is the category label to which each permeability gradient belongs, thereby realizing the classification of the area within the macroblock. The specific steps include randomly initializing the centroid, assigning the grid points to the nearest centroid, and updating the centroid position until convergence. After the K-means algorithm clustering is completed, the cluster labels are reordered into high permeability zones, medium permeability zones, and low permeability zones according to the permeability gradient. Finally, the middle block division result of each macroblock is output. Each middle block contains a continuous grid area, and its grid points have similar permeability change trends. Define the first The permeability field of a macroblock is , the calculated The permeability gradient of a macroblock is , is the gradient calculation; The macroblocks are divided into sub-blocks, defining the The middle block is , In a specific implementation, the high permeability zone is middle block 1, the medium permeability zone is middle block 2, and the low permeability zone is middle block 3.

[0023] The process of micro-block division is as follows: obtain all grid cells in each medium block, and evenly group them by grid index. units, one unit is a microblock; The blocks are divided into units, defining the first The middle block is , .

[0024] The multi-level block division (macroblock, medium block, microblock) of the reservoir parameter field is mainly a hierarchical management strategy designed to adapt to different computing needs. In the present invention, the macroblock is the smallest unit directly operated by the algorithm, which is used to decompose the global high-dimensional parameter field into local sub-problems. The computing tasks of different macroblocks can be independently assigned to different computing nodes to achieve distributed processing. The microblock is the smallest division unit of the parameter field, corresponding to the grid of the reservoir, and the parameter field update actually acts on the microblock level. The medium block is an intermediate level between the macroblock and the microblock. When the computing tasks in the macroblock are unbalanced (such as some macroblocks with high parameter sensitivity and large amount of calculation), the macroblock can be further divided into medium blocks to optimize task allocation. In addition, in a memory-constrained environment, the medium block can be used as the smallest unit for data loading. In this case, the medium block is equivalent to a macroblock. To simplify the description, the smallest unit directly operated in the subsequent steps is referred to as a macroblock.

[0025] Step 1.3: Use Johnson–Lindenstrauss Lemma (JL Lemma, a high-dimensional data dimensionality reduction theory) 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; It is a normal distribution, which describes the probability distribution characteristics of random variables; 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 by random projection; Step 1.4: Initialize the momentum and second-order moment of each macroblock at the same time for subsequent adaptive parameter update. 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.

[0026] Step 2: Design a time-adaptive sampling mechanism to dynamically select key time windows based on time sensitivity and optimize computing resource allocation. The specific process is as follows: 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, and its value range is ; For the The parameter vector of samples describes the parameter configuration of the model; is the mean vector of all sample parameters; For at time point The model response function maps the parameters to the corresponding response values; 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; 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, giving priority to the observation data at highly sensitive moments.

[0027] 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.

[0028] Step 3: Calculate the hierarchical random low-rank gain; the specific process is: 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 The low-dimensional parameter vector under samples is the representation after the original high-dimensional parameters are compressed into a low-dimensional space through the Gaussian random projection matrix; For the The macroblock is The original high-dimensional parameter vector under samples is used to describe the geological parameter configuration of the macroblock; 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 maps low-dimensional parameters to corresponding response values; 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 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.

[0029] Step 4: Design asynchronous momentum acceleration and adaptive learning rate update strategy, calculate momentum and second-order moment estimation by weighted average of historical gradient directions, and improve parameter update stability and convergence speed; the specific process is as follows: Step 4.1: Calculate momentum by weighted average of historical gradient directions to reduce oscillation of parameter updates. 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, with a value range of [0,1), which is used to control the weight of historical momentum in the current momentum calculation; Step 4.2: Dynamically adjust the learning rate based on the second-order moment estimation, and use a larger step size for highly sensitive areas; the specific process is: For each macroblock, first update the second-order moment estimate: ; in, , Respectively Macroblock sequence The second moment estimate at the iteration; is the second-order moment estimation attenuation coefficient, with a value range of [0,1), which is used to control the weight of the historical second-order moment estimation in the current estimation; Then calculate the adaptive learning rate: ; in, For the Macroblock Adaptive learning rate at iterations; It is the learning rate scaling factor, which controls the overall learning rate; is a numerical stability constant used to avoid division by zero and ensure the stability of calculation; Step 4.3: Each macroblock updates parameters asynchronously based on the local reconstructed Kalman gain matrix and the observed residuals without waiting for global synchronization.

[0030] For each macroblock, update the parameters in parallel: ; in, , Respectively Macroblock sequence 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.

[0031] Step 5: Globally and synchronously 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: Step 5.1: Aggregate the updated parameter fields of each macroblock into global parameters, and eliminate the boundary differences of the overlapping areas through weighted averaging to ensure spatial continuity. The formula is: ; in, For the Global parameters at iterations; Step 5.2: Introduce geological rules such as non-negative constraints (such as non-negative permeability) and phase control range restrictions to correct possible non-physical solutions. Apply geophysical constraints to the parameters, including non-negative constraints and phase control range restrictions, which are: The non-negativity constraints are: (If the permeability is non-negative). 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.

[0032] Step 6: When the convergence condition is met, the iteration is terminated and the parameter field that meets the convergence condition is output; the specific process is as follows: Step 6.1: Monitor the rate of decrease 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 is terminated.

[0033] The condition for iterative termination is: the observed residual norm is less than the preset threshold, that is, , or the maximum number of iterations is reached. Observation data are the reservoir production dynamic data actually measured; is the model response function, mapping parameters to corresponding response values; 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.

[0034] In order to demonstrate the feasibility and superiority of the present invention, the following specific examples are given.

[0035] The embodiment of the present invention adopts the classic channel type reservoir model (EGG model), and the model structure is as follows: Figure 2As shown in the figure, the model grid system is 60×60×1, with a total grid number of 3600, simulating the distribution of heterogeneous channel sand bodies. The reservoir contains 4 production wells (numbered as production well 1 to production well 4) and 8 water injection wells (numbered as water injection well 1 to water injection well 8), with a production cycle of 750 days, divided into 25 time points (30 days per step). The target inversion parameter is the permeability field, and the initial prior model set (100 realizations) is generated through geological modeling. The observation data include the oil production (WOPR) and water production (WWPR) of each production well, with a total of 200 observation points (2 observations per production well per 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 permeability of the non-channel area (low permeability zone) is 1-100 mD.

[0036] Initialization and macroblock stratification were performed according to the process of step 1; first, the parameter field, observation data and hyperparameters were loaded, and the number of iterations was set to 4, the learning rate was set to 1, and the momentum decay factor was set to 0.9. Then, the reservoir was divided into 12 macroblocks according to the characteristics of the stratigraphic sedimentary phase (river channel, levee, floodplain), and each macroblock covered a continuous geological unit. K-means algorithm clustering was performed within the macroblock according to the permeability gradient, and the macroblock was further subdivided into mesoblocks, finally forming a grid parameter field at the microblock level. Finally, a Gaussian random projection matrix was generated for each macroblock to compress the original parameter space to a low dimension.

[0037] The time sensitivity calculation and sampling were performed according to the process of step 2; first, the parameter-response covariance matrix of each time point was calculated, and the sensitivity of each time point was evaluated by the Frobenius norm. Then, a probability distribution was generated based on the sensitivity, and 5 highly sensitive time windows (each group contained 5 consecutive time points) were sampled without replacement from the 25 time points, covering the stage of drastic changes in oil / water production rate.

[0038] The hierarchical low-rank gain calculation is performed according to the process of step 3; first, the parameters of each macroblock are low-rank projected; then, the low-dimensional response covariance and the low-dimensional Kalman gain matrix are calculated in the low-dimensional space; finally, the low-dimensional Kalman gain matrix is ​​restored to the original space by inverse mapping.

[0039] Perform asynchronous momentum acceleration update according to the process in step 4. For each macroblock, first, update the momentum term using the weighted average of historical gradient directions. Then, update the second-order moment and update the parameters in combination with the adaptive learning rate. The adaptive learning rate is dynamically adjusted based on the second-order moment estimate, and the step size in highly sensitive areas is increased by 30%.

[0040] Perform global synchronization and constraints according to the process in step 5; first, perform parameter field aggregation, that is, merge the update results of each macroblock, eliminate boundary differences through weighted averaging, and ensure the spatial continuity of the permeability field. Then, perform physical constraints, that is, impose non-negative constraints (for example, permeability in the river area ≥ 0) and phase control range (10mD≤permeability in the river area≤500mD), and correct non-physical solutions.

[0041] Iteration termination and result verification are performed according to the process of step 6. First, convergence judgment is performed by calculating the decrease rate of the monitoring objective function (observation residual norm), and the preset threshold (1e-4) is reached after 4 iterations. Then, simulation verification is performed 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%.

[0042] 100 groups of samples and 40 groups of samples of oil and water well production dynamic data were selected for production data fitting experiments, and the fitting results were Figure 3 and Figure 4 shown. Figure 3 and Figure 4 In the figure, the initial a priori model simulates the production dynamic data curve (light gray). After 4 iterations, the fitted generation dynamic data curve (dark gray) closely wraps the actual observed generation dynamic data curve (black). Figure 4 Compared with more samples Figure 3 , the iterative fitting production data matches the actual observed data well, and the fitting error is smaller.

[0043] The convergence of the objective function for 100 groups of samples and 40 groups of samples are as follows: Figure 5 and Figure 6 As shown, Figure 5 and Figure 6 It shows that the historical fitting objective function value decreases with the increase of iteration number and finally converges. Figure 6 Compared with more samples Figure 5 , the objective function is easier to converge and has a smaller value.

[0044] The above results verify the efficiency and stability of the method in large-scale reservoir history matching. The computational complexity is significantly reduced through macroblock stratification and low-rank covariance compression, and the convergence speed is further improved by time adaptive sampling and momentum acceleration, providing a feasible solution for real-time fitting of tens of millions of grid models.

[0045] 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 technicians in this technical field within the essential scope 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