A method for reservoir history matching based on information-constrained stein transport

Abstract

The application provides a reservoir history matching method based on information constraint Stein transport, and relates to the technical fields of reservoir numerical simulation and dynamic data assimilation, and specifically comprises the following steps: obtaining a prior set, actual observation data, an observation error standard deviation, observation grouping, a state variable local block and a composite mapping; constructing a state anomaly matrix and a predicted observation anomaly matrix; performing redundant correction on the diagonal observation error variance of the observation; constructing a local observation whitening weighting matrix and an information constraint observation projection matrix; constructing a Stein transport update direction, and updating the set members of the state variable local block; merging the update results, and performing the composite mapping once after the completion of each backfire layer until the output of a posterior set and a posterior predicted observation set. The technical scheme of the application overcomes the problems in the prior art that the effective independent information amount of each group of observation cannot be identified, the redundant observation cannot be compressed into a few independent information modes, and the set cannot be updated using Stein transport.

Description

Technical Field

[0001] This invention relates to the field of reservoir numerical simulation and dynamic data assimilation technology, specifically to a reservoir history fitting method based on information-constrained Stein transport. Background Technology

[0002] Ensemble Smoother Multiple Data Assimilation (ESMDA) is a commonly used automatic assimilation method in reservoir history fitting. This type of method constructs a prior set to represent reservoir parameter uncertainties, gradually introduces observational information into the model ensemble through multi-level tempering assimilation, and uses linear updates in the form of Kalman gain to correct ensemble members. Localized ESMDA further utilizes spatial distance or sample correlation to truncate the update relationship, mitigating spurious correlations caused by finite ensemble size.

[0003] However, in real reservoir history fitting, injection-production dynamics observations typically exhibit dense time series. The oil production, water production, water cut, and gas production of the same production well, or the actual water injection rate and bottom hole pressure of the same injection well, are highly correlated between adjacent time points. Increasing the density of an injection-production dynamics curve from 15 time points to 60 time points does not necessarily result in a fourfold increase in independent information. However, in practical engineering, the observation error covariance matrix is ​​usually only available in diagonal form. When traditional ESMDA or localized ESMDA processes such dense time series, it easily approximates highly correlated time points as independent observations, leading to redundant and repeated counts.

[0004] Redundant observation re-counting leads to the assimilation operator repeatedly injecting similar information into the posterior set, causing a continuous decrease in the posterior set variance and even set collapse. Set collapse significantly underestimates posterior uncertainty and may lead to overfitting of the model to dense observations with correlated errors. While localization methods can suppress spurious correlations caused by finite sets, they cannot identify the effective independent information within the same time series, and therefore cannot fundamentally solve the problem of re-counting dense time observations.

[0005] Existing observation error diagnosis methods, such as error estimation methods based on prior innovation or Desroziers' ideas, typically rely on the statistical characteristics of observation residuals. These methods are suitable for checking whether the observation error variance is set too large or too small, but not for determining how many independent information directions the assimilation operator can provide during the posterior update phase. Set collapse caused by redundant observations mainly occurs during the assimilation update process, not during the prior residual statistics phase; therefore, simply relying on the magnitude of prior innovation is insufficient to accurately identify the effective independent information content of dense observations.

[0006] Therefore, there is a need for a reservoir history fitting method based on information-constrained Stein transport that can identify the effective independent information content of each set of observations using response correlation structure and posterior influence spectrum, compress dense and redundant observations into a few independent information patterns, update the set using Stein transport with kernel exclusion term in the compressed observation space, and recover the necessary variance in the data null space of the compressed observation operator. Summary of the Invention

[0007] The main objective of this invention is to provide a reservoir history fitting method based on information-constrained Stein transport, in order to solve the problems in the prior art that cannot use response correlation structure and posterior influence spectrum to identify the effective independent information content of each set of observations, cannot compress dense and redundant observations into a few independent information patterns, cannot use Stein transport with kernel exclusion term to update the set in the compressed observation space, and cannot recover the necessary variance in the data null space of the compressed observation operator.

[0008] To achieve the above objectives, this invention provides a reservoir history fitting method based on information-constrained Stein transport, specifically including the following steps: S1 obtains the prior set of reservoir state variables to be updated, actual observation data, standard deviation of observation error, observation grouping, local blocks of state variables, and composite mapping between reservoir model and observation operator.

[0009] S2, perform a composite mapping on the prior set to obtain the predicted observation set, and construct the state anomaly matrix and the predicted observation anomaly matrix based on the prior set and the predicted observation set, respectively.

[0010] S3, based on the predicted observation anomaly matrix, estimates the number of effective independent observations for each observation group according to the observation group, and performs redundancy correction on the diagonal observation error variance of the observations based on the number of effective independent observations.

[0011] S4. For each local block of state variable, select observations related to the local block of state variable, construct a local observation whitening weighted matrix, and construct an information-constrained observation projection matrix based on the local observation whitening weighted matrix, compressing the local observation residual matrix and the local predicted observation anomaly matrix into an independent information pattern space.

[0012] S5. Construct the Stein transport update direction in the independent information pattern space, and update the set members of the local block of state variables according to the Stein transport update direction.

[0013] S6. Construct a data null space based on the compressed prediction observation anomaly matrix. Perform variance recovery on the updated set members in the data null space, merge the update results of all local blocks of state variables, and run a composite mapping once after each tempering layer is completed until the posterior set and the posterior prediction observation set are output.

[0014] Furthermore, the specific steps of S1 are as follows: S1.1, Determine the reservoir state variables to be updated. , for The state variables are: permeability, porosity, fault multiplier, capillary pressure parameter, relative permeability parameter, and well control quantity.

[0015] S1.2, Generate a priori set based on geological understanding, well logging interpretation, seismic inversion results, dynamic understanding, or existing historical fitting results. Represented as: ; in, For the first A set of members; This represents the total number of members in the set.

[0016] S1.3, Obtain actual observation data and set the standard deviation of observation error. Actual observation data includes: oil production, water production, water cut, gas production, water injection volume, water injection rate, bottom hole pressure, pressure recovery interpretation results, and seismic attributes.

[0017] S1.4 Observation groups are divided according to well number, data type, seismic attribute window or connected seismic region, and local blocks of state variables are divided according to grid block, well control area, flow unit, layer segment or fault block.

[0018] S1.5, establish a composite mapping between the reservoir model and the observation operator.

[0019] Furthermore, the specific steps of S2 are as follows: S2.1, Perform a composite mapping on each set member: ; in, For the first A set of members; For the first Predicted observation vectors for each set member; It is a composite mapping function; Obtain the predicted observation set for: ; in, For the first The predicted observation vectors of each set member.

[0020] S2.2, using the mean-removing method and... Scaled projection matrix Construct the state anomaly matrix and the predicted observation anomaly matrix : ; .

[0021] Furthermore, the specific steps of S3 include: S3.1, for the first For each observation group, extract the corresponding prediction observation anomaly submatrix from the prediction observation anomaly matrix. : ; in, For the set of observation indexes; It is the set of real numbers; For the first Number of observations in each observation group; To predict the observation anomaly matrix from a global perspective Selecting row indexes belongs to The submatrix obtained by taking all rows and retaining all N columns has a dimension of . ; Index for observation groups.

[0022] S3.2, using the predicted observation anomaly submatrix Constructing the whitening prediction observation anomaly matrix : ; in, This is the operator for constructing diagonal matrices; For the first Standard deviation of observation error for each observation data; Index for observation data.

[0023] S3.3, will Normalize each row to the unit norm, and we get... And calculate the correlation matrix of the observed response. : ; in, For transpose; To predict the whitening observation anomaly matrix The matrix obtained by normalizing each non-zero row according to the L2 norm; The observed response correlation matrix is ​​shrunk, and the shrunk observed response correlation matrix is ​​calculated. : ; in, The shrinkage coefficient; It is an identity matrix.

[0024] S3.4, for Find the eigenvalues ​​and calculate the posterior influence strength. : ; in, for eigenvalues; Index for source members.

[0025] Furthermore, S3 also includes: S3.5, Calculate the number of valid independent observations based on the posterior influence strength. And then cut it: ; ; in, To obtain the minimum value.

[0026] S3.6, for the first Each observation group defines the redundancy expansion factor. for: ; in, To obtain the maximum value.

[0027] S3.7, Correct the diagonal observation error variance for the observation data within this observation group: ; in, For the first The observation error variance after correction for each observation data; For the first Redundancy inflation factor of the group to which each observation data belongs; For the first The original error variance of each observation data; Indicates the first The observation group to which each observation data belongs; Corrected diagonal observation error matrix for: ; in, For the first The original error variance of the observation data.

[0028] Furthermore, the specific steps of S4 are as follows: S4.1, for each local block of state variables, select a set of observation indices related to the local block of state variables based on sample correlation strength, physical distance, well control relationship, flow connectivity relationship, or geological barrier relationship. Select the local block that affects the state variable Observation index set : ; in, The correlation truncation threshold; For the local block of state variables and the first The strength of the sample correlation between individual observation data; To take the absolute value.

[0029] S4.2, according to Determine localization weights : .

[0030] S4.3, based on localization weights and the corrected standard deviation of observation error Constructing a local observation whitening weighted matrix : .

[0031] S4.4 uses the local observation whitening weighted matrix to process local predicted observation anomalies, calculates the number of local effective independent observations, and determines the local compressed dimension.

[0032] Furthermore, the specific steps of S4.4 are as follows: S4.4.1, for any local block of state variables First, extract the local blocks related to the state variables. Related predictive observation anomaly submatrix : ; in, To predict the observation anomaly matrix from a global perspective Selecting an index set The local submatrix obtained by taking all the corresponding rows and retaining all columns has a dimension of . .

[0033] S4.4.2 Constructing a local whitening prediction observation anomaly matrix using a local observation whitening weighted matrix. : .

[0034] S4.4.3, based on The row correlation matrix and posterior influence spectrum are used to calculate the number of locally effective independent observations. ,according to Take the local compression dimension : ; in, This is the floor symbol.

[0035] S4.4.4, retain the former eigenvectors Define information-constrained observation projection matrix : .

[0036] S4.4.5, for any predicted observation vector The local compression residual is : ; in, This is the actual observation vector composed of all historical observations.

[0037] S4.4.6, using the information-constrained observation projection matrix, the local predicted observation anomaly submatrix is ​​projected onto the independent information pattern space to obtain the locally compressed predicted observation anomaly matrix. : .

[0038] Furthermore, the specific steps of S5 are as follows: S5.1, construct a data attraction term based on the local compression residuals and the local compression predicted observation anomaly matrix, for the There are one tempering layer, and the tempering increment is... Local blocks of state variables The Middle Data attraction items of a set member for: ; in, This is a local block state anomaly matrix; This is a transpose.

[0039] S5.2 Construct a priori anchoring terms based on the deviations between the current set members and the prior set members. : ; in, For prior anchoring strength; During the current tempering update phase, the first The set members are in the local block of the state variable. State subvectors on; For the first The set members are in the local block of the state variable. The fixed reference state subvector on, This represents the prior state.

[0040] S5.3 standardizes the state variables within the local block of the state variables, assuming the first... The standard deviation of each state variable in the current set is , define the first Standardized state variables of each set member for: ; in, Local block of state variables for the current tempering phase The set mean vector; This is the index for the global state variable.

[0041] S5.4 Construct the Stein kernel matrix based on the standardized set membership distances. : ; in, It is a natural exponential function; For bandwidth parameters; It is an L2 norm.

[0042] S5.5 Calculate the Stein transport velocity based on the data attraction term, the prior anchoring term, and the kernel repulsion term in the Stein kernel matrix. ; ; in, The metric matrix corresponding to the standardized distance. ; For the first In the first tempering stage The set members are in the local block of the state variable. The state subvectors on.

[0043] S5.6 Select the update step size based on the linearized prediction residual of the current set, and use the update step size and Stein transport speed to obtain the temporary update set of the local block of state variables.

[0044] Furthermore, the specific steps of S6 are as follows: S6.1, using local compression to predict the observed anomaly matrix to calculate the information-constrained linear Gaussian target transformation. : ; in, for The identity matrix; Calculate the target variance trace : ; in, Let be the trace of the matrix.

[0045] S6.2, Suppose the temporarily updated state variable local block anomaly is... Its variance trace for: .

[0046] S6.3, if Then for Perform singular value decomposition: ; in, This is the left singular vector orthogonal matrix of the locally compressed predicted observation anomaly matrix; It is an orthogonal matrix of right singular vectors; It is a rectangular diagonal matrix composed of non-negative singular values; This represents the minimum variance retention ratio for dimensionless quantities.

[0047] Furthermore, S6 also includes: S6.4, Preserve the right singular vector corresponding to the data activity subspace. Define the null space projection of data in the set member space. And magnify the set anomaly along the data null space projection direction: ; in, It is a vector consisting entirely of 1s.

[0048] S6.5 If a state variable belongs to multiple state variable local blocks, the update results of multiple state variable local blocks are merged according to the following formula: ; in, Weights for merging local blocks of state variables; For state variables, after Stein transport update and data null space variance recovery are completed for local blocks of state variables; To complete the local transport, variance recovery, and weighted merging of overlapping blocks in the current tempered layer, the first... The first set member A global state variable.

[0049] S6.6, after merging all local blocks of state variables, the updated set is obtained. And run a numerical simulation: ; in, For the first The posterior predictive observation set of the tempered layer output; For numerical simulation functions; right Repeat steps S2.2 to S6.6, obtaining an updated posterior set in each round. and its posterior predictive observation set After the loop ends, the final posterior set is output. and posterior prediction observation set The total number of orthogonal runs is: ; in, This represents the number of tempering layers.

[0050] The present invention has the following beneficial effects: This invention estimates the effective independent observations of an observation group using posterior influence spectrum, automatically identifying redundant time points in dense production curves and avoiding directly treating dense sampling points as a large number of independent observations. It eliminates the need to construct a complete observation error correlation matrix, correcting the observation information content using a redundancy expansion factor even under diagonal observation error conditions. It compresses locally correlated observations into independent information patterns using a local observation whitening weighting matrix and an information-constrained observation projection matrix, reducing redundant constraints on ensemble updates. It employs Stein transport instead of Kalman gain as the primary update mechanism, maintaining the dispersion of ensemble members through a kernel exclusion term, thus reducing the risk of ensemble collapse. It recovers the ensemble variance in the null space of the compressed observation operator, improving posterior ensemble reliability while minimizing disruption to historical fitting effects. Each tempering layer only requires one full ensemble forward model, maintaining the same total number of forward model calls as the classic ESMDA. Attached Figure Description

[0051] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In the drawings: Figure 1 The overall flowchart of the reservoir history fitting method based on information-constrained Stein transport of the present invention is shown.

[0052] Figure 2The diagram shows a comparison of the posterior marginal Wasserstein distance in three types of redundant observation verification experiments.

[0053] Figure 3 The diagram shows a comparison of the proportions of posterior variance in three types of redundant observation verification experiments of the present invention.

[0054] Figure 4 The diagram shows a comparison of the predicted normalized root mean square error in three types of redundant observation verification experiments of the present invention.

[0055] Figure 5 The diagram shows the original observation mismatch comparison in the three types of redundant observation verification experiments of this invention.

[0056] Figure 6 The diagram shows a comparison of the 90% edge coverage rate of the present invention in three types of redundant observation verification experiments.

[0057] Figure 7 The diagram shows the initial set, the fitted set, and the matching results of historical observations of the actual water injection rate of water injection well 2 in the reservoir numerical simulation experiment of this invention.

[0058] Figure 8 The figure shows the initial set, the fitted set, and the matching results of historical observations of the oil production rate of production well 8 in the reservoir numerical simulation experiment of this invention.

[0059] Figure 9 The diagram shows the convergence results of the objective function of all set members in the initial batch and four updated batches in the reservoir numerical simulation experiment of this invention. Detailed Implementation

[0060] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0061] like Figure 1 The reservoir history fitting method based on information-constrained Stein transport, as shown, specifically includes the following steps: S1 obtains the prior set of reservoir state variables to be updated, actual observation data, standard deviation of observation error, observation grouping, local blocks of state variables, and composite mapping between reservoir model and observation operator.

[0062] S2, perform a composite mapping on the prior set to obtain the predicted observation set, and construct the state anomaly matrix and the predicted observation anomaly matrix based on the prior set and the predicted observation set, respectively.

[0063] S3, based on the predicted observation anomaly matrix, estimates the number of effective independent observations for each observation group according to the observation group, and performs redundancy correction on the diagonal observation error variance of the observations based on the number of effective independent observations.

[0064] S4. For each local block of state variable, select observations related to the local block of state variable, construct a local observation whitening weighted matrix, and construct an information-constrained observation projection matrix based on the local observation whitening weighted matrix, compressing the local observation residual matrix and the local predicted observation anomaly matrix into an independent information pattern space.

[0065] S5. Construct the Stein transport update direction in the independent information pattern space. The Stein transport update direction includes data attraction terms, prior anchor terms, and kernel repulsion terms. Update the set members of the local block of state variables according to the Stein transport update direction.

[0066] S6. Construct a data null space based on the compressed prediction observation anomaly matrix. Perform variance recovery on the updated set members in the data null space, merge the update results of all local blocks of state variables, and run a composite mapping once after each tempering layer is completed until the posterior set and the posterior prediction observation set are output.

[0067] Specifically, the steps of S1 are as follows: S1.1, Determine the reservoir state variables to be updated. , for The state variables are: permeability, porosity, fault multiplier, capillary pressure parameter, relative permeability parameter, and well control quantity.

[0068] S1.2, a priori sets are generated based on geological understanding, well logging interpretation, seismic inversion results, dynamic understanding, or existing historical fitting results. These priori sets represent the uncertainties in reservoir geological parameters, fluid parameters, or well control quantities. Represented as: ; in, For the first A set of members; This represents the total number of members in the set.

[0069] S1.3, acquire actual observation data and set the standard deviation of observation error. The actual observation data includes oil production, water production, water cut, gas production, water injection volume, water injection rate, bottom hole pressure, pressure recovery interpretation results and seismic attributes.

[0070] S1.4 Observations are grouped according to well number, data type, seismic attribute window or connected seismic region, and local blocks of state variables are divided according to grid block, well control area, flow unit, layer segment or fault block. The purpose of setting local blocks of state variables is to ensure that each local Stein transport update only uses the observation mode related to the local state, thereby reducing the interference of distant irrelevant observations on the local state update.

[0071] S1.5 Establish a composite mapping between the reservoir model and the observation operator. The composite mapping is used to map any reservoir state variable into a predicted observation vector with the same dimension as the actual observation data.

[0072] Specifically, the steps of S2 are as follows: S2.1, Perform a composite mapping on each set member: ; in, For the first A set of members; For the first Predicted observation vectors for each set member; It is a composite mapping function; Obtain the predicted observation set for: ; in, For the first The predicted observation vectors of each set member.

[0073] S2.2, using the mean-removing method and... Scaled projection matrix Construct the state anomaly matrix and the predicted observation anomaly matrix : ; ; in, Describe the direction of change of the current set in the state space. Describes the direction of change of the current set in the observation response space.

[0074] Specifically, the steps of S3 are as follows: S3.1, for the first For each observation group, extract the corresponding prediction observation anomaly submatrix from the prediction observation anomaly matrix. : ; in, For the set of observation indexes; It is the set of real numbers; For the first Number of observations in each observation group; To predict the observation anomaly matrix from a global perspective Selecting row indexes belongs to The submatrix obtained by extracting all rows and retaining all columns has a dimension of . ; Index for observation groups.

[0075] S3.2, using the predicted observation anomaly submatrix Constructing the whitening prediction observation anomaly matrix : ; in, This is the operator for constructing diagonal matrices; For the first Standard deviation of observation error for each observation data; Index for observation data; If the first The observed data remain almost unchanged within the set, i.e. No. If the row norm is very small, the observed data are not identifiable to the current set and can be removed from the effective independent observations estimate or given very small weights.

[0076] S3.3, will Normalize each row to the unit norm, and we get... And calculate the correlation matrix of the observed response. : ; in, For transpose; To predict the whitening observation anomaly matrix The matrix obtained by normalizing each non-zero row according to the L2 norm; The observed response correlation matrix is ​​shrunk, and the shrunk observed response correlation matrix is ​​calculated. : ; in, The shrinkage coefficient; It is an identity matrix.

[0077] S3.4, in an optional implementation, ,right Find the eigenvalues ​​and calculate the posterior influence strength. : ; in, for eigenvalues; Index for source members.

[0078] S3.5, Calculate the number of valid independent observations based on the posterior influence strength. And then cut it: ; ; in, To minimize this value, the effective independent observation number is determined based on the response correlation structure of the predicted observation anomalies and its posterior influence spectrum, without using the magnitude of prior innovation residuals for estimation. This effective independent observation number characterizes the number of independent information directions that the observation data in the current set can provide. For the intensive injection-production dynamic curve of the same well, when the number of time sampling points increases but the number of independent patterns included in the curve shape does not increase accordingly, the effective independent observation number remains within a relatively stable range; when the response differences among observation points within the observation group are significant, the effective independent observation number increases with the increase in the number of identifiable patterns.

[0079] S3.6, determine the redundancy inflation factor based on the number of observations and the number of effective independent observations in the observation group. For the first... Each observation group defines the redundancy expansion factor. for: ; in, To obtain the maximum value.

[0080] S3.7, using the redundancy expansion factor to correct the diagonal observation error variance of the observation data within this observation group: Correct the diagonal observation error variance for the observation data within this observation group: ; in, For the first The observation error variance after correction for each observation data; For the first Redundancy inflation factor of the group to which each observation data belongs; For the first The original error variance of each observation data; Indicates the first The observation group to which each observation data belongs.

[0081] Redundancy correction for diagonal observation error variance is determined based on the ratio between the number of observations and the number of effective independent observations. Let the first... The number of observations in each observation group is The number of effective independent observations is Then, by correcting the diagonal observation error variance of the observation data within the observation group through the redundancy inflation factor, the total information contribution of the observation group to the posterior update is equal to that of the observation group. Matching each independent observation, this process preserves the form of the diagonal observation error matrix, eliminating the need to explicitly construct the complete observation error correlation matrix.

[0082] Corrected diagonal observation error matrix for: ; in, For the first The original error variance of the observation data.

[0083] Step S3 does not aim to accurately estimate the correlation coefficient of each observation error, but rather ensures that the total information content of each redundant time series for the posterior set is approximately [missing information]. An independent observation, rather than Each independent observation.

[0084] Specifically, the steps of S4 are as follows: S4.1, for each local block of state variable, select a set of observation indices related to that local block of state variable based on sample correlation strength, physical distance, well control relationship, flow connectivity relationship, or geological barrier relationship. Select the local block that affects the state variable Observation index set : ; in, The correlation truncation threshold can be set to [value missing]. ; For the local block of state variables and the first The strength of sample correlation between individual observation data The representative variable within the local block of the state variable can be compared with the first... The sample correlation coefficient of an observation data can also be the mean or maximum value of the absolute values ​​of the correlation of variables within a local block of the state variable. To take the absolute value; In practical reservoir applications, it can also be determined based on physical distance, well control relationships, or flow connectivity. .

[0085] S4.2, according to Determine localization weights : .

[0086] S4.3, based on localization weights and the corrected standard deviation of observation error Constructing a local observation whitening weighted matrix ; ; The local observation whitening weighting matrix is ​​used to simultaneously achieve observation error scale normalization, redundancy information reduction, and local correlation weighting. The local observation whitening weighting matrix is ​​essentially a diagonal scaling matrix used to achieve three functions simultaneously: first, utilizing... First, the observation error scale is normalized; second, the repetitive effects of redundant observations are reduced by using the error standard deviation after redundancy correction; third, through... Reduce the impact of observation data that are weakly correlated with local blocks of the current state variable.

[0087] S4.4 uses the local observation whitening weighted matrix to process local predicted observation anomalies, calculates the number of local effective independent observations, and determines the local compressed dimension.

[0088] Specifically, the steps in S4.4 are as follows: S4.4.1, for any local block of state variables First, extract the local blocks related to the state variables. Block-related predictive observation anomaly submatrix : ; in, To predict the observation anomaly matrix from a global perspective Selecting an index set The local submatrix obtained by taking all the corresponding rows and retaining all columns has a dimension of . .

[0089] S4.4.2 Constructing a local whitening prediction observation anomaly matrix using a local observation whitening weighted matrix. : .

[0090] S4.4.3, based on The row correlation matrix and posterior influence spectrum are used to calculate the number of locally effective independent observations. ,according to Take the local compression dimension : ; in, This is the floor symbol.

[0091] S4.4.4, retain the former eigenvectors Define information-constrained observation projection matrix : ; Information-constrained observation projection matrix Compress dense time series into If there are only about 10 valid models out of 60 observation points, then the subsequent assimilation will only be carried out in a residual space of about 10 dimensions, but these 10 models still retain the main fitting information of the entire production curve.

[0092] S4.4.5, for any predicted observation vector The local compression residual is : ; in, This is the actual observation vector composed of all historical observations.

[0093] S4.4.6, using the information-constrained observation projection matrix, the local predicted observation anomaly submatrix is ​​projected onto the independent information pattern space to obtain the locally compressed predicted observation anomaly matrix. : .

[0094] Specifically, the steps of S5 are as follows: S5.1, a data attraction term is constructed based on the local compression residuals and the local compression prediction observation anomaly matrix. The data attraction term is used to reduce the compression residuals. For the first... There are one tempering layer, and the tempering increment is... Local blocks of state variables The Middle Data attraction items of a set member for: ; in, This is a local block state anomaly matrix; For transpose; The data attraction term is the likelihood gradient under the precondition of set covariance, using... Structure, but not formed It also does not use the Kalman gain of ESMDA, so it is not an ESMDA update, but a data score estimate in the Stein stream.

[0095] S5.2 Construct a priori anchoring terms based on the deviations between the current set members and the prior set members. : ; in, The prior anchoring strength is used to constrain the update results from deviating excessively from prior geological knowledge. During the current tempering update phase, the first The set members are in the local block of the state variable. State subvectors on; For the first The set members are in the local block of the state variable. The fixed reference state subvector on, This represents the prior state.

[0096] S5.3 standardizes the state variables within the local block of the state variables, assuming the first... The standard deviation of each state variable in the current set is , define the first Standardized state variables of each set member for: ; in, Local block of state variables for the current tempering phase The set mean vector; This is the index for the global state variable.

[0097] S5.4 Construct the Stein kernel matrix based on the standardized set membership distances. : ; in, It is a natural exponential function; For bandwidth parameters; It is an L2 norm.

[0098] S5.5 Calculate the Stein transport velocity based on the data attraction term, the prior anchoring term, and the kernel repulsion term in the Stein kernel matrix. ; ; in, The metric matrix corresponding to the standardized distance. ; For the first In the first tempering stage The set members are in the state variable block. State subvectors on; The first term of the Stein transport rate pushes the set members in a direction that reduces the compression residual, while the second term is the nuclear repulsion term, which is used to maintain the dispersion of the set members.

[0099] S5.6 Select the update step size based on the linearized prediction residual of the current set, and use the update step size and Stein transport speed to obtain the temporary update set of the local block of state variables.

[0100] During local updates, the information-constrained observation projection matrix projects dense observations related to local blocks of state variables onto a few key information directions, allowing subsequent Stein transport updates to apply data attraction only to the compressed residual space. This process preserves key fitting information while reducing redundant time points from repeatedly compressing the ensemble variance. When a local block of a state variable has a large number of visible observations but few effective independent information patterns, the projection matrix can compress high-dimensional observation residuals into low-dimensional independent information patterns, thereby reducing the ill-conditioned nature of local updates.

[0101] The Stein transport update consists of a data attraction term, a priori anchoring term, and a kernel repulsion term. The data attraction term reduces compression residuals, causing ensemble members to move in a direction consistent with historical data. The priori anchoring term constrains the update results from deviating excessively from prior geological understanding and the original model ensemble. The kernel repulsion term maintains the dispersion among ensemble members, reducing the risk of multiple ensemble members shrinking to the same position. This update process does not require explicit formation of the Kalman gain in ESMDA, nor does it require perturbing observations; therefore, the update mechanism differs from traditional Kalman-type ensemble smoothing updates.

[0102] Specifically, the steps of S6 are as follows: The Stein kernel repulsion term is used to mitigate ensemble shrinkage. When observations are dense and residual attraction is strong, the variance of the temporarily updated ensemble is further compared with the variance of the information constraint target. If the former is lower than a preset ratio, singular value decomposition is performed on the compressed predicted observation anomalies to separate the data activity direction and the approximate data null space direction from the ensemble membership space, and the ensemble anomaly amplitude is recovered along the approximate data null space direction. Therefore, this invention recovers the variance in the approximate data null space of the compressed observation operator.

[0103] S6.1, using local compression to predict the observed anomaly matrix to calculate the information-constrained linear Gaussian target transformation. : ; in, for The identity matrix; Calculate the target variance trace : ; in, Let be the trace of the matrix.

[0104] S6.2, Suppose the temporarily updated state variable local block anomaly is... Its variance trace for: .

[0105] S6.3, if Then for Perform singular value decomposition: ; in, This is the left singular vector orthogonal matrix of the locally compressed predicted observation anomaly matrix; It is an orthogonal matrix of right singular vectors; It is a rectangular diagonal matrix composed of non-negative singular values; This represents the minimum variance retention ratio for dimensionless quantities.

[0106] S6.4, Preserve the right singular vector corresponding to the data activity subspace. Define the null space projection of data in the set member space. And magnify the set anomaly along the data null space projection direction: ; in, It is a vector consisting entirely of 1s; Data null space variance recovery is used to address the problem of local variance insufficiency that may still occur under dense and redundant observations. When the locally updated ensemble variance is lower than the reasonable level given by the linear objective of the information constraint, the anomalous amplitude is recovered along the first-order data null space of the compression observation operator. Since this recovery direction is located in the approximate null space of the compression observation operator, the first-order impact on the compression residual is small when recovering the variance, and it can compensate for ensemble discrepancies without compromising the historical fitting effect as much as possible.

[0107] S6.5 If a state variable belongs to multiple local blocks of state variables, the update results of multiple local blocks are merged according to the following formula: ; in, Weights for merging local blocks of state variables; For state variables, after Stein transport update and data null space variance recovery are completed for local blocks of state variables; To complete the local transport, variance recovery, and weighted merging of overlapping blocks in the current tempered layer, the first... The first set member A global state variable.

[0108] S6.6, after merging all local blocks of state variables, the updated set is obtained. And run a numerical simulation: ; in, For the first The posterior prediction observations of each tempered layer output are combined; For numerical simulation functions; right Repeat steps S2.2 to S6.6, obtaining the updated set in each round. and its posterior predictive observation set After the loop ends, the final posterior set is output. and posterior prediction observation set The total number of orthogonal runs is: ; in, This represents the number of tempering layers.

[0109] To verify the method provided by this invention, three types of redundant observation verification experiments were conducted, including: linear dense redundant observation experiments, two-dimensional moderate nonlinear redundant observation experiments, and high-dimensional surrogate experiments.

[0110] Linear Dense Redundant Observation Experiment: This experiment uses a 50-dimensional one-dimensional Gaussian random field and sets up 60 densely time-sampled linearly smoothed observations. The actual observation errors are time-dependent, but only diagonal observation errors are used in the assimilation calculation. Validation results show that the method of this invention yields an average effective independent observation count of 10.05 and an average redundancy inflation factor of 5.973, indicating that the 60 dense observation points are identified as approximately 10 effective independent information patterns. The posterior quality distance of the method of this invention is 0.08225, which is better than the 0.09002 of the ensemble smoothing method and the 0.1255 of the localized ensemble smoothing method. The posterior variance ratio of the method of this invention is 0.8819, which is higher than the 0.7851 of the ensemble smoothing method, indicating that the method of this invention can reduce excessive shrinkage of the ensemble variance while maintaining the fitting ability.

[0111] Two-Dimensional Moderate Nonlinear Redundant Observation Experiment: This experiment uses a two-dimensional analytical moderately monotonic nonlinear observation function with 50 redundant time points. Validation results show that the posterior quality distance of the proposed method is 0.04513, a reduction of approximately 70.73% compared to the optimal baseline method's 0.1542; the prediction normalization error of the proposed method is 0.1774, superior to the ensemble smoothing method's 0.3278 and the localized ensemble smoothing method's 1.2405; the posterior variance ratio of the proposed method is 0.8809, significantly higher than the ensemble smoothing method's 0.2788. These results demonstrate that, under conditions of moderate nonlinearity and temporal redundancy in observations, the proposed method can simultaneously improve posterior quality, prediction accuracy, and variance preservation.

[0112] High-dimensional surrogate experiment: This experiment uses a 100-dimensional high-dimensional field, including 14 effective principal components, and sets up 3 groups of dense well response observations, totaling 90 observation points. The observed responses are moderately monotonically nonlinear. Validation results show that the posterior quality distance of the proposed method is 0.1497, which is better than the 0.2112 of the ensemble smoothing method and the 0.1705 of the localized ensemble smoothing method; the posterior variance ratio of the proposed method is 0.9075, which is higher than the 0.8118 of the ensemble smoothing method and is closer to a reasonable posterior variance level. These results demonstrate that even with a high-dimensional state space and multiple groups of dense well responses, the proposed method can still effectively reduce the impact of redundant observation counts and maintain posterior ensemble reliability.

[0113] like Figure 2 As shown, in the three types of redundant observation verification experiments, the posterior quality distance of the method of the present invention is lower than the best results of the set smoothing method and the localized set smoothing method, indicating that the posterior set obtained by the present invention is closer to the reference posterior as a whole.

[0114] like Figure 3 As shown, in the three types of redundant observation verification experiments, the posterior variance ratio of the method of the present invention remained around 0.88 to 0.91, indicating that it can maintain the reliability of the posterior set variance well.

[0115] like Figure 4 As shown, in the two-dimensional moderately nonlinear redundant observation experiment, the method of the present invention significantly reduces the prediction error; in the linear dense redundant observation experiment and the high-dimensional surrogate experiment, the method of the present invention maintains a prediction error similar to that of the baseline method.

[0116] like Figure 5 As shown, the original observation mismatch of the method of the present invention is not always the lowest. This is because the ensemble smoothing method may overfit dense observations with correlated errors, while the present invention emphasizes avoiding ensemble collapse while fitting the observations.

[0117] like Figure 6 As shown, the method of the present invention has a significant improvement over the over-shrinking set smoothing method in terms of 90% coverage, indicating that its posterior uncertainty expression is more reliable.

[0118] Figures 2 to 6 The verification results demonstrate that this invention, through a combined mechanism of effective independent observation estimation, local observation whitening weighting, information-constrained observation compression, Stein kernel exclusion, and data null space variance recovery, obtains more reliable posterior sets and more robust history fitting results under the same forward modeling budget. To verify the method provided by this invention, reservoir numerical simulation experiments were conducted.

[0119] This experiment applies the method of this invention to automatic history fitting of two-dimensional heterogeneous reservoirs. The ECLIPSE reservoir numerical simulator is called to complete the oil-water two-phase flow calculation for each set member. This experiment uses the actual water injection rate of four water injection wells and the oil production rate of nine production wells as fitting observations. The experiment focuses on explaining the complete engineering process of reservoir model establishment, prior set parameterization, injection-production dynamic grouping, information constraint transport update, set dispersion diagnosis, and result output.

[0120] A basic reservoir model was established. The reservoir model grid has a dimension of 60×60×1, with a total of 3600 effective grids; the plane size of a single grid is 20m×20m, the thickness is 10m, the porosity is fixed at 0.2, and the state variable to be updated is the horizontal permeability of each grid. Four water injection wells are located at grid positions (17,17), (17,43), (43,17), and (43,43), respectively; nine production wells are located at (4,4), (4,30), (4,56), (30,4), (30,30), (30,56), (56,4), (56,30), and (56,56), respectively. The water injection wells adopt a target water injection rate of 1500 cubic meters / day and a bottom hole pressure upper limit of 410 bar; when the pressure upper limit constraint is activated, the actual water injection rate varies with the reservoir permeability. The production wells are controlled with a bottom hole pressure of 310 bar. The historical period was divided into 50 time steps, each lasting 30 days, for a total historical period of 1500 days. This experiment selected the actual water injection rate of four injection wells and the oil production rate of nine production wells as the subjects for joint historical data fitting.

[0121] A priori set is constructed and parameterized. 1000 candidate permeability fields with a preset spatial correlation structure are generated, ranging from 30 millidarcy to 3500 millidarcy. The first permeability field is selected as the reference true model and removed from the prior training samples. A principal component analysis model is established using the remaining 999 permeability fields, reducing the 3600-dimensional grid permeability to 100-dimensional principal component coefficients. During history fitting, 100 initial set members following a standard normal distribution are generated in the 100-dimensional principal component space, and 100 initial permeability fields are obtained through inverse principal component transformation. Thus, all 100-dimensional principal component coefficients are used as a user-specified state variable block; in other implementations, multiple local state variable blocks can be further divided according to well control zones, flow units, or layers.

[0122] Observational data and observation groups were constructed. A reservoir numerical simulator was run once on the reference real model, simultaneously extracting the actual water injection rates of four injection wells and the oil production rates of nine production wells. Each well included 50 consecutive time points, resulting in 13 injection-production dynamic time series, totaling 650 observation points. Thirteen observation groups were formed according to well number and response type, each group containing 50 consecutive time points; the first four groups represented the actual water injection rates of the injection wells, and the latter nine groups represented the oil production rates of the production wells. The standard deviation of the observation error was set to 10% of the absolute value of the corresponding reference response, with a lower limit of 0.5 cubic meters per day for absolute error. The corresponding Gaussian error was superimposed on the reference response to obtain the historical observation vector. The diagonal observation error matrix form was maintained during assimilation calculations; a complete 650×650 observation error correlation matrix was not constructed.

[0123] Set the ICST-RST parameters and complete the initial forward modeling. Number of set members. Set to 100, number of tempering layers Set to 4, the coefficient for preserving the variance of the null space. The values ​​were set to 0.95, the Stein kernel repulsion coefficient to 0.04, the redundancy observation variance inflation index to 1.0, and the minimum standard deviation retention ratio for each state set to 0.70. First, the reservoir numerical simulator was run once on the 100 initial set members to obtain the initial predicted observation set, and a 100-dimensional state anomaly matrix and a 650-dimensional predicted observation anomaly matrix were constructed.

[0124] The number of effective independent observations was estimated group by group, and an information-constrained projection was constructed. For each of the 13 observation groups, a predicted observation anomaly submatrix was extracted. After observation error whitening, row normalization, and eigenvalue decomposition of the response correlation matrix, the posterior influence spectrum was calculated. The average number of effective independent observations per group in the first to fourth update layers were 9.4244, 10.6356, 12.4667, and 12.9756, respectively, with corresponding average redundancy inflation factors of 5.8657, 4.8842, 4.0799, and 3.9631. The average number of effective independent observations across the four layers was 11.3756, indicating that the 50 dense time points in each group mainly represent approximately 9 to 13 identifiable information directions in the current set. The total rank of the information-constrained projection in each update layer was 118, 140, 167, and 172, respectively, thus compressing the original 650-dimensional observation residuals into the main information pattern space.

[0125] Stein transport and data null space variance recovery are performed. Within each updated layer, the near-end data attraction direction of the ensemble subspace is solved based on the compressed observation residuals and compressed predicted observation anomalies, and a radial basis function kernel repulsion term is superimposed. The variance recovery objective uses the likelihood score consistent with the current layer's tempering coefficient, and distinguishes between significant and weak information directions based on the posterior influence of each singular direction; when the variance of the weak information direction is lower than a preset proportion, the anomaly amplitude is recovered along the approximate data null space. When dense observations make the ensemble member space numerically full-rank and the updated total standard deviation is lower than 0.70 of the previous layer, minimum dispersion protection is activated, and the total standard deviation is restored to this proportion. Subsequently, the 100-dimensional principal component coefficients are inversely transformed into a 3600-dimensional permeability field, and the next batch of reservoir numerical simulations is run on 100 updated members.

[0126] The forward modeling budget is controlled and posterior results are output. History fitting includes one initial ensemble forward modeling batch and four updated forward modeling batches, each batch containing 100 ensemble members. Therefore, the reservoir numerical simulation calls during the assimilation phase are 500, the same as the classic ESMDA forward modeling budget under the same ensemble size and four-layer tempering setting. A reference real model is run again to uniformly generate historical observations for 13 wells; therefore, this experiment actually performs 501 reservoir numerical simulations. After completing the fourth updated layer, the posterior principal component set, posterior permeability set, posterior actual water injection rate set for the four water injection wells, posterior oil production rate set for the nine production wells, the number of effective independent observations for each layer, the redundancy expansion factor, the ensemble dispersion diagnostic quantity, and the objective function are output.

[0127] Evaluation of the dynamic fitting results of injection and sampling: All ensemble members were evaluated using a normalized objective function including data mismatch and principal component regularization terms. The mean of the initial batch objective function was 125.6584; the mean objective functions after the first to fourth update layers were 15.8078, 2.6597, 1.4026, and 1.1267 respectively, with a median of 1.1268 for the final objective function. The mean of pure data mismatch, standardized by observation error, decreased from 125.5813 to 1.0482, a reduction of 99.1653%. The root mean square error of the initial ensemble mean response relative to historical observations was 207.1377 cubic meters per day, while the root mean square error of the final ensemble mean response was 50.9868 cubic meters per day, a reduction of 75.3851%. The root mean square error of the posterior mean permeability field relative to the reference true field decreased from 412.1001 millidarcy for the prior mean field to 306.5555 millidarcy, a decrease of 25.6114%.

[0128] In this experiment, the ensemble size, number of tempered layers, observation error ratio, variance preservation coefficient, kernel exclusion coefficient, redundant observation variance inflation index, and minimum standard deviation preservation ratio for each layer were all optional settings. In other reservoirs, these parameters can be adjusted based on the ensemble size, injection-production dynamic observation density, forward modeling cost, and posterior variance diagnosis results. When the posterior ensemble dispersion is lower than the engineering requirements, the minimum standard deviation preservation ratio or kernel exclusion coefficient for each layer can be increased, and the extrapolation capability of the posterior model can be tested using retained prediction period data.

[0129] like Figure 7 As shown, the fitting results of the actual water injection rate of injection well 2 over 50 historical time steps are displayed separately. The gray curve cluster represents the response of the 100 initial ensemble members, the blue curve cluster represents the response of the 100 posterior ensemble members after four ICST-RST updates, and the red discrete observation points represent noisy historical observations. No lines are drawn between the observation points. The blue filled area is not drawn in the figure, nor is a separate ensemble mean line drawn; the visually dark blue band is formed by the superposition of multiple semi-transparent posterior member curves. The root mean square error of the ensemble mean response of injection well 2 relative to the historical observations decreased from 370.2777 cubic meters / day in the initial ensemble to 20.9055 cubic meters / day in the posterior ensemble, a reduction of 94.3541%, indicating that the change level of the actual water injection rate of this well has been effectively corrected.

[0130] like Figure 8 As shown, the oil production rate fitting results of production well 8 over 50 historical time steps are displayed separately, with plotting rules consistent with... Figure 7 Consistent. Historical observations of production well 8 initially increased from approximately 100 cubic meters per day to approximately 250 cubic meters per day, then gradually decreased to approximately 20 cubic meters per day at the end of the historical period; the posterior ensemble can reproduce the trend of initial increase followed by decrease and the intermediate decline process. The root mean square error of the mean response of this well relative to historical observations decreased from 98.7695 cubic meters per day to 16.7692 cubic meters per day, a reduction of 83.0219%. Figure 7 and Figure 8Only injection well 2 and production well 8 were selected as representative wells for display, without changing the joint assimilation object and update process. The actual water injection rate of all 4 injection wells and the oil production rate of all 9 production wells were still included in the historical fitting, with each well including 50 time points, forming a total of 13 observation groups and 650 observation points. The mean standard deviation of the 100-dimensional principal component set decreased from 0.9968 to 0.3452, with a retention rate of 34.63%; the effective rank of the set decreased from 48.7116 to 22.0774, with a retention rate of 45.32%; the grid mean standard deviation of the final posterior permeability set was 59.8156 millidarcy, and the mean standard deviation of the final well dynamic response was 7.2931 cubic meters / day, which is 5.34% of the initial value of 136.5682 cubic meters / day. The above results indicate that the joint observation of the 13 wells resulted in a significant concentration of the posterior response set, but the trajectories of the 100 posterior members were still distinguishable, and the effective rank of the set did not degenerate to 1. In engineering applications, the effective rank, set standard deviation, and data from the retained prediction period should be combined to further evaluate the posterior uncertainty.

[0131] like Figure 9 As shown, the horizontal axis sequentially arranges all 500 set members of the initial batch and the four updated batches, while the vertical axis represents the normalized objective function on a logarithmic scale. The objective function point cloud gradually decreases from the highly discrete and mismatched state of the initial batch, with the mean objective function value of each batch decreasing from 125.6584 to 15.8078, 2.6597, 1.4026, and 1.1267 respectively. Figure 7 and Figure 8 As shown by the 100 blue posterior member curves, the 34.63% principal component standard deviation retention rate, and the 45.32% effective rank retention rate, this invention can complete the information compression and history fitting of dense injection and acquisition dynamics while maintaining the predetermined assimilation forward modeling budget, and retain the multi-member posterior representation through non-zero posterior effective rank retention. For significant ensemble contraction caused by joint observations, the minimum standard deviation retention rate or nuclear exclusion coefficient can be further adjusted according to the posterior dispersion diagnosis results.

[0132] 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 reservoir history fitting method based on information-constrained Stein transport, characterized in that, Specifically, the steps include the following: S1, obtain the prior set of reservoir state variables to be updated, actual observation data, standard deviation of observation error, observation grouping, local blocks of state variables, and composite mapping of reservoir model and observation operator; S2, perform a composite mapping on the prior set to obtain the predicted observation set, and construct the state anomaly matrix and the predicted observation anomaly matrix based on the prior set and the predicted observation set, respectively; S3, based on the predicted observation anomaly matrix, estimate the effective independent observation number of each observation group according to the observation group, and perform redundancy correction on the diagonal observation error variance of the observation according to the effective independent observation number. The effective independent observation number is determined according to the response correlation structure and posterior influence spectrum of the predicted observation anomaly, and the redundancy correction of the diagonal observation error variance is determined according to the ratio between the number of observations and the effective independent observation number. S4. For each local block of state variable, select observations related to the local block of state variable, construct a local observation whitening weighting matrix, and construct an information-constrained observation projection matrix based on the local observation whitening weighting matrix. The local observation whitening weighting matrix is ​​used to process local predicted observation anomalies, calculate the number of local effective independent observations and determine the local compression dimension. The information-constrained observation projection matrix is ​​used to compress the dense time series into the local compression dimension of independent information directions, and compress the local observation residual matrix and the local predicted observation anomaly matrix into the independent information pattern space. S5, construct the Stein transport update direction in the independent information pattern space, and update the set members of the local block of state variables according to the Stein transport update direction; S6. Construct a data null space based on the compressed prediction observation anomaly matrix. Perform variance recovery on the updated set members in the data null space, merge the update results of all local blocks of state variables, and run a composite mapping once after each tempering layer is completed until the posterior set and the posterior prediction observation set are output.

2. The reservoir history fitting method based on information-constrained Stein transport according to claim 1, characterized in that, The specific steps of S1 are as follows: S1.1, Determine the reservoir state variables to be updated. , for The state variables are: permeability, porosity, fault multiplier, capillary pressure parameter, relative permeability parameter, and well control quantity. S1.2, Generate a priori set based on geological understanding, well logging interpretation, seismic inversion results, dynamic understanding, or existing historical fitting results. Represented as: ； in, For the first A set of members; The total number of members in the set; S1.3, Obtain actual observation data and set the standard deviation of observation error. The actual observation data includes: oil production, water production, water cut, gas production, water injection volume, water injection rate, bottom hole pressure, pressure recovery interpretation results, and seismic attributes. S1.4, divide the observation groups according to well number, data type, seismic attribute window or connected seismic region, and divide the local blocks of state variables according to grid block, well control area, flow unit, layer segment or fault block; S1.5, establish a composite mapping between the reservoir model and the observation operator.

3. The reservoir history fitting method based on information-constrained Stein transport according to claim 1, characterized in that, The specific steps of S2 are as follows: S2.1, Perform a composite mapping on each set member: ； in, For the first A set of members; For the first Predicted observation vectors for each set member; It is a composite mapping function; Obtain the predicted observation set for: ； in, For the first Predicted observation vectors for each set member; S2.2, using the mean-removing method and... Scaled projection matrix Construct the state anomaly matrix and the predicted observation anomaly matrix : ； 。 4. The reservoir history fitting method based on information-constrained Stein transport according to claim 1, characterized in that, The specific steps of S3 include: S3.1, for the first For each observation group, extract the corresponding prediction observation anomaly submatrix from the prediction observation anomaly matrix. : ； in, For the set of observation indexes; It is the set of real numbers; For the first Number of observations in each observation group; To predict the observation anomaly matrix from a global perspective Selecting row indexes belongs to The submatrix obtained by extracting all rows and retaining all columns has a dimension of . ; Index for observation groups; S3.2, using the predicted observation anomaly submatrix Constructing the whitening prediction observation anomaly matrix : ； in, This is the operator for constructing diagonal matrices; For the first Standard deviation of observation error for each observation data; Index for observation data; S3.3, will Normalize each row to the unit norm, and we get... And calculate the correlation matrix of the observed response. : ； in, For transpose; To predict the whitening observation anomaly matrix The matrix obtained by normalizing each non-zero row according to the L2 norm; The observed response correlation matrix is ​​shrunk, and the shrunk observed response correlation matrix is ​​calculated. : ； in, The shrinkage coefficient; It is the identity matrix; S3.4, for Find the eigenvalues ​​and calculate the posterior influence strength. : ； in, for eigenvalues; Index for source members.

5. The reservoir history fitting method based on information-constrained Stein transport according to claim 4, characterized in that, S3 also includes the following steps: S3.5, Calculate the number of valid independent observations based on the posterior influence strength. And then cut it: ； ； in, To obtain the minimum value; S3.6, for the first Each observation group defines the redundancy expansion factor. for: ； in, To obtain the maximum value; S3.7, Correct the diagonal observation error variance for the observation data within this observation group: ； in, For the first The observation error variance after correction for each observation data; For the first Redundancy inflation factor of the group to which each observation data belongs; For the first The original error variance of each observation data; Indicates the first The observation group to which each observation data belongs; Corrected diagonal observation error matrix for: ； in, For the first The original error variance of the observation data.

6. The reservoir history fitting method based on information-constrained Stein transport according to claim 1, characterized in that, The specific steps of S4 are as follows: S4.1, for each local block of state variables, select a set of observation indices related to the local block of state variables based on sample correlation strength, physical distance, well control relationship, flow connectivity relationship, or geological barrier relationship. Select the local block that affects the state variable Observation index set : ； in, The correlation truncation threshold; For the local block of state variables and the first The strength of the sample correlation between individual observation data; To take the absolute value; S4.2, according to Determine localization weights : ； S4.3, based on localization weights and the corrected standard deviation of observation error Constructing a local observation whitening weighted matrix : ； S4.4 uses the local observation whitening weighted matrix to process local predicted observation anomalies, calculates the number of local effective independent observations, and determines the local compressed dimension.

7. The reservoir history fitting method based on information-constrained Stein transport according to claim 6, characterized in that, The specific steps of S4.4 are as follows: S4.4.1, for any local block of state variables First, extract the local blocks related to the state variables. Related predictive observation anomaly submatrix : ； in, To predict the observation anomaly matrix from a global perspective Selecting an index set The local submatrix obtained by taking all the corresponding rows and retaining all columns has a dimension of . ; S4.4.2 Constructing a local whitening prediction observation anomaly matrix using a local observation whitening weighted matrix. : ； S4.4.3, based on The row correlation matrix and posterior influence spectrum are used to calculate the number of locally effective independent observations. ,according to Take the local compression dimension : ； in, The rounding up symbol; S4.4.4, retain the former eigenvectors Define information-constrained observation projection matrix : ； S4.4.5, for any predicted observation vector The local compression residual is : ； in, The actual observation vector is composed of all historical observations; S4.4.6, using the information-constrained observation projection matrix, the local predicted observation anomaly submatrix is ​​projected onto the independent information pattern space to obtain the locally compressed predicted observation anomaly matrix. : 。 8. The reservoir history fitting method based on information-constrained Stein transport according to claim 1, characterized in that, The specific steps of S5 include: S5.1, construct a data attraction term based on the local compression residuals and the local compression predicted observation anomaly matrix, for the There are one tempering layer, and the tempering increment is... Local blocks of state variables The Middle Data attraction items of a set member for: ； in, This is a local block state anomaly matrix; For transpose; S5.2 Construct a priori anchoring terms based on the deviations between the current set members and the prior set members. : ； in, For prior anchoring strength; During the current tempering update phase, the first The set members are in the local block of the state variable. State subvectors on; For the first The set members are in the local block of the state variable. The fixed reference state subvector on, Indicates the prior state; S5.3 standardizes the state variables within the local block of the state variables, assuming the first... The standard deviation of each global state variable in the current set is , define the first Standardized state variables of each set member for: ； in, Local block of state variables for the current tempering phase The set mean vector; Index for global state variables; S5.4 Construct the Stein kernel matrix based on the standardized set membership distances. : ； in, It is a natural exponential function; For bandwidth parameters; It is an L2 norm; S5.5 Calculate the Stein transport velocity based on the data attraction term, the prior anchoring term, and the kernel repulsion term in the Stein kernel matrix. ; ； in, The metric matrix corresponding to the standardized distance. ; For the first In the first tempering stage The set members are in the local block of the state variable. State subvectors on; S5.6 Select the update step size based on the linearized prediction residual of the current set, and use the update step size and Stein transport speed to obtain the temporary update set of the local block of state variables.

9. The reservoir history fitting method based on information-constrained Stein transport according to claim 1, characterized in that, The specific steps of S6 include: S6.1, using local compression to predict the observed anomaly matrix to calculate the information-constrained linear Gaussian target transformation. : ； in, for The identity matrix; Calculate the target variance trace : ； in, The trace of the matrix; S6.2, Suppose the temporarily updated state variable local block anomaly is... Its variance trace for: ； S6.3, if Then for Perform singular value decomposition: ； in, This is the left singular vector orthogonal matrix of the locally compressed predicted observation anomaly matrix; It is an orthogonal matrix of right singular vectors; It is a rectangular diagonal matrix composed of non-negative singular values; This represents the minimum variance retention ratio for dimensionless quantities.

10. The reservoir history fitting method based on information-constrained Stein transport according to claim 9, characterized in that, S6 also includes the following steps: S6.4, Preserve the right singular vector corresponding to the data activity subspace. Define the null space projection of data in the set member space. And magnify the set anomaly along the data null space projection direction: ； in, It is a vector consisting entirely of 1s; S6.5 If a state variable belongs to multiple state variable local blocks, the update results of multiple state variable local blocks are merged according to the following formula: ； in, Weights for merging local blocks of state variables; For state variables, after Stein transport update and data null space variance recovery are completed for local blocks of state variables; To complete the local transport, variance recovery, and weighted merging of overlapping blocks in the current tempered layer, the first... The first set member One global state variable; S6.6, after merging all local blocks of state variables, the updated set is obtained. And run a numerical simulation: ； in, For the first The posterior prediction observations of each tempered layer output are combined; For numerical simulation functions; right Repeat steps S2.2 to S6.6, obtaining an updated posterior set in each round. and its posterior predictive observation set After the loop ends, the final posterior set is output. and posterior prediction observation set The total number of orthogonal runs is: ； in, This represents the number of tempering layers.

Images