A prestack seismic stochastic inversion method based on frequency-domain co-simulation

Abstract

The present invention belongs to the technical field of geophysical oil and gas exploration, and particularly relates to a prestack seismic stochastic inversion method based on frequency-domain co-simulation. First, binary joint probability distributions of the P-wave to S-wave velocity ratio with the P-wave modulus and the P-wave to S-wave velocity ratio with density are established using well logging data. An initial P-wave to S-wave velocity ratio model is established using FFT-MA simulation, and initial P-wave modulus and density models are established using FFT-MA co-simulation. Then, an objective function formula is established according to Bayes' theorem. Next, an iterative update algorithm is used to generate a new P-wave to S-wave velocity ratio model, which is used as soft data. An iterative update algorithm is used to generate new P-wave modulus and density models, and the Metropolis criterion is used to determine whether to retain the new models. This step is repeated until the maximum number of iterations is reached or the objective function value is less than the threshold value. The present invention combines the fast Fourier moving average co-simulation method with the iterative update method, improving the accuracy of the inversion parameter simulation and the convergence and accuracy in the inversion iteration process.

Description

Technical Field

[0001] The present invention belongs to the technical field of geophysical oil and gas exploration, and particularly relates to a prestack seismic stochastic inversion method based on frequency-domain co-simulation. Background Art

[0002] Seismic inversion can use seismic data to predict relevant information such as reservoir elastic parameters, physical properties, lithology, porosity, water saturation, etc., laying a foundation for oil and gas resource evaluation. The inversion problem is usually a non-linear problem. Due to the limitations of the seismic inversion method itself, there are large errors in solving the inversion problem using conventional linear inversion methods in many cases, and the resolution of the inversion results is difficult to meet the requirements of thin reservoir prediction. Compared with the linear inversion method, stochastic inversion does not require obtaining an analytical solution to the inversion problem, and can effectively utilize the high-frequency information in well logging data to improve the vertical resolution of the inversion results, providing strong technical support for identifying thin reservoirs.

[0003] Stochastic inversion methods usually regard the inversion problem as an iterative optimization problem. They use stochastic simulation to extract samples from the prior distribution and adopt optimization algorithms (such as Metropolis-Hastings, simulated annealing) to optimize the samples, and finally obtain inversion results that converge to the posterior probability distribution. Stochastic simulation is the core of stochastic inversion, which restricts the accuracy and efficiency of inversion. According to the simulation method, it is mainly divided into two categories. One is the sequential simulation method, mainly including sequential Gaussian simulation, direct sequential simulation, etc. The sequential simulation method, based on well logging data, analyzes the spatial distribution trend of parameters and uses the weighted average method to calculate the simulation mean and variance at each grid point, thereby establishing a high-resolution reservoir parameter model. However, this method has a large amount of calculation and it is difficult to establish a reliable reservoir parameter model in the absence of well data. Another stochastic simulation method is the frequency-domain simulation method, such as fast Fourier transform moving average (FFT-MA) simulation, stochastic medium modeling, etc. This method can establish a reservoir parameter model reflecting the seismic spatial structure characteristics by estimating the covariance and autocorrelation matrix from seismic data. Compared with the traditional sequential Gaussian simulation, this method has higher computational efficiency. However, it is difficult for the frequency-domain simulation method to effectively realize co-simulation between multiple parameters, which restricts its application in prestack AVA inversion.

[0004] Optimization algorithms can optimize the samples generated by stochastic simulation and make them gradually converge to the posterior probability distribution through continuous iteration. Markov chain Monte Carlo (MCMC) algorithms and their improved algorithms, genetic algorithms, simulated annealing and other algorithms are commonly used optimization methods in stochastic inversion. However, these algorithms are all global optimization algorithms with low computational efficiency. Summary of the Invention

[0005] In view of the various deficiencies of the prior art, the inventor has researched and designed a prestack seismic stochastic inversion method based on frequency-domain co-simulation in long-term practice, which combines the FFT-MA co-simulation and the iterative update method to improve the accuracy and computational efficiency of prestack stochastic inversion.

[0006] The stochastic inversion method of the present invention includes the following steps:

[0007] Step 1: Using well logging data, establish the binary joint probability distributions of the P-wave to S-wave velocity ratio with the P-wave modulus and the P-wave to S-wave velocity ratio with density, respectively.

[0008] Step 2: Use FFT-MA simulation to establish an initial P-wave to S-wave velocity ratio model.

[0009] Step 3: Based on the P-wave to S-wave velocity ratio model, combine the joint probability distributions of the P-wave to S-wave velocity ratio with the P-wave modulus and the P-wave to S-wave velocity ratio with density, and calculate the value ranges of the P-wave modulus and density.

[0010] Step 4: Using the value ranges of the P-wave modulus and density as constraints, use FFT-MA co-simulation to establish an initial P-wave modulus and density model.

[0011] Step 5: According to Bayes' theorem, establish the objective function formula and calculate the initial objective function value.

[0012] Step 6: Based on the existing P-wave to S-wave velocity ratio model, use the iterative update algorithm to generate a new P-wave to S-wave velocity ratio model.

[0013] Step 7: Based on the new P-wave to S-wave velocity ratio model, calculate the value ranges of the P-wave modulus and density at each grid point.

[0014] Step 8: Based on the existing P-wave modulus and density model, using the new P-wave to S-wave velocity ratio model as soft data, use the iterative update algorithm to generate a new P-wave modulus and density model.

[0015] Step 9: Calculate the objective function value corresponding to the new model, compare it with the objective function value corresponding to the original model, and use the Metropolis criterion to determine whether to retain the new model.

[0016] Step 10: Repeat Steps 6 to 9 until the maximum number of iterations is reached or the objective function value is less than the threshold value.

[0017] Further, in Step 1, the binary joint probability distributions of the P-wave to S-wave velocity ratio with the P-wave modulus and the P-wave to S-wave velocity ratio with density are established respectively, and the specific method is as follows:

[0018] According to the well logging P-wave modulus, P-wave to S-wave velocity ratio, and density curves, establish the P-wave to S-wave velocity ratio and the P-wave modulus The joint probability distribution of two variables , the ratio of P-wave velocity to S-wave velocity and density The joint probability distribution of two variables .

[0019] Furthermore, in step two, an initial ratio model of P-wave velocity to S-wave velocity is established by using FFT-MA simulation. The specific method is as follows:

[0020] Based on the FFT-MA simulation method, the simulation model can be expressed in the following form:

[0021] Formula 1

[0022] Formula 2

[0023] In the formula, is the mean of the simulation model, is the standard deviation of the logging data, represents the conjugate root of the covariance function, represents Gaussian white noise with a mean of 0 and a variance of 1.

[0024] Furthermore, in step three, the value ranges of the P-wave modulus and density are calculated. The specific method is as follows:

[0025] Taking the ratio of P-wave velocity to S-wave velocity as the secondary variable and the P-wave modulus and density as the primary variables. Based on the simulation results of the ratio of P-wave velocity to S-wave velocity , an interval is defined. Therefore, the conditional probability distribution of the primary variable can be expressed as:

[0026] Formula 3

[0027] where represents the length of the given secondary variable interval. The Gaussian distribution fitting is performed on the conditional probability distribution , that is

[0028] Formula 4

[0029] In the formula, is the mean of the Gaussian distribution, is the variance. Therefore, when the simulation result of the secondary variable is , the value range of the primary variable is .

[0030] Further, in Step 4, an initial P-wave modulus and density model is established using FFT-MA co-simulation. Taking the P-to-S wave velocity ratio as the secondary variable, and the P-wave modulus and density as the primary variables respectively, the following calculations are carried out separately:

[0031] Sub-step 4.1: First, simulate the secondary variable through FFT-MA simulation to obtain the simulation result of the secondary variable , and then establish the following co-simulation formula to obtain the simulation result of the primary variable , realizing the co-simulation of multiple inversion parameters:

[0032] Equation 5

[0033] Equation 6

[0034] In the formula, is the mean of the simulation result of the primary variable, and are the variances of the primary variable and the secondary variable respectively, is the correlation coefficient between the primary variable and the secondary variable;

[0035] Sub-step 4.2: According to the value range of the primary variable being , at this time the final co-simulation result can be written as:

[0036] Equation 7

[0037] Further, in Step 5, calculate the initial objective function value, including the following sub-steps:

[0038] Sub-step 5.1: Based on the Zeoppritz approximation and the seismic convolution model, establish the forward relationship between the inversion parameters and the seismic records. Based on the Aki-Richards approximation, derive the seismic reflection coefficient equation based on the P-wave modulus, P-to-S wave velocity ratio, and density:

[0039] Equation 8

[0040] In the formula, is the incident angle; , , respectively represent the average values of the P-wave modulus, P-to-S wave velocity ratio, and density underground; , , respectively represent the changes in the P-wave modulus, P-to-S wave velocity ratio, and density on both sides of the formation interface;

[0041] Seismic reflection coefficient and seismic data The mathematical relationship can be characterized by using the convolution formula:

[0042] Formula 9

[0043] Wherein, represents the wavelet matrix;

[0044] By combining Formula 8 and Formula 9, the forward seismic relationship between elastic parameters such as P-wave modulus, P-S wave velocity ratio, and density and seismic records can be established;

[0045] Sub-step 5.2: According to Bayes' theorem, the posterior probability distribution can be written as the product of the prior distribution and the likelihood function , that is:

[0046] Formula 10

[0047] Assume that both the prior distribution and the likelihood function satisfy the Gaussian distribution, that is:

[0048] Formula 11

[0049] Formula 12

[0050] Therefore, according to Formulas 10, 11, and 12, the objective function can be written as:

[0051] Formula 13

[0052] Wherein, represents the forward equation, that is, Formulas 8 and 9; represents the covariance matrix of seismic noise, represents the covariance matrix of inversion parameters; represents the prior smoothing constraint model, which is usually obtained by interpolating and extrapolating well logging data; represents the number of sampling points.

[0053] Furthermore, in Step 6, an iterative update algorithm is used to generate a new P-S wave velocity ratio model. The specific method is as follows:

[0054] According to Formulas 1 and 7, the new model can be expressed as:

[0055] Formula 14

[0056] Wherein, is the weight coefficient; is the current model. For the secondary variable P-S wave velocity ratio, For the main variables of longitudinal wave modulus and density, 。

[0057] Furthermore, in step nine, the Metropolis criterion is used to determine whether to retain the new model, including the following sub-steps:

[0058] Sub-step 9.1, combine the inversion likelihood function and the prior constraint model , and establish the following objective function:

[0059] Equation 15

[0060] In the formula, represents the forward equation, that is, Equation 8 and Equation 9; represents the covariance matrix of seismic noise, represents the covariance matrix of inversion parameters;

[0061] Sub-step 9.2, in order to avoid the inversion objective value falling into local extrema, the Metropolis criterion is used to optimize the generated new model, and the following acceptance probability expression is established:

[0062] Equation 16

[0063] In the formula, represents the inversion objective function, is an adjustable parameter. Generate a random number between 0 and 1. When the acceptance probability is greater than the random number, update the current optimal solution. When the acceptance probability is less than the random number, retain the original optimal solution.

[0064] The beneficial effects of the present invention are:

[0065] By combining the fast Fourier sliding average covariance simulation method with the iterative update method, compared with the traditional fast Fourier sliding average simulation, the fast Fourier sliding average covariance simulation improves the accuracy of inversion parameter simulation by fully utilizing the correlation between inversion parameters; while the iterative update method can improve the stability in the inversion iteration process, effectively improving the convergence and accuracy of inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 is the two-dimensional P-S wave velocity ratio, longitudinal wave modulus and density model of the present invention.

[0067] Figure 2 is the partial stacked seismic data corresponding to the two-dimensional model of the present invention.

[0068] Figure 3 is the two-dimensional joint probability distribution between the P-S wave velocity ratio and the longitudinal wave modulus and the conditional probability distribution of the longitudinal wave modulus of the present invention.

[0069] Figure 4 It is the FFT-MA simulation result of the P-wave to S-wave velocity ratio of the present invention.

[0070] Figure 5 It is the comparison between the FFT-MA simulation result and the FFT-MA co-simulation result of the P-wave modulus of the present invention.

[0071] Figure 6 It is the comparison between the relative error of the FFT-MA simulation and the relative error of the FFT-MA co-simulation of the P-wave modulus of the present invention.

[0072] Figure 7 It is the inversion result of the two-dimensional model of the present invention based on the iterative update algorithm.

[0073] Figure 8 It is the inversion result of the two-dimensional model of the present invention based on the traditional stochastic inversion strategy.

[0074] Figure 9 It is the single-trace comparison of the inversion results of the present invention based on the iterative update algorithm and the traditional stochastic inversion strategy.

[0075] Figure 10 It is the variation of the objective function of different inversion methods of the present invention with the number of iterations.

[0076] Figure 11 It is the partially stacked seismic data of the actual near, middle, and far angles of the present invention.

[0077] Figure 12 It is the simulation results of the P-wave to S-wave velocity ratio, P-wave modulus, and density of the actual two-dimensional survey line of the present invention.

[0078] Figure 13 It is the inversion results of the P-wave to S-wave velocity ratio, P-wave modulus, and density of the actual two-dimensional survey line of the present invention.

[0079] Figure 14 It is the comparison between the inversion result of the well-side trace at the position of Well B of the present invention and the actual logging curve. Detailed implementation manners

[0080] The present invention will be further described in detail below in conjunction with the accompanying drawings and implementation examples. By describing these implementation examples in sufficient detail, those skilled in the art can understand and practice the present invention. Without departing from the gist and scope of the present invention, logical, implementation, and other changes can be made to the implementation. Therefore, the following detailed description should not be construed in a limiting sense, and the scope of the present invention is only defined by the claims.

[0081] Aiming at the deviation problem existing in the measurement of incident short-wave radiation on the current buoy, the present invention proposes a pre-stack seismic stochastic inversion method based on frequency-domain co-simulation, which specifically includes the following steps:

[0082] Step 1: Using well logging data, establish the bivariate joint probability distributions of the P-wave to S-wave velocity ratio versus the P-wave modulus and the P-wave to S-wave velocity ratio versus density, respectively.

[0083] Based on the well logging P-wave modulus, P-wave to S-wave velocity ratio, and density curves, establish the bivariate joint probability distribution of the P-wave to S-wave velocity ratio versus the P-wave modulus and the bivariate joint probability distribution of the P-wave to S-wave velocity ratio versus density .

[0084] Step 2: Use FFT-MA simulation to establish an initial P-wave to S-wave velocity ratio model.

[0085] Stochastic inversion is an iterative optimization problem, and stochastic simulation is the core of stochastic inversion, which can establish a series of model samples. Previous studies have found that compared with traditional geostatistical simulation methods (such as sequential Gaussian simulation), the FFT-MA simulation method has higher computational efficiency.

[0086] Based on the FFT-MA simulation method, the simulation model can be characterized in the following form:

[0087] Equation 1

[0088] Equation 2

[0089] In the formula, is the mean of the simulation model, is the standard deviation of the well logging data, represents the conjugate root of the covariance function, represents Gaussian white noise with a mean of 0 and a variance of 1.

[0090] Step 3: Based on the P-wave to S-wave velocity ratio model, combined with the joint probability distributions of the P-wave to S-wave velocity ratio versus the P-wave modulus and the P-wave to S-wave velocity ratio versus density, calculate the value ranges of the P-wave modulus and density.

[0091] Take the P-wave to S-wave velocity ratio as the secondary variable and the P-wave modulus and density as the primary variables. Based on the simulation results of the P-wave to S-wave velocity ratio , define an interval , so the conditional probability distribution of the primary variable can be characterized as:

[0092] Equation 3

[0093] where represents the length of a given secondary variable interval. For the conditional probability distribution perform Gaussian distribution fitting, that is

[0094] Formula 4

[0095] wherein, is the mean of the Gaussian distribution, is the variance. Therefore, when the simulation result of the secondary variable is , the value range of the primary variable is .

[0096] Step 4: Using the value ranges of the P-wave modulus and density as constraints, establish an initial P-wave modulus and density model by FFT-MA co-simulation.

[0097] Pre-stack AVA inversion usually requires simultaneous inversion of three elastic parameters (such as P-wave modulus, P-to-S wave velocity ratio, and density), and there is a certain correlation between the inversion parameters. However, the FFT-MA simulation does not consider the correlation between the inversion parameters, resulting in an increase in simulation uncertainty and a corresponding impact on the simulation accuracy. Therefore, we propose an FFT-MA co-simulation method to achieve the co-simulation of multiple inversion parameters.

[0098] Taking the co-simulation of the P-to-S wave velocity ratio and the P-wave modulus as an example, the P-to-S wave velocity ratio is used as the secondary variable, and the P-wave modulus is used as the primary variable. First, simulate the secondary variable through FFT-MA simulation to obtain the simulation result of the secondary variable ; then, the following co-simulation formula can be established to obtain the simulation result of the primary variable (P-wave modulus or density) :

[0099] Formula 5

[0100] Formula 6

[0101] wherein, is the mean of the simulation result of the primary variable, and are the variances of the primary variable and the secondary variable respectively, is the correlation coefficient between the primary variable and the secondary variable.

[0102] However, the simulation result obtained from Formula 5 is only applicable to the case where the secondary variable and the primary variable satisfy linear correlation. According to the joint probability distribution, when the simulation result of the secondary variable is , the value range of the primary variable is . At this time, the final co-simulation result can be written as:

[0103] Formula 7

[0104] Step 5: According to Bayes' theorem, establish the objective function formula and calculate the initial objective function value.

[0105] The forward seismic relationship is the basis for constructing the objective function, mainly including the reflection coefficient equation and the seismic convolution model, which describe the mathematical relationship between reservoir elastic parameters and seismic data. Based on the Zeoppritz approximation and the seismic convolution model, establish the forward relationship between the inversion parameters and the seismic records. Based on the Aki-Richards approximation, derive the seismic reflection coefficient equation based on the P-wave modulus, P-S wave velocity ratio, and density:

[0106] Formula 8

[0107] In the formula, is the incident angle; , , respectively represent the average values of the P-wave modulus, P-S wave velocity ratio, and density underground; , , respectively represent the changes in the P-wave modulus, P-S wave velocity ratio, and density on both sides of the formation interface;

[0108] The seismic reflection coefficient and the seismic data The mathematical relationship can be characterized by the convolution formula:

[0109] Formula 9

[0110] In the formula, represents the wavelet matrix;

[0111] Combining Formula 8 and Formula 9, the forward seismic relationship between elastic parameters such as P-wave modulus, P-S wave velocity ratio, and density and seismic records can be established;

[0112] According to Bayes' theorem, the posterior probability distribution can be written as the product of the prior distribution and the likelihood function , that is:

[0113] Formula 10

[0114] Assume that both the prior distribution and the likelihood function satisfy the Gaussian distribution, that is:

[0115] Formula 11

[0116] Formula 12

[0117] Therefore, according to Formulas 10, 11, and 12, the objective function can be written as:

[0118] Formula 13

[0119] In the formula, represents the forward equation, that is, Formulas 8 and 9; represents the covariance matrix of seismic noise, represents the covariance matrix of inversion parameters; represents the prior smoothing constraint model, which is usually obtained by interpolating and extrapolating well logging data; represents the number of sampling points.

[0120] Step Six, based on the existing P-wave to S-wave velocity ratio model, use the iterative update algorithm to generate a new P-wave to S-wave velocity ratio model.

[0121] Traditional stochastic inversion usually uses stochastic simulation to generate new samples. However, stochastic simulation has strong randomness, which restricts the computational efficiency and accuracy of inversion. To improve the convergence and stability of stochastic inversion, the present invention proposes an iterative update method.

[0122] According to Formulas 1 and 7, the new model can be expressed as:

[0123] Formula 14

[0124] In the formula, is the weight coefficient; is the current model. For the secondary variable of P-wave to S-wave velocity ratio, . For the primary variables of P-wave modulus and density, .

[0125] Step Seven, based on the new P-wave to S-wave velocity ratio model, calculate the value range of P-wave modulus and density at each grid point.

[0126] Step Eight, based on the existing P-wave modulus and density model, use the new P-wave to S-wave velocity ratio model as soft data and use the iterative update algorithm to generate a new P-wave modulus and density model. Specifically, according to Formula 14 in Step Six, the primary variable P-wave modulus and density model can be obtained, that is .

[0127] Step Nine, calculate the objective function value corresponding to the new model and compare it with the objective function value corresponding to the original model, and use the Metropolis criterion to determine whether to retain the new model. Specifically, it includes the following sub-steps:

[0128] Sub-step 9.1: Combine the inversion likelihood function with the prior constraint model , and establish the following objective function:

[0129] Equation 15

[0130] In the formula, represents the forward equation, namely Equation 8 and Equation 9; represents the covariance matrix of seismic noise, represents the covariance matrix of inversion parameters.

[0131] Sub-step 9.2: To avoid the inversion objective value falling into a local extreme, use the Metropolis criterion to optimize the generated new model and establish the following acceptance probability expression:

[0132] Equation 16

[0133] In the formula, represents the inversion objective function, is an adjustable parameter. Generate a random number between 0 and 1. When the acceptance probability is greater than the random number, update the current optimal solution. When the acceptance probability is less than the random number, retain the original optimal solution.

[0134] Step Ten: Repeat Step Six to Step Nine until the maximum number of iterations is reached or the objective function value is less than the threshold value.

[0135] To test the feasibility of the method, the inventor established a two-dimensional P-wave modulus, P-to-S wave velocity ratio, and density model for inversion testing. The CDP range is 1 - 100. Based on the seismic forward relationship, the central angles of some stacked data are set to 7°, 15°, and 24°. Use a 25Hz Ricker wavelet to establish the angular partial stacked data of the two-dimensional model and use it as the observed seismic data. Figure 1 Respectively show the two-dimensional models of P-wave modulus, P-to-S wave velocity ratio, and density. Figure 2 Shows the partially stacked seismic data.

[0136] First, based on the two-dimensional model, establish the binary joint distribution between inversion parameters. Given the simulation result of the P-to-S wave velocity ratio, that is, calculate the conditional probability distribution of the P-wave modulus. This process is as Figure 3 shown. It can be found that the conditional value range of the P-wave modulus at this time is much smaller than the global value range of the P-wave modulus.

[0137] Based on the partially stacked seismic data, estimate the variogram of the model. Use FFT-MA simulation to obtain the simulation result of the P-to-S wave velocity ratio, as Figure 4 shown. On this basis, use FFT-MA co-simulation to obtain the simulation results of the P-wave modulus and density.Figure 5 respectively show the simulated results of the P-wave to S-wave velocity ratio obtained by FFT-MA simulation, the simulated results of the P-wave modulus obtained by FFT-MA simulation, and the co-simulated results of the P-wave modulus obtained by FFT-MA co-simulation. Figure 6 respectively show the relative errors of the simulated results of the P-wave modulus obtained by FFT-MA simulation and FFT-MA co-simulation. According to Figure 5 and Figure 6 it can be found that the relative error of the FFT-MA co-simulated results is lower than that of the FFT-MA simulated results. Therefore, it can be concluded that when simultaneously simulating multiple parameters, the model established by FFT-MA co-simulation is more accurate than the traditional FFT-MA simulation.

[0138] Based on the FFT-MA co-simulation, the iterative update method and the Metropolis criterion are combined to optimize the simulated samples until the objective function reaches the threshold. Figure 7 is the inversion result based on the iterative update method. The traditional random method directly uses stochastic simulation to obtain new model samples. Figure 8 shows the inversion result based on the traditional stochastic inversion method. Figure 9 is the comparison between the inversion results of the 50th trace of the traditional stochastic method and the iterative update method and the true model. From Figure 7 , Figure 8 and Figure 9 it can be seen that the inversion result obtained by the iterative update method is more accurate than that of the traditional stochastic method.

[0139] Figure 10 is the change of the inversion objective function value with the number of iterations. Compared with the traditional stochastic method (4000 iterations), the iterative update method has a faster convergence speed (3000 iterations). The inversion method based on the iterative update method is superior to the traditional stochastic method in terms of simulation effect and computational efficiency. The two-dimensional model test verifies the reliability and stability of this method.

[0140] In order to further verify the feasibility of the method, this method is applied to the estimation of elastic parameters of an actual tight sandstone reservoir. Through reservoir petrophysical analysis, it is found that the reservoir shows low VP / VS values and low P-wave modulus. Therefore, the reservoir and non-reservoir can be distinguished by inverting the P-wave modulus and VP / VS values. In order to improve the computational efficiency of the inversion, the angle summation is performed on the diagonal gathers to obtain partially angle-stacked seismic data, as Figure 11 shown.

[0141] The method proposed in the present invention is used to carry out the simulation and inversion of reservoir parameters. Figure 12 are respectively the simulated results of the P-wave modulus, the P-wave to S-wave velocity ratio, and density. Figure 13The inversion results of the longitudinal wave modulus, VP / VS, and density are shown. The black curves in the profile are the corresponding well curves. It can be found that the inversion results are consistent with the well data and have higher vertical resolution. In addition, this method can accurately identify the thin reservoir between 1.62 and 1.65 s at Well A. Figure 14 This is a comparison between the inversion result of the well-side trace of this method and the traditional deterministic inversion result. The blue line in the figure is the logging data and the observed seismic data, the red dashed line is the inversion result of this method, and the green dashed line is the conventional deterministic inversion result. It can be seen that compared with the traditional deterministic inversion method, the inversion result of the optimized stochastic inversion method is usually consistent with the well data and has higher resolution.

[0142] The present invention has been described in detail above. The above description is only a preferred embodiment of the present invention, and it cannot limit the scope of the implementation of the present invention. That is, all equivalent changes and modifications made according to the scope of this application should still fall within the scope covered by the present invention.

Claims

1. A prestack seismic stochastic inversion method based on frequency domain co-simulation, characterized in that: The following steps are involved: Step 1: Using logging data, establish the binary joint probability distribution of the ratio of longitudinal and transverse wave velocity to longitudinal wave modulus, and the ratio of longitudinal and transverse wave velocity to density respectively; Step 2: Use FFT-MA simulation to establish the initial longitudinal and transverse wave velocity ratio model; Step 3, based on the P-wave velocity ratio model, the value range of the P-wave modulus and density is calculated by combining the joint probability distribution of the P-wave velocity ratio and the P-wave modulus, and the P-wave velocity ratio and the density; Step 4: Using the range of longitudinal wave modulus and density as constraints, the initial longitudinal wave modulus and density model is established using FFT-MA co-simulation; Step 5: According to Bayes' theorem, establish the objective function formula and calculate the initial objective function value; Step 6, based on the existing longitudinal and transverse wave velocity ratio model, using an iterative update algorithm to generate a new longitudinal and transverse wave velocity ratio model; Step 7: Based on the new longitudinal and transverse wave velocity ratio model, the value range of the longitudinal wave modulus and density at each grid point is calculated; Step 8: Based on the existing longitudinal wave modulus and density model, the new longitudinal and transverse wave velocity ratio model is used as soft data, and a new longitudinal wave modulus and density model is generated by using an iterative update algorithm; Step nine, calculate the objective function value corresponding to the new model, compare it with the objective function value corresponding to the original model, and use the Metropolis criterion to determine whether to retain the new model; Step 10: Repeat steps 6 to 9 until the maximum number of iterations is reached or the objective function value is less than the threshold value.

2. The method according to claim 1, characterized in that In step 1, the binary joint probability distribution of the longitudinal and transverse wave velocity ratio and the longitudinal wave modulus, and the longitudinal and transverse wave velocity ratio and density are established respectively. The specific method is: According to the logging P-wave modulus, P-wave velocity ratio and density curve, the binary joint probability distribution F(τ,M) of P-wave velocity ratio τ and P-wave modulus M and the binary joint probability distribution F(τ,ρ) of P-wave velocity ratio τ and density ρ are established respectively.

3. The method according to claim 2, characterized in that In step 2, the initial longitudinal and transverse wave velocity ratio model is established by using FFT-MA simulation. The specific method is as follows: Based on the FFT-MA simulation method, the simulation model m p It can be represented as follows: m p =m0+δ m ·f Formula 1 f=g*z Formula 2 Where m0 is the mean of the simulation model; δ m is the standard deviation of the logging data; g represents the conjugate root of the covariance function; z represents Gaussian white noise, whose mean is 0 and variance is 1.

4. The method according to claim 3, characterized in that In step 3, the range of longitudinal wave modulus and density is calculated. The specific method is: The longitudinal and transverse wave velocity ratio is taken as the secondary variable, and the longitudinal wave modulus and density are taken as the main variables; in the simulation result of the longitudinal and transverse wave velocity ratio m sec Based on this, we define an interval [m sec -0.5ψ,m sec +0.5ψ], so the conditional probability distribution of the main variable p(m pri |m sec ) can be characterized as: p(m pri |m sec )≈p(m pri |m sec -0.5ψ≤m sec ≤m sec +0.5ψ) Formula 3 Where ψ represents the length of the given secondary variable interval; for the conditional probability distribution p(m pri |m sec ) to fit the Gaussian distribution, that is p(m pri |m sec )≈N(m g ,ε g ) Formula 4 In the formula, m g is the mean of the Gaussian distribution, ε g is the variance; therefore, when the secondary variable simulation result is m sec When the value range of the main variable is [m g -ε g ,m g +ε g ].

5. The method according to claim 4, characterized in that In step 4, the initial P-wave modulus and density model is established by using FFT-MA co-simulation, taking the P-wave velocity ratio as the secondary variable and the P-wave modulus and density as the primary variables, respectively, and performing the following calculations: Sub-step 4.1, first simulate the secondary variable through FFT-MA simulation to obtain the simulation result m of the secondary variable sec , and then establish the following co-simulation formula to obtain the simulation result m of the main variable pri , to achieve collaborative simulation of multiple inversion parameters: In the formula, is the mean of the simulation results of the main variable, C pri With C sec are the variances of the main variable and the secondary variable, respectively, and λ is the correlation coefficient between the main variable and the secondary variable; Sub-step 4.2, according to the value range of the main variable [m g -ε g ,m g +ε g ], at this time the final co-simulation result m cosim It can be written as:

6. The method according to claim 5, characterized in that In step 5, the initial objective function value is calculated, including the following sub-steps: Sub-step 5.1: Based on the Zeoppritz approximation and the seismic convolution model, establish the forward relationship between the inversion parameters and the seismic records; based on the Aki-Richards approximation, derive the seismic reflection coefficient equation based on the P-wave modulus, P-wave velocity ratio and density: Where θ is the incident angle; Represent the average values ​​of the P-wave modulus, P-wave velocity ratio and density of the underground respectively; ΔM, Δμ, Δρ represent the changes of the P-wave modulus, P-wave velocity ratio and density on both sides of the stratum interface respectively; The mathematical relationship between the seismic reflection coefficient R(t,θ) and the seismic data d(t,θ) can be characterized by the convolution formula: Where W represents the wavelet matrix; Combining Formula 8 with Formula 9, the seismic forward modeling relationship between the P-wave modulus, P-wave velocity ratio, and density elastic parameters and seismic records can be established; Sub-step 5.2: According to Bayes’ theorem, the posterior probability distribution p(m|d) can be written as the product of the prior distribution p(m) and the likelihood function p(m|d), that is: p(m|d)∝p(m)·p(m|d) Formula 10 Assume that both the prior distribution and the likelihood function satisfy the Gaussian distribution, that is: Therefore, according to formulas 10, 11 and 12, the objective function can be written as: J(m) = [d - A(m)] T (C e ) -1 [d - A(m)] + (m - m′) T (C m ) -1 (m - m′) Formula 13 Where A represents the forward equation, i.e., Formula 8 and Formula 9; C e represents the covariance matrix of seismic noise, C m represents the covariance matrix of the inversion parameters; m′ represents the prior smooth constraint model, which is usually obtained by interpolation and extrapolation of logging data; N represents the number of sampling points.

7. The method according to claim 6, characterized in that In step six, an iterative update algorithm is used to generate a new longitudinal and transverse wave velocity ratio model. The specific method is as follows: According to Formula 1 and Formula 7, the new model m * It can be expressed as: m * =ηm cur +(1-η)m s Formula 14 Where η is the weight coefficient; m cur is the current model; for the secondary variable longitudinal and transverse wave velocity ratio, m s =m sim ; For the main variables longitudinal wave modulus and density, m s =m cosim .

8. The method according to claim 7, characterized in that In step nine, the Metropolis criterion is used to determine whether to retain the new model, including the following sub-steps: Sub-step 9.1, combining the inversion likelihood function with the prior constraint model m′, establish the following objective function: J(m) = [d - A(m)] T (C e ) -1 [d - A(m)] + (m - m′) T (C m ) -1 (m - m′) Formula 15 Where A represents the forward equation, i.e., Formula 8 and Formula 9; C e represents the covariance matrix of seismic noise, C m represents the covariance matrix of the inversion parameters; Sub-step 9.2, in order to avoid the inversion target value falling into the local extreme value, the Metropolis criterion is used to optimize the generated new model and establish the following acceptance probability expression: Where J represents the inversion objective function and T is an adjustable parameter. A random number between 0 and 1 is generated. When the acceptance probability is greater than the random number, the current optimal solution is updated. When the acceptance probability is less than the random number, the original optimal solution is retained.

Images