Space-time correlation adaptive MCMC wave impedance inversion method

Abstract

The invention relates to the technical field of seismic exploration, and relates to an oil and gas seismic impedance inversion method, in particular to a space-time correlation adaptive MCMC wave impedance inversion method. According to the method, the problems of instability, low calculation efficiency, high multiplicity of solutions and insufficient uncertainty quantization of a traditional random inversion method in the prior art are solved. The method comprises the following steps: constructing a space-time correlation covariance matrix; constructing a target function of the fitting degree; performing iterative inversion; and post-processing an inversion result. According to the method, an inversion strategy based on a Bayesian framework is adopted, space-time correlation covariance is constructed as a prior constraint, spatial correlation and stratum vertical continuity describing impedance parameters are displayed, a physical model is combined to ensure that a prediction result meets the matching performance of observed seismic data, meanwhile, deep fusion with a delay rejection self-adaptive MCMC algorithm is carried out, and the prediction accuracy is improved. And collaborative innovation of physical model constraint and statistical optimization in seismic wave impedance inversion is realized. And the requirement of high-precision seismic inversion is met.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration, and relates to an oil and gas seismic impedance inversion method, in particular to a spatio-temporally correlated adaptive MCMC wave impedance inversion method. Background Art

[0002] In recent years, due to the fact that the focus of oil and gas exploration and development at home and abroad has gradually shifted to deep formations, the identification, tracking, and quantitative evaluation of reservoirs have become an urgent need for oilfield development. In the process of actual oil and gas exploration, seismic inversion technology has been increasingly widely used in seismic reservoir prediction. However, in the face of increasingly complex geological targets, traditional geophysical inversion methods are difficult to meet the accuracy requirements of modern exploration.

[0003] At present, the equation-solving algorithms for seismic inversion are mainly divided into two categories. The first major category is the linear and quasi-linear solution based on gradient calculation. The second major category is the non-linear algorithm gradually evolved based on the biological evolution theory, also known as the stochastic solution algorithm. The limitations of non-linearity are as follows: First, it requires multiple iterative calculations of the forward and inverse processes, and the computational workload is huge. When dealing with large-scale seismic data, it faces limitations in computational efficiency and storage space. Second, the convergence of non-linear inversion is difficult to guarantee. Especially under complex geological conditions, the inversion process may not converge or the convergence speed is very slow, and the inversion result has extremely strong oscillation characteristics and is difficult to show the true result. In the development process of the field of stochastic seismic parameter inversion, the inversion method based on Bayesian theory once occupied an important position. This method relies on Bayesian theory to construct a seismic inversion objective function, transforms seismic parameter information into a conditional probability distribution, and solves the objective function through various optimization algorithms, and finally realizes the reliability evaluation of seismic parameters. However, the application of Bayesian theory in seismic inversion was once limited by the complexity of the denominator integral calculation. However, with the rapid progress of technology, especially the improvement of computer technology and the introduction of the Markov chain Monte Carlo (MCMC) method, the originally complex numerical solution process of the denominator integral has become more efficient and convenient. The MCMC method effectively avoids the complex marginal integral calculation in the Bayesian formula through its unique algorithm design, and significantly promotes the practical application of Bayesian theory. The core of this method is to construct a Markov chain whose stationary distribution is consistent with the target posterior distribution, and through repeated iteration until a stationary state is reached, so as to obtain samples of the posterior distribution.

[0004] However, the integration of well logging information and the high- and low-frequency components provided by the model in the inversion results introduces the problem of non-uniqueness in inversion. At the same time, its effect is also restricted by the number of wells and the well pattern distribution. To reduce the non-uniqueness of inversion, a prior information based on geological knowledge is usually introduced during the inversion process. However, the existing method for imposing prior information constraints directly calculates the error with the inversion results, resulting in a decrease in the accuracy of the inversion results.

[0005] In summary, the current research on seismic stochastic impedance inversion methods has the following problems:

[0006] 1. Difficult to converge: The convergence is difficult to guarantee. Especially under complex geological conditions, the inversion process may not converge or the convergence speed is very slow.

[0007] 2. The problem of non-uniqueness is prominent: Due to the non-uniqueness of the inversion problem, the results may have multiple interpretations, and the inversion results have strong oscillation characteristics, making it difficult to show the true results and increasing the difficulty of geological interpretation.

[0008] 3. Low computational efficiency: Multiple iterations of forward and inversion processes are required, and the computational amount is huge. When dealing with large-scale seismic data, it faces limitations in computational efficiency and storage space.

[0009] 4. High dependence on prior information: Stochastic inversion usually needs to rely on prior information to constrain the sampling range. However, if the prior information is inaccurate or incomplete, the inversion results may deviate.

[0010] There is an urgent need for a new type of seismic stochastic impedance inversion method that can solve the above problems. Summary of the Invention

[0011] The present invention proposes a spatio-temporally correlated adaptive MCMC wave impedance inversion method, which solves the problems of instability, low computational efficiency, strong non-uniqueness, and insufficient uncertainty quantification existing in the traditional stochastic inversion methods in the prior art.

[0012] The technical solution of the present invention is realized as follows: A spatio-temporal related adaptive MCMC wave impedance inversion method, comprising the following steps: S1: Construct a spatio-temporal related covariance matrix: Taking the vertical continuity of the geological sequence as the core prior, using the Toeplitz structure combined with the squared exponential kernel function, and introducing a regularization term at the same time; S2: Construct an objective function for the fitting degree between the impedance model and the observed data; Combining physical models and statistical priors through multiple constraint terms; S3: Iterative inversion: By repeating the iteration and controlling the maximum number of iterations through the inversion residual, an optimal impedance inversion structure is obtained; The iterative content includes: An adaptive mechanism, specifically: Through the exponentially weighted strategy moving average dynamics, fusing historical sampling information, updating the covariance matrix every certain number of iterations, so that the proposal distribution gradually approaches the actual geometric form of the posterior parameter space; S4: Post-processing of the inversion result: Calculate the posterior standard deviation, quantify the uncertainty of the inversion result, and provide a risk reference for exploration decisions; At the same time, check the convergence of the Markov chain through autocorrelation analysis to ensure the reliability of the inversion result.

[0013] A further technical solution, the sampling process of the adaptive mechanism in step S3 is specifically to improve the sampling efficiency through the rejection jump mechanism DRAM by two-stage dynamic adjustment: S31: In the first stage, generate global exploratory samples based on the current covariance and the M-H criterion; S32: When the proposal is rejected, in the second stage, shrink the current covariance to generate a local refined search step size, generate sub-optimal proposals within the local neighborhood, and combine physical constraints to achieve refined search with controllable exploration risks.

[0014] A preferred technical solution, the adaptive mechanism in step S3 is specifically: Considering that the wave impedance is blocky, use adaptive window smoothing filtering to constrain the seismic wave impedance inversion result.

[0015] A further technical solution, the construction of the covariance in step S1 is specifically:

[0016] S11: Data preparation and processing: Construct a real acoustic impedance model containing 350 layers to simulate the layered structure of the underground formation; According to the acoustic impedance model, calculate the reflection coefficient at the adjacent formation interfaces;

[0017]

[0018] Where x represents the well logging impedance parameter and r represents the reflection coefficient; Through the convolution of the seismic wavelet and the reflection coefficient:

[0019] s = w × r #(2)

[0020] Where \(w\) represents the Ricker wavelet, \(\times\) represents the convolution operator, and \(s\) represents the synthetic seismic record; it is used to generate the synthetic seismic record; simulate the actual seismic data and add noise to be closer to the real observation;

[0021] s noise = s+\(\epsilon\), \(N\sim(0,\sigma\) 2 I)#(3)

[0022] Where \(s\) noise is the seismic record after adding noise, \(\sigma\) 2 is used to control the noise intensity, and \(\epsilon\) is the noise term;

[0023] S12: Construction of the initial covariance: According to the correlation hypothesis, construct the initial spatio-temporal correlation covariance matrix through the space of the formation, and define the covariance decay relationship between parameters based on the Gaussian kernel function;

[0024]

[0025] Where \(d\) represents the layer spacing, \(l\) represents the correlation length, \(\delta\) 2 represents the variance, and \(C\) is the one-dimensional covariance column vector constructed at this time; then through the Toeplitz matrix transformation, expand the one-dimensional covariance vector into a symmetric positive definite covariance matrix, so that the diagonal elements represent the variance of the parameters themselves, and the non-diagonal elements decay exponentially with the increase of the layer spacing; finally, add a regularization term to ensure the numerical stability of the matrix for subsequent Gaussian kernel decomposition;

[0026]

[0027] Convert the vector \(C\) into a Toeplitz matrix, generate a symmetric covariance matrix, and add a regularization term as follows:

[0028] \(C_0 = Toeplitz(C)+\epsilon I\)#(6)

[0029] Where \(C_0\) is the initial covariance, \(i\) and \(j\) are the row and column indices of the matrix \(C\), \(\epsilon\) is the regularization coefficient, and \(I\) is the identity matrix of the same type as \(C_0\).

[0030] For a further technical solution, the construction of the objective function in step S2 is specifically: The objective function is the core in Bayesian inversion, which is used to quantify the matching degree between the impedance in the model and the observed data, and at the same time incorporate prior constraints; First, construct a difference matrix:

[0031]

[0032] After that, an objective function for the fitting degree between the impedance model and the observed data is constructed, and the sparsity of the reflection coefficient is promoted by the L1 norm, implicitly assuming geological conditions; and the L2 norm constrains the gradient of the impedance model, making the inversion result smooth, suppressing high-frequency noise and pulling the inversion result towards the initial guess model to avoid the model deviating from the reasonable range due to data noise;

[0033]

[0034] where θ is the current impedance value, β i＝1 , 2 is the coefficient of the constraint term, and ss(θ) is the objective function at the current impedance value.

[0035] A further technical solution is that the specific content of S3 is:

[0036] S31: The first proposal:

[0037] θ t+1 = θ t + τ, τ ~ N(0, C t )#(9)

[0038] where θ t is the current model, C t is the current covariance matrix, and θ t+1 is the candidate model generated; Calculate the acceptance probability according to the M-H criterion:

[0039]

[0040] Since the proposal distribution is a uniform distribution, it is symmetric. At the same time, based on Bayesian theory, the above formula is further transformed into:

[0041]

[0042] where α1 is the acceptance probability, is the scaling factor;

[0043] S32: If the above candidate value is not accepted, enter the second proposal:

[0044] θ′ t+1 = θ t + τ′, τ′ ~ N(0, γ m C t )#(12)

[0045] where θ t is the current model, C t is the current covariance matrix, γ m is the shrinkage factor of the covariance, and θ′ t+1 is the candidate model generated; Similarly, we get:

[0046]

[0047] At this time, the second candidate value is selected with probability α2. If it is accepted, then θ′ t+1 = θ t , and if it is not accepted, θ is continued to be maintained t and the next iteration is carried out.

[0048] For a further technical solution, specifically in the step S3, the process of adaptively adjusting the covariance is as follows:

[0049] When the number of iterations exceeds a preset threshold, at fixed intervals, impedance samples of a recent window are extracted from the historical sampling results, and its empirical covariance matrix is calculated;

[0050]

[0051] Through the above formula transformation, the covariance formula for the (j + 1)-th time can be obtained:

[0052]

[0053] When updating, the exponential weighted moving average method is adopted to perform weighted fusion on the historical covariance matrix and the current window covariance; among them, the historical weight coefficient controls the retention ratio of old information and gives a higher weight to long-term statistical characteristics; the current window covariance reflects the recent sampling trend and enhances the adaptability to non-stationary posterior distributions; to avoid numerical problems caused by insufficient samples or high-dimensional singularity of the covariance matrix, a small amount of diagonal regularization term is added to it after each update to force the positive definiteness of the matrix; at the same time, the matrix ill-condition is detected through eigenvalue decomposition, and the diagonal elements are further adjusted if necessary to ensure the physical validity of the proposed distribution; considering that the covariance may mutate and the result may be greatly affected, the exponential weighted moving average is introduced:

[0054]

[0055] In the above formulas (14)(15)(16), C j represents the covariance matrix updated at the j-th time, t0 is the initial threshold, s d represents the scaling factor, ε represents the regularization parameter, represents the update weight, I d is the identity matrix of the same type as C0; through repeated iteration and controlling the maximum number of iterations through the inversion residual, the optimal impedance inversion structure is obtained; then, the seismic wave impedance inversion result is constrained by the adaptive window smoothing filter (formula 17) to obtain a more accurate inversion result; among them, the formula of EPS is as follows:

[0056]

[0057] where x1 is the predicted impedance of the i-th layer, ω is the half-window width, and ω k is the weight coefficient; since the moving average method is adopted, k = -ω...ω is used to calculate the value of ω k ; finally, substituting into (17) to obtain the finally predicted impedance value.

[0058] For a further technical solution, the specific determination of the correctness of the inversion result in step S4 is as follows: The standard deviation of relative error and the autocorrelation coefficient are dual indicators for evaluating the quality of MCMC inversion; for the calculation of the standard deviation, the relative error needs to be calculated first:

[0059]

[0060] Then, the standard deviation of relative error is calculated through the standard deviation formula;

[0061]

[0062] In the formula, Er represents the relative error, σ represents the standard deviation of relative error, x represents the true impedance value, and x smooth represents the impedance value finally predicted.

[0063] A spatio-temporally correlated adaptive MCMC wave impedance inversion method disclosed by the present invention has the following beneficial effects:

[0064] (1) By introducing a spatio-temporally correlated covariance prior display to describe the spatial correlation of impedance parameters and the vertical continuity of strata, it not only reflects the spatial characteristics of the layered geological structure, significantly improves the physical rationality and numerical robustness of the inversion result, and the prior constraint on the parameter space restricts the fluctuation range of parameters during the inversion process, thereby avoiding the inversion result from falling into local extrema or unstable states; in seismic inversion, impedance parameters usually have strong spatial correlation, and introducing a covariance matrix can effectively describe this correlation, making the inversion process more stable.

[0065] (2) In the construction of the objective function for the fitting degree between the impedance model and the observed data, by adding multiple constraint terms to combine the physical model and statistical prior, it ensures the high accuracy and geological rationality of the inversion result, reduces the non-uniqueness, and at the same time improves the resolution of the interlayer interface, meeting the requirements of fine reservoir characterization under complex formation conditions.

[0066] (3) By using the adaptive rejection jump method to generate candidate points, it significantly shortens the convergence time in the high-dimensional parameter space during inversion, while maintaining the ergodicity of the Markov chain, ensuring the global optimality of the inversion result.

[0067] (4) Check the convergence of the Markov chain through the posterior standard deviation and autocorrelation analysis to ensure the reliability of the inversion results and meet the dual requirements of high-precision seismic exploration and risk control. Description of the Drawings

[0068] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0069] Figure 1 : Flowchart of the inversion method of the present invention;

[0070] Figure 2 : Noiseless post-stack seismic data and initial impedance parameter model input in the embodiment of the present invention;

[0071] Figure 3 : Noisy post-stack seismic data and initial impedance parameter model input in the embodiment of the present invention;

[0072] Figure 4 : Impedance, its stochastic realizations and uncertainty information obtained by the spatio-temporal correlation-based MCMC seismic inversion method in the noiseless case of the embodiment of the present invention;

[0073] Figure 5 : Impedance, its stochastic realizations and uncertainty information obtained by the spatio-temporal correlation-based MCMC seismic inversion method in the case of noise (signal-to-noise ratio of 5) in the embodiment of the present invention. Detailed Embodiments

[0074] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, not all of them. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention. Detailed Embodiment 1

[0076] As Figure 1 shown in the specific flowchart of the inversion algorithm of the present invention:

[0077] 101. Construct the prior: Generate a spatio-temporal correlation covariance matrix through a hidden Markov directed graph, and use a regularized squared exponential Toeplitz covariance matrix so that the covariance matrix can reflect the spatio-temporal correlation of the layered geological structure.

[0078] 102. Objective function design: By combining multiple constraint terms with physical models and statistical priors, including seismic data fitting terms, initial model prior constraint terms, and sparse constraint terms, to ensure the high precision and geological rationality of the inversion results.

[0079] 103. DR: In the first stage, global exploratory samples are generated based on the current covariance and the M-H criterion. When the proposal is rejected, in the second stage, local refined search step sizes are generated by shrinking the covariance, using the delayed rejection strategy.

[0080] 104. AM: Use an adaptive mechanism to dynamically fuse historical sampling information through exponentially weighted moving average, and update the covariance matrix every certain number of iterations, so that the proposal distribution gradually approaches the actual structure of the posterior parameter space.

[0081] 105. Through the repeated iteration of steps (103) and (104), and controlling the maximum number of iterations by the inversion residual, the optimal impedance inversion structure is obtained. Then, an adaptive window smoothing filter is used to constrain the seismic wave impedance inversion results to obtain more accurate inversion results.

[0082] 106. Post-processing of the inversion results: Calculate the posterior standard deviation to quantify the uncertainty of the inversion results and provide a risk reference for exploration decisions; at the same time, check the convergence of the Markov chain through autocorrelation analysis to ensure the reliability of the inversion results.

[0083] This embodiment specifically adopts the following working steps to implement the above technical solution: Based on the spatio-temporal correlated covariance matrix of the hidden Markov directed graph model, that is, using the Toeplitz structure combined with the squared exponential kernel function to construct the prior covariance for subsequent sample generation, constructing an objective function for the fitting degree of the impedance model with constraints and observed data as the optimization object, making point-taking judgments in two stages under the DR strategy, adaptively updating the covariance matrix during the iteration, continuously iterating to stabilize the Markov chain and obtain the predicted impedance value, using an adaptive window smoothing filter (EPS) to constrain the obtained impedance, calculating the posterior standard, and performing autocorrelation analysis to ensure the correctness of the inversion.

[0084] The inversion strategy of this embodiment is based on the Bayesian framework. By introducing the hidden Markov chain model to construct the spatio-temporal correlated covariance as a prior constraint, it shows the spatial correlation of impedance parameters and the vertical continuity of strata, combines the physical model to ensure that the prediction results meet the matching of the observed seismic data, and at the same time is deeply integrated with the delayed rejection adaptive MCMC algorithm, realizing the collaborative innovation of physical model constraints and statistical optimization in seismic wave impedance inversion; meeting the requirements of high-precision seismic inversion. Specific Embodiment Two

[0086] Based on the technology of Specific Embodiment 1, a spatio-temporal related adaptive MCMC wave impedance inversion method is detailed, including the following steps:

[0087] S1: Construct a spatio-temporal related covariance matrix based on the hidden Markov chain model: Taking the vertical continuity of the geological sequence as the core prior, using the Toeplitz structure combined with the squared exponential kernel function, and introducing a regularization term at the same time, so that the covariance matrix can not only reflect the spatio-temporal correlation of the layered geological structure, but also avoid the sampling failure caused by the ill-conditioned matrix, and display the stationary spatial correlation of the impedance parameters; The specific construction of the covariance is as follows:

[0088] S11: Data preparation and processing: Assume that the seismic wavelet is known before inversion. First, construct a true acoustic impedance model containing 350 layers to simulate the layered structure of the underground formation; By expanding the impedance values of 15 coarse layers into 350 finer layers, the thickness differences of different geological layers are reflected; Then, according to the acoustic impedance model, calculate the reflection coefficients at the adjacent formation interfaces;

[0089]

[0090] In the formula, x represents the well logging impedance parameter, and r represents the reflection coefficient; Through the convolution of the seismic wavelet and the reflection coefficient:

[0091] s = w × r#(2)

[0092] In the formula, w represents the Ricker wavelet, × represents the convolution operator, and s represents the synthetic seismic record; It is used to generate the synthetic seismic record; Simulate the actual seismic data and add noise to be closer to the real observation;

[0093] s noise = s + ∈, N~(0, σ 2 I)#(3)

[0094] In the formula, s noise is the seismic record after adding noise, σ 2 is used to control the noise intensity, and ∈ is the noise term;

[0095] S12: Construction of the initial covariance: According to the correlation hypothesis, construct an initial spatio-temporal related covariance matrix through the space of the formation, and define the covariance decay relationship between parameters based on the Gaussian kernel function;

[0096]

[0097] In the formula, d represents the layer spacing, l represents the correlation length, δ 2Denote the variance as, and \(C\) is the one-dimensional covariance column vector constructed at this time. Then, through Toeplitz matrix transformation, the one-dimensional covariance vector is extended to a symmetric positive definite covariance matrix, where its diagonal elements represent the variances of the parameters themselves, and the off-diagonal elements decay exponentially as the layer spacing increases. Finally, a regularization term is added to ensure the numerical stability of the matrix for subsequent Gaussian kernel decomposition;

[0098]

[0099] Convert the vector \(C\) into a Toeplitz matrix to generate a symmetric covariance matrix, and add a regularization term as follows:

[0100] \(C_0 = \text{Toeplitz}(C)+\varepsilon I\) #(6)

[0101] In the formula, \(C_0\) is the initial covariance, \(i\) and \(j\) are the row and column indices of the matrix \(C\), \(\varepsilon\) is the regularization coefficient, and \(I\) is the identity matrix of the same type as \(C_0\). This design enables the proposed distribution to balance local exploration and global correlation, providing an initial perturbation pattern driven by geological significance for MCMC sampling. This design enables the proposed distribution to balance local exploration and global correlation, providing an initial perturbation pattern driven by geological significance for MCMC sampling.

[0102] S2: Construction of the objective function for the fitting degree between the impedance model and the observed data: As the optimization objective of the inversion, this method combines the physical model and statistical prior through multiple constraint terms, including the seismic data fitting term, the initial model prior constraint term, and the sparse constraint term, to ensure the high precision and geological rationality of the inversion results. The objective function is the core in Bayesian inversion, used to quantify the matching degree between the impedance in the model and the observed data, while incorporating prior constraints. First, construct the difference matrix:

[0103]

[0104] After that, construct the objective function for the fitting degree between the impedance model and the observed data, and promote the sparsity of the reflection coefficient through the L1 norm, implicitly assuming geology; and the L2 norm constrains the gradient of the impedance model, making the inversion results smooth, suppressing high-frequency noise, and pulling the inversion results towards the initial guess model to avoid the model deviating from the reasonable range due to data noise;

[0105]

[0106] In the formula, \(\theta\) is the current impedance value, \(\beta\) i＝1 , 2 is the coefficient of the constraint term, and \(ss(\theta)\) is the objective function at the current impedance value.

[0107] S3: Iterative inversion: By repeating the iteration and controlling the maximum number of iterations through the inversion residual, an optimal impedance inversion structure is obtained; the iterative content includes: an adaptive mechanism, specifically: by means of an exponentially weighted strategy moving average dynamics, fusing historical sampling information, updating the covariance matrix every certain number of iterations, so that the proposal distribution gradually approaches the actual geometric form of the posterior parameter space; the sampling process of the adaptive mechanism is specifically to improve the sampling efficiency through a two-stage dynamic adjustment of the rejection jump mechanism DRAM:

[0108] S31: In the first stage, generate global exploratory samples based on the current covariance and the M-H criterion;

[0109] The first proposal:

[0110] θ t+1 = θ t + τ, τ ~ N(0, C t )#(9)

[0111] where θ t is the current model, C t is the current covariance matrix, and θ t+1 is the candidate model generated;

[0112] Calculate the acceptance probability according to the M-H criterion:

[0113]

[0114] Since the proposal distribution is a uniform distribution, it is symmetric, and based on Bayesian theory, the above formula is further transformed into:

[0115]

[0116] In the formula, α1 is the acceptance probability, is the scaling factor;

[0117] S32: When the proposal is rejected; in the second stage, by shrinking the current covariance, generate a local refined search step size, generate a sub-optimal proposal in the local neighborhood, and combine physical constraints to achieve refined search with controllable exploration risk, while maintaining statistical independence from the initial proposal; during this process, the candidate model needs to meet the preset physical constraint conditions (such as impedance value range) to ensure the rationality of the inversion result; this mechanism, through a multi-scale perturbation strategy, not only retains the ability to jump out of local extrema but also enhances the detailed search of the local area, thus significantly improving the sampling efficiency and convergence speed while maintaining the statistical properties of the Markov chain and ergodicity.

[0118] If the above candidate value is not accepted, then enter the second proposal:

[0119] θ t′ +1 = θ t + τ ′ , τ ′ ~ N(0, γ m C t )#(12)

[0120] where θ t is the current model, C t is the current covariance matrix, γ m is the shrinkage factor of the covariance, and θ t ′ +1 is the generated candidate model;

[0121] Similarly to the above, we get:

[0122]

[0123] At this time, we choose whether to accept the second candidate value with probability α2. If accepted, then θ t ′ +1 = θ t , if not accepted, then continue to keep θ t and perform the next iteration.

[0124] Among them, the specific process of adaptively adjusting the covariance is:

[0125] When the number of iterations exceeds the preset threshold, at fixed intervals, impedance samples in the recent window are extracted from the historical sampling results, and its empirical covariance matrix is calculated;

[0126]

[0127] Through the above formula transformation, the covariance formula for the (j + 1)-th time can be obtained:

[0128]

[0129] When updating, the exponential weighted moving average method is adopted to perform weighted fusion on the historical covariance matrix and the current window covariance; among them, the historical weight coefficient controls the retention ratio of old information, and gives higher weight to long-term statistical characteristics; the current window covariance reflects the recent sampling trend and enhances the adaptability to non-stationary posterior distributions; to avoid numerical problems caused by insufficient samples or high-dimensional singularity of the covariance matrix, a small amount of diagonal regularization term is added to it after each update to enforce matrix positive definiteness; at the same time, the matrix ill-condition is detected through eigenvalue decomposition, and the diagonal elements are further adjusted if necessary to ensure the physical validity of the proposal distribution;

[0130] Considering that the covariance may mutate and greatly affect the results, the exponential weighted moving average is introduced:

[0131]

[0132] In the above equations (14), (15), and (16), C j represents the covariance matrix of the j-th update, t0 is the initial threshold, and s d represents the scaling factor, and ε represents the regularization parameter. represents the update weight, and I d is the identity matrix of the same type as C0;

[0133] Through repeated iteration and controlling the maximum number of iterations by the inversion residual, the optimal impedance inversion structure is obtained;

[0134] After that, the seismic wave impedance inversion result is constrained by adaptive window smoothing filtering (Equation 17) to obtain a more accurate inversion result;

[0135] The formula for EPS is as follows:

[0136]

[0137] where x1 is the predicted impedance of the i-th layer, ω is the half-window width, and ω k is the weight coefficient; since the moving average method is adopted, k = -ω...ω is used to calculate the value of ω k ; finally, substituting into (17) to obtain the finally predicted impedance value.

[0138] S4: Post-processing of the inversion result: Calculate the posterior standard deviation to quantify the uncertainty of the inversion result and provide a risk reference for exploration decision-making; at the same time, check the convergence of the Markov chain through autocorrelation analysis to ensure the reliability of the inversion result. Specifically, to determine the correctness of the inversion result:

[0139] The relative error standard deviation and the autocorrelation coefficient are dual indicators for evaluating the quality of MCMC inversion; for the calculation of the standard deviation, the relative error needs to be calculated first:

[0140]

[0141] Then, the relative error standard deviation is calculated through the standard deviation formula;

[0142]

[0143] In the formula, Er represents the relative error, σ represents the relative error standard deviation, x represents the true impedance value, and x smooth represents the impedance value finally predicted.

[0144] In this embodiment (1), by introducing a spatio-temporal related covariance prior display, the spatial correlation of impedance parameters and the vertical continuity of strata are described, which not only reflects the spatial characteristics of layered geological structures, significantly improves the physical rationality and numerical robustness of the inversion results for the prior constraints on the parameter space, restricts the fluctuation range of parameters during the inversion process, and thus avoids the inversion results from falling into local extrema or unstable states; in seismic inversion, impedance parameters usually have strong spatial correlation, and introducing a covariance matrix can effectively describe this correlation, making the inversion process more stable.

[0145] (2) In the construction of the objective function for the fitting degree between the impedance model and the observed data, by adding multiple constraint terms to combine the physical model and statistical prior, the high-precision and geological rationality of the inversion results are ensured, the multi-solution problem is reduced, and at the same time, the resolution of the interlayer interface is improved to meet the requirements of fine reservoir characterization under complex formation conditions.

[0146] (3) By using the method of adaptive rejection jump to generate candidate points, the convergence time of the high-dimensional parameter space is significantly shortened during inversion, while maintaining the ergodicity of the Markov chain to ensure the global optimality of the inversion results.

[0147] (4) By checking the convergence of the Markov chain through posterior standard deviation and autocorrelation analysis, the reliability of the inversion results is ensured, meeting the dual requirements of high-precision seismic exploration and risk control. Specific Embodiment Three

[0149] Based on the above embodiment, considering that the wave impedance is blocky, the adaptive mechanism uses adaptive window smoothing filtering to constrain the seismic wave impedance inversion results to obtain the optimal impedance solution, and finally evaluates the quality of MCMC inversion through the standard deviation of relative error and autocorrelation coefficient.

[0150] Certainly, without departing from the spirit and essence of the present invention, those skilled in the art should be able to make various corresponding changes and deformations according to the present invention, but these corresponding changes and deformations should fall within the protection scope of the appended claims of the present invention.

Claims

1. A spatio-temporal related adaptive MCMC wave impedance inversion method, characterized in that: It includes the following steps: S1: Construct a spatio-temporal covariance matrix: Taking the vertical continuity of the geological sequence as the core prior, using the Toeplitz structure combined with the squared exponential kernel function, and introducing a regularization term at the same time; S2: Construct an objective function for the fitting degree between the impedance model and the observed data: Combining the physical model and the statistical prior through multiple constraint terms; S3: Iterative inversion: Obtain the optimal impedance inversion structure by repeating iterations and controlling the maximum number of iterations through the inversion residual; The iterative content includes: an adaptive mechanism, specifically: through the exponentially weighted strategy moving average dynamics, fusing historical sampling information, updating the covariance matrix every certain number of iterations, so that the proposal distribution gradually approaches the actual geometric shape of the posterior parameter space; S4: Post-processing of the inversion result: Calculate the posterior standard deviation, quantify the uncertainty of the inversion result, and provide a risk reference for exploration decision-making; At the same time, check the convergence of the Markov chain through autocorrelation analysis to ensure the reliability of the inversion result.

2. The adaptive MCMC wave impedance inversion method related to space-time according to claim 1, wherein: The sampling process of the adaptive mechanism in step S3 is specifically to improve the sampling efficiency through the two-stage dynamic adjustment of the rejection jump mechanism DRAM: S31: In the first stage, generate global exploratory samples based on the current covariance and the M-H criterion; S32: When the proposal is rejected, in the second stage, shrink the current covariance, generate a local refined search step size, generate sub-optimal proposals within the local neighborhood, and combine physical constraints to achieve refined search with controllable exploration risks.

3. An adaptive MCMC wave impedance inversion method related to space-time according to claim 1 or 2, characterized in that: The adaptive mechanism in step S3 is specifically: Considering that the wave impedance is blocky, use adaptive window smoothing filtering to constrain the seismic wave impedance inversion result.

4. An adaptive MCMC wave impedance inversion method related to space-time according to claim 3, characterized in that: The construction of the covariance in step S1 is specifically: S11: Data preparation and processing: Construct a real acoustic impedance model containing 350 layers to simulate the layered structure of the underground formation; According to the acoustic impedance model, calculate the reflection coefficient at the adjacent formation interface; In the formula, x represents the logging impedance parameter, and r represents the reflection coefficient; Through the convolution of the seismic wavelet and the reflection coefficient: s = w × r #(2) In the formula, w represents the Ricker wavelet, × represents the convolution operator, and s represents the synthetic seismic record; Used to generate synthetic seismic records, simulate actual seismic data, and add noise to be closer to real observations; s noise = s + ∈, N~(0, σ 2 I)#(3) where s noise is the seismic record after adding noise, σ 2 is used to control the noise intensity, and ∈ is the noise term; S12: Construction of the initial covariance: According to the correlation hypothesis, construct an initial spatio-temporal covariance matrix through the space of the formation, and define the covariance decay relationship between parameters based on the Gaussian kernel function; where d represents the layer spacing, l represents the correlation length, δ 2 represents the variance, and C is the one-dimensional covariance column vector constructed at this time; Then, through the Toeplitz matrix transformation, expand the one-dimensional covariance vector into a symmetric positive definite covariance matrix, so that its diagonal elements represent the variance of the parameters themselves, and the non-diagonal elements decay exponentially with the increase of the layer spacing; Finally, add a regularization term to ensure the numerical stability of the matrix for subsequent Gaussian kernel decomposition; Convert the vector C into a Toeplitz matrix, generate a symmetric covariance matrix, and add a regularization term as follows: C0 = Toeplitz(C) + εI #(6) In the formula, C0 is the initial covariance, i and j are the row and column indices of the matrix C, ε is the regularization coefficient, and I is the identity matrix of the same type as C0.

5. The spatio-temporal correlation-based adaptive MCMC wave impedance inversion method according to claim 4, characterized in that: The construction of the objective function in step S2 is specifically as follows: The objective function is the core in Bayesian inversion, which is used to quantify the matching degree between the impedance in the model and the observed data, and at the same time incorporate prior constraints. First, construct a difference matrix: After that, construct the objective function for the fitting degree between the impedance model and the observed data, and promote the sparsity of the reflection coefficient through the L1 norm, implicitly assuming geology; and use the L2 norm to constrain the gradient of the impedance model, making the inversion result smooth, suppressing high-frequency noise, and pulling the inversion result towards the initial guess model to avoid the model deviating from the reasonable range due to data noise; where θ is the current impedance value, β i＝1 , 2 is the coefficient of the constraint term, and ss(θ) is the objective function at the current impedance value.

6. The spatio-temporal correlation-based adaptive MCMC wave impedance inversion method according to claim 5, wherein: The specific content of S3 is as follows: S31: The first proposal: θ t+1 = θ t + τ, τ ~ N(0, C t )#(9) where θ t is the current model, C t is the current covariance matrix, and θ t+1 generates a candidate model; Calculate the acceptance probability according to the M-H criterion: Since the proposal distribution is a uniform distribution, it is symmetric. At the same time, based on Bayesian theory, the above formula is further transformed into: where α1 is the acceptance probability, is the scaling factor; S32: If the above candidate value is not accepted, enter the second proposal: θ′ t+1 = θ t + τ′, τ′ ~ N(0, γ m C t )#(12) where θ t is the current model, C t is the current covariance matrix, γ m is the reduction factor of the covariance, θ′ t+1 is the generated candidate model; Similarly to the above, we get: At this time, the second candidate value is selected with probability α2. If it is accepted, then θ′ t+1 = θ t , and if it is not accepted, θ is continued to be maintained t and the next iteration is carried out.

7. An adaptive MCMC wave impedance inversion method related to space-time according to claim 6, characterized in that: The specific content of the adaptive adjustment of the covariance process in step S3 is as follows: When the number of iterations exceeds the preset threshold, at fixed intervals, extract the impedance samples in the recent window from the historical sampling results and calculate their empirical covariance matrix: Through the above transformation, the covariance formula for the (j + 1)-th time can be obtained: When updating, use the exponential weighted moving average method to fuse the historical covariance matrix and the current window covariance; among them, the historical weight coefficient controls the retention ratio of old information and assigns a higher weight to the long-term statistical characteristics; the current window covariance reflects the recent sampling trend and enhances the adaptability to the non-stationary posterior distribution; to avoid numerical problems caused by insufficient samples or high-dimensional singularity of the covariance matrix, add a small diagonal regularization term to it after each update to force the matrix to be positive definite; at the same time, detect the ill-conditioning of the matrix through eigenvalue decomposition and further adjust the diagonal elements if necessary to ensure the physical validity of the proposal distribution; Considering that the covariance may mutate and greatly affect the result, introduce the exponential weighted moving average: In the above equations (14), (15), and (16), C j represents the covariance matrix updated at the j-th time, t0 is the initial threshold, and s d represents the scaling factor, ε represents the regularization parameter, represents the updated weight, and I d is the identity matrix of the same type as C0; Obtain the optimal impedance inversion structure by repeating the iteration and controlling the maximum number of iterations through the inversion residual; After that, use adaptive window smoothing filtering (Equation 17) to constrain the seismic wave impedance inversion result to obtain a more accurate inversion result; The formula for EPS is as follows: where x1 is the predicted impedance of the i-th layer, ω is the half-window width, and ω k is the weight coefficient; since the moving average method is adopted, thus is used to calculate the value of ω k ; finally, substituting into (17) to obtain the finally predicted impedance value.

8. A spatio-temporal related adaptive MCMC wave impedance inversion method according to claim 7, characterized in that: The specific content of step S4 to determine the correctness of the inversion result is as follows: The standard deviation of the relative error and the autocorrelation coefficient are dual indicators for evaluating the quality of MCMC inversion; for the calculation of the standard deviation, first calculate the relative error: Then calculate the standard deviation of the relative error through the standard deviation formula; where Er represents the relative error, σ represents the standard deviation of the relative error, x represents the true impedance value, and x smooth represents the impedance value finally predicted.

Images