VTI medium anisotropy parameter pre-stack seismic synchronous inversion method and system

Abstract

The application discloses a VTI medium anisotropy parameter prestack seismic synchronous inversion method and system. The method comprises the following steps: constructing a target functional based on simultaneous inversion of VTI medium anisotropy coupling parameters and reflectivity of elastic parameters; inputting the angle gather seismic data, seismic wavelet and prior model into the target functional to calculate the reflectivity of the four parameters under different regularization conditions; estimating the absolute values of the VTI medium anisotropy coupling parameters and the elastic parameters according to the reflectivity of the four parameters, and then calculating the VTI medium anisotropy parameters. The application directly estimates the anisotropy parameters of the VTI medium through the prestack synchronous inversion results of the four parameters, effectively simplifying the inversion calculation process; meanwhile, the absolute values of the anisotropy coupling parameters and the elastic parameters are estimated in the frequency domain, effectively guaranteeing the stability and fidelity of the estimated anisotropy coupling parameters and the elastic parameters and the finally estimated anisotropy parameters.

Description

Technical Field

[0001] The present invention belongs to the technical field of pre-stack seismic inversion for seismic exploration, and more particularly, relates to a pre-stack seismic synchronous inversion method and system for anisotropic parameters of VTI media. Background Art

[0002] Prestack AVO inversion of isotropic media can theoretically directly determine multiple elastic parameters closely related to lithology and even hydrocarbon content, such as P-wave and S-wave velocities and density. However, it is not suitable for inverting elastic parameters in subsurface formations with low-angle natural fractures such as bedding fractures and horizontal fractures. Therefore, the VTI medium reflection coefficient equation is used to invert elastic and anisotropic parameters. Although the VTI medium reflection and transmission approximate formula derived by Rüger (1997, 1998) is widely used, this approximate formula requires simultaneous inversion of five unknown parameters, and the coefficient of the anisotropic parameter term is relatively small, making direct simultaneous inversion highly ill-posed.

[0003] To this end, Zhang (2019) proposed a simplified Rüger approximation formula with three terms. The elastic parameters related to the anisotropic parameters were obtained using isotropic and anisotropic split-step prestack inversion methods, respectively. The anisotropic parameters for predicting low-angle natural fractures were then estimated. This method effectively reduced the ill-posedness of the five-parameter direct inversion and achieved an estimation of the Thomsen anisotropy parameter ε, which is significantly correlated with the development of low-angle fractures. However, this method uses a split-step inversion, which is equivalent to performing two independent rounds of AVO prestack three-parameter inversion: an isotropic three-parameter low-angle prestack AVO inversion and a VTI medium three-parameter prestack AVO inversion. This increases the complexity, computational effort, and uncertainty of the inversion workflow. Any inversion round that is not performed well may affect the reliability of the final anisotropic parameter ε estimate. In particular, the stability and accuracy of the inversion results of the isotropic three-parameter low-angle prestack AVO inversion are subject to significant uncertainty due to the lack of medium and large angle information. Furthermore, this method uses the reflection coefficient integration method to calculate the absolute values ​​of the reflectivities of various elastic parameters obtained from reflection, which is subject to error accumulation and computational distortion. Therefore, for seismic data with a low signal-to-noise ratio, a large regularization factor is often used to ensure the stability of the inversion results. This significantly reduces the resolution of the inversion results and loses much detailed information.

[0004] Therefore, how to simplify the inversion workflow of anisotropic parameters of VTI media and improve the stability and precision of the calculation results of the absolute values ​​of elastic parameters and anisotropic parameters is of great practical significance.

[0005] The information disclosed in the background technology section of the present invention is only intended to deepen the understanding of the general background technology of the present invention, and should not be regarded as an admission or any form of suggestion that the information constitutes the prior art already known to those skilled in the art. Summary of the Invention

[0006] The purpose of the present invention is to propose a method for directly estimating the anisotropic parameters of VTI media through the results of four-parameter prestack synchronous inversion, which effectively simplifies the inversion workflow of the anisotropic parameters of VTI media and improves the stability and fidelity of the calculation results of the absolute values ​​of elastic parameters and anisotropic parameters.

[0007] To achieve the above objectives, the present invention proposes a pre-stack seismic synchronous inversion method and system for anisotropic parameters of VTI media.

[0008] According to a first aspect of the present invention, a method for prestack seismic synchronous inversion of anisotropic parameters of VTI media is proposed, comprising:

[0009] Obtaining a priori models of reflectivity from angle gather seismic data, seismic wavelets, and VTI medium anisotropic coupling parameters and elastic parameters;

[0010] Construct a target functional for the simultaneous inversion of the reflectivity of the anisotropic coupling parameters and elastic parameters of VTI media based on angle gather seismic data;

[0011] Inputting the angle gather seismic data, the seismic wavelet and the prior model into the target functional, and calculating the reflectivity of the four anisotropic coupling parameters and elastic parameters of the VTI medium under different regularization conditions through an optimization algorithm;

[0012] estimating the absolute values ​​of the anisotropic coupling parameter and the elastic parameter of the VTI medium based on the conversion relationship between the reflectivity of the anisotropic coupling parameter and the elastic parameter and the absolute value thereof in the frequency domain and the reflectivity;

[0013] The anisotropic parameters of the VTI medium are calculated according to the coupling relationship between the anisotropic coupling parameters and the elastic parameters and the absolute values.

[0014] Optionally, the target functional for inverting the reflectivity of the anisotropic coupling parameters and elastic parameters of the VTI medium specifically includes:

[0015] Based on the theoretical calculation approximate formula of the angle domain reflection coefficient of VTI media and the convolution model, a group of linear equations for the angle gather seismic data and the anisotropic coupling parameter reflectivity and elastic parameter reflectivity are established. Combined with the prior model of the anisotropic coupling parameter reflectivity and elastic parameter reflectivity of the VTI medium, a target functional for simultaneously inverting the reflectivity of the underground anisotropic coupling parameters and elastic parameters based on the angle gather seismic data is constructed.

[0016] Optionally, it is characterized in that the expression of the theoretical calculation approximate formula of the angle domain reflection coefficient of the VTI medium is:

[0017] R VTI (θ)=r1+tan 2 θr2-8gsin 2 θr3+(1-4gsin 2 θ)r4;

[0018] Among them, R VTI (θ) is the angle domain reflection coefficient of the VTI medium, θ is the incident angle of the seismic wave, r1 is the reflectivity of the elastic coefficient A1, A1=V P0 , r2 is the reflectivity of the anisotropic coupling parameter A2, A2=V P0 e ε , r3 anisotropic coupling parameter A3 reflectivity, r4 is the reflectivity of the elastic coefficient A4, A4=ρ,V P0 is the longitudinal wave velocity, V s0 is the shear wave velocity, ρ is the density, Δ is the difference in medium parameters above and below the formation interface, ε and δ are the Thomsen anisotropy parameters of the VTI medium, g is the square of the ratio of the average shear wave to the average compressional wave velocity,

[0019] Optionally, the linear equations are expressed as:

[0020] d(t,θ)=w(t)*R VTI (t,θ)=w(t)*(r1(t)+tan 2 θr2(t)-8gsin 2 θr3(t)+(1-4gsin 2 θ)r4(t));

[0021] Among them, R VTI (t,θ) is the angle domain reflection coefficient, d(t,θ) is the angle gather seismic data, t is the vertical two-way travel time of the seismic data, w(t) is the seismic wavelet, and * is the convolution operation.

[0022] Optionally, the target functional is:

[0023] J=||WLr-d||2+μ||D x (r-r0)|| p ;

[0024] Assume that there are M discrete angle gather seismic data channels, each with N sample points. The angle gather seismic data constitute a one-dimensional vector d of length N·M. The four anisotropic coupling parameters and elastic parameters to be inverted, namely, the reflectivities r1(t), r2(t), r3(t), and r4(t), constitute a one-dimensional vector r of length 4N, where r = [r1 r2 r3 r4] T , J is the target functional for simultaneously inverting the anisotropic coupling parameters and elastic parameter reflectivity of VTI media based on angle gather seismic data, W is the N·M×N·M wavelet matrix composed of seismic wavelets w(t), L is the coefficient of the four reflectivities 1, tan 2 θ, -8gsin 2 θ, 1-4gsin 2 θ forms a coefficient matrix of N·M×4N, p is the norm constraint of the reflectivity regularization term, μ is the regularization parameter, μ>0, D x is the regularization operator, x is the differential order, and r0 is the prior model of r.

[0025] Optionally, the different regularization conditions include:

[0026] p is set to 2 to represent the reflectivity regularization term as a 2-norm constraint, and a smooth reflectivity r solution is obtained;

[0027] Taking p as 1 represents that the reflectivity regularization term is constrained by 1 norm, and a sparse edge-preserving reflectivity r solution is obtained.

[0028] Optionally, the calculating of the reflections of the four anisotropic coupling parameters and elastic parameters of the VTI medium under different regularization conditions by using an optimization algorithm specifically includes:

[0029] When p=2, the target functional J is minimized and the damped least squares method is used to solve the reflectivity r of the four anisotropic coupling parameters and the elastic parameter:

[0030]

[0031] When p=1, the target functional J is minimized and the iterative reweighted least squares method is used to solve the reflectivity r of the four anisotropic coupling parameters and the elastic parameter:

[0032]

[0033] Among them, k is the kth iteration of the setting, it max The maximum number of iterations set for iterative solution, Q k is the reweighting matrix of the kth iteration, η is a small positive constant, r' is the reflectivity without low-frequency components, r' obtained at the kth iteration (k)The jth element of the vector after the differential operation is summed with η and the reciprocal is taken as Q k The jth element of the diagonal, L T is the transposed matrix of L, W T is the transposed matrix of W, D x T D x The transposed matrix of .

[0034] Optionally, the expression for estimating the absolute value of the anisotropic coupling parameter and the elastic parameter of the VTI medium is:

[0035]

[0036] Where i = 1, 2, 3 or 4, are the absolute values ​​of the anisotropic coupling parameters and elastic parameters within the estimated effective frequency band b(f), b(f) is the filter with an asymmetric structure, FT -1 (·) represents the inverse Fourier transform, i is the imaginary unit, A i0 A i (t), a is the regularization factor, which is a small value greater than zero, e is a natural constant, T(f) is the Tikhonov smoothing operator as a regularization term, T(f) = 1 / |f| n , n>1, r i ′(f) is the reflectivity spectrum of the i-th parameter in the four-parameter reflectivity without low-frequency components.

[0037] Optionally, the expression for calculating the anisotropy parameter of the VTI medium is:

[0038]

[0039] in, is the final estimated anisotropy parameter of the VTI medium, that is, the estimated value of ε, is the absolute value of the first parameter within the estimated effective frequency band b(f), that is, V P0 The estimated value of is the absolute value of the second parameter within the estimated effective frequency band b(f), that is, V P0 e ε Estimated value of .

[0040] According to a second aspect of the present invention, a system for pre-stack seismic synchronous inversion of anisotropic parameters of VTI media is proposed, which is used to execute the pre-stack seismic synchronous inversion method of anisotropic parameters of VTI media according to any one of the first aspects, comprising:

[0041] Acquisition module, used to obtain the angle gather seismic data, seismic wavelet and reflection prior model of VTI medium anisotropic coupling parameters and elastic parameters;

[0042] A construction module is used to construct a target functional for simultaneously inverting the anisotropic coupling parameters and elastic parameter reflectivity of VTI media based on angle gather seismic data;

[0043] An input and calculation module, configured to input the angle gather seismic data, the seismic wavelet, and the prior model into the target functional, and calculate the reflectivity of the four anisotropic coupling parameters and elastic parameters of the VTI medium under different regularization conditions through an optimization algorithm;

[0044] an estimation module, configured to estimate the absolute values ​​of the anisotropic coupling parameter and the elastic parameter of the VTI medium based on a conversion relationship between the reflectivity of the anisotropic coupling parameter and the elastic parameter and their absolute values ​​in the frequency domain and the reflectivity;

[0045] The calculation module is used to calculate the anisotropic parameters of the VTI medium according to the coupling relationship between the anisotropic coupling parameters and the elastic parameters and the absolute value.

[0046] The beneficial effects of the present invention are as follows: the present invention directly estimates the anisotropic parameters of the VTI medium through the pre-stack synchronous inversion results of four parameters (two elastic parameters and two anisotropic coupling parameters), without the need to sort small-angle and medium-large angle gathers, and performs small-angle pre-stack isotropic three-parameter inversion and medium-large angle pre-stack anisotropic three-parameter inversion in steps, effectively simplifying the inversion calculation process and reducing the uncertainty and computational complexity of the step-by-step inversion process. At the same time, the present invention uses a frequency domain method to estimate the absolute values ​​of the anisotropic coupling parameters and elastic parameters, avoiding the error accumulation and calculation distortion problems caused by the traditional reflectivity integration method, eliminating the influence of the distorted low-frequency and high-frequency information components, and effectively ensuring the stability and fidelity of the estimated anisotropic coupling parameters and elastic parameters and their final estimated anisotropic parameters. In addition, the present invention only requires the inversion of four parameters, which is fewer than the five-parameter inversion parameters to be determined, further reducing the inversion ambiguity.

[0047] The system of the present invention has other features and advantages that will be apparent from or will be described in detail in the accompanying drawings and subsequent detailed description incorporated herein, which together serve to explain the specific principles of the invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] The above and other objects, features and advantages of the present invention will become more apparent through a more detailed description of exemplary embodiments of the present invention with reference to the accompanying drawings, in which like reference numerals generally represent like components.

[0049] Figure 1 A flow chart showing the steps of the pre-stack seismic synchronous inversion method for VTI medium anisotropy parameters according to the present invention.

[0050] Figure 2 The anisotropy parameters δ and ε, the longitudinal wave velocity V according to Example 1 of the present invention are shown. p0 , shear wave velocity V S0 Schematic diagram of the theoretical curves of and density ρ.

[0051] Figure 3 a and Figure 3 b shows a schematic diagram of angle gather seismic data according to Example 1 of the present invention.

[0052] Figure 4 A schematic diagram showing the steps of estimating anisotropic parameters of a VTI medium according to embodiment 1 of the present invention.

[0053] Figure 5 a and Figure 5 b is a schematic diagram showing the comparison results of anisotropic parameters of angle gather seismic data inversion according to Example 1 of the present invention with the real model and the prior model.

[0054] Figure 6 A schematic diagram of a pre-stack seismic synchronous inversion system for anisotropic parameters of VTI media according to embodiment 2 of the present invention is shown. DETAILED DESCRIPTION

[0055] The present invention will now be described in more detail with reference to the accompanying drawings. While preferred embodiments of the present invention are shown in the accompanying drawings, it should be understood that the present invention may be implemented in various forms and is not limited to the embodiments set forth herein. Rather, these embodiments are provided to make the present invention more thorough and complete and to fully convey the scope of the present invention to those skilled in the art.

[0056] like Figure 1 As shown, a pre-stack seismic synchronous inversion method for anisotropic parameters of VTI media according to the present invention includes:

[0057] Obtaining a priori models of reflectivity from angle gather seismic data, seismic wavelets, and VTI medium anisotropic coupling parameters and elastic parameters;

[0058] Construct a target functional for the simultaneous inversion of the reflectivity of the anisotropic coupling parameters and elastic parameters of VTI media based on angle gather seismic data;

[0059] The angle gather seismic data, seismic wavelet and prior model are input into the target functional, and the reflectivity of four VTI media anisotropic coupling parameters and elastic parameters under different regularization conditions is calculated by the optimization algorithm.

[0060] The absolute values ​​of anisotropic coupling parameters and elastic parameters of VTI medium are estimated based on the conversion relationship between the reflectivity of anisotropic coupling parameters and elastic parameters and their absolute values ​​in the frequency domain and the reflectivity;

[0061] The anisotropic parameters of VTI media are calculated based on the coupling relationship and absolute values ​​between anisotropic coupling parameters and elastic parameters.

[0062] Specifically, the present invention uses a theoretical calculation approximate formula between the angle-domain reflection coefficient of a VTI medium and the reflectivity of anisotropic coupling parameters and elastic parameters, and establishes a linear equation group between angle gather seismic data and the reflectivity of anisotropic coupling parameters and elastic parameters based on a convolution model. Under the constraints of an a priori reflection model of the anisotropic coupling parameters and elastic parameters of the VTI medium, a target functional is constructed to simultaneously invert the reflectivity of the anisotropic coupling parameters and elastic parameters of the underground VTI medium. Then, the absolute values ​​of the anisotropic coupling parameters and elastic parameters of the VTI medium are estimated based on the conversion relationship between the reflectivity of the anisotropic coupling parameters and elastic parameters and their absolute values ​​in the frequency domain and the reflectivity. Finally, the absolute values ​​of the anisotropic coupling parameters and elastic parameters of the underground VTI medium are estimated based on the conversion relationship between the reflectivity of the anisotropic coupling parameters and elastic parameters and their absolute values ​​in the frequency domain, and the anisotropic parameters of the VTI medium are calculated based on the coupling relationship between the anisotropic coupling parameters and elastic parameters.

[0063] In one example, constructing a target functional for inverting the anisotropic coupling parameters and elastic parameter reflectivity of a VTI medium specifically includes:

[0064] Based on the approximate formula for the theoretical calculation of the angle-domain reflection coefficient of VTI media and the convolution model, a set of linear equations for the angle gather seismic data and the anisotropic coupling parameter reflectivity and elastic parameter reflectivity are established. Combined with the prior models of the anisotropic coupling parameter reflectivity and elastic parameter reflectivity of VTI media, a target functional is constructed for simultaneously inverting the reflectivity of the underground anisotropic coupling parameter and elastic parameter based on the angle gather seismic data.

[0065] Specifically, the theoretical calculation approximate formula of the angle domain reflection coefficient of VTI media and the theoretical calculation approximate formula between the anisotropic coupling parameter reflectivity and the elastic parameter reflectivity are used, and a group of linear equations between the angle gather seismic data and the anisotropic coupling parameter reflectivity and the elastic parameter reflectivity are established based on the convolution model. Then, under the prior model constraints of the anisotropic coupling parameter reflectivity and the elastic parameter reflectivity of the VTI medium, a target functional is constructed to simultaneously invert the reflectivity of the anisotropic coupling parameters and elastic parameters of the underground VTI medium.

[0066] In one example, the theoretical approximate formula for calculating the angle domain reflection coefficient of a VTI medium is expressed as follows:

[0067] R VTI (θ)=r1+tan 2 θr2-8gsin 2 θr3+(1-4gsin 2 θ)r4;

[0068] Among them, R VTI (θ) is the angle domain reflection coefficient, θ is the incident angle of the seismic wave, r1 is the reflectivity of the elastic coefficient A1, A1=V P0 , r2 is the reflectivity of the anisotropic coupling parameter A2, A2=V P0 e ε , r3 anisotropic coupling parameter A3 reflectivity, r4 is the reflectivity of the elastic coefficient A4, A4=ρ,V P0 is the longitudinal wave velocity, V s0 is the shear wave velocity, ρ is the density, Δ is the difference in medium parameters above and below the formation interface (lower layer minus upper layer), ε and δ are the Thomsen anisotropy parameters of the VTI medium, e is a natural constant, g is the square of the ratio of the average shear wave velocity to the average compressional wave velocity,

[0069] Specifically, the Aki-Richards isotropic medium angle domain reflection coefficient approximate formula is used to replace the isotropic term in the VTI medium angle domain reflection coefficient approximate formula proposed by Rüger, and the trigonometric function relationship sin is applied. 2 θtan 2 θ=tan 2 θ-sin 2 θ, sec 2 θ=tan 2 θ+1, we can obtain the new approximate formula of the VTI medium angle domain reflection coefficient:

[0070]

[0071] In order to construct the theoretical calculation approximate formula between the angle domain reflection coefficient of VTI medium and the anisotropic coupling parameter and elastic parameter reflectivity, further approximation is made to each reflectivity term.

[0072] because

[0073]

[0074] Assuming that the elastic parameters and anisotropic parameters of adjacent strata in the VTI medium are continuously differentiable over time and the differences are small, we have Furthermore, there Then formula (2) can be rewritten as:

[0075]

[0076] Similarly, we can get:

[0077]

[0078] Combining equations (3) and (4), equation (1) can be rewritten as:

[0079]

[0080] make

[0081] Then we have:

[0082] R VTI (θ)=r1+tan 2 θr2-8gsin 2 θr3+(1-4gsin 2 θ)r4 (7);

[0083] It can be seen from formula (7) that the formula has a standard reflectivity form. The reflection coefficient of the VTI medium is composed of four parameters (elastic parameter A1 = V P0 , A4 = ρ; anisotropic coupling parameter A2 = V P0 e ε ,

[0084] ) reflectivity r i (i=1,2,3,4) weighted sum is obtained.

[0085] In one example, the convolution model is expressed as:

[0086] d(t,θ)=w(t)*R VTI (t,θ)=w(t)*(r1(t)+tan 2 θr2(t)-8gsin 2 θr3(t)+(1-4gsin2 θ)r4(t));

[0087] Among them, R VTI (t,θ) is the angle domain reflection coefficient, d(t,θ) is the angle gather seismic data, t is the vertical two-way travel time of the seismic data, w(t) is the seismic wavelet, and * is the convolution operation.

[0088] Specifically, in the time domain, assuming that the angle domain reflection coefficient R VTI (t,θ) and the angle gather seismic data d(t,θ) satisfy the convolution model relationship, then:

[0089] d(t,θ)=w(t)*R VTI (t,θ)=w(t)*(r1(t)+tan 2 θr2(t)-8gsin 2 θr3(t)+(1-4gsin 2 θ)r4(t)) (8);

[0090] In one example, the target functional is:

[0091] J=||WLr-d||2+μ||D x (r-r0)|| p ;

[0092] Assume that there are M discrete angle gather seismic data channels, each with N sample points. The angle gather seismic data constitute a one-dimensional vector d of length N·M. The four anisotropic coupling parameters and elastic parameters to be inverted, namely, the reflectivities r1(t), r2(t), r3(t), and r4(t), constitute a one-dimensional vector r of length 4N, where r = [r1 r2 r3 r4] T , J is the target functional for simultaneously inverting the anisotropic coupling parameters and elastic parameter reflectivity of VTI media based on angle gather seismic data, W is the N·M×N·M wavelet matrix composed of seismic wavelets w(t), L is the coefficient of the four reflectivities 1, tan 2 θ, -8gsin 2 θ, 1-4gsin 2 θ forms a coefficient matrix of N·M×4N, p is the norm constraint of the reflectivity regularization term, μ is the regularization parameter, μ>0, D x is the regularization operator, x is the differential order, and r0 is the prior model of r.

[0093] Specifically, suppose there are M channels of discrete angle gather seismic data, and each channel has N sample points. Then, the angle gather seismic data can form a one-dimensional vector d with a length of N·M, that is, the i-th angle gather seismic data sample point of the j-th channel is recorded as d(t i ,θ j), i = 1, 2, ..., N, j = 1, 2, ..., M; the four reflectivities to be inverted (r1(t), r2(t), r3(t), r4(t)) can form a one-dimensional vector r = [r1 r2 r3 r4] with a length of 4N T Its prior model is r0, which is usually calculated using A from logging data. i (t) Smoothed low-frequency model A i0 Calculated according to formula (6); the coefficients of the four reflectances (1, tan 2 θ, -8gsin 2 θ, 1-4gsin 2 θ) can form an N·M×4N coefficient matrix L; the seismic wavelet w(t) can form an N·M×N·M wavelet matrix W;

[0094] Based on the linear inversion theory, the target functional J for inverting the anisotropic coupling parameters and elastic parameter reflectivity r of VTI media can be established:

[0095] J=||WLr-d || 2+μ||D x (r-r0)|| p (9);

[0096] Where μ is the regularization parameter, μ>0, which is related to the noise level of the angle gather seismic data d. The heavier the noise, the larger μ is to ensure the stability of the solution; D x is the regularization operator, usually a differential operator, x is the differential order, usually 0, 1 or 2, x 0 represents D x is the unit matrix I, x takes 1 to represent D x is the first-order differential operator D1; x is 2 to represent D x is the second-order differential operator D2;

[0097]

[0098] In one example, different regularization conditions include:

[0099] p is set to 2 to represent the reflectivity regularization term as a 2-norm constraint, and a smooth reflectivity r solution is obtained;

[0100] Taking p as 1 represents that the reflectivity regularization term is constrained by 1 norm, and a sparse edge-preserving reflectivity r solution is obtained.

[0101] Specifically, p is 2, which means that the reflectivity regularization term is constrained by the 2-norm, and a smooth reflectivity r solution can be obtained. p is 1, which means that the reflectivity regularization term is constrained by the 1-norm, and a sparse edge-preserving reflectivity r solution can be obtained. The p value can be selected according to actual needs.

[0102] In one example, the reflectivities of four VTI media anisotropic coupling parameters and elastic parameters under different regularization conditions include:

[0103] When p = 2, the target functional J is minimized and the damped least squares method is used to solve the reflectivity r of the four anisotropic coupling parameters and the elastic parameter:

[0104]

[0105] When p = 1, the target functional J is minimized and the iterative reweighted least squares method is used to solve the reflectivity r of the four anisotropic coupling parameters and the elastic parameter:

[0106]

[0107] Among them, k represents the kth iteration, it max The maximum number of iterations set for iterative solution, Q k is the reweighting matrix of the kth iteration, η is a small positive constant, r' is the reflectivity without low-frequency components, r' obtained at the kth iteration (k) The jth element of the vector after the differential operation is summed with η and the reciprocal is taken as Q k The jth element of the diagonal, L T is the transposed matrix of L, W T is the transposed matrix of W, D x T D x The transposed matrix of .

[0108] Specifically, when p = 2, the target functional J is minimized and the damped least squares method is used to solve the four reflectance r items:

[0109]

[0110] When p = 1, the target functional J is minimized and the iterative reweighted least squares method is used to solve the four reflectance r:

[0111]

[0112] Among them, k represents the kth iteration, it max The maximum number of iterations set for iterative solution, Q k is the reweighting matrix of the kth iteration, and η is a small positive constant.

[0113] According to Equation (10) or Equation (11), the anisotropic coupling parameters of the VTI medium and the four reflectivities r or r′ (r′ generally has no low-frequency components) of the elastic parameters under two different regularization conditions can be inverted. The optimization algorithms used are the damped least squares method and the iterative reweighted least squares method, respectively. Other optimization algorithms can be selected according to actual conditions.

[0114] In one example, the expressions for estimating the absolute values ​​of the anisotropic coupling parameters and elastic parameters of the VTI medium are:

[0115]

[0116] Where i = 1, 2, 3 or 4, are the absolute values ​​of the anisotropic coupling parameters and elastic parameters within the estimated effective frequency band b(f), b(f) is the filter with an asymmetric structure, FT -1 (·) represents the inverse Fourier transform, i is the imaginary unit, A i0 A i (t), a is the regularization factor, which is a small value greater than zero, e is a natural constant, T(f) is the Tikhonov smoothing operator as a regularization term, T(f) = 1 / |f| n , r i ′(f) is the reflectivity spectrum of the i-th parameter in the four-parameter reflectivity without low-frequency components.

[0117] Specifically, according to the reflectivity trace integral theory, the anisotropic coupling parameter, the absolute value of the elastic parameter and its reflectivity in the present invention satisfy the following relationship:

[0118]

[0119] Where i = 1, 2, 3, 4;

[0120] Obviously, when the reflectivity r i In the presence of noise or systematic errors, The existence of error accumulation will cause instability or distortion in the absolute value calculated by equation (12).

[0121] Assuming the formation anisotropic coupling parameters and elastic parameters A i (t) is continuously differentiable with the two-way propagation time t of the seismic wave and the numerical difference between two adjacent time points is small. The four parameters reflectivity r can be derived by taking the natural logarithm derivative of formula (11) i The approximate relationship is:

[0122]

[0123] Performing Fourier transform on equation (12), we have:

[0124]

[0125] Where FT(·) represents Fourier transform, is the relative value of anisotropic coupling parameter and elastic parameter, i is the imaginary unit, A i0 Usually A i (t) average value, when A i (t) When low-frequency information is missing or missing a lot, the A calculated using logging data is usually i (t) Smoothed low-frequency model A i0 (t), which corresponds to the reflection prior model r0 in equation (9).

[0126] Considering that when f = 0, there are singular values ​​in Equation (14), therefore, the regularization parameter I(f) = 1 / |f| is introduced. p ,

[0127]

[0128] Then we can introduce the Tikhonov smoothing operator T(f)=1 / |f| n As a regularization term, in order to ensure that Equation (3) is stable in the extreme case f→0, n should be greater than 1. Here, n can be 2. a is the regularization factor, which is a small value greater than zero. Although lnA is obtained by Equation (15), i (t) spectrum, but for noisy seismic data, the inverted r i (f) Below the lowest cutoff frequency f l and above the highest cutoff frequency f h The spectrum components are also unreliable. The frequency boosting process will amplify the influence of high and low frequency noise. -1 (·) obtained ln[A i (t) / A i0 ] will be contaminated. Therefore, filtering operation needs to be added before the inverse operation, and data within the effective frequency band should be used to estimate lnA. i (t) / A i0 The Tukey window function can achieve this goal better, so the asymmetric filter b(f) is designed as follows:

[0129]

[0130] Among them, the frequency band range of the filter is B = f h -f l, λ1,λ2∈(0,1) is the ratio of the cosine rim length of the bandpass filter to the filter length. Usually, since the low- and medium-frequency information of seismic inversion is relatively reliable, the value of λ1 is small and smaller than λ2.

[0131] Assume that the relative values ​​of anisotropic coupling parameters and elastic parameters after b(f) filtering are Then we have:

[0132]

[0133] Take the inverse Fourier transform of equation (17) and then take the natural exponential, then we have:

[0134]

[0135] in, The absolute values ​​of anisotropic coupling parameters and elastic parameters within the estimated effective frequency band b(f) are usually lower than A i (t), while its stability and fidelity are higher than A i (t). Obviously, when the seismic data itself is noise-free, b(f) = 1, then Equivalent to A i (t).

[0136] Obviously, Equation (18) avoids the error accumulation problem caused by the reflectivity integral in Equation (11), and eliminates the influence of the distorted low-frequency and high-frequency information components. The estimated Stability and fidelity can be effectively guaranteed.

[0137] In addition, considering that when f = 0, there are singular values ​​in Equation (14), although the introduction of regularization terms can ensure the stability of numerical calculations, it is impossible to avoid the r in Equation (15). i The distortion problem after dividing the value near zero frequency by f. Therefore, in the actual calculation of formula (18), the reflectivity r without low-frequency components inverted by formula (10) or (11) can be selected. i ′ is replaced by A calculated with logging data. i (t) Smoothed low-frequency model A i0 (t) is used to compensate for this. Therefore, Equation (18) can be replaced by Equation (19), which can further improve the estimation accuracy of the absolute values ​​of the anisotropic coupling parameters and elastic parameters.

[0138]

[0139] In one example, the expression for calculating the anisotropy parameter of the VTI medium is:

[0140]

[0141] in, is the final estimated anisotropy parameter of the VTI medium, that is, the estimated value of ε, is the absolute value of the first parameter within the estimated effective frequency band b(f), that is, V P0 The estimated value of is the absolute value of the second parameter within the estimated effective frequency band b(f), that is, V P0 e ε Estimated value of .

[0142] Specifically, due to Then we have:

[0143]

[0144] in, This is the final estimated anisotropy parameter of the VTI medium, which can be used to predict the development degree of low-angle fractures in shale layers under certain geological conditions.

[0145] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, but they are not intended to limit the present invention. It should be noted that, unless there is a conflict, the embodiments and features in the embodiments of the present invention may be combined with each other.

[0146] Example 1

[0147] This embodiment provides a pre-stack seismic synchronous inversion method for anisotropic parameters of VTI media, including:

[0148] Obtain the priori model of reflectivity of angle gather seismic data, seismic wavelet and VTI medium anisotropic coupling parameters and elastic parameters; use the well logging data of the marine shale gas target layer in a certain exploration area in the Sichuan Basin to conduct theoretical simulation seismic data testing, and calculate the anisotropic parameters δ and ε, the P-wave velocity V p0 , shear wave velocity V S0 The theoretical curve (real model) of and density ρ is as follows Figure 2 As shown; according to the theoretical curve, the seismic data of the angle gather is simulated, as shown Figure 3 (a) and Figure 3 (b), where Figure 3 (a) is a schematic diagram of noise-free angle gather seismic data. Figure 3 (b) Schematic diagram of seismic data with 10% random noise-free angle gathers, and the prior model of reflectivity of seismic wavelet and VTI medium anisotropic coupling parameters and elastic parameters;

[0149] The theoretical calculation approximate formula of the VTI medium angle domain reflection coefficient and the theoretical calculation approximate formula between the anisotropic coupling parameter reflectivity and the elastic parameter reflectivity are as follows:

[0150] R VTI(θ)=r1+tan 2 θr2-8gsin 2 θr3+(1-4gsin 2 θ)r4, based on the convolution model, a linear equation system is established between the angle gather seismic data and the anisotropic coupling parameter reflectivity and elastic parameter reflectivity: d(t,θ)=w(t)*R VTI (t,θ)=w(t)*(r1(t)+tan 2 θr2(t)-8gsin 2 θr3(t)+(1-4gsin 2 θ)r4(t)), then under the constraints of the prior model of the anisotropic coupling parameter reflectivity and elastic parameter reflectivity of the VTI medium, a target functional is constructed to simultaneously invert the anisotropic coupling parameter and elastic parameter reflectivity of the underground VTI medium: J = ||WLr-d||2+μ||D x (r-r0)|| p ; Among them, R VTI (θ) is the angle domain reflection coefficient, θ is the incident angle of the seismic wave, r1 is the reflectivity of the elastic coefficient A1, A1=V P0 , r2 is the reflectivity of the anisotropic coupling parameter A2, A2=V P0 e ε , r3 anisotropic coupling parameter A3 reflectivity, r4 is the reflectivity of the elastic coefficient A4, A4=ρ,V P0 is the longitudinal wave velocity, V s0 is the shear wave velocity, ρ is the density, Δ is the difference in medium parameters above and below the formation interface (lower layer minus upper layer), ε and δ are the Thomsen anisotropy parameters of the VTI medium, e is a natural constant, g is the square of the ratio of the average shear wave velocity to the average compressional wave velocity, R VTI (t,θ) is the angle domain reflection coefficient of the VTI medium, d(t,θ) is the angle gather seismic data, t is the vertical two-way travel time of the seismic data, w(t) is the seismic wavelet, and * is the convolution operation. Assume that there are M discrete angle gather seismic data, each with N sample points. The angle gather seismic data constitutes a one-dimensional vector d of length N·M. The reflectivity of the four anisotropic coupling parameters and elastic parameters to be inverted (r1(t), r2(t), r3(t), r4(t)) constitutes a one-dimensional vector r of length 4N, r=[r1 r2 r3 r4] T, J is the target functional for simultaneously inverting the anisotropic coupling parameters and elastic parameter reflectivity of VTI media based on angle gather seismic data, W is the N·M×N·M wavelet matrix composed of seismic wavelets w(t), L is the coefficient of the four reflectivities (1, tan 2 θ, -8gsin 2 θ, 1-4gsin 2 θ) constitutes a coefficient matrix of N·M×4N, p is the norm constraint of the reflectivity regularization term, μ is the regularization parameter, μ>0, which is related to the noise level of the angle gather seismic data d. The heavier the noise, the larger μ is to ensure the stability of the solution, D x is the regularization operator, x is the differential order, usually 0, 1 or 2, x 0 represents D x is the unit matrix I, x takes 1 to represent D x is the first-order differential operator D1; x is 2 to represent D x is the prior model of the second-order differential operator D2 r0 is r;

[0151]

[0152] like Figure 4 As shown, the angle gather seismic data, seismic wavelet and prior model are input into the target functional, and the reflectivity of the four anisotropic coupling parameters and elastic parameters of the VTI medium under different regularization conditions is calculated by the optimization algorithm; the different regularization conditions include: p is taken as 2 to represent the reflectivity regularization term as a 2-norm constraint, and a smooth reflectivity r solution is obtained; p is taken as 1 to represent the reflectivity regularization term as a 1-norm constraint, and a sparse edge-preserving reflectivity r solution is obtained; set p to 2; calculate the reflectivity of the four anisotropic coupling parameters and elastic parameters of the VTI medium when p = 2 or p = 1;

[0153] When p = 2, set the values ​​of x and μ according to the actual situation, take the minimum value of the target functional J, and use the damped least squares method to solve the four reflectivity r:

[0154]

[0155] When p=1, set x, μ, η and it according to the actual situation max The value of , takes the minimum value of the target functional J, and uses the iterative reweighted least squares method to solve the four reflectance r:

[0156]

[0157] Among them, k represents the kth iteration, it max The maximum number of iterations set for iterative solution, Q kis the reweighting matrix of the kth iteration, η is a small positive constant. According to the different values ​​of p, the four reflectivities r or r′ (r′ is the four-parameter reflectivity without low-frequency components) of the anisotropic coupling parameters and elastic parameters of the VTI medium under two different regularization conditions can be inverted. r' obtained at the kth iteration (k) The jth element of the vector after the differential operation is summed with η and the reciprocal is taken as Q k The jth element of the diagonal, L T is the transposed matrix of L, W T is the transposed matrix of W, D x T D x The transposed matrix of

[0158] Input filter b(f), regularization factor a, low-frequency model A i0 And reflectivity r, according to the expression:

[0159]

[0160] Estimate the absolute values ​​of the anisotropic coupling parameters and elastic parameters of the VTI medium, where i = 1, 2, 3 or 4, are the absolute values ​​of the anisotropic coupling parameters and elastic parameters within the estimated effective frequency band b(f), b(f) is the filter with an asymmetric structure, FT -1 (·) represents the inverse Fourier transform, i is the imaginary unit, A i0 A i (t), a is the regularization factor, which is a small value greater than zero, e is a natural constant, T(f) is the Tikhonov smoothing operator as a regularization term, T(f) = 1 / |f| n , n>1, r i ′(f) is the reflectivity spectrum of the i-th parameter in the four-parameter reflectivity without low-frequency components.

[0161] According to the absolute values ​​and expressions of the anisotropic coupling parameters and elastic parameters of VTI media:

[0162] Calculate the anisotropic parameters of VTI medium; This is the final estimated anisotropy parameter of the VTI medium, that is, the estimated value of ε. is the absolute value of the first parameter within the estimated effective frequency band b(f), that is, V P0 The estimated value of is the absolute value of the second parameter within the estimated effective frequency band b(f), that is, V P0 e ε estimated value of;

[0163] The comparison results of anisotropy parameters inverted from angle gather seismic data with the real model and prior model are as follows: Figure 5 (a) and Figure 5 (b), where Figure 5 The anisotropy parameters in (b) are based on Figure 3 (b) is obtained by inverting the angle gather seismic data containing 10% random noise; when the theoretically simulated angle gather seismic data is noise-free, the inversion result is highly consistent with the true model, and the average relative error is less than 2.8%, indicating that the method proposed in the present invention can estimate the anisotropy parameter ε with high precision; when the theoretically simulated angle gather seismic data contains 10% random noise, the degree of consistency between the inversion result and the true model is slightly lower than that of the noise-free inversion result, but is still high as a whole, and the average relative error is still less than 4.5%, indicating that the method proposed in the present invention has certain noise resistance and can obtain relatively stable and accurate anisotropy parameter ε from data containing a certain amount of noise.

[0164] Example 2

[0165] like Figure 6 As shown, this embodiment provides a VTI medium anisotropic parameter pre-stack seismic synchronous inversion system, which is used to execute the VTI medium anisotropic parameter pre-stack seismic synchronous inversion method described in Example 1, including:

[0166] Acquisition module, used to obtain the angle gather seismic data, seismic wavelet and reflection prior model of VTI medium anisotropic coupling parameters and elastic parameters;

[0167] A construction module is used to construct a target functional for simultaneously inverting the anisotropic coupling parameters and elastic parameter reflectivity of VTI media based on angle gather seismic data;

[0168] The input and calculation module is used to input the angle gather seismic data, seismic wavelet and prior model into the target functional, and calculate the reflectivity of the four anisotropic coupling parameters and elastic parameters of the VTI medium under different regularization conditions through the optimization algorithm;

[0169] An estimation module, for estimating the absolute values ​​of the anisotropic coupling parameters and elastic parameters of the VTI medium based on the conversion relationship between the reflectivity of the anisotropic coupling parameters and elastic parameters and their absolute values ​​in the frequency domain and the reflectivity;

[0170] The calculation module is used to calculate the anisotropic parameters of the VTI medium according to the coupling relationship and absolute values ​​between the anisotropic coupling parameters and the elastic parameters.

[0171] While various embodiments of the present invention have been described above, the above description is intended to be illustrative, not exhaustive, and not limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments.

Claims

1. A pre-stack seismic synchronous inversion method for anisotropic parameters of VTI media, characterized by: include: Acquiring angle gather seismic data, seismic wavelets, and a priori model, wherein the priori model is a priori model of reflectivity of anisotropic coupling parameters and elastic parameters of a VTI medium; Construct a target functional for the simultaneous inversion of the reflectivity of the anisotropic coupling parameters and elastic parameters of VTI media based on angle gather seismic data; Inputting the angle gather seismic data, the seismic wavelet and the prior model into the target functional, and calculating the reflectivity of four VTI medium anisotropic coupling parameters and elastic parameters under different regularization conditions through an optimization algorithm; estimating the absolute values ​​of the anisotropic coupling parameter and the elastic parameter of the VTI medium based on the conversion relationship between the reflectivity of the anisotropic coupling parameter and the elastic parameter and the absolute value thereof in the frequency domain and the reflectivity; The anisotropic parameters of the VTI medium are calculated according to the coupling relationship between the anisotropic coupling parameters and the elastic parameters and the absolute values.

2. The method for prestack seismic synchronous inversion of VTI medium anisotropy parameters according to claim 1, characterized in that: The target functional for simultaneously inverting the reflectivity of the anisotropic coupling parameters and elastic parameters of the VTI medium based on the angle gather seismic data specifically includes: Based on the theoretical calculation approximate formula of the angle domain reflection coefficient of VTI media and the convolution model, a group of linear equations for the angle gather seismic data and the anisotropic coupling parameter reflectivity and elastic parameter reflectivity are established. Combined with the prior model of the anisotropic coupling parameter reflectivity and elastic parameter reflectivity of the VTI medium, a target functional for simultaneously inverting the reflectivity of the underground anisotropic coupling parameters and elastic parameters based on the angle gather seismic data is constructed.

3. The method for pre-stack seismic synchronous inversion of VTI medium anisotropy parameters according to claim 2, characterized in that: The theoretical calculation approximate formula of the VTI medium angle domain reflection coefficient is expressed as follows: ； in, is the angular domain reflection coefficient of the VTI medium, is the incident angle of the seismic wave, is the elastic parameter The reflectivity, = , = , is the anisotropic coupling parameter The reflectivity, = , , Anisotropic coupling parameters The reflectivity, , , is the elastic parameter The reflectivity, , , is the longitudinal wave velocity, is the shear wave velocity, is the density, is the difference in medium parameters above and below the formation interface, and is the Thomsen anisotropy parameter of the VTI medium, is a natural constant, is the square of the ratio of the average velocity of transverse waves to that of longitudinal waves, .

4. The method for pre-stack seismic synchronous inversion of VTI medium anisotropy parameters according to claim 3, characterized in that: The expression of the linear equation system is: ； in, is the angle domain reflection coefficient, is the angle gather seismic data, is the vertical two-way travel time of seismic data, is the seismic wavelet, is the convolution operation.

5. The method for pre-stack seismic synchronous inversion of VTI medium anisotropic parameters according to claim 4, characterized in that: The target functional is: ； Among them, suppose that the total number of discrete angle gather seismic data is Dao, each Dao has sample points, the length of the angle gather seismic data is One-dimensional vector of , the reflectivity of the four anisotropic coupling parameters and elastic parameters to be inverted 、 、 and The length of the component is 4 N One-dimensional vector of , , It is a target functional for simultaneously inverting the anisotropic coupling parameters and elastic parameter reflectivity of VTI media based on angle gather seismic data. Seismic wavelet composed of wavelet matrix, is the coefficient 1 of the 4-term reflectivity, , , constitute The coefficient matrix of is the norm constraint of the reflectivity regularization term, is the regularization parameter, >0, is the regularization operator, is the differential order, for The prior model.

6. The method for pre-stack seismic synchronous inversion of VTI medium anisotropic parameters according to claim 5, characterized in that: The different regularization conditions include: Taking 2 to represent the reflectivity regularization term is a 2-norm constraint, we get the smooth reflectivity untie; Taking 1 to represent the reflectivity regularization term is 1 norm constraint, we get the sparse edge-preserving reflectivity untie.

7. The method for pre-stack seismic synchronous inversion of VTI medium anisotropic parameters according to claim 6, characterized in that: The calculation of the reflection of the four anisotropic coupling parameters and elastic parameters of the VTI medium under different regularization conditions by the optimization algorithm specifically includes: when When the target functional is taken Minimum, using the damped least squares method to solve the reflectivity of the four VTI medium anisotropic coupling parameters and elastic parameters : ； when When the target functional is taken The minimum value is obtained by using the iterative reweighted least squares method to solve the reflectivity of the four anisotropic coupling parameters and elastic parameters of the VTI medium. : ； in, k For the setting k iterations, it max The maximum number of iterations set for iterative solution, Q k For the k The reweighting matrix of the iteration, η is a small positive constant set, is the reflectivity of the four parameters without low-frequency components, For the k The iterative The first vector after differentiation j elements, and η Sum and take the reciprocal as Q k The diagonal j elements, for The transposed matrix of for The transposed matrix of for The transposed matrix of .

8. The method for pre-stack seismic synchronous inversion of VTI medium anisotropic parameters according to claim 7, characterized in that: The expressions for estimating the absolute values ​​of the anisotropic coupling parameters and elastic parameters of the VTI medium are: ； in, i =1, 2, 3 or 4, is the estimated effective frequency band b ( f ) in the range of the absolute values ​​of the anisotropic coupling parameters and elastic parameters, b ( f ) is an asymmetric filter, FT -1 ( ) represents the inverse Fourier transform, i is the imaginary unit, , for The average value of a is the regularization factor, which is a small value greater than zero. is a natural constant, is the Tikhonov smoothing operator as a regularization term, , n>1, The first of the four parameters of reflectivity without low-frequency components i Reflectance spectrum of the item parameters.

9. The method for pre-stack seismic synchronous inversion of VTI medium anisotropic parameters according to claim 8, characterized in that: The expression for calculating the anisotropic parameter of the VTI medium is: ； in, is the final estimated anisotropic parameter of the VTI medium, namely The estimated value of is the estimated effective frequency band b ( f ) within the range of the first parameter absolute value, that is The estimated value of is the estimated effective frequency band b ( f ) within the range of the second parameter absolute value, that is Estimated value of .

10. A pre-stack seismic synchronous inversion system for VTI media anisotropic parameters, characterized by: An acquisition and input module is used to acquire angle gather seismic data, seismic wavelets and a priori models, wherein the priori models are a priori models of reflectivity of anisotropic coupling parameters and elastic parameters of VTI media; A construction module for constructing a target functional for simultaneously inverting the reflectivity of anisotropic coupling parameters and elastic parameters of VTI media based on angle gather seismic data; An input and calculation module, configured to input the angle gather seismic data, the seismic wavelet, and the prior model into the target functional, and calculate the reflectivity of four VTI medium anisotropic coupling parameters and elastic parameters under different regularization conditions through an optimization algorithm; an estimation module, configured to estimate the absolute values ​​of the anisotropic coupling parameter and the elastic parameter of the VTI medium based on a conversion relationship between the reflectivity of the anisotropic coupling parameter and the elastic parameter and their absolute values ​​in the frequency domain and the reflectivity; The calculation module is used to calculate the anisotropic parameters of the VTI medium according to the coupling relationship between the anisotropic coupling parameters and the elastic parameters and the absolute value.

Images