A method for enhancing weak earthquake signals based on compressed sensing and statistical learning

Abstract

The invention discloses a method for enhancing weak seismic signals based on compressed sensing and statistical learning. The invention belongs to the field of seismic signal processing. The traditional curvelet domain denoising method suppresses random noise by thresholding the amplitude of the curvelet coefficient, while weak signals and random noise have the same weak amplitude characteristics in the curvelet domain, so it is difficult to distinguish weak seismic signals from random noise. The invention applies compressed sampling to noisy seismic data. After compressed sampling, the curvelet coefficients of the effective seismic signal have a significant amplitude reduction, while random noise does not. Based on this characteristic, a sensitive coefficient matrix is ​​designed, which is a statistical superposition matrix obtained by calculating the ratio of the original coefficient to the sampled coefficient after multiple subsampling analysis and outlier removal. In the curvelet domain, the original weak signal is effectively enhanced through the constraints of the sensitive coefficient matrix and the L1 norm matrix, and the noise coefficient in the sensitive coefficient matrix is ​​significantly lower than the effective signal, so that the method can effectively suppress noise while enhancing the weak signal. Both synthetic and actual seismic data tests show that the method is superior to the traditional curvelet domain method in terms of weak seismic signal enhancement. Compared with current similar technologies, this method can significantly enhance deep weak signals and improve the overall signal-to-noise ratio. The effect is far superior to traditional algorithms, and it provides a new idea for weak signal enhancement and noise suppression in deep seismic data.

Description

Technical Field

[0001] The present invention relates to the field of exploration geophysics and signal processing, and in particular to a seismic signal denoising and weak signal enhancement method combining a curvelet transform algorithm with a compressed sensing and statistical algorithm. Background Art

[0002] With the continuous deepening of seismic exploration, the research direction has shifted from conventional to unconventional, from medium and shallow layers to deep layers, and from simple structures to complex structures. The Tarim region in western my country is rich in oil and gas resources and has great exploration potential, but most of the resources are stored in deep or ultra-deep underground layers. The underground loess area in this area has serious absorption attenuation, and is affected by the undulating surface and the underlying high-speed layer. Various interference waves are developed and the effective reflection is weak. The quality of seismic signals received through seismic exploration is difficult to guarantee, which seriously affects the accuracy of full waveform inversion and prestack migration imaging, and brings great difficulties to the subsequent imaging and interpretation of geophysical structures. In order to cope with the processing and interpretation needs of seismic data of deep, steep structures and complex geological bodies in this area, it is necessary to further improve the resolution and signal-to-noise ratio of reflection seismic exploration data.

[0003] Existing methods for enhancing weak seismic signals are mainly from two aspects: one is to compensate for the attenuation amplitude, and the other is to suppress random noise. Traditional denoising methods are usually based on the characteristic differences between effective signals and noise in different transform domains, and design appropriate thresholds for filtering or other processing to suppress noise components. Although such methods can improve the quality of seismic records to a certain extent, in the denoising process, it is often necessary to adjust the corresponding parameters according to the actual noise characteristics to achieve the optimal denoising effect. Since the effective weak signals and noise interference frequency bands in the deep layers of seismic data are small and difficult to distinguish, further, traditional denoising methods often do not provide sufficient protection for weak signals, and even misjudge weak signals as noise removal, especially when the signal-to-noise ratio in the deep layers is extremely low, the quality of seismic data cannot be guaranteed, and it is difficult to meet imaging and interpretation requirements. Summary of the invention

[0004] In order to solve the problems existing in the prior art, the present invention provides a method for enhancing weak seismic signals based on compressed sensing and statistical learning, which can effectively enhance weak signals while suppressing random noise.

[0005] The original noisy data is compressed and sampled using the observation matrix. The influence of multiple subsampling on the curvelet coefficient is studied in the curvelet domain. That is, the amplitude of the curvelet coefficient of the weak seismic signal is significantly reduced after compression sampling, while the influence of compression sampling on the curvelet coefficient of random noise is very small. According to this characteristic, combined with statistical methods, a sensitivity matrix is ​​designed, which is the absolute ratio of the curvelet coefficient before and after compression sampling. In the curvelet domain, the sensitivity coefficient can better distinguish weak seismic signals from random noise by thresholding the amplitude of the curvelet coefficient. The results of the application of actual seismic data show that the proposed method can effectively enhance weak seismic signals while denoising, and the processing effect is significantly better than the traditional curvelet domain denoising method, which effectively improves the signal-to-noise ratio and resolution of seismic data, especially deep seismic data.

[0006] The specific technical solution adopted by the present invention is:

[0007] First, the original seismic data containing noise is subjected to discrete curvelet transformation to obtain the original curvelet coefficients. Secondly, the original seismic data is compressed and sampled using a measurement matrix (e.g., a Gaussian random matrix), and the curvelet coefficients after compression sampling are calculated; the ratio of the original curvelet coefficients to the corresponding elements of the curvelet coefficients after subsampling is calculated to obtain a ratio matrix; the outliers of the ratio matrix are calculated according to statistical methods, and the outliers are defined as elements that are greater than the upper quartile (75%) by more than 1.5 interquartile differences or less than the lower quartile (25%) by more than 1.5 interquartile differences. Set the size of the logic matrix equal to the curvelet coefficients, and set the corresponding positions of the logic matrix to 1 according to the outliers found, to represent the positions of the curvelet coefficient elements with large changes in the curvelet coefficients after subsampling. According to analysis, the effective seismic signal is sensitive to subsampling and insensitive to noise, that is, the corresponding positions of the logic matrix indicate that the curvelet domain transformation of the effective signal is large. Thirdly, in order to eliminate the randomness of single sampling, the number of sub-sampling is set to 50-100, and the logical matrices obtained each time are superimposed, and the sensitivity coefficient matrix is ​​obtained after iteration; then, the discrete Fourier transform of the unit matrix with the same size as the original data is calculated to obtain the norm matrix of the curvelet coefficients; finally, the obtained sensitivity matrix and norm matrix are constrained on the original curvelet coefficients, and the discrete curvelet inverse transform is implemented to obtain the processing results after denoising and weak signal enhancement.

[0008] Compared with the existing methods, the present invention implements multiple sub-sampling of the original data based on the compressed sensing theory, and obtains the designed sensitivity matrix through statistical analysis, which involves fewer parameters that need to be set manually; the sensitivity matrix is ​​designed according to the sensitivity of the effective signal and noise to sub-sampling. In theory, the maximum value is the set number of iterations. The closer to the upper limit value, the more it indicates a valid signal, and the smaller the value, the greater the possibility of noise. Simulation and actual data processing results show that the present invention can effectively enhance weak signals while removing noise, especially in deep data, showing excellent weak signal enhancement capabilities. BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 is a flow chart of the method of the present invention;

[0010] Figure 2 Content diagram of earthquake weak signal enhancement method based on compressed sensing and statistical learning;

[0011] Figure 3 This is a schematic diagram of the model of earthquake weak signal enhancement method based on compressed sensing and statistical learning;

[0012] Figure 4 The comparison diagram of the height value transformation before and after the reconstruction of a certain curvelet coefficient tuple by the designed sensitivity matrix and different methods;

[0013] (a) Contour distribution of the original curvelet coefficient C{2}{2}; (b) Sensitivity matrix σ N {2}{2}Contour distribution;

[0014] (c) The CCS{2}{2} contour distribution after reconstruction by our method; (d) CCS{2}{2} after BCD threshold processing BCD {2}{2}Contour distribution;

[0015] Figure 5 This is a comparison chart of the actual data processing effect;

[0016] (a) Original noisy seismic profile; (b) BCD denoising result; (c) Processing result of this method;

[0017] Figure 6 for Figure 5 The deep local magnified comparison picture shown;

[0018] (a) Original noisy seismic data; (b) BCD method; (c) the proposed method;

[0019] Figure 7 The following is a spectrum diagram comparing the noise suppression results of actual seismic data using different methods. DETAILED DESCRIPTION

[0020] The present invention discloses a method for enhancing weak seismic signals based on compressed sensing and statistical learning, and belongs to the field of seismic data processing. The present invention designs a sensitive matrix for characterizing signal validity through compressed sensing theory and statistical analysis theory. The sensitivity of different curve domain coefficients is obtained by multiple subsampling of original data, and the original curve coefficients are constrained by the sensitive matrix, so that noise can be suppressed and weak signals can be enhanced. Weak signal enhancement and noise suppression of deep low signal-to-noise ratio data are achieved, and the quality of seismic data in complex work areas is effectively improved.

[0021] The following introduces the specific implementation steps of the present invention and makes further specific descriptions in combination with examples and drawings. According to the following description, the content, characteristics, advantages and effects of the present invention will become more clear.

[0022] Figure 1 The method flow chart of the present invention is as follows:

[0023] Step 101: Select seismic data. If it is post-stack data, outliers can be removed, data normalization and conventional filtering can be performed; the number of data channels M and the number of sampling points N are saved as scale and angle calculation parameters for discrete curvelet inverse transform in step 102;

[0024] Step 102: Perform a curvelet forward transform on the noisy seismic data to obtain original curvelet coefficients c(j, l, k).

[0025] The actual collected seismic data contains effective signals and noise, and the composition of seismic data can be described by the following expression:

[0026] f(t)=o(t)+n(t)

[0027] Where f(t) represents the noisy seismic data, o(t) represents the original seismic data without noise, and n(t) represents the noise. Where x is the spatial domain factor, j is the scale factor, j = 1, 2, ..., r. The total scale number r is determined by the number of seismic traces M and the number of sampling points N of the original seismic data f(t).

[0028]

[0029] in, This is a round-up operation. Fourier transform of With the Fourier window function U j (ω)Satisfaction ω is the frequency domain factor, and all curve functions can be obtained through Translation and rotation are obtained. Introduce the displacement factor sequence k=(k1,k2)∈Z 2 , rotation angle and location Where 0≤θ l <2π, l is the angle factor, and l=0,1,…, is a floor operation, then the curve function It can be expressed as:

[0030]

[0031] In the formula, the rotation matrix

[0032]

[0033] For the original seismic signal f(t), the original curvelet coefficient c is obtained by curvelet forward transformation in the real domain: j,l,k ,

[0034]

[0035] In the formula, is the Fourier transform of f(t).

[0036] Step 103: Select a sparse matrix to subsample the original profile

[0037] In the theory of compressed sensing, the most widely used is the random Gaussian measurement matrix, which is designed as follows: construct an MxN matrix so that Each element in is independently subject to a Gaussian distribution with a mean of 0 and a variance of 1 / M, that is:

[0038]

[0039] The measurement matrix has strong randomness. It can be proved that when the number of measurements of the random Gaussian measurement matrix M ≥ cKlog N / K, it will satisfy the finite isometry property (RIP condition) with great probability. In the process of compressed sensing, the random Gaussian measurement matrix is ​​widely used mainly because it is irrelevant to most orthogonal bases or orthogonal dictionaries, and the number of measurements required for accurate reconstruction is relatively small.

[0040] If the number of sampling points of the original data is M>N, set the subsampling matrix (G represents the Gaussian measurement matrix); otherwise, Transpose.

[0041] Get the seismic data after subsampling:

[0042]

[0043] Step 104: Perform discrete curvelet transformation on the sub-sampled seismic profile to obtain curvelet coefficients

[0044]

[0045] Step 105: Construct sensitivity coefficient matrix

[0046] Compare the curvelet coefficients before subsampling and the coefficients after subsampling to get the ratio matrix δ

[0047]

[0048] The outliers of the ratio matrix are calculated according to the statistical method. The outliers are defined as elements that are greater than 1.5 interquartile differences from the upper quartile (75%) or less than 1.5 interquartile differences from the lower quartile (25%). The sensitivity coefficient matrix σ is set to a logical matrix whose size is equal to the curvelet coefficient. The corresponding position of σ is set to 1 according to the outliers found, indicating the position of the curvelet coefficient element with large curvelet coefficient transformation after subsampling.

[0049] Step 105: Obtain the sensitivity matrix through statistical analysis of the curvelet coefficients after multiple subsampling

[0050] In order to eliminate the randomness of single sampling, the number of sub-sampling is set to 50-100, and the logic matrix obtained each time is superimposed. After iteration, the sensitivity coefficient matrix σ is obtained. N , the maximum value of a single element is N, and the minimum value is 0; the larger the value, the greater the influence of subsampling on its curve domain coefficient, representing the effective signal, and the smaller the value, the smaller the influence of subsampling, representing the random noise contained in the data.

[0051] Step 106: Obtain the de-noised seismic profile using the sensitivity matrix and inverse curvelet transform

[0052] Reconstructing the curvelet coefficients using the sensitivity coefficient matrix

[0053]

[0054] Where E is the L1 norm of the seismic data. The new curvelet coefficients are subjected to inverse curvelet transform to obtain the seismic data after denoising and weak signal enhancement. The inverse curvelet transform formula is as follows:

[0055]

[0056] Finally, the present invention evaluates the processing results by comparing the spectrum diagrams.

[0057] The implementation process of the present invention is demonstrated below in conjunction with specific embodiments. Figure 3 This is a schematic diagram of a model of specific implementation steps provided by the present invention.

[0058] Step 1: In order to demonstrate the denoising and weak signal enhancement effects of the present invention, actual seismic data from a certain work area in the west of China is used as the object, and the effectiveness of the present invention is further illustrated by comparing it with the Bayesian adaptive threshold denoising method (BCD method).

[0059] Step 2: Set the scale and angle of the curvelet transform according to the number of seismic data channels and sampling points, and perform the curvelet forward transform to obtain the original curvelet coefficients;

[0060] Step 3: Sub-sample the original data using a random Gaussian measurement matrix;

[0061] Step 4: Perform a curvelet forward transform on the sub-sampled data to obtain the sub-sampled curvelet coefficients;

[0062] Step 5: Calculate the ratio of the elements corresponding to the original coefficient and the sub-sampling coefficient to obtain the ratio matrix. According to the quartile statistical analysis, obtain the sensitive logic matrix of this sub-sampling;

[0063] Step 6: Set the number of sub-sampling times N=50, repeat steps 3-5, and superimpose the sensitive logic matrix obtained by each sub-sampling to obtain a sensitive matrix; Figure 4 The sensitivity matrix obtained through multiple iterations and the height value transformation comparison diagram of a certain curvelet coefficient tuple before and after reconstruction by different methods show that the curvelet coefficient reconstructed by the proposed method is not for a fixed frequency band or value, but constrains the noise through the full frequency band of the sensitivity matrix. The details of Figure (4-b) are richer, and the curvelet coefficient for weak signals is clearly expressed; compared with the proposed method ( Figure 4 -c) and the BCD method ( Figure 4 -d) The reconstructed curvelet coefficients. It can be seen that the protection of weak signal coefficients in this paper is significantly better than that of the BCD method;

[0064] Step 7: Constrain the original curvelet coefficients according to the sensitivity matrix to obtain the reconstructed curvelet coefficients, and perform inverse curvelet transform to obtain the processed seismic data.

[0065] Step 8: Compare the processing results of the present invention and the curvelet domain BCD method, and globally compare Figure 5 , deep local contrast Figure 6 The proposed method has obvious denoising effect in the deep layer of 3000-5000 points, and the denoising and weak signal enhancement capabilities in the deep layer below 5000 points are obviously better. The comparison results show that the proposed method can effectively improve the quality of low signal-to-noise ratio data in deep layers in complex areas.

[0066] Step 9: Calculate the spectrum of the entire seismic data and compare them. Comparison of the global spectrum of the original seismic data with the BCD method and the proposed method Figure 7 As shown in the figure, the BCD method has a certain suppression effect on noise, but due to the threshold processing method, the signal energy in the effective frequency band is also lost. The proposed method has a significant compensation effect on seismic data within the effective frequency band, and the high-frequency noise suppression effect is significantly better than the BCD method. The comparison results show that this method can effectively suppress noise signals while enhancing weak signals.

[0067] Through the above process, it can be seen that the present invention provides a method for enhancing weak seismic signals based on compressed sensing and statistical learning, and designs a sensitive matrix that characterizes the effectiveness of signals through compressed sensing theory and statistical analysis theory. Through the reconstruction of the curvelet coefficients, noise suppression and formation weak signal enhancement are achieved. The comparison results show that the method proposed by the present invention significantly improves the signal-to-noise ratio, and its effect far exceeds that of traditional algorithms.

Claims

1. A method for enhancing weak earthquake signals based on compressed sensing and statistical learning, characterized in that: Performing a forward curvelet transform on the original seismic data containing noise to obtain the original curvelet coefficients; sub-sampling the original seismic data using a measurement matrix and calculating the curvelet coefficients after sub-sampling, the measurement matrix includes a Gaussian random matrix and a local Fourier transform matrix; calculating the ratio of the original curvelet coefficients to the corresponding elements of the curvelet coefficients after sub-sampling to obtain a ratio matrix; calculating the outliers of the ratio matrix using a statistical method; Set the size of the logic matrix to be equal to the curvelet coefficient, and set the corresponding position of the logic matrix to 1 according to the outlier value found; in order to eliminate the randomness of single sampling, set the number of sub-sampling to 50-100, superimpose the logic matrix obtained each time, and obtain the sensitivity coefficient matrix after iteration; calculate the discrete Fourier transform of the unit matrix of the same size as the original data to obtain the L1 norm matrix of the curvelet coefficient; constrain the original curvelet coefficient with the obtained sensitivity coefficient matrix and norm matrix, and implement the inverse curvelet transform to obtain the processing results after denoising and weak signal enhancement; The specific implementation steps of this method are as follows: Step 101: Select seismic data, save the data channel number M and sampling point number N: Step 102: Performing a curvelet forward transform on the noisy seismic data to obtain original curvelet coefficients c(j, l, k); The actual collected seismic data contains effective signals and noise. The following expression is used to describe the composition of seismic data: f(t)=o(t)+n(t); Among them, f(t) represents the noisy seismic data, o(t) represents the original seismic data without noise, and n(t) represents the noise; the mother curve wave is defined as Where x is the spatial domain factor, j is the scale factor, j = 1, 2, ..., r; the total scale number r is determined by the number of seismic traces M and the number of sampling points N of the original seismic data f(t); in, It is the rounding operation; Fourier transform of With the Fourier window function U j (ω)Satisfaction ω is the frequency domain factor, and all curve functions are passed through Translation and rotation are obtained; the displacement factor sequence k=(k1,k2)∈Z is introduced 2 , rotation angle and location Where 0≤θ l <2π, l is the angle factor, and l=0,1,…, is a floor operation, then the curve function It is expressed as: Where the rotation matrix is: For the original seismic signal f(t), the original curvelet coefficient c is obtained by curvelet forward transformation: j,l,k ; In the formula, is the Fourier transform of f(t); Step 103: Select a measurement matrix to subsample the original profile; If M>N, set the measurement matrix G represents a Gaussian random matrix, otherwise, Perform transposition; Get the seismic data after subsampling: Step 104: performing a forward curvelet transformation on the sub-sampled seismic profile to obtain a curvelet coefficient; Step 105: construct a sensitivity coefficient matrix; Compare the curvelet coefficients before subsampling and the coefficients after subsampling to obtain the ratio matrix δ; The outliers of the ratio matrix were calculated according to the statistical method, and the outliers were defined as the elements that were greater than 75% of the upper quartile by more than 1.5 interquartile differences or less than 25% of the lower quartile by more than 1.5 interquartile differences; Set δ(j,l,k) to a logical matrix whose size is equal to the curvelet coefficient. Set the corresponding position of δ(j,l,k) to 1 according to the outlier value found, indicating the position of the curvelet coefficient element with large curvelet coefficient transformation after subsampling. The sensitivity coefficient matrix is ​​obtained by statistical analysis of curvelet coefficients after multiple subsampling; Set the number of subsampling times to 50-100, superimpose the logic matrix obtained each time, and obtain the sensitivity coefficient matrix σ after iteration. N , the maximum value of a single element is N, and the minimum value is 0; Step 106: Obtain the de-noised seismic profile using the sensitivity coefficient matrix and inverse curvelet transform; Reconstruct the curvelet coefficients using the sensitivity coefficient matrix: Where E is the L1 norm of the seismic data; the new curvelet coefficients are subjected to inverse curvelet transform to obtain the seismic data after denoising and weak signal enhancement. The inverse curvelet transform formula is as follows: Finally, the processing results are evaluated by comparing the spectrum diagrams.

2. The method for enhancing weak earthquake signals based on compressed sensing and statistical learning according to claim 1, characterized in that: The observation matrix that satisfies the compressed sensing theory is used as the measurement matrix, and the local Fourier transform matrix is ​​more suitable for sub-sampling of seismic data in the domain.

Images