A high-resolution direct inversion method for elastic impedance without well constraints in the frequency domain

Abstract

The present invention discloses a frequency-domain, well-free, high-resolution direct inversion method for elastic impedance, comprising the following steps: obtaining a relative elastic impedance expression and an elastic impedance low-frequency trend model; obtaining a well-free, frequency-domain absolute elastic impedance direct inversion calculation formula and a seismic wavelet energy correction factor based on the relative elastic impedance expression and the elastic impedance low-frequency trend model; and obtaining a frequency-domain, high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction based on the well-free, frequency-domain absolute elastic impedance direct inversion calculation formula and the seismic wavelet energy correction factor. The present invention constructs a frequency-domain, high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction. The obtained elastic impedance has high accuracy and is suitable for exploration in early-stage exploration areas with few or no wells. The calculation is simple, avoiding the strong dependence of well logging-constrained seismic inversion on well information and the poor stability and low accuracy of traditional direct inversion.

Description

Technical Field

[0001] The present invention belongs to the seismic inversion technology in the field of seismic exploration and relates to a high-resolution elastic impedance direct inversion method in the frequency domain without well constraints. Background Art

[0002] Seismic inversion can extract elastic parameters closely related to formation lithology identification, lithofacies analysis, and reservoir description from seismic data. The commonly used post-stack acoustic impedance inversion method has inherent advantages in inversion accuracy and stability, but its ability to mine oil and gas information is limited due to the lack of utilization of angle-domain AVO information. Although pre-stack AVO inversion can theoretically directly obtain multiple elastic parameters closely related to lithology and even oil and gas content, such as P-wave and S-wave velocities, density, etc., it is affected by factors such as uneven near-, mid-, and far-range coverage of pre-stack seismic data, large quality differences, lack of or distortion of high-angle reflection information, and the accuracy of obtaining the Zoepritz equation coefficient matrix. As a result, the application results of the pre-stack three-parameter simultaneous inversion in the exploration of complex lithologic oil and gas reservoirs often fall far short of the theoretical results.

[0003] To reconcile the computational practicality of stacked acoustic inversion with the information mining advantages of prestack AVO inversion, Connolly (1999) proposed the definition of elastic impedance (EI). This definition is a function of P-wave velocity, S-wave velocity, density, and angle of incidence. Its value varies with the angle of incidence θ, reflecting the relationship between the amplitude of the reflected wave at nonzero offset (nonzero angle of incidence θ) and the change in elastic impedance of the strata on both sides of the interface. Its mathematical form is consistent with acoustic impedance, which is a special case of elastic impedance at normal incidence. Because elastic impedance inversion can be flexibly performed using stacked data in the angle domain, it can naturally extract seismic AVO information, often obtaining elastic impedance information that highlights lithologic or fluid anomalies (Wang et al., 2005). This definition allows elastic impedance inversion to be performed based on the principles and methods of post-stack seismic inversion. By optimizing the angular domain stacking data at different angles to invert the elastic impedance, the required elastic parameters (such as P- and S-wave velocities, P- and S-wave impedances, density, Lame constant, shear modulus, and Poisson's ratio) are then calculated from the elastic impedance at different angles (Cao et al., 2006). Therefore, the accuracy of elastic parameter determination depends on the precision (stability and resolution) of the elastic impedance inversion results.

[0004] Elastic impedance inversion is mainly divided into two categories: direct inversion and model-based iterative inversion. Direct inversion is mainly divided into trace integration and recursive inversion. The former is limited by the effective seismic bandwidth and has low resolution. It cannot adapt to the needs of thin-layer interpretation and cannot obtain the absolute wave impedance of the formation, making it difficult to directly use it for lithology identification. The latter is affected by the large cumulative error when recursively recursing the reflection coefficient, and the reliability of the inversion results is difficult to guarantee. Model-based iterative inversion requires the use of a large amount of well logging and geological data to establish a relatively accurate initial model (including some high-frequency and low-frequency information missing from the seismic data itself). In areas with few wells, the application effect is difficult to guarantee. When the inversion parameters and initial model settings are not ideal, it is easy to fall into local optimality, and the inversion results cannot reflect the actual underground geological conditions.

[0005] Ivan et al. (2012) proposed a one-step direct inversion method for seismic acoustic impedance. This seismic inversion framework overcomes the drawbacks of traditional direct inversion and theoretically proposes a method for directly deriving high-resolution acoustic impedance from seismic data. However, due to the lack of consideration of the impact of low- and high-frequency signals distorted by actual seismic data on the inversion results (i.e., the lack of consideration of the optimal utilization of effective high- and low-frequency information in post-stack seismic data to obtain robust high-resolution acoustic impedance), the inversion results contain low- and high-frequency anomalies, significantly reducing the robustness and interpretability of the inversion results. Furthermore, the problem of deriving high-resolution absolute acoustic impedance when wavelet energy correction is not possible through well-seismic calibration is not considered. Therefore, directly applying Ivan's inversion framework to elastic impedance inversion also presents numerous problems.

[0006] It is well known that in early exploration areas with few or no wells, it is difficult to perform well-seismic calibration to correct the wavelet energy estimated by statistical methods and establish a high-precision initial model. Therefore, how to effectively perform high-resolution direct inversion of elastic impedance without well constraints to obtain high-precision elastic impedance is of great practical significance.

[0007] Therefore, a high-resolution direct inversion method of elastic impedance without well constraints is particularly needed to obtain elastic impedance with higher accuracy. Summary of the Invention

[0008] The purpose of the present invention is to propose a high-resolution direct inversion method of elastic impedance without well constraints, which can obtain elastic impedance with higher accuracy.

[0009] In view of this, the present invention provides a frequency domain high-resolution direct elastic impedance inversion method without well constraints, which at least solves the problem in the prior art that direct elastic impedance inversion for work areas with few or no wells cannot obtain high-precision elastic impedance.

[0010] The present invention provides a frequency domain high-resolution direct inversion method for elastic impedance without well constraints, comprising the following steps: obtaining a relative elastic impedance expression and an elastic impedance low-frequency trend model; obtaining a frequency domain absolute elastic impedance direct inversion calculation formula without well constraints and a seismic wavelet energy correction factor based on the relative elastic impedance expression and the elastic impedance low-frequency trend model; and obtaining a frequency domain high-resolution direct inversion calculation formula for elastic impedance based on low-frequency trend model energy correction based on the frequency domain absolute elastic impedance direct inversion calculation formula without well constraints and the seismic wavelet energy correction factor.

[0011] Optionally, obtaining the relative elastic impedance expression includes: obtaining a basic relationship between the frequency domain seismic trace and the relative elastic impedance; introducing a regularization term into the basic relationship between the frequency domain seismic trace and the relative elastic impedance; adding a filter to the basic relationship between the frequency domain seismic trace and the relative elastic impedance after the regularization term is introduced to obtain a filtered frequency domain relative elastic impedance expression; performing an inverse Fourier transform on the filtered frequency domain relative impedance expression to obtain a relative elastic impedance expression.

[0012] Optionally, the basic relationship between the frequency domain seismic trace and the relative elastic impedance is:

[0013]

[0014] Wherein, EI(t,θ) is the absolute elastic impedance, InEI(t,θ) is the natural logarithm of EI(t,θ), which is the relative elastic impedance, t is the propagation time of the seismic reflection wave with an incident angle of θ in the formation, t∈[T1,T2], [T1,T2] is the time period of the formation elastic impedance to be inverted, θ is the incident angle, S(f,θ) is the frequency spectrum of the angle domain stacked seismic trace s(t,θ) with an incident angle of θ, R(f,θ) is the frequency spectrum of the angle domain reflection coefficient r(t,θ), W(f,θ) is the frequency spectrum of the seismic wavelet w(t,θ) of the angle domain stacked seismic trace s(t,θ), f is the frequency, i is the imaginary unit, FT(·) is the Fourier transform.

[0015] Optionally, the basic relationship between the frequency domain seismic trace and the relative elastic impedance after the regularization term is introduced is:

[0016]

[0017] Where W(f,θ) * is the complex conjugate of W(f,θ), I(f) is the regularization term, and μ is the regularization factor.

[0018] Optionally, the frequency domain relative elastic impedance expression after filtering is:

[0019]

[0020] Among them, FT b [In EI(t,θ)] is the relative elastic impedance in the frequency domain after filtering, and the filter b(f) is:

[0021]

[0022] f l and f h is the frequency parameter, λ1,λ2∈(0,1) is the ratio of the filter length to the cosine rim length, and B is the frequency band.

[0023] Optionally, the relative elastic impedance expression is:

[0024]

[0025] Among them, FT- 1 is the inverse Fourier transform.

[0026] Optionally, the elastic impedance low-frequency trend model is:

[0027]

[0028] Among them, V p0 is the low-frequency trend model of longitudinal wave velocity V p The average value of (t), is the low-frequency trend model of shear wave velocity V s (t), ρ0 is the average value of the density low-frequency trend model ρ(t), K = (Vs1 + Vs2) / (Vp1 + Vp2) is the proportional factor of the P-wave and S-wave velocities of the adjacent layers above and below the formation time interface, Vs1 is the S-wave velocity of the upper layer, Vs2 is the S-wave velocity of the lower layer, Vp1 is the P-wave velocity of the upper layer, Vp2 is the P-wave velocity of the lower layer, and θ is the incident angle.

[0029] Optionally, the frequency domain absolute elastic impedance direct inversion calculation formula without well constraints is:

[0030]

[0031] where EI'(t,θ) is the absolute elastic impedance without correction of the energy of the seismic wavelet w(t,θ).

[0032] Optionally, the seismic wavelet w(t, θ) energy correction factor is obtained by: based on the relative elastic impedance expression, obtaining the maximum and minimum values of the relative elastic impedance that is not corrected for the seismic wavelet energy within the time period of the formation elastic impedance to be inverted; based on the elastic impedance low-frequency trend model, obtaining the relative elastic impedance of the maximum and minimum values of the low-frequency trend model of the actual underground elastic impedance within the time period of the formation elastic impedance to be inverted; and calculating the seismic wavelet energy correction factor using the following formula:

[0033]

[0034] Where a is the wavelet energy correction factor, In EI′ max and In EI′ min In EI is the maximum and minimum relative elastic impedance without correction of seismic wavelet energy. 0max and In EI 0min is the maximum value EI of the low-frequency trend model of the true elastic impedance of the underground 0max and minimum EI 0min The relative elastic impedance.

[0035] Optionally, the frequency-domain high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction is:

[0036]

[0037] Where EI(t,θ)” is the high-resolution absolute elastic impedance after correction of the seismic wavelet energy based on the low-frequency trend model and is the final result of the high-resolution direct inversion of elastic impedance without well constraints in the frequency domain.

[0038] The beneficial effects of the present invention are as follows: the frequency domain high-resolution elastic impedance direct inversion method without well constraints of the present invention constructs a frequency domain high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction, the obtained elastic impedance has high accuracy, does not require high abundance of prior information, does not require fine calibration of well seismic data to establish a high-precision initial model and extract energy-corrected wavelets, is suitable for exploration areas with no or few wells in the early stages of exploration, has simple calculations, avoids the strong dependence of logging-constrained seismic inversion on well information and the problems of poor stability and low accuracy of traditional direct inversion, and achieves a good balance between the stability and resolution of the results.

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

[0040] 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, wherein like reference numerals generally represent like components throughout the exemplary embodiments of the present invention.

[0041] Figure 1 The flowchart of a frequency domain high-resolution direct inversion method of elastic impedance without well constraints according to one embodiment of the present invention is shown.

[0042] Figure 2 A detailed flow chart of a frequency-domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown.

[0043] Figure 3 A Tukey bandpass filter diagram of an asymmetric structure of a frequency domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown.

[0044] Figure 4a An angle domain stacked imaging seismic trace diagram of a frequency domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown.

[0045] Figure 4b The present invention shows a low-frequency trend model of elastic impedance of a frequency-domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention.

[0046] Figure 4c FIG2 shows an absolute elastic impedance inversion result diagram of a frequency domain high-resolution elastic impedance direct inversion method without well constraints according to an embodiment of the present invention.

[0047] Figure 5a A seismic wavelet diagram before low-frequency trend model energy correction is shown for a frequency-domain high-resolution elastic impedance direct inversion method without well constraints according to an embodiment of the present invention.

[0048] Figure 5b A seismic wavelet diagram after energy correction based on a low-frequency trend model is shown for a frequency-domain high-resolution elastic impedance direct inversion method without well constraints according to an embodiment of the present invention. DETAILED DESCRIPTION

[0049] The preferred embodiments of the present invention will be described in more detail below. Although the preferred embodiments of the present invention are described below, it should be understood that the present invention can be implemented in various forms and should not be limited to the embodiments set forth herein.

[0050] The present invention provides a frequency domain high-resolution elastic impedance direct inversion method without well constraints, comprising the following steps: obtaining a relative elastic impedance expression and an elastic impedance low-frequency trend model; obtaining a frequency domain absolute elastic impedance direct inversion calculation formula without well constraints and a seismic wavelet energy correction factor based on the relative elastic impedance expression and the elastic impedance low-frequency trend model; and obtaining a frequency domain high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction based on the frequency domain absolute elastic impedance direct inversion calculation formula without well constraints and the seismic wavelet energy correction factor.

[0051] Specifically, a relative elastic impedance expression is obtained, and the expression is used to count the maximum and minimum values of the relative elastic impedance without correction of the seismic wavelet energy. The elastic impedance low-frequency trend model is input to count the relative elastic impedance of the minimum and maximum values of the low-frequency component of the true underground elastic impedance. The maximum and minimum values of the relative elastic impedance without correction of the seismic wavelet energy and the relative elastic impedance of the minimum and maximum values of the low-frequency component of the true underground elastic impedance are substituted into the seismic wavelet energy correction factor formula to calculate the energy correction factor. The obtained energy correction factor is then added to the frequency domain absolute elastic impedance direct inversion calculation formula without well constraints to obtain the frequency domain high-resolution elastic impedance direct inversion calculation formula based on the energy correction of the low-frequency trend model, and then the final absolute elastic impedance is calculated.

[0052] According to an exemplary embodiment, a frequency-domain high-resolution elastic impedance direct inversion method without well constraints constructs a frequency-domain high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction. The obtained elastic impedance has high accuracy and does not require high prior information abundance. It does not require fine well-seismic calibration to establish a high-precision initial model and extract the energy-corrected wavelet. It overcomes the problem that elastic impedance inversion methods based on high-precision initial model iteration or strong logging constraints are not suitable for insufficient prior information in the early stages of exploration. It is suitable for exploration areas with no or few wells in the early stages of exploration. The calculation is simple, avoiding the strong dependence of logging-constrained seismic inversion on well information and the poor stability and low accuracy of traditional direct inversion. The stability and resolution of the results are well balanced.

[0053] As an optional solution, obtaining the relative elastic impedance expression includes: obtaining a basic relationship between frequency domain seismic traces and relative elastic impedance; introducing a regularization term into the basic relationship between frequency domain seismic traces and relative elastic impedance; adding a filter into the basic relationship between frequency domain seismic traces and relative elastic impedance after the regularization term is introduced to obtain a filtered frequency domain relative elastic impedance expression; performing an inverse Fourier transform on the filtered frequency domain relative elastic impedance expression to obtain the relative elastic impedance expression.

[0054] Specifically, first, the input angle domain stacking trace and the seismic wavelet estimated by the statistical method are converted into frequency domain data using discrete Fourier transform to obtain the basic relationship between the frequency domain seismic trace and the relative elastic impedance. Then, a regularization term is introduced and a filter is added to obtain the filtered frequency domain relative elastic impedance expression. The filtered frequency domain relative elastic impedance expression is subjected to inverse Fourier transform to obtain the relative elastic impedance.

[0055] As an optional solution, the basic relationship between frequency domain seismic traces and relative elastic impedance is:

[0056]

[0057] Wherein, EI(t,θ) is the absolute elastic impedance, InEI(t,θ) is the natural logarithm of EI(t,θ), which is the relative elastic impedance, t is the propagation time of the seismic reflection wave with an incident angle of θ in the formation, t∈[T1,T2], [T1,T2] is the time period of the formation elastic impedance to be inverted, θ is the incident angle, S(f,θ) is the frequency spectrum of the angle domain stacked seismic trace s(t,θ) with an incident angle of θ, R(f,θ) is the frequency spectrum of the angle domain reflection coefficient r(t,θ), W(f,θ) is the frequency spectrum of the seismic wavelet w(t,θ) of the angle domain stacked seismic trace s(t,θ), f is the frequency, i is the imaginary unit, FT(·) is the Fourier transform.

[0058] Specifically, assuming that the absolute elastic impedance EI(t,θ) of the formation is continuously differentiable with the seismic wave propagation time t, the formation angle domain reflection coefficient r(t,θ) can be approximately equivalent to:

[0059] r(t,θ)≈d[In EI(t,θ)] / 2dt (1)

[0060] Where In EI(t,θ) is usually defined as the relative elastic impedance of the formation, and θ is the incident angle of the seismic wave. Obviously, the frequency domain expression of Equation (1) is R(f,θ)≈if / 2FT[In EI(t,θ)]. Combined with the frequency domain expression of the seismic convolution model S(f,θ)=R(f,θ)W(f,θ), the basic relationship between the frequency domain seismic trace and the relative elastic impedance can be derived as follows:

[0061]

[0062] Where S(f, θ) is the frequency spectrum of the angle-domain stacked seismic trace s(t, θ) with an incident angle of θ, R(f, θ) is the frequency spectrum of the formation angle-domain reflection coefficient r(t, θ), W(f, θ) is the frequency spectrum of the seismic wavelet w(t, θ) of the angle-domain stacked seismic data s(t, θ), f is the frequency, i is the imaginary unit, FT(·) is the Fourier transform. Obviously, the process of spectrum division by wavelet in Equation (2) has a frequency-increasing effect, which improves the resolution of elastic impedance.

[0063] As an optional solution, the basic relationship between the frequency domain seismic trace and the relative elastic impedance after the regularization term is introduced is:

[0064]

[0065] Where W(f,θ) * is the complex conjugate of W(f,θ), I(f) is the regularization term, and μ is the regularization factor.

[0066] Specifically, considering that the spectrum of actual seismic data is of limited bandwidth, below the lowest cutoff low frequency f l (SNR spectrum is less than 1) and above the highest cutoff frequency f h The spectral components of the signal-to-noise ratio spectrum (SNR spectrum is less than 1) are usually distorted. The seismic wavelet estimated by the average statistical results of the seismic data amplitude spectrum also has similar properties. Therefore, in order to stably solve Equation (2), a regularization term must be introduced. Then Equation (2) can be rewritten as:

[0067]

[0068] Where W(f,θ) * is the complex conjugate of W(f, θ), I(f)>0 is the regularization term, and μ>0 is the regularization factor. Considering the frequency term f in the denominator of Equation (3), it is obvious that Equation (3) is stable for high-frequency components, but unstable for low-frequency components. Therefore, the design of the regularization term should be targeted at low-frequency conditions. Since the previous article assumes that the formation wave impedance Z(t) is continuously differentiable with time and depth, in order to ensure the smoothness of the result, the Tikhonov smoothing operator I(f)=1 / |f| can be introduced. p As a regularization term, in order to ensure that Equation (3) is stable in the extreme case f→0, p should be greater than 1, and here p can be 2.

[0069] As an alternative, the frequency-domain relative elastic impedance after filtering can be expressed as:

[0070]

[0071] Where the filter b(f) is:

[0072]

[0073] f l and f h is the frequency parameter, λ1,λ2∈(0,1) is the ratio of the cosine rim length of the bandpass filter to the filter length, and B is the frequency band range.

[0074] Specifically, although the high-resolution spectrum of the relative elastic impedance is obtained by equation (3), it is lower than the lowest cutoff low frequency f l and above the highest cutoff frequency f h The spectrum components of the frequency spectrum are also unreliable. The frequency increase process will amplify the influence of high and low frequency noise. Through the Fourier inverse operation FT- 1 The relative elastic impedance obtained by (·) will be contaminated. Therefore, a filtering operation needs to be added before the inverse operation to optimize the effective low-frequency and high-frequency information to increase the frequency while suppressing the influence of the distortion spectrum. The Tukey window function can better achieve this goal. Therefore, an asymmetric filter b(f) is designed as follows:

[0075]

[0076] Among them, the frequency band of the filter is B=f h -f l In actual work, f l and f h The value of can usually be determined by frequency scanning or signal-to-noise ratio spectrum estimation. λ1,λ2∈(0,1) is the ratio of the cosine rim length of the bandpass filter to the filter length. Generally, since mid- and low-frequency information plays a controlling role in the overall inversion grid, in order to reduce the loss of mid- and low-frequency energy, λ1 is set to a small value and less than λ2. Therefore, the frequency domain relative elastic impedance after filtering is:

[0077]

[0078] As an alternative, the relative elastic impedance expression is:

[0079]

[0080] Among them, FT- 1 is the inverse Fourier transform.

[0081] Specifically, the filtered frequency domain relative elastic impedance is subjected to inverse Fourier transform to obtain the time domain relative elastic impedance InEI(t,θ):

[0082]

[0083] As an alternative, the elastic impedance low-frequency trend model is:

[0084]

[0085] Among them, V p0 is the low-frequency trend model of longitudinal wave velocity V p The average value of (t), is the low-frequency trend model of shear wave velocity Vs (t), ρ0 is the average value of the density low-frequency trend model ρ(t), K = (Vs1 + Vs2) / (Vp1 + Vp2) is the proportional factor of the P-wave and S-wave velocities of the adjacent layers above and below the formation time interface, Vs1 is the S-wave velocity of the upper layer, Vs2 is the S-wave velocity of the lower layer, Vp1 is the P-wave velocity of the upper layer, Vp2 is the P-wave velocity of the lower layer, and θ is the incident angle.

[0086] As an alternative, the direct inversion calculation formula for absolute elastic impedance in the frequency domain without well constraints is:

[0087]

[0088] where EI'(t,θ) is the absolute elastic impedance without correction of the energy of the seismic wavelet w(t,θ).

[0089] Specifically, according to the absolute elastic impedance calculation formula: EI(t,θ)=EI0(t,θ)e InEI(t,θ) , the direct inversion formula for absolute elastic impedance in the frequency domain without well constraints can be derived:

[0090]

[0091] in, is the elastic impedance low-frequency trend model, and K = (Vs1 + Vs2) / (Vp1 + Vp2) is the proportional factor for the P-wave and S-wave velocities of the adjacent layers above and below the formation time interface t. Clearly, the elastic impedance low-frequency trend model for this exploration area can be easily constructed by processing seismic data from neighboring wells or the current exploration area under the same sedimentary background to obtain an interval velocity model and empirical rock physics formulas.

[0092] As an optional solution, the following steps are used to obtain the seismic wavelet energy correction factor: based on the relative elastic impedance expression, the maximum and minimum values of the relative elastic impedance that are not corrected for the seismic wavelet energy are obtained within the time period of the formation elastic impedance to be inverted; based on the elastic impedance low-frequency trend model, the relative elastic impedance of the maximum and minimum values of the low-frequency trend model of the actual underground elastic impedance is obtained within the time period of the formation elastic impedance to be inverted; and the seismic wavelet energy correction factor is calculated using the following formula:

[0093]

[0094] Where a is the seismic wavelet energy correction factor, In EI′ max and In EI′ min In EI is the maximum and minimum relative elastic impedance without correction of seismic wavelet energy. 0max and In EI 0min The maximum value EI of the low-frequency trend model of the true elastic impedance of the underground 0max and minimum EI0min The relative elastic impedance.

[0095] Specifically, when there is no well or the well cannot be effectively calibrated, the seismic wavelet w(t, θ) energy calculated by using the average amplitude spectrum of the seismic data cannot be calibrated. Assuming that w(t, θ) is a times the seismic wavelet energy after fine calibration of the real wavelet or the well logging synthetic record, the seismic reflection coefficient in the frequency domain or time domain is obtained. or That is the true reflection coefficient a of the underground -1 times. At this time, the relative elastic impedance (denoted as In EI′(t,θ)) obtained by directly applying Equation (6) is not energy equivalent to the relative elastic impedance (denoted as In EI(t,θ)) obtained by the true underground reflection coefficient. The absolute elastic impedance value (i.e., EI′(t,θ)) obtained by the relative elastic impedance In EI′(t,θ) and the true elastic impedance low-frequency trend model EI0(t,θ) through Equation (7) will deviate greatly from the true value (denoted as EI(t,θ)) and lose its physical interpretability.

[0096] Therefore, it is necessary to numerically scale the relative elastic impedance energy. Existing technology can easily obtain the low-frequency velocity trend through seismic data processing or obtain the absolute elastic impedance low-frequency background trend model EI0(t,θ) at any angle (generally less than 10Hz, which is approximately the low-frequency component of the real elastic impedance) from the virtual wells in the adjacent area. Energy correction is performed on the relative elastic impedance within the value range (within the time period of the formation elastic impedance to be inverted). The correction method is derived as follows:

[0097] Let the relative elastic impedance of the real underground reflection coefficient In EI(t,θ), In EI′(t,θ) be the relative elastic impedance to be corrected, and a be the wavelet energy correction factor. 0min ,EI 0max ] is the low-frequency component of the underground real elastic impedance, and its value range is usually basically consistent with the range of change of the underground real impedance EI(t,θ), that is, EI min (t,θ)≈EI 0min ,EI max (t,θ)≈EI 0max , then according to the inversion theory of integrated elastic impedance We know that there always exists a t1, t2 such that:

[0098]

[0099] Then, defined by the relative elastic impedance, available:

[0100]

[0101] From formula (9), we can get:

[0102]

[0103] Similarly, the reflection coefficient obtained by using the statistical wavelet before energy correction and the relative elastic impedance In EI′(t,θ) calculated by channel integration also satisfy:

[0104]

[0105] From formula (11), we can get:

[0106]

[0107] Where t1 and t2 are the time positions of the minimum and maximum elastic impedance, respectively, and r(t,θ) is the reflection coefficient. The relative elastic impedance to be corrected for energy is obtained by equation (6) as In EI′(t,θ), and its minimum and maximum values are statistically calculated as In EI′ min and In EI′ max , a can be approximately calculated from formula (10) and formula (12), that is:

[0108]

[0109] It should be noted here that when inverting the entire two-dimensional or three-dimensional angle domain seismic stacking data, the average values of the minimum and maximum values of InEI0(t,θ) and InEI′(t,θ) in the two-dimensional or three-dimensional space should be estimated respectively, and then substituted into formula (13) to calculate a, ensuring that the entire seismic data body uses the same energy correction value, retaining the spatial relative change relationship of the seismic data anomaly, facilitating comparative interpretation and analysis, and preventing the low-frequency trend model from affecting the inversion results due to local distortion. If it is known that the deviations of the minimum and maximum values of each trace of the low-frequency trend elastic impedance model from the minimum and maximum values of the actual underground elastic impedance are small, the spatial variation a can be calculated to carry out correction. In addition, when the established low-frequency trend model is not sufficient to estimate EI more accurately, 0min and EI 0max The EI of the stratum to be inverted can be estimated based on the wells in the work area or adjacent areas, the test data of outcrop rock samples in the work area, the processing velocity model and other results. 0min and EI 0max , as a reference value for relative impedance energy correction, and then repeatedly test and correct EI based on experience 0min and EI 0max , until the absolute elastic impedance inversion result falls within the empirical value range of the exploration area.

[0110] As an alternative, the frequency-domain high-resolution direct inversion calculation formula of elastic impedance based on energy correction of the low-frequency trend model is:

[0111]

[0112] Where EI(t,θ)” is the high-resolution absolute elastic impedance after correction of the seismic wavelet energy based on the low-frequency trend model and is the final result of the high-resolution direct inversion of elastic impedance without well constraints in the frequency domain.

[0113] Specifically, the frequency-domain high-resolution direct inversion calculation formula of elastic impedance based on energy correction of the low-frequency trend model incorporates a robust frequency-boosting term for optimizing the effective spectrum and a wavelet energy correction term based on the low-frequency trend model. This can obtain absolute elastic impedance with a high signal-to-noise ratio and resolution, making it possible to use pre-stack information to identify favorable targets in the early stages of exploration and to extract lithologic and fluid-sensitive elastic parameters based on elastic impedance. It is a concise and efficient method for obtaining elastic impedance with high precision and is highly practical.

[0114] Example 1

[0115] Figure 1 A flow chart of a frequency domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown. Figure 2 A detailed flow chart of a frequency-domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown.

[0116] Figure 3 A Tukey bandpass filter diagram of an asymmetric structure of a frequency domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown. Figure 4a An angle domain stacked imaging seismic trace diagram of a frequency domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention is shown. Figure 4b The present invention shows a low-frequency trend model of elastic impedance of a frequency-domain high-resolution direct inversion method of elastic impedance without well constraints according to an embodiment of the present invention. Figure 4c FIG2 shows an absolute elastic impedance inversion result diagram of a frequency domain high-resolution elastic impedance direct inversion method without well constraints according to an embodiment of the present invention. Figure 5a A seismic wavelet diagram before low-frequency trend model energy correction is shown for a frequency-domain high-resolution elastic impedance direct inversion method without well constraints according to an embodiment of the present invention. Figure 5b A seismic wavelet diagram after energy correction based on a low-frequency trend model is shown for a frequency-domain high-resolution elastic impedance direct inversion method without well constraints according to an embodiment of the present invention.

[0117] Combine Figure 1 、 Figure 2 、 Figure 3 、 Figure 4a 、 Figure 4b 、 Figure 4c 、 Figure 5a and Figure 5b As shown in Figure 2, the frequency domain high-resolution direct inversion method for elastic impedance without well constraints includes:

[0118] Step 1: Obtain the relative elastic impedance expression and elastic impedance low-frequency trend model;

[0119] Among them, obtaining the relative elastic impedance expression includes: obtaining the basic relationship between the frequency domain seismic trace and the relative elastic impedance; introducing a regularization term into the basic relationship between the frequency domain seismic trace and the relative elastic impedance; adding a filter to the basic relationship between the frequency domain seismic trace and the relative elastic impedance after the regularization term is introduced to obtain the filtered frequency domain relative elastic impedance expression; performing an inverse Fourier transform on the filtered frequency domain relative impedance expression to obtain the relative elastic impedance expression.

[0120] Among them, the basic relationship between frequency domain seismic trace and relative elastic impedance is:

[0121]

[0122] Wherein, EI(t,θ) is the absolute elastic impedance, InEI(t,θ) is the natural logarithm of EI(t,θ), which is the relative elastic impedance, t is the propagation time of the seismic reflection wave with an incident angle of θ in the formation, t∈[T1,T2], [T1,T2] is the time period of the formation elastic impedance to be inverted, θ is the incident angle, S(f,θ) is the frequency spectrum of the angle domain stacked seismic trace s(t,θ) with an incident angle of θ, R(f,θ) is the frequency spectrum of the angle domain reflection coefficient r(t,θ), W(f,θ) is the frequency spectrum of the seismic wavelet w(t,θ) of the angle domain stacked seismic trace s(t,θ), f is the frequency, i is the imaginary unit, FT(·) is the Fourier transform.

[0123] Among them, the basic relationship between the frequency domain seismic trace and the relative elastic impedance after the regularization term is introduced is:

[0124]

[0125] Where W(f,θ) * is the complex conjugate of W(f,θ), I(f) is the regularization term, and μ is the regularization factor.

[0126] The frequency domain relative elastic impedance expression after filtering is:

[0127]

[0128] Where the filter b(f) is:

[0129]

[0130] f l and f h is the frequency parameter, λ1,λ2∈(0,1) is the ratio of the cosine rim length of the bandpass filter to the filter length, and B is the frequency band range.

[0131] The relative elastic impedance expression is:

[0132]

[0133] Among them, FT -1 is the inverse Fourier transform.

[0134] Among them, the low-frequency trend model of elastic impedance is:

[0135]

[0136] Among them, V p0 is the low-frequency trend model of longitudinal wave velocity V p The average value of (t), is the low-frequency trend model of shear wave velocity V s (t), ρ0 is the average value of the density low-frequency trend model ρ(t), K = (Vs1 + Vs2) / (Vp1 + Vp2) is the proportional factor of the P-wave and S-wave velocities of the adjacent layers above and below the formation time interface, Vs1 is the S-wave velocity of the upper layer, Vs2 is the S-wave velocity of the lower layer, Vp1 is the P-wave velocity of the upper layer, Vp2 is the P-wave velocity of the lower layer, and θ is the incident angle.

[0137] Step 2: Based on the relative elastic impedance expression and the elastic impedance low-frequency trend model, the frequency-domain absolute elastic impedance direct inversion formula and the seismic wavelet energy correction factor without well constraints are obtained;

[0138] The direct inversion calculation formula for absolute elastic impedance in the frequency domain without well constraints is:

[0139]

[0140] where EI'(t,θ) is the absolute elastic impedance without correction of the energy of the seismic wavelet w(t,θ).

[0141] The seismic wavelet energy correction factor is obtained by the following steps: based on the relative elastic impedance expression, the maximum and minimum values of the relative elastic impedance that are not corrected for the seismic wavelet energy are obtained within the time period of the formation elastic impedance to be inverted; based on the elastic impedance low-frequency trend model, the relative elastic impedance of the maximum and minimum values of the low-frequency trend model of the actual underground elastic impedance is obtained within the time period of the formation elastic impedance to be inverted; and the seismic wavelet energy correction factor is calculated using the following formula:

[0142]

[0143] Where a is the seismic wavelet energy correction factor, In EI′ max and In EI′ min In EI is the maximum and minimum relative elastic impedance without correction of seismic wavelet energy. 0max and In EI 0min is the relative elastic impedance of the maximum and minimum values of the low-frequency trend model of the true elastic impedance of the subsurface.

[0144] Step 3: Based on the frequency domain absolute elastic impedance direct inversion formula without well constraints and the seismic wavelet energy correction factor, the frequency domain high-resolution elastic impedance direct inversion formula based on the low-frequency trend model energy correction is obtained.

[0145] The calculation formula for direct inversion of high-resolution elastic impedance in the frequency domain based on energy correction of the low-frequency trend model is:

[0146]

[0147] where EI(t,θ)” is the high-resolution absolute elastic impedance after correction of the seismic wavelet energy based on the low-frequency trend model.

[0148] The experiment was conducted using actual seismic data from a marine shale gas target layer in an exploration area in the Sichuan Basin. Figure 4a For angle domain stacking imaging seismic traces, Figure 4b The low-frequency trend elastic impedance model reflects that the elastic impedance of the shale layer (between T1 and T3) has a low value trend and the limestone layer (underlying T2) has a high value trend. Figure 4c The absolute elastic impedance is obtained by the method of the present invention. Between T2 and T3, there is a set of stable, continuous but thin low-value elastic impedance strata, representing high-quality shale intervals. Between T1 and T2 is a set of thick mud shale layers containing multiple sets of relatively high-resistance interlayers, reflecting a geological depositional environment different from that of the high-quality shale between T2 and T3. Figure 5b The energy of the wavelet after energy correction is increased to 10 compared with the wavelet before correction. 4The order of magnitude ensures that the elastic impedance inversion results fall close to the numerical range of the low-frequency trend model and have physically interpretable meanings.

[0149] 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 high-resolution direct inversion method for elastic impedance in the frequency domain without well constraints, characterized by: include: Obtain the relative elastic impedance expression and elastic impedance low-frequency trend model; Based on the relative elastic impedance expression and the elastic impedance low-frequency trend model, a frequency domain absolute elastic impedance direct inversion calculation formula without well constraints and a seismic wavelet energy correction factor are obtained; Based on the frequency domain absolute elastic impedance direct inversion calculation formula without well constraints and the seismic wavelet energy correction factor, a frequency domain high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction is obtained.

2. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 1 is characterized in that: The obtaining of the relative elastic impedance expression comprises: Obtain the basic relationship between frequency domain seismic traces and relative elastic impedance; Introducing a regularization term into the basic relationship between the frequency domain seismic trace and the relative elastic impedance; A filter is added to the basic relationship between frequency domain seismic traces and relative elastic impedance after the regularization term is introduced to obtain the filtered frequency domain relative elastic impedance expression; Performing inverse Fourier transform on the filtered frequency-domain relative elastic impedance expression to obtain the relative elastic impedance expression.

3. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 2, characterized in that: The basic relationship between the frequency domain seismic trace and the relative elastic impedance is: Wherein, EI(t,θ) is the absolute elastic impedance, InEI(t,θ) is the natural logarithm of EI(t,θ), which is the relative elastic impedance, t is the propagation time of the seismic reflection wave with an incident angle of θ in the formation, t∈[T1,T2], [T1,T2] is the time period of the formation elastic impedance to be inverted, θ is the incident angle, S(f,θ) is the frequency spectrum of the angle domain stacked seismic trace s(t,θ) with an incident angle of θ, R(f,θ) is the frequency spectrum of the angle domain reflection coefficient r(t,θ), W(f,θ) is the frequency spectrum of the seismic wavelet w(t,θ) of the angle domain stacked seismic trace s(t,θ), f is the frequency, i is the imaginary unit, FT(·) is the Fourier transform.

4. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 3 is characterized in that: The basic relationship between the frequency domain seismic trace and the relative elastic impedance after the regularization term is introduced is: Where W(f,θ) * is the complex conjugate of W(f,θ), I(f) is the regularization term, and μ is the regularization factor.

5. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 4 is characterized in that: The frequency domain relative elastic impedance expression after filtering is: Among them, FT b [InEI(t,θ)] is the relative elastic impedance in the frequency domain after filtering, and the filter b(f) is: f l and f h is the frequency parameter, λ1,λ2∈(0,1) is the ratio of the filter length to the cosine rim length, and B is the frequency band.

6. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 5, characterized in that: The relative elastic impedance expression is: Among them, FT -1 is the inverse Fourier transform.

7. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 1, characterized in that: The elastic impedance low-frequency trend model is: Among them, V p0 is the low-frequency trend model of longitudinal wave velocity V p The average value of (t), is the low-frequency trend model of shear wave velocity V s (t), ρ0 is the average value of the density low-frequency trend model ρ(t), K = (Vs1 + Vs2) / (Vp1 + Vp2) is the proportional factor of the P-wave and S-wave velocities of the adjacent layers above and below the formation time interface, Vs1 is the S-wave velocity of the upper layer, Vs2 is the S-wave velocity of the lower layer, Vp1 is the P-wave velocity of the upper layer, Vp2 is the P-wave velocity of the lower layer, and θ is the incident angle.

8. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 6 or 7, characterized in that: The direct inversion calculation formula of absolute elastic impedance in the frequency domain without well constraints is: where EI′(t,θ) is the absolute elastic impedance without correction of the energy of the seismic wavelet w(t,θ).

9. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 8, characterized in that: The seismic wavelet energy correction factor is obtained using the following steps: Based on the relative elastic impedance expression, obtaining the maximum and minimum values of the relative elastic impedance that is not corrected for seismic wavelet energy within the time period of the formation elastic impedance to be inverted; Based on the elastic impedance low-frequency trend model, obtaining relative elastic impedance of the maximum and minimum values of the low-frequency trend model of the actual underground elastic impedance within the time period of the formation elastic impedance to be inverted; The seismic wavelet energy correction factor is calculated using the following formula: Where a is the seismic wavelet energy correction factor, In EI′ max and In EI′ min In EI is the maximum and minimum relative elastic impedance without correction of seismic wavelet energy. 0max and In EI 0min is the maximum value EI of the low-frequency trend model of the true elastic impedance of the underground 0max and minimum EI 0min The relative elastic impedance.

10. The frequency domain high-resolution direct inversion method for elastic impedance without well constraints according to claim 9, characterized in that: The frequency domain high-resolution elastic impedance direct inversion calculation formula based on low-frequency trend model energy correction is: where EI(t,θ)” is the high-resolution absolute elastic impedance after correction of the seismic wavelet energy based on the low-frequency trend model.

Images