Earthquake surface wave suppression method based on EWT-Curvelet

Abstract

The invention discloses a seismic surface wave suppression method based on EWT-Curvelet, and belongs to the technical field of seismic exploration data processing. The method comprises the following steps: firstly, performing empirical wavelet transform (EWT) on post-stack seismic data channel by channel, extracting a frequency localized modal component through adaptive spectrum segmentation, and constructing a surface wave noise model; then converting the model to a curvelet domain, performing physical constraint filtering at a specific scale and in a nearly horizontal direction, and combining pixel connectivity analysis to suppress a continuous noise region; and finally, reconstructing and optimizing noise through inverse transformation and removing the noise from the original data. The method has the advantages that the modal aliasing problem is solved through channel-by-channel adaptive spectrum segmentation, a surface wave energy dense region is accurately locked by curvelet domain scale-direction combined filtering, isolated effective signals are protected in combination with connectivity analysis, high-precision suppression of surface waves is efficiently achieved in a non-iterative mode, meanwhile, the structural integrity of reflected waves is kept, and high-precision suppression of the surface waves is achieved. The seismic section signal-to-noise ratio and the geological interpretation precision are remarkably improved, and the method is suitable for land seismic exploration data processing.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration data processing, and in particular to a seismic surface wave suppression method based on empirical wavelet transform (EWT) and curvelet transform. Background Art

[0002] In the field of seismic exploration, interference from surface waves (also known as ground roll) has always been a core problem restricting data quality. This type of noise typically manifests as low-frequency (2-25Hz), high-amplitude coherent energy that propagates along the shallow surface with strong horizontal consistency. Its energy intensity can be more than 10 times that of the effective reflected wave, forming a large-scale "strip-like" contamination zone in the time and space domain, seriously obscuring the effective reflection signals in the middle and deep layers. Especially in areas with complex shallow geological structures (such as deserts and piedmont zones), the frequency and propagation direction of surface waves and reflected waves overlap significantly, making it difficult to achieve effective separation using traditional filtering methods.

[0003] Early surface roll suppression methods primarily relied on signal spectra or apparent velocity differences. For example, high-pass filtering directly removes low-frequency components by setting a fixed cutoff frequency (e.g., 15 Hz). However, this indiscriminately eliminates all signals within this frequency band, permanently losing weak shallow reflections. Frequency-wavenumber (fk) domain filtering attempts to suppress surface rolls by exploiting their elliptical distribution in the wavenumber domain. However, when subsurface structures are tilted, the reflected wave events in the fk domain also exhibit low-frequency, high-wavenumber characteristics, and their overlap with the surface roll distribution can exceed 40%, resulting in a significant number of valid signals being misclassified as noise. While the Radon transform can separate linear noise using time differences, it lacks adaptability to surface rolls with non-hyperbolic morphologies, and its separation accuracy plummets in energy-ambiguous regions of the velocity spectrum (e.g., near-offset).

[0004] To improve adaptability, multi-channel signal decomposition technology has been introduced into this field in recent years. Multi-channel variational mode decomposition (MVMD) decomposes the seismic profile into bandwidth-limited modal components through global optimization, concentrating the surface wave energy on the first few low-frequency modes. However, this method is extremely sensitive to the selection of initial values ​​and regularization parameters: when the initial value deviates from the optimal solution by 10%, the modal aliasing rate can reach 15%-20%, causing part of the reflected wave energy to leak into the surface wave mode. In addition, the "inter-channel stationarity" assumption implicit in MVMD fails in non-uniform geological areas such as faults and pinch-outs. Measured data show that when the waveform difference between adjacent channels exceeds 30%, the surface wave residual rate rises to more than 25%.

[0005] Deep learning-driven implicit neural representation (INR) offers a new approach to noise modeling. This method utilizes a multi-layer perceptron (MLP) to map spatial coordinates to signal values, reconstructing clean signals from noisy data through unsupervised learning. Although INR performs well in suppressing random noise, it has inherent flaws in modeling the physical properties of surface waves. First, the network struggles to learn the characteristics of noise with clear directionality, resulting in an error of over 25% in directional suppression of surface waves. Second, the training process relies on hyperparameter tuning, requiring more than 5,000 iterations (taking >30 minutes) for a single profile, which cannot meet the real-time processing requirements of exploration sites. Third, the model suffers from poor interpretability and cannot provide the physical basis for noise separation required by geologists.

[0006] The core contradiction facing current technologies is that physical-driven methods (such as MVMD) have difficulty adapting to the strong non-stationarity of complex geological scenarios, while data-driven methods (such as INR) lack the ability to accurately model the directionality and spatial continuity of surface rolls. This contradiction is particularly acute in the spectral overlap region (5-15Hz), and existing methods cannot simultaneously meet the industrial requirements of "high fidelity, strong adaptability, and high efficiency." Therefore, it is urgent to develop an innovative solution that integrates adaptive signal decomposition and physically constrained directional filtering to fundamentally solve the problems of modal aliasing, signal damage, and computational efficiency in surface roll suppression. Summary of the Invention

[0007] In view of the defects of the prior art, the present invention provides a seismic surface wave suppression method based on EWT-Curvelet.

[0008] In order to achieve the above object of the invention, the technical solution adopted by the present invention is as follows:

[0009] A seismic surface wave suppression method based on EWT-Curvelet, characterized by comprising the following steps:

[0010] Step 1, adaptive EWT decomposition on a per-trace basis: Perform empirical wavelet transform (EWT) on the post-stack seismic data, adaptively determine the decomposition scale through spectral analysis, and extract multiple frequency-localized intrinsic mode functions (IMFs);

[0011] Step 2, surface wave model construction: select the IMF containing significant surface wave components from the IMF obtained in step 1 to construct the surface wave noise model;

[0012] Step 3, Curvelet Domain Directional Filtering: The surface roll noise model is zero-filled and then converted to the curvelet domain. Multi-scale analysis is performed at scales containing surface roll components and in the near-horizontal direction of the curvelet transform. Combined with pixel connectivity screening, regions of surface roll energy with high spatial continuity are suppressed.

[0013] Step 4: Signal reconstruction and noise removal: Perform inverse transform on the filtered curvelet coefficients to reconstruct the spatiotemporal noise model, subtract the noise model from the original seismic data, and obtain the output data of surface roll suppression.

[0014] Furthermore, the EWT decomposition in step 1 includes:

[0015] Perform Gaussian filtering and smoothing on the seismic signal spectrum;

[0016] Determine the cutoff point based on the smoothed spectrum and construct an empirical wavelet filter bank;

[0017] The IMF components of different frequency bands in the signal are extracted through the filter bank.

[0018] Furthermore, the basis for selecting the low-order IMF in step 2 is:

[0019] The spectrum energy of IMF is concentrated in the low frequency band (0-15Hz);

[0020] IMF presents horizontally consistent waveform characteristics in the time and space domain.

[0021] Furthermore, the pixel connectivity screening in step 3 is specifically as follows:

[0022] Identify 8-connected regions in the curvelet coefficients;

[0023] Coefficient suppression is performed only on continuous regions with an area greater than 200 pixels.

[0024] Furthermore, the near-horizontal direction in step 3 is defined as a coefficient subset with direction indices of 0°±15° and 180°±15° in the curvelet domain.

[0025] Furthermore, the EWT in step 1 can be replaced by any of the following signal decomposition methods:

[0026] Empirical mode decomposition (EMD); variational mode decomposition (VMD); time-varying filtering (TVF).

[0027] Furthermore, the curvelet transform in step 3 can be replaced by shearlet transform or wavelet packet transform.

[0028] Compared with the prior art, the advantages of the present invention are:

[0029] 1. Through a channel-by-channel adaptive spectrum segmentation mechanism, EWT dynamically constructs a filter bank that matches the signal's spectral characteristics, completely avoiding the modal aliasing problem caused by a fixed number of decomposition layers. This ensures that the surface roll energy is strictly constrained within physically clear low-frequency modes, providing a pure input source for subsequent noise modeling.

[0030] 2. The curvelet transform performs directional filtering within the physically constrained space of coarse-scale and near-horizontal directions. Combined with pixel-level connectivity analysis, it intelligently distinguishes continuous noise from isolated reflection points, preserving the edge details of complex geological structures such as faults and river channels while eliminating the risk of inadvertent damage to tilted events encountered with traditional methods.

[0031] 3. EWT's non-iterative frequency band segmentation eliminates the burden of parameter tuning and achieves an order of magnitude improvement in computational speed compared to variational mode decomposition (VMD). At the same time, the trace-by-trace processing strategy endows the method with natural robustness to drastic inter-trace variations, maintaining stable performance in diverse geological scenarios from deserts to mountains. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 10 IMF component maps obtained by decomposing the synthetic seismic traces through EWT in an embodiment of the present invention;

[0033] Figure 2 is a third-order IMF diagram containing surface roll components according to an embodiment of the present invention;

[0034] Figure 3 : is an initial surface wave model diagram constructed by superposition of third-order IMFs in an embodiment of the present invention;

[0035] Figure 4 8 scale component images after curvelet transformation of the surface wave model according to the embodiment of the present invention;

[0036] Figure 5 This is a 6-8 scale surface roll suppression image containing surface roll components according to an embodiment of the present invention;

[0037] Figure 6 This is a graph of synthetic data after surface roll suppression according to an embodiment of the present invention;

[0038] Figure 7 It is a flow chart of an embodiment of the present invention. DETAILED DESCRIPTION

[0039] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and examples.

[0040] like Figure 7 As shown, the present invention provides a seismic surface wave suppression method based on EWT-Curvelet, comprising the following steps:

[0041] (1) Channel-by-channel adaptive EWT decomposition

[0042] Execute the following for each signal of the stacked seismic data:

[0043] The spectrum is Gaussian smoothed (standard deviation σ = 0.05 × sampling rate) to eliminate spectrum burrs;

[0044] Detect local maximum points and divide the frequency bands based on the minimum points ( Figure 1 );

[0045] An empirical wavelet filter bank is constructed to extract the frequency-localized IMF components.

[0046] (2) Construction of surface wave model

[0047] Selection criteria: Select low-frequency IMFs (usually the first three orders) whose energy is concentrated in the range of 0-15Hz, and whose spatiotemporal waveforms satisfy horizontal continuity ( Figure 2 );

[0048] Model generation: superimpose selected IMFs to construct the initial surface wave model ( Figure 3 ).

[0049] (3) Curvelet domain directional filtering

[0050] Zero padding: fill the boundary of the surface wave model with zeros to an integer power of 2 to avoid boundary effects;

[0051] Curvelet transform: decomposed into 8 scales ( Figure 4 ), each scale contains 24 directional sub-bands;

[0052] Coefficient screening:

[0053] Target scale: 6-8 scale (corresponding to coarse-scale features with wavelengths > 100 m);

[0054] Target direction: 0°±15°, 180°±15° (horizontal propagation direction);

[0055] Connectivity analysis: only suppress the coefficients of 8 connected regions with an area greater than 200 pixels ( Figure 5 ).

[0056] (4) Signal reconstruction and noise removal

[0057] Perform inverse curvelet transform on the filtered coefficients to reconstruct the optimized surface wave model;

[0058] Subtract the model from the original data and output the surface wave suppression result ( Figure 6 ).

[0059] Example 1: Synthetic Data Testing

[0060] Data parameters:

[0061] Synthetic record: 256 channels × 512 sampling points, main frequency 30Hz reflection wave + 10Hz surface wave;

[0062] Surface wave energy: The amplitude is three times that of the reflected wave, and the frequency overlap band is 5-15Hz.

[0063] Implementation steps:

[0064] (1) EWT decomposition: 12 IMFs are extracted per channel, with a Gaussian filter of σ = 2 Hz (sampling rate 40 Hz);

[0065] (2) Select IMF1-3 to construct the surface wave model (spectral peak 8Hz);

[0066] (3) Curvelet transform: scale 6-8 (wavelength 120-200m), direction angle 175°-185°;

[0067] (4) Suppress connected regions with an area greater than 200 pixels (accounting for 85% of the target coefficient).

[0068] result:

[0069] Surface wave energy attenuation is 96.3%, and reflected wave amplitude loss is less than 2%;

[0070] PSNR=23.1dB, 42 seconds (CPU: i7-11800H).

[0071] Example 2: Measured land earthquake data

[0072] Data source: Post-stack profile of an oil field (1024 traces × 2048 points);

[0073] Parameter adjustment:

[0074] To address the strong inter-channel variations, median filtering preprocessing is added to the EWT cutoff point detection;

[0075] The connected area threshold was raised to 300 pixels (to avoid suppressing tilted events).

[0076] Effect:

[0077] After surface wave suppression, the clarity of shallow fault imaging is improved by 40% (assessed by geological experts);

[0078] Processing speed: 3.2 minutes / section (vs. 28 minutes for INR).

[0079] Alternatives

[0080] EWT alternative: VMD decomposition (penalty coefficient α = 2000, mode number K = 8), but the PSNR dropped to 20.8dB;

[0081] Curvelet replacement: Shearlet transform (scale 5, direction 16), direction error causes the reflection wave loss to increase by 5%.

[0082] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the implementation methods of the present invention, and it should be understood that the scope of protection of the present invention is not limited to such specific descriptions and embodiments. Those skilled in the art can make various other specific variations and combinations based on the technical teachings disclosed in the present invention without departing from the essence of the present invention, and such variations and combinations are still within the scope of protection of the present invention.

Claims

1. A seismic surface wave suppression method based on EWT-Curvelet, characterized in that: The following steps are involved: Step 1, adaptive EWT decomposition on a per-trace basis: perform empirical wavelet transform (EWT) on the post-stack seismic data, adaptively determine the decomposition scale through spectral analysis, and extract multiple frequency-localized intrinsic mode functions (IMFs); Step 2: Surface wave model construction: Select low-order IMFs containing significant surface wave components from the IMFs obtained in step 1 to construct a surface wave noise model; Step 3, Curvelet Domain Directional Filtering: The surface roll noise model is zero-filled and then converted to the curvelet domain. Multi-scale analysis is performed at the high-frequency scale and near-horizontal direction containing the surface roll components of the curvelet transform. Combined with pixel connectivity screening, the surface roll energy areas with high spatial continuity are suppressed. Step 4: Signal reconstruction and noise removal: Perform inverse transform on the filtered curvelet coefficients to reconstruct the spatiotemporal noise model, subtract the noise model from the original seismic data, and obtain the output data of surface roll suppression.

2. The seismic surface wave suppression method according to claim 1, characterized in that: The EWT decomposition described in step 1 includes: Perform Gaussian filtering and smoothing on the seismic signal spectrum; Determine the cutoff point based on the smoothed spectrum and construct an empirical wavelet filter bank; The IMF components of different frequency bands in the signal are extracted through the filter bank.

3. The seismic surface wave suppression method according to claim 1, characterized in that: The basis for selecting the low-order IMF in step 2 is: The spectrum energy of IMF is concentrated in the low frequency band; IMF presents horizontally consistent waveform characteristics in the time and space domain.

4. The seismic surface wave suppression method according to claim 1, wherein: The pixel connectivity screening in step 3 is specifically as follows: Identify 8-connected regions in the curvelet coefficients; Coefficient suppression is performed only on continuous regions with an area greater than 200 pixels.

5. The seismic surface wave suppression method according to claim 1, characterized in that: The near-horizontal direction in step 3 is defined as a coefficient subset with direction indices of 0°±15° and 180°±15° in the curvelet domain.

6. The seismic surface wave suppression method according to claim 1, characterized in that: The EWT in step 1 can be replaced by any of the following signal decomposition methods: Empirical mode decomposition (EMD); variational mode decomposition (VMD); time-varying filtering (TVF).

7. The seismic surface wave suppression method according to claim 1, characterized in that: The curvelet transform in step 3 can be replaced by shearlet transform or wavelet transform.

Images