CN117687086A - Earthquake high-resolution processing method and system, storage medium and electronic equipment - Google Patents
Earthquake high-resolution processing method and system, storage medium and electronic equipment Download PDFInfo
- Publication number
- CN117687086A CN117687086A CN202311677336.7A CN202311677336A CN117687086A CN 117687086 A CN117687086 A CN 117687086A CN 202311677336 A CN202311677336 A CN 202311677336A CN 117687086 A CN117687086 A CN 117687086A
- Authority
- CN
- China
- Prior art keywords
- spectrum
- frequency
- time
- reflection coefficient
- seismic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003672 processing method Methods 0.000 title claims abstract description 12
- 238000001228 spectrum Methods 0.000 claims abstract description 215
- 238000000034 method Methods 0.000 claims abstract description 51
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 17
- 238000012545 processing Methods 0.000 claims abstract description 12
- 230000009466 transformation Effects 0.000 claims description 41
- 230000006870 function Effects 0.000 claims description 28
- 230000001052 transient effect Effects 0.000 claims description 21
- 230000003595 spectral effect Effects 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 9
- 230000002123 temporal effect Effects 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 4
- 230000002035 prolonged effect Effects 0.000 claims description 3
- 238000004590 computer program Methods 0.000 claims description 2
- 102000008297 Nuclear Matrix-Associated Proteins Human genes 0.000 claims 1
- 108010035916 Nuclear Matrix-Associated Proteins Proteins 0.000 claims 1
- 210000000299 nuclear matrix Anatomy 0.000 claims 1
- 238000004458 analytical method Methods 0.000 abstract description 7
- 230000000694 effects Effects 0.000 abstract description 7
- 238000011084 recovery Methods 0.000 abstract description 7
- 239000010410 layer Substances 0.000 description 29
- 238000005516 engineering process Methods 0.000 description 10
- 230000015572 biosynthetic process Effects 0.000 description 8
- 238000013213 extrapolation Methods 0.000 description 8
- 230000002087 whitening effect Effects 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 6
- 238000001914 filtration Methods 0.000 description 6
- 238000010521 absorption reaction Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 239000002131 composite material Substances 0.000 description 2
- 238000013524 data verification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 238000000985 reflectance spectrum Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000013502 data validation Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 239000002344 surface layer Substances 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a seismic high-resolution processing method and system, a storage medium and electronic equipment, which are applied to the field of seismic data processing and aim at the problems that the spectrum prolongation effect is poor and thin-layer seismic is difficult to identify in the prior art; the invention is based on a time-frequency analysis method, deconvolution, harmonic decomposition and base tracking inversion are carried out on the time spectrum of the seismic data, thereby realizing high-frequency reconstruction of the time spectrum, improving the resolution of the seismic data and forming a new spectrum prolongation method. Compared with the prior spectrum recovery method, the method has better spectrum prolongation effect on a complex stratum and higher resolution and signal-to-noise ratio.
Description
Technical Field
The invention belongs to the field of seismic data processing, and particularly relates to a thin-layer seismic identification technology.
Background
With the deep depth of oil and gas seismic exploration, thin-layer seismic identification has gained widespread attention in the oil and gas industry. A thin layer refers to a layer of seismic waves interfering with reflected waves on its top and bottom surfaces, which are often difficult to identify on a seismic profile due to the resolution limitations of the seismic data. The earliest studies of thin layers by Widess et al defined the ultimate resolution of the thin layer as less than one eighth of the wavelength of the seismic wavelet. In consideration of the influence of noise, wavelets and other factors, the geological geophysical expert gate defines the ultimate resolution of a thin layer as a quarter of the wavelet wavelength in practical application. In this sense, a thin layer refers to a formation that is less than a quarter of a seismic wavelet wavelength thick and is not readily discernable in a seismic profile. Therefore, in the process of studying thin layers, improving the resolution of seismic data becomes a non-negligible problem. Studies have shown that the thinner the target layer, the greater the period of the layer notch, requiring a higher seismic data frequency and a wider seismic data frequency band to resolve the thin layer on the seismic profile, i.e., an efficient expansion of the data frequency bandwidth. With the increase of the main frequency and the increase of the frequency bandwidth of the seismic data, the resolution of the seismic data is reasonably improved, so that the broadband seismic data has very important role in thin layer identification.
Because of the geofiltering effect, especially the strong absorption effect of the surface layer, the seismic data collected by the seismic exploration has the characteristics of low frequency and narrow frequency band, and seriously influences the accuracy of thin layer identification. In order to identify the thin layer, the frequency bandwidth of the seismic data is expanded while the main frequency of the seismic data is improved, namely, the high-frequency part information is reasonably compensated while the low-frequency part information is kept. There are many methods for expanding the bandwidth at present, and the methods are mainly classified into 4 types: spectral whitening, deconvolution, absorption compensation, and spectral restoration techniques.
The spectrum whitening method does not change the phase spectrum, but whitens the seismic data amplitude spectrum in a specific frequency band, and compensates the whole frequency. Chen Chuanren and the like combine wavelet transformation with spectral whitening, so that the local frequency characteristics of the signals can be analyzed, wang Jidi shows a Hilbert-Huang transformation-based spectral whitening method, and the local frequency characteristics of the signals in the time domain and the frequency domain can be enhanced simultaneously. Although the spectral whitening method can enhance the high frequency part of the data, the influence of the phase on the resolution of the seismic wavelet is ignored, and the spatial relation of the amplitude of each frequency can be destroyed.
The deconvolution method is based on the seismic deconvolution model, and the reflection coefficient sequence is restored by compressing the seismic wavelet, so that the time resolution of the seismic data is improved, and the method comprises the methods of pulse deconvolution, homomorphic deconvolution, prediction deconvolution, spectrum simulation deconvolution and the like. Robinson provides assumption of convolution model and predictive deconvolution algorithm, and lays the theory and application foundation of seismic deconvolution. Oppenheim proposes homomorphism deconvolution to separate seismic wavelets and reflection coefficients in the complex spectrum by nonlinear filtering. Zhao Bo and the like propose a spectrum simulation deconvolution method, in which the wavelet amplitude spectrum is assumed to be much smoother than the reflection coefficient amplitude spectrum, and the amplitude spectrum of the seismic record can be fitted by using a polynomial, so that the wavelet amplitude spectrum is obtained, the frequency bandwidth of the wavelet amplitude spectrum is widened, and the purpose of improving the resolution of the seismic data is achieved. However, these convolution theory-based methods do not extract effective information fundamentally, and it is difficult to ensure the fidelity and signal-to-noise ratio after processing, and the reflection is not necessarily the true situation of the stratum.
Common methods for absorption compensation techniques are anti-Q filtering and spectral whitening. The anti-Q filtering can compensate the amplitude loss, frequency and phase attenuation of the seismic waves by estimating the stratum quality factor Q so as to improve the resolution of the seismic data. Hale proposes an inverse Q filter algorithm at the earliest, wang proposes a stable and effective inverse Q filter algorithm, which can compensate amplitude attenuation and phase distortion at the same time, and popularize the algorithm to the condition that the Q value changes with time or depth. Wang Shoudong the inverse Q filtering problem is reduced to an inversion problem, and attenuation compensation is achieved by inversion using a regularization method. Although the anti-Q filtering method can compensate the high frequency part to a certain extent, the Q value of the quality factor is generally difficult to accurately calculate, which may lead to inaccurate anti-Q filtering results.
Methods for improving the resolution of the seismic data, such as spectral whitening, deconvolution and absorption compensation, are limited in the cut-off bandwidth of the seismic data, and various mathematical algorithm transformation is utilized to enhance the high-frequency information, so that the obtained high-resolution result has great human randomness and is difficult to reflect the real situation of the stratum. The spectral recovery method predicts reconstructed high frequency components from low frequency components of the seismic data, rather than simply emphasizing the high frequency components. For example, the spectrum recovery sparse layer reflection coefficient inversion technology proposed by Castagna recovers and reconstructs high-frequency information of data from original seismic data, so that the seismic resolution reaches one eighth or even one sixteenth wavelength, and compared with the previous method, the resolution and the fidelity are obviously improved. Harmonic extrapolation carries out harmonic decomposition on the reflection coefficient sequence spectrum in the frequency band of the original seismic data, and carries out spectrum prolongation according to the periodicity of the sine and cosine functions. These spectrum restoration techniques, while capable of widening the frequency band of seismic data, still have the following problems:
(1) Most of frequency spectrum recovery methods including harmonic extrapolation are to analyze frequency spectrums in a frequency domain, however, frequency spectrum characteristics of the whole time domain are obtained by Fourier transformation, and moments corresponding to frequency variation characteristics of signals cannot be described;
(2) The periodicity of the reflectance spectrum is difficult to adequately characterize in the seismic data band in the face of the seismic response created by many thin layer tuning. More importantly, for complex stratum, the reflection coefficient sequence may be a random sequence which is closely arranged, a plurality of periodic functions are overlapped on the corresponding reflection coefficient frequency spectrum, the structure is very complex, and the uncertainty of inversion of the reflection coefficient of the frequency spectrum recovery stratum is enhanced.
Identifying the thin layer needs to widen the frequency band of the seismic data and improve the resolution of the seismic data. In recent years, various high-resolution processing methods are endless, including deconvolution algorithms, absorption compensation techniques and spectrum restoration techniques. However, most frequency expansion methods are limited in the frequency bandwidth of the high signal-to-noise ratio of the seismic data, so that more effective information is not fundamentally increased, the parameter selection is unsuitable, and the fidelity and the signal-to-noise ratio of the seismic data after high-resolution processing are difficult to ensure. Although the spectrum restoration technique represented by harmonic extrapolation can calculate effective high-frequency components from low-frequency components of data, and has better resolution and fidelity, the spectrum prolongation effect is generally poor when the stratum is complex due to the limitation of fourier transform.
Disclosure of Invention
In order to solve the technical problems, the invention provides a seismic high-resolution processing method for instantaneous spectrum prolongation of a time-frequency domain reflection coefficient, which carries out harmonic decomposition and spectrum prolongation on a time-frequency spectrum with a simpler structure on the reflection coefficient sequence, recovers high-frequency information while maintaining low-frequency information, and obviously improves the accuracy of spectrum prolongation.
The invention adopts the technical scheme that: a seismic high-resolution processing method for temporal spectrum prolongation of time-frequency domain reflection coefficient includes:
s1, inputting seismic data x (t) and performing generalized S transformation; the time spectrum GST (τ, f) is obtained:
wherein f represents frequencies from 0 to Nyquist frequency, t is time position, τ represents center time of the Gaussian window, and λ and p are parameters;
s2, in the time-frequency domain, taking each time point t of the time spectrum GST (τ, f) i The corresponding amplitude spectrum, denoted GST (t i F) is carried out; amplitude spectrum GST (t) i Dividing f) by the seismic spectrum W (f) to obtain the instantaneous spectrum R of the reflection coefficient sequence in the amplitude spectrum bandwidth of each time point i (f)；
S3, transient spectrum R of reflection coefficient sequence i (f) Harmonic decomposition is carried out to obtain a sine and cosine function weighted sum:
s4, broadening an instantaneous spectrum of the reflection coefficient sequence by using a base tracking inversion method;
s5, multiplying the transient spectrum of the reflection coefficient sequence after frequency expansion by the wavelet spectrum with the main frequency half of the bandwidth after the expansion, and carrying out frequency expansion on the amplitude spectrum at each moment; obtaining a time spectrum GST_high (tau, f) with enhanced high-frequency information;
s6, performing generalized S inverse transformation on the time spectrum GST_high (tau, f) with the enhanced high-frequency information to obtain high-resolution seismic data in a time domain.
W (f) is specifically estimated from the seismic data or a Rake wavelet of the same zero phase as the primary frequency of the seismic signal is selected.
The invention has the beneficial effects that: the invention provides a seismic high-resolution processing technology for prolonging the transient spectrum of a time-frequency domain reflection coefficient based on generalized S transformation and harmonic decomposition theory, which can improve the resolution of seismic data. Experimental results of synthetic data and actual data demonstrate that: the reflection coefficient spectrum is subjected to harmonic decomposition and base tracking inversion on the basis of the time spectrum, so that the seismic data with higher main frequency can be obtained, the frequency band is effectively widened, the frequency band is prolonged while the effective frequency band inner components of the original seismic data are maintained, the resolution of the seismic data is obviously improved, and favorable seismic data conditions are provided for thin-layer identification and thin-layer oil and gas exploration.
Drawings
FIG. 1 is a two-layer reflectance model;
FIG. 2 is a flow chart of the method of the present invention;
FIG. 3 is a graph showing the comparison of resolvable effects before and after spectrum continuation based on generalized S-transform of synthetic data;
wherein, (a) the primary data has a primary frequency of 30 Hz; (b) The main frequency is increased to 60Hz based on a generalized S transformation spectrum continuation method; (c) Is synthetic seismic data obtained by convolution of a reflection coefficient model and 30Hz Rake wavelets; (d) The method is a result of improving the main frequency of the synthesized data to 60Hz based on a generalized S transformation spectrum continuation method;
FIG. 4 is a graph showing the comparison of the 5 th data before and after the prolongation of the spectrum based on generalized S transformation;
(a) Is the original data section; (b) is a profile based on generalized S-transform spectrum prolongation; (c) Comparing the amplitude spectrum of the data before and after the generalized S-transform spectrum prolongation of the data of the 5 th channel;
FIG. 5 is a spectrum comparison chart of generalized S transformation before and after the spectrum prolongation of the data in the 5 th channel;
wherein, (a) is the spectrum of the 5 th channel of original data when in generalized S transformation; (b) a time spectrum after extending the dominant frequency to 60Hz;
FIG. 6 is a graph of the comparison result of the generalized S-transform spectrum prolongation based on the actual data and the harmonic extrapolation based;
wherein (a) is the original data profile; (b) is a profile based on generalized S-transform spectrum prolongation; (c) Carrying out spectrum continuation on actual data for generalized S transformation before and after data amplitude spectrum comparison;
FIG. 7 shows a seismic high resolution processing system with temporal spectrum extension of time-frequency domain reflection coefficients.
Detailed Description
To facilitate understanding of the technical content of the present invention by those skilled in the art, the following techniques are first described:
1. reflection coefficient model
For a two-layer reflectance formation model as shown in FIG. 1, the reflectance in the time domain is expressed as:
g(t)＝r 1 δ(t-t 1 )+r 2 δ(t-t 1 -T) (1)
wherein: r is (r) 1 Is the top reflection coefficient; r is (r) 2 Is the bottom reflection coefficient; delta (t) is a unit pulse function; t is the time position, t 1 Is the time position of the top reflection coefficient, and T is the time thickness of the formation. Locating the analysis point at the midpoint of the formation, the reflection coefficient of that point can be expressed as:
reflection coefficient e of top-bottom interface 1 And r 2 When the amplitudes are equal and the signs are the same, this layer is called even pulse pair, when r 1 And r 2 When the amplitudes are equal and the signs are opposite, this layer is referred to as an odd pulse pair.
Fourier transforming equation (2) can obtain the spectrums of the even pulse pair and the odd pulse pair as follows:
II(f)＝2rcos(πΔtf) (3)
I I (f)＝i2rsin(πΔtf) (4)
wherein: i is an imaginary unit; r is the reflection coefficient; Δt is the time thickness of the formation; f includes all frequencies within the range of 0 to nyquist frequency. Any sequence of reflection coefficients can be considered as a superposition of even and odd pulse pairs at the analysis point.
In general, the smaller the formation thickness, the greater the period of the notch. The thin layer is only distinguishable when at least one notch period is included in the seismic data band, so that the thinner the formation, the larger the notch period, and the wider the seismic data band required to distinguish this formation.
2. Harmonic extrapolation
In the convolution model, the seismic record in the time domain can be represented as:
s(t)＝w(t)*r(t) (5)
wherein s (t) is a seismic record, w (t) is a seismic wavelet, and r (t) is a reflection coefficient sequence.
The seismic record is represented in the frequency domain as:
S(f)＝W(f)×R(f) (6)
wherein: s (f) is a seismic record spectrum, W (f) is a seismic wavelet spectrum, and R (f) is a reflectance sequence spectrum.
The harmonic extrapolation method first needs deconvolution of the spectrum in the original seismic data bandwidth, namely dividing the seismic data spectrum by the seismic wavelet spectrum, so as to obtain the reflection coefficient sequence spectrum in the seismic data band. The real part of the deconvolution result is a superposition of cosine functions of different periodicity and the imaginary part is a superposition of sine functions of different periodicity. And then according to the periodicity of the sine and cosine functions, decomposing the reflection coefficient sequence frequency spectrum in the frequency band into a form of a weighted sum of a plurality of sine and cosine functions through harmonic decomposition. The reflection coefficient sequence with 2n+1 sampling points corresponds to k=n+1 pulse pairs, and the spectrum of any one reflection coefficient pair can be expressed as the spectrum of even pulse pairs and odd pulse pairs:
G(f,n)＝2r e ·cos(2π·n·dt·f)+i2r o (n)·sin(2π·n·dt·f)) (7)
wherein: dt is the sampling rate, n is half the time thickness of the reflection coefficient pair, r e And r o Coefficients of the spectrum of the even pulse pair and the spectrum of the odd pulse pair, respectively. Thus, the spectrum over any sequence of reflection coefficients from 0 to nyquist frequencies within the analysis window can be expressed as:
converting equation (8) into a matrix form:
S d ＝φa+ε (9)
wherein: s is S d Representing the deconvolution of the resulting reflectance sequence spectrum, phi being a matrix of nuclei consisting of sinusoidal elements, a being a matrix comprising r e And r o Epsilon is the prediction residual. The dimension of the kernel matrix is m×n, M is the number of frequencies determined by the available bandwidth and the frequency sampling rate, and K is the number of pulse pairs.
In order to improve inversion stability, the invention solves the formula (9) by using a base tracking inversion method, and decomposes the data spectrum in the original bandwidth into sparse superposition of stratum responses to obtain an output solution r e (n) and r o (n) corresponding to the real and imaginary parts of the reflection coefficient spectrum, respectively. According to the sine and cosine function periodicity, the reconstructed reflection coefficient spectrum is widened to be outside the bandwidth of the original data, and corresponding K=1 to N+1 (N is the Nyquist frequency) in the summation formula are summed to contain all possible frequency periods, so that the frequency in the reflection coefficient sequence from 0 to the Nyquist frequency is obtained.
The base-trace inversion method solves the coefficients in equation (5) by minimizing the l2 norm of the error term and the l1 norm of the solution at the same time:
wherein: lambda is a regularization parameter for sparsity constraint. Increasing the value of λ yields a higher sparsity result.
The spectrum of the reflection coefficient sequence after the spectrum is widened is multiplied by the wavelet spectrum with higher main frequency and wider frequency band, and the Fourier inverse transformation is carried out to obtain the seismic data with wider frequency band and higher resolution.
3. Time-frequency analysis-generalized S-transform
The S-transform is improved on the basis of short-time Fourier transform and wavelet transform, and a Gaussian window function with width changing inversely with frequency is used while the phase information of the signal is maintained, so that the time-frequency resolution of frequency self-adaptive change is realized. The S-transform of the non-stationary signal h (t) is defined as:
where f represents all frequencies within 0 to nyquist frequency, t is a time variable, τ represents the center time of the gaussian window, ω (t) represents the gaussian window function:
standard deviation is defined as the inverse of frequency:
|f| represents the absolute value of f;
in order to overcome the defect of fixed variation trend of the S-transform Gaussian window function, the S-transform window function is improved. The time window should be narrower because of the strong signal variation in the high frequency part and relatively wider because of the relatively smaller time period and relatively smooth signal variation in the low frequency part. Therefore, in order to more flexibly adjust the window function, parameters mu and p are introduced into the Gaussian window function, so that the variation trend of the window function omega (t) along with the frequency scale f can be adjusted, and the generalized S transformation is obtained. The standard deviation is given by:
the expression of the generalized S transform is as follows:
the window function of the generalized S transformation has nonlinear variation characteristics along with the transformation of frequency, and has better flexibility and higher resolution, and the time-frequency focusing performance is superior to that of the S transformation.
Example 1
Compared with the traditional Fourier transform, the time-frequency analysis is used as a non-stationary signal processing method, and can describe the instantaneous change characteristics and rules of the frequency spectrum energy of the signal in time. On the time spectrum, the sparsity and regularity of the reflection coefficient spectrum obtained by calculation are better. The time-frequency analysis method mainly comprises short-time Fourier transformation, continuous wavelet transformation, S transformation, generalized S transformation and the like. The generalized S transformation can flexibly adjust the width and attenuation trend of the window function, and has better time-frequency focusing property. The generalized S transformation is combined with harmonic decomposition, and the seismic spectrum prolongation is carried out on the time-frequency domain of the seismic data, so that higher resolution can be obtained. The invention provides a seismic high-resolution processing technology for instantaneous spectrum prolongation of a time-frequency domain reflection coefficient based on a time-frequency analysis and harmonic extrapolation spectrum recovery technology. The technology of the invention converts the object with spectrum recovery from complex Fourier transform spectrum into transient spectrum with simpler structure, and carries out harmonic decomposition and spectrum prolongation on the amplitude spectrum corresponding to each moment of the transient spectrum, thereby realizing effective prediction and compensation on high-frequency information of the seismic data, keeping the low-frequency information unchanged, achieving the purpose of improving the resolution of the seismic data and laying a data foundation for thin-layer identification.
The invention firstly carries out generalized S transformation on the seismic signal to obtain a time spectrum with high time-frequency focusing property, carries out time-varying wavelet deconvolution and harmonic decomposition on an amplitude spectrum corresponding to each moment in the time spectrum, decomposes a reflection coefficient sequence transient spectrum in a high signal-to-noise ratio frequency band into a form of a sine-cosine function weighted sum, expands the reflection coefficient sequence transient spectrum to be out of a limited frequency band by utilizing base tracking inversion, multiplies the reflection coefficient sequence transient spectrum with a wavelet spectrum with higher main frequency, and carries out generalized S inverse transformation to obtain a broadband high-resolution seismic signal, so that seismic interpreters can better identify the spatial distribution of a thin layer. As shown in fig. 2, the following procedure is included:
step1, inputting a seismic data, and performing generalized S transformation
First, the seismic signal x (t) is subjected to a generalized S-transform according to equation (16):
where f represents all frequencies within the range of 0 to the nyquist frequency and τ represents the center time of the gaussian window. Adjusting parameters p and μ can adjust the width and attenuation trend of the window function, and appropriate values of p and μ are selected to optimize the focusing and resolution of the spectrum during the generalized S-transform, resulting in a time spectrum GST (τ, f).
Step2 for t i Deconvolution of time-of-day amplitude spectrum to obtain sequence spectrum of reflection coefficient in frequency band
In the time-frequency domain, the generalized S-transform time spectrum GST (τ, f) is taken for each time point t i The corresponding amplitude spectrum, denoted GST (t i F). Amplitude spectrum GST (t) i Dividing f) by the spectrum W (f) (which can be estimated from the seismic data or can be a Rake wavelet with zero phase same as the main frequency of the seismic signal), the instantaneous spectrum R of the reflection coefficient sequence in the amplitude spectrum bandwidth of each time point can be obtained i (f)。
Step3, carrying out harmonic decomposition on the transient spectrum of the reflection coefficient sequence to obtain a sine and cosine function weighted sum
Since the spectrum of any one reflection coefficient pair can be represented by the spectrum of even pulse pair and odd pulse pair, the reflection coefficient sequence spectrum in the bandwidth is decomposed into the form of weighted sum of sine and cosine functions through harmonic decomposition:
step4, broadening the reflectance spectrum using the basis pursuit inversion method
Converting equation (17) into a matrix form:
S d ＝φa+ε (18)
wherein: s is S d Representing the deconvolution of the resulting reflectance sequence spectrum, phi being a matrix of nuclei consisting of sinusoidal elements, a being a matrix comprising r e And r o Epsilon is the prediction residual. L for base tracking method 1 Norm instead of L 0 The norm solves the optimization problem, describing equation (18) as a sparse solution under constraint:
wherein: first term of constraint d -Φa|| 2 L representing a vector obtained by subtracting the deconvolution of the reflection coefficient sequence spectrum from the reconstructed reflection coefficient spectrum 2 Norm, second term a|| 1 L representing vector a 1 Norms. The value conforming to the constraint of the formula (19) is the reflection coefficient vector a obtained by the basis tracking inversion. λ is a regularization parameter, typically ranging from 0 to 1, where increasing λ produces a more sparsely structured structure, reducing λ may amplify inversion noise, typically taking a value of 0.01.
Decomposing the transient spectrum of the reflection coefficient sequence into sparse superposition of stratum response by a basis tracking inversion method, and outputting a weighting coefficient r of a sine-cosine function e (n) and r o (n) real and imaginary parts of the instantaneous spectrum, respectively. Calculating the transient spectrum of the reflection coefficient sequence from 0 to the Nyquist frequency through the formula (17), and widening the transient spectrum of the reflection coefficient sequence in the bandwidth to be outside the original bandwidth.
Step 5. The transient spectrum of the extended reflection coefficient sequence is multiplied by the wavelet spectrum of a certain dominant frequency
Multiplying the extended reflection coefficient sequence spectrum with a wavelet spectrum with a dominant frequency half of the extended frequency bandwidth, so as to obtain a new spectrum i The amplitude spectrum at the moment is frequency-extended.
Repeating steps from Step2 to Step5, realizing frequency expansion of the instantaneous amplitude spectrum at each moment, namely realizing frequency expansion of the whole time spectrum, and obtaining the time spectrum GST_high (tau, f) with enhanced high-frequency information.
GST high(τ,f) ＝{R 1 (f)·W high(f) ，…，R i (f)·W_high(f)，…R N (f)·W_high(f)} (20)
Wherein R is i (f) For t calculated in Step4 i The instantaneous spectrum of the reflection coefficient sequence from 0 to the Nyquist frequency at the moment, i=1, …, N, N is the sampling point number; w_high (f) is a wavelet spectrum with a dominant frequency that is half the bandwidth after the rubbing.
Step6, performing generalized S inverse transformation on the frequency spectrum after frequency expansion
Because the generalized S transformation is lossless and reversible, the generalized S inverse transformation is realized on the time spectrum after frequency expansion, and the high-resolution seismic data in the time domain are obtained:
f (t) is the broadband seismic signal after spectrum extension.
And repeating Step1 to Step6 aiming at each seismic data of the research work area, so as to realize spectrum prolongation of the whole seismic data volume, and finally obtaining the high-resolution seismic data volume. The technology of the invention is not only suitable for improving the resolution of post-stack seismic data, but also can improve the resolution of pre-stack seismic data.
The technical effects of the invention are explained by experiments as follows:
synthetic data validation
The method of the invention is validated by establishing a wedge model. Fig. 3 (a) and 3 (b) are a velocity model and a reflectance model, respectively, used in synthesizing seismic data. FIG. 3 (c) is synthetic seismic data obtained by convolving a reflectance model with a 30Hz Rake wavelet, with a target layer thickness of 0ms on the left side of the model, 30ms on the far right side, and 1ms sampling interval. Fig. 3 (d) is a result of the spectrum continuation method based on the generalized S transform after the dominant frequency of the synthesized data is raised to 60 Hz. The comparison of fig. 3 (c) and fig. 3 (d) shows that only the time thickness of the target stratum in the original synthetic data is greater than 7ms, and the seismic high-resolution data extended by the transient spectrum of the time-frequency domain reflection coefficient can resolve the stratum with the time thickness of 3 ms. Fig. 4 (a) and fig. 4 (b) show cross-sectional views of the original data and the spectrum extended data, respectively, and compared with the velocity model and the reflection coefficient model, it can be found that the high-resolution cross-section after the transient spectrum extension of the time-frequency domain reflection coefficient not only improves the dominant frequency, but also has better amplitude preservation. Fig. 4 (c) compares the fourier amplitude spectra before and after the frequency extension of the 5 th composite data, and can observe that the effective bandwidth of the composite data amplitude spectrum is obviously increased after the instantaneous spectrum extension of the reflection coefficient of the time-frequency domain, and the components in the effective bandwidth of the original seismic data are maintained while the high-frequency components are enhanced.
Fig. 5 compares the generalized S-transform time spectrum before and after the spectrum extension of the 5 th data, fig. 5 (a) is the generalized S-transform time spectrum of the original data, and fig. 5 (b) is the generalized S-transform time spectrum after the main frequency is extended to 60 Hz. It can be seen that the information of the high frequency part of the time spectrum after frequency expansion is obviously enhanced, and the low frequency information is not damaged, thus proving the effectiveness of the technology.
Actual data verification
The method of the invention is subjected to actual data verification. Fig. 6 compares the result graph of the generalized S-transform spectrum extension and the harmonic extrapolation, fig. 6 (a) is an original cross section, fig. 6 (b) is a cross section based on the generalized S-transform spectrum extension, and fig. 6 (c) shows a spectrum comparison graph before and after the generalized S-transform spectrum extension of the 300 th actual data. After the generalized S transformation is carried out for spectrum prolongation, the main frequency of the seismic data is improved, the frequency bandwidth is increased, the resolution of the whole section is increased, the thin layer which is difficult to identify originally is clearer, and the obtained section has better amplitude-preserving property and signal to noise ratio. The frequency band of the seismic data is widened to high frequency, the high frequency components are obviously enhanced, and the components in the effective bandwidth of the original seismic data are not damaged.
The experimental results of the synthetic data and the actual data prove that the high-resolution processing technology of the earthquake provided by the invention for prolonging the transient spectrum of the time-frequency domain reflection coefficient improves the resolution of the earthquake data, widens the frequency band of the earthquake data, enhances the high-frequency component and simultaneously well maintains the component in the effective bandwidth of the original earthquake data, and has higher signal-to-noise ratio and amplitude retention compared with the prior method technology.
Example 2
The invention also provides a seismic high-resolution processing system for the temporal spectrum prolongation of the time-frequency domain reflection coefficient, as shown in fig. 7, comprising: the device comprises a generalized S transformation module, a reflection coefficient sequence instantaneous spectrum extraction module, a harmonic decomposition module, a widening module, a high-frequency information enhancement module and a generalized S inverse transformation module; the generalized S transformation module is input into a channel of seismic data x (t) and outputs a time spectrum GST (tau, f) corresponding to the x (t); the instantaneous spectrum extraction module of the reflection coefficient sequence inputs the time spectrum GST (tau, f) and outputs the instantaneous spectrum R of the reflection coefficient sequence i (f) The method comprises the steps of carrying out a first treatment on the surface of the The harmonic decomposition module inputs as the instantaneous spectrum R of the reflection coefficient sequence i (f) Output is the instantaneous spectrum R of the reflection coefficient sequence i (f) A sine and cosine function weighted sum representation of (a); the broadening module inputs as the instantaneous spectrum R of the reflection coefficient sequence i (f) The sine and cosine functions of (a) are weighted and expressed, and output is a reflection coefficient sequence spectrum after frequency expansion; the input of the high-frequency information enhancement module is a reflection coefficient sequence spectrum after frequency expansion, and the output is a time spectrum GST_high (tau, f) after high-frequency information enhancement; the generalized S inverse transformation module is input into a time spectrum GST_high (tau, f) with enhanced high-frequency information, and outputs a broadband seismic signal with prolonged frequency spectrum.
Example 3
The present embodiment also provides a computer readable storage medium having a computer program stored thereon, which when executed by a processor performs the steps of the seismic high resolution processing method for temporal spectrum continuation of time-frequency domain reflection coefficients described above.
Example 4
The embodiment also provides an electronic device, including: the system comprises a processor, a memory and a bus, wherein the memory stores machine-readable instructions executable by the processor, the processor and the memory are communicated through the bus when the electronic equipment runs, and the machine-readable instructions are executed by the processor to execute the steps of the earthquake high-resolution processing method for prolonging the transient spectrum of the time-frequency domain reflection coefficient.
Those of ordinary skill in the art will recognize that the embodiments described herein are for the purpose of aiding the reader in understanding the principles of the present invention and should be understood that the scope of the invention is not limited to such specific statements and embodiments. Various modifications and variations of the present invention will be apparent to those skilled in the art. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention should be included in the scope of the claims of the present invention.
Claims (7)
1. A seismic high-resolution processing method for temporal spectrum prolongation of time-frequency domain reflection coefficient is characterized by comprising the following steps:
s1, inputting seismic data x (t) and performing generalized S transformation; the time spectrum GST (τ, f) is obtained:
wherein f represents the frequency from 0 to the Nyquist frequency, t is the time position, τ represents the center time of the Gaussian window, and μ, p are parameters;
s2, in the time-frequency domain, taking each time point t of the time spectrum GST (τ, f) i The corresponding amplitude spectrum, denoted GST (t i F) is carried out; amplitude spectrum GST (t) i Dividing f) by the seismic spectrum W (f) to obtain the instantaneous spectrum R of the reflection coefficient sequence in the amplitude spectrum bandwidth of each time point i (f)；
S3, transient spectrum R of reflection coefficient sequence i (f) Harmonic decomposition is carried out to obtain a sine and cosine function weighted sum:
where dt is the sampling rate, n is half the time thickness of the reflection coefficient pair, r e (n) and r o (n) coefficients of the spectrum of the even pulse pair and the spectrum of the odd pulse pair, respectively;
s4, broadening an instantaneous spectrum of the reflection coefficient sequence by using a base tracking inversion method;
s5, multiplying the transient spectrum of the reflection coefficient sequence after frequency expansion by the wavelet spectrum with the main frequency half of the bandwidth after the expansion, and carrying out frequency expansion on the amplitude spectrum at each moment; obtaining a time spectrum GST_high (tau, f) with enhanced high-frequency information;
s6, performing generalized S inverse transformation on the time spectrum GST_high (tau, f) with the enhanced high-frequency information to obtain high-resolution seismic data in a time domain.
2. The method of claim 1, wherein in step S2, W (f) is estimated from the seismic data, or a rake wavelet with zero phase identical to the main frequency of the seismic signal is selected.
3. The method for seismic high-resolution processing of temporal spectrum extension of reflection coefficients in time-frequency domain according to claim 1, wherein step S4 comprises the following steps:
s41, transient spectrum R of reflection coefficient sequence i (f) Is converted into a matrix form:
S d ＝φa+ε
wherein S is d Representing a reflection coefficient sequence spectrum obtained by deconvolution, wherein phi is a nuclear matrix composed of sine elements, a is a coefficient vector, and epsilon is a prediction residual;
s42, use L 1 Norm instead of L 0 The matrix form of step S41 is specifically described as a sparse solution problem under constraint conditions:
where λ is the regularization parameter, the first term of the constraint is S d -Φa|| 2 Representing the deconvoluted reflectance sequence spectrum S d L of vector obtained by subtracting the reconstructed reflection coefficient spectrum phia 2 Norm, second term a|| 1 L representing vector a 1 A norm;
s43, a value conforming to the constraint condition of the sparse solution problem in the step S42 is the reflection sparse sequence a obtained by the basis tracking inversion.
4. A method of seismic high resolution processing for temporal-frequency domain reflection coefficient transient spectrum prolongation according to claim 3, wherein the high resolution seismic data in the time domain are represented as:
f (t) is the broadband seismic signal after spectrum extension.
5. A seismic high-resolution processing system for temporal spectrum prolongation of time-frequency domain reflection coefficients, comprising: the device comprises a generalized S transformation module, a reflection coefficient sequence instantaneous spectrum extraction module, a harmonic decomposition module, a widening module, a high-frequency information enhancement module and a generalized S inverse transformation module; the generalized S transformation module is input into a channel of seismic data x (t) and outputs a time spectrum GST (tau, f) corresponding to the x (t); the instantaneous spectrum extraction module of the reflection coefficient sequence inputs the time spectrum GST (tau, f) and outputs the instantaneous spectrum R of the reflection coefficient sequence i (f) The method comprises the steps of carrying out a first treatment on the surface of the The harmonic decomposition module inputs as the instantaneous spectrum R of the reflection coefficient sequence i (f) Output is the instantaneous spectrum R of the reflection coefficient sequence i (f) A sine and cosine function weighted sum representation of (a); the broadening module inputs as the instantaneous spectrum R of the reflection coefficient sequence i (f) The sine and cosine functions of (a) are weighted and expressed, and output is a reflection coefficient sequence spectrum after frequency expansion; the input of the high-frequency information enhancement module is after frequency expansionThe reflection coefficient sequence spectrum of (2) is output as a time spectrum GST_high (tau, f) after the high-frequency information is enhanced; the generalized S inverse transformation module is input into a time spectrum GST_high (tau, f) with enhanced high-frequency information, and outputs a broadband seismic signal with prolonged frequency spectrum.
6. A computer readable storage medium, characterized in that the computer readable storage medium has stored thereon a computer program which, when executed by a processor, performs the steps of the seismic high resolution processing method of temporal spectral extension of reflection coefficients in the time-frequency domain according to any of claims 1 to 4.
7. An electronic device, comprising: a processor, a memory and a bus, said memory storing machine readable instructions executable by said processor, said processor in communication with said memory via the bus when the electronic device is running, said machine readable instructions when executed by said processor performing the steps of the seismic high resolution processing method of temporal spectral extension of reflection coefficients in the time-frequency domain as defined in any one of claims 1 to 4.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311677336.7A CN117687086A (en) | 2023-12-08 | 2023-12-08 | Earthquake high-resolution processing method and system, storage medium and electronic equipment |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311677336.7A CN117687086A (en) | 2023-12-08 | 2023-12-08 | Earthquake high-resolution processing method and system, storage medium and electronic equipment |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117687086A true CN117687086A (en) | 2024-03-12 |
Family
ID=90131155
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311677336.7A Pending CN117687086A (en) | 2023-12-08 | 2023-12-08 | Earthquake high-resolution processing method and system, storage medium and electronic equipment |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117687086A (en) |
-
2023
- 2023-12-08 CN CN202311677336.7A patent/CN117687086A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Sun et al. | Cross-correlation analysis and time delay estimation of a homologous micro-seismic signal based on the Hilbert–Huang transform | |
| CN109031422A (en) | A kind of seismic signal noise suppressing method based on CEEMDAN and Savitzky-Golay filtering | |
| Chen et al. | Enhancing seismic reflections using empirical mode decomposition in the flattened domain | |
| CN106680876B (en) | A kind of seismic data joint denoising method | |
| CN111427088A (en) | Seismic data low-frequency compensation method for identifying thin mutual reservoir | |
| CN117111155B (en) | Microseism data denoising method based on integrated framework | |
| CN111474582B (en) | Precise S transformation method for generating high-precision time frequency spectrum | |
| CN117687086A (en) | Earthquake high-resolution processing method and system, storage medium and electronic equipment | |
| Liu et al. | Adaptive time-reassigned synchrosqueezing transform for seismic random noise suppression | |
| CN110749923A (en) | Deconvolution method for improving resolution based on norm equation | |
| Wu et al. | Iterative deblending based on the modified singular spectrum analysis | |
| Kravchenko et al. | New analytical WA-systems of Kravchenko functions | |
| Karslı et al. | Post-stack high-resolution deconvolution using Cauchy norm regularization with FX filter weighting | |
| Nose-Filho et al. | Algorithms for sparse multichannel blind deconvolution | |
| CN111856559B (en) | Multi-channel seismic spectrum inversion method and system based on sparse Bayes learning theory | |
| Pan et al. | The Interplay of Framelet Transform and l p Quasi-Norm to Interpolate Seismic Data | |
| Weishi et al. | Statistical denoising of signals in the S-transform domain | |
| Chen et al. | Nonstationary spectral inversion of seismic data | |
| CN112925013A (en) | Seismic data high-resolution processing method based on full-band continuation fidelity | |
| CN112327354A (en) | Method and device for improving low-frequency weak signal, electronic equipment and readable medium | |
| Zheng et al. | Time-Frequency Domain Deconvolution based on Synchrosqueezing Generalized S Transform | |
| CN113419275B (en) | High-resolution seismic processing method based on sparse dictionary learning | |
| He et al. | Seismic Random Noise Simultaneous Attenuation in the Time–Frequency Domain Using Lp-Variation and γ Norm Constraint | |
| CN113359187B (en) | Wavelet sidelobe elimination method for seismic data | |
| Li et al. | Enhancement of the seismic data resolution through Q-compensated denoising based on dictionary learning |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |