Seismic data denoising method based on Contourlet transformation

Abstract

The invention provides a seismic data denoising method based on Contourlet transformation, and the method is characterized in that the method comprises the following steps: 1, carrying out the loading and preprocessing of seismic data; step 2, denoising is carried out by adopting multi-scale cyclic translation Contourlet transformation; step 3, adopting improved non-local mean filtering; and 4, processing and outputting data. According to the method, a cyclic translation algorithm and non-local mean filtering (FNLM) are combined, so that the denoising capability and weak signal identification capability of seismic signals are effectively improved.

Description

Technical Field

[0001] The present invention provides a seismic data denoising method based on Contourlet transform, belonging to the technical field of seismic exploration. Background Technique

[0002] In the fields of seismic exploration, image processing, medical image analysis, etc., signals are usually interfered by various noises, resulting in a decline in data quality, which in turn affects subsequent analysis and interpretation. How to effectively remove noise while maximizing the retention of the detailed features of the original signal, especially the enhancement of weak signals, is the core issue in this field.

[0003] Currently, the methods for removing random noise from seismic data are mainly divided into space-time domain methods and transform domain methods:

[0004] 1. Space-time domain filtering method:

[0005] It mainly includes median filtering, mean filtering, Fourier transform filtering, etc., and denoises relying on the local statistical characteristics of the signal. These methods are simple to calculate, but it is difficult to distinguish effective signals from random noise, and signal loss is likely to occur.

[0006] 2. Transform domain denoising method:

[0007] Representative methods include wavelet transform and empirical mode decomposition (EMD).

[0008] These methods use the different scale characteristics of the signal for denoising and can better preserve the signal details, but there are still limitations in terms of directionality and high-resolution processing.

[0009] Denoising method based on wavelet transform: Wavelet transform is a multi-scale analysis tool and is widely used in signal denoising. The disadvantages are:

[0010] Insufficient direction information: Wavelet transform mainly decomposes the signal along the horizontal and vertical directions and lacks the ability to express oblique features.

[0011] Pseudo-Gibbs phenomenon: Unnatural fluctuations may be introduced at signal mutation points.

[0012] Denoising method based on empirical mode decomposition (EMD). It is a time-frequency analysis method suitable for non-stationary signals and is mainly used to remove low-frequency trends and high-frequency noises. The disadvantages are:

[0013] Mode mixing problem: Some IMF components may contain both effective signals and noises at the same time, making it difficult to separate.

[0014] High computational complexity: EMD requires iterative calculation, with a high computational cost, making it difficult to apply to large-scale data processing.

[0015] The denoising effect is affected by signal characteristics: For some types of seismic signals, EMD may not be able to effectively remove noise. Summary of the Invention

[0016] The object of the present invention is to optimize the seismic data denoising performance, improve the signal fidelity and weak signal recognition ability on the basis of the Contourlet transform, in combination with the cyclic translation algorithm and non-local means filtering (FNLM). This scheme can enhance seismic event axes, boundary features, and weak signals while suppressing random noise, achieving high-quality signal recovery.

[0017] The specific technical solution is as follows:

[0018] A method for denoising seismic data based on the Contourlet transform, comprising the following steps:

[0019] Step 1: Loading and preprocessing of seismic data;

[0020] (1) Data parsing and initialization

[0021] Input parameters: filename (path of the SEGY file with noise), sigma_contourlet (Contourlet threshold parameter), sigma_fnlm (non-local means smoothing coefficient);

[0022] Parse the SEGY file through the ReadSegy function to obtain the seismic data matrix seis and its header information;

[0023] Dynamically construct the processing matrix: Automatically calculate the optimal processing size n = 2^nextpow2(max(high, longth)) according to the input data size [high, longth].

[0024] (2) Data central embedding

[0025] Adopt a symmetric padding strategy:

[0026] start_row = floor((n - high) / 2) + 1;

[0027] start_col = floor((n - longth) / 2) + 1;

[0028] Data_single(start_row:start_row+high-1,start_col:start_col+longth-1) = seis;

[0029] Step 2: Use multi-scale circular shift Contourlet transform for denoising;

[0030] (1) Parameter system construction

[0031] Decomposition structure: pfilt = 'pkva' (dual-tree complex wavelet filter) / dfilt = 'pkva' (directional filter bank);

[0032] Decomposition levels: nlevs = [4, 4, 3] (the number of pyramid decomposition levels decreases gradually);

[0033] Noise modeling: rho = 3 (noise intensity coefficient), sigma = sig / rho (dynamic noise standard deviation).

[0034] (2) Circular shift optimization

[0035] Perform 4 translation iterations (num_translations = 4, the number of loops can be adjusted: 2 - 8 (even)), and the translation amount t ∈ [1, 4] each time

[0036] Core processing flow:

[0037] Translate data → Contourlet decomposition → Hard threshold function processing → Inverse transform reconstruction → Translation anti-correction

[0038] The threshold calculation model uses hard threshold denoising:

[0039] cth = sigma_contourlet * sigma * sqrt(nvar)

[0040] where nvar is the Contourlet sub-band noise variance

[0041] Step 3: Use improved non-local means filtering;

[0042] (1) Data normalization processing

[0043] Dynamic range compression:

[0044] normalized_data = (cumulative_result - data_min) / (data_max - data_min)

[0045] Retain the data distribution characteristics and avoid gradient disappearance

[0046] (2) FNLM Parameter Configuration

[0047] Smoothness control: sigma_fnlm;

[0048] Similarity measurement window: default is 8×8 search domain;

[0049] Weighted coefficient calculation:

[0050] w_ij = exp(-||I(i)-I(j)||2 / (h2))

[0051] where h = sigma_fnlm * sqrt(max(high, longth)).

[0052] Step Four: Data Processing and Output;

[0053] (1) Spatial Domain Restoration

[0054] Size restoration: Cropped to the original data size;

[0055] Edge compensation: Weighted average processing for the filled area;

[0056] (2) SEGY Format Output

[0057] File naming rule: [Prefix]_[File Name]_contourlet.segy

[0058] Metadata inheritance: Completely retain the original SegyTraceHeaders information.

[0059] The specific algorithm in the above method is:

[0060] 1. Contourlet transform for multi-scale and multi-directional decomposition

[0061] (1) Laplacian pyramid (LP) for multi-scale decomposition

[0062] The Laplacian pyramid is used to perform multi-scale decomposition on the input signal f(x,y). The decomposition steps are as follows:

[0063] a. Gaussian pyramid (Gaussian Pyramid)

[0064] Low-frequency component L j (x,y) can be achieved through the Gaussian blur filter G(x,y) and downsampling:

[0065] G j (x,y) = (G * G j-1 )(x,y)↓2

[0066] where, G * G j-1Indicates a two-dimensional convolution, and ↓2 indicates downsampling.

[0067] b. Laplacian pyramid to calculate the high-frequency residual

[0068] L j (x, y) = G j-1 (x, y) - Upsample(G j (x, y))

[0069] where Upsample represents interpolation magnification to restore to the scale of G j-1 .

[0070] c. Low-frequency part G J (x, y) is used as the final low-frequency information, and the high-frequency part L j (x, y) contains edge and noise information.

[0071] (2) Directional filter bank (DFB) for direction decomposition

[0072] Apply the directional filter bank (DFB) to the high-frequency part L j (x, y) obtained from the LP decomposition to enhance the direction information:

[0073] a. Directional filtering operation

[0074] Let L j (x, y) be the high-frequency residual. After the DFB transformation, each direction component can be expressed as:

[0075]

[0076] where H k (x, y) is the directional filter corresponding to k directions (such as horizontal, vertical, diagonal, etc.).

[0077] b. Signal reconstruction

[0078]

[0079] where k is the number of direction components, usually k = 8, 16 to provide more refined direction information.

[0080] 2. Circular translation algorithm

[0081] In the threshold denoising process, if the transform lacks translational invariance, pseudo-Gibbs phenomena will occur at positions adjacent to the discontinuity points of the effective signal, resulting in visual distortion. This is especially obvious for the edge details in the image. The circular translation algorithm is proposed to suppress this phenomenon generated during the wavelet transform threshold denoising process.

[0082] The Contourlet transform itself also has this property. During the attenuation of random noise in seismic data, "circular translation - threshold denoising - inverse circular translation" is performed in the row and column directions. In this way, threshold denoising of the signal after each translation will cause the pseudo - Gibbs phenomenon to appear in different places, and different denoising results will be obtained for each sub - translation process in the row and column directions:

[0083]

[0084] Linear accumulation of all denoising results will obtain a denoising result that suppresses the pseudo - Gibbs phenomenon:

[0085]

[0086] Among them, K1 and K2 respectively represent the maximum translation amounts in the row and column directions, S is the circular translation operator, i, j and -i, -j are the translation amounts in the row and column directions respectively, CT is the Contourlet transform operator, CT -1 is the inverse Contourlet transform operator, and ∧ is the threshold operator.

[0087] 3. Fast Non - Local Means Filtering (FNLM)

[0088] It is used to remove residual noise and enhance signal details simultaneously.

[0089] (1) Similar block search

[0090] a. Calculate the Euclidean distance between two signal blocks P(x, y) and P(x′, y′):

[0091]

[0092] Among them, ω(i, j) is the weight.

[0093] b. Calculate the weight:

[0094]

[0095] Among them, h is the filtering parameter.

[0096] (2) Calculate the denoised signal

[0097] Optimized calculation of the mean value using the integral image:

[0098] I(x, y) = I(x - 1, y)+I(x, y - 1)-I(x - 1, y - 1)+S(x, y)

[0099] Among them, I(x, y) is the integral image, which reduces the computational complexity.

[0100] As a more advanced transformation method, the Contourlet transform combines the Laplacian pyramid (LP) and the directional filter bank (DFB), providing richer directional information than the traditional wavelet transform and being able to more accurately describe the directional characteristics of complex signals. The present invention combines the cyclic translation algorithm, non-local means filtering (FNLM), and directional smoothing diffusion (DSD) to effectively improve the denoising ability and weak signal recognition ability of seismic signals.

[0101] The technical effects of the technical solution of the present invention are as follows:

[0102] 1. Improve the ability to express directional information:

[0103] Compared with the wavelet transform, the Contourlet transform can more finely represent complex directional features and improve the signal recognition ability.

[0104] 2. Reduce the pseudo-Gibbs phenomenon:

[0105] Combined with the cyclic translation algorithm, through the superposition of multiple translations - transforms - inverse transforms, the denoising robustness is improved, the sensitivity of the Contourlet transform to translation is eliminated, the translation invariance of the transform is improved, the edge artifacts are reduced, and the pseudo-Gibbs effect is reduced.

[0106] 3. Enhance the ability to retain weak signals:

[0107] Through non-local means filtering, the loss of useful signals is avoided and the signal quality is improved.

[0108] 4. Optimize the computational efficiency:

[0109] Adopt the fast FNLM algorithm to improve the processing efficiency of large-scale data and reduce the computational cost.

[0110] 5. Improve the denoising stability:

[0111] The combination of the cyclic translation algorithm and the fast FNLM algorithm further optimizes the signal quality, enabling it to achieve good results in different application scenarios. Brief description of the drawings

[0112] Figure 1 is the technical flow chart of the present invention;

[0113] Figure 2 is the original Marmousi model of the embodiment;

[0114] Figure 3a is the original noisy seismic data of the embodiment;

[0115] Figure 3b is the symmetrically padded noisy seismic data of the embodiment;

[0116] Figure 4a Effect diagram after the first cycle translation of the embodiment;

[0117] Figure 4b Effect diagram after the second cycle translation of the embodiment;

[0118] Figure 4c Effect diagram after the third cycle translation of the embodiment;

[0119] Figure 4d Effect diagram after the fourth cycle translation of the embodiment;

[0120] Figure 4e Total effect diagram after the cycle translation algorithm of the embodiment:

[0121] Figure 5 Effect diagram after the embodiment performs fast non-local means filtering (FNLM);

[0122] Figure 6a Profile diagram after adding random noise to the embodiment;

[0123] Figure 6b Profile diagram after denoising by Contourlet transform;

[0124] Figure 6c Profile diagram after denoising by Contourlet transform + cycle translation algorithm;

[0125] Figure 6d Profile diagram after denoising by the technology of the present invention in the embodiment;

[0126] Figure 7 Comparison of local profiles after denoising simulated post-stack data using different methods:

[0127] Figure 7a Profile diagram after adding random noise to the embodiment;

[0128] Figure 7b Comparison of local profiles after denoising by Contourlet transform;

[0129] Figure 7c Comparison of local profiles after denoising by Contourlet transform + cycle translation algorithm;

[0130] Figure 7d Comparison of local profiles after denoising by the technology of the present invention. Detailed implementation manners

[0131] The present invention is based on the Contourlet transform, combines the cycle translation algorithm and non-local means filtering (FNLM) to perform denoising operations on seismic data. To make the purpose, technical solution and advantages of the present invention clearer, the following technical roadmap is attached as Figure 1 and uses the classic Marmousi model (Figure 2 ) As a test model for simulating post-stack seismic data, embodiments are enumerated to further elaborate the present invention. The main steps are as follows:

[0132] A seismic data denoising method based on Contourlet transform includes the following steps:

[0133] Step 1: Seismic data loading and preprocessing

[0134] 1. Data parsing and initialization

[0135] Input parameters: filename (path of SEGY file with noise), sigma_contourlet (Contourlet threshold parameter), sigma_fnlm (non-local means smoothing coefficient);

[0136] Parse the SEGY file through the ReadSegy function to obtain the seismic data matrix seis and its header information;

[0137] Dynamically construct the processing matrix: Automatically calculate the optimal processing size n = 2^nextpow2(max(high, longth)) according to the input data size [high, longth].

[0138] 2. Data centered embedding

[0139] As Figure 3a and Figure 3b , adopt the symmetric padding strategy:

[0140] start_row = floor((n - high) / 2)+1;

[0141] start_col = floor((n - longth) / 2)+1;

[0142] Data_single(start_row:start_row + high - 1, start_col:start_col + longth - 1) = seis;

[0143] Step 2: Multi-scale circular shift Contourlet transform denoising

[0144] 1. Parameter system construction

[0145] Decomposition structure: pfilt = 'pkva' (dual-tree complex wavelet filter) / dfilt = 'pkva' (directional filter bank);

[0146] Decomposition levels: nlevs = [4, 4, 3] (the number of pyramid decomposition layers decreases gradually);

[0147] Noise modeling: rho = 3 (noise intensity coefficient), sigma = sig / rho (dynamic noise standard deviation).

[0148] 2. Circular translation optimization, such as Figures 4a to 4e ;

[0149] Perform 4 translation iterations (num_translations = 4, the number of loops can be adjusted: 2 - 8 (even numbers)), and the translation amount t ∈ [1, 4] each time;

[0150] Core processing flow:

[0151] Translate data → Contourlet decomposition → Hard threshold function processing → Inverse transform reconstruction → Translation anti - correction

[0152] The threshold calculation model uses hard threshold denoising:

[0153] cth = sigma_contourlet * sigma * sqrt(nvar)

[0154] where nvar is the Contourlet sub - band noise variance

[0155] Step 3: Improved non - local means filtering, such as Figure 5 ;

[0156] 1. Data normalization processing

[0157] Dynamic range compression:

[0158] normalized_data = (cumulative_result - data_min) / (data_max - data_min)

[0159] Retain the data distribution characteristics and avoid gradient disappearance

[0160] 2. FNLM parameter configuration

[0161] Smoothness control: sigma_fnlm (can be modified)

[0162] Similarity measurement window: default 8×8 search domain (can be modified)

[0163] Weighted coefficient calculation:

[0164] w_ij = exp(-||I(i)-I(j)||2 / (h2))

[0165] where h = sigma_fnlm * sqrt(max(high, longth)).

[0166] Step 4: Data Processing and Output

[0167] 1. Spatial Domain Restoration

[0168] Size Restoration: Crop to the size of the original data;

[0169] Edge Compensation: Perform weighted average processing on the filled area.

[0170] 2. SEGY Format Output

[0171] File Naming Rule: [Prefix]_[File Name]_contourlet.segy;

[0172] Metadata Inheritance: Completely retain the original SegyTraceHeaders information.

[0173] To make the purpose and actual effects of the present invention clearer and more understandable, the effects of each link of the present invention technology are compared.

[0174] Figures 6a - 6d

[0175] Figures 7a - 7d Shows the effect comparison diagrams of adding random noise to the original model, denoising by Contourlet transform, denoising by Contourlet transform + cyclic translation algorithm, and denoising by the technology of the present invention. It can be seen that to a certain extent, all three methods can suppress random noise, but the pseudo-Gibbs phenomenon of the Contourlet transform is serious; the Contourlet transform + cyclic translation algorithm removes the pseudo-Gibbs effect, but the noise removal quality is still not good; the denoising technology of the present invention removes the pseudo-Gibbs effect and effectively removes random noise at the same time, with better denoising quality. Figure 7b Figure 7c Figure 7d

[0176]

[0177]

[0178] In addition to the method of the present invention based on Contourlet transform + cyclic translation algorithm + non-local means filtering (FNLM), the following several alternative solutions can also be considered to achieve similar denoising goals:

[0177] 1. Joint Denoising Method Combining Multiple Transforms

[0178] Adopt Contourlet transform + Wavelet Transform, use Contourlet transform for multi-directional analysis, and combine with wavelet transform for multi-scale decomposition to enhance the representation ability of different frequency components of the signal.

[0179] This method can more comprehensively retain the structural characteristics of seismic signals, reduce the deficiency of wavelet transform in extracting direction information, and at the same time reduce the impact brought by the computational complexity of Contourlet transform.

[0180] 2. Adaptive denoising method based on deep learning

[0181] Use U-Net, Transformer or Denoising CNN (DnCNN) for end-to-end seismic data denoising.

[0182] Train a data-driven model, train the network with a large amount of noisy and denoised data pairs so that it can adapt to different types of noise and reduce the need for manual parameter adjustment.

[0183] This method can learn complex signal features and noise distribution patterns, but requires a large amount of training data and it is difficult to explain the denoising mechanism of the model.

[0184] 3. Adaptive denoising method based on BM3D (Block Matching 3D Filtering)

[0185] The BM3D (block matching three-dimensional denoising) method uses block matching to find similar signal blocks and performs denoising in the three-dimensional transform domain (such as DCT transform).

[0186] This method performs well in the field of image processing, can maintain the edge and texture features of the signal while reducing noise, and also has certain potential in seismic data denoising.

[0187] However, this method has a high computational complexity, a slow processing speed for large-scale data, and is difficult to be applied to real-time or high-efficiency data processing scenarios.

[0188] 4. Noise suppression method combined with statistical modeling

[0189] Use Bayesian Estimation or adaptive Kalman filtering for noise modeling.

[0190] These methods dynamically adjust the filtering parameters according to the statistical characteristics of the noise to achieve adaptive denoising for different types of noise.

[0191] It is applicable to the cases where the noise characteristics are known or can be modeled, but when the noise is complex or non-stationary, the denoising effect may be limited.

[0192] 5. Method of combining empirical mode decomposition (EMD) and adaptive filtering

[0193] EMD (Empirical Mode Decomposition) can decompose a signal into a series of intrinsic mode functions (IMFs) and suppress noise for different IMF components.

[0194] Combined with adaptive filtering techniques (such as Wiener filtering, RLS filtering, etc.), the denoising effect can be further optimized.

[0195] This method is applicable to non-stationary signals, but there may be mode mixing problems, resulting in limited noise separation effect.

Claims

1. A seismic data denoising method based on Contourlet transform, characterized in that, It includes the following steps: Step 1: Earthquake data loading and preprocessing; Step 2: Denoising using multi-scale cyclic translation Contourlet transform; Step 3: Using improved non-local means filtering; Step 4: Data processing and output.

2. The seismic data denoising method based on Contourlet transform according to claim 1, characterized in that, Step 1 includes the following processes: (1) Data parsing and initialization Input parameters: filename, sigma_contourlet, sigma_fnlm; Among them, filename is the path of the noisy SEGY file to be input, sigma_contourlet is the Contourlet threshold parameter, and sigma_fnlm is the non-local means smoothing coefficient; Parse the SEGY file through the ReadSegy function to obtain the seismic data matrix seis and its header information; Dynamically construct the processing matrix: Automatically calculate the optimal processing size according to the input data size [high, longth] n = 2^nextpow2(max(high, longth)); (2) Data centered embedding Adopt the symmetric padding strategy: start_row = floor((n - high) / 2) + 1; start_col = floor((n - longth) / 2) + 1; Data_single(start_row:start_row + high - 1, start_col:start_col + longth - 1) = seis.

3. A seismic data denoising method based on Contourlet transform according to claim 1, characterized in that, Step 2 includes the following processes: (1) Parameter system construction Decomposition structure: pfilt = 'pkva', dfilt = 'pkva'; Decomposition levels: nlevs = [4, 4, 3]; Noise modeling: rho = 3, sigma = sig / rho; Among them, pfilt and dfilt are the Laplacian pyramid filter and the directional filter bank, and can be replaced by several filters such as "9-7", "5-3", "Burt", "pkva"; nlevs is the number of layers specified for the directional filter bank DFB decomposition at each pyramid level, and the number of layers can be modified; rho is the noise intensity coefficient, generally default; sigma is the dynamic noise standard deviation; (2) Cyclic translation optimization Perform 4 translation iterations: num_translations = 4, where the number of cycles can be adjusted: 2 - 8, and it needs to be an even number. Each translation amount t ∈ [1, 4] Core processing flow: Translate data → Contourlet decomposition → Hard threshold function processing → Inverse transform reconstruction → Translation anti-correction The threshold calculation model uses hard threshold denoising: cth = sigma_contourlet * sigma * sqrt(nvar) where nvar is the Contourlet subband noise variance.

4. A seismic data denoising method based on Contourlet transform according to claim 1, characterized in that, Step 3 includes the following processes: [[ID= ​ ​ Preserve the data distribution characteristics and avoid gradient disappearance; (2) FNLM parameter configuration Smoothness control: sigma_fnlm; Similarity measurement window: default 8×8 search domain; Among them, both the sigma_fnlm and the similarity measurement window parameters can be modified according to the seismic data situation Weighted coefficient calculation: w_ij = exp(-||I(i)-I(j)||2 / (h2)) Where h = sigma_fnlm * sqrt(max(high, longth)).

5. A seismic data denoising method based on Contourlet transform according to claim 1, characterized in that, Step 4 includes the following processes: (1) Spatial domain restoration Size reduction: Crop to the original data size; Edge compensation: Perform weighted average processing on the filled area; (2) SEGY format output File naming rule: [prefix]_[filename]_contourlet.segy; Metadata inheritance: Completely retain the original SegyTraceHeaders information.

Images