Stratum Q value estimation method and system based on time-frequency domain two-dimensional Legendre polynomial decomposition

Abstract

The invention belongs to the technical field of oil-gas exploration seismic data processing, and particularly discloses a stratum Q value estimation method and system based on time-frequency domain two-dimensional Legendre polynomial decomposition. The method comprises the following steps: firstly, calculating a time-varying logarithmic amplitude spectrum of seismic data based on short-time Fourier transform, and then mapping the time-varying logarithmic amplitude spectrum of the seismic data into a two-dimensional square integrable function space spanned by a Legendre polynomial; and performing two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of the seismic data to estimate a time-varying wavelet amplitude spectrum, and finally calculating a stratum equivalent Q value by using the time-varying wavelet amplitude spectrum. According to the method, the Legendre polynomial decomposition is used for estimating the time-varying wavelet amplitude spectrum, and compared with a method for decomposing the amplitude spectrum into Fourier series and a traditional polynomial, the method has the advantages of being fast in convergence, high in approximation capability, capable of reducing matrix morbidity and the like, can eliminate influences of factors such as local strong reflection and waveform coupling, and improves stability and precision of Q value estimation.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data processing in 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 becomes narrower. This absorption characteristic of the formation for seismic wave energy is an inherent characteristic of the formation medium, called the quality factor, which is commonly 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 into 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 lithology and physical properties 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 and other processes.

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

[0004] At present, 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 accuracy of Q value estimation 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: A method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, comprising 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.

[0008] 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 system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, which adopts the following technical solutions: A system for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain, comprising the following modules: 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; 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; 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; And a formation equivalent Q value estimation module, configured to calculate the formation equivalent Q value using the estimated time-varying wavelet amplitude spectrum.

[0009] 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; The computer device includes a memory and one or more processors. An executable code is stored in the memory. When the processor executes the executable code, the steps of the method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain as described above are implemented.

[0010] 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 provides a computer-readable storage medium; A program is stored on the computer-readable storage medium. When the program is executed by a processor, the steps of the method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain as described above are implemented.

[0011] The present invention has the following advantages: As described above, the present invention relates to 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. Compared with the forms of basis functions such as Fourier series expansion and traditional polynomials, the method of the present invention represents the time-varying wavelet amplitude spectrum in the form of two-dimensional Legendre polynomial expansion, which has the advantages of fast convergence, strong approximation ability, reduction of Gibbs phenomenon, and reduction of 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. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 It is a flowchart of the method for estimating formation Q value by two-dimensional Legendre polynomial decomposition in the time-frequency domain in an embodiment of the present invention.

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

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

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

[0016] Figure 5 It is Figure 2 the waveform diagram of the synthetic seismic record with Q attenuation obtained by convolving the time-varying wavelet and Figure 3 the reflection coefficient.

[0017] Figure 6 It is Figure 5 the waveform diagram of the result after T compensation.

[0018] Figure 7 It is Figure 3Gray-scale display map of the time-varying amplitude spectrum of the reflection coefficient sequence.

[0019] Figure 8 is Figure 2 Gray-scale display map of the normalized time-varying amplitude spectrum of the time-varying wavelet.

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

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

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

[0023] Figure 12 is Figure 8 and Figure 10 Comparison graph of the theoretical value and the estimated value of the instantaneous wavelet at 200 ms and 600 ms in

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

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

[0026] Figure 15 Waveform display map of the synthetic seismic record with Q attenuation in the example of the present invention.

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

[0028] Figure 17 is in the example of the present invention Figure 16 compensation result and Figure 14 Difference waveform display map between the synthetic seismic record without attenuation.

[0029] Figure 18 Gray-scale display map of the post-stack actual seismic data in the example of the present invention.

[0030] Figure 19 Gray-scale display map of the time-varying amplitude spectrum obtained by performing short-time Fourier transform on the 200th trace of the post-stack actual seismic data.

[0031] Figure 20 This is a grayscale display plot of the time-varying wavelet amplitude spectrum estimated by using the method of the present invention for the 200th post-stack actual seismic data.

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

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

[0034] The present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners: Example 1 This Example 1 describes a method for estimating the formation Q-value by two-dimensional Legendre polynomial decomposition in the time-frequency domain to achieve the purpose of stable and high-precision estimation of the formation equivalent Q-value for 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, reduction of Gibbs phenomenon, reduction of matrix ill-conditioning, etc. 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, etc. on the Q-value estimation result, thereby improving the Q-value estimation accuracy and stability.

[0035] 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 example specifically includes the following steps: Step 1. Calculate the time-varying logarithmic amplitude spectrum of the seismic data based on the short-time Fourier transform.

[0036] In Step 1, first perform the short-time Fourier transform on the seismic data to obtain the time-frequency spectrum of the seismic data, then perform modulus processing on the time-frequency spectrum of the seismic data to obtain the time-varying amplitude spectrum of the seismic data, and finally perform logarithmic processing on 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: Perform the short-time Fourier transform on the seismic data to obtain the time-frequency spectrum of the seismic data denoted as , then: (1) 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.

[0037] Perform modulus processing on the time-frequency spectrum of seismic data to obtain the time-varying amplitude spectrum of seismic data, denoted as , that is: (2) where represents the modulus processing of a complex number.

[0038] Perform logarithmic processing on the time-varying amplitude spectrum of seismic data to obtain the time-varying logarithmic amplitude spectrum of seismic data, which is calculated using the following formula: (3) where is a non-negative small constant, The role of is to avoid the situation where taking the logarithm of zero values in the data results in negative infinity.

[0039] 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 adopts a Gaussian function, and its expression is: (4) where is the standard deviation, and the standard deviation determines the width of the window.

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

[0041] In step 2, the time variable and frequency variable of the time-varying logarithmic amplitude spectrum of 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: 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 into the form of orthogonal decomposition of Legendre polynomials. If the logarithmic amplitude spectrum of the time-varying wavelet is expressed in the form of two-dimensional Legendre polynomial orthogonal decomposition, then the discrete time variable, that is, the time translation parameter and the frequency variable, that is, the frequency need to be mapped to through translation and scaling processing.

[0042] Let the time-varying logarithmic amplitude spectrum of seismic data The discrete-time variable, i.e., the time-shift parameter has a minimum value of 0 ms, and the time variable has a maximum value of , then any time variable within the time range can be mapped to the two-dimensional square-integrable function space spanned by Legendre polynomials through formula (5): on: (5) where represents the value of the time variable after mapping from the time variable in the two-dimensional real space, i.e., the time-shift parameter to the two-dimensional square-integrable function space .

[0043] When estimating the logarithmic amplitude spectrum of the time-varying wavelet, let the time-varying logarithmic amplitude spectrum of the seismic data whose frequency variable is the frequency has a minimum value of 0 Hz, and the frequency has a maximum value of , then any frequency variable within the frequency range can be mapped to the two-dimensional square-integrable function space spanned by Legendre polynomials through formula (6): on: (6) where represents the value of the frequency variable after mapping from the frequency variable in the two-dimensional real space to the two-dimensional square-integrable function space .

[0044] 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.

[0045] In Step 3, represent the logarithmic amplitude spectrum of the time-varying wavelet in the form of two-dimensional Legendre polynomials, and fit the time-varying logarithmic amplitude spectrum of the seismic data in the least squares sense to obtain the coefficients of each term of the two-dimensional Legendre polynomials. Substitute the coefficients into the two-dimensional Legendre polynomials to obtain the fitted logarithmic amplitude spectrum of the 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: Step 3.1. Orthogonal decomposition of the Legendre polynomials of the logarithmic amplitude spectrum of the time-varying wavelet.

[0046] In the square-integrable function space In it, any square-integrable function can be expanded into the form of orthogonal decomposition of Legendre polynomials. Therefore, the logarithmic amplitude spectrum of the time-varying wavelet is written as the form of orthogonal decomposition of two-dimensional Legendre polynomials as shown in formula (7): (7) Wherein, represents the logarithmic amplitude spectrum of the time-varying wavelet, is the Legendre polynomial of order is the Legendre polynomial of order represents the Legendre polynomial coefficient, , , and are the orders of orthogonal decomposition of Legendre polynomials in the time and frequency directions respectively.

[0047] The mathematical expression of is:(8) Wherein, , is the floor operation.

[0048] The mathematical expression of is:(9) Wherein, .

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

[0050] On the premise of the preset Legendre polynomial orders and , the logarithmic amplitude spectrum of the time-varying wavelet shown in formula (7) is used to approximate the logarithmic amplitude spectrum of the seismic data. By least squares fitting, the best Legendre polynomial coefficients are found, and the coefficients of each Legendre polynomial are calculated.

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

[0052] Substitute the calculated values of the coefficients of each order of Legendre polynomials 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. Its calculation formula is: (10) Substitute equations (5), (6), (8), and (9) into equation (10) to calculate the time-varying wavelet amplitude spectrum.

[0053] In addition, in step 3 of the formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain of this embodiment, when calculating the Legendre polynomial coefficients The specific process is as follows:[[]]END]] Solve for the two-dimensional Legendre polynomial coefficients by least squares fitting to minimize the following sum of squared residuals: (11) where represents the error function; take the partial derivative of the error function with respect to each Legendre polynomial coefficient : (12) where , .

[0054] Simplify equation (12) to obtain: (13) For the convenience of writing: Denote as ; Denote as .

[0055] Then let: .

[0056] , .

[0057] Rewrite equation (13) into the matrix form shown in equation (14): (14) where is a matrix of calculated from the Legendre polynomials, is the coefficient matrix of the to be solved, is the vector.

[0058] Solve equation (14) to obtain the coefficient matrix of the Legendre polynomial as: (15).

[0059] The coefficients of each Legendre polynomial are obtained from the calculation formula (15). .

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

[0061] In Step 4, 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 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 slope of the linear equation is solved using the least squares method, and the equivalent Q value is further calculated as the final result for output.

[0062] Specifically, based on obtaining the time-varying wavelet amplitude spectrum of the 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 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: (16) where is the instantaneous wavelet amplitude spectrum estimated at the th trace data at the moment, , is the total number of traces.

[0063] 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: (17) where corresponds to the equivalent Q value describing the average absorption effect of the formation.

[0064] Rewrite formula (17) into a linear equation as shown in formula (18): (18) 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, Indicates 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: (19) After further simplification, the formula for solving the formation equivalent Q value as shown in formula (20) is obtained: (20) 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.

[0065] 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, reducing matrix ill-conditioning, etc. 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.

[0066] To verify the effectiveness of the method proposed in the present invention, the present invention also gives two specific experiments: 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.

[0067] Experiment 1 Figure 2 Uses the Ricker wavelet with a main 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. The Ricker wavelet with a main frequency of 45 Hz is convolved with this reflection coefficient sequence to obtain a stationary synthetic seismic record as shown in Figure 4 . Using the time-varying wavelet of Figure 2 and the reflection coefficient sequence of Figure 3 to perform time-varying convolution processing, a synthetic seismic record with Q attenuation as shown in Figure 5 is obtained.

[0068] Figures 3 to 6 is a comparison chart of the synthetic seismic record before and after Q attenuation and the T compensation result. FromFigure 5 and Figure 4 It can be seen from the comparison that after the attenuation of the formation Q value, the amplitude of the synthetic seismogram gradually decreases, and the seismic wave energy is lost. For Figure 5 the attenuated synthetic seismogram, T compensation processing is carried out, and the result is as Figure 6 shown. Theoretically, T compensation processing only compensates the amplitude of seismic waves and does not change the bandwidth of the data. From Figure 6 it can be seen that after T compensation, the amplitude of the attenuated synthetic seismogram is effectively restored. Further comparison of 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 Figure 6 at 700 to 900 ms in the figure. In other words, Q attenuation not only reduces the amplitude of seismic waves, but also changes the frequency band of seismic data and reduces the resolution of 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, stable and high-precision estimation of the formation Q value is very necessary.

[0069] Figure 7 shows 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 normalization 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 seriously. Figure 9 is Figure 5 the time-varying amplitude spectrum obtained by short-time Fourier transform of the synthetic record with Q attenuation of Figure 9 . And for the convenience of comparing the frequency band changes in the time direction, energy gain processing is carried out in the time direction. It can be seen from

[0070] 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 . The parameters selected for two-dimensional Legendre polynomial decomposition are respectively and From the figure, it can be seen 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 serious. The estimated time-varying wavelet amplitude spectrum is very stable and is not affected by factors such as local strong reflection and waveform coupling of seismic records. Figure 10 The estimated value of Figure 8 The theoretical value of , we can see that the two are basically consistent, and further compare the difference between the estimated value and the theoretical value of the time-varying wavelet amplitude spectrum, such as Figure 11 As shown in the figure, it can be seen that the difference between the two is not obvious.

[0071] Further analysis of the estimation accuracy of the time-varying wavelet amplitude spectrum, Figure 12 The comparison diagrams of the estimated and theoretical values ​​of the instantaneous wavelet amplitude spectrum at 200ms and 600ms are given respectively. Figure 12 It can be seen that the estimated values ​​of the main frequencies of the instantaneous wavelet amplitude spectrum at 200ms and 600ms are basically consistent with the theoretical values, and the bandwidth of the estimated value of the instantaneous wavelet amplitude spectrum is also basically consistent with the theoretical value. Compared with the bandwidth of the estimated value of the instantaneous wavelet amplitude spectrum at 200ms, the bandwidth of the estimated value of the instantaneous wavelet amplitude spectrum at 600ms is closer to the theoretical value. The error of the estimated value of the instantaneous wavelet amplitude spectrum at 200ms mainly appears at the high-frequency end, while the error of the estimated value of the instantaneous wavelet amplitude spectrum at 600ms mainly appears in the frequency band range of the weak amplitude region. In short, the accuracy of the estimated values ​​of the instantaneous wavelet amplitude spectrum at 200ms and 600ms is relatively high. Comprehensive comparison Figure 9 , Figure 10 , Figure 11 and Figure 12 It can be seen 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 coefficient and waveform coupling on the instantaneous wavelet amplitude spectrum estimation. The estimation result is relatively stable and the instantaneous wavelet estimation accuracy is relatively high.

[0072] Figure 13 Given the use Figure 10 Comparison of the formation equivalent Q value calculated by 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 the following Figure 9 The estimation result is very stable without the influence of the local strong reflection factor near 820ms.

[0073] Figures 14 to 17 This is a comparison chart of Q compensation processing using the estimated formation equivalent Q value, where Figure 14 For a smooth synthetic seismic record, Figure 15 The synthetic seismic record after Q attenuation under the same reflection sequence conditions is used to estimate the Q value of the data using the method of the present invention. 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 Figure Figure 16 and Figure 14 shows the difference between them. It can be seen from the figure that the waveform difference on the synthetic seismic record caused by 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.

[0074] Experiment 2 Figure 18 is a grayscale display map of the post-stack actual data of a certain exploration area. It can be seen from the figure 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. It can be seen from the figure 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.

[0075] The formation equivalent Q value is calculated by combining the instantaneous wavelet amplitude spectrum with the travel time of the source wavelet. 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 the 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 the corresponding formation equivalent Q value should be relatively small.

[0076] Figure 19 Figure Figure 18 is the time-varying amplitude spectrum obtained by short-time Fourier transform processing of the 200th post-stack actual seismic data in Figure 20 Figure Figure 19 and Figure 20 is the time-varying wavelet amplitude spectrum estimated by the method of the present invention for this seismic data. Comparing

[0077] Since the post-stack data is generally the result of pre-stack resolution enhancement 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 wavelet 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 bandwidth.

[0078] 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.

[0079] The formation equivalent Q value of the area surrounded by the black dotted 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 dotted line frame 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

[0080] Example 2 This Example 2 describes a formation Q value estimation system based on 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 based on two-dimensional Legendre polynomial decomposition in the above Example 1.

[0081] The formation Q value estimation system based on 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.

[0082] 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.

[0083] An amplitude spectrum estimation module, which is used to estimate the time-varying wavelet amplitude spectrum by performing two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data in a two-dimensional square-integrable function space.

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

[0085] 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 described in the corresponding steps of the method in the above-mentioned Embodiment 1, and will not be elaborated here.

[0086] Embodiment 3 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.

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

[0088] Embodiment 4 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.

[0089] Of course, the above description is only a preferred embodiment 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 formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain, characterized in that: The steps include: Step 1. Calculate the time-varying logarithmic amplitude spectrum of seismic data based on short-time Fourier transform; Step 2. Map the time-varying logarithmic amplitude spectrum of seismic data into a 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 the seismic data to estimate the time-varying wavelet amplitude spectrum; Step 4. Calculate the equivalent Q value of the formation using the estimated time-varying wavelet amplitude spectrum.

2. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain according to claim 1 is characterized in that: The step 1 is specifically as follows: Seismic data Perform short-time Fourier transform and obtain the time-frequency spectrum of seismic data as : （1） in, is the observation time of earthquake record, is the time shift parameter, is the frequency, is the imaginary unit, represents the kernel function of the short-time Fourier transform; The time-varying amplitude spectrum of the seismic data is obtained by modulo processing the time-varying amplitude spectrum of the seismic data, which is recorded as : （2） in, Indicates the modulo processing of complex numbers; Time-varying amplitude spectrum of seismic data Take logarithmic processing to obtain the time-varying logarithmic amplitude spectrum of seismic data : （3） in, is a non-negative constant.

3. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain according to claim 2 is characterized in that: In step 1, the kernel function of the short-time Fourier transform Using Gaussian function, its expression is: （4） in, 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 step 2 is specifically as follows: Let the time-varying logarithmic amplitude spectrum of the seismic data be The discrete time variable of the time shift parameter The minimum value of the time variable is 0ms. The maximum value of ,but Any time variable in range can be mapped to the two-dimensional square integrable function space spanned by Legendre polynomials through formula (5): superior: （5） in, Represented by the time variable in two-dimensional real space Mapping to the two-dimensional space of square integrable functions The time variable value after Let the time-varying logarithmic amplitude spectrum of the seismic data be The frequency variable is the frequency The minimum value of the frequency is 0Hz. The maximum value of ,but Any frequency variable in the range can be mapped to the two-dimensional square integrable function space spanned by Legendre polynomials through formula (6): superior: （6） in, Represented by the frequency variable in the two-dimensional real number space Mapping to the two-dimensional space of square integrable functions The frequency variable value after .

5. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain according to claim 4 is characterized in that: The step 3 is specifically as follows: The logarithmic amplitude spectrum of the time-varying wavelet is written as the orthogonal decomposition of the two-dimensional Legendre polynomials shown in formula (7): （7） in, represents the logarithmic amplitude spectrum of the time-varying wavelet, for Legendre polynomials of degree, for The Legendre polynomials of degree, represents the Legendre polynomial coefficients, , , and The orders of the orthogonal decomposition of Legendre polynomials in the time and frequency directions respectively; The mathematical expression is: （8） in, , This is a round down operation; The mathematical expression is: （9） in, ; In the preset Legendre polynomial order and Based on the premise of , the mathematical expression of the logarithmic amplitude spectrum of the time-varying wavelet shown in formula (7) is used to approximate the logarithmic amplitude spectrum of the seismic data ; Solving the coefficients of two-dimensional Legendre polynomials using least squares fitting , each Legendre polynomial coefficient calculated is Substituting into formula (7), we get the logarithmic amplitude spectrum of the time-varying wavelet: , and Pick Exponential processing to obtain the time-varying wavelet amplitude spectrum , and its calculation formula is: （10） Substituting formula (5), formula (6), formula (8) and formula (9) into formula (10), the time-varying wavelet amplitude spectrum is calculated.

6. The formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain according to claim 5, characterized in that: In step 3, the Legendre polynomial coefficients are calculated The specific process is: Solving for 2D Legendre Polynomial Coefficients by Least Squares Fitting , so that the residual sum of squares is minimized as shown in formula (11): （11） in, represents the error function; the error function For each Legendre polynomial coefficient Find the partial derivative: （12） in, , ; Simplifying formula (12) yields: （13） Will Recorded as ; Will Recorded as ; Then let: ； ， ； Rewrite formula (13) into the matrix form shown in formula (14): （14） in, is calculated using Legendre polynomials The matrix of It is waiting for The coefficient matrix of yes A vector of Solving formula (14) yields the coefficient matrix of Legendre polynomials: （15） Calculate formula (15) to get the coefficients of each Legendre polynomial .

7. 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 step 4 is specifically as follows: The instantaneous wavelet amplitude spectrum with the widest frequency band at the first arrival time of the reflection wave in the time-varying wavelet amplitude spectrum is taken as the source wavelet amplitude spectrum; Perform multi-channel statistical averaging on the multi-channel data and record the first arrival time of the reflection wave as , the road number is To express it, the amplitude spectrum of the source wavelet after multi-channel statistics is It is expressed as: （16） in, is the frequency, For the Estimated from the data The instantaneous wavelet amplitude spectrum at time , , is the total number of channels; by As the initial time for estimating the equivalent Q value of the formation, using the Q value definition formula, the wave propagation time is The amplitude spectrum of the source wavelet after It is expressed as: （17） in, is the equivalent Q value describing the average absorption effect of the formation; Rewrite formula (17) into the linear equation shown in formula (18): （18） The frequency band range selected for calculating the equivalent Q value of the formation is recorded as ,in is the total number of sample points of the selected frequency, represents the starting frequency for equivalent Q value estimation, represents the cutoff frequency used for equivalent Q value estimation, and the slope of the linear equation is solved using the least squares method. ,Right now: （19） Then we get the formula for solving the equivalent Q value of the formation as shown in formula (20): （20） The calculated time-varying wavelet amplitude spectrum and source wavelet amplitude spectrum Substitute it into formula (20) to calculate the equivalent Q value of the formation that changes with time, and output the calculated equivalent Q value of the formation as the final estimation result.

8. A formation Q value estimation system based on two-dimensional Legendre polynomial decomposition in time-frequency domain, characterized in that: Includes the following modules: A time-varying logarithmic amplitude spectrum calculation module, used for calculating the time-varying logarithmic amplitude spectrum of seismic data based on short-time Fourier transform; A mapping processing module, used for mapping the time-varying logarithmic amplitude spectrum of seismic data into a two-dimensional square integrable function space spanned by Legendre polynomials; An amplitude spectrum estimation module is used to estimate the time-varying wavelet amplitude spectrum by performing two-dimensional Legendre polynomial decomposition on the time-varying logarithmic amplitude spectrum of seismic data in a two-dimensional square integrable function space; and a formation equivalent Q value estimation module, which is used to calculate the formation equivalent Q value by using the estimated time-varying wavelet amplitude spectrum.

9. A computer device comprising a memory and one or more processors, wherein the memory stores executable code, characterized in that: When the processor executes the executable code, the steps of the formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in time-frequency domain as described in any one of claims 1 to 7 are implemented.

10. A computer-readable storage medium having a program stored thereon, characterized in that: When the program is executed by a processor, the steps of the formation Q value estimation method based on two-dimensional Legendre polynomial decomposition in the time-frequency domain as described in any one of claims 1 to 7 are implemented.

Images