A stepwise data assimilation method for set subspaces in nonlinear inverse problems

Abstract

This invention discloses a stepwise data assimilation method for a set subspace in a nonlinear inverse problem, relating to the field of data processing technology. The invention constructs an initial prior physical set based on the physical state variables to be inverted and optimized, and historical observation data. It extracts the static subspace basis anomaly matrix, initializes the latent variable set and particle weights, calculates fractional-step incremental reweighting, evaluates the current likelihood mismatch penalty using predicted data, updates and normalizes the particle weights in the logarithmic domain, obtains the effective sample number, performs system resampling operations in conjunction with a preset resampling tolerance coefficient, eliminates low-weight particles and replicates high-weight particles, synchronously updates the latent variable set and predicted data, executes a dynamic MCMC mutation loop to generate proposed latent state vectors and affinely maps them to a high-dimensional physical space, and updates the latent variable set by evaluating the annealing target energy within the latent variable subspace until all fractional steps are traversed. Finally, it outputs the latent variable set and maps it back to the physical space.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, specifically to a stepwise data assimilation method for ensemble subspaces for nonlinear inverse problems, applicable to automatic history fitting of oil and gas reservoirs, geological carbon sequestration parameter inversion, and groundwater parameter estimation. Background Technology

[0002] Integrating complex numerical models controlled by nonlinear partial differential equations with sparse, noisy dynamic observation data is a core task of current predictive computational science and a key approach to achieving efficient oil and gas field development, safe monitoring of geological carbon sequestration, and groundwater resource management. This type of problem, inferring the implicit physical state of a system from limited external observations, is known as Bayesian inversion or data assimilation. With the widespread adoption of high-resolution numerical grids, the physical dimensions of unknown parameters (such as the spatially heterogeneous porosity and permeability fields of reservoir rocks) have jumped to an extremely large scale of millions or even tens of millions of dimensions. Simultaneously, due to the complex transport mechanisms of multiphase fluids underground, the system's observation operators typically exhibit extremely strong nonlinearity.

[0003] However, current methods for solving high-dimensional inverse problems and data assimilation have significant limitations. Currently, industry heavily relies on set-based approximation solvers, such as set smoothers with multiple data assimilation and set stochastic maximum likelihood methods. These methods approximate gradient calculations by utilizing the cross-covariance of set samples, avoiding explicit computation of the large Jacobian matrix and significantly improving inversion efficiency. However, these mainstream set smoothing methods heavily rely on the linear Gaussian assumption in their underlying mathematical derivations. When the physical forward modeling operator exhibits strong nonlinear characteristics, the true posterior energy landscape becomes highly nonconvex and exhibits strong ambiguity (i.e., multimodal distribution). Faced with a multimodal posterior distribution, current gradient-based set smoothers reveal mathematical fragility. First, greedy local gradient descent makes them prone to getting trapped in shallow local minima and deadlocking. Second, linear combinations of spatial gradients in opposite directions cancel each other out, forcing the entire set to collapse to the arithmetic mean of the individual solutions, leading to severe mode collapse and a significant underestimation of prediction uncertainty. On the other hand, the pure Bayesian unbiased sampling algorithm has extremely low exploration efficiency when facing high-dimensional parameter spaces, is prone to the curse of dimensionality and weight degradation, and requires the use of an expensive forward numerical simulator for each state evaluation, resulting in high computational costs and making it difficult to meet the needs of practical engineering applications. Summary of the Invention

[0004] This invention aims to address the problems of traditional ensemble data assimilation methods being prone to local minima deadlock, mode collapse, and high computational costs in strongly nonlinear multi-peak inversion. It proposes an ensemble subspace stepwise data assimilation method for nonlinear inverse problems, which is beneficial for accurate inversion to obtain multiple geophysical fields and provides a basis for automatic history fitting of oil and gas reservoirs, inversion of geological carbon sequestration parameters, and estimation of groundwater parameters.

[0005] The present invention adopts the following technical solution: A stepwise data assimilation method for set subspaces in nonlinear inverse problems includes the following steps: Step 1: Obtain the physical state variables to be inverted and optimized and the historical observation data, construct the initial prior physical set, extract the static subspace basis anomaly matrix, and initialize the set of latent variables and the uniform importance weights of all particles; Step 2: Perform fractional-step incremental reweighting, calculate the temperature increment under the current fractional step, evaluate the current likelihood mismatch penalty using the cached prediction data, and complete the stable update and normalization of all particle weights in the logarithmic domain. Step 3: Obtain the number of valid samples in the current prior physical set. When the number of valid samples is lower than the preset resampling tolerance ratio, perform system resampling operation, eliminate low-weight particles and copy high-weight particles, and synchronously update the latent variable set and prediction data. Step 4: Execute a derivativeless dynamic MCMC mutation loop in the latent variable subspace, adaptively adjust the transfer kernel parameters according to the annealing stage of the current fractional step, generate the proposed latent state vector and affine map it to the high-dimensional physical space, call the forward numerical simulator to evaluate the annealing target energy in the latent variable subspace, and execute Metropolis-Hastings to update the latent variable set. Step 5: Execute iterative loop control, repeat steps 2 to 4 until all preset fractional steps are traversed; Step 6: Map the final accepted set of latent variables back to the physical space and output it to obtain the unweighted physical set inversion solution.

[0006] Preferably, in step 1.1, the set size is set. Constructing a system based on observation data that includes The initial prior physical set of each member , ,in, The dimensions of the grid attribute variables after spatial discretization. Given the set of real numbers; calculate the column mean vector of the initial prior physical set. , ,in, , It is a column vector consisting entirely of 1s; Step 1.2, Set the central operator matrix , ,in, , for An identity matrix of order 1. It is the transpose matrix; based on the centralization operator matrix Determine the scaling projection operator matrix , Extracting the static subspace basis anomaly matrix , , ; Step 1.3, Initialize the set of implicit variables , , Used to ensure the reconstruction of equations Strictly valid; Step 1.4: Obtain the column vector of actual historical observation data. , ,in, The observation data dimension is used to construct the observation error covariance matrix based on the tolerance softening mechanism. , The observation error covariance matrix For a diagonal matrix, the first line on its main diagonal is... The observation variance corresponding to each element is The standard deviation is , ,in, It is a function with maximum value. This represents the relative error tolerance ratio. The first column vector of actual observed historical data One portion, This is an absolute minimum error. Step 1.5: Initialize the predicted observation column vectors and uniform importance weights for all particles; The formula for calculating the predicted observation column vector is as follows: ,in, For particle indexing, For the first The predicted observation column vector of each particle. , For the first hidden variable in the set The hidden state column vector corresponding to the column, It is a nonlinear forward modeling operator; the formula for calculating the uniform importance weight is: ,in, For the first Uniform importance weights for each particle.

[0007] Preferably, in step 2, a nonlinear convex temperature scheduling sequence is obtained, which includes multiple fractional temperatures at fractional steps. ,in, The fractional step number, For the first Fractional temperature of each fractional step, Given the total number of fractional steps; calculate the fractional temperature increment for the current fractional step. , ,in, For the first Fractional temperature of each fractional step; The current likelihood mismatch penalty for each particle is evaluated based on the predicted observation column vector and the observation error covariance matrix, calculated using the following formula: ,in, For particle indexing, For the first The current likelihood mismatch penalty for each particle For the first The predicted observation column vector of each particle. , This is a column vector of actual historical observation data. , For the dimensions of the observation data, For the set of real numbers, It is the inverse of the observation error covariance matrix; In the logarithmic domain, the unnormalized logarithmic weights are updated. The formula for calculating the unnormalized logarithmic weights is: ,in, It is the natural logarithm function. For the first The fractional step Unnormalized logarithmic weights of individual particles For the first The fractional step Unnormalized logarithmic weights of each particle; Extracting the global maximum of unnormalized logarithmic weights , ,in, This is the maximum value function; a Log-Sum-Exp offset is performed to prevent floating-point underflow, and the result is updated to the [value]. The normalized weights of each particle at each fractional step are calculated using the following formula: ,in, It is a natural exponential function. For the summation index, The size of the set.

[0008] Preferably, in step 3, the effective sample number, which reflects the number of effective independent particles in the current prior physics set, is calculated based on the normalized weights of all particles in the current prior physics set. , ,in, For particle indexing, For the size of the set, This is the sum of squares of the normalized weights for all particles; when At that time, among them, Using a preset resampling tolerance ratio, a system resampling operation with replacement is performed on the current prior physical set using the cumulative distribution function. The index of each surviving particle is used to overwrite the set of latent variables. Forward prediction vector predicts observation column vector The likelihood mismatch penalty term for each particle is cached, and the weights of each particle are reset to a uniform distribution. The i-th particle in the current prior physics set is then... The weights of each particle are set to , .

[0009] Preferably, in step 4, according to the first Fractional temperature at each fractional step The annealing phase interval performs Markov chain Monte Carlo transfer kernels, adaptively adjusting and setting the number of mutation cycles. Differential evolution probability random walk probability Preconditional Crank-Nicolson autoregressive probability ; Independent execution In each MCMC mutation cycle, a uniformly random number is drawn during the calculation. Generate the proposed hidden state vector, where, , for Uniform distribution; like Then, subspace differential evolution transition is performed, and the formula for calculating the proposed hidden state vector is: ; like If a local random walk is performed, the formula for calculating the proposed hidden state vector is: ; like Then, preconditional Crank-Nicolson autoregression is performed, and the formula for calculating the proposed hidden state vector is: ; In the formula, For the first The proposed hidden state vector of each particle; For the first hidden variable in the set The hidden state column vector corresponding to the column; , , These are the step size scaling factors corresponding to the difference evolution operator, random walk operator, and autoregressive operator, respectively. To satisfy the autoregressive scaling factor defined in the real field, a mandatory constraint is applied. ; , All are sets of latent variables The hidden state vectors of two different particles randomly selected from the sample; It is a perturbation scale factor for differential evolution, used to provide the perturbation capability of forced spatial stripping when particles are highly overlapping; , , All are independent white noise vectors that follow a standard multivariate Gaussian distribution. Specifically, it is a vector with a mean of zero and a covariance of... The multidimensional Gaussian distribution, in which, for An identity matrix of order 1.

[0010] Preferably, in step 4, the annealing target energy assessment and the Metropolis-Hastings acceptance decision process are calculated within the latent variable subspace. If the mutation mechanism is a preconditional Crank-Nicolson autoregression, since its proposed distribution and Gaussian prior reference measure jointly satisfy the detailed stationarity condition, the proportions of the prior density terms strictly cancel each other out in the Metropolis-Hastings acceptance rate calculation, and the current likelihood energy is calculated. and global proposal energy , , ,in, For the first The current hidden state vector of each particle. For the first Each particle corresponds to The predicted observation column vector; This is a column vector of actual historical observation data. , For the dimensions of the observation data, It is the set of real numbers; It is the transpose matrix; It is the inverse of the observation error covariance matrix; If the mutation mechanism is another non-preconditioned Crank-Nicolson autoregressive mutation mechanism, then calculate the global current energy including the prior. and global proposal energy , , ,in, For the first The proposed hidden state vector of each particle; For the first Each particle corresponds to The predicted observation column vector; , The inner product of the current hidden state vector and the proposed hidden state vector is equivalent to the Mahalanobis prior penalty in the high-dimensional physical space under subspace constraints. Draw uniformly random numbers for acceptance judgment. , ,in, for A uniform distribution; if it satisfies If the proposal is accepted, the current state remains unchanged; otherwise, the current state is maintained. Afterward, regardless of acceptance or rejection, the next MCMC mutation is performed within the current fractional step, repeating until the preset target is reached within the current fractional step. After the first mutation, proceed to the next fractional step to continue the dynamic MCMC mutation loop.

[0011] Preferably, in step 5, to prevent spatial exploration stagnation caused by the clones of high-weight particles during differential evolution transitions, an adaptive perturbation expansion mechanism is introduced. When the system resampling operation is triggered in the current fractional step, the white noise vector is dynamically amplified. The corresponding spatial perturbation scale factor forcibly breaks the spatial overlap phenomenon of hidden state particles; if the current fractional step does not trigger the system resampling operation, the differential evolution perturbation scale factor is used.

[0012] The present invention has the following beneficial effects: (1) This invention proposes a stepwise data assimilation method for set subspaces for nonlinear inverse problems. By projecting the subspace, a low-rank anomalous basis spanned by the initial prior set is extracted. The tens of millions-dimensional physical search space is projected into an extremely low-dimensional latent variable subspace. In the latent variable subspace, the Mahalanobis prior energy assessment, which originally involved the inversion of the high-dimensional covariance matrix, is equivalently transformed into the Euclidean L2 norm squared calculation of the latent state vector. After this fidelity-preserving dimensionality reduction process, the initial prior physical set is subjected to an optimal linear approximation (i.e. Gaussianization) that preserves the statistical properties of second-order spatial variation. While overcoming the problems of dimensionality curse and memory explosion, this invention provides a reference measurement basis for the subsequent exploration of multi-peak likelihood landscapes using the dynamic Markov chain Monte Carlo method.

[0013] (2) This invention proposes a stepwise data assimilation method for set subspaces of nonlinear inverse problems. By performing derivative-free differential evolution Markov chain Monte Carlo mutation in the latent space, a dynamic Markov chain Monte Carlo stage control mechanism is introduced. In the early stage of annealing, high-frequency differential evolution gives particles the ability to cross local minima over a wide range. In the later stage of annealing, the probability of the precondition Crank-Nicolson proposal is increased to perform detailed local optimization, thus preserving the multi-solution characteristics of geological parameters.

[0014] (3) This invention proposes a stepwise data assimilation method for ensemble subspaces for nonlinear inverse problems, which decomposes steep nonlinear likelihood functions into smooth fractional steps through incremental likelihood reweighting. At the same time, by introducing a tolerance softening mechanism and an early stopping monitoring mechanism, it effectively prevents the initial weight collapse due to large deviations, while also avoiding the fragmentation of the geophysical field caused by overfitting in the later stages. Attached Figure Description

[0015] Figure 1 This is a flowchart of a stepwise data assimilation method for set subspaces in nonlinear inverse problems according to the present invention.

[0016] Figure 2 This is a schematic diagram of the permeability field of a heterogeneous river channel-type three-dimensional reservoir model.

[0017] Figure 3 The figure shows the dynamic data of oil production in the reservoir when using large-scale ensembles for automatic historical data fitting. In the figure, (a) is the historical data fitting curve of oil production of production well 1, (b) is the historical data fitting curve of oil production of production well 2, (c) is the historical data fitting curve of oil production of production well 3, and (d) is the historical data fitting curve of oil production of production well 4.

[0018] Figure 4 The figure shows the dynamic data of water production in the reservoir when using large-scale ensembles for automatic historical data fitting. In the figure, (a) is the historical data fitting curve of water production of production well 1, (b) is the historical data fitting curve of water production of production well 2, (c) is the historical data fitting curve of water production of production well 3, and (d) is the historical data fitting curve of water production of production well 4.

[0019] Figure 5 This is a graph showing the change of the objective function when using a large-scale ensemble for automatic reservoir history fitting.

[0020] Figure 6The figure shows the dynamic data of oil production in the reservoir when using a small-scale set for automatic historical data fitting. In the figure, (a) is the historical data fitting curve of oil production of production well 1, (b) is the historical data fitting curve of oil production of production well 2, (c) is the historical data fitting curve of oil production of production well 3, and (d) is the historical data fitting curve of oil production of production well 4.

[0021] Figure 7 The figure shows the dynamic data of reservoir water production when using small-scale ensembles for automatic historical data fitting. In the figure, (a) is the historical data fitting curve of water production of production well 1, (b) is the historical data fitting curve of water production of production well 2, (c) is the historical data fitting curve of water production of production well 3, and (d) is the historical data fitting curve of water production of production well 4.

[0022] Figure 8 This is a graph showing the change of the objective function when using a small-scale set for automatic reservoir history fitting.

[0023] Figure 3 , Figure 4 , Figure 6 , Figure 7 In the diagram, gray dots represent the initial prior physical set, blue dots represent the output physical posterior set, and red curves represent the actual reservoir observation data. Detailed Implementation

[0024] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0025] Example 1 This invention proposes a stepwise data assimilation method for set subspaces in nonlinear inverse problems, such as... Figure 1 As shown, this method is applicable to automatic historical data fitting for oil and gas reservoirs, geological carbon sequestration parameter inversion, and groundwater parameter estimation. Specifically, it includes the following steps: Step 1: Obtain the physical state variables and historical observation data to be inverted and optimized, construct the initial prior physical set, extract the static subspace basis anomaly matrix, and initialize the set of latent variables and the uniform importance weights of all particles, including the following sub-steps: Step 1.1, Set the size of the collection To ensure the diversity of geological priors, a prior library containing multiple candidate geological models is loaded, and then, according to the set random extraction mode, samples are blindly drawn from the prior library without replacement. Each member constructs an initial prior physics set. , ,in, The dimensions of the grid attribute variables after spatial discretization. It is the set of real numbers.

[0026] In this embodiment, for grid attribute variables (such as permeability fields) with strict non-negative physical boundary constraints, it is necessary to perform nonlinear mapping transformations such as logarithmic transformation on them before constructing the initial prior physical set, and project them into an unbounded real space to ensure the legality of the physical quantities throughout the dimensionality reduction, inversion and reconstruction.

[0027] Next, calculate the column mean vector of the initial prior physics set. , ,in, , It is a column vector consisting entirely of 1s.

[0028] Step 1.2, Set the central operator matrix , ,in, , for An identity matrix of order 1. It is the transpose matrix; based on the centralization operator matrix Determine the scaling projection operator matrix , Extracting the static subspace basis anomaly matrix , , It is used to mathematically implicitly encode the variation function and connectivity features of the real physical space.

[0029] Step 1.3: Initialize the set of hidden variables in the low-dimensional memory space. , , Used to ensure the reconstruction of equations Strictly valid.

[0030] Step 1.4: Obtain the column vector of actual historical observation data. , ,in, The dimension is the observation data. To mitigate the rapid collapse of weights in the early stages of strongly nonlinear inversion, a tolerance softening mechanism is introduced to construct the observation error covariance matrix. , The observation error covariance matrix For a diagonal matrix, the first line on its main diagonal is... The observation variance corresponding to each element is The standard deviation is , ,in, It is a function with maximum value. This represents the relative error tolerance ratio. The first column vector of actual observed historical data One portion, To mitigate weight collapse during the initial stage of strong nonlinearity, this embodiment uses an absolute minimum error. Set to 5.0, relative error tolerance ratio Set it to 0.1.

[0031] Step 1.5: Initialize the predicted observation column vectors and uniform importance weights for all particles.

[0032] A forward simulation is performed on the physical states corresponding to all hidden vectors to obtain the predicted observation column vectors for all particles. Specifically, the formula for calculating the predicted observation column vectors is as follows: ,in, For particle indexing, For the first The predicted observation column vector of each particle. , For the first hidden variable in the set The hidden state column vector corresponding to the column, It is a nonlinear forward modeling operator; the formula for calculating the uniform importance weight is: ,in, For the first Uniform importance weights for each particle.

[0033] Step 2: Perform fractional-step incremental reweighting, calculate the temperature increment under the current fractional step, evaluate the current likelihood mismatch penalty using the cached prediction data, and complete the stable update and normalization of all particle weights in the logarithmic domain.

[0034] Specifically, a nonlinear convex temperature scheduling sequence is obtained, which includes multiple fractional temperatures at fractional steps. ,in, The fractional step number, For the first Fractional temperature of each fractional step, , The total number of steps is fractional. To control the nonlinear exponent of the heating rate in the early stage of annealing, ; Calculate the fractional temperature increment of the current fractional step. , ,in, For the first Fractional temperature of each fractional step.

[0035] The current likelihood mismatch penalty for each particle is evaluated based on the predicted observation column vector and the observation error covariance matrix, calculated using the following formula: ,in, For particle indexing, For the first The current likelihood mismatch penalty for each particle For the first The predicted observation column vector of each particle. , This is a column vector of actual historical observation data. , For the dimensions of the observation data, For the set of real numbers, It is the inverse of the observation error covariance matrix.

[0036] In the logarithmic domain, the unnormalized logarithmic weights are updated. The formula for calculating the unnormalized logarithmic weights is: ,in, It is the natural logarithm function. For the first The fractional step Unnormalized logarithmic weights of individual particles For the first The fractional step Unnormalized logarithmic weights of each particle.

[0037] Extracting the global maximum of unnormalized logarithmic weights , ,in, This is the maximum value function; to prevent floating-point underflow, this global maximum value is used to perform a shift operation (i.e., Log-Sum-Exp offset) on the exponent term, updating to obtain the nth value. The normalized weights of each particle at each fractional step are calculated using the following formula: ,in, It is a natural exponential function. For the summation index, The size of the set.

[0038] Step 3: Obtain the number of valid samples in the current prior physical set. When the number of valid samples is lower than the preset resampling tolerance coefficient, perform system resampling operation, eliminate low-weight particles and copy high-weight particles, and synchronously update the latent variable set and prediction data.

[0039] In this embodiment, the normalized weights of all particles in the current prior physics set are used to calculate the number of effective samples reflecting the number of effective independent particles in the current prior physics set. , ,in, For particle indexing, For the size of the set, This is the sum of squares of the normalized weights for all particles.

[0040] when In this embodiment, the resampling tolerance ratio is... The preset value is 0.6, indicating that a minority of particles have already dominated the weights. At this point, resampling of the cumulative distribution function of the current prior physics set is triggered, through sampling with replacement. The index of each surviving particle is used to overwrite the set of latent variables. Forward prediction vector predicts observation column vector The likelihood mismatch penalty term for each particle is cached, and the weights of each particle are reset to a uniform distribution. The i-th particle in the current prior physics set is then... The weights of each particle are set to , This eliminates low-weight particles and replicates high-weight particles.

[0041] Step 4: Execute a derivativeless dynamic MCMC (Markov Chain Monte Carlo) mutation loop in the latent variable subspace. Adaptively adjust the transition kernel parameters according to the current annealing stage of the fractional step, generate the proposed latent state vector and affine map it to the high-dimensional physical space, call the forward numerical simulator to evaluate the annealing target energy in the latent variable subspace, and execute Metropolis-Hastings to update the latent variable set.

[0042] In this embodiment, after resampling, a derivative-free dynamic MCMC mutation loop is executed in the latent variable subspace. During the dynamic MCMC mutation loop, according to the first... Fractional temperature at each fractional step The annealing phase interval performs Markov chain Monte Carlo transfer kernels, adaptively adjusting and setting the number of mutation cycles. Differential evolution probability random walk probability Preconditional Crank-Nicolson autoregressive probability .

[0043] Independent execution In each MCMC mutation cycle, a uniformly random number is drawn during the calculation. Generate the proposed hidden state vector, where, , for The uniform distribution.

[0044] like Then, subspace differential evolution transition is performed, and the formula for calculating the proposed hidden state vector is: ; like If a local random walk is performed, the formula for calculating the proposed hidden state vector is: ; like Then, preconditional Crank-Nicolson autoregression is performed, and the formula for calculating the proposed hidden state vector is: ; In the formula, For the first The proposed hidden state vector of each particle; For the first hidden variable in the set The hidden state column vector corresponding to the column; , , These are the step size scaling factors corresponding to the difference evolution operator, random walk operator, and autoregressive operator, respectively. To satisfy the autoregressive scaling factor defined in the real field, a mandatory constraint is applied. ; , All are sets of latent variables The hidden state vectors of two different particles randomly selected from the sample; It is a perturbation scale factor for differential evolution, used to provide the perturbation capability of forced spatial stripping when particles are highly overlapping; , , All are independent white noise vectors that follow a standard multivariate Gaussian distribution. Specifically, it is a vector with a mean of zero and a covariance of... The multidimensional Gaussian distribution, in which, for An identity matrix of order 1.

[0045] In this embodiment, to maintain stable and minimal space exploration, the differential evolution perturbation scale factor can be adjusted. The value is fixed at 0.01, but the differential evolution perturbation scaling factor can be dynamically adjusted based on the number of valid samples in the current fractional step. The value is determined by the following process: if resampling is not triggered, the basic perturbation scale is used; if resampling is triggered, the differential evolution perturbation scale factor is amplified in a stepwise manner. The value of this value is used to forcibly break the spatial overlap phenomenon of hidden state particles.

[0046] Meanwhile, the adaptive adjustment strategy in this embodiment is specifically set as follows: In the early stages of annealing at lower temperatures (i.e.) (At that time), set a lower number of mutation cycles. This gives rise to a higher probability of differential evolution. It is used to dominate the global large-step exploration; as the annealing temperature increases, it enters the mid-annealing stage (i.e. (At that time), gradually increase the number of mutation cycles. At the same time, reduce the probability of differential evolution. And improve the preconditioning Crank-Nicolson autoregressive probability. Balancing global and local considerations, and strengthening the fine-grained local development in the region approximating the true solution; with fractional temperature The temperature rises to the late annealing stage (i.e., (Time), set the highest frequency mutation cycle number. Significantly increases the preconditioning Crank-Nicolson autoregressive probability Deep local fine-tuning is performed in local regions that approximate the true solution, and the step size scaling factor corresponding to the synchronously decreasing differential evolution operator, random walk operator, and autoregressive operator is adjusted to maintain a reasonable acceptance rate for mutation proposals.

[0047] The annealing target energy assessment and Metropolis-Hastings acceptance decision process are calculated within the latent variable subspace. If the mutation mechanism is a preconditional Crank-Nicolson autoregression, since its proposed distribution and Gaussian prior reference measure jointly satisfy the detailed stationarity condition, the proportions of the prior density terms strictly cancel each other out in the Metropolis-Hastings acceptance rate calculation. The current likelihood energy is then calculated. and global proposal energy , , ,in, For the first The current hidden state vector of each particle. For the first Each particle corresponds to The predicted observation column vector; This is a column vector of actual historical observation data. , For the dimensions of the observation data, It is the set of real numbers; It is the transpose matrix; It is the inverse of the observation error covariance matrix.

[0048] If the mutation mechanism is another non-preconditioned Crank-Nicolson autoregressive mutation mechanism, then calculate the global current energy including the prior. and global proposal energy , , ,in, For the first The proposed hidden state vector of each particle; For the first Each particle corresponds to The predicted observation column vector; , The inner product of the current hidden state vector and the proposed hidden state vector is approximately equivalent to the Mahalanobis prior penalty in the high-dimensional physical space under subspace constraints.

[0049] Draw uniformly random numbers for acceptance judgment. , ,in, for A uniform distribution; if it satisfies If the proposal is accepted, it is accepted; otherwise, the current state remains unchanged. Afterward, regardless of acceptance or rejection, the next MCMC mutation within the current fractional step is initiated. This process is repeated until the set parameters are accumulated within the current fractional step. The mutation proceeds to the next fractional step and continues the dynamic MCMC mutation cycle.

[0050] Step 5: Execute iterative loop control, repeat steps 2 to 4 until all preset fractional steps are traversed.

[0051] In this embodiment, after the MCMC mutation calculation of a single fractional step is completed, the system calculates the relative decrease rate of the average objective function value of the entire set in real time. If the relative decrease rate is detected to be lower than a preset threshold in the continuous state and the current average objective function value has approached or fallen below the preset safe physical tolerance limit, the system determines that the optimal balance has been reached between the data fit and the prior physical constraints. In order to prevent geological field fragmentation caused by overfitting, the system will trigger an early stop command to exit the loop in advance. If the early stop is not triggered, the fractional step is updated and the dynamic MCMC mutation loop calculation is continued in the next fractional step.

[0052] Step 6: Map the final accepted set of latent variables back to the physical space and output it to obtain the unweighted physical set inversion solution.

[0053] In this embodiment, after the dynamic MCMC mutation cycle calculation terminates, the final set of latent variables will be accepted and retained. Through the affine mapping formula Accurate affine reconstruction back to the higher-dimensional physical space, in which... Given a physical posterior set, the output is an approximate physical posterior set controlled by regularization (i.e., an unweighted large-scale high-dimensional multi-peak physical set inversion solution). This set effectively preserves the multi-solution characteristics of geological parameters and effectively avoids mode collapse.

[0054] Example 2 To demonstrate the feasibility and superiority of the method of this invention in solving large-scale nonlinear inverse problems, the method of this invention is applied to the automatic history fitting and high-dimensional full-field permeability inversion of a classic heterogeneous channel-type three-dimensional reservoir model.

[0055] In this embodiment, the permeability field of the heterogeneous river channel-type three-dimensional reservoir model is as follows: Figure 2 As shown, the physical grid system of the heterogeneous channel-type 3D reservoir model is set to 60×60×1, with a total of 3600 spatial grids planned for the entire area. The heterogeneous channel-type 3D reservoir model includes one water injection well and four production wells. The water injection well is located at the center of the heterogeneous channel-type 3D reservoir model, while the production wells are evenly distributed around the water injection well at the four corners of the heterogeneous channel-type 3D reservoir model. The observation data of the heterogeneous channel-type 3D reservoir model includes a total of 50 time steps, with each time step corresponding to a duration of 30 days. A total of eight dynamic fitting indicators are used to monitor the oil production and water production of the four production wells, forming a complex nonlinear observation response vector with a total dimension of 400.

[0056] To verify the adaptive processing capability of the method of the present invention for sets of different sizes and the intelligent advantages of the early stopping mechanism, this embodiment sets different set sizes. Large set of 100 and set size For a small set of 40, two parallel control groups were used, and the total number of fractional steps was set. It is 16.

[0057] In this embodiment, high-dimensional physical prior initialization and subspace extraction are performed according to the process in step 1. First, 500 models are loaded from the candidate geological model library as the total prior library. An initial prior physical set matrix is ​​constructed by randomly selecting 100 and 40 members respectively. Then, the static subspace basis anomaly matrix is ​​extracted, mapping the high-dimensional physical space to the low-dimensional latent variable subspace, and the latent variable set and observation error covariance matrix are initialized.

[0058] Then, perform fractional-step annealing and likelihood increment reweighting according to step 2. Calculate the fractional temperature increment of the current fractional step based on the set nonlinear convex temperature scheduling sequence, reuse the cached prediction data to evaluate the current likelihood mismatch penalty, and complete the stable update and normalization of the weights in the logarithmic domain.

[0059] Then, follow the process in step 3 to perform dynamic ESS monitoring and system resampling. Using the normalized probability weights, calculate the effective sample number of the effective independent particle count in the current prior physics set in real time. When the system detects that a small number of particles dominate the weights (i.e., below the resampling tolerance ratio), it triggers cumulative distribution function resampling, eliminating inferior particles and cloning high-probability superior particles, while resetting the weights to a uniform distribution.

[0060] Then, following the process in step 4, a low-dimensional MCMC mutation loop with dynamic parameter tuning is performed. In this embodiment, after resampling, derivative-free mutation is performed in the latent variable subspace. The system adaptively schedules parameters according to the current annealing temperature; that is, in the early stage of annealing, a larger differential evolution probability is used to drive particles to perform global large-step transitions; in the middle and late stages, the preconditioned Crank-Nicolson autoregressive probability and the number of mutation loops are automatically increased, and in-depth fine-tuning is performed in the local region that approximates the true solution.

[0061] Following the process in step 5, iterative loop control and early stopping mechanism judgment are performed. After a single fractional-step mutation, the system calculates the relative decrease rate of the average objective function value of the entire set in real time. If the decrease rate is lower than the set minimum threshold and the objective function value has met the physical tolerance limit, an early stopping instruction is triggered to exit the loop in advance to prevent geological field fragmentation; otherwise, the fractional steps are incremented to enter the next round of MCMC mutation loop iteration.

[0062] Finally, following the process in step 6, the final deconstruction and physical space mapping are performed. After the loop terminates or an early stop is triggered, the final set of latent variables is precisely reconstructed back to the high-dimensional physical space through affine mapping, and the physical posterior set is output to obtain the unweighted physical set inversion solution.

[0063] Furthermore, comparing the automatic reservoir history fitting results using large-scale and small-scale sets in this embodiment, we obtain the following:

[0064] For large-scale datasets, when using the method of this invention for automatic reservoir history fitting, the production dynamic data of each production well is as follows: Figure 3 and Figure 4 As shown, analysis Figure 3 and Figure 4 It can be seen that as the annealing temperature increases, the system adaptively executes dynamic MCMC parameter scheduling, the number of effective samples remains within a healthy range, and the blue envelope band fits the actual red observation data well. Combined with... Figure 5 As shown in the graph of the objective function variation of large-scale sets, large-scale sets have extremely strong global exploration and fitting capabilities in complex posterior spaces.

[0065] For small-scale datasets, the production dynamics data of each production well during automatic reservoir history fitting using the method of this invention are as follows: Figure 6 and Figure 7 As shown, the analysis is as follows Figure 6 and Figure 7 It can be seen that even using only 40 samples and terminating early, the predicted envelope after inversion still closely matches the historical curve, combined with... Figure 8The graph showing the objective function changes for a small set reveals that its stepwise descent remains stable. While the objective function decreases rapidly in the early stages of annealing, the rate of decline slows significantly in the later stages. By the 12th fractional step, the system detects that the average objective function has met the physical tolerance limit, and the relative descent rate has dropped to an extremely low level. At this point, the system determines that continued assimilation would lead to overfitting, and the early stopping mechanism successfully intervenes and forcibly terminates the loop.

[0066] The above comparison fully verifies the efficiency of the method of this invention in handling high-dimensional, strongly nonlinear inverse problems. At the same time, the method of this invention successfully maintains the validity of the posterior set by utilizing the dynamic MCMC mechanism, avoids mode collapse, and better identifies the risk of overfitting by combining with the early stopping mechanism. While taking into account the computational cost, it protects the physical meaning of the inverted geological field and successfully solves the problems of traditional set subspace step-by-step data assimilation methods being prone to local minimum deadlock, mode collapse, and high computational cost of standard sampling methods in strongly nonlinear multi-peak inversion.

[0067] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A stepwise data assimilation method for set subspaces in nonlinear inverse problems, characterized in that, Includes the following steps: Step 1: Obtain the physical state variables to be inverted and optimized and the historical observation data, construct the initial prior physical set, extract the static subspace basis anomaly matrix, and initialize the set of latent variables and the uniform importance weights of all particles; Step 2: Perform fractional-step incremental reweighting, calculate the temperature increment under the current fractional step, evaluate the current likelihood mismatch penalty using the cached prediction data, and complete the stable update and normalization of all particle weights in the logarithmic domain. Step 3: Obtain the number of valid samples in the current prior physical set. When the number of valid samples is lower than the preset resampling tolerance ratio, perform system resampling operation, eliminate low-weight particles and copy high-weight particles, and synchronously update the latent variable set and prediction data. Step 4: Execute a derivativeless dynamic MCMC mutation loop in the latent variable subspace, adaptively adjust the transfer kernel parameters according to the annealing stage of the current fractional step, generate the proposed latent state vector and affine map it to the high-dimensional physical space, call the forward numerical simulator to evaluate the annealing target energy in the latent variable subspace, and execute Metropolis-Hastings to update the latent variable set. Step 5: Execute iterative loop control, repeat steps 2 to 4 until all preset fractional steps are traversed; Step 6: Map the final accepted set of latent variables back to the physical space and output it to obtain the unweighted physical set inversion solution.

2. The method for stepwise data assimilation of set subspaces for nonlinear inverse problems according to claim 1, characterized in that, Step 1.1, Set the size of the collection Constructing a system based on observation data that includes The initial prior physical set of each member , ,in, The dimensions of the grid attribute variables after spatial discretization. Given the set of real numbers; calculate the column mean vector of the initial prior physical set. , ,in, , It is a column vector consisting entirely of 1s; Step 1.2, Set the central operator matrix , ,in, , for An identity matrix of order 1. It is a transpose matrix; based on the centralization operator matrix Determine the scaling projection operator matrix , Extracting the static subspace basis anomaly matrix , , ; Step 1.3, Initialize the set of implicit variables , , Used to ensure the reconstruction of equations Strictly valid; Step 1.4: Obtain the column vector of actual historical observation data. , ,in, The observation data dimension is used to construct the observation error covariance matrix based on the tolerance softening mechanism. , The observation error covariance matrix For a diagonal matrix, the first line on its main diagonal is... The observation variance corresponding to each element is The standard deviation is , ,in, It is a function with maximum value. This represents the relative error tolerance ratio. The first column vector of actual observed historical data One portion, This is an absolute minimum error. Step 1.5: Initialize the predicted observation column vectors and uniform importance weights for all particles; The formula for calculating the predicted observation column vector is as follows: ,in, For particle indexing, For the first The predicted observation column vector of each particle. , The first in the set of latent variables The hidden state column vector corresponding to the column, It is a nonlinear forward modeling operator; the formula for calculating the uniform importance weight is: ,in, For the first Uniform importance weights for each particle.

3. The method for stepwise data assimilation of set subspaces for nonlinear inverse problems according to claim 1, characterized in that, In step 2, a nonlinear convex temperature scheduling sequence is obtained, which includes multiple fractional temperatures at fractional steps. ,in, The fractional step number, For the first Fractional temperature of each fractional step, Given the total number of fractional steps; calculate the fractional temperature increment for the current fractional step. , ,in, For the first Fractional temperature of each fractional step; The current likelihood mismatch penalty for each particle is evaluated based on the predicted observation column vector and the observation error covariance matrix, calculated using the following formula: ,in, For particle indexing, For the first The current likelihood mismatch penalty for each particle For the first The predicted observation column vector of each particle. , This is a column vector of actual historical observation data. , For the dimensions of the observation data, For the set of real numbers, It is the inverse of the observation error covariance matrix; In the logarithmic domain, the unnormalized logarithmic weights are updated. The formula for calculating the unnormalized logarithmic weights is: ,in, It is the natural logarithm function. For the first The fractional step Unnormalized logarithmic weights of individual particles For the first The fractional step Unnormalized logarithmic weights of each particle; Extracting the global maximum of the unnormalized logarithmic weights , ,in, This is the maximum value function; a Log-Sum-Exp offset is performed to prevent floating-point underflow, and the result is updated to the [value]. The normalized weights of each particle at each fractional step are calculated using the following formula: ,in, It is a natural exponential function. For the summation index, The size of the set.

4. The method for stepwise data assimilation of set subspaces for nonlinear inverse problems according to claim 1, characterized in that, In step 3, the number of effective samples reflecting the number of effective independent particles in the current prior physics set is calculated based on the normalized weights of all particles in the current prior physics set. , ,in, For particle indexing, For the size of the set, This is the sum of squares of the normalized weights for all particles; when At that time, among them, Using a preset resampling tolerance ratio, a system resampling operation with replacement is performed on the current prior physical set using the cumulative distribution function. The index of each surviving particle is used to overwrite the set of latent variables. Forward prediction vector predicts observation column vector The likelihood mismatch penalty term for each particle is cached, and the weights of each particle are reset to a uniform distribution. The i-th particle in the current prior physics set is then... The weights of each particle are set to , .

5. The method for stepwise data assimilation of set subspaces for nonlinear inverse problems according to claim 1, characterized in that, In step 4, according to the first Fractional temperature at each fractional step The annealing phase interval performs Markov chain Monte Carlo transfer kernels, adaptively adjusting and setting the number of mutation cycles. Differential evolution probability random walk probability Preconditional Crank-Nicolson autoregressive probability ; Independent execution In each MCMC mutation cycle, a uniformly random number is drawn during the calculation. Generate the proposed hidden state vector, where, , for Uniform distribution; like Then, subspace differential evolution transition is performed, and the formula for calculating the proposed hidden state vector is: ； like If a local random walk is performed, the formula for calculating the proposed hidden state vector is: ； like Then, preconditional Crank-Nicolson autoregression is performed, and the formula for calculating the proposed hidden state vector is: ； In the formula, For the first The proposed hidden state vector of each particle; The first in the set of latent variables The hidden state column vector corresponding to the column; , , These are the step size scaling factors corresponding to the difference evolution operator, random walk operator, and autoregressive operator, respectively. To satisfy the autoregressive scaling factor defined in the real field, a mandatory constraint is applied. ; , All are sets of latent variables The hidden state vectors of two different particles randomly selected from the sample; It is a perturbation scale factor for differential evolution, used to provide the perturbation capability of forced spatial stripping when particles are highly overlapping; , , All are independent white noise vectors that follow a standard multivariate Gaussian distribution. Specifically, it is a vector with a mean of zero and a covariance of... The multidimensional Gaussian distribution, in which, for An identity matrix of order 1.

6. The method for stepwise data assimilation of set subspaces for nonlinear inverse problems according to claim 5, characterized in that, In step 4, the annealing target energy assessment and the Metropolis-Hastings acceptance decision process are calculated within the latent variable subspace. If the mutation mechanism is a preconditional Crank-Nicolson autoregression, since its proposed distribution and Gaussian prior reference measure jointly satisfy the detailed stationarity condition, the proportions of the prior density terms strictly cancel each other out in the Metropolis-Hastings acceptance rate calculation. The current likelihood energy is then calculated. and global proposal energy , , ,in, For the first The current hidden state vector of each particle. For the first Each particle corresponds to The predicted observation column vector; This is a column vector of actual historical observation data. , For the dimensions of the observation data, It is the set of real numbers; It is the transpose matrix; It is the inverse of the observation error covariance matrix; If the mutation mechanism is another non-preconditioned Crank-Nicolson autoregressive mutation mechanism, then calculate the global current energy including the prior. and global proposal energy , , ,in, For the first The proposed hidden state vector of each particle; For the first Each particle corresponds to The predicted observation column vector; , The inner product of the current hidden state vector and the proposed hidden state vector is equivalent to the Mahalanobis prior penalty in the high-dimensional physical space under subspace constraints. Draw uniformly random numbers for acceptance judgment. , ,in, for A uniform distribution; if it satisfies If the proposal is accepted, the current state remains unchanged; otherwise, the current state is maintained. Afterward, regardless of acceptance or rejection, the next MCMC mutation is performed within the current fractional step, repeating until the preset target is reached within the current fractional step. After the first mutation, proceed to the next fractional step to continue the dynamic MCMC mutation loop.

7. The method for stepwise data assimilation of set subspaces for nonlinear inverse problems according to claim 6, characterized in that, In step 5, to prevent spatial exploration stagnation caused by the clones of high-weighted particles during differential evolution transitions, an adaptive perturbation expansion mechanism is introduced. When the system resampling operation is triggered in the current fractional step, the white noise vector is dynamically amplified. The corresponding spatial perturbation scale factor forcibly breaks the spatial overlap phenomenon of hidden state particles; if the current fractional step does not trigger the system resampling operation, the differential evolution perturbation scale factor is used.

Images