A method, device and electronic device for de-noising seismic signals

Abstract

The present application proposes a method, device and electronic device for denoising seismic signals. The embedding dimension is determined by calculating the power spectrum density of the seismic signal. For the selection of the denoising order k, the cumulative energy contribution rate of the signal is first calculated, and the denoising order is determined by using an adaptive threshold based on the noise variance. Then, the signal components are screened by using the kurtosis value to remove the components dominated by Gaussian noise to obtain a preliminary denoised signal. Then, the preliminary denoised signal is subjected to wavelet decomposition, and the signal is screened by using the kurtosis value to remove the components dominated by Gaussian noise. In the signal-dominated components, the wavelet coefficients are processed by a preset evaluation model and a preset wavelet estimator, and a denoised signal is obtained by inverse transformation. The method described in the present application can efficiently and accurately denoise seismic signals.

Description

Technical Field

[0001] The present application relates to the technical field of seismic data processing, and in particular to a seismic signal denoising method, device and electronic equipment. Background Art

[0002] Seismic signal denoising has always been a hot topic in the field of seismic research. In recent years, wavelet threshold denoising has been studied in depth by many scholars because it can effectively remove noise in the same frequency band. However, this method currently has the following problems:

[0003] (1) The commonly used wavelet transform method is discrete wavelet transform. Although it has high computational efficiency, it will cause aliasing of the wavelet coefficients of the signal during the downsampling process. After the wavelet coefficients are processed by the threshold method, the contraction factor near the singular point of the signal has a serious pseudo-Gibbs phenomenon, which makes the signal discontinuous or oscillating after reconstruction. Although continuous wavelet transform methods such as continuous wavelet transform and synchronous compression wavelet transform can solve the pseudo-Gibbs phenomenon, the transform will produce too many redundant coefficients, resulting in low overall computational efficiency of the algorithm.

[0004] (2) The selection of thresholds largely determines the denoising effect of the algorithm. Traditional item-by-item processing methods, such as soft / hard thresholds, maximum and minimum thresholds, unbiased risk estimation thresholds, heuristic thresholds, etc., will produce singular points. In addition, seismic signals are continuous within a limited frequency band in a limited time. The item-by-item processing method is easy to remove effective signals, reducing the denoising effect.

[0005] Therefore, how to efficiently and accurately denoise seismic signals has become an urgent problem that needs to be solved. Summary of the invention

[0006] The present application provides a seismic signal denoising method, device and electronic equipment, which can achieve efficient and accurate denoising of seismic signals.

[0007] In a first aspect, an embodiment of the present application provides a method for denoising a seismic signal, the method comprising: acquiring a seismic signal; determining a trajectory matrix corresponding to an embedding dimension based on a power spectral density and an embedding dimension calculation formula, the power spectral density being determined according to the seismic signal; determining an energy contribution rate based on the trajectory matrix, singular spectrum analysis and an energy contribution rate calculation formula; determining a denoising order based on the energy contribution rate and an adaptive threshold adjustment mechanism; determining a signal component based on the denoising order and a signal component calculation formula; determining a kurtosis value corresponding to the signal component based on the signal component and a kurtosis value calculation formula; determining preliminary denoising data based on the kurtosis value and a screening formula, and determining denoised data by processing the preliminary denoising data through a complex wavelet block threshold.

[0008] In a possible implementation, the processing of the preliminary denoised data through a complex wavelet block threshold to determine the denoised data includes: determining multi-scale wavelet coefficients through dual-tree complex wavelet transform processing based on the preliminary denoised data; determining a valid signal through the screening formula based on the multi-scale wavelet coefficients; determining an optimal block size and a target threshold of the valid signal through a preset evaluation model based on the valid signal; and determining the denoised data based on the optimal block size and the target threshold, the multi-scale wavelet coefficients and a preset wavelet estimator.

[0009] In a possible implementation, the embedding dimension calculation formula is:

[0010]

[0011] Among them, fmax represents the frequency corresponding to the maximum peak when the power spectrum density of the seismic signal is calculated, Fs represents the sampling frequency of the signal, a represents the adjustment factor, d represents the embedding dimension, and n represents the signal length.

[0012] In a possible implementation, the energy contribution rate is determined based on the trajectory matrix, singular spectrum analysis and energy contribution rate calculation formula, including: performing autocorrelation analysis on the trajectory matrix to obtain a covariance matrix; performing singular spectrum decomposition on the covariance matrix to determine the modified singular values ​​corresponding to the covariance matrix; and determining the energy contribution rate based on the modified singular values ​​using the energy contribution rate calculation formula.

[0013] In a possible implementation, the energy contribution rate calculation formula is:

[0014]

[0015] Among them, η represents the energy contribution rate, λ' i represents the modified singular value corresponding to the covariance matrix, and k represents the denoising order.

[0016] In a possible implementation, the kurtosis value calculation formula is:

[0017]

[0018] Among them, μ 4 represents the fourth-order center distance, σ represents the standard deviation of noise; x represents the time series corresponding to each signal component after singular spectrum decomposition, and kurt represents the kurtosis value.

[0019] In a possible implementation, the screening formula is:

[0020]

[0021] Among them, a represents the adjustment factor and 24 / n represents the variance.

[0022] In a possible implementation, the preset evaluation model is:

[0023]

[0024] Among them, σ represents the standard deviation of noise, L represents the block size, λ represents the threshold level, Wy represents the wavelet coefficient of the noisy signal, N represents the data length, λ F =2L*ln N.

[0025] In a second aspect, an embodiment of the present application provides a seismic signal denoising device, the device comprising: an acquisition module for acquiring seismic signals; a determination module for determining a trajectory matrix corresponding to an embedding dimension based on a power spectrum density and an embedding dimension calculation formula, wherein the power spectrum density is determined according to the seismic signal; the determination module is also used to determine an energy contribution rate based on the trajectory matrix, singular spectrum analysis and an energy contribution rate calculation formula; the determination module is also used to determine a denoising order based on the energy contribution rate and an adaptive threshold adjustment mechanism; the determination module is also used to determine a signal component based on the denoising order and a signal component calculation formula; the determination module is also used to determine a kurtosis value corresponding to the signal component based on the signal component and a kurtosis value calculation formula; the determination module is also used to determine preliminary denoising data based on the kurtosis value and a screening formula, and determine denoised data by processing the preliminary denoising data through a complex wavelet block threshold.

[0026] In a third aspect, an embodiment of the present application provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the method described in the first aspect or any one of the implementation methods thereof is implemented.

[0027] In a fourth aspect, an embodiment of the present application provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the method described in the first aspect or any one of the implementation methods thereof is implemented.

[0028] The present application obtains seismic signals and then uses improved singular spectrum analysis and complex wavelet block threshold for joint denoising. First, the singular spectrum analysis method is improved. In view of the problem of empirical selection of embedded dimensions in the original algorithm, power spectrum density is used for adaptive adjustment; for the selection of denoising order, the cumulative energy contribution rate of the signal is first calculated, and the denoising order is determined by an adaptive threshold based on noise variance, and then the signal components are screened by kurtosis value to remove the components dominated by Gaussian noise to obtain preliminary denoised data. In view of the problem that the improved singular spectrum analysis cannot suppress the same-band noise, an algorithm based on complex wavelet block threshold is used to remove the same-band noise. In this method, double-tree complex wavelet transform is first used for wavelet decomposition, and then the kurtosis value is used to screen the signal to remove the components dominated by Gaussian noise. In the signal-dominated components, the optimal block size and target threshold are obtained by a preset evaluation model (stein's unbiased risk estimation), and then the preset wavelet estimator (wavelet estimator based on James-Stein shrinkage criterion) is substituted to process the wavelet coefficients, and the denoised signal is obtained after inverse transformation. By comparing with other existing methods, the method provided in this application can efficiently and accurately denoise seismic signals, thereby improving the quality of seismic signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0030] Figure 1 A schematic diagram of a flow chart of a seismic signal denoising method provided in one embodiment of the present application;

[0031] Figure 2 The time domain and frequency domain diagrams of the seismic signal and the noisy signal with a noise variance of 0.5 provided in this application;

[0032] Figure 3 This is the SSA decomposition result diagram provided by this application when the kurtosis value is not added;

[0033] Figure 4 The SSA decomposition result diagram when adding the kurtosis value provided by this application;

[0034] Figure 5 A comparison chart of different denoising methods provided for this application;

[0035] Figure 6 A comparison chart of the residuals of different denoising methods provided in this application;

[0036] Figure 7 Time domain and time-frequency domain diagrams of different signals provided for this application;

[0037] Figure 8 A signal amplitude diagram after smoothing and threshold processing provided by the present application;

[0038] Fig. 9 A comparison chart of denoising results with and without adding kurtosis value provided in this application;

[0039] Fig.10 This is the signal diagram after SSA+BT-DTCWT joint denoising provided in this application;

[0040] Fig.11 Performance diagram of various denoising methods under different input noise variances provided by this application;

[0041] Fig.12 A structural block diagram of a seismic signal denoising device provided in one embodiment of the present application;

[0042] Fig.13 A schematic diagram of the structure of an electronic device provided in one embodiment of the present application. DETAILED DESCRIPTION

[0043] The following will be combined with the implementation methods of this application to clearly and completely describe the technical solutions in the implementation methods of this application. Obviously, the described implementation methods are only part of the implementation methods of this application, not all of the implementation methods. Based on the implementation methods in this application, all other implementation methods obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0044] After the seismic signal is detected, it is necessary to further use a denoising algorithm to process the noisy signal, especially the signal with low signal-to-noise ratio, to provide high-quality data for subsequent data analysis. The noise processed in this application is the most common random noise in seismic records, including environmental noise, secondary biological noise, and system noise. It is characterized by unpredictable apparent velocity, unfixed propagation direction, and uncertain frequency distribution range. In some cases, seismic signals, especially P-wave signals, will be covered by noise, which greatly affects the quality of the data. The noise exists all the time and is spread throughout the entire spectrum of the data, and there is a lot of low-frequency noise.

[0045] There are many methods for denoising seismic signals:

[0046] 1. Methods based on signal decomposition. The signal decomposition method aims to decompose a complex signal into several simpler components for easy analysis and processing. Currently, there are mainly methods based on modal decomposition and methods based on matrix rank reduction. The modal decomposition method is an adaptive signal analysis method. Its advantage is that it can extract and select the intrinsic mode functions (IMFs) existing in the signal according to different signal characteristics, so as to further analyze, process and interpret it.

[0047] 2. Methods based on time-frequency domain transformation. In the field of signal processing, time-frequency domain transformation denoising technology is an effective method that utilizes the difference in amplitude between the transformed coefficients of the signal and the noise. In the transform domain, the coefficients of the signal usually have a large amplitude, while the noise coefficients are relatively small. By applying hard threshold or soft threshold technology, noise can be effectively suppressed while retaining the important features of the signal. The key to this method is to select a suitable transform basis and threshold processing strategy, as well as the subsequent inverse transform process to ensure that the denoised signal retains its original information content. Sparse transform is an important form of transform domain denoising. Its core idea is to represent the signal as a linear combination of a series of transform coefficients so that most of these coefficients are zero, thereby achieving sparse representation of the signal. Wavelet transform is a typical example of sparse transform. It can capture the local characteristics of the signal through multi-scale analysis and is very suitable for processing non-stationary signals.

[0048] When using wavelet threshold denoising, the wavelet transform method and the selection of threshold have a great influence on the denoising result. Currently, the commonly used wavelet transforms are continuous wavelet transform (CWT) and discrete wavelet transform (DWT). Subsequently, a series of transform methods such as synchronized compressed wavelet transform (SS-CWT) and empirical wavelet transform (EWT) have gradually developed. There are three commonly used denoising methods: modulus maximum denoising, scale correlation denoising and threshold denoising. Among them, the amplitude-based soft and hard threshold denoising proposed by Donoho is the most widely used.

[0049] 4. Based on the filter denoising method, the filters here include analog filters and digital filters, such as high-pass, band-pass, low-pass and band-stop filters, which are designed according to the frequency characteristics of seismic signals and noise to remove noise in a specific frequency range. The key point of this method is to determine the selection of the filter and its corresponding frequency range. The selection of the filter is one of the core issues of filter processing. For different types of research content, the filters selected by various research institutions are also different. From the perspective of causality, the filtering methods can be divided into causal and non-causal filtering. Among them, causal filtering uses unidirectional filtering, which will change the phase spectrum of the data; while non-causal filtering will not change the phase spectrum of the data. It filters in the forward and backward directions in the entire time domain, and in the process of reverse filtering, the phase change introduced by the first filtering is offset, so that the data does not produce phase changes as a whole.

[0050] 3. Denoising method based on deep learning. This method utilizes the powerful learning ability of deep neural networks and achieves signal denoising by learning from a large amount of noisy and non-noisy data. This method does not need to rely on accurate noise models or other prior information, but automatically learns how to recover clean signals from noisy data in a data-driven way. This method is divided into four categories according to the structure of the network model: autoencoder, convolutional neural network, recurrent neural network and generative adversarial network. It can learn denoising tasks through end-to-end training. First, it requires a large number of noisy signals, corresponding pure signals and noisy signals as data sets. At present, the research is still immature and there are some shortcomings: (1) There is a lack of a large number of high signal-to-noise ratio signals as training sets; (2) The generalization ability is poor. Due to the regional characteristics of seismic signals, the overall denoising effect is general.

[0051] In order to efficiently and accurately remove noise from seismic signals, the method used in this application combines improved singular spectrum analysis and complex wavelet block threshold technology, aiming to effectively reduce the noise of seismic signals. This method first implements adaptive bandpass filtering through singular spectrum analysis, and then removes noise in the same frequency band as the signal through complex wavelet block threshold technology, thereby improving the signal-to-noise ratio of the signal and achieving more accurate seismic data processing.

[0052] Figure 1 A schematic flow chart of a seismic signal denoising method provided in one embodiment of the present application, the method comprising:

[0053] S110, obtaining seismic signals.

[0054] The earthquake and noise data used in this application are all raw data recorded by seismographs downloaded from the Global Seismic Network. After obtaining the acceleration recording data, the instrument response should be removed first. This is because the raw data output by the modern seismic recording system has been converted from acceleration to voltage, and then from voltage to count. The physical vibration of the surface is convolved with the instrument response of the sensor to obtain an electrical signal that is easy to transmit. The value is in counts and has no actual physical meaning. In order to obtain the physical quantities (displacement, velocity and acceleration) of the actual vibration of the surface, it is necessary to perform a deconvolution operation on the raw data, that is, to remove the instrument response.

[0055] The instrument response of the sensor is not an ideal constant. It will drop sharply from an extremely high gain before Fs / 2. This causes the data amplitude of the high-frequency band to increase sharply relative to the other bands when performing the instrument response removal operation, and may even submerge the seismic signal in the high-frequency noise. Therefore, when performing instrument response, it is necessary to set the range of the bandpass filter in advance to remove high-frequency signals. This application mainly processes seismic recording data with a sampling frequency of 200Hz. The instrument response begins to drop sharply from 80Hz. Therefore, it is necessary to set the high-frequency cutoff frequency of the bandpass filter below 80Hz.

[0056] Therefore, the present application obtains seismic signals by removing instrument responses based on data recorded by seismographs downloaded from the global seismic network.

[0057] S120, based on the power spectrum density and the embedding dimension calculation formula, determining the trajectory matrix corresponding to the embedding dimension, the power spectrum density is determined according to the seismic signal.

[0058] This application uses an improved singular spectrum analysis (SSA) method to preliminarily process seismic signals and obtain noise reduction data. The core of SSA is to construct the time series signal into a trajectory matrix through Takens embedding theorem, and perform singular value decomposition on its covariance matrix to obtain the corresponding eigenvalues ​​and eigenvectors, and then reconstruct the eigenvectors corresponding to some eigenvalues ​​to obtain multiple signal components, thereby realizing signal decomposition and removing some noise.

[0059] The choice of embedding dimension has a significant impact on the results and performance of signal decomposition. Too high a dimension will not only increase the calculation time, but also lead to over-decomposition of the signal; too low a dimension may not be able to effectively separate the different components in the signal, and the decomposition performance is average. In order to obtain better signal decomposition effects and avoid empirical selection of embedding dimensions, this application calculates the power spectral density (PSD) of the seismic signal X to obtain the energy distribution of each frequency component in the signal. According to the analysis results of PSD, the embedding dimension is adaptively selected to achieve reconstruction of the trajectory matrix and ensure that SSA can fully decompose the seismic signal.

[0060] For a known time series signal x = {x1, x2, x3, ···, xn}, i.e., a seismic signal, where n represents the length of the signal. First, the one-dimensional time series is constructed into a multidimensional time series matrix by the Takens embedding theorem, which is expressed as:

[0061]

[0062] Where d is the embedding dimension, τ is the delay time, and the window length m = n-(d-1)·τ. In the process of constructing the trajectory matrix X, there are two main parameters d and τ. This application adopts the method of selecting the embedding dimension and the delay time. The expression for calculating the embedding dimension is as follows:

[0063]

[0064] In the formula, fmax is the frequency corresponding to the maximum peak obtained by calculating the PSD of the original signal, Fs is the sampling frequency of the signal, and a is the adjustment factor. According to the characteristics of the seismic signal, the adjustment factor is tested multiple times and finally a=1.8 is selected; the delay time τ is usually 1.

[0065] S130, determining the energy contribution rate based on the trajectory matrix, singular spectrum analysis and energy contribution rate calculation formula, and determining the noise reduction order based on the energy contribution rate and the adaptive threshold adjustment mechanism.

[0066] In one possible implementation, the energy contribution rate is determined based on the trajectory matrix, singular spectrum analysis and the energy contribution rate calculation formula, including: performing autocorrelation analysis on the trajectory matrix to obtain the covariance matrix; performing singular spectrum decomposition on the covariance matrix to determine the modified singular values ​​corresponding to the covariance matrix; and determining the energy contribution rate based on the modified singular values ​​through the energy contribution rate calculation formula.

[0067] After completing the construction of the trajectory matrix X, the dimension of the matrix is ​​m×d. To further construct the required matrix expression, the matrix X is autocorrelated and the covariance matrix A is obtained:

[0068] A=XT X(3) At this time, the matrix dimension of A is d×d, and then the covariance matrix A is decomposed into a singular spectrum:

[0069] A=UΛV T (4)

[0070] In the formula, U and V are the eigenvector matrices of matrix A, Λ=diag(λ1,…,λd), where diag(·) represents a diagonal matrix, λ1,…,λd are the corresponding eigenvalues ​​of A, and λ1≥λ2≥···≥λd>0.

[0071] Ideally, the covariance matrix A can be regarded as the superposition of the signal matrix and the noise matrix, that is:

[0072] A=S+N(5)

[0073] In the formula, S represents the signal-dominated part of the matrix A, and N represents the noise-dominated part. The following assumptions are made: (1) STN = 0, indicating that the useful signal and noise are independent of each other and have no correlation; (2) is the noise power, I d is the identity matrix.

[0074] Therefore, the singular spectrum decomposition of the covariance matrix A and the signal matrix S is expressed as:

[0075]

[0076]

[0077] in, The following relationship exists between the block matrices:

[0078]

[0079] Λ 2 =σ W I d-k (9)

[0080]

[0081] From the above formula, it can be seen that the larger denoising order k eigenvalues ​​contain most of the effective signal and a small part of the noise signal, and the smaller dk eigenvalues ​​behind are mainly noise. Therefore, the value of k is extremely important for signal denoising and will directly affect the denoising performance of the algorithm. This application determines the denoising order k by calculating the energy contribution rate η.

[0082] Often the noise is randomly distributed on the singular values. To reduce its influence on the denoising order k, all the modified singular values ​​are determined by referring to formula (8):

[0083]

[0084] Then the corrected singular values ​​are superimposed to calculate the energy contribution rate η:

[0085]

[0086] In formula (12), the cumulative energy contribution rate increases with the increase of the noise reduction order k. When it exceeds a certain threshold, the corresponding k is the noise reduction order. Because the threshold has a greater impact on the signal, if it is set to a fixed threshold, a large amount of noise will be retained when the signal-to-noise ratio of the signal is low, and it is easy to eliminate the effective signal when the signal-to-noise ratio of the signal is high. In order to adapt to the change of the signal-to-noise ratio, an adaptive threshold adjustment mechanism is introduced. Considering that before the earthquake arrives, the signal is mainly composed of noise components. Therefore, this application uses the data before the earthquake arrives to estimate the noise variance. First, roughly estimate the earthquake arrival time t, then estimate the noise variance of the data before the earthquake arrives, and then set the threshold of the energy contribution rate.

[0087]

[0088] thr=1-a·σ (14)

[0089] Where a is the adjustment factor. When the contribution rate η increases to exceed the threshold, the noise reduction order k can be obtained.

[0090] S140, determining the signal component based on the noise reduction order and the signal component calculation formula, and determining the kurtosis value corresponding to the signal component based on the signal component and the kurtosis value calculation formula.

[0091] Specifically, the coefficient matrix S is constructed using the eigenvector matrix Q corresponding to the first k eigenvalues ​​and the original trajectory matrix X:

[0092]

[0093] Furthermore, the reconstruction matrix Z is obtained through the eigenvector matrix Q and the coefficient matrix S:

[0094] Z i =Q i S i (16)

[0095] In the formula, Z i (i=1,2,…,k) represents a single-component reconstruction matrix, the matrix size is m*d, and the elements in the matrix are defined as z ij , and then force the matrix size to be converted. The formula is as follows:

[0096]

[0097] Then, the reconstruction matrix can be transformed into a one-dimensional time series with a length of n by the diagonal averaging method, and then k groups of initial single-component signals can be obtained, namely y 1 ,y 2 ,…,y k The formula for diagonal averaging is as follows:

[0098]

[0099] However, after analyzing the obtained k groups of signal components, it was found that some strong noise still existed. Therefore, this application uses the kurtosis value to screen the decomposed signal components, aiming to remove some signal components with strong amplitude and strong Gaussianity. The kurtosis value is a statistic that describes the distribution morphology of the signal, which can quantitatively evaluate the Gaussianity of the time series. The kurtosis value calculation formula is as follows:

[0100]

[0101] In the formula, μ 4 represents the fourth-order center distance, σ represents the noise standard deviation, x represents the time series corresponding to each signal component after singular spectrum decomposition, and kurt represents the kurtosis value.

[0102] S150, based on the kurtosis value and the screening formula, determining preliminary denoised data, processing the preliminary denoised data by complex wavelet block thresholding, and determining denoised data.

[0103] Specifically, the screening formula is:

[0104]

[0105] In the formula, a represents the confidence percentage value of the adjustment factor, which is generally taken as 90%. If the kurtosis value calculation result of the signal is within this interval, it means that the signal presents Gaussianity. Therefore, if x is mainly noise (mainly Gaussian noise), kurt satisfies formula (20) and approaches 0. If x is mainly useful signal, kurt will increase and does not satisfy formula (20). By screening the kurtosis value, some signal components dominated by noise can be removed to determine the preliminary noise reduction data.

[0106] Furthermore, the preliminary denoised data is processed by complex wavelet block thresholding (BT-DTCWT) to determine the denoised data.

[0107] At present, there are many methods for calculating wavelet thresholds, such as: maximum and minimum thresholds, unbiased risk estimation thresholds, heuristic thresholds, etc. At the same time, there are also three methods for estimating thresholds: (1) independent of noise variance, using a fixed threshold; (2) using the noise variance at the first scale; (3) using the noise variance at different scales. The core of this algorithm is to distinguish between signal and noise by setting a suitable threshold and remove the noise component from the signal. However, this threshold processing method compares the threshold of each wavelet coefficient without considering the influence of the surrounding wavelet coefficients on it, which leads to the easy appearance of isolated wavelet coefficients and the elimination of effective information near strong wavelet coefficients.

[0108] The block threshold is an improved wavelet threshold denoising technique. Instead of processing each wavelet coefficient one by one, it groups all the wavelet coefficients at a single scale into wavelet blocks and processes all the coefficients in the block uniformly. The advantages of this method are:

[0109] (1) Reduce isolated coefficients. By processing the wavelet coefficients in a region at the same time, the occurrence of isolated wavelet coefficients can be reduced, thus avoiding unnecessary signal distortion during the denoising process;

[0110] (2) Retain valid signals. The block threshold method can better retain the local characteristics of the signal because it takes into account the mutual influence between adjacent wavelet coefficients and reduces the possibility of valid signals being mistakenly eliminated.

[0111] (3) Adaptive estimation. By adjusting the block size and threshold level, the block threshold method can achieve adaptive parameter estimation for different signals.

[0112] The threshold processing method adopted in this application is a wavelet estimator based on block threshold, which can adaptively select the block size and threshold level at a single scale through a data-driven method. This data-driven method obtains the optimal block size and target threshold by presetting the evaluation model, that is, solving the minimum Stein's unbiased risk estimate, and allows them to have different values ​​at different scales. Compared with the block threshold estimator with a fixed block size, this method has significant advantages.

[0113] In one possible implementation, preliminary denoising data is processed through complex wavelet block thresholding (BT-DTCWT) to determine denoised data, including: determining multi-scale wavelet coefficients through dual-tree complex wavelet transform processing based on the preliminary denoising data; determining a valid signal through a screening formula based on the multi-scale wavelet coefficients; determining an optimal block size and a target threshold of the valid signal through a preset evaluation model based on the valid signal; and determining denoised data based on the optimal block size and the target threshold, the multi-scale wavelet coefficients and a preset wavelet estimator.

[0114] Specifically, based on the preliminary denoised data, the preliminary denoised data is decomposed into multi-scale wavelet coefficients through dual-tree complex wavelet transform (DTCWT), and the real and imaginary parts of the wavelet coefficients are threshold denoised respectively. The wavelet coefficients of the noisy signal can be expressed as:

[0115]

[0116] In the formula, represents the wavelet coefficients of the noisy signal, represents the wavelet coefficients of the ideal signal, represents the wavelet coefficient of the noise term, and σ represents the noise standard deviation of the noisy signal.

[0117] Furthermore, the kurtosis value is introduced to screen the wavelet coefficients of each scale. When the kurtosis value of the wavelet coefficients at a certain scale satisfies equation (20), it is determined to be the noise-dominated component, and the wavelet coefficients of the entire scale are set to zero. The advantage of this is that the wavelet transform itself has a high frequency resolution for the low-frequency end of the signal, making it easier to remove low-frequency noise while retaining the effective signal.

[0118] The goal of the block threshold based wavelet estimator is to estimate the wavelet coefficients of the ideal signal from the wavelet coefficients of the measured data while calculating the minimum mean square error of the estimator by adjusting the parameters:

[0119]

[0120] In the formula, It represents the wavelet coefficients after being processed by the wavelet estimator, that is, the denoised data.

[0121] Under scale a, assuming that the length of each wavelet block is L ≥ 1, there are a total of m = n / L non-overlapping wavelet blocks. If m contains a decimal, it is rounded up and zero padding is used for wavelet blocks with insufficient data.

[0122] They represent the measured value, ideal signal, and noise coefficients of the b-th wavelet block respectively.

[0123] This application uses a wavelet estimator based on the James-Stein shrinkage criterion, that is, a preset wavelet estimator, to estimate the wavelet coefficients of the signal on each wavelet block:

[0124]

[0125] in, Represents the energy of each wavelet block, which is the sum of the squares of the wavelet coefficients in the bth wavelet block:

[0126]

[0127] The noise standard deviation σ is calculated by the following formula:

[0128]

[0129] It can be seen from formula (23) that, in addition to the noise standard deviation σ, the selection of block size L and threshold level λ largely determines the performance of the final preset wavelet estimator. Therefore, this application first sets a grid search range, uses Stein's unbiased risk estimation, that is, a preset evaluation model, to calculate the minimum mean square error, and finally obtains the optimal block size L and threshold λ by solving the minimum Stein's unbiased risk estimation. Assume that the signal estimator is:

[0130]

[0131] Substituting formula (27) into formula (23), we can obtain:

[0132]

[0133] where g is R L -R L The function of 2 represents the noise variance. According to Stein's theory, g is differentiable. Therefore, the minimum mean square error of the wavelet block can be calculated by the following formula:

[0134]

[0135] Substitute equation (29) into equation (30):

[0136]

[0137] The overall minimum mean square error calculation can be expressed as the superposition of the minimum mean square errors of each wavelet block:

[0138]

[0139] Add the search range of L and λ, and obtain the optimal block size L at this scale by solving the minimized Stein's unbiased risk estimate * and target threshold λ * :

[0140]

[0141] Among them, σ represents the standard deviation of noise, L represents the block size, λ represents the threshold level, Wy represents the wavelet coefficient of the noisy signal, N represents the data length, λ F =2L*ln N.

[0142] The preliminary denoised data and denoised data were evaluated through simulation analysis, and the specific results are as follows:

[0143] (1) Evaluate preliminary noise reduction data

[0144] In order to effectively evaluate the performance of the preliminary denoised data, quantitative evaluation is required through objective evaluation indicators. This application uses signal-to-noise ratio, root mean square error, and eigenvalue peak as performance evaluation indicators of the denoising algorithm.

[0145] The signal-to-noise ratio is the ratio of signal power to noise power, expressed in decibels (dB). The calculation formula is as follows:

[0146]

[0147] In the formula, x0 represents the original signal, x represents the signal after noise reduction, and n represents the signal length. The higher the signal-to-noise ratio, the better the performance of the noise reduction algorithm.

[0148] The smaller the RMS error is, the better the signal denoising effect is. Its mathematical expression is as follows:

[0149]

[0150] In addition, different noise reduction algorithms will produce a certain degree of loss after processing the signal. This application uses the characteristic peak value of the signal in the time domain to measure the degree of signal loss. The calculation formula is as follows:

[0151] y max =max(y i )(39)

[0152] When conducting simulation experiments, we first superimpose a magnitude 4 earthquake signal with a sampling frequency of 200 Hz and a high signal-to-noise ratio with real noise at different amplitude ratios to simulate noisy seismic signals with different signal-to-noise ratios, thereby verifying the effectiveness of the denoising algorithm on real seismic signals.

[0153] Figure 2 The time domain and frequency domain diagrams of the seismic signal and the noisy signal with a noise variance of 0.5 provided in this application are as follows: Figure 2 (a) and (b) show the time domain and frequency domain diagrams of a magnitude 4 earthquake signal with a very high signal-to-noise ratio after standard deviation normalization. Figure 2 In (a), the earthquake's P-wave and S-wave signals as well as the arrival time of the P-wave can be clearly observed. Figure 2 From (b), we can see that the frequency band of the seismic signal is in the low-frequency region and is continuous and has a clustering characteristic.

[0154] By superimposing noise with a mean of 0 and a variance of 0.5 on the seismic signal, a noisy seismic signal is obtained, whose time domain and frequency domain are as follows: Figure 2 As shown in (c) and (d), the noise almost drowns the P-wave signal, making it impossible to effectively detect the arrival of an earthquake using methods such as STA / LTA. By observing the frequency domain of the high signal-to-noise ratio signal and the noisy signal, it can be found that the noise is randomly distributed in the entire sampling frequency band, and there is a higher noise in the low frequency band.

[0155] When using the SSA algorithm to denoise the seismic signal, in order to improve the stability of the algorithm and prevent interference from high-frequency strong noise, the present application first pre-processes the signal through a Butterworth low-pass filter to remove high-frequency noise above 40 Hz.

[0156] For the denoising order, the energy contribution rate is chosen instead of the singular difference spectrum to determine the denoising order because the singular difference spectrum determines the denoising order by selecting the position of the maximum value after the singular value difference. This method only retains the first few larger singular values, which may easily lead to some valid signals being discarded as noise. By using the energy contribution rate, the required denoising order can be selected by controlling the threshold size, thereby retaining the valid signal.

[0157] However, the threshold setting of energy contribution rate is also extremely important and has a great impact on the denoising result. When a fixed threshold is used to determine the denoising order, the algorithm's anti-interference ability will be poor and it will not be able to adapt to data with different signal-to-noise ratios.

[0158] In order to improve the anti-interference ability of the algorithm, this application uses the prior information of the data to adjust the threshold of the energy contribution rate through the noise variance to achieve adaptive threshold adjustment. When the signal-to-noise ratio of the noisy signal is low, the energy proportion of the effective signal is reduced, but the noise variance will increase relatively, and the corresponding threshold will decrease, thereby eliminating more noise components; when the signal-to-noise ratio of the noisy signal is high, the energy proportion of the signal is high, and the noise variance of the data will decrease, and the corresponding threshold will increase, so that more useful components can be retained.

[0159] In order to remove noise as much as possible, the Monte Carlo algorithm is used here to take the signal-to-noise ratio after denoising as the evaluation index. The signal-to-noise ratio of 195 seismic signals with high signal-to-noise ratio under different noise variances is counted. The change of signal-to-noise ratio with the adjustment factor is finally obtained, and the adjustment factor a in formula (14) is 0.2.

[0160] The results of denoising using fixed threshold and adaptive threshold are shown in Table 1. Here, the fixed threshold is set to 95%. When the signal-to-noise ratio is high, the two methods are not much different. However, as the noise variance increases, the signal-to-noise ratio of the signal decreases, and the performance of the adaptive threshold is better.

[0161] Table 1 Signal-to-noise ratio (dB) of fixed threshold and adaptive threshold denoising when input noise variance is different

[0162]

[0163] Figure 3 The SSA decomposition result without adding kurtosis value provided in this application, Figure 3 (a) is the time domain, (b) is the frequency domain, Figure 3 To filter the signal using only the energy contribution rate. The signal is decomposed as Figure 3 As shown, for the convenience of observation, the correlation coefficients of the decomposed signal components are calculated and the components with similarity exceeding 90% are merged. Figure 3 In the above figure, SSC represents the retained singular spectral component, and Res represents the noise to be removed. It can be seen that although this method removes most of the noise, it retains the SSC5 component in the low-frequency band, which is mainly low-frequency noise. This is because the SSA algorithm distinguishes between effective signals and noise by the size of the corresponding singular value of the signal. However, when the noise in a certain frequency band has a strong amplitude, the corresponding singular value will also be large, and the SSA algorithm will also retain it.

[0164] against Figure 3 The low-frequency noise components with strong amplitude and the possible high-frequency noise components that appear in the image are calculated according to formula (19), and the threshold interval is set by formula (20), thereby eliminating the components dominated by Gaussian noise. The SSA decomposition result when adding the kurtosis value is as follows: Figure 4 As shown, Figure 4 (a) is the time domain, (b) is the frequency domain, Res contains Figure 3 The low-frequency noise SSC5 is not removed in the image, which means that the algorithm successfully removes the Gaussian noise with strong amplitude and improves the signal quality.

[0165] After data preprocessing, this method is compared with the variational mode decomposition (VMD) method. Figure 5 As shown, Figure 5 (a) is the time domain, and (b) is the frequency domain. Variational mode decomposition is also a signal decomposition method that can effectively extract the dominant frequency band of the signal. It is worth noting that the denoising effect of variational mode decomposition is greatly affected by parameters. Here, the modal component is set to 5, and the center frequency is randomly initialized. After the decomposition is completed, the kurtosis value is also used to screen the signal components.

[0166] The results show that both methods can achieve adaptive bandpass filtering of seismic signals and remove high-frequency and low-frequency noise. However, in order to prevent frequency aliasing, variational mode decomposition uses an iterative method to continuously update the modal components of the signal until the tolerance of the convergence criterion set in advance is reached, which leads to its low computational efficiency. For the same data, SSA only takes 0.18s to complete the calculation, while VMD takes 4.68s, which seriously affects the processing efficiency of the denoising algorithm. In addition, the embedding dimension of SSA is very high, resulting in the number of singular spectral components being much larger than the modal components of VMD, which makes the resolution of SSA for signal decomposition higher than VMD, and can achieve more detailed signal processing. Therefore, compared with VMD, the SSA algorithm has better denoising effect and higher computational efficiency.

[0167] In addition, in order to reduce the noise level, VMD performs Wiener filtering on each modal component when processing data. The seismic signal has a wide dominant frequency band, which results in a tendency to retain high-amplitude frequency components when decomposing the signal using VMD, while the signals around this frequency are attenuated to varying degrees, reducing the quality of the denoised signal. Figure 6 is the residual comparison of different denoising methods, Figure 6 (a) is the time domain, and (b) is the residual term after the noisy signal is removed after VMD and SSA processing in the frequency domain. From the time domain, we can see that the residual after VMD denoising contains some seismic signals, which indicates that VMD may also remove some useful signals while removing noise. Further observation of the spectrum corresponding to the residual shows that the residual corresponding to the VMD algorithm has some signals between 5Hz-20Hz, and this frequency band is the dominant frequency band of the seismic signal. Combined with the frequency domain of the original signal Figure 2 (b), it can be seen that VMD easily eliminates valid signals.

[0168] In order to test the robustness of the denoising algorithm, noise with different variances is input into the signal with high signal-to-noise ratio, noisy data with different signal-to-noise ratios are simulated, and the performance of the SSA algorithm and the VMD algorithm are compared and analyzed. The results are shown in Table 2. The performance of SSA denoising is better than that of VMD overall, especially when the data signal-to-noise ratio is high.

[0169] Table 2 Signal-to-noise ratio (dB) after VMD and SSA denoising under different noise variances

[0170]

[0171] Although VMD and SSA can both perform signal denoising, they have certain limitations when dealing with noise in the same frequency band. This is because these two methods are essentially based on the frequency characteristics of the signal for decomposition and reconstruction, and when separating signals and noise, they mainly rely on the frequency distribution and variation of the signal.

[0172] In actual seismic signal processing, the composition of noise is relatively complex, including environmental noise, instrument noise, etc. These noises may be in the same frequency band as the signal and are randomly distributed. When the signal and noise overlap in frequency, traditional frequency domain-based denoising methods are difficult to effectively separate, resulting in poor final denoising effect.

[0173] Figure 7 are the time domain and time-frequency domain diagrams of different signals, such as Figure 7 As shown in the figure, the time-frequency analysis of the original signal (seismic signal), noisy signal and SSA denoised signal by CWT reveals the changing characteristics of the signal at different times and frequencies, which can provide a deeper understanding of the characteristics of the signal and the performance of the denoising algorithm.

[0174] The time-frequency diagram of the original signal is as follows Figure 7 As shown in (b), we know that seismic signals are non-stationary signals, and their energy is sparse and concentrated in a limited frequency band within a period of time. The characteristics of this signal require the denoising algorithm to accurately locate the energy concentration area of ​​the signal in the time-frequency domain and effectively separate the effective signal from the noise. Figure 7 (d) shows that the random noise is distributed throughout the entire area, especially in the large-scale (low-frequency) area where the amplitude is very high, and part of the noise is also mixed with the seismic signal.

[0175] Perform time-frequency analysis on the signal after SSA denoising, such as Figure 7 (f). This algorithm basically removes all high-frequency noise and most of the low-frequency noise. However, when facing the same-band noise, SSA will not work. In order to further improve the denoising effect, it is necessary to use a method based on the time-frequency domain for denoising.

[0176] (2) Evaluate the denoised data

[0177] The generation of simulation signals is the same as Section 3.3.3. Considering the sampling frequency and frequency band of the seismic signal, the number of wavelet decomposition layers is determined in the following way:

[0178]

[0179] Where Fs represents the sampling frequency, f 0 Represents the high frequency cutoff frequency of the approximate coefficient. The sampling frequency Fs of the data used in this application is 200Hz. Let the frequency range of the approximate coefficient in the wavelet decomposition be 0-0.5Hz. The number of wavelet decomposition layers is calculated by formula (3-40) to be 7 layers. In this way, the low frequency end of the seismic signal can be fully decomposed, which is convenient for removing low frequency noise in the data. Then set the DTCWT filter to "nearsym5_7" and the filter length to 16.

[0180] Before performing wavelet threshold denoising, it is necessary to make a rough estimate of the arrival time of the seismic signal to facilitate the subsequent solution of the noise variance of the wavelet coefficients at the corresponding scale. Figure 8 is the signal amplitude after smoothing and threshold processing, Figure 8 (a) is smoothing processing, Figure 8 (b) is threshold processing, such as Figure 8 As shown in (a), the absolute value of the seismic signal is first taken and smoothed, and then the data with amplitude less than 0.2 times the maximum amplitude is set to zero, and the Figure 8 (b) Assuming that the time point when the first amplitude is not zero is the time of the earthquake, the data before this point are all noise data, from which the noise variance of the data can be roughly calculated.

[0181] Then, DTCWT is used to perform wavelet decomposition on the data, and the real and imaginary parts of the wavelet coefficients are subjected to wavelet block threshold denoising respectively. In order to eliminate the noise-dominated components, the kurtosis value is introduced according to formula (19) to calculate the Gaussianity of the wavelet coefficients at each scale, and the threshold is calculated according to formula (20). If the kurtosis value of a wavelet scale is lower than the threshold, it means that the signal at this scale is dominated by noise, and its wavelet coefficient is set to zero. If the kurtosis value of a wavelet scale is higher than the threshold, it means that the wavelet coefficient at this scale is dominated by the signal, and block threshold denoising is performed on it.

[0182] The comparison of denoising results with and without adding kurtosis value is shown below Fig. 9 As shown, Fig. 9 (a) and (b) are the BT-DTCWT denoising results without adding kurtosis value screening. Although it effectively removes the noise in the same frequency band as the signal, the low-frequency noise cannot be completely removed. Fig. 9 (c) and (d) show that low-frequency noise components can be effectively removed by adding kurtosis value screening.

[0183] However, the problem with this method is that the frequency resolution of the wavelet decomposition in the high-frequency part is low, and the seismic signal is relatively sparse. Noise items are easily retained in the denoised signal, especially high-frequency noise in the same time period as the seismic signal. This noise will affect the amplitude of the signal in the time domain.

[0184] Therefore, this application finally uses SSA and BT-DTCWT for joint denoising. First, adaptive bandpass filtering is performed through the improved SSA. Since DTCWT has a high frequency resolution at the low-frequency end, the denoising effect of BT-DTCWT at the low-frequency end is better than that of SSA. Therefore, when using the improved SSA denoising, the low-frequency noise components with a center frequency of less than 10Hz are retained, and only the high-frequency noise is removed. BT-DTCWT is then used to remove low-frequency noise and noise in the same frequency band. The signal after SSA+BT-DTCWT joint denoising is as follows. Fig.10 As shown, Fig.10(a) is the time domain, Fig.10 (b) Time-frequency domain. It can be seen from (b) that this method can effectively remove high-frequency noise and low-frequency noise, and can also remove noise in the same frequency band, greatly improving the quality of the signal.

[0185] Finally, the denoising method of this application is compared and analyzed with other methods. Under different input noise variances, this application compares the performance in terms of signal-to-noise ratio, root mean square error, peak absolute error (PAE), etc. Fig.11 Performance of various denoising methods under different input noise variances.

[0186] exist Fig.11 As can be seen from (a), both BT-DTCWT and SSA+BT-DTCWT achieve good results. Compared with Sure-DWT, a commonly used wavelet threshold denoising algorithm, the signal-to-noise ratio is improved by 4 dB. Fig.11 In (b), BT-DTCWT and SSA+BT-DTCWT have smaller RMS errors and perform better than other algorithms. Fig.11 As shown in (c), when calculating the PAE of a signal, facing signals with different input noise variances, SSA+BT-DTCWT can achieve a smaller and more stable error, and its performance is better than that of the BT-DTCWT algorithm alone.

[0187] The solution provided by this application is: first, the singular spectrum analysis method is improved. In view of the problem of empirical selection of the embedded dimension in the original algorithm, PSD is used for adaptive adjustment; for the selection of the denoising order k, the cumulative energy contribution rate of the signal is first calculated, and the denoising order is determined by using an adaptive threshold based on the noise variance, and then the signal components are screened by the kurtosis value to remove the components dominated by Gaussian noise to obtain the final denoised signal. In view of the problem that the improved singular spectrum analysis cannot suppress the same-band noise, an algorithm based on the complex wavelet block threshold is used to remove the same-band noise. First, DTCWT is used for wavelet decomposition, and then the kurtosis value is used to screen the signal to remove the components dominated by Gaussian noise. In the signal-dominated components, the optimal block size L and threshold λ are obtained by calculating the minimum stein's unbiased risk estimate, and then the wavelet coefficients are processed by a wavelet estimator based on the James-Stein shrinkage criterion. The denoised signal is obtained after inverse transformation, which is efficient and accurate in denoising the seismic signal compared with the prior art.

[0188] Fig.12 This is a structural block diagram of a seismic signal denoising device provided in one embodiment of the present application. For the sake of convenience, only the parts related to the embodiment of the present application are shown. Fig.12The seismic signal denoising device 1200 includes an acquisition module 1201 and a determination module 1202 .

[0189] In one implementation, the apparatus 500 may be used to implement the above Figure 1 For example, the acquisition module 501 is used to implement S110, and the determination module 1202 is used to implement S120 and S150.

[0190] The solution provided by this application is: first, the singular spectrum analysis method is improved. In view of the problem of empirical selection of the embedded dimension in the original algorithm, PSD is used for adaptive adjustment; for the selection of the denoising order k, the cumulative energy contribution rate of the signal is first calculated, and the denoising order is determined by using an adaptive threshold based on the noise variance, and then the signal components are screened by the kurtosis value to remove the components dominated by Gaussian noise to obtain the final denoised signal. In view of the problem that the improved singular spectrum analysis cannot suppress the same-band noise, an algorithm based on the complex wavelet block threshold is used to remove the same-band noise. First, DTCWT is used for wavelet decomposition, and then the kurtosis value is used to screen the signal to remove the components dominated by Gaussian noise. In the signal-dominated components, the optimal block size L and threshold λ are obtained by calculating the minimum stein's unbiased risk estimate, and then the wavelet coefficients are processed by a wavelet estimator based on the James-Stein shrinkage criterion. The denoised signal is obtained after inverse transformation, which is efficient and accurate in denoising the seismic signal compared with the prior art.

[0191] Fig.13 This is a schematic diagram of the structure of an electronic device provided by an embodiment of the present application. Fig.13 As shown, the electronic device 13 of this embodiment includes: at least one processor 130 ( Fig.13 Only one is shown in the figure) a processor, a memory 131, and a computer program 132 stored in the memory 131 and executable on the at least one processor 130, and when the processor 130 executes the computer program 132, the steps in any of the above method embodiments are implemented.

[0192] The electronic device 13 may be a computing device such as a desktop computer, a notebook, a PDA, or a cloud server. The electronic device may include, but is not limited to, a processor 130 and a memory 131. Those skilled in the art will appreciate that Fig.13 It is only an example of the electronic device 13 and does not constitute a limitation on the electronic device 13. It may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, it may also include input and output devices, network access devices, etc.

[0193] The processor 130 may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor may be a microprocessor or any conventional processor, etc.

[0194] In some embodiments, the memory 131 may be an internal storage unit of the electronic device 13, such as a hard disk or memory of the electronic device 13. In other embodiments, the memory 131 may also be an external storage device of the electronic device 13, such as a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card, etc. equipped on the electronic device 13. Further, the memory 131 may also include both an internal storage unit of the electronic device 13 and an external storage device. The memory 131 is used to store an operating system, an application program, a boot loader (BootLoader), data, and other programs, such as the program code of the computer program. The memory 131 may also be used to temporarily store data that has been output or is to be output.

[0195] The above description is only a preferred implementation mode of the present application and is not intended to limit the present application. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present application should be included in the protection scope of the present application.

Claims

1. A seismic signal denoising method, characterized in that: The method comprises: Acquire seismic signals; Determine a trajectory matrix corresponding to the embedding dimension based on a power spectrum density and an embedding dimension calculation formula, wherein the power spectrum density is determined according to the seismic signal; Based on the trajectory matrix, singular spectrum analysis and energy contribution rate calculation formula, an energy contribution rate is determined, and based on the energy contribution rate and an adaptive threshold adjustment mechanism, a noise reduction order is determined; Determine the signal component based on the noise reduction order and the signal component calculation formula, and determine the kurtosis value corresponding to the signal component based on the signal component and the kurtosis value calculation formula; Based on the kurtosis value and the screening formula, preliminary denoised data is determined, and the preliminary denoised data is processed by complex wavelet block threshold to determine denoised data; The step of processing the preliminary denoised data by complex wavelet block threshold to determine denoised data comprises: Based on the preliminary denoised data, determining multi-scale wavelet coefficients through dual-tree complex wavelet transform processing; Based on the multi-scale wavelet coefficients, determining a valid signal through the screening formula; Based on the valid signal, determining an optimal block size and a target threshold of the valid signal through a preset evaluation model; Determining denoised data based on the optimal block size and target threshold, the multi-scale wavelet coefficients and a preset wavelet estimator; The kurtosis value calculation formula is: ； Among them, μ4 represents the fourth-order center distance, Characterizes the noise standard deviation; x Characterizes the time series corresponding to each signal component after singular spectrum decomposition, kurt characterizes the kurtosis value, n Characterizes the signal length.

2. The method according to claim 1, characterized in that The embedding dimension calculation formula is: ； in, f max Characterizes and calculates the power spectrum density of the seismic signal to obtain the frequency corresponding to the maximum peak. F Characterizes the sampling frequency of the signal, a Characterize regulatory factors, d Represents the embedding dimension, n Characterizes the signal length.

3. The method according to claim 1, characterized in that The determining of the energy contribution rate based on the trajectory matrix, singular spectrum analysis and energy contribution rate calculation formula includes: Performing autocorrelation analysis on the trajectory matrix to obtain a covariance matrix; Performing singular spectrum decomposition on the covariance matrix to determine modified singular values ​​corresponding to the covariance matrix; The energy contribution rate is determined based on the modified singular value through the energy contribution rate calculation formula.

4. The method according to claim 3, characterized in that The energy contribution rate calculation formula is: ； Among them, η represents the energy contribution rate, represents the modified singular value corresponding to the covariance matrix, k represents the denoising order, and d represents the embedding dimension.

5. The method according to claim 1, characterized in that The screening formula is: ； in, a Characterize regulatory factors, 24 / n Characterize the variance, n Characterizes the signal length.

6. A seismic signal denoising device, characterized in that: The device comprises: An acquisition module, used for acquiring seismic signals; A determination module, used to determine a trajectory matrix corresponding to the embedding dimension based on a power spectrum density and an embedding dimension calculation formula, wherein the power spectrum density is determined according to the seismic signal; The determination module is further used to determine the energy contribution rate based on the trajectory matrix, singular spectrum analysis and energy contribution rate calculation formula, and determine the noise reduction order based on the energy contribution rate and the adaptive threshold adjustment mechanism; The determination module is further used to determine the signal component based on the noise reduction order and the signal component calculation formula, and determine the kurtosis value corresponding to the signal component based on the signal component and the kurtosis value calculation formula; The determination module is further used to determine preliminary denoised data based on the kurtosis value and the screening formula, and to process the preliminary denoised data through complex wavelet block threshold to determine denoised data; The step of processing the preliminary denoised data by complex wavelet block threshold to determine denoised data comprises: Based on the preliminary denoised data, determining multi-scale wavelet coefficients through dual-tree complex wavelet transform processing; Based on the multi-scale wavelet coefficients, determining a valid signal through the screening formula; Based on the valid signal, determining an optimal block size and a target threshold of the valid signal through a preset evaluation model; Determining denoised data based on the optimal block size and target threshold, the multi-scale wavelet coefficients and a preset wavelet estimator; The kurtosis value calculation formula is: ； Among them, μ4 represents the fourth-order center distance, Characterizes the noise standard deviation; x Characterizes the time series corresponding to each signal component after singular spectrum decomposition, kurt characterizes the kurtosis value, n Characterizes the signal length.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the method according to any one of claims 1 to 5 is implemented.

Images