A high-resolution seismic processing method and system, storage medium, and electronic device.

Abstract

This invention discloses a high-resolution seismic processing method and system, storage medium, and electronic device, applied to the field of seismic data processing. Addressing the problem that existing techniques suffer from poor spectral extension effects and difficulty in identifying thin-layer earthquakes, this invention, based on time-frequency analysis, performs deconvolution, harmonic decomposition, and basis pursuit inversion on the time-frequency spectrum of seismic data. This achieves high-frequency reconstruction of the time-frequency spectrum, improving the resolution of seismic data and forming a novel spectral extension method. Compared to previous spectral recovery methods, this method achieves better spectral extension effects on more complex strata, with higher resolution and signal-to-noise ratio.

Description

Technical Field

[0001] This invention belongs to the field of earthquake data processing, and specifically relates to a thin-layer earthquake identification technology. Background Technology

[0002] With the deepening of oil and gas seismic exploration, thin-layer seismic identification has received widespread attention from the oil and gas industry. A thin layer refers to a stratum where reflected seismic waves interfere on its top and bottom surfaces. Due to the limitations of seismic data resolution, it is usually difficult to identify on seismic profiles. Widess et al. were among the first to study thin layers, defining the limiting resolution of a thin layer as less than one-eighth of the seismic wavelet wavelength. Considering the influence of noise and wavelet factors, in practical applications, geological and geophysical experts define the limiting resolution of a thin layer as one-quarter of the wavelet wavelength. In this sense, a thin layer refers to a stratum with a thickness less than one-quarter of the seismic wavelet wavelength and difficult to distinguish on a seismic profile. Therefore, improving the resolution of seismic data has become an indispensable issue in the study of thin layers. Research shows that the thinner the target layer, the larger the notch period. To distinguish this thin layer on a seismic profile, a higher seismic data frequency and a wider seismic data bandwidth are required, i.e., an effective expansion of the data frequency bandwidth is needed. As the dominant frequency and bandwidth of seismic data increase, the resolution of seismic data is reasonably improved. Therefore, broadband seismic data plays a very important role in thin-layer identification.

[0003] Due to the filtering effect of the earth, especially the strong absorption effect of the surface layer, seismic data acquired through seismic exploration is characterized by low frequency and narrow bandwidth, which seriously affects the accuracy of thin-layer identification. To identify thin layers, it is necessary to increase the dominant frequency of the seismic data while expanding its bandwidth; that is, to reasonably compensate for the high-frequency information while preserving the low-frequency information. Currently, there are many methods for bandwidth expansion, mainly divided into four categories: spectral whitening, deconvolution, absorption compensation, and spectral restoration techniques.

[0004] Spectral whitening methods do not alter the phase spectrum; instead, they whiten the amplitude spectrum of seismic data within a specific frequency band, compensating for overall frequency distribution. Chen Chuanren et al. combined wavelet transform with spectral whitening to analyze local frequency characteristics of the signal. Wang Ji proposed a spectral whitening method based on Hilbert-Huang transform, which can simultaneously enhance local signal characteristics in both the time and frequency domains. While spectral whitening methods can enhance the high-frequency components of the data, they neglect the influence of phase on seismic wavelet resolution, thus disrupting the spatial relationships between amplitudes at different frequencies.

[0005] Deconvolution methods, based on seismic convolution models, improve the temporal resolution of seismic data by compressing seismic wavelets to recover the reflection coefficient sequence. These methods include impulse deconvolution, homomorphic deconvolution, predictive deconvolution, and spectral simulation deconvolution. Robinson proposed the assumptions of the convolution model and the predictive deconvolution algorithm, laying the theoretical and applied foundation for seismic deconvolution. Oppenheim proposed homomorphic deconvolution, which separates the seismic wavelet and reflection coefficient in the complex spectrum using nonlinear filtering. Zhao et al. proposed the spectral simulation deconvolution method, which assumes that the wavelet amplitude spectrum is much smoother than the reflection coefficient amplitude spectrum. It uses polynomial fitting to the amplitude spectrum of the seismic record to obtain the wavelet amplitude spectrum and broaden its bandwidth, thereby improving the resolution of seismic data. However, these methods based on convolution theory do not fundamentally extract effective information, making it difficult to guarantee the fidelity and signal-to-noise ratio after processing, and they may not reflect the true state of the strata.

[0006] Commonly used absorption compensation techniques include inverse Q-filtering and spectral whitening. Inverse Q-filtering, by estimating the formation quality factor Q, can compensate for amplitude loss and frequency and phase attenuation of seismic waves, thereby improving the resolution of seismic data. Hale first proposed the inverse Q-filtering algorithm, and Wang proposed a stable and effective algorithm that can simultaneously compensate for amplitude attenuation and phase distortion, extending the algorithm to cases where the Q value varies with time or depth. Wang Shoudong reduced the inverse Q-filtering problem to an inversion problem, using regularization methods to achieve attenuation compensation through inversion. Although the inverse Q-filtering method can compensate for high-frequency components to some extent, the quality factor Q is often difficult to obtain accurately, leading to inaccurate inverse Q-filtering results.

[0007] Methods for improving seismic data resolution, such as spectral whitening, deconvolution, and absorption compensation, are limited to the cutoff bandwidth of the seismic data. While various mathematical algorithms enhance high-frequency information, the resulting high-resolution results are highly arbitrary and fail to reflect the true geological conditions. Spectral restoration methods predict and reconstruct high-frequency components from low-frequency components of seismic data, rather than simply enhancing them. For example, Castagna's spectral restoration technique for sparse layer reflection coefficient inversion recovers and reconstructs high-frequency information from the original seismic data, achieving a resolution of one-eighth or even one-sixteenth of the wavelength, significantly improving resolution and fidelity compared to previous methods. Harmonic extrapolation decomposes the spectrum of the reflection coefficient sequence within the original seismic data's frequency band and extends the spectrum based on the periodicity of sine and cosine functions. Although these spectral restoration techniques can broaden the frequency band of seismic data, they still have the following problems:

[0008] (1) Most spectrum recovery methods, including harmonic extrapolation, analyze the spectrum in the frequency domain. However, Fourier transform obtains the spectrum characteristics of the entire time domain, which cannot characterize the time corresponding to the frequency change characteristics of the signal.

[0009] (2) Faced with seismic responses formed by many thin-layer tuning, the periodicity of the reflection coefficient spectrum is difficult to fully characterize in the seismic data frequency band. More importantly, for complex strata, the reflection coefficient sequence may be a closely arranged random sequence, and the corresponding reflection coefficient spectrum is superimposed with many periodic functions, resulting in a very complex structure, which increases the uncertainty of spectrum recovery of strata reflection coefficient inversion.

[0010] Identifying thin layers requires broadening the frequency band of seismic data and improving its resolution. In recent years, various high-resolution processing methods have emerged, including deconvolution algorithms, absorption compensation techniques, and spectral reconstruction techniques. However, most spectral reconstruction methods are limited to the high signal-to-noise ratio frequency bandwidth of seismic data, failing to fundamentally add more effective information. Furthermore, inappropriate parameter selection makes it difficult to guarantee the fidelity and signal-to-noise ratio of the processed seismic data. While spectral reconstruction techniques, such as harmonic extrapolation, can calculate effective high-frequency components from low-frequency components of the data, offering better resolution and fidelity, the limitations of Fourier transform often result in poor spectral extrapolation performance in complex geological formations. Summary of the Invention

[0011] To address the aforementioned technical problems, this invention proposes a high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients. This method performs harmonic decomposition and spectral extension on the simpler time-frequency spectrum of the reflection coefficient sequence, restoring high-frequency information while preserving low-frequency information, thus significantly improving the accuracy of spectral extension.

[0012] The technical solution adopted in this invention is: a high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients, comprising:

[0013] S1. Input a seismic data point x(t) and perform a generalized S-transform; obtain the time spectrum GST(τ,f):

[0014]

[0015] Where f represents the frequency from 0 to the Nyquist frequency, t is the time position, τ represents the center time of the Gaussian window, and λ and p are parameters;

[0016] S2. In the time-frequency domain, take the time spectrum GST(τ,f) at each time point t. i The corresponding amplitude spectrum is denoted as GST(t) i ,f); GST(t) amplitude spectrum iDividing f by the seismic wavelet spectrum W(f) yields the instantaneous spectrum R of the reflection coefficient sequence within the amplitude spectral bandwidth at each time point. i (f);

[0017] S3, the instantaneous spectrum R of the reflection coefficient sequence i (f) Obtain the weighted sum of sine and cosine functions by performing harmonic decomposition:

[0018]

[0019] S4. Use the basis pursuit inversion method to broaden the instantaneous spectrum of the reflection coefficient sequence;

[0020] S5. Multiply the instantaneous spectrum of the reflection coefficient sequence after frequency extension with the sub-wavelength spectrum whose main frequency is half the bandwidth of the extended frequency, and extend the amplitude spectrum at each time point; to obtain the time spectrum GST_high(τ,f) after high frequency information enhancement.

[0021] S6. Perform a generalized inverse S-transform on the enhanced time spectrum GST_high(τ,f) to obtain high-resolution seismic data in the time domain.

[0022] W(f) is specifically estimated from earthquake data, or a zero-phase Reich wavelet with the same dominant frequency as the earthquake signal is selected.

[0023] The beneficial effects of this invention: Based on the generalized S-transform and harmonic decomposition theory, this invention proposes a high-resolution seismic processing technique for instantaneous spectral extension of the reflection coefficient in the time-frequency domain, which can improve the resolution of seismic data. Experimental results using synthetic and real data demonstrate that harmonic decomposition and basis tracing inversion of the reflection coefficient spectrum based on the time-frequency spectrum can yield seismic data with higher dominant frequencies, effectively broadening the bandwidth. While maintaining the components within the effective bandwidth of the original seismic data, the bandwidth is extended, significantly improving the resolution of seismic data and providing favorable seismic data conditions for thin-layer identification and thin-layer oil and gas exploration. Attached Figure Description

[0024] Figure 1 A two-layer reflection coefficient model;

[0025] Figure 2 This is a flowchart of the method of the present invention;

[0026] Figure 3 A comparison of the resolvable effects of synthesized data before and after generalized S-transform spectral extension;

[0027] Among them, (a) is the original data with a dominant frequency of 30Hz; (b) the dominant frequency is increased to 60Hz by the spectral extension method based on the generalized S-transform; (c) is the synthetic seismic data obtained by convolving the reflection coefficient model with the 30Hz Ricker wavelet; and (d) is the result after the dominant frequency of the synthetic data is increased to 60Hz by the spectral extension method based on the generalized S-transform.

[0028] Figure 4 This is a comparison chart of the data in track 5 before and after generalized S-transform spectral extension;

[0029] (a) is the original data profile; (b) is the profile after generalized S-transform spectral extension; (c) is a comparison of the amplitude spectra of the data before and after generalized S-transform spectral extension of the 5th channel.

[0030] Figure 5 A comparison of the generalized S-transform spectrum before and after the spectral extension of the fifth data channel;

[0031] Wherein, (a) is the time spectrum of the generalized S-transform of the original data of the 5th channel; (b) is the time spectrum after the main frequency is extended to 60Hz;

[0032] Figure 6 The figure shows a comparison between the generalized S-transform spectral extension and the harmonic extrapolation method based on actual data.

[0033] Among them, (a) is the original data profile; (b) is the profile after spectral extension based on generalized S-transform; and (c) is a comparison of the amplitude spectra of the actual data before and after spectral extension by generalized S-transform.

[0034] Figure 7 This is a high-resolution seismic processing system for instantaneous spectral extension of time-frequency domain reflection coefficients. Detailed Implementation

[0035] To facilitate understanding of the technical content of this invention by those skilled in the art, the following techniques will be described first:

[0036] 1. Reflection Coefficient Model

[0037] For example Figure 1 The two-layer stratigraphic model shown has the reflection coefficients expressed in the time domain as follows:

[0038] g(t)=r1δ(t-t1)+r2δ(t-t1-T) (1)

[0039] Where: r1 is the top reflection coefficient; r2 is the bottom reflection coefficient; δ(t) is the unit impulse function; t is the time position, t1 is the time position of the top reflection coefficient, and T is the time thickness of the formation. If the analysis point is located at the midpoint of the formation, the reflection coefficient at that point can be expressed as:

[0040]

[0041] When the reflection coefficients e1 and r2 at the top and bottom interfaces have equal amplitudes and the same sign, this layer is called an even pulse pair. When the amplitudes r1 and r2 are equal but opposite in sign, this layer is called an odd pulse pair.

[0042] Performing a Fourier transform on equation (2) yields the spectra of even pulse pairs and odd pulse pairs as follows:

[0043] II(f)=2rcos(πΔtf) (3)

[0044] I I (f)=i2rsin(πΔtf) (4)

[0045] Where: i is the imaginary unit; r is the reflection coefficient; Δt is the temporal thickness of the formation; f includes all frequencies from 0 to the Nyquist frequency. Any reflection coefficient sequence can be considered as a superposition of even pulse pairs and odd pulse pairs at the analysis point.

[0046] Generally, the thinner the stratum, the longer the notch period. A thin stratum can only be resolved when the seismic data band includes at least one notch period. Therefore, the thinner the stratum, the longer the notch period, and the wider the seismic data band required to resolve the stratum.

[0047] 2. Harmonic extrapolation method

[0048] In the convolution model, seismic records in the time domain can be represented as:

[0049] s(t)=w(t)*r(t) (5)

[0050] Where s(t) is the earthquake record, w(t) is the seismic wavelet, and r(t) is the reflection coefficient sequence.

[0051] Seismic records are represented in the frequency domain as follows:

[0052] S(f)=W(f)×R(f) (6)

[0053] Where: S(f) is the seismic record spectrum, W(f) is the seismic wavelet spectrum, and R(f) is the reflection coefficient sequence spectrum.

[0054] The harmonic extrapolation method first requires deconvolution of the spectrum within the original seismic data bandwidth, i.e., dividing the seismic data spectrum by the seismic wavelet spectrum to obtain the reflection coefficient sequence spectrum within the seismic data bandwidth. The real part of the deconvolution result is a superposition of cosine functions of different periods, and the imaginary part is a superposition of sine functions of different periods. Then, based on the periodicity of the sine and cosine functions, harmonic decomposition is used to decompose the reflection coefficient sequence spectrum within the bandwidth into a weighted sum of several sine and cosine functions. A reflection coefficient sequence with 2N+1 sampling points corresponds to K = N+1 pulse pairs. The spectrum of any reflection coefficient pair can be represented by the spectra of even pulse pairs and odd pulse pairs as follows:

[0055] G(f,n)=2r e ·cos(2π·n·dt·f)+i2r o (n)·sin(2π·n·dt·f)) (7)

[0056] Where: dt is the sampling rate, n is half the time thickness of the reflection coefficient pair, and r e and r o These are the coefficients of the spectrum of even pulse pairs and the spectrum of odd pulse pairs, respectively. Therefore, the spectrum of any reflection coefficient sequence from 0 to the Nyquist frequency within the analysis window can be expressed as:

[0057]

[0058] Convert formula (8) into matrix form:

[0059] S d =φa+ε (9)

[0060] Wherein: S d The spectrum of the reflection coefficient sequence obtained by deconvolution is represented by φ, which is a kernel matrix composed of sine elements, and a is a matrix containing r e and r o The coefficient vector is ε, where ε is the prediction residual. The kernel matrix has a dimension of M×N, where M is the number of frequencies determined by the available bandwidth and frequency sampling rate, and K is the number of pulse pairs.

[0061] To improve the stability of the inversion, this invention uses the basis pursuit inversion method to solve equation (9), decomposing the data spectrum within the original bandwidth into a sparse superposition of the formation response, and obtaining the output solution r. e (n) and r o (n), corresponding to the real and imaginary parts of the reflection coefficient spectrum, respectively. Based on the periodicity of the sine and cosine functions, the reconstructed reflection coefficient spectrum is broadened beyond the bandwidth of the original data. The summation formula is then used to sum the values ​​from K=1 to N+1 (where N is the Nyquist frequency) to include all possible frequency periods, thus obtaining the frequencies of this reflection coefficient sequence from 0 to the Nyquist frequency.

[0062] The basis pursuit inversion method solves for the coefficients in equation (5) by simultaneously minimizing the l2 norm of the error term and the l1 norm of the solution:

[0063]

[0064] Where λ is the regularization parameter used to impose sparsity constraints. Increasing the value of λ will produce results with higher sparsity.

[0065] By multiplying the spectrum of the reflection coefficient sequence after spectral broadening by the sub-wavelength spectrum with a higher dominant frequency and wider bandwidth, and performing an inverse Fourier transform, we can obtain seismic data with a wider bandwidth and higher resolution.

[0066] 3. Time-Frequency Analysis – Generalized S-Transform

[0067] The S-transform improves upon the short-time Fourier transform and wavelet transform by using a Gaussian window function whose width varies inversely with frequency, while preserving the signal's phase information. This provides adaptive time-frequency resolution. The S-transform of a non-stationary signal h(t) is defined as follows:

[0068]

[0069] Where f represents all frequencies from 0 to the Nyquist frequency, t is the time variable, τ represents the center time of the Gaussian window, and ω(t) represents the Gaussian window function:

[0070]

[0071] Standard deviation is defined as the reciprocal of frequency:

[0072]

[0073] |f| represents the absolute value of f;

[0074] To overcome the limitation of the Gaussian window function in the S-transform having a fixed variation trend, the window function of the S-transform is improved. Since the signal changes drastically in the high-frequency range with a relatively short time period, the time window should be narrower, while in the low-frequency range, the signal changes relatively smoothly with a relatively long time period, requiring a wider time window. Therefore, to allow for more flexible adjustment of the window function, parameters μ and p are introduced into the Gaussian window function, enabling the window function ω(t) to vary with the frequency scale f, resulting in the generalized S-transform. The formula for the standard deviation is as follows:

[0075]

[0076] The expression for the generalized S-transform is as follows:

[0077]

[0078] The window function of the generalized S-transform exhibits nonlinear variation with frequency, offering greater flexibility and higher resolution, and superior time-frequency focusing compared to the S-transform.

[0079] Example 1

[0080] Compared to the traditional Fourier transform, time-frequency analysis, as a non-stationary signal processing method, can describe the instantaneous changes and patterns of the spectral energy of a signal over time. The calculated reflection coefficient spectrum exhibits better sparsity and regularity in the time-frequency domain. Time-frequency analysis methods mainly include short-time Fourier transform, continuous wavelet transform, S-transform, and generalized S-transform. Among them, the generalized S-transform can flexibly adjust the width and attenuation trend of the window function, resulting in better time-frequency focusing. Combining the generalized S-transform with harmonic decomposition and performing seismic spectral extension in the time-frequency domain of seismic data can achieve higher resolution. This invention, based on time-frequency analysis and harmonic extrapolation spectral recovery techniques, proposes a high-resolution seismic processing technique for instantaneous spectral extension of reflection coefficients in the time-frequency domain. This invention transforms the object of spectrum recovery from the complex Fourier transform spectrum to the simpler instantaneous spectrum. Harmonic decomposition and spectral extension are performed on the amplitude spectrum corresponding to each moment of the instantaneous spectrum, thereby enabling effective prediction and compensation of high-frequency information in seismic data while maintaining low-frequency information, thus improving the resolution of seismic data and laying a data foundation for thin-layer identification.

[0081] This invention first performs a generalized S-transform on the seismic signal to obtain a time-frequency focused spectrum. Then, time-varying wavelet deconvolution and harmonic decomposition are performed on the amplitude spectrum corresponding to each moment in the time-frequency spectrum. This decomposes the instantaneous spectrum of the reflection coefficient sequence within the high signal-to-noise ratio band into a weighted sum of sine and cosine functions. Basis pursuit inversion is then used to extend the instantaneous spectrum of the reflection coefficient sequence beyond a finite frequency band. Multiplying this spectrum with the wavelet spectrum of a higher dominant frequency and performing a generalized inverse S-transform yields a wideband, high-resolution seismic signal, enabling seismic interpreters to better identify the spatial distribution of thin layers. Figure 2 As shown, the process includes the following:

[0082] Step 1: Input one seismic data point and perform a generalized S-transform.

[0083] First, perform a generalized S-transform on the seismic signal x(t) according to formula (16):

[0084]

[0085] Where f represents all frequencies from 0 to the Nyquist frequency, and τ represents the center time of the Gaussian window. Adjusting the parameters p and μ can adjust the width and attenuation trend of the window function. Choosing appropriate values ​​of p and μ optimizes the focusing and resolution of the generalized S-transform time spectrum, resulting in the time spectrum GST(τ,f).

[0086] Step 2: For t i Deconvolution of the amplitude spectrum at time intervals yields the reflection coefficient sequence spectrum within the frequency band.

[0087] In the time-frequency domain, take the generalized S-transform time spectrum GST(τ,f) at each time point t. i The corresponding amplitude spectrum is denoted as GST(t) i ,f). The amplitude spectrum GST(t) i Dividing W(f) by the seismic wavelet spectrum W(f) (which can be estimated from seismic data or by selecting a zero-phase Ricker wavelet with the same dominant frequency as the seismic signal) yields the instantaneous spectrum R of the reflection coefficient sequence within the amplitude spectral bandwidth at each time point. i (f).

[0088] Step 3: Perform harmonic decomposition on the instantaneous spectrum of the reflection coefficient sequence to obtain the weighted sum of sine and cosine functions.

[0089] Since the spectrum of any pair of reflection coefficients can be represented by the spectra of even pulse pairs and odd pulse pairs, the spectrum of the reflection coefficient sequence within the bandwidth can be decomposed into a weighted sum of sine and cosine functions through harmonic decomposition:

[0090]

[0091] Step 4: Broaden the reflection coefficient spectrum using the basis pursuit inversion method

[0092] Convert formula (17) into matrix form:

[0093] S d =φa+ε (18)

[0094] Wherein: S d The spectrum of the reflection coefficient sequence obtained by deconvolution is represented by φ, which is a kernel matrix composed of sine elements, and a is a matrix containing r e and r o The reflection coefficient vector, ε is the prediction residual. The basis pursuit method uses the L1 norm instead of the L0 norm to solve the optimization problem, and describes equation (18) as a sparse solution problem under constraints:

[0095]

[0096] Where: the first term of the constraint condition ||S d-Φa||2 represents the L2 norm of the vector obtained by subtracting the spectrum of the reflection coefficient sequence obtained by deconvolution from the spectrum of the reconstructed reflection coefficients, and the second term ||a||1 represents the L1 norm of vector a. The value that satisfies the constraint condition of equation (19) is the reflection coefficient vector a obtained by basis pursuit inversion. λ is a regularization parameter, which is generally between 0 and 1. Increasing λ within this range will produce a more sparsity structure, while decreasing λ may amplify the inversion noise. It is usually taken as 0.01.

[0097] The instantaneous spectrum of the reflection coefficient sequence is decomposed into a sparse superposition of formation responses using the basis pursuit inversion method, and the weighting coefficients r of the sine and cosine functions are output. e (n) and r o (n), corresponding to the real and imaginary parts of the instantaneous spectrum, respectively. The instantaneous spectrum of the reflection coefficient sequence from 0 to the Nyquist frequency is calculated by equation (17), thus broadening the instantaneous spectrum of the reflection coefficient sequence within the bandwidth to outside the original bandwidth.

[0098] Step 5: Multiply the instantaneous spectrum of the frequency-spread reflection coefficient sequence by the sub-wavelength spectrum of a certain dominant frequency.

[0099] Multiplying the spectrum of the extended reflection coefficient sequence with the spectrum of the sub-wavelength whose main frequency is half the bandwidth of the extended frequency allows us to obtain the t i The amplitude spectrum at time t is frequency-extended.

[0100] Repeat steps 2 to 5 to achieve frequency extension of the instantaneous amplitude spectrum at each moment, which means achieving frequency extension of the entire time spectrum and obtaining the time spectrum GST_high(τ,f) after high-frequency information enhancement.

[0101] GST high(τ,f) ={R1(f)·W high(f) , ..., R i (f)·W_high(f), …R N (f)·W_high(f)} (20)

[0102] Among them, R i (f) is the t calculated in Step 4. i The instantaneous spectrum of the reflection coefficient sequence from time 0 to the Nyquist frequency, i = 1, ..., N, where N is the number of sampling points; W_high(f) is the sub-wave spectrum with the main frequency being half the bandwidth after the frequency extension.

[0103] Step 6: Perform a generalized inverse S-transform on the frequency-extensioned time-spectrum.

[0104] Since the generalized S-transform is lossless and reversible, performing the generalized inverse S-transform on the frequency-spreaded time spectrum yields high-resolution seismic data in the time domain:

[0105]

[0106] f(t) is the broadband seismic signal after spectral extension.

[0107] Steps 1 through 6 are repeated for each seismic data track in the study area to achieve spectral extension of the entire seismic data volume, ultimately obtaining a high-resolution seismic data volume. This invention is suitable not only for improving the resolution of post-stack seismic data but also for improving the resolution of pre-stack seismic data.

[0108] The technical effects of the present invention are illustrated by the following experiments:

[0109] Synthetic data validation

[0110] A wedge-shaped model was established to verify the effectiveness of the method of the present invention. Figure 3 (a) and Figure 3 (b) shows the velocity model and reflection coefficient model used when synthesizing seismic data. Figure 3 (c) is the synthetic seismic data obtained by convolving the reflection coefficient model with the 30Hz Ricker wavelet. The target layer thickness on the left side of the model is 0ms, the thickness on the rightmost side is 30ms, and the sampling interval is 1ms. Figure 3 (d) shows the result after the frequency of the synthesized data is increased to 60Hz by the spectral extension method based on the generalized S-transform. Figure 3 (c) and Figure 3 (d) Comparative analysis shows that only the target strata with a time thickness greater than 7ms can be distinguished in the original synthetic data, while the seismic high-resolution data with instantaneous spectrum extension of time-frequency domain reflection coefficient can distinguish strata with a time thickness of 3ms. Figure 4 (a) and Figure 4 (b) The profiles of the original data and the data after spectral extension are shown respectively. By comparing them with the velocity model and the reflection coefficient model, it can be found that the high-resolution profile of the instantaneous spectral extension of the time-frequency domain reflection coefficient not only improves the main frequency, but also has good amplitude preservation. Figure 4 (c) By comparing the Fourier amplitude spectra of the fifth synthetic data before and after frequency extension, it can be observed that the effective bandwidth of the synthetic data amplitude spectrum after the instantaneous spectrum extension of the time-frequency domain reflection coefficient is significantly increased, and the components within the effective bandwidth of the original seismic data are maintained while enhancing the high-frequency components.

[0111] Figure 5 The generalized S-transform time spectrum before and after the spectral extension of the fifth data channel was compared. Figure 5 (a) is the spectrum of the generalized S-transform of the original data. Figure 5 (b) is the generalized S-transform time spectrum after extending the main frequency to 60Hz. It can be seen that the high-frequency information of the extended time spectrum is significantly enhanced, while the low-frequency information is not damaged, proving the effectiveness of the invention.

[0112] Actual data verification

[0113] The method of the present invention was verified with actual data. Figure 6 The results based on the generalized S-transform spectral extension and the harmonic extrapolation method were compared. Figure 6 (a) is the original cross-section. Figure 6 (b) is a profile based on the spectral extension of the generalized S-transform. Figure 6 (c) This section shows a comparison of the spectra before and after the generalized S-transform spectral extension of the 300th actual data. After spectral extension using the generalized S-transform, the dominant frequency of the seismic data is increased, the bandwidth is widened, and the overall resolution of the profile is improved. This makes previously difficult-to-identify thin layers clearer, and the resulting profile has better amplitude preservation and signal-to-noise ratio. The seismic data bandwidth is widened towards higher frequencies, and the high-frequency components are significantly enhanced, while the components within the effective bandwidth of the original seismic data remain undamaged.

[0114] Experimental results using synthetic and real data demonstrate that the seismic high-resolution processing technique based on the instantaneous spectral extension of the time-frequency domain reflection coefficient proposed in this invention improves the resolution of seismic data, broadens the seismic data bandwidth, and effectively preserves the components within the effective bandwidth of the original seismic data while enhancing high-frequency components. Compared with previous methods and techniques, it has a higher signal-to-noise ratio and amplitude preservation.

[0115] Example 2

[0116] This invention also provides a high-resolution seismic processing system for instantaneous spectral extension of time-frequency domain reflection coefficients, such as... Figure 7 As shown, it includes: a generalized S-transform module, a reflection coefficient sequence instantaneous spectrum extraction module, a harmonic decomposition module, a broadening module, a high-frequency information enhancement module, and a generalized inverse S-transform module; the generalized S-transform module takes a seismic data point x(t) as input and outputs the time spectrum GST(τ,f) corresponding to x(t); the reflection coefficient sequence instantaneous spectrum extraction module takes the time spectrum GST(τ,f) as input and outputs the reflection coefficient sequence instantaneous spectrum R0. i (f); The input to the harmonic decomposition module is the instantaneous spectrum R of the reflection coefficient sequence. i (f) The output is the instantaneous spectrum R of the reflection coefficient sequence. i (f) is expressed as a weighted sum of sine and cosine functions; the input to the broadened module is the instantaneous spectrum R of the reflection coefficient sequence. i The weighted sum of sine and cosine functions (f) is used as the input to the frequency-extended reflection coefficient sequence spectrum. The high-frequency information enhancement module takes the frequency-extended reflection coefficient sequence spectrum as its input and outputs the frequency-extended time spectrum GST_high(τ,f). The generalized S-inverse transform module takes the frequency-extended time spectrum GST_high(τ,f) as its input and outputs the broadband seismic signal with extended spectrum.

[0117] Example 3

[0118] This embodiment also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, performs the steps of the above-described high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients.

[0119] Example 4

[0120] This embodiment also provides an electronic device, including: a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the electronic device is running, the processor communicates with the memory via the bus. When the machine-readable instructions are executed by the processor, the steps of the above-described high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients are performed.

[0121] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the invention should be included within the scope of the claims of the invention.

Claims

1. A high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients, characterized in that, include: S1. Input a seismic data point x(t) and perform a generalized S-transform; The time spectrum GST(τ,f) is obtained: Where f represents the frequency from 0 to the Nyquist frequency, t is the time position, τ represents the center time of the Gaussian window, and μ and p are parameters; S2, in the time-frequency domain, take the time-frequency spectrum GST(τ, f) of each time point t i The corresponding amplitude spectrum is denoted as GST(t i ,f); divide the amplitude spectrum GST(t i ,f) by the seismic wavelet spectrum W(f) to obtain the reflection coefficient sequence instantaneous spectrum R i (f) within the amplitude spectrum bandwidth of each time point. S3, the reflection coefficient sequence instantaneous spectrum R i (f) performing a harmonic decomposition to obtain a weighted sum of cosine and sine functions: Where dt is the sampling rate, n is half the time thickness of the reflection coefficient pair, and r e (n) and r o (n) are the coefficients of the spectrum of even pulse pairs and the spectrum of odd pulse pairs, respectively; S4. Use the basis pursuit inversion method to broaden the instantaneous spectrum of the reflection coefficient sequence; S5. Multiply the instantaneous spectrum of the reflection coefficient sequence after frequency extension with the sub-wavelength spectrum whose main frequency is half the bandwidth of the extended frequency, and extend the amplitude spectrum at each time point; to obtain the time spectrum GST_high(τ,f) after high frequency information enhancement. S6. Perform a generalized inverse S-transform on the enhanced time spectrum GST_high(τ,f) to obtain high-resolution seismic data in the time domain.

2. The high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients according to claim 1, characterized in that, In step S2, W(f) is specifically estimated from the seismic data, or a zero-phase Reich wavelet with the same dominant frequency as the seismic signal is selected.

3. The high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients according to claim 1, characterized in that, Step S4 specifically includes the following sub-steps: S41, the instantaneous spectrum R of the reflection coefficient sequence i The weighted sum of sine and cosine functions in (f) can be converted into matrix form: S d =φa+e Among them, S d φ represents the spectrum of the reflection coefficient sequence obtained by deconvolution, φ is the kernel matrix composed of sinusoidal elements, a is the coefficient vector, and ε is the prediction residual; S42. Solve the optimization problem by replacing the L0 norm with the L1 norm. Specifically, describe step S41 in matrix form as a sparse solution problem under constraints: Where λ is the regularization parameter, and the first term of the constraint is ||S d -Φa||2 represents the spectrum S of the reflection coefficient sequence obtained by deconvolution. d The L2 norm of the vector obtained by subtracting the reconstructed reflection coefficient spectrum Φa, where the second term ||a||1 represents the L1 norm of vector a; S43, the value that meets the constraints of the sparse problem in step S42, is the reflection sparse sequence a obtained by basis pursuit inversion.

4. The high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients according to claim 3, characterized in that, High-resolution seismic data in the time domain is represented as follows: f(t) is the broadband seismic signal after spectral extension.

5. A high-resolution seismic processing system for instantaneous spectral extension of time-frequency domain reflection coefficients, characterized in that, include: The system includes a generalized S-transform module, a reflection coefficient sequence instantaneous spectrum extraction module, a harmonic decomposition module, a broadening module, a high-frequency information enhancement module, and a generalized inverse S-transform module. The generalized S-transform module takes a single seismic data point x(t) as input and outputs the time-frequency spectrum GST(τ,f) corresponding to x(t). The reflection coefficient sequence instantaneous spectrum extraction module takes the time-frequency spectrum GST(τ,f) as input and outputs the reflection coefficient sequence instantaneous spectrum R0. i (f); The input to the harmonic decomposition module is the instantaneous spectrum R of the reflection coefficient sequence. i (f) The output is the instantaneous spectrum R of the reflection coefficient sequence. i (f) is expressed as a weighted sum of sine and cosine functions; the input to the broadened module is the instantaneous spectrum R of the reflection coefficient sequence. i The weighted sum of sine and cosine functions (f) is used as the input to the frequency-extended reflection coefficient sequence spectrum. The high-frequency information enhancement module takes the frequency-extended reflection coefficient sequence spectrum as its input and outputs the frequency-extended time spectrum GST_high(τ,f). The generalized S-inverse transform module takes the frequency-extended time spectrum GST_high(τ,f) as its input and outputs the broadband seismic signal with extended spectrum.

6. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, performs the steps of the high-resolution seismic processing method for instantaneous spectral extension of time-frequency domain reflection coefficients as described in any one of claims 1 to 4.

7. An electronic device, characterized in that, include: The device includes a processor, a memory, and a bus. The memory stores machine-readable instructions executable by the processor. When the electronic device is running, the processor communicates with the memory via the bus. When the machine-readable instructions are executed by the processor, they perform the steps of the seismic high-resolution processing method for instantaneous spectral extension of time-frequency domain reflection coefficients as described in any one of claims 1 to 4.

Images