High-resolution seismic data synthesis method, device, electronic equipment and medium
By constructing guided filtering and matching pursuit wavelet decomposition and extrapolating low-frequency and high-frequency wavelets, the problems of noise pollution and insufficient low-frequency information in seismic data processing in existing technologies are solved, the synthesis of high-resolution seismic data is achieved, and the signal-to-noise ratio and bandwidth are improved.
Patent Information
- Application Number
- CN202210182568.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-02-25
- Publication Date
- 2025-09-26
- Estimated Expiration
- 2042-02-25
AI Technical Summary
Existing seismic data processing methods suffer from severe noise pollution when improving high-frequency information and insufficient expansion of low-frequency information, resulting in low signal-to-noise ratio of seismic profiles and insufficient vertical and horizontal resolution.
By constructing guided filtering and matching pursuit wavelet decomposition, the main frequency, low frequency and high frequency ranges of the original seismic data are determined, smoothing filtering and wavelet decomposition are performed, low frequency and high frequency wavelets are extrapolated, and high-resolution seismic data are synthesized.
While maintaining the original signal-to-noise ratio, the resolution of seismic data is improved, the relative bandwidth is expanded, the ability to depict stratigraphic details is enhanced, and the vertical and lateral resolutions are improved.
Smart Images

Figure CN116699679B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the field of geophysical exploration, and more particularly to a method, device, electronic equipment and medium for synthesizing high-resolution seismic data. Background Art
[0002] With the increasing exploration and development of oil and gas fields, thin interbeds and small targets have become crucial components of upstream reserve growth. High-quality seismic data are crucial for discovering these hidden reservoirs. Therefore, improving the resolution of seismic data has been a hot topic of research for scholars both domestically and internationally. Deconvolution was first proposed by the MIT Geophysical Analysis Group in the 1960s. Based on the Robinson model, the inverse operator was first calculated and then used to determine the reflection coefficient series. Subsequently, through continuous effort and experimentation, various deconvolution methods were developed to meet diverse needs. Peacock and Treitel pioneered predictive deconvolution in 1969. In 1975, Burg deconvolution was proposed, building on Burg's 1967 maximum entropy analysis method. That same year, Makhoul proposed the linear predictive deconvolution algorithm, building on predictive deconvolution. In 1992, Xu Boxun et al., building on previous research, proposed the adaptive predictive deconvolution algorithm. In 1994, Wang Chengshu proposed multi-channel predictive deconvolution. In 1971, Ulrych proposed homomorphic deconvolution, which circumvents the assumptions of minimum phase of seismic wavelets and white noise of reflection coefficients, allowing for the simultaneous extraction of seismic wavelets and reflection coefficients. In 2003, Gao Shaowu et al. improved homomorphic deconvolution by employing the L1 norm, the Parsimony criterion, and minimum entropy deconvolution to select the optimal wavelet. This improved method significantly improved both processing effectiveness and efficiency. In 1995, Ling Yun et al. developed zero-phase homomorphic deconvolution based on homomorphic theory. This processing method can significantly improve the vertical resolution of seismic profiles.
[0003] Margrave first proposed time-varying deconvolution in 1998, and subsequently in 2011, he used non-stationary deconvolution to estimate reflection coefficients using the Gabor transform, breaking through traditional time-domain or frequency-domain deconvolution methods. Conventional deconvolution methods use the least-squares matching principle to design an anti-wavelet filter, which allows seismic data to output reflection coefficient pulses after passing through the filter, thereby eliminating the band-limited filtering characteristics of the seismic wavelet. This least-squares deconvolution method assumes that the seismic wavelet is a minimum-phase wavelet and the reflection coefficient is a white noise sequence. However, these assumptions are difficult to hold in most cases and cannot accurately characterize the reflection characteristics of thin layers. In 2011, Gary et al. proposed Gabor deconvolution, generalizing the stable convolution model to non-stationary convolution. This method represents the seismic record as three components: an attenuation function, a source wavelet, and a reflection coefficient. Similar to the Fourier decomposition of stable convolution, this method approximates the decomposition of the seismic trace in the time-frequency domain as the Gabor transform of the reflection coefficient multiplied by the attenuation function, which is then multiplied by the Fourier transform of the source wavelet. Under the assumptions of white noise reflection coefficients and minimum phase wavelets in the time-frequency domain, the attenuation function and source wavelet are estimated directly from the time-frequency spectrum of the seismic trace. The time-frequency spectrum of the reflection coefficient is obtained through mathematical operations, and then the inverse transform is used to obtain the reflection coefficient in the time domain. This algorithm leverages the concept of inverse Q filtering and applies hyperbolic smoothing to directly estimate the attenuation function from the time-frequency spectrum of the seismic trace, addressing the nonstationarity of seismic records and the instability of the inverse Q filtering algorithm. However, this algorithm uses the Gabor transform, which involves the issue of Gaussian window length truncation. Furthermore, the hyperbolic smoothing assumes a uniformly attenuating medium, so the actual data often produce some artifacts.
[0004] Deconvolution techniques, based on white noise reflection coefficients and the assumption of known wavelet phases, aim to achieve time-domain wavelet compression. This typically elevates high-frequency information, but high-frequency information in actual seismic signals is heavily contaminated by noise. This elevating of high-frequency signals also increases noise, resulting in low signal-to-noise ratios (SNRs) and discontinuous events in the horizontal direction of the processed seismic profile. Furthermore, low-frequency information in the seismic signal is not expanded, resulting in poor vertical relative amplitude relationships and limited relative bandwidth expansion of the processed seismic signal.
[0005] Therefore, it is necessary to develop a high-resolution seismic data synthesis method, device, electronic equipment and medium.
[0006] The information disclosed in the background technology section of the present invention is only intended to deepen the understanding of the general background technology of the present invention, and should not be regarded as an admission or any form of suggestion that the information constitutes the prior art already known to those skilled in the art. Summary of the Invention
[0007] The present invention proposes a high-resolution seismic data synthesis method, device, electronic device and medium, which can improve seismic resolution by constructing guided filtering and matching pursuit wavelet decomposition, while maintaining the original seismic signal-to-noise ratio, and extrapolating high- and low-frequency wavelets.
[0008] In a first aspect, an embodiment of the present disclosure provides a method for synthesizing high-resolution seismic data, comprising:
[0009] Determine the main frequency, low frequency and high frequency ranges of the original seismic data;
[0010] Performing structure-guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio;
[0011] performing matching pursuit wavelet decomposition on the filtered seismic data to obtain a plurality of seismic wavelets;
[0012] Extrapolating, based on the dominant frequencies of the plurality of seismic wavelets, to obtain low-frequency wavelets and high-frequency wavelets that match the dominant frequency, low-frequency range, and high-frequency range of the original seismic data;
[0013] The main frequencies of the multiple seismic wavelets, the low-frequency wavelets and the high-frequency wavelets are superimposed to synthesize high-resolution seismic data.
[0014] Preferably, determining the main frequency, low frequency and high frequency ranges of the original seismic data includes:
[0015] In the time domain, a short-time window Fourier transform is performed on each trace of original seismic data, and the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges.
[0016] Preferably, performing structure-guided smoothing filtering on the raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio comprises:
[0017] Determine tensor diffusion expressions;
[0018] Calculate the structure tensor and then calculate the diffusion tensor;
[0019] Substituting the diffusion tensor into the tensor diffusion expression and performing discretization to obtain a filtering iteration formula;
[0020] The original seismic data is subjected to structure-guided smoothing filtering by the filtering iterative formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0021] Preferably, the filtering iteration formula is:
[0022]
[0023] Among them, u k and u k+1They represent the filtering results of the original seismic data at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration, div is the divergence operator, is the gradient operator, and D represents the diffusion tensor.
[0024] Preferably, matching pursuit wavelet decomposition is performed using formula (2):
[0025]
[0026] Where, f(t) is the filtered seismic data, a n is the seismic wavelet w of the nth iteration n The amplitude of (t), R (N) f is the residual signal energy after N iterations, and N is the total number of iterations.
[0027] Preferably, extrapolating the main frequencies of the plurality of seismic wavelets to obtain low-frequency wavelets and high-frequency wavelets that match the main frequency, low-frequency, and high-frequency ranges of the original seismic data comprises:
[0028] Initialize the extrapolation parameters of each seismic wavelet separately;
[0029] Obtaining the corresponding low-frequency wavelet and high-frequency wavelet according to the main frequency of the seismic wavelet and the corresponding extrapolation parameter;
[0030] The main frequencies of the multiple seismic wavelets, the spectra of the low-frequency wavelets and the high-frequency wavelets are superimposed, compared with the spectrum of the original seismic data, and the extrapolation parameters are adjusted to obtain the low-frequency wavelet and the high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency ranges of the original seismic data.
[0031] Preferably, the low-frequency wavelet is:
[0032]
[0033] The high-frequency wavelet is:
[0034] ω″=αω0 (4)
[0035] Among them, ω′ is the low-frequency wavelet, ω0 is the main frequency of the seismic wavelet, α is the extrapolation parameter, and ω″ is the high-frequency wavelet.
[0036] As a specific implementation of the embodiment of the present disclosure,
[0037] In a second aspect, the present disclosure also provides a high-resolution seismic data synthesis device, comprising:
[0038] Frequency determination module, which determines the main frequency, low frequency and high frequency range of the original seismic data;
[0039] A filtering module, performing structure-guided smoothing filtering on the raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio;
[0040] a wavelet decomposition module, performing matching pursuit wavelet decomposition on the filtered seismic data to obtain a plurality of seismic wavelets;
[0041] An extrapolation module, which extrapolates, based on the dominant frequencies of the plurality of seismic wavelets, to obtain low-frequency and high-frequency wavelets that match the dominant frequency, low-frequency and high-frequency ranges of the original seismic data;
[0042] The synthesis module superimposes the main frequencies of the plurality of seismic wavelets with the low-frequency and high-frequency wavelets obtained by extrapolation to synthesize high-resolution seismic data.
[0043] Preferably, determining the main frequency, low frequency and high frequency ranges of the original seismic data includes:
[0044] In the time domain, a short-time window Fourier transform is performed on each trace of original seismic data, and the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges.
[0045] Preferably, performing structure-guided smoothing filtering on the raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio comprises:
[0046] Determine tensor diffusion expressions;
[0047] Calculate the structure tensor and then calculate the diffusion tensor;
[0048] Substituting the diffusion tensor into the tensor diffusion expression and performing discretization to obtain a filtering iteration formula;
[0049] The original seismic data is subjected to structure-guided smoothing filtering by the filtering iterative formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0050] Preferably, the filtering iteration formula is:
[0051]
[0052] Among them, u k and u k+1 They represent the filtering results of the original seismic data at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration, div is the divergence operator, is the gradient operator, and D represents the diffusion tensor.
[0053] Preferably, matching pursuit wavelet decomposition is performed using formula (2):
[0054]
[0055] Where, f(t) is the filtered seismic data, a n is the seismic wavelet w of the nth iteration n The amplitude of (t), R (N) f is the residual signal energy after N iterations, and N is the total number of iterations.
[0056] Preferably, extrapolating the main frequencies of the plurality of seismic wavelets to obtain low-frequency wavelets and high-frequency wavelets that match the main frequency, low-frequency, and high-frequency ranges of the original seismic data comprises:
[0057] Initialize the extrapolation parameters of each seismic wavelet separately;
[0058] Obtaining the corresponding low-frequency wavelet and high-frequency wavelet according to the main frequency of the seismic wavelet and the corresponding extrapolation parameter;
[0059] The main frequencies of the multiple seismic wavelets, the spectra of the low-frequency wavelets and the high-frequency wavelets are superimposed, compared with the spectrum of the original seismic data, and the extrapolation parameters are adjusted to obtain the low-frequency wavelet and the high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency ranges of the original seismic data.
[0060] Preferably, the low-frequency wavelet is:
[0061]
[0062] The high-frequency wavelet is:
[0063] ω″=αω0 (4)
[0064] Among them, ω′ is the low-frequency wavelet, ω0 is the main frequency of the seismic wavelet, α is the extrapolation parameter, and ω″ is the high-frequency wavelet.
[0065] In a third aspect, an embodiment of the present disclosure further provides an electronic device, the electronic device comprising:
[0066] a memory storing executable instructions;
[0067] A processor runs the executable instructions in the memory to implement the high-resolution seismic data synthesis method.
[0068] In a fourth aspect, an embodiment of the present disclosure further provides a computer-readable storage medium, which stores a computer program, and when the computer program is executed by a processor, the high-resolution seismic data synthesis method is implemented.
[0069] The methods and apparatus of the present invention have other features and advantages that will be apparent from or will be described in detail in the accompanying drawings and subsequent detailed descriptions incorporated herein, which together serve to explain the specific principles of the invention. BRIEF DESCRIPTION OF THE DRAWINGS
[0070] The above and other objects, features and advantages of the present invention will become more apparent through a more detailed description of exemplary embodiments of the present invention with reference to the accompanying drawings, wherein like reference numerals generally represent like components throughout the exemplary embodiments of the present invention.
[0071] Figure 1 A flow chart showing the steps of a high-resolution seismic data synthesis method according to one embodiment of the present invention.
[0072] Figure 2 A schematic diagram illustrating a seismic section of raw seismic data according to an embodiment of the present invention is shown.
[0073] Figure 3 A schematic diagram illustrating a seismic cross section of high-resolution seismic data according to one embodiment of the present invention.
[0074] Figure 4 A block diagram of a high-resolution seismic data synthesis device according to an embodiment of the present invention is shown.
[0075] Description of reference numerals:
[0076] 201. Frequency determination module; 202. Filtering module; 203. Wavelet decomposition module; 204. Extrapolation module; 205. Synthesis module. DETAILED DESCRIPTION
[0077] The preferred embodiments of the present invention will be described in more detail below. Although the preferred embodiments of the present invention are described below, it should be understood that the present invention can be implemented in various forms and should not be limited to the embodiments set forth herein.
[0078] The present invention provides a high-resolution seismic data synthesis method, comprising:
[0079] Determine the main frequency, low frequency and high frequency ranges of the original seismic data;
[0080] Perform structurally guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio;
[0081] Perform matching pursuit wavelet decomposition on filtered seismic data to obtain multiple seismic wavelets;
[0082] According to the main frequencies of multiple seismic wavelets, low-frequency and high-frequency wavelets that match the main frequency, low-frequency and high-frequency ranges of the original seismic data are obtained by extrapolation;
[0083] The main frequencies of multiple seismic wavelets are superimposed with the low-frequency and high-frequency wavelets obtained by extrapolation to synthesize high-resolution seismic data.
[0084] In one example, determining the main frequency, low frequency, and high frequency ranges of the raw seismic data includes:
[0085] In the time domain, a short-time window Fourier transform is performed on each trace of original seismic data, and the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges.
[0086] In one example, performing structure-guided smoothing filtering on raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio includes:
[0087] Determine tensor diffusion expressions;
[0088] Calculate the structure tensor and then calculate the diffusion tensor;
[0089] Substitute the diffusion tensor into the tensor diffusion expression and discretize it to obtain the filtering iteration formula;
[0090] The original seismic data is subjected to structure-guided smoothing filtering through a filtering iterative formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0091] In one example, the filtering iteration formula is:
[0092]
[0093] Among them, u k and u k+1 They represent the filtering results of the original seismic data at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration, div is the divergence operator, is the gradient operator, and D represents the diffusion tensor.
[0094] In one example, matching pursuit wavelet decomposition is performed using formula (2):
[0095]
[0096] Where, f(t) is the filtered seismic data, a n is the seismic wavelet w of the nth iteration n The amplitude of (t), R (N) f is the residual signal energy after N iterations, N is the total number of iterations, and is also the number of seismic wavelets.
[0097] In one example, extrapolating the main frequencies of the plurality of seismic wavelets to obtain low-frequency wavelets and high-frequency wavelets that match the main frequency, low-frequency, and high-frequency ranges of the original seismic data includes:
[0098] Initialize the extrapolation parameters of each seismic wavelet separately;
[0099] According to the main frequency of the seismic wavelet and the corresponding extrapolation parameters, the corresponding low-frequency wavelet and high-frequency wavelet are obtained;
[0100] The spectra of the main frequency, low-frequency wavelet and high-frequency wavelet of multiple seismic wavelets are superimposed, compared with the spectrum of the original seismic data, and the extrapolation parameters are adjusted to obtain the low-frequency wavelet and high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency range of the original seismic data.
[0101] In one example, the low frequency wavelet is:
[0102]
[0103] The high-frequency wavelet is:
[0104] ω″=αω0 (4)
[0105] Among them, ω′ is the low-frequency wavelet, ω0 is the main frequency of the seismic wavelet, α is the extrapolation parameter, and ω″ is the high-frequency wavelet.
[0106] Specifically, in the time domain, a short-time window Fourier transform is performed on each piece of original seismic data, and then the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges, which are used for the subsequent extrapolation of high and low frequency wavelets.
[0107] Anisotropic diffusion is used for smoothing filtering to improve the signal-to-noise ratio of the input signal, which is beneficial for the matching pursuit wavelet decomposition of the seismic trace in the later stage. The image of the original seismic data is used as the initial condition, and the diffused image is obtained by solving the partial differential equation about time. In the diffusion equation, the local structural information (faults, pinch-outs, etc.) is obtained by introducing the structural tensor. The diffusion tensor is designed based on this structural information, that is, different diffusion coefficients are used in different directions to achieve the effect of denoising while protecting the edge. The filtering of the anisotropic diffusion equation is realized by the tensor diffusion expression - formula (5):
[0108]
[0109] Among them, u(x,y,0) is the image of the original seismic data, div is the divergence operator, is the gradient operator, D represents the diffusion tensor, and its elements are designed based on the local structure information of the image extracted by the structure tensor S.
[0110] To calculate the diffusion tensor, first calculate the structure tensor (Structure Tensor), which is the first-order partial differential information u of the image x ,u y Another form of expression, which provides a matrix field information, so that each point in the image corresponds to a 3 × 3 (in the case of three-dimensional) real symmetric matrix, such as formula (6):
[0111]
[0112] Among them, G σ Represents the Gaussian kernel with σ as parameter, which avoids the influence of noise when estimating the gradient. ρ Convolution can take surrounding information into account and avoid the disadvantage of canceling out edges with the same direction but opposite signs when orienting the edges.
[0113] The diffusion tensor uses the same eigenvectors as the structure tensor. The diffusion coefficients μ1, μ2, and μ3 are designed based on the local features extracted from the structure tensor. The diffusion behaviors in the directions of υ1, υ2, and υ3 are independently controlled to achieve the effect of reducing noise while enhancing edges. υ1, υ2, and υ3 are the eigenvectors of the structure tensor matrix, obtained by solving the eigenvalues and eigenvectors. Here, the diffusion tensor D is designed as formula (7):
[0114]
[0115] in, It is usually a very small positive number, and in this example, α = 0.001. This design implies the following mechanism: along the υ1 direction, that is, the direction parallel to the gradient, or the direction with the fastest rate of change, the diffusion coefficient is very small and the edge is protected; along the υ2 direction, the diffusion coefficient is adjusted according to the coherence: the greater the coherence, the closer the diffusion coefficient is to the maximum diffusion coefficient of 1. On the contrary, if there is no significant coherence in the image locally, this may be the end of the reflection in the seismic image. At this time, μ2 = α, and the diffusion will also be very slow, thus protecting geological structures such as faults and pinch-outs. Therefore, the elements of the diffusion tensor can be calculated as d 12 =μ1υ 11 υ 12 +μ2υ 21 υ 22 +μ3υ 31 υ 32 , d 13 =μ1υ 11 υ 13 +μ2υ 21 υ 23 +μ3υ31 υ 33 , d 23 =μ1υ 12 υ 13 +μ2υ 22 υ 23 +μ3υ 32 υ 33 .
[0116] Discretize the tensor diffusion expression and obtain the differential filtering iterative formula:
[0117]
[0118] Among them, u k and u k+1 They represent the filtering results of the original image at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration (0.1≤Δt≤0.2 is the best). The noise profile is obtained by subtracting the original seismic data from the filtered seismic data. The number of iterations is determined according to the characteristics of the noise profile, usually 5-8 times.
[0119] In the iterative process, the differential is replaced by the available difference, and the filtering iteration formula is formula (1). The original seismic data is subjected to a construction-guided smoothing filter using the filtering iteration formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0120] The filtered seismic data is decomposed by matching pursuit wavelet using formula (2) to obtain multiple seismic wavelets. The seismic wavelets can use Morlet wavelet, as shown in formula (9), but are not limited to this wavelet:
[0121]
[0122] Among them, ω0 is the main frequency.
[0123] Initialize the extrapolation parameters of each seismic wavelet separately; calculate the corresponding low-frequency wavelet and high-frequency wavelet by formula (3) and formula (4) respectively according to the main frequency of the seismic wavelet and the corresponding extrapolation parameters; superimpose the frequency spectra of the main frequency, low-frequency wavelet and high-frequency wavelet of multiple seismic wavelets, compare them with the frequency spectrum of the original seismic data, adjust the extrapolation parameters, and obtain the low-frequency wavelet and high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency range of the original seismic data.
[0124] The main frequencies of multiple seismic wavelets are superimposed with their corresponding low- and high-frequency wavelets to synthesize high-resolution seismic data. Because the extrapolated wavelets are noise-free, the signal-to-noise ratio of the resulting high-resolution profile is not reduced, while the relative bandwidth of the processed seismic signal is expanded.
[0125] The present invention also provides a high-resolution seismic data synthesis device, comprising:
[0126] Frequency determination module, which determines the main frequency, low frequency and high frequency range of the original seismic data;
[0127] The filtering module performs structure-guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio;
[0128] The wavelet decomposition module performs matching pursuit wavelet decomposition on the filtered seismic data to obtain multiple seismic wavelets;
[0129] An extrapolation module, based on the dominant frequencies of multiple seismic wavelets, extrapolates to obtain low-frequency and high-frequency wavelets that match the dominant frequency, low-frequency and high-frequency ranges of the original seismic data;
[0130] The synthesis module superimposes the main frequencies of multiple seismic wavelets with the low-frequency and high-frequency wavelets obtained by extrapolation to synthesize high-resolution seismic data.
[0131] In one example, determining the main frequency, low frequency, and high frequency ranges of the raw seismic data includes:
[0132] In the time domain, a short-time window Fourier transform is performed on each trace of original seismic data, and the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges.
[0133] In one example, performing structure-guided smoothing filtering on raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio includes:
[0134] Determine tensor diffusion expressions;
[0135] Calculate the structure tensor and then calculate the diffusion tensor;
[0136] Substitute the diffusion tensor into the tensor diffusion expression and discretize it to obtain the filtering iteration formula;
[0137] The original seismic data is subjected to structure-guided smoothing filtering through a filtering iterative formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0138] In one example, the filtering iteration formula is:
[0139]
[0140] Among them, u k and u k+1 They represent the filtering results of the original seismic data at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration, div is the divergence operator, is the gradient operator, and D represents the diffusion tensor.
[0141] In one example, matching pursuit wavelet decomposition is performed using formula (2):
[0142]
[0143] Where, f(t) is the filtered seismic data, a n is the seismic wavelet w of the nth iteration n The amplitude of (t), R (N) f is the residual signal energy after N iterations, N is the total number of iterations, and is also the number of seismic wavelets.
[0144] In one example, extrapolating the main frequencies of the plurality of seismic wavelets to obtain low-frequency wavelets and high-frequency wavelets that match the main frequency, low-frequency, and high-frequency ranges of the original seismic data includes:
[0145] Initialize the extrapolation parameters of each seismic wavelet separately;
[0146] According to the main frequency of the seismic wavelet and the corresponding extrapolation parameters, the corresponding low-frequency wavelet and high-frequency wavelet are obtained;
[0147] The spectra of the main frequency, low-frequency wavelet and high-frequency wavelet of multiple seismic wavelets are superimposed, compared with the spectrum of the original seismic data, and the extrapolation parameters are adjusted to obtain the low-frequency wavelet and high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency range of the original seismic data.
[0148] In one example, the low frequency wavelet is:
[0149]
[0150] The high-frequency wavelet is:
[0151] ω″=αω0 (4)
[0152] Among them, ω′ is the low-frequency wavelet, ω0 is the main frequency of the seismic wavelet, α is the extrapolation parameter, and ω″ is the high-frequency wavelet.
[0153] Specifically, in the time domain, a short-time window Fourier transform is performed on each piece of original seismic data, and then the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges, which are used for the subsequent extrapolation of high and low frequency wavelets.
[0154] The anisotropic diffusion process is used for smoothing filtering to improve the signal-to-noise ratio of the input signal, which is beneficial to the matching pursuit wavelet decomposition of the seismic trace in the later stage. The image of the original seismic data is used as the initial condition, and the diffused image is obtained by solving the partial differential equation about time. In the diffusion equation, the local structural information (faults, pinch-outs, etc.) is obtained by introducing the structural tensor. The diffusion tensor is designed based on this structural information, that is, different diffusion coefficients are used in different directions to achieve the effect of denoising while protecting the edge. The filtering of the anisotropic diffusion equation is realized by the tensor diffusion expression - formula (5),
[0155]
[0156] Among them, u(x,y,0) is the image of the original seismic data, div is the divergence operator, is the gradient operator, D represents the diffusion tensor, and its elements are designed based on the local structure information of the image extracted by the structure tensor S.
[0157] To calculate the diffusion tensor, first calculate the structure tensor (Structure Tensor), which is the first-order partial differential information u of the image x ,u y Another form of expression, which provides a matrix field information, so that each point in the image corresponds to a 3 × 3 (in the case of three-dimensional) real symmetric matrix, such as formula (6):
[0158]
[0159] Among them, G σ Represents the Gaussian kernel with σ as parameter, which avoids the influence of noise when estimating the gradient. ρ Convolution can take surrounding information into account and avoid the disadvantage of canceling out edges with the same direction but opposite signs when orienting the edges.
[0160] The diffusion tensor uses the same eigenvectors as the structure tensor. The diffusion coefficients μ1, μ2, and μ3 are designed based on the local features extracted from the structure tensor. The diffusion behavior in the directions of υ1, υ2, and υ3 is independently controlled, thereby achieving the effect of reducing noise while enhancing edges. υ1, υ2, and υ3 are the eigenvectors of the structure tensor matrix, obtained by solving the eigenvalues and eigenvectors. Here, the diffusion tensor D is designed as formula (7):
[0161]
[0162] in, It is usually a very small positive number, and in this example, α = 0.001. This design implies the following mechanism: along the υ1 direction, that is, the direction parallel to the gradient, or the direction with the fastest rate of change, the diffusion coefficient is very small and the edge is protected; along the υ2 direction, the diffusion coefficient is adjusted according to the coherence: the greater the coherence, the closer the diffusion coefficient is to the maximum diffusion coefficient of 1. On the contrary, if there is no significant coherence in the image locally, this may be the end of the reflection in the seismic image. At this time, μ2 = α, and the diffusion will also be very slow, thus protecting geological structures such as faults and pinch-outs. Therefore, the elements of the diffusion tensor can be calculated as d 12 =μ1υ 11 υ 12 +μ2υ 21 υ 22 +μ3υ 31 υ 32 , d 13 =μ1υ 11 υ 13 +μ2υ 21 υ 23 +μ3υ 31 υ 33 , d 23 =μ1υ 12 υ 13 +μ2υ 22 υ 23 +μ3υ 32 υ 33 .
[0163] Discretize the tensor diffusion expression and obtain the differential filtering iterative formula:
[0164]
[0165] Among them, u k and u k+1 The values represent the filtering results of the original image at kΔt and (k+1)Δt, respectively. Δt is the diffusion time for one iteration (0.1≤Δt≤0.2 is the best choice). The noise profile is obtained by subtracting the original seismic data from the filtered data. The number of iterations is determined based on the characteristics of the noise profile, typically 5-8.
[0166] In the iterative process, the differential is replaced by the available difference, and the filtering iteration formula is formula (1). The original seismic data is subjected to a construction-guided smoothing filter using the filtering iteration formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0167] The filtered seismic data is decomposed by matching pursuit wavelet using formula (2) to obtain multiple seismic wavelets. The seismic wavelets can use Morlet wavelet, as shown in formula (9), but are not limited to this wavelet:
[0168]
[0169] Among them, ω0 is the main frequency.
[0170] Initialize the extrapolation parameters of each seismic wavelet separately; calculate the corresponding low-frequency wavelet and high-frequency wavelet by formula (3) and formula (4) respectively according to the main frequency of the seismic wavelet and the corresponding extrapolation parameters; superimpose the frequency spectra of the main frequency, low-frequency wavelet and high-frequency wavelet of multiple seismic wavelets, compare them with the frequency spectrum of the original seismic data, adjust the extrapolation parameters, and obtain the low-frequency wavelet and high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency range of the original seismic data.
[0171] The main frequencies of multiple seismic wavelets are superimposed with their corresponding low- and high-frequency wavelets to synthesize high-resolution seismic data. Because the extrapolated wavelets are noise-free, the signal-to-noise ratio of the resulting high-resolution profile remains unchanged, while the relative bandwidth of the processed seismic signal is expanded.
[0172] The present invention also provides an electronic device, which includes: a memory storing executable instructions; and a processor running the executable instructions in the memory to implement the above-mentioned high-resolution seismic data synthesis method.
[0173] The present invention also provides a computer-readable storage medium storing a computer program, which implements the above-mentioned high-resolution seismic data synthesis method when executed by a processor.
[0174] To facilitate understanding of the solutions and effects of the embodiments of the present invention, four specific application examples are given below. Those skilled in the art should understand that these examples are only for facilitating understanding of the present invention, and any specific details thereof are not intended to limit the present invention in any way.
[0175] Example 1
[0176] Figure 1 A flow chart showing the steps of a high-resolution seismic data synthesis method according to one embodiment of the present invention.
[0177] like Figure 1As shown, the high-resolution seismic data synthesis method includes: step 101, determining the main frequency, low frequency and high frequency range of the original seismic data; step 102, performing structure-guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio; step 103, performing matching pursuit wavelet decomposition on the filtered seismic data to obtain multiple seismic wavelets; step 104, extrapolating the main frequencies of the multiple seismic wavelets to obtain low-frequency and high-frequency wavelets that match the main frequency, low frequency and high frequency range of the original seismic data; step 105, superimposing the main frequencies of the multiple seismic wavelets with the low-frequency and high-frequency wavelets obtained by extrapolation to synthesize high-resolution seismic data.
[0178] Figure 2 A schematic diagram illustrating a seismic section of raw seismic data according to an embodiment of the present invention is shown.
[0179] Figure 3 A schematic diagram illustrating a seismic cross section of high-resolution seismic data according to one embodiment of the present invention.
[0180] contrast Figure 2 、 Figure 3 It can be seen that the processed seismic resolution is significantly improved, and the signal-to-noise ratio is high. In terms of details, the strong trough at the bottom of the seismic profile corresponds to the coal seams of the Shanxi and Taiyuan Formations, overlying the reservoirs of the Shihezi Formation. On the pre-processed seismic profile, two layers are vaguely visible on the left, merging into a single layer to the right, indicating that the bottom boundary of the Shihezi Formation is unclear from the top boundary of the Shanxi Formation. After processing, the two sets of coaxial lines are separated and have good continuity, which can better track the bottom boundary of the lower Shihezi Formation.
[0181] It can be seen that the processed seismic profile not only has improved resolution in the vertical direction and continuous event axes in the horizontal direction, but also has no reduction in signal-to-noise ratio and maintains a good original relative amplitude relationship; the energy in the high-frequency band is compensated and the energy in the low-frequency band is increased, thereby widening the relative effective frequency band of the seismic data, better depicting the stratigraphic details, and effectively improving the resolution of the seismic data.
[0182] Example 2
[0183] Figure 4 A block diagram of a high-resolution seismic data synthesis device according to an embodiment of the present invention is shown.
[0184] like Figure 4 As shown, the high-resolution seismic data synthesis device comprises:
[0185] Frequency determination module 201, determines the main frequency, low frequency and high frequency range of the original seismic data;
[0186] The filtering module 202 performs structure-guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio;
[0187] The wavelet decomposition module 203 performs matching pursuit wavelet decomposition on the filtered seismic data to obtain multiple seismic wavelets;
[0188] An extrapolation module 204, based on the dominant frequencies of the plurality of seismic wavelets, extrapolates to obtain low-frequency and high-frequency wavelets that match the dominant frequency, low-frequency and high-frequency ranges of the original seismic data;
[0189] The synthesis module 205 superimposes the main frequencies of the multiple seismic wavelets with the low-frequency and high-frequency wavelets obtained by extrapolation to synthesize high-resolution seismic data.
[0190] As an optional method, determining the main frequency, low frequency and high frequency ranges of the original seismic data includes:
[0191] In the time domain, a short-time window Fourier transform is performed on each trace of original seismic data, and the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges.
[0192] As an optional solution, structurally guided smoothing filtering is performed on the raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio. The following methods are used:
[0193] Determine tensor diffusion expressions;
[0194] Calculate the structure tensor and then calculate the diffusion tensor;
[0195] Substitute the diffusion tensor into the tensor diffusion expression and discretize it to obtain the filtering iteration formula;
[0196] The original seismic data is subjected to structure-guided smoothing filtering through a filtering iterative formula to obtain filtered seismic data with a high signal-to-noise ratio.
[0197] As an optional solution, the filtering iteration formula is:
[0198]
[0199] Among them, u k and u k+1 They represent the filtering results of the original seismic data at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration, div is the divergence operator, is the gradient operator, and D represents the diffusion tensor.
[0200] As an alternative, matching pursuit wavelet decomposition is performed using formula (2):
[0201]
[0202] Where, f(t) is the filtered seismic data, a n is the seismic wavelet w of the nth iteration n The amplitude of (t), R (N) f is the residual signal energy after N iterations, and N is the total number of iterations.
[0203] As an optional solution, based on the dominant frequencies of multiple seismic wavelets, low-frequency wavelets and high-frequency wavelets that match the dominant frequency, low-frequency, and high-frequency ranges of the original seismic data are obtained by extrapolation, including:
[0204] Initialize the extrapolation parameters of each seismic wavelet separately;
[0205] According to the main frequency of the seismic wavelet and the corresponding extrapolation parameters, the corresponding low-frequency wavelet and high-frequency wavelet are obtained;
[0206] The spectra of the main frequency, low-frequency wavelet and high-frequency wavelet of multiple seismic wavelets are superimposed, compared with the spectrum of the original seismic data, and the extrapolation parameters are adjusted to obtain the low-frequency wavelet and high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency range of the original seismic data.
[0207] As an alternative, the low-frequency wavelet is:
[0208]
[0209] The high-frequency wavelet is:
[0210] ω″=αω0 (4)
[0211] Among them, ω′ is the low-frequency wavelet, ω0 is the main frequency of the seismic wavelet, α is the extrapolation parameter, and ω″ is the high-frequency wavelet.
[0212] Example 3
[0213] The present disclosure provides an electronic device comprising: a memory storing executable instructions; and a processor executing the executable instructions in the memory to implement the above-mentioned high-resolution seismic data synthesis method.
[0214] An electronic device according to an embodiment of the present disclosure includes a memory and a processor.
[0215] The memory is used to store non-transitory computer-readable instructions. Specifically, the memory may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. The volatile memory may include, for example, random access memory (RAM) and / or cache memory. The non-volatile memory may include, for example, read-only memory (ROM), a hard disk, flash memory, etc.
[0216] The processor may be a central processing unit (CPU) or other form of processing unit having data processing capability and / or instruction execution capability, and may control other components in the electronic device to perform desired functions. In one embodiment of the present disclosure, the processor is used to execute the computer-readable instructions stored in the memory.
[0217] Those skilled in the art should understand that in order to solve the technical problem of how to obtain a good user experience, this embodiment may also include well-known structures such as a communication bus and an interface, and these well-known structures should also be included in the scope of protection of this disclosure.
[0218] For detailed description of this embodiment, please refer to the corresponding description in the aforementioned embodiments, which will not be repeated here.
[0219] Example 4
[0220] An embodiment of the present disclosure provides a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, the high-resolution seismic data synthesis method is implemented.
[0221] According to an embodiment of the present disclosure, a computer-readable storage medium stores non-transitory computer-readable instructions, which, when executed by a processor, execute all or part of the steps of the aforementioned methods of the embodiments of the present disclosure.
[0222] The above-mentioned computer-readable storage media include, but are not limited to, optical storage media (e.g., CD-ROMs and DVDs), magneto-optical storage media (e.g., MOs), magnetic storage media (e.g., magnetic tapes or mobile hard disks), media with built-in rewritable non-volatile memory (e.g., memory cards), and media with built-in ROM (e.g., ROM cartridges).
[0223] Those skilled in the art should understand that the above description of the embodiments of the present invention is only for the purpose of illustrative purposes only to illustrate the beneficial effects of the embodiments of the present invention, and is not intended to limit the embodiments of the present invention to any given examples.
[0224] While various embodiments of the present invention have been described above, the above description is intended to be illustrative, not exhaustive, and not limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments.
Claims
1. A high-resolution seismic data synthesis method, characterized in that: include: Determine the main frequency, low frequency and high frequency ranges of the original seismic data; Performing structure-guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio; performing matching pursuit wavelet decomposition on the filtered seismic data to obtain a plurality of seismic wavelets; Extrapolating, based on the dominant frequencies of the plurality of seismic wavelets, to obtain low-frequency wavelets and high-frequency wavelets that match the dominant frequency, low-frequency range, and high-frequency range of the original seismic data; The main frequencies of the multiple seismic wavelets, the low-frequency wavelets and the high-frequency wavelets are superimposed to synthesize high-resolution seismic data.
2. The high-resolution seismic data synthesis method according to claim 1, wherein: Determine the main frequency, low frequency and high frequency ranges of the original seismic data including: In the time domain, a short-time window Fourier transform is performed on each trace of original seismic data, and the frequency information of the statistical seismic signal is superimposed to obtain the main frequency, low frequency and high frequency ranges.
3. The high-resolution seismic data synthesis method according to claim 1, wherein: Performing structure-guided smoothing filtering on the original seismic data to obtain filtered seismic data with a high signal-to-noise ratio includes: Determine tensor diffusion expressions; Calculate the structure tensor and then calculate the diffusion tensor; Substituting the diffusion tensor into the tensor diffusion expression and performing discretization to obtain a filtering iteration formula; The original seismic data is subjected to structure-guided smoothing filtering by the filtering iterative formula to obtain filtered seismic data with a high signal-to-noise ratio.
4. The high-resolution seismic data synthesis method according to claim 3, wherein: The filtering iteration formula is: Among them, u k and u k+1 They represent the filtering results of the original seismic data at kΔt and (k+1)Δt, respectively. Δt is the diffusion time of one iteration, div is the divergence operator, is the gradient operator, and D represents the diffusion tensor.
5. The high-resolution seismic data synthesis method according to claim 1, wherein: Perform matching pursuit wavelet decomposition using formula (2): Where, f(t) is the filtered seismic data, a n is the seismic wavelet w of the nth iteration n The amplitude of (t), R (N) f is the residual signal energy after N iterations, and N is the total number of iterations.
6. The high-resolution seismic data synthesis method according to claim 1, wherein: Extrapolating, based on the dominant frequencies of the plurality of seismic wavelets, to obtain low-frequency wavelets and high-frequency wavelets that match the dominant frequency, low-frequency range, and high-frequency range of the original seismic data comprises: Initialize the extrapolation parameters of each seismic wavelet separately; Obtaining the corresponding low-frequency wavelet and high-frequency wavelet according to the main frequency of the seismic wavelet and the corresponding extrapolation parameter; The main frequencies of the multiple seismic wavelets, the spectra of the low-frequency wavelets and the high-frequency wavelets are superimposed, compared with the spectrum of the original seismic data, and the extrapolation parameters are adjusted to obtain the low-frequency wavelet and the high-frequency wavelet corresponding to each seismic wavelet that matches the main frequency, low-frequency and high-frequency ranges of the original seismic data.
7. The high-resolution seismic data synthesis method according to claim 6, wherein: The low-frequency wavelet is: The high-frequency wavelet is: ω″=αω0 (4) Among them, ω′ is the low-frequency wavelet, ω0 is the main frequency of the seismic wavelet, α is the extrapolation parameter, and ω″ is the high-frequency wavelet.
8. A high-resolution seismic data synthesis device, characterized in that: include: Frequency determination module, which determines the main frequency, low frequency and high frequency range of the original seismic data; A filtering module, performing structure-guided smoothing filtering on the raw seismic data to obtain filtered seismic data with a high signal-to-noise ratio; a wavelet decomposition module, performing matching pursuit wavelet decomposition on the filtered seismic data to obtain a plurality of seismic wavelets; An extrapolation module, which extrapolates, based on the dominant frequencies of the plurality of seismic wavelets, to obtain low-frequency and high-frequency wavelets that match the dominant frequency, low-frequency and high-frequency ranges of the original seismic data; The synthesis module superimposes the main frequencies of the plurality of seismic wavelets with the low-frequency and high-frequency wavelets obtained by extrapolation to synthesize high-resolution seismic data.
9. An electronic device, characterized in that: The electronic device comprises: a memory storing executable instructions; A processor, wherein the processor runs the executable instructions in the memory to implement the high-resolution seismic data synthesis method according to any one of claims 1 to 7.
10. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the high-resolution seismic data synthesis method according to any one of claims 1 to 7.
Citation Information
Patent Citations
Seismic record broadband expanding method
CN104122589A
BPFE method and device for high-resolution processing of seismic data
CN114035225A