Pre-stack waveform inversion method based on viscoelastic layered medium

Abstract

The invention discloses a pre-stack waveform inversion method based on a viscoelastic layered medium, which relates to the field of geophysical exploration of petroleum, and comprises the following steps: carrying out forward modeling based on a viscoelastic reflection matrix method, and accurately simulating the propagation characteristics of seismic waves in the viscoelastic layered medium, thereby effectively improving the inversion prediction precision. According to the method, the elastic characteristic and the attenuation characteristic of the stratified medium are considered at the same time in the seismic wave propagation process, the P-wave speed, the S-wave speed, the density, the P-wave quality factor and the S-wave quality factor can be inverted at the same time, the precision of the inversion result is improved through fine simulation of the underground reservoir medium, and the method has the advantages of being high in practicability and the like. The multi-parameter inversion method has the advantages that strong support is provided for reservoir prediction, a Bayesian inversion framework is adopted, a multi-parameter inversion process is optimized by fusing prior information, robustness and noise immunity of the inversion method are improved, and stability and accuracy of the multi-parameter inversion method in actual seismic data application are guaranteed.

Description

Technical Field

[0001] The present invention belongs to the field of petroleum geophysical exploration, and particularly relates to a prestack waveform inversion method based on viscoelastic layered media. Background Technique

[0002] Reservoir prediction is one of the important contents of reservoir description and is the key work for carrying out the analysis of the development potential of oilfields. The underground reservoir media usually show a layered structure, and there are different elastic and attenuation characteristics between layers. These inherent attenuation effects are closely related to the properties of reservoir fluids and cannot be ignored in prestack waveform inversion. Accurately predicting the physical properties of underground reservoirs, especially the P-wave velocity, S-wave velocity, density, and quality factor Q, is the basis for effective reservoir prediction and oil and gas exploration.

[0003] Currently, in the field of oil and gas exploration, common reservoir prediction methods are mostly based on seismic inversion techniques, such as the common AVA inversion method. However, these methods mainly limit the attenuation between layers to the interface for reflection characteristic analysis and rarely consider the attenuation effect of seismic waves during propagation in the medium layer. In a complex geological environment, this treatment method of ignoring propagation attenuation will lead to a decrease in the accuracy of the inversion result and affect the accurate prediction of the reservoir. When dealing with layered media with inherent attenuation characteristics, the propagation characteristics of seismic waves are affected not only by the elastic characteristics of the medium but also by the quality factor Q of the medium. Therefore, traditional inversion methods based on elastic media often cannot accurately reflect the true properties of the medium, especially in areas with significant attenuation effects. Summary of the Invention

[0004] The purpose of the present invention is to carry out prestack waveform inversion by introducing the viscoelastic reflection matrix method, considering both the elastic and attenuation characteristics of layered media, realizing the fine simulation of underground reservoir media, improving the accuracy of the inversion result, and providing a more comprehensive reservoir description for oil and gas exploration.

[0005] To solve the above technical problems, the present invention adopts the following technical solution: A prestack waveform inversion method based on viscoelastic layered media, comprising:

[0006] S1: Conduct forward modeling based on the viscoelastic reflection matrix method, and use the plane wave theory to characterize the wave propagation law under the assumption of a layered homogeneous medium through the reflection matrix method, which is represented by the wave propagation matrix M;

[0007]

[0008] In the above formula, A is the eigenvector matrix, E is the phase shift matrix of the up-going and down-going waves in the layer, the subscripts 1 and n represent the top and bottom media respectively, the subscript k represents any intermediate layer medium, and the frequency dispersion of the formation quality factor Q for the P-wave and S-wave velocities is expressed by the following formula:

[0009]

[0010] In the above formula, ω represents frequency, V P (ω) and V S (ω) represent the complex velocities of P-wave and S-wave respectively, V P0 and V S0 represent the P-wave and S-wave velocities at the reference frequency ω c respectively, Q P and Q S represent the quality factors of P-wave and S-wave respectively. The eigenvector matrix A and phase shift matrix E are related to the P-wave velocity V P (ω) and S-wave velocity V S (ω). Substitute the complex velocities V P (ω) and V S (ω) in equations (2) and (3) into the eigenvector matrix A and phase shift matrix E in equation (1) to obtain the wave propagation matrix of the attenuating layered medium including the influence of the quality factor Q. Finally, obtain the frequency-domain reflection coefficient R(ω,θ) through the wave propagation matrix, multiply it by the frequency-domain wavelet W(ω,θ), and obtain the prestack seismic data through the inverse Fourier transform:

[0011]

[0012] In the above formula, d pp represents the angle gather in the time domain, θ represents the angle, i represents the imaginary unit, t represents time, G represents the non-linear forward operator, and m is m = [V P , V S , ρ, Q P , Q S T , where ρ is the medium density.

[0013] S2: Construct the objective function for prestack waveform inversion, and construct the objective functional that integrates Bayesian theory:

[0014]

[0015] In the above formula, J represents the objective function, d obs represents the observed seismic data, η represents the regularization parameter of the prior information, μ represents the expectation of the model, represents the inverse of the covariance matrix of the model parameter m.

[0016] Furthermore, the wave propagation matrix M in step S1 is represented by a block matrix:

[0017]

[0018] ​where the subscripts U and D represent the up-going wave and the down-going wave respectively, and the reflection and transmission coefficient matrices of the layered medium can be expressed by the block matrices M UU 、M UD 、M DU and M DD as follows:

[0019]

[0020] R represents the reflection coefficient matrix and T represents the transmission coefficient matrix.

[0021] Furthermore, the P-wave quality factor Q P and the S-wave quality factor Q S are obtained from empirical formulas:

[0022] Q P = 10.76·V P 2 ; (8)

[0023] Q S = 10.76·V S 2 . (9)

[0024] Furthermore, in step S2, prestack waveform inversion is completed in the Bayesian framework to construct the synthetic seismic record d syn :

[0025] d syn = G(m) + n; (10)

[0026] In the above formula, d syn represents the synthetic seismic record, G represents the non-linear forward operator, m represents the model parameters, and n is the random noise.

[0027] Furthermore, after the observed data d obs is given, the posterior distribution of the model parameters can be expressed as:

[0028]

[0029] In the above formula, P(d|m) is the likelihood function, P(m) is the prior distribution, and P(d) is the marginal probability density, which is a constant.

[0030] Furthermore, assuming that the noise is Gaussian noise and independent of each other, the likelihood function is expressed as:

[0031]

[0032] In the above formula, P0 is the normalization constant, is the covariance of the noise, and d obs represents the observed seismic record.

[0033] Furthermore, the prior distribution P(m) is a Gaussian distribution:

[0034]

[0035] In the above formula, μ is the expectation of the model, represents the inverse of the covariance matrix of the model parameter m, with a size of 5N * 5N, where N represents the data length of the elastic parameter.

[0036] Furthermore, the covariance matrix C m is expressed as:

[0037] C m = Kron(Cov{V P , V S , ρ, Q P , Q S}, I); (14)

[0038] In the above formula, Kron represents the Kronecker product, Cov represents the solution of the covariance matrix, I is the identity matrix, and the matrix size is N * N.

[0039] Furthermore, the posterior distribution P(m|d) is expressed as:

[0040]

[0041] In the above formula, η represents the regularization parameter of the prior information.

[0042] Furthermore, maximizing the posterior distribution P(m|d) means minimizing the objective function J(m).

[0043]

[0044] When J > error threshold or the number of iterations < maximum number of iterations, update the model parameter m, where:

[0045] Δm = -H(m k ) -1 g(m k ); (17)

[0046] m k+1 = m k + Δm; (18)

[0047] In the above formula, g and H respectively represent the gradient operator and Hessian matrix for obtaining the parameters. The updated model parameter m k+1 is reused to calculate the synthetic seismic record d syn .

[0048] Beneficial effects: A prestack waveform inversion method based on viscoelastic layered media disclosed by the present invention accurately simulates the propagation characteristics of seismic waves in attenuating layered media, and more comprehensively considers the elastic and attenuation characteristics of the media compared with traditional methods. By integrating the Bayesian inversion framework and introducing prior information, the stability and robustness of the inversion process are significantly improved. When facing multi-parameter prediction of attenuating media, it can effectively cope with the influence of the initial model and the interference of random noise, providing reliable support for reservoir prediction. Description of the Drawings

[0049] The drawings are used to provide a further understanding of the present invention, and constitute a part of the specification. They are used together with the embodiments of the present invention to explain the present invention, and do not constitute a limitation to the present invention.

[0050] In the drawings:

[0051] Figure 1 is the prestack waveform inversion flowchart of the present invention based on the viscoelastic reflection matrix method (VERMM);

[0052] Figure 2 is a schematic diagram comparing the forward modeling results of the present invention based on the Aki approximation formula, the reflectivity method in the wavenumber domain (RMM), and VERMM;

[0053] Figure 3 are the model parameters, observed data, and inversion results of the present invention;

[0054] Figure 4 are the VERMM inversion results of different initial models of the present invention;

[0055] Figure 5 are the VERMM inversion results of adding different random noises to the observed data of the present invention;

[0056] Figure 6 is the comparison of the two-dimensional inversion profile of the actual seismic data of the present invention and the inversion results of the well positions. Detailed Embodiments

[0057] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. The following text is only used to describe the implementation manner of a prestack waveform inversion method based on viscoelastic layered media of the present invention, and does not strictly limit the scope of protection specifically requested by the present invention.

[0058] In addition, the technical solutions between the various embodiments can be combined with each other, but it must be based on the ability of those of ordinary skill in the art to implement. When the combination of technical solutions results in contradictions or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection required by the present invention.

[0059] Example 1: As Figure 1 shown, the present invention discloses a prestack waveform inversion method based on viscoelastic layered media, including:

[0060] S1: Forward modeling is carried out based on the viscoelastic reflection matrix method. Using the plane wave theory, the wave propagation law under the assumption of layered homogeneous media is characterized by the reflection matrix method and represented by the wave propagation matrix M;

[0061]

[0062] In the above formula, A is the eigenvector matrix, E is the phase shift matrix of the up-going wave and down-going wave within the layer. The subscripts 1 and n represent the top and bottom layer media respectively, and the subscript k represents any intermediate layer media. The frequency dispersion of the formation quality factor Q for the P-wave and S-wave velocities is expressed by the following formula:

[0063]

[0064] In the above formula, ω represents the frequency, V P (ω) and V S (ω) represent the complex velocities of the P-wave and S-wave respectively, V P0 and V S0 represent the P-wave and S-wave velocities at the reference frequency ω c respectively, Q P and Q S represent the P-wave and S-wave quality factors respectively. The eigenvector matrix A and the phase shift matrix E are related to the P-wave velocity V P (ω), S-wave velocity V S (ω). Substitute the complex velocities V P (ω) and V S (ω) in equations (2) and (3) into the eigenvector matrix A and the phase shift matrix E in equation (1) to obtain the wave propagation matrix of the attenuating layered media including the influence of the quality factor Q. Finally, obtain the frequency-domain reflection coefficient R(ω,θ) through the wave propagation matrix, multiply it by the frequency-domain wavelet W(ω,θ), and obtain the prestack seismic data through the inverse Fourier transform:

[0065]

[0066] In the above formula, d pp represents the time-domain angle gather, θ represents the angle, i represents the imaginary unit, t represents the time, G represents the non-linear forward operator, and m is m = [V P , V S , ρ, Q P , Q S T , where ρ is the medium density.

[0067] ​S2: Construct the objective function for prestack waveform inversion and construct the objective functional that integrates Bayesian theory:

[0068]

[0069] In the above formula, J represents the objective function, d obs represents the observed seismic data, η represents the regularization parameter of the prior information, μ represents the expectation of the model, represents the inverse of the covariance matrix of the model parameter m.

[0070] In the first embodiment, the wave propagation matrix M in the step S1 is represented by a block matrix:

[0071]

[0072] where the subscripts U and D represent the up-going wave and the down-going wave respectively, and the reflection and transmission coefficient matrices of the layered medium can be represented by the block matrices M UU , M UD , M DU and M DD as follows:

[0073]

[0074] R represents the reflection coefficient matrix and T represents the transmission coefficient matrix.

[0075] In the first embodiment, the P-wave quality factor Q P and the S-wave quality factor Q S are obtained from empirical formulas:

[0076] Q P = 10.76·V P 2 ; (8)

[0077] Q S = 10.76·V S 2 ; (9)

[0078] In the first embodiment, in the step S2, prestack waveform inversion is completed under the Bayesian framework to construct the synthetic seismic record d syn :

[0079] d syn = G(m) + n; (10)

[0080] In the above formula, d syn represents the synthetic seismic record, G represents the non-linear forward operator, m represents the model parameter, and n is the random noise.

[0081] In the first embodiment, given the observed data d obsAfter that, the posterior distribution of the model parameters can be expressed as:

[0082]

[0083] In the above formula, P(d|m) is the likelihood function, P(m) is the prior distribution, and P(d) is the marginal probability density, which is a constant.

[0084] In the first embodiment, the noise is Gaussian noise and independent of each other, and the likelihood function is expressed as:

[0085]

[0086] In the formula, P0 is the normalization constant, is the covariance of the noise, and d obs represents the observed seismic record.

[0087] In the first embodiment, the prior distribution P(m) is a Gaussian distribution:

[0088]

[0089] In the above formula, μ is the expectation of the model, represents the inverse of the covariance matrix of the model parameter m, with a size of 5N*5N, and N represents the data length of the elastic parameter.

[0090] In the first embodiment, the covariance matrix C m is expressed as:

[0091] C m = Kron(Cov{V P ,V S ,ρ,Q P ,Q S},I); (14)

[0092] In the above formula, Kron represents the Kronecker product, Cov represents the solution of the covariance matrix, I is the identity matrix, and the matrix size is N*N.

[0093] In the first embodiment, the posterior distribution P(m|d) is expressed as:

[0094]

[0095] In the above formula, η represents the regularization parameter of the prior information.

[0096] In the first embodiment, maximizing the posterior distribution P(m|d) is equivalent to minimizing the objective function J(m),

[0097] That is:

[0098] When J > error threshold or the number of iterations < maximum number of iterations, update the model parameter m, where

[0099] Δm = -H(m k ) -1 g(m k ); (17)

[0100] m k+1 =m k +Δm; (18)

[0101] In the above formula, g and H respectively represent the gradient operator and Hessian matrix for obtaining parameters. By reusing the updated model parameter m k+1 to calculate the synthetic seismic record d syn .

[0102] Example 2: The present invention conducts a comprehensive model comparison test on the prestack waveform inversion method based on VERMM. The viscoelastic reflection matrix method (VERMM) is respectively used for forward modeling with the reflection matrix method (RMM) and the conventional Aki method, and VERMM inversion, RMM inversion, and AVA inversion are respectively performed, where the Aki forward modeling corresponds to the AVA inversion.

[0103] As Figure 2 shown in a, a model of physical parameters including the longitudinal wave velocity V p , the transverse wave velocity V s , density, the longitudinal wave attenuation Q p and the transverse wave attenuation Q S is constructed. As Figure 2 shown in b, based on the longitudinal wave velocity V p , the transverse wave velocity V s , density, the longitudinal wave attenuation Q p and the transverse wave attenuation Q S , the forward modeling result is obtained. As Figure 3 shown in c, the prestack waveform inversion result based on VERMM of the present invention is compared with the traditional AVA inversion result and RMM inversion result. When facing an attenuating medium and when the medium quality factor is relatively high, the inversion method based on VERMM can more accurately predict the elastic parameters of the reservoir. In particular, it can invert the quality factor Q. The inversion results of the traditional AVA inversion and RMM inversion have relatively large errors in the region with a relatively high medium quality factor Q. The inversion result based on VERMM significantly reduces the errors of the traditional methods in the attenuating medium.

[0104] Example 3: As Figure 4As shown in the figure, in order to further verify the stability and robustness of the method, we tested the influence on the inversion results by changing the smoothness of the initial model and setting the smooth window sizes of smooth-50, smooth-70, and smooth-100 respectively. The inversion results of VERMM have a higher degree of fit with the true model, indicating that the prestack waveform inversion method based on VERMM still maintains high accuracy under different sizes of smooth windows, showing that the prestack waveform inversion method based on VERMM has good robustness to the initial model.

[0105] Example 4: As Figure 5 shown in the figure, the influence test on the inversion results was carried out by simulating different degrees of random noise. The noise intensity was measured by the signal-to-noise ratio (SNR). Noises with SNR = 10 and SNR = 5 were added to the observed data respectively. Although the noise was increased, the inversion method based on VERMM could still effectively predict the model parameters, demonstrating that the method has a certain anti-noise ability. The model data test fully proves the superiority of the prestack waveform inversion method based on VERMM in the attenuating medium, and it has strong stability, robustness, and anti-noise ability, providing strong support for practical applications.

[0106] Example 5: As Figure 6 shown in Figure a, the prestack waveform inversion method of the present invention was applied to the actual seismic data of an oil and gas exploration area. After data processing and inversion, the P-wave and S-wave quality factor Q profiles of this area were successfully predicted. The inversion results not only provided an accurate estimate of the elastic parameters but also provided a more detailed description of the attribute characteristics of the underground reservoir. As Figure 6 shown in Figure b, the comparison between the inversion results at Well A and the well curve data was presented. The inverted P-wave and S-wave quality factor Q profiles were consistent with the Q-value data measured by the well curve.

[0107] The actual application shows that the prestack waveform inversion method based on VERMM can be successfully applied in the real attenuating medium and provide high-precision reservoir prediction results, providing important technical support for further oil and gas reservoir characterization and development decision-making.

[0108] The above only describes the present invention and its implementation manners. This description is not restrictive. What is shown in the drawings is only one of the implementation manners of the present invention. The actual structure is not limited thereto. Therefore, if those of ordinary skill in the art are inspired by it and design similar structural modes and embodiments without creative efforts without departing from the spirit of the present invention, they shall fall within the protection scope of the present invention.

Claims

1. A prestack waveform inversion method based on viscoelastic layered media, characterized in that, Including: S1: Forward modeling is carried out based on the viscoelastic reflection matrix method. The wave propagation law under the assumption of a layered homogeneous medium is characterized by the reflection matrix method using the plane wave theory. The wave propagation matrix M is expressed as: In the above formula, A(V P ,V S ,ρ,Q P ,Q S ) is the eigenvector matrix, and E(V P ,V S ,Q P ,Q S ) is the phase shift matrix of the up-going wave and down-going wave within the layer, where V p and V s represent the longitudinal wave velocity and the shear wave velocity respectively, Q p and Q s represent the quality factor of the P-wave and the quality factor of the S-wave respectively, where ρ is the medium density, the subscripts 1 and n represent the top layer and bottom layer media respectively, the subscript k represents any intermediate layer medium, and the frequency dispersion of the formation quality factor Q with respect to the P-wave and S-wave velocities is expressed by the following formula: In the above formula, ω represents frequency, V P (ω) and V S (ω) represent the complex velocities of P-wave and S-wave respectively, V P0 and V S0 represent the P-wave and S-wave velocities at the reference frequency ω c respectively. Substitute V P (ω) and V S (ω) into Equation (1) to obtain the wave propagation matrix that includes the influence of the quality factor Q. Construct the frequency-domain reflection coefficient R(ω, θ) through the wave propagation matrix, multiply R(ω, θ) by the frequency-domain wavelet W(ω, θ), and obtain the prestack seismic data through the inverse Fourier transform: In the above formula, d pp represents the angle gather in the time domain, θ represents the angle, i represents the imaginary unit, t represents the time, G represents the nonlinear forward operator, and m = [V P , V S , ρ, Q P , Q S T ;​ S2: Construction of the objective function for prestack waveform inversion, constructing an objective functional that integrates Bayesian theory: In the above formula, J represents the objective function, d obs represents the observed seismic data, η represents the regularization parameter of the prior information, μ represents the expectation of the model, represents the inverse of the covariance matrix of the model parameter m.

2. The prestack waveform inversion method based on a viscoelastic layered medium according to claim 1, wherein: The wave propagation matrix M described in step S1 is represented by a block matrix: where the subscripts U and D represent the up-going wave and the down-going wave respectively, and the reflection and transmission coefficient matrix of the layered medium can be expressed by the block matrices M UU , M UD , M DU and M DD as follows: R represents the reflection coefficient matrix, and T represents the transmission coefficient matrix.

3. The prestack waveform inversion method based on a viscoelastic layered medium according to claim 1, wherein: P-wave quality factor Q P and S-wave quality factor Q S obtained from empirical formulas: Q P = 10.76·V P 2 ; (8) Q S = 10.76·V S 2 。 (9) 4. A prestack waveform inversion method based on a viscoelastic layered medium according to claim 1, characterized in that: In step S2, prestack waveform inversion is completed under the Bayesian framework to construct a synthetic seismic record d syn : d syn = G(m) + n; (10) In the above formula, G represents the nonlinear forward operator, m represents the model parameter, and n is the random noise.

5. A prestack waveform inversion method based on a viscoelastic layered medium according to claim 1, characterized in that: The posterior distribution of the model parameter m is expressed as: In the above formula, P(d|m) is the likelihood function, P(m) is the prior distribution, and P(d) is the marginal probability density, which is a constant.

6. A prestack waveform inversion method based on a viscoelastic layered medium according to claim 5, characterized in that: The noise is Gaussian noise and independent of each other. The likelihood function is expressed as: In the above formula, P0 is the normalization constant, is the covariance of the noise, and d obs represents the observed seismic record.

7. A prestack waveform inversion method based on a viscoelastic layered medium according to claim 5, wherein: The prior distribution P(m) is a Gaussian distribution: In the above formula, μ is the expectation of the model, represents the inverse of the covariance matrix of the model parameter m, with a size of 5N * 5N, where N represents the data length of the elastic parameter.

8. The prestack waveform inversion method based on a viscoelastic layered medium according to claim 7, characterized in that: Covariance matrix C m denote C m = Kron(Cov{V P ,V S ,ρ,Q P ,Q S},I); (14) In the above formula, Kron represents the Kronecker product, Cov represents the solution of the covariance matrix, I is the identity matrix, and the matrix size is N*N.

9. A prestack waveform inversion method based on a viscoelastic layered medium according to claim 8, characterized in that: The posterior distribution P(m|d) is expressed as: In the above formula, η represents the regularization parameter of the prior information.

10. A prestack waveform inversion method based on a viscoelastic layered medium according to claim 9, characterized in that: Maximizing the posterior distribution P(m|d) means minimizing the objective function J(m). When J > error threshold or the number of iterations < maximum number of iterations, update the model parameter m, where: Δm = -H(m k ) -1 g(m k ); (17) m k+1 = m k + Δm; (18) In the above formula, g and H represent the gradient operator and Hessian matrix for obtaining parameters respectively, and the updated model parameter m k+1 is reused to calculate the synthetic seismic record d syn .

Images