CN111175827A - A High-Performance Time-Frequency Domain Filtering Method for Enhanced Seismic Exploration Signals - Google Patents
A High-Performance Time-Frequency Domain Filtering Method for Enhanced Seismic Exploration Signals Download PDFInfo
- Publication number
- CN111175827A CN111175827A CN202010130916.4A CN202010130916A CN111175827A CN 111175827 A CN111175827 A CN 111175827A CN 202010130916 A CN202010130916 A CN 202010130916A CN 111175827 A CN111175827 A CN 111175827A
- Authority
- CN
- China
- Prior art keywords
- frequency
- time
- filtering
- noise
- signal
- 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.)
- Granted
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 101
- 238000000034 method Methods 0.000 title claims abstract description 84
- 238000009826 distribution Methods 0.000 claims abstract description 44
- 230000001360 synchronised effect Effects 0.000 claims abstract description 43
- 238000009792 diffusion process Methods 0.000 claims abstract description 33
- 230000002708 enhancing effect Effects 0.000 claims abstract description 11
- 238000004220 aggregation Methods 0.000 claims abstract description 7
- 238000001125 extrusion Methods 0.000 claims description 48
- 238000001228 spectrum Methods 0.000 claims description 19
- 230000008569 process Effects 0.000 claims description 13
- 230000000694 effects Effects 0.000 claims description 12
- 230000009466 transformation Effects 0.000 claims description 11
- 230000008707 rearrangement Effects 0.000 claims description 6
- 230000008030 elimination Effects 0.000 claims description 4
- 238000003379 elimination reaction Methods 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 239000000654 additive Substances 0.000 claims description 3
- 230000000996 additive effect Effects 0.000 claims description 3
- 230000003247 decreasing effect Effects 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims 1
- 230000009467 reduction Effects 0.000 abstract description 8
- 230000002776 aggregation Effects 0.000 abstract description 3
- 230000001629 suppression Effects 0.000 abstract description 3
- 238000004321 preservation Methods 0.000 abstract description 2
- 238000012545 processing Methods 0.000 description 21
- 238000010586 diagram Methods 0.000 description 12
- 238000002474 experimental method Methods 0.000 description 6
- 230000014759 maintenance of location Effects 0.000 description 6
- 238000003672 processing method Methods 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 4
- 238000004445 quantitative analysis Methods 0.000 description 4
- 230000006872 improvement Effects 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 238000004451 qualitative analysis Methods 0.000 description 3
- 238000010183 spectrum analysis Methods 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000004215 lattice model Methods 0.000 description 2
- 239000000047 product Substances 0.000 description 2
- 238000009825 accumulation Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000011056 performance test Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/324—Filtering
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
一种增强地震勘探信号的高性能时频域滤波方法，先对含噪的地震信号做同步挤压小波变换得到高分辨率的时频平面，再在所得的时频平面上采用Catte模型进行各向异性扩散滤波，然后对滤波后的时频分布做同步挤压小波重构，最后得到有效信号的时域估计即为最终的滤波结果，并对滤波结果进行分析；本发明充分利用了同步挤压时频分布高分辨、高聚集特性和Catte模型各向异性扩散滤波能够兼顾有效信号特征保持与随机噪声消减的特性，提出了一种增强地震勘探信号的高性能时频域滤波方法；经验证，该方法在地震勘探信号的增强以及随机噪声的压制方面具有非常优越的性能。
A high-performance time-frequency domain filtering method for enhancing seismic exploration signals. First, synchronous squeezing wavelet transform is performed on noisy seismic signals to obtain a high-resolution time-frequency plane, and then the Catte model is used on the obtained time-frequency plane to perform various Anisotropic diffusion filtering is performed, and then synchronous squeezing wavelet reconstruction is performed on the filtered time-frequency distribution. Finally, the time domain estimation of the effective signal is the final filtering result, and the filtering result is analyzed. The present invention makes full use of synchronous squeezing. The high resolution and high aggregation characteristics of the pressure time-frequency distribution and the anisotropic diffusion filtering of the Catte model can take into account the characteristics of effective signal feature preservation and random noise reduction. A high-performance time-frequency domain filtering method for enhancing seismic exploration signals is proposed. , the method has very superior performance in the enhancement of seismic exploration signal and the suppression of random noise.
Description
Technical Field
The invention relates to a filtering technology, in particular to a high-performance time-frequency domain filtering method for enhancing seismic exploration signals by anisotropic diffusion filtering in a high-resolution time-frequency domain.
Background
In order to obtain a seismic record with clearer event information, the acquired seismic exploration data must be processed. For seismic data acquired from an easily explorable and exploitable area, noise interference is less and the situation is not very complex, a certain processing requirement can be met by applying a common filtering method, but if a processing result with very clear information of a phase axis is obtained, a more advanced filtering method still needs to be found. With the increase of depth and difficulty of exploration and development, the noise condition in the collected seismic exploration data is more complex and serious. Random noise is one of the important interference types affecting the quality of seismic data, and can obscure, truncate and even annihilate effective in-phase axis information. If the noise interference can be well suppressed, the effective in-phase axis can be greatly enhanced, so that the quality of seismic exploration data is greatly improved.
To date, there is no method for enhancing seismic exploration signals by applying anisotropic diffusion filtering in a high-resolution time-frequency transform domain, namely a synchronous extrusion transform domain. In the existing research result reports, there are various time-frequency domain filtering methods (such as wavelet filtering, S-transform filtering, time-frequency peak filtering (TFPF), etc.) for processing seismic data, and some image processing methods (such as gaussian filtering, anisotropic diffusion filtering, filtering based on curvelet and contourlet transform, etc.) for processing seismic images. Although these studies have achieved significant results, there is still room for improvement. When the traditional time-frequency filtering method is used for reducing random noise of seismic exploration, the filtering effect is difficult to improve to a certain degree due to the limitation of each time-frequency distribution (such as a time window selection problem). When the image processing method is used, if the seismic data are converted into pictures with a certain format for processing, the conversion problem among data with different formats is involved when a single-channel signal is extracted for waveform and frequency spectrum contrast analysis, and is often difficult to solve; more methods need to process the seismic data by zero filling to form a matrix with the number of elements in rows and columns being 2 to the power of an integer, and after the processing is finished, a matrix corresponding to the original data is intercepted from a filtering result, so that the whole processing process is not fast and convenient.
Currently, among many time-frequency filtering methods, time-frequency peak filtering (TFPF) has very good filtering performance. The method can carry out unbiased estimation on signals which linearly change along with time under the condition of Gaussian white noise (the signal-to-noise ratio can be as low as-9 dB). The implementation process comprises the steps of firstly carrying out frequency modulation on a noise-containing signal to obtain an analytic signal of an original signal, then carrying out time-frequency transformation on the analytic signal by utilizing Wigner-Ville distribution (PWVD) or a windowing form thereof, namely pseudo-Wigner-Ville distribution (PWVD), and then carrying out peak value search on a time-frequency plane, wherein an obtained peak value set is an estimated value of an effective signal. The application of the method in seismic exploration signal processing is very mature, and remarkable effects are achieved. As the seismic exploration signals have nonlinearity and non-stationarity, PWVD is suitable for being adopted in the TFPF method, so that the method has certain limitations: the window length is chosen in a trade-off between retention of the effective signal and reduction of random noise. Therefore, the filtering effect of the method is limited. The existing filtering method based on synchronous extrusion time-frequency transformation mainly performs region interception (such as band-pass filtering) or division (such as threshold division) on a time-frequency plane, and then performs synchronous extrusion reconstruction on the time-frequency distribution of a reserved part to obtain an estimation value of an effective signal. The method can effectively remove the interference of different frequencies with the reserved part, but the enhancement effect on the reserved part is limited, so that the filtering effect still has a space for improvement.
The anisotropic diffusion filtering in the image processing method is applied to enhance the seismic image, and the filtering processing is only carried out along the gradient direction of the effective homophase axis, while the homophase axes in the seismic exploration data are various and have different directions, so that the concentrated filtering cannot be carried out in a specific area, and the achieved effect is limited to a certain extent. When a considerable number of image processing methods, such as a filtering method based on curvelet transform, are used to process an image, the number of elements in rows and columns of a matrix needs to be reduced to an integer power of 2, and either blocking and intercepting data or padding zero and completing numbers can cause the processing process to be not concise. In fact, many image processing methods have a good processing effect on picture files with a certain format, if seismic data are converted into a picture format for filtering processing, and then in the analysis of filtering results, a single-channel signal needs to be extracted for waveform comparison and frequency spectrum comparison, some problems related to data format conversion exist, and in most cases, the single-channel signal cannot be converted into the original seismic data format for analysis, so that the method is very inconvenient.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention aims to provide a method for performing anisotropic diffusion filtering in a high-resolution time-frequency domain, namely a high-performance filtering method for enhancing seismic exploration signals by adopting Catte model filtering in a synchronous extrusion time-frequency domain, so that effective signals are better enhanced and noise is better reduced to supplement each other, seismic data with higher definition of a same-phase axis is obtained, and a good premise is provided for subsequent data interpretation and other work.
In order to achieve the purpose, the technical scheme of the invention is as follows:
a high performance time-frequency domain filtering method for enhancing seismic exploration signals, comprising the steps of:
step one, performing synchronous extrusion wavelet transformation on a noise-containing seismic signal to obtain high-resolution time-frequency distribution:
setting the noise-containing signal as:
s(t)＝x(t)+n(t) (1)
x (t) is the effective signal, n (t) is additive random noise;
and (3) performing continuous wavelet transform on s (t), and referring to the formula (2), wherein the continuous wavelet transform comprises the following steps:
where a is the scale factor, b is the translation factor, ψ (t) is the wavelet basis function,
frequency domain table of equation (2) as a conjugate function of ψ (t) according to Pasteval's theoremThe expression is as follows:
where ξ is the angular frequency,
is the Fourier transform of s (t), phi (a ξ) is phi (t), the convolution form of the time domain is transformed into the product form of the frequency domain by the formula, and the instantaneous frequency, namely the instantaneous frequency, can be preliminarily estimated by differentiating the wavelet coefficients
Where i is the imaginary unit; this converts the time-scale plane (b, a) into a time-frequency plane (b, ω)s(a, b)), both the frequency ω and the scale a are discretized, thus only at discrete points ai{ai-ai-1＝(Δa)iCalculating wavelet coefficient Ws(a, b); the implementation process of frequency rearrangement in synchronous extrusion wavelet transformation comprises the following steps: assuming that the sampling frequency of the signal is cf and the number of scales taken during wavelet transform is N, the discrete frequency value on the rearranged time-frequency spectrum is ωl＝l*cf/N,l∈[1,N](ii) a Synchronous extrusion value T of wavelet coefficients(ωlB) by extrusion of any center frequency ωlNearby section
Is obtained by a value of (i.e. a
Wherein, akIs a discrete scale, and Δ ak＝ak-ak-1，Ts(ωlB) is the high-resolution and high-aggregation time-frequency distribution after extrusion rearrangement;
step two, performing Catte model anisotropic diffusion filtering on the high-resolution time-frequency plane, wherein the Catte model is as follows:
where U (x, y, t) represents the pixel value at time (x, y) t, U0(x, y) represents the original image, ▽ () is the divergence operator, ▽ is the gradient operator,. GσIs a gaussian kernel function with a scale factor sigma, C is a diffusion coefficient, and a non-negative monotonically decreasing function is usually chosen. Perona and Malik propose two diffusion coefficient functions:
and
where k is a constant, an appropriate integer |. G can be selected depending on the object to be processedσU (x, y, t) | is the gradient mode of the constant coefficient thermal diffusion equation with the scale factor σ;
the time frequency distribution T obtained by the formula (5) iss(ωlAnd b) regarding the image as an image, performing anisotropic diffusion filtering on the image, and substituting the image into the formula (6) to obtain:
thus, the filtered time-frequency plane is obtained, the energy gathering point of the effective signal is obviously enhanced, the random noise is suppressed to a great extent, and particularly the high-frequency random noise is greatly reduced; if the random noise in the noise-containing signal is strong, a part of noise at the intermediate frequency still remains on the filtered time-frequency plane, so that in order to obtain a better noise elimination effect, time-frequency domain band-pass filtering is needed;
step three, performing synchronous extrusion wavelet reconstruction on the filtered time-frequency distribution to obtain a filtered time-domain waveform:
the time-frequency distribution after the filter of the Catte model can be obtained by the formula (7) and is recorded as Ts′(ωlAnd b), adopting a synchronous extrusion reconstruction formula as follows:
wherein,
re represents the operation of the real part; and x (b) is a time domain waveform reconstructed after the time-frequency surface filtering, namely the estimation of the effective signal x (t) in the formula (1).
The invention firstly provides a method for enhancing seismic exploration signals by adopting a Catte anisotropic diffusion filtering model on a synchronous extrusion time-frequency transform domain, and also provides a method for powerfully suppressing random noise of seismic exploration. The synchronous extrusion wavelet transform is a time-frequency analysis method with quite high time-frequency resolution and time-frequency aggregation at present, the energy of effective signals in the time-frequency distribution obtained through the transform is highly aggregated together, and random noise energy is scattered in a time-frequency plane. And the obtained time-frequency plane is regarded as an image, and the image is subjected to Catte model anisotropic diffusion filtering, so that the energy accumulation point of the effective signal is obviously enhanced, and meanwhile, the random noise energy is reduced to a great extent. And (4) for the filtered time-frequency distribution, estimating the effective signal by adopting a synchronous extrusion wavelet transform reconstruction process.
The method fully utilizes the high resolution and high aggregation characteristics of synchronous extrusion time-frequency transformation and the capability of the Catte model anisotropic diffusion filtering for considering both characteristic retention and noise elimination, so as to provide a more advanced filtering technology for the enhancement of effective signals and the reduction of random noise in seismic exploration data. The advantages of the method provided by the invention are verified by testing the simulated seismic signals and the artificially synthesized seismic records, and a quite ideal result is obtained by subsequent filtering processing on actual seismic data.
Drawings
FIG. 1 is a flow chart of a time-frequency domain anisotropic diffusion filtering method for synchronous squeeze wavelet transform in the present invention.
FIG. 2 is a diagram of a single Ricker wavelet waveform, including both noise-free and noise-added waveforms. FIG. 2(a) is a Ricker wavelet without noise; FIG. 2(b) is a Ricker wavelet with a signal-to-noise ratio of about 15 dB.
FIG. 3 is a single wavelet synchronous extrusion time-frequency distribution diagram, which includes a noise-free and noise-containing wavelet time-frequency diagram and a time-frequency distribution diagram after anisotropic diffusion filtering of a synchronous extrusion time-frequency domain Catte model. FIG. 3(a) is a diagram of a noise-free Ricker wavelet synchronous extrusion time-frequency distribution; FIG. 3(b) is a distribution diagram of the Ricker wavelet synchronous extrusion time-frequency with a signal-to-noise ratio of about 15 dB; FIG. 3(c) is a time-frequency distribution diagram of anisotropic diffusion filtering of time-frequency domain of noise-containing Ricker wavelet synchronous extrusion.
FIG. 4 is a diagram of two-component Ricker wavelet waveforms, including noise-free and noise-added waveforms. FIG. 4(a) is a two-component Ricker wavelet without noise; FIG. 4(b) is a two-component Ricker wavelet with a signal-to-noise ratio of about 5 dB.
FIG. 5 is a two-component Ricker wavelet synchronous extrusion time-frequency distribution diagram, which includes a noise-free and noise-containing wavelet time-frequency diagram and a time-frequency distribution after anisotropic diffusion band-pass filtering of a synchronous extrusion time-frequency domain Catte model. FIG. 5(a) is a diagram of a noise-free two-component Ricker wavelet synchronous extrusion time-frequency distribution; FIG. 5(b) is a two-component Ricker wavelet synchronous extrusion time-frequency distribution graph with a signal-to-noise ratio of about 5 dB; FIG. 5(c) is a time-frequency distribution diagram of anisotropic diffusion band-pass filtering of time-frequency domain of synchronous extrusion of noisy two-component Ricker wavelets.
FIG. 6 is a comparison graph of Time Frequency Peak Filtering (TFPF) and the filtering waveform of the present invention for two component Ricker wavelets and a partial enlarged view thereof. FIG. 6(a) is a comparison of a noisy waveform, a non-noisy waveform, a TFPF waveform, and a new method filtered waveform; FIG. 6(b) is a comparison graph of four waveforms of the first sub-dominant signal segment; FIG. 6(c) is a comparison of four waveforms of the second sub-dominant signal segment; FIG. 6(d) is a comparison graph of four waveforms of the first wavelet non-dominant signal segment; FIG. 6(e) is a comparison of four waveforms for the second wavelet non-dominant signal segment.
FIG. 7 is a graph of a comparison of the frequency spectrum of a TFPF and the method of the present invention after filtering two component Ricker wavelets and a close-up view thereof. FIG. 7(a) is a graph of a comparison of the frequency spectrum of a noisy signal, a non-noisy signal, a TFPF signal, and a new method filtered signal; FIG. 7(b) is a comparison graph of the middle and low frequency and middle frequency partial frequency spectrums; FIG. 7(c) is a graph comparing the spectra of the high frequency part.
FIG. 8 is a graph of the results of a filtering experiment of TFPF and the method of the present invention for artificially synthesized seismic records having a signal-to-noise ratio of about-4 dB. FIG. 8(a) is a synthetic seismic record without noise; FIG. 8(b) is a synthetic seismic record with a signal-to-noise ratio of about-4 dB; FIG. 8(c) is a filtered record of the method of the present invention; FIG. 8(d) is a TFPF record; FIG. 8(e) is a filtered difference record for the method of the present invention; fig. 8(f) is the TFPF difference record.
FIG. 9 is a graph of the results of a filtering experiment of TFPF and the method of the present invention on artificially synthesized seismic records with a signal-to-noise ratio of-12 dB. FIG. 9(a) is a synthetic seismic record without noise; FIG. 9(b) is a synthetic seismic record with a signal-to-noise ratio of about-12 dB; FIG. 9(c) is a filtered record of the method of the present invention; FIG. 9(d) is a TFPF record; FIG. 9(e) is a comparison of the waveform of trace 23; FIG. 9(f) is a comparison of the first wavelet waveform; FIG. 9(g) is a second wavelet waveform comparison.
FIG. 10 is a graph of the results of filtering an actual seismic recording by TFPF and the method of the present invention. FIG. 10(a) is an original record; FIG. 10(b) is a filtered record of the method of the present invention; fig. 10(c) is a TFPF record.
Fig. 11 is a time-frequency distribution graph of a single-channel signal in an actual seismic record after synchronous extrusion time-frequency domain lattice model anisotropic diffusion filtering. FIG. 11(a) is the synchronous extrusion time-frequency distribution of the 49 th signal; FIG. 11(b) is a time-frequency domain distribution of time-frequency domain anisotropic diffusion filtering for the 49 th channel signal; FIG. 11(c) is the synchronous extrusion time-frequency distribution of the 67 th signal; FIG. 11(d) is the time-frequency domain distribution of the time-frequency domain anisotropic diffusion filter for the 67 th channel signal.
FIG. 12 is a graph of a comparison of the spectrum of the 49 th signal in an actual seismic record and its filtered results. FIG. 12(a) is a comparison of the frequency spectra before and after filtering of the whole channel signal; fig. 12(b) is a partially enlarged view of the spectral contrast of the channel signal.
Detailed Description
The technical scheme of the invention is described in detail in the following with reference to the accompanying drawings.
Referring to fig. 1, a method of enhancing high performance time-frequency domain filtering of seismic survey signals, comprising the steps of:
step one, performing synchronous extrusion wavelet transformation on a noise-containing seismic signal to obtain high-resolution time-frequency distribution:
setting the noise-containing signal as:
s(t)＝x(t)+n(t) (1)
x (t) is the effective signal, n (t) is additive random noise;
and (3) performing continuous wavelet transform on s (t), and referring to the formula (2), wherein the continuous wavelet transform comprises the following steps:
where a is the scale factor, b is the translation factor, ψ (t) is the wavelet basis function,
a conjugate function of ψ (t); according to the paseuler theorem, the frequency domain expression of equation (2) is:
where ξ is the angular frequency,
is the Fourier transform of s (t), phi (a ξ) is phi (t), the convolution form of the time domain is transformed into the product form of the frequency domain by the formula, and the instantaneous frequency can be preliminarily estimated by differentiating the wavelet coefficients, i.e. the instantaneous frequency is estimated
Where i is the imaginary unit; this converts the time-scale plane (b, a) into a time-frequency plane (b, ω)s(a, b)), both the frequency ω and the scale a are discretized, thus only at discrete points ai{ai-ai-1＝(Δa)iCalculating wavelet coefficient Ws(a, b); the implementation process of frequency rearrangement in synchronous extrusion wavelet transformation comprises the following steps: assuming a signalThe sampling frequency is cf, the number of scales taken during wavelet transformation is N, and the discrete frequency value is omega on the rearranged time frequency spectruml＝l*cf/N,l∈[1,N](ii) a Synchronous extrusion value T of wavelet coefficients(ωlB) by extrusion of any center frequency ωlNearby section
Is obtained by a value of (i.e. a
Wherein, akIs a discrete scale, and Δ ak＝ak-ak-1，Ts(ωlAnd b) is the high-resolution and high-aggregation time-frequency distribution after extrusion rearrangement.
The resulting time-frequency distribution is shown in fig. 3 and 5. FIGS. 3(a) and (b) are the simultaneous extrusion time-frequency distributions of a single Ricker wavelet (with a dominant frequency of 25Hz) without noise (shown in FIG. 2 (a)) and with noise (shown in FIG. 2 (b)), respectively; FIGS. 5(a) and (b) are the simultaneous extrusion time-frequency distributions of two-component Ricker wavelets (with 20Hz and 15Hz main frequencies, respectively) without noise (shown in FIG. 4 (a)) and with noise (shown in FIG. 4 (b)); it can be seen that no matter whether the random noise is weak or strong, the energy concentration of the effective signal on the time-frequency plane is high, the frequency component curve is fine, and the time-frequency curve is not clear due to interference of a certain degree only when the noise is strong. Comparing fig. 3(b) and fig. 5(b), it can be seen that random noise is more spread in the high frequency part in the latter.
Step two, performing Catte model anisotropic diffusion filtering on the high-resolution time-frequency plane, wherein the Catte model is as follows:
where U (x, y, t) represents the pixel value at time (x, y) t, U0(x, y) represents the original image, ▽ () is the divergence operator, ▽ is the gradient operator,. GσIs a Gaussian kernel function with a scale factor of sigma, C is diffusionCoefficients, usually a non-negative monotonically decreasing function is chosen; perona and Malik propose two diffusion coefficient functions:
and
where k is a constant, an appropriate integer |. G can be selected depending on the object to be processedσU (x, y, t) | is the gradient mode of the constant coefficient thermal diffusion equation with the scale factor σ;
the time frequency distribution T obtained by the formula (5) iss(ωlAnd b) regarding the image as an image, performing anisotropic diffusion filtering on the image, and substituting the image into the formula (6) to obtain:
thus, the filtered time-frequency plane is obtained, the energy gathering point of the effective signal is obviously enhanced, the random noise is suppressed to a great extent, and particularly the high-frequency random noise is greatly reduced; if the random noise in the noisy signal is strong, a part of noise at the intermediate frequency still remains on the filtered time-frequency plane, so that in order to obtain a better noise elimination effect, time-frequency domain band-pass filtering is needed.
As shown in fig. 3(c), the time-frequency distribution after the 15dB single Ricker wavelet synchronous extrusion time-frequency domain cable model is filtered is relatively clean, and no further processing is needed; the two-component Ricker wavelet shown in fig. 4(b) has a relatively low signal-to-noise ratio, and the random noise energy on the synchronous extrusion time-frequency plane (shown in fig. 5 (b)) is relatively strong, so that the filter of the Catte model on the time-frequency plane needs to intercept the enhanced effective signal energy distribution by means of band-pass filtering, as shown in fig. 5 (c).
Step three, performing synchronous extrusion wavelet reconstruction on the filtered time-frequency distribution to obtain a filtered time-domain waveform:
the time-frequency distribution after the filter of the Catte model can be obtained by the formula (7) and is recorded as Ts′(ωlAnd b), adopting a synchronous extrusion reconstruction formula as follows:
wherein,
re represents the operation of the real part; and x (b) is a time domain waveform reconstructed after the time-frequency surface filtering, namely the estimation of the effective signal x (t) in the formula (1). Fig. 6 shows a comparison of waveforms processed by the filtering method proposed by the present invention and the TFPF method for two-component noisy wavelets. It can be seen that the method of the present invention is more advantageous in maintaining the effective signal (peaks, troughs) and superior to TFPF in random noise suppression.
And (3) performance test of the filtering method:
firstly, the step needs to start with a filtering experiment for artificially synthesizing a seismic data model, analyze a filtering result and then process actual seismic data.
Firstly, processing a single-channel noise-added analog seismic signal according to the steps of I, II and III, and then processing a noise-added synthetic event seismic record. The filtering process for the synthesized record is also performed according to the processing method of the single-channel signal. In general, Ricker wavelets are used to simulate seismic reflection signals, and the mathematical expression is as follows:
x(t)＝[1-2(πf0t)2]exp[-(πf0t)2](9)
the filtering experiment is carried out under different signal-to-noise ratios, so that the method performance can be tested more comprehensively. For single-channel analog seismic signals and synthetic records containing the same-phase axis, random noises with different intensities can be added, so that the signal-to-noise ratio is taken from high to low.
In the filtering experiment process, random noise is added to make the signal-to-noise ratio of two component Ricker wavelets about 5dB, so that the signal-to-noise ratio of two in-phase axes artificially synthesized and recorded is about-4 dB and-12 dB. FIG. 6(a) - (e) are a comparison of waveforms before and after two component Ricker wavelet filtering and a partial enlarged view thereof; FIGS. 8(a) - (f) show two noise-free in-phase axis synthetic recordings and their noisy recordings (-4dB), as well as the recordings and difference recordings filtered using the inventive method and the TFPF method, respectively; fig. 9(a) - (d) show two in-phase axis synthetic recordings without noise and their noisy recordings (-12dB), as well as the recordings filtered using the inventive method and the TFPF method, respectively.
It should be noted here that although there are some differences in the form between the actual seismic signal and the simulated seismic signal, the filtering method designed by the present invention has good robustness and universality, so that the effect consistent with the experimental conclusion of the simulated seismic signal can be obtained when the actual seismic data is processed. FIGS. 10(b), (c) are the results of processing an actual seismic record (shown in FIG. 10 (a)) using the filtering method and TFPF method, respectively, according to the present invention.
And secondly, performing qualitative and quantitative analysis on the seismic exploration signals and the recorded filtering results.
After the analog seismic signals and the artificially synthesized seismic records are filtered, the filtering results can be qualitatively and quantitatively analyzed due to the fact that noise-free conditions exist in the analog data experiment as reference standards. However, since there is no noise-free condition for reference in the actual seismic data processing, qualitative analysis is usually mainly performed.
The qualitative analysis mainly comprises modes of waveform comparison, spectrum analysis, difference making and the like; the quantitative analysis is mainly to calculate the signal-to-noise ratio and the mean square error before and after filtering. And comparing the waveforms before and after filtering with the ideal waveforms, wherein the closer the waveform is to the ideal condition, the better the filtering effect is. The waveform pairs are shown in fig. 6(a) - (e) and fig. 9(e) - (g). As can be seen from these figures, the waveform filtered by the method of the present invention is closer to a noise-free waveform, indicating that the retention of the effective signal and the reduction of the random noise are better than those of the TFPF method.
The most common method of spectrum analysis is to use fourier transform for spectrum analysis, and mainly observe the reduction of middle and high frequency random noise (flatness of spectrum) and the retention of effective signals at low frequency (retention of spectrum amplitude). For the spectrum pairs such as fig. 7(a) - (c) and fig. 12(a), (b), it can be seen that the signal obtained by filtering with the method of the present invention has a very flat spectrum in the high frequency part, while the signal obtained by TFPF has a jitter phenomenon in the spectrum, which indicates that the former has a stronger suppression capability for random noise. While observing the amplitude preservation of the low frequency portion of the spectrum (as shown in fig. 12 (a)), it can be seen that TFPF remains slightly inferior.
The differencing method is to subtract the filtered record from the original noisy record to obtain a difference record, which is mainly a filtered random noise and a small portion of the lost effective in-phase axis. If the less effective in-phase axis is lost in the difference record, the stronger the filtering method can maintain the effective signal. The difference comparisons are made as in fig. 8(e), (f), and it can be seen that the effective in-phase axis of loss is more pronounced in the TFPF difference recording. FIGS. 9(c) and (d) show the sharpness and continuity of the in-phase axes after the filter, and when the SNR is very low, it can be seen that the in-phase axes in the TFPF-processed records are discontinuous (shown by the rectangular box in the figure), while the in-phase axes in the TFPF-processed records are more clearly continuous.
FIG. 10 is a comparison of the results of actual seismic records processed by two methods, and it is clear that the records processed by the method of the present invention have high definition of effective phase axes, good continuity, and high random noise reduction, while the phase axes in the records processed by TFPF are still not clear enough. Since the actual recording is not too noisy (the random noise energy is weak as seen by the time-frequency distribution shown in fig. 11(a) - (d)), no further processing is required after filtering by the simultaneous extrusion time-frequency domain lattice model. From the comparison of the single-channel signal frequency spectrums shown in fig. 12(a) and (b), it can be seen that the method of the present invention has great advantages in suppressing the high-frequency random noise, while the TFPF method has a relatively limited capability in suppressing the high-frequency random noise, and the obtained signal frequency spectrum is not flat enough up to the high-frequency part, and has a small jitter.
In the quantitative analysis, the improvement degree of the filtered signal-to-noise ratio can be obtained by calculating the signal-to-noise ratio and comparing the signal-to-noise ratio with the original signal-to-noise ratio; and calculating the mean square error and comparing the mean square error with the original mean square error to obtain the reduction degree of the mean square error after filtering. The mathematical formula for calculating the signal-to-noise ratio and the mean square error is as follows (taking two-dimensional seismic records as an example):
in the formula, x (i, j) is a record without noise,
for the record after filtering, N is the number of sampling points, and T is the number of tracks.
Table 1 shows the quantitative analysis of the artificially synthesized records of different signal-to-noise ratios and the filtering results thereof, i.e., the signal-to-noise ratios and the mean square errors before and after filtering are calculated by the formulas (10) and (11).
TABLE 1 Artificial synthesis record of different SNR and evaluation index of filtering result
The method provided by the invention not only makes full use of the high-resolution and high-aggregation characteristics of synchronous extrusion time-frequency transformation, namely, the energy of effective signals can be highly aggregated together, the instantaneous frequency curve is very fine, but also makes use of the capability of the Catte model anisotropic diffusion filtering to consider the retention of the characteristics of the effective signals and the reduction of random noise. The method is easy to implement, high in operation speed and strong in parameter robustness, and is suitable for seismic exploration data processing under different noise intensities. Therefore, the seismic exploration signal enhancement method has good popularization and application prospect and value.
Those skilled in the art will appreciate that the deeper and deeper the oil and gas reservoir resources are buried, the more complex the exploration environment is, the more difficult the exploration becomes, and the more difficult the identification of valid signals in the collected seismic exploration data becomes. In order to provide high-quality interpretation materials for data interpretation work, a filtering method with excellent performance is required to recover and extract effective signals in a data processing link and eliminate various noises, and random noises are distributed in the whole data and are one of the main existing noise types. In view of this, a high-performance filtering method capable of significantly enhancing effective signals and greatly reducing random noise is required to process seismic exploration data, so as to lay a solid foundation for the subsequent data interpretation link.
Claims (1)
1. A method of enhancing high performance time-frequency domain filtering of seismic survey signals, comprising the steps of:
step one, performing synchronous extrusion wavelet transformation on a noise-containing seismic signal to obtain high-resolution time-frequency distribution:
setting the noise-containing signal as:
s(t)＝x(t)+n(t) (1)
x (t) is the effective signal, n (t) is additive random noise;
and (3) performing continuous wavelet transform on s (t), and referring to the formula (2), wherein the continuous wavelet transform comprises the following steps:
where a is the scale factor, b is the translation factor, ψ (t) is the wavelet basis function,
a conjugate function of ψ (t); according to the paseuler theorem, the frequency domain expression of equation (2) is:
where ξ is the angular frequency,
is the Fourier transform of s (t), phi (a ξ) is phi (t), and the time will be obtained by the formulaThe convolution form of the inter-domain is transformed into the product form of the frequency domain, and the instantaneous frequency can be estimated initially by derivation of the wavelet coefficients, i.e.
Where i is the imaginary unit; this converts the time-scale plane (b, a) into a time-frequency plane (b, ω)s(a, b)), both the frequency ω and the scale a are discretized, thus only at discrete points ai{ai-ai-1＝(Δa)iCalculating wavelet coefficient Ws(a, b); the implementation process of frequency rearrangement in synchronous extrusion wavelet transformation comprises the following steps: assuming that the sampling frequency of the signal is cf and the number of scales taken during wavelet transform is N, the discrete frequency value on the rearranged time-frequency spectrum is ωl＝l*cf/N,l∈[1,N](ii) a Synchronous extrusion value T of wavelet coefficients(ωlB) by extrusion of any center frequency ωlNearby section
Is obtained by a value of (i.e. a
Wherein, akIs a discrete scale, and Δ ak＝ak-ak-1，Ts(ωlB) is the high-resolution and high-aggregation time-frequency distribution after extrusion rearrangement;
step two, performing Catte model anisotropic diffusion filtering on the high-resolution time-frequency plane, wherein the Catte model is as follows:
where U (x, y, t) represents the pixel value at time (x, y) t, U0(x, y) represents an original image,
is a divergence operator, and is a function of the divergence,
is a gradient operator; gσThe method is characterized in that the method is a Gaussian kernel function with a scale factor of sigma, C is a diffusion coefficient, and a non-negative monotonically decreasing function is usually selected; perona and Malik propose two diffusion coefficient functions:
and
in the formula, k is a constant, and a proper integer can be selected according to a processed object;
is the gradient mode of the constant coefficient thermal diffusion equation with the scale factor of sigma;
the time frequency distribution T obtained by the formula (5) iss(ωlAnd b) regarding the image as an image, performing anisotropic diffusion filtering on the image, and substituting the image into the formula (6) to obtain:
thus, the filtered time-frequency plane is obtained, the energy gathering point of the effective signal is obviously enhanced, the random noise is suppressed to a great extent, and particularly the high-frequency random noise is greatly reduced; if the random noise in the noise-containing signal is strong, a part of noise at the intermediate frequency still remains on the filtered time-frequency plane, so that in order to obtain a better noise elimination effect, time-frequency domain band-pass filtering is needed;
step three, performing synchronous extrusion wavelet reconstruction on the filtered time-frequency distribution to obtain a filtered time-domain waveform:
the time-frequency distribution after the filter of the Catte model can be obtained by the formula (7) and is recorded as Ts′(ωlAnd b), adopting a synchronous extrusion reconstruction formula as follows:
wherein,
re represents the operation of the real part; and x (b) is a time domain waveform reconstructed after the time-frequency surface filtering, namely the estimation of the effective signal x (t) in the formula (1).
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010130916.4A CN111175827B (en) | 2020-02-28 | 2020-02-28 | A High-Performance Time-Frequency Domain Filtering Method for Enhanced Seismic Exploration Signals |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010130916.4A CN111175827B (en) | 2020-02-28 | 2020-02-28 | A High-Performance Time-Frequency Domain Filtering Method for Enhanced Seismic Exploration Signals |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111175827A true CN111175827A (en) | 2020-05-19 |
| CN111175827B CN111175827B (en) | 2022-11-01 |
Family
ID=70653300
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010130916.4A Active CN111175827B (en) | 2020-02-28 | 2020-02-28 | A High-Performance Time-Frequency Domain Filtering Method for Enhanced Seismic Exploration Signals |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111175827B (en) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115079254A (en) * | 2022-05-19 | 2022-09-20 | 西安石油大学 | VSP up-down traveling wave separation method based on Radon domain AWT |
| CN116593831A (en) * | 2023-07-19 | 2023-08-15 | 西安交通大学 | A cable defect location method, equipment and medium |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090292475A1 (en) * | 2007-05-25 | 2009-11-26 | Aftab Alam | Time-Space Varying Spectra for Seismic Processing |
| US20100166330A1 (en) * | 2008-11-07 | 2010-07-01 | Thyagarajan Kadayam S | Systems and Methods of Using Spatial/Spectral/Temporal Imaging for Hidden or Buried Explosive Detection |
| CN103926616A (en) * | 2014-04-11 | 2014-07-16 | 中国海洋石油总公司 | Multi-scale anisotropic diffusion filtering method based on pre-stack CRP trace sets |
| CN106772588A (en) * | 2017-02-24 | 2017-05-31 | 西南石油大学 | A kind of frequency domain self-adaptation nonlinear seismic imaging filtering method |
| CN107991707A (en) * | 2017-12-05 | 2018-05-04 | 吉林大学 | A kind of borehole microseismic first break picking method based on kurtosis characteristic in shear let domains |
| CN109669213A (en) * | 2019-02-25 | 2019-04-23 | 中国石油化工股份有限公司 | Frequency dividing diffusing filter tomography intensifying method based on optimization Morlet small echo |
-
2020
- 2020-02-28 CN CN202010130916.4A patent/CN111175827B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090292475A1 (en) * | 2007-05-25 | 2009-11-26 | Aftab Alam | Time-Space Varying Spectra for Seismic Processing |
| US20100166330A1 (en) * | 2008-11-07 | 2010-07-01 | Thyagarajan Kadayam S | Systems and Methods of Using Spatial/Spectral/Temporal Imaging for Hidden or Buried Explosive Detection |
| CN103926616A (en) * | 2014-04-11 | 2014-07-16 | 中国海洋石油总公司 | Multi-scale anisotropic diffusion filtering method based on pre-stack CRP trace sets |
| CN106772588A (en) * | 2017-02-24 | 2017-05-31 | 西南石油大学 | A kind of frequency domain self-adaptation nonlinear seismic imaging filtering method |
| CN107991707A (en) * | 2017-12-05 | 2018-05-04 | 吉林大学 | A kind of borehole microseismic first break picking method based on kurtosis characteristic in shear let domains |
| CN109669213A (en) * | 2019-02-25 | 2019-04-23 | 中国石油化工股份有限公司 | Frequency dividing diffusing filter tomography intensifying method based on optimization Morlet small echo |
Non-Patent Citations (5)
| Title |
|---|
| YAN-PING LIU，等: "THE ENHANCEMENT OF SEISMIC SIGNAL BY THE FILTERING METHOD BASED ON SYNCHROSQUEEZING TRANSFORM", 《IEEE》 * |
| 刘晗等: "利用同步挤压变换检测微地震信号", 《中国科技论文》 * |
| 张建中，等: "基于同步挤压变换的水合物储层地震信号时频分析", 《海洋地质前沿》 * |
| 张琪，等: "基于同步挤压变换的地震勘探信号时频特性分析", 《勘探开发》 * |
| 闫安菊，等: "改进的各向异性扩散滤波方法压制地震数据噪声", 《断块油气田》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115079254A (en) * | 2022-05-19 | 2022-09-20 | 西安石油大学 | VSP up-down traveling wave separation method based on Radon domain AWT |
| CN116593831A (en) * | 2023-07-19 | 2023-08-15 | 西安交通大学 | A cable defect location method, equipment and medium |
| CN116593831B (en) * | 2023-07-19 | 2023-11-07 | 西安交通大学 | Cable defect positioning method, device and medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111175827B (en) | 2022-11-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gómez et al. | A simple method inspired by empirical mode decomposition for denoising seismic data | |
| US11079418B2 (en) | Methods for extending frequency transforms to resolve features in the spatio-temporal domain | |
| CN110174702B (en) | Method and system for recovering low-frequency weak signals of marine seismic data | |
| CN106405645B (en) | Frequency processing method is opened up in a kind of controllable earthquake of signal-to-noise ratio based on data quality analysis | |
| CN103995289A (en) | Time-varying mixed-phase seismic wavelet extraction method based on time-frequency spectrum simulation | |
| CN111505716A (en) | Seismic time-frequency analysis method for extracting generalized Chirplet transform based on time synchronization | |
| Czarnecki | The instantaneous frequency rate spectrogram | |
| Ma et al. | Noise reduction for desert seismic data using spectral kurtosis adaptive bandpass filter | |
| CN102183787A (en) | Method for improving seismic data resolution based on seismographic record varitron wave model | |
| CN106680876B (en) | A kind of seismic data joint denoising method | |
| CN107132579A (en) | A kind of attenuation of seismic wave compensation method for protecting earth formation | |
| CN107179550A (en) | A kind of seismic signal zero phase deconvolution method of data-driven | |
| CN111175827A (en) | A High-Performance Time-Frequency Domain Filtering Method for Enhanced Seismic Exploration Signals | |
| CN110531420A (en) | A Non-destructive Separation Method of Industrial Interference Noise in Seismic Data | |
| CN102096101A (en) | Method and device for extracting hybrid-phase seismic wavelets | |
| CN113608259B (en) | Seismic thin layer detection method based on ICEEMDAN constraint generalized S transformation | |
| CN106125134A (en) | Based on the geological data signal-noise ratio computation method of window during hyperbolic | |
| Liu et al. | Random noise reduction using SVD in the frequency domain | |
| CN108646296A (en) | Desert seismic signal noise reduction methods based on Adaptive spectra kurtosis filter | |
| CN106772588A (en) | A kind of frequency domain self-adaptation nonlinear seismic imaging filtering method | |
| Oropeza et al. | Multifrequency singular spectrum analysis | |
| CN112925013B (en) | Seismic data high-resolution processing method based on full-band continuation fidelity | |
| CN109212608A (en) | Borehole microseismic signal antinoise method based on 3D shearlet transformation | |
| Zhang et al. | Seismic random noise attenuation and signal-preserving by multiple directional time-frequency peak filtering | |
| CN114791628A (en) | Seismic result data absorption compensation quantitative analysis and evaluation method |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |