Method and system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in time-frequency domain

Abstract

The present invention belongs to the technical field of seismic data processing for oil and gas exploration, and specifically discloses a method and system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain. The method of the present invention first calculates the time-varying logarithmic amplitude spectrum of seismic data based on the short-time Fourier transform, and then maps the time-varying logarithmic amplitude spectrum of seismic data into a two-dimensional square-integrable function space spanned by Legendre polynomials. In this space, the time-varying logarithmic amplitude spectrum of seismic data is decomposed by two-dimensional Legendre polynomials to estimate the time-varying wavelet amplitude spectrum, and finally the equivalent formation Q value is calculated using the time-varying wavelet amplitude spectrum. The method of the present invention uses Legendre polynomial decomposition to estimate the time-varying wavelet amplitude spectrum. Compared with the methods of decomposing into Fourier series and traditional polynomials, this method has the advantages of fast convergence, strong approximation ability, and reduction of matrix ill-conditioning. It can also eliminate the influence of factors such as local strong reflections and waveform coupling, and improve the stability and accuracy of Q value estimation.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data processing for oil and gas exploration, and particularly relates to a method and system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain. Background Art

[0002] Affected by the non-perfect elastic factors of the formation, the energy of the observed seismic wave signal gradually decays with time, the main frequency decreases and the frequency band gradually narrows. This absorption characteristic of the formation for seismic wave energy is an inherent characteristic of the formation medium, called the quality factor, and is usually represented by the Q value. The Q value of the formation is usually divided into the layer Q value and the formation equivalent Q value, which correspond to different scales and calculation methods, and the two can be converted to each other. The layer Q value is used to describe the absorption attenuation of a single formation, and is suitable for well data analysis and research on the lithological and physical characteristics of local formations. The formation equivalent Q value is used for large-scale average attenuation and is applicable to seismic wave propagation and Q compensation processing.

[0003] The Q values of formations with different lithologies, physical properties, and their pore fillings with different fluids are different. Precise Q values can be used for energy compensation of seismic data, improving the imaging quality of deep formations, enhancing the resolution of seismic data, oil and gas reservoir prediction, lithological feature identification, etc. In the hydrocarbon-bearing area, the absorption attenuation of the reservoir for seismic waves is more obvious, which provides a basis and approach for hydrocarbon detection based on seismic data. Therefore, estimating the stable and high-precision formation Q value has clear practical significance.

[0004] Currently, there are many Q value estimation methods, such as the spectral ratio method, frequency shift method, centroid frequency method, waveform matching method, etc. These methods mostly directly estimate the formation Q value based on seismic data. However, affected by factors such as the signal-to-noise ratio of the observed data, multiple waves, waveform coupling caused by complex reflection coefficient combinations, local strong reflections, local energy focusing or scattering caused by special structures, etc., the spectrum of the reflected wave will be distorted, interfering with the accurate extraction of the amplitude attenuation characteristics, resulting in low accuracy and poor stability of these methods for directly estimating the formation Q value using seismic data.

[0005] From the perspective of the convolution model, the essence of the attenuation of seismic wave energy caused by the formation Q value is to cause the time variation of the seismic wavelet, resulting in the energy and frequency band of the seismic wavelet changing with time. To improve the accuracy of formation Q value estimation, it is necessary to eliminate the influence of the above factors affecting the Q value estimation accuracy as much as possible. Therefore, obtaining a time-varying wavelet or the amplitude spectrum of the time-varying wavelet that changes stably with time is the key to improving the stability and accuracy of Q value estimation. Summary of the Invention

[0006] The object of the present invention is to propose a method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain. This method uses two-dimensional Legendre polynomials to fit a stable time-varying wavelet amplitude spectrum from the time-varying amplitude spectrum of seismic records, and further calculates the formation equivalent Q value from the time-varying wavelet amplitude spectrum, which can improve the stability and accuracy of Q value estimation.

[0007] To achieve the above object, the present invention adopts the following technical solutions:

[0008] The method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain includes the following steps:

[0009] Step 1. Calculate the time-varying logarithmic amplitude spectrum of seismic data based on the short-time Fourier transform;

[0010] Step 2. Map the time-varying logarithmic amplitude spectrum of seismic data into the two-dimensional square-integrable function space spanned by Legendre polynomials;

[0011] Step 3. In the two-dimensional square-integrable function space, perform two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data to estimate the time-varying wavelet amplitude spectrum;

[0012] Step 4. Calculate the formation equivalent Q value using the estimated time-varying wavelet amplitude spectrum.

[0013] In addition, based on the above method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, the present invention also proposes a corresponding system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, which adopts the following technical solutions:

[0014] The system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain includes the following modules:

[0015] A time-varying logarithmic amplitude spectrum calculation module, configured to calculate the time-varying logarithmic amplitude spectrum of seismic data based on the short-time Fourier transform;

[0016] A mapping processing module, configured to map the time-varying logarithmic amplitude spectrum of seismic data into the two-dimensional square-integrable function space spanned by Legendre polynomials;

[0017] An amplitude spectrum estimation module, configured to perform two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data in the two-dimensional square-integrable function space to estimate the time-varying wavelet amplitude spectrum;

[0018] And a formation equivalent Q value estimation module, configured to calculate the formation equivalent Q value using the estimated time-varying wavelet amplitude spectrum.

[0019] In addition, based on the above method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, the present invention also proposes a computer device;

[0020] The computer device includes a memory and one or more processors. Executable code is stored in the memory. When the processor executes the executable code, the steps of the formation Q-value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain as described above are implemented.

[0021] In addition, based on the formation Q-value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain as described above, the present invention also provides a computer-readable storage medium;

[0022] A program is stored on the computer-readable storage medium. When the program is executed by the processor, the steps of the formation Q-value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain as described above are implemented.

[0023] The present invention has the following advantages:

[0024] As described above, the present invention relates to a formation Q-value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain. This method uses two-dimensional Legendre polynomials to fit a stable time-varying wavelet amplitude spectrum from the time-varying amplitude spectrum of seismic records. Compared with the forms of expanding into Fourier series, traditional polynomials and other basis functions, the method of the present invention represents the time-varying wavelet amplitude spectrum in the form of two-dimensional Legendre polynomial expansion, having the advantages of fast convergence, strong approximation ability, reducing Gibbs phenomenon, and reducing matrix ill-conditioning. The method of the present invention also further calculates the formation equivalent Q-value by estimating the time-varying wavelet amplitude spectrum. Compared with the method of directly calculating the formation equivalent Q-value from seismic records, the method of the present invention can eliminate the influence of factors such as reflection coefficients, especially local strong reflections and waveform coupling, and has high estimation accuracy and good stability for the formation equivalent Q-value. Description of the Drawings

[0025] Figure 1 It is a flowchart of the formation Q-value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain in an embodiment of the present invention.

[0026] Figure 2 It is a grayscale display diagram of a time-varying wavelet after constant Q attenuation of a Ricker wavelet with a main frequency of 45 Hz generated in an example of the present invention.

[0027] Figure 3 It is a simulated reflection coefficient sequence diagram generated in an example of the present invention.

[0028] Figure 4 For the Ricker wavelet of 45 Hz and Figure 3 It is a waveform diagram of a synthetic seismic record without attenuation obtained by convolving the reflection coefficient sequence.

[0029] Figure 5 For Figure 2 The time-varying wavelet of Figure 3Waveform diagram of synthetic seismic record of Q attenuation obtained by convolution of reflection coefficient.

[0030] Figure 6 is Figure 5 Waveform diagram of the result after T compensation.

[0031] Figure 7 is Figure 3 Gray-scale display diagram of time-varying amplitude spectrum of reflection coefficient sequence.

[0032] Figure 8 is Figure 2 Gray-scale display diagram of normalized time-varying amplitude spectrum of time-varying wavelet.

[0033] Figure 9 is Figure 5 Gray-scale display diagram of time-varying amplitude spectrum of synthetic seismic record of Q attenuation after energy gain processing in time direction.

[0034] Figure 10 is from Figure 5 Gray-scale display diagram of normalized time-varying amplitude spectrum of time-varying wavelet estimated by using the seismic record after attenuation of with the method of the present invention.

[0035] Figure 11 is Figure 10 Gray-scale display diagram of error of time-varying wavelet amplitude spectrum estimated by using the method of the present invention of .

[0036] Figure 12 is Figure 8 and Figure 10 Comparison diagram of theoretical value and estimated value of instantaneous wavelet at 200 ms and 600 ms in and .

[0037] Figure 13 is from Figure 5 Comparison diagram of estimated formation equivalent Q value and theoretical value by using the seismic record after attenuation of with the method of the present invention.

[0038] Figure 14 Waveform display diagram of synthetic seismic record without attenuation in the example of the present invention.

[0039] Figure 15 Waveform display diagram of synthetic seismic record of Q attenuation in the example of the present invention.

[0040] Figure 16 Waveform display diagram of the result after Q compensation by using the Q value estimated by the method of the present invention in the example of the present invention.

[0041] Figure 17 is in the example of the present invention Figure 16 compensation result and Figure 14 difference waveform display diagram of synthetic seismic record without attenuation.

[0042] Figure 18 This is the grayscale display map of the post-stack actual seismic data in the embodiments of the present invention.

[0043] Figure 19 This is the grayscale display map of the time-varying amplitude spectrum obtained by performing short-time Fourier transform on the 200th post-stack actual seismic trace.

[0044] Figure 20 This is the grayscale display map of the time-varying wavelet amplitude spectrum estimated by using the method of the present invention from the 200th post-stack actual seismic trace.

[0045] Figure 21 This is the amplitude spectrum map of the source wavelet estimated by using the method of the present invention from the post-stack actual seismic data.

[0046] Figure 22 This is the grayscale display map of the result of the formation equivalent Q-value estimated by using the method of the present invention from the post-stack actual seismic data. Detailed implementation manners

[0047] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners:

[0048] Embodiment 1

[0049] This Embodiment 1 describes a method for estimating the formation Q-value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, so as to achieve the purpose of stably and highly accurately estimating the formation equivalent Q-value of seismic data. The method of the present invention uses two-dimensional Legendre polynomials to fit a stable time-varying wavelet amplitude spectrum from the time-varying amplitude spectrum of the seismic record, and further calculates the formation equivalent Q-value from the time-varying wavelet amplitude spectrum. Compared with the forms of basis functions such as expanding into Fourier series and traditional polynomials, the method based on two-dimensional Legendre polynomials has the advantages of fast convergence, strong approximation ability, reducing Gibbs phenomenon, and reducing matrix ill-conditioning. Compared with the method of directly calculating the formation equivalent Q-value from the seismic record, the method of the present invention can eliminate the influence of reflection coefficients, especially local strong reflections and waveform coupling, on the Q-value estimation result, thereby improving the accuracy and stability of the Q-value estimation.

[0050] As Figure 1 shown, the method for estimating the formation Q-value by two-dimensional Legendre polynomial decomposition in the time-frequency domain in this embodiment specifically includes the following steps:

[0051] Step 1. Calculate the time-varying logarithmic amplitude spectrum of the seismic data based on short-time Fourier transform.

[0052] In step 1, first perform a short-time Fourier transform on the seismic data to obtain the time-frequency spectrum of the seismic data, then take the modulus of the time-frequency spectrum of the seismic data to obtain the time-varying amplitude spectrum of the seismic data, and finally take the logarithm of the time-varying amplitude spectrum of the seismic data to obtain the time-varying logarithmic amplitude spectrum of the seismic data. The specific steps are as follows:

[0053] Perform a short-time Fourier transform on the seismic data to obtain the time-frequency spectrum of the seismic data, denoted as , then:

[0054] (1)

[0055] where is the observation time of the seismic record, is the time shift parameter, is the frequency, is the imaginary unit, represents the kernel function of the short-time Fourier transform.

[0056] Take the modulus of the time-frequency spectrum of the seismic data to obtain the time-varying amplitude spectrum of the seismic data, denoted as , that is:

[0057] (2)

[0058] where represents taking the modulus of a complex number.

[0059] Take the logarithm of the time-varying amplitude spectrum of the seismic data to obtain the time-varying logarithmic amplitude spectrum of the seismic data, calculated using the following formula:

[0060] (3)

[0061] where is a non-negative small constant, the role of which is to avoid the situation where taking the logarithm of a zero value in the data results in negative infinity.

[0062] In addition, in step 1 of the method for estimating the formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain of this embodiment, the kernel function of the short-time Fourier transform uses a Gaussian function, and its expression is:

[0063] (4)

[0064] where is the standard deviation, and the standard deviation determines the width of the window.

[0065] Step 2. Map the time-varying logarithmic amplitude spectrum of the seismic data into the two-dimensional square-integrable function space spanned by Legendre polynomials.

[0066] In Step 2, the time variable and frequency variable of the time-varying logarithmic amplitude spectrum of the seismic data in the two-dimensional real space are mapped into the two-dimensional square-integrable function space spanned by Legendre polynomials through translation and scaling processing. The specific steps are as follows:

[0067] In the square-integrable function space , Legendre polynomials form an orthogonal basis, which can span an orthogonal space. Any square-integrable function in this orthogonal space can be expanded in the form of an orthogonal decomposition of Legendre polynomials. If the logarithmic amplitude spectrum of the time-varying wavelet is expressed in the form of an orthogonal decomposition of two-dimensional Legendre polynomials, then the discrete time variable, i.e., the time translation parameter and the frequency variable, i.e., the frequency need to be mapped onto through translation and scaling processing.

[0068] Let the discrete time variable of the time-varying logarithmic amplitude spectrum of the seismic data, i.e., the time translation parameter , have a minimum value of 0 ms, and the maximum value of the time variable be . Then any time variable within the time range can be mapped onto the two-dimensional square-integrable function space spanned by Legendre polynomials through formula (5):

[0069] (5)

[0070] where represents the value of the time variable after mapping the time variable of the two-dimensional real space, i.e., the time translation parameter , onto the two-dimensional square-integrable function space .

[0071] When estimating the logarithmic amplitude spectrum of the time-varying wavelet, let the minimum value of the frequency variable of the time-varying logarithmic amplitude spectrum of the seismic data, i.e., the frequency , be 0 Hz, and the maximum value of the frequency be . Then any frequency variable within the frequency range can be mapped onto the two-dimensional square-integrable function space spanned by Legendre polynomials through formula (6):

[0072] (6)

[0073] Among them, represents the value of the frequency variable after being mapped from the frequency variable in the two-dimensional real number space to the two-dimensional square-integrable function space .

[0074] Step 3. In the two-dimensional square-integrable function space, perform two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of the seismic data to estimate the time-varying wavelet amplitude spectrum.

[0075] In Step 3, represent the logarithmic amplitude spectrum of the time-varying wavelet in the form of a two-dimensional Legendre polynomial, and fit the time-varying logarithmic amplitude spectrum of the seismic data in the least squares sense, so as to obtain the coefficients of each term of the two-dimensional Legendre polynomial. Substitute the coefficients into the two-dimensional Legendre polynomial to obtain the logarithmic amplitude spectrum of the fitted time-varying wavelet, and then perform exponential processing on the processing result to obtain the time-varying wavelet amplitude spectrum. The specific steps are as follows:

[0076] Step 3.1. Orthogonal decomposition of the Legendre polynomial of the logarithmic amplitude spectrum of the time-varying wavelet.

[0077] In the square-integrable function space , any square-integrable function can be expanded into the form of orthogonal decomposition of the Legendre polynomial. Therefore, write the logarithmic amplitude spectrum of the time-varying wavelet in the form of orthogonal decomposition of the two-dimensional Legendre polynomial as shown in formula (7):

[0078] (7)

[0079] Among them, represents the logarithmic amplitude spectrum of the time-varying wavelet, is the -th Legendre polynomial, is the -th Legendre polynomial, represents the Legendre polynomial coefficient, , , and are the orders of orthogonal decomposition of the Legendre polynomial in the time and frequency directions respectively.

[0080] The mathematical expression of

[0081] is:

[0082] Among them, , is the floor operation.

[0083] The mathematical expression is:

[0084] (9)

[0085] Wherein, .

[0086] Step 3.2. Solve the coefficients of the two-dimensional Legendre polynomial by least squares fitting.

[0087] On the premise of the preset Legendre polynomial order and , use the mathematical expression of the logarithmic amplitude spectrum of the time-varying wavelet shown in formula (7) to approximate the logarithmic amplitude spectrum of the seismic data . By least squares fitting, find the optimal Legendre polynomial coefficients , and calculate the coefficients of each Legendre polynomial .

[0088] Step 3.3. Calculate the amplitude spectrum of the time-varying wavelet.

[0089] Substitute the calculated values of the coefficients of each order of the Legendre polynomial into formula (7) to obtain the logarithmic amplitude spectrum of the time-varying wavelet , and take exponential processing to obtain the amplitude spectrum of the time-varying wavelet , and its calculation formula is:

[0090] (10)

[0091] Substitute formula (5), formula (6), formula (8) and formula (9) into formula (10) to calculate the amplitude spectrum of the time-varying wavelet.

[0092] In addition, in step 3 of the method for estimating the formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain of this embodiment, the process of calculating the Legendre polynomial coefficients is specifically as follows:

[0093] Solve the coefficients of the two-dimensional Legendre polynomial by least squares fitting to minimize the following sum of squared residuals:

[0094] (11)

[0095] Wherein, represents the error function; take the partial derivative of the error function with respect to each Legendre polynomial coefficient :

[0096] (12)

[0097] Among them, , .

[0098] Simplifying formula (12) gives:

[0099] (13)

[0100] For the sake of convenient writing:

[0101] Let be denoted as ;

[0102] Let be denoted as .

[0103] Then let:

[0104] .

[0105] , .

[0106] Rewrite formula (13) into the matrix form as shown in formula (14):

[0107] (14)

[0108] Among them, is the matrix of calculated from Legendre polynomials, is the coefficient matrix of to be solved, is vector of.

[0109] Solving formula (14) gives the coefficient matrix of Legendre polynomials as:

[0110] (15).

[0111] Calculating formula (15) gives the coefficients of each Legendre polynomial .

[0112] Step 4. Calculate the formation equivalent Q value using the estimated time-varying wavelet amplitude spectrum.

[0113] In step 4, the instantaneous wavelet amplitude spectrum with the widest frequency band at the first arrival time of the reflected wave in the time-varying wavelet amplitude spectrum is used as the source wavelet amplitude spectrum. If it is multi-channel data, it needs to be statistically averaged for multi-channel data and then used as the amplitude spectrum of the source wavelet. Then, combined with the Q-value definition formula, a linear equation for estimating the equivalent Q-value is obtained. Next, the frequency band range for calculating the equivalent Q-value estimate is preset, and the least squares method is used to solve the slope of the linear equation, and further calculate the equivalent Q-value as the final result for output.

[0114] Specifically, based on obtaining the time-varying wavelet amplitude spectrum of seismic data in step 3, the following method is used to calculate the formation equivalent Q-value. The instantaneous wavelet amplitude spectrum with the widest frequency band near the initial time in the time-varying wavelet amplitude spectrum is used as the source wavelet amplitude spectrum. If the instantaneous wavelet amplitude spectrum with the widest frequency band near the initial time in the time-varying wavelet amplitude spectrum is multi-channel data, the following multi-channel statistical weighting process needs to be performed. Denote the first arrival time of the reflected wave as , and the trace number is denoted by , then the source wavelet after multi-channel statistics is expressed as:

[0115] (16)

[0116] Where is the instantaneous wavelet amplitude spectrum estimated from the th trace data at the moment, , is the total number of traces.

[0117] Taking as the initial time for estimating the formation equivalent Q-value, using the Q-value definition formula, the amplitude spectrum after the wave propagation time of is expressed as:

[0118] (17)

[0119] Where corresponds to the equivalent Q-value describing the average absorption effect of the formation.

[0120] Rewrite formula (17) into a linear equation as shown in formula (18):

[0121] (18)

[0122] Select the frequency band range for calculating the formation equivalent Q-value and denote it as , where is the total number of sample points of the selected frequency, represents the starting frequency for equivalent Q-value estimation, Denote the cut-off frequency for equivalent Q-value estimation. The slope of the linear equation can be solved using the least squares method , that is:

[0123] (19)

[0124] After further simplification, the formula for solving the formation equivalent Q-value as shown in formula (20) is obtained:

[0125] (20)

[0126] Substitute the calculated time-varying wavelet amplitude spectrum and the source wavelet into formula (20), then the time-varying formation equivalent Q-value can be calculated, and the calculated formation equivalent Q-value is output as the final result of the formation equivalent Q-value estimation.

[0127] The method of the present invention uses two-dimensional Legendre polynomials to fit a stable time-varying wavelet amplitude spectrum from the time-varying amplitude spectrum of seismic records. Compared with the forms of basis functions such as expanding into Fourier series and traditional polynomials, the method based on two-dimensional Legendre polynomials has the advantages of fast convergence, strong approximation ability, reducing Gibbs phenomenon, and reducing matrix ill-conditioning. In addition, this method further calculates the formation equivalent Q-value by estimating the time-varying wavelet of seismic records. Compared with the method of directly calculating the formation Q-value from seismic records, it can eliminate the influence of reflection coefficients, especially local strong reflections and waveform coupling, etc. on the Q-value estimation result, and improve the stability and accuracy of Q-value estimation.

[0128] To verify the effectiveness of the method proposed in the present invention, the present invention also gives two specific experiments:

[0129] In specific experiment 1 and experiment 2, the formation Q-value estimation method of two-dimensional Legendre polynomial decomposition in the time-frequency domain of the present invention is described by combining theoretical synthetic seismic records and post-stack actual seismic data respectively.

[0130] Experiment 1

[0131] Figure 2 uses the Ricker wavelet with a dominant frequency of 45 Hz, that is, the Ricker wavelet, as the source wavelet. When the formation equivalent Q-value is 40, the time-varying wavelet calculated using the Q-value definition formula. Figure 3 is the generated Gaussian-Bernoulli reflection coefficient sequence. Convolve the Ricker wavelet with a dominant frequency of 45 Hz with this reflection coefficient sequence to obtain a stationary synthetic seismic record as shown in Figure 4 . Use the time-varying wavelet of Figure 2 to perform time-varying convolution processing with the reflection coefficient sequence of Figure 3 to obtain a synthetic seismic record with Q attenuation as shown in Figure 5 .

[0132] Figures 3 to 6 It is a comparison chart of synthetic seismograms before and after Q attenuation and the results of T compensation. It can be seen from the comparison between Figure 5 and Figure 4 that after the Q value of the formation decays, the amplitude of the synthetic seismogram gradually decreases, and the energy of the seismic wave is lost. For the synthetic seismogram with the attenuation of Figure 5 , T compensation processing is carried out, and the results are as shown in Figure 6 . Theoretically speaking, T compensation processing only compensates the amplitude of the seismic wave and does not change the bandwidth of the data. It can be seen from Figure 6 that after T compensation, the amplitude of the attenuated synthetic seismogram is effectively restored. Further comparison between Figure 6 and Figure 4 shows that although the amplitude of the synthetic seismogram is effectively restored after T compensation processing, the resolution is still relatively low, such as at 700 to 900 ms in Figure 6 . In other words, Q attenuation not only reduces the amplitude of the seismic wave, but also changes the frequency band of the seismic data and reduces the resolution of the seismic data. To compensate or eliminate the effects of amplitude reduction and bandwidth narrowing caused by the absorption of seismic wave energy by the formation, it is necessary to accurately estimate the Q value of the formation for Q compensation processing. Therefore, it is very necessary to stably and accurately estimate the Q value of the formation.

[0133] Figure 7 gives Figure 3 the time-varying amplitude spectrum obtained by short-time Fourier transform of the reflection coefficient sequence. It can be seen from the figure that the overall trend of the time-varying amplitude spectrum of the reflection coefficient sequence changing with time and frequency is relatively stable, but there is strong amplitude energy in the reflection coefficient locally in time. Figure 8 is Figure 2 the normalized display of the time-varying amplitude spectrum obtained by short-time Fourier transform of the time-varying wavelet. It can be seen from Figure 8 that after the action of the Q value, the main frequency of the seismic wavelet gradually decreases with time, and the high-frequency energy decays severely. Figure 9 is Figure 5 the time-varying amplitude spectrum obtained by short-time Fourier transform of the synthetic record with Q attenuation, and for the convenience of comparing the frequency band changes in the time direction, energy gain processing is performed in the time direction. It can be seen from Figure 9 that the oscillation of the energy strength change of the time-varying amplitude spectrum of the synthetic seismogram in time is very serious, which greatly affects the stability and accuracy of Q value estimation by traditional methods.

[0134] Figure 10 is the time-varying wavelet amplitude spectrum estimated by the method of the present invention from the synthetic seismogram after Q attenuation of Figure 5 , where the parameters selected for two-dimensional Legendre polynomial decomposition are respectively and 。It can be seen from the figure that the estimated time-varying wavelet amplitude spectrum changes slowly with time, the frequency band gradually narrows with time, the main frequency gradually decreases with time, and the high-frequency attenuation is severe. The estimated time-varying wavelet amplitude spectrum is very stable and is not affected by factors such as local strong reflections and waveform coupling in the seismic record. Comparing Figure 10 the estimated value with Figure 8 the theoretical value, it can be seen that the two are basically exactly the same. Further comparing the difference between the estimated value and the theoretical value of the time-varying wavelet amplitude spectrum, as shown in Figure 11 , it can be seen from the figure that the difference between the two is not obvious.

[0135] Further analyzing the estimation accuracy of the time-varying wavelet amplitude spectrum, Figure 12 comparison diagrams of the estimated value and the theoretical value of the instantaneous wavelet amplitude spectrum at 200 ms and 600 ms are respectively given. It can be seen from Figure 12 that the estimated values of the main frequencies of the instantaneous wavelet amplitude spectra at 200 ms and 600 ms are basically exactly the same as the theoretical values, and the bandwidths of the estimated values of the instantaneous wavelet amplitude spectra are also basically the same as the theoretical values. Compared with the bandwidth of the estimated value of the instantaneous wavelet amplitude spectrum at 200 ms, the bandwidth of the estimated value of the instantaneous wavelet amplitude spectrum at 600 ms is closer to the theoretical value. The error of the estimated value of the instantaneous wavelet amplitude spectrum at 200 ms mainly appears at the high-frequency end, while the error of the estimated value of the instantaneous wavelet amplitude spectrum at 600 ms mainly appears in the frequency band range of the weak amplitude area. In short, the estimation accuracies of the instantaneous wavelet amplitude spectra at 200 ms and 600 ms are relatively high. Comprehensive comparison of Figure 9 , Figure 10 , Figure 11 and Figure 12 shows that the method of the present invention uses two-dimensional Legendre polynomial fitting to estimate the instantaneous wavelet amplitude spectrum, which can eliminate the influence of local strong reflection coefficients and waveform coupling on the estimation of the instantaneous wavelet amplitude spectrum, the estimation result is relatively stable, and the estimation accuracy of the instantaneous wavelet is relatively high.

[0136] Figure 13 shows a comparison diagram of the formation equivalent Q value result calculated from the estimated time-varying wavelet amplitude spectrum and the theoretical value. It can be seen from the figure that the estimated result is very close to the theoretical value, the error between the two is small, and it is not affected by local strong reflection factors near 820 ms as in Figure 10 , and the estimated result is very stable. Figure 9 Figures 14 to 17

[0137] Figures 14 to 17 is a comparison diagram for Q compensation processing using the estimated formation equivalent Q value, where Figure 14 is a stationary synthetic seismic record, Figure 15 is the synthetic seismic record after Q decay under the condition of the same reflection sequence. The Q value of this data is estimated using the method of the present invention, and the estimated Q is used for Figure 15Perform Q compensation processing to restore the lost amplitude and frequency band. The result of the Q compensation processing is as Figure 16 shown. Comparing Figure 16 the seismic record after Q compensation with Figure 14 the stationary synthetic seismic record, it can be seen that the two are basically completely consistent. Figure 17 gives Figure 16 and Figure 14 the difference. From the figure, it can be seen that the waveform difference on the synthetic seismic record due to the Q estimation error is very small, which also indicates that the Q value estimation result of the method of the present invention has high accuracy and stability.

[0138] Experiment 2

[0139] Figure 18 is the gray-scale display map of the post-stack actual data of a certain exploration area. From the figure, it can be seen that its signal-to-noise ratio is relatively high. Generally speaking, the reflection wave isophase axis above 2.5 seconds is relatively thin, indicating that the resolution of the data is relatively high. Below 2.5 seconds, the reflection wave isophase axis of the data is relatively thick, indicating that the resolution of the data below 2.5 seconds is relatively low. Further analyze the resolution of the sedimentary area enclosed by the black dotted line in the range of 2.0 to 3.0 s in the figure. From the figure, it can be seen that compared with the surrounding wall rock, the reflection wave isophase axis of the area enclosed by the black dotted line is relatively thin, indicating that the resolution of the data is relatively high, while the reflection wave isophase axis of the surrounding wall rock is relatively thick and the resolution of the data is relatively low.

[0140] The formation equivalent Q value is calculated by combining the instantaneous wavelet amplitude spectrum and the source wavelet amplitude spectrum with the travel time. The relatively high resolution of the area enclosed by the black dotted line to a certain extent indicates that the main frequency of this area is relatively high and the frequency band width is relatively wide, and its corresponding formation equivalent Q value should be relatively large. While the resolution of the surrounding wall rock is relatively low, it can be inferred that the main frequency of the reflection wave isophase axis of the wall rock is relatively low and the frequency band width is relatively narrow, and its corresponding formation equivalent Q value should be relatively small.

[0141] Figure 19 is Figure 18 the time-varying amplitude spectrum obtained by short-time Fourier transform processing of the 200th post-stack actual seismic data in Figure 20 is the time-varying wavelet amplitude spectrum estimated for this seismic data by using the method of the present invention. Comparing Figure 19 and Figure 20 it can be seen that the main frequency and frequency band width of the two are basically the same, and the time-varying wavelet amplitude spectrum eliminates the influence of factors such as local strong reflections and waveform coupling in the time-varying amplitude spectrum of the seismic data, indicating that the time-varying wavelet amplitude spectrum estimated by the method of the present invention has high accuracy and stability, and can accurately reflect the changes of the main frequency and frequency band width of the seismic data with time, which is the premise of accurate Q value estimation.

[0142] Since the post-stack data is generally the result of pre-stack resolution improvement processing and migration, it is difficult to estimate the source wavelet in a strict sense. In the present invention, the result of multi-channel statistical averaging of the instantaneous wavelets estimated from the actual seismic data within the range of 500 ms to 600 ms is used as the source wavelet, and its amplitude spectrum is as Figure 21 shown. It can be seen from Figure 21 that the main frequency of the amplitude spectrum of the estimated source wavelet is about 38 Hz, and it has a certain frequency band width.

[0143] Using the amplitude spectrum of the instantaneous wavelet estimated from the post-stack actual seismic data by the method of the present invention, and combining the travel time and the estimated source wavelet to estimate the formation equivalent Q value, the estimation result is as Figure 22 shown. It can be seen from Figure 22 that as time increases, the overall Q value of the estimated seismic data shows a decreasing trend with time, and there are also obvious differences in the Q values of different geological sedimentary regions.

[0144] The formation equivalent Q value of the area surrounded by the black dashed line in Figure 18 analyzed from the seismic section mentioned above is much higher than the formation equivalent Q value of the surrounding wall rock, which is consistent with the Figure 22 estimated result of the formation equivalent Q value in Figure 22 The formation equivalent Q value of the area surrounded by the black dashed line box is much higher than the formation equivalent Q value of the surrounding wall rock. It can also be seen from the comparison and analysis of the formation equivalent Q value estimation result in Figure 22 and the Figure 18 section details in

[0145] Example 2

[0146] This Example 2 describes a formation Q value estimation system for two-dimensional Legendre polynomial decomposition in the time-frequency domain. This system is based on the same inventive concept as the formation Q value estimation method for two-dimensional Legendre polynomial decomposition in the above Example 1.

[0147] The formation Q value estimation system for two-dimensional Legendre polynomial decomposition in the time-frequency domain includes the following modules:

[0148] A time-varying logarithmic amplitude spectrum calculation module, which is used to calculate the time-varying logarithmic amplitude spectrum of seismic data based on the short-time Fourier transform.

[0149] A mapping processing module, configured to map the time-varying logarithmic amplitude spectrum of seismic data into a two-dimensional square-integrable function space spanned by Legendre polynomials.

[0150] An amplitude spectrum estimation module, configured to perform two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data in the two-dimensional square-integrable function space to estimate the time-varying wavelet amplitude spectrum.

[0151] And a formation equivalent Q value estimation module, configured to calculate the formation equivalent Q value by using the estimated time-varying wavelet amplitude spectrum.

[0152] It should be noted that in the formation Q value estimation system with two-dimensional Legendre polynomial decomposition in the time-frequency domain, the implementation processes of the functions and roles of each functional module are specifically detailed in the corresponding steps of the method in the above-mentioned Embodiment 1, and will not be elaborated here.

[0153] Embodiment 3

[0154] This Embodiment 3 describes a computer device, which includes a memory and one or more processors. An executable code is stored in the memory, and when the processor executes the executable code, it is used to implement the steps of the formation Q value estimation method with two-dimensional Legendre polynomial decomposition in the above-mentioned Embodiment 1.

[0155] In this embodiment, the computer device is any device or apparatus with data processing capabilities, which will not be elaborated here.

[0156] Embodiment 4

[0157] This Embodiment 4 describes a computer-readable storage medium, on which a program is stored. When the program is executed by a processor, it is used to implement the steps of the formation Q value estimation method with two-dimensional Legendre polynomial decomposition in Embodiment 1. The computer-readable storage medium can be an internal storage unit of any device or apparatus with data processing capabilities, such as a hard disk or a memory, or an external storage device of any device with data processing capabilities, such as a plug-in hard disk, a Smart Media Card (SMC), an SD card, a Flash Card, etc. equipped on the device.

[0158] Of course, the above description is only the preferred embodiments of the present invention. The present invention is not limited to listing the above embodiments. It should be noted that all equivalent substitutions and obvious deformation forms made by any person skilled in the art under the teaching of this specification fall within the substantial scope of this specification and should be protected by the present invention.

Claims

1. A method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in time-frequency domain, characterized in that It includes the following steps: Step 1. Calculate the time-varying logarithmic amplitude spectrum of seismic data based on the short-time Fourier transform; Step 2. Map the time-varying logarithmic amplitude spectrum of seismic data into the two-dimensional square-integrable function space spanned by Legendre polynomials; Step 3. In the two-dimensional square-integrable function space, perform two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data to estimate the time-varying wavelet amplitude spectrum; Step 4. Calculate the formation equivalent Q value using the estimated time-varying wavelet amplitude spectrum; The specific content of step 3 is as follows: Write the logarithmic amplitude spectrum of the time-varying wavelet in the form of two-dimensional Legendre polynomial orthogonal decomposition shown in formula (7): Among them, G(τ,f) represents the logarithmic amplitude spectrum of the time-varying wavelet, and P m (T(τ)) is the Legendre polynomial of degree m, and P n (F(f)) is the Legendre polynomial of degree n, and a mn represents the Legendre polynomial coefficient, m ∈ [0,M], n ∈ [0,N], and M and N are the orders of the orthogonal decomposition of the Legendre polynomial in the time and frequency directions respectively; P m (T(τ)) has the following mathematical expression: wherein, is a floor operation; P n (F(f)) has the following mathematical expression: Among them, On the premise of the preset orders M and N of Legendre polynomial orthogonal decomposition, use the mathematical expression of the logarithmic amplitude spectrum of the time-varying wavelet shown in formula (7) to approximate the time-varying logarithmic amplitude spectrum B(τ, f) of seismic data; Solve for the coefficients a of the two-dimensional Legendre polynomial using least squares fitting mn , and substitute each calculated coefficient a of the Legendre polynomial mn into formula (7) to obtain the logarithmic amplitude spectrum G(τ, f) of the time-varying wavelet, and take the exponential of G(τ, f) to obtain the amplitude spectrum W(τ, f) of the time-varying wavelet. Its calculation formula is as follows: Substitute formula (5), formula (6), formula (8) and formula (9) into formula (10) to calculate the time-varying wavelet amplitude spectrum.

2. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain according to claim 1, characterized in that, The specific content of step 1 is as follows: Perform short-time Fourier transform on the seismic data x(t) to obtain the time-frequency spectrum of the seismic data denoted as P(τ, f): where t is the observation time of the seismic record, τ is the time shift parameter, f is the frequency, j is the imaginary unit, and w(t - τ) represents the kernel function of the short-time Fourier transform; Take the modulus of the time-frequency spectrum of the seismic data to obtain the time-varying amplitude spectrum of the seismic data denoted as A(τ, f): A(τ, f) = |P(τ, f)| (2) where |·| represents taking the modulus of a complex number; Take the logarithm of the time-varying amplitude spectrum A(τ, f) of the seismic data to obtain the time-varying logarithmic amplitude spectrum B(τ, f) of the seismic data: B(τ, f) = ln(A(τ, f) + α) (3) where α is a non-negative constant.

3. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain according to claim 2, characterized in that, In step 1, the kernel function w(t - τ) of the short-time Fourier transform uses a Gaussian function, and its expression is: where σ is the standard deviation.

4. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain according to claim 1, characterized in that The specific content of step 2 is as follows: Let the minimum value of the discrete-time variable of the time-varying logarithmic amplitude spectrum B(τ, f) of seismic data, i.e., the time-shift parameter τ, be 0 ms, and the maximum value of the time variable τ be τ max , then any time variable τ within [0, τ max can be mapped to the two-dimensional square-integrable function space L 2 ([-1, 1]) by formula (5): Among them, T(τ) represents the value of the time variable after mapping the time variable τ in the two-dimensional real number space to the two-dimensional square-integrable function space L 2 ([-1, 1]); Let the minimum value of the frequency variable of the time-varying logarithmic amplitude spectrum B(τ, f) of seismic data, i.e., the frequency f, be 0 Hz, and the maximum value of the frequency f be f max , then any frequency variable f within the range [0, f max can be mapped to the two-dimensional square-integrable function space L 2 ([-1, 1]) by formula (6): Among them, F(f) represents the value of the frequency variable after mapping the frequency variable f in the two-dimensional real number space to the two-dimensional square-integrable function space L 2 ([-1,1]).

5. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain according to claim 4, characterized in that, In step 3, the process of calculating the Legendre polynomial coefficient a mn is specifically as follows: Solving for the coefficients a of the two-dimensional Legendre polynomial by least squares fitting mn , such that the sum of the squares of the residuals as shown in formula (11) is minimized: where E represents the error function; taking the partial derivative of the error function E with respect to each Legendre polynomial coefficient a pq : where p ∈ [0, M], q ∈ [0, N]; Simplify formula (12) to get: Respectively, is denoted as Σ τ,f P m (T)P n (F)P p (T)P q (F), Then let: Rewrite formula (13) into the matrix form shown in formula (14): CA = D (14) where C is a (M + 1)(N + 1) × (M + 1)(N + 1) matrix calculated from Legendre polynomials, A is a (M + 1)(N + 1) × 1 coefficient matrix to be solved, and D is a (M + 1)(N + 1) × 1 vector; Solve formula (14) to obtain the coefficient matrix of Legendre polynomials as: A = C -1 D (15) Obtain each Legendre polynomial coefficient a using calculation formula (15). mn .

6. The formation Q-value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain according to claim 1, characterized in that The specific content of step 4 is as follows: Take the instantaneous wavelet amplitude spectrum with the widest frequency band at the arrival time of the reflected wave in the time-varying wavelet amplitude spectrum as the source wavelet amplitude spectrum; Perform multi-channel statistical averaging on multi-channel data. Denote the arrival time of the reflected wave as t1 and the trace number as s. Then the source wavelet amplitude spectrum W(t1, f) after multi-channel statistics is expressed as: where f is the frequency, W s (t1, f) is the instantaneous wavelet amplitude spectrum at time t1 estimated from the data of the s-th trace, s ∈ [1, l], and l is the total number of traces; Take t1 as the initial time for estimating the formation equivalent Q value. Using the Q value definition formula, the source wavelet amplitude spectrum W(t1 + Δt, f) after the wave propagation time of Δt is expressed as: where Q is the equivalent Q value describing the average absorption effect of the formation; Rewrite Equation (17) into a linear equation as shown in Equation (18): The frequency band range selected for calculating the equivalent Q value of the formation is denoted as [f1, f v , where v is the total number of sample points of the selected frequency, f1 represents the starting frequency for equivalent Q value estimation, and f v represents the cut-off frequency for equivalent Q value estimation. The slope of the linear equation is solved using the least squares method That is: Furthermore, obtain the formula for solving the equivalent Q value of the formation as shown in Equation (20): Substitute the calculated time-varying wavelet amplitude spectrum W(τ,f) and the source wavelet amplitude spectrum W(t1,f) into Equation (20), calculate the time-varying equivalent Q value of the formation, and output the calculated equivalent Q value of the formation as the final estimation result.

7. A formation Q value estimation system for two-dimensional Legendre polynomial decomposition in the time-frequency domain for implementing the formation Q value estimation method of two-dimensional Legendre polynomial decomposition in the time-frequency domain as described in Claim 1, characterized in that The formation Q value estimation system for two-dimensional Legendre polynomial decomposition in the time-frequency domain includes the following modules: A time-varying logarithmic amplitude spectrum calculation module for calculating the time-varying logarithmic amplitude spectrum of seismic data based on the short-time Fourier transform; A mapping processing module for mapping the time-varying logarithmic amplitude spectrum of seismic data into the two-dimensional square-integrable function space spanned by Legendre polynomials; An amplitude spectrum estimation module for performing two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data in the two-dimensional square-integrable function space to estimate the time-varying wavelet amplitude spectrum; And a formation equivalent Q value estimation module for calculating the formation equivalent Q value using the estimated time-varying wavelet amplitude spectrum.

8. A computer device, comprising a memory and one or more processors, wherein executable code is stored in the memory, characterized in that, When the processor executes the executable code, it implements the steps of the formation Q value estimation method of two-dimensional Legendre polynomial decomposition in the time-frequency domain as described in any one of Claims 1 to 6.

9. A computer-readable storage medium having a program stored thereon, characterized in that, When the program is executed by the processor, it implements the steps of the formation Q value estimation method of two-dimensional Legendre polynomial decomposition in the time-frequency domain as described in any one of Claims 1 to 6.

Images