Seismic signal denoising method and device

Abstract

The present invention discloses a method and device for denoising seismic signals. The method comprises: decomposing the seismic signal into multiple finite bandwidth modes by a preset variational modal decomposition method to obtain a first group of eigenmodal components; extracting independent source signals from the first group of eigenmodal components by a preset independent component analysis method to obtain a first group of independent source components and a mixing matrix; performing noise judgment and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components; and reconstructing the second group of independent source components to obtain a denoised seismic signal. The present invention first decomposes the seismic signal into different modal components from high frequency to low frequency, then extracts the independent source signal of the seismic signal from it, and denoises the independent source signal, so that the effective signal will not be damaged while the noise is removed, thereby achieving a good separation between the effective signal and random noise, and significantly improving the signal-to-noise ratio of the seismic signal.

Description

Technical Field

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

[0002] The collection of seismic signals is affected by a variety of factors, including vibration interference, equipment failure, and data loss during transmission. These influences are reflected in the seismic signal as random noise, reducing its quality and increasing the difficulty of subsequent processing and analysis. Currently, the main methods for suppressing random noise in seismic signals include those based on filtering theory, wavelet domain transforms, matrix theory, and signal decomposition theory. All of these methods can denoise seismic signals to a certain extent. Summary of the Invention

[0003] The inventors have discovered that existing methods have achieved certain results in seismic signal denoising, but are still imperfect and have room for improvement. By improving existing denoising methods, better denoising effects can be achieved. To at least partially address the technical problems of the prior art, the inventors have devised the present invention, which provides the following technical solutions through specific implementation methods:

[0004] In a first aspect, an embodiment of the present invention provides a method for denoising a seismic signal, the method comprising the following steps:

[0005] The seismic signal is decomposed into multiple finite bandwidth modes using a preset variational mode decomposition method and converted into a time domain signal to obtain the first set of eigenmode components;

[0006] Extracting independent source signals from the first group of intrinsic mode components using a preset independent component analysis method to obtain a separation matrix of the seismic signal, and obtaining a first group of independent source components and a mixing matrix of the seismic signal based on the separation matrix;

[0007] performing noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components;

[0008] The second group of independent source components is multiplied by the mixing matrix to construct a second group of eigenmode components; and the second group of eigenmode components is summed to obtain a denoised seismic signal.

[0009] Furthermore, performing noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components specifically includes:

[0010] For each independent source component in the first set of independent source components:

[0011] The time domain signal of the target independent source component is uniformly sampled at a preset time interval to obtain a certain number of sampling point signal values;

[0012] Determining a noise threshold of the target independent source component according to the signal values ​​of the certain number of sampling points;

[0013] Comparing the signal values ​​of the certain number of sampling points with the noise threshold respectively, determining whether the sampling point signal corresponding to each sampling point signal value is noise, and performing amplitude processing on the signal of each sampling point according to the determination result to obtain a quantized signal value of each sampling point;

[0014] A new independent source component of the target independent source component is determined according to the obtained quantized signal values ​​of each sampling point.

[0015] Furthermore, the step of comparing the signal values ​​of the certain number of sampling points with the noise threshold, determining whether the sampling point signal corresponding to each sampling point signal value is noise, and performing amplitude processing on the signal of each sampling point according to the determination result to obtain the quantized signal value of each sampling point includes:

[0016] The signal values ​​of the certain number of sampling points are compared with the noise threshold value respectively by the following formula (1), and it is judged whether the sampling point signal corresponding to each sampling point signal value is noise. Then, the amplitude of each sampling point signal is processed according to the judgment result to obtain the quantized signal value of each sampling point:

[0017]

[0018] Among them, y′ i represents the i-th independent source component in the first group of independent source components, y″ i represents y′ i The new independent source component after noise judgment and amplitude processing, j represents the jth sampling point in the preset number of sampling points of the independent source component, y′ i (j) represents the signal value of the jth sampling point of the i-th independent source component, t i is the threshold of the i-th independent source component in the first group of independent source components, a is a given minimum value used to suppress the amplitude of the noise signal, and sign is the sign function.

[0019] Furthermore, determining the noise threshold of the target independent source component based on the signal values ​​of the certain number of sampling points includes:

[0020] Substitute the signal values ​​of the certain number of sampling points into the following formula (2) to obtain the noise threshold of the target independent source component:

[0021]

[0022] Among them, t i is the threshold of the ith independent source component in the first group of independent source components, σ i represents the mean square error of the i-th independent source component, N refers to the component length, that is, the number of preset sampling points, median is a function for finding the median, abs is a function for finding the absolute value, y′ i represents the i-th independent source component in the first group of independent source components, y′ i (j) represents the signal value of the jth sampling point among the N sampling points of the i-th independent source component.

[0023] Furthermore, the seismic signal is decomposed into multiple finite bandwidth modes by a preset variational modal decomposition method and converted into a time domain signal, including:

[0024] Substituting the signal value of the seismic signal into the following formula (3) to obtain K finite bandwidth modes, where K is the preset number of finite bandwidth modes, and solving formula (3) by the alternating direction multiplier algorithm to obtain K eigenmodal components converted into time domain signals, the K eigenmodal components constituting the first group of eigenmodal components:

[0025]

[0026] Among them, t is the time variable of the seismic signal, u k (t) represents the kth finite bandwidth mode of the seismic signal, and k = 1, 2, ···, K, u k n+1 (t) represents the result of the n+1th iteration, arg min represents the value of the independent variable when the function takes the minimum value, is the derivative with respect to time t, δ(t) is the unit impulse function, j is the complex unit, * is the convolution operator, The meaning of is to transform each finite bandwidth mode u k (t) becomes an analytical signal, making the real-valued signal u k (t) is transformed into a complex value, represents the square of the second norm, ω k (t) is the corresponding eigenmode component u k (t) is the frequency center, and α is the penalty factor, x(t) represents K finite bandwidth modes u k (t), and λ(t) is the Lagrange multiplication operator.

[0027] Furthermore, the extracting independent source signals from the first group of intrinsic mode components by a preset independent component analysis method includes:

[0028] The first group of eigenmode components is substituted into the following formula (4), and according to the preset convergence condition, the optimal solution of formula (4) is iterated by the alternating direction multiplier algorithm to obtain all row vectors of the separation matrix. The separation matrix of the seismic signal is obtained by adding all the row vectors:

[0029]

[0030] Where z represents the set of the first group of eigenmode components, w represents a row vector of the separation matrix W, and w i and w i+1 denote the i-th and i+1-th row vectors of the separation matrix, respectively. E(·) denotes the expected value of the random variable. g(·) denotes the derivative of the function G(·). g'(·) denotes the derivative of g(·). The value of the function G(·) satisfies the following formula. T denotes the transpose of the vector. represents the square of the two-norm;

[0031]

[0032] Here, x represents a random variable.

[0033] Furthermore, the preset convergence condition includes: when the following formula is satisfied or when the number of iterations reaches a preset number of iterations, the iteration is terminated;

[0034] ||w i+1 -w i ||<ε ε>0

[0035] Among them, ε is the preset accuracy threshold.

[0036] Furthermore, before extracting independent source signals from the first group of intrinsic mode components using a preset independent component analysis method, the method further includes:

[0037] Centralizing the first group of eigenmode components.

[0038] Furthermore, before extracting independent source signals from the first group of intrinsic mode components using a preset independent component analysis method, the method further includes:

[0039] A whitening process is performed on the first group of eigenmode components.

[0040] In a second aspect, an embodiment of the present invention provides a seismic signal denoising device, the seismic signal denoising device comprising:

[0041] A modal decomposition module is used to decompose the seismic signal into multiple finite bandwidth modes using a preset variational modal decomposition method, and convert it into a time domain signal to obtain a first set of eigenmodal components;

[0042] an independent source extraction module, configured to extract independent source signals from the first group of intrinsic mode components using a preset independent component analysis method to obtain a separation matrix of the seismic signal, and obtain a first group of independent source components and a mixing matrix of the seismic signal based on the separation matrix;

[0043] a noise suppression module, configured to perform noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components;

[0044] A signal reconstruction module is used to multiply the second group of independent source components by the mixing matrix to construct a second group of eigenmode components; and sum the second group of eigenmode components to obtain a denoised seismic signal.

[0045] Furthermore, the noise suppression module includes:

[0046] a noise threshold calculation module, configured to uniformly sample the time domain signal of the target independent source component at a preset time interval for each independent source component in the first group of independent source components, obtain a certain number of sampling point signal values, and determine the noise threshold of the target independent source component based on the certain number of sampling point signal values;

[0047] a signal amplitude suppression module, configured to compare the signal values ​​of the certain number of sampling points with the noise threshold, determine whether the sampling point signal corresponding to each sampling point signal value is noise, perform amplitude processing on the signal of each sampling point according to the determination result, obtain a quantized signal value of each sampling point, and determine a new independent source component of the target independent source component according to the obtained quantized signal values ​​of each sampling point.

[0048] Furthermore, the seismic signal denoising device also includes a preprocessing module, which is used to perform centering and whitening on the first group of intrinsic mode components before the independent source extraction module processes the first group of intrinsic mode components, and input the first group of intrinsic mode components that have undergone the centering and whitening processing into the independent source extraction module.

[0049] In a third aspect, an embodiment of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the seismic signal denoising method as described in any of the above solutions.

[0050] The beneficial effects of the above technical solutions provided by the embodiments of the present invention include at least:

[0051] The embodiment of the present invention first decomposes the seismic signal into different bandwidth modes through a preset variational modal decomposition method, and converts it into a time domain signal to obtain a first group of eigenmodal components, so as to reduce the residual noise in each mode and reduce the redundant modes in the seismic signal; then, fixed-point independent component analysis is performed on the first group of eigenmodal components to obtain independent source components of the seismic signal, so as to facilitate noise judgment on each independent source component; amplitude processing is performed on each independent source component according to the noise judgment result to reduce the amplitude of the noise signal; finally, each independent source component that has undergone amplitude processing is reconstructed to obtain a denoised seismic signal, thereby achieving effective denoising of the seismic signal.

[0052] Other features and advantages of the present invention will be described in the following description, and in part will become apparent from the description, or will be understood by practicing the present invention. The purposes and other advantages of the present invention can be realized and obtained by the structures particularly pointed out in the written description, claims, and drawings.

[0053] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0055] Figure 1 This is a flow chart of a seismic signal denoising method in an embodiment of the present invention;

[0056] Figure 2a A preset effective seismic signal image in an embodiment of the present invention;

[0057] Figure 2b Yes Figure 2a A synthetic seismic signal image after adding noise to the seismic signal in;

[0058] Figure 2c yes Figure 2b A denoised seismic signal image is obtained after the synthetic seismic signal in the embodiment of the present invention is processed in step S1;

[0059] Figure 2d yes Figure 2b A denoised seismic signal image obtained by processing the synthetic seismic signal in the embodiment of the present invention through the complete denoising process;

[0060] Figure 3 In an embodiment of the present invention, a first group of eigenmode components is obtained by decomposing a preset synthetic seismic signal;

[0061] Figure 4 yes Figure 3 The second group of eigenmode components is obtained by extracting the first group of eigenmode components in and reconstructing them after processing the independent source signals;

[0062] Figure 5 This is an effect diagram of the preset synthetic seismic signal after denoising in an embodiment of the present invention, wherein Figure (a) is a preset effective seismic signal image, Figure (b) is a synthetic seismic signal image after adding noise to the seismic signal in Figure (a), Figure (c) is a denoised seismic signal image obtained by processing the synthetic seismic signal in Figure (b), and Figure (d) is a noise signal image obtained by processing the synthetic seismic signal in Figure (b). DETAILED DESCRIPTION

[0063] Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0064] Example 1

[0065] The embodiment of the present invention provides a method for denoising a seismic signal. Figure 1 As shown, the following steps are included:

[0066] S1. Decompose the seismic signal into multiple finite bandwidth modes using a preset variational mode decomposition method, and convert it into a time domain signal to obtain the first set of eigenmode components.

[0067] The variational modal decomposition method is a type of signal decomposition estimation method that can decompose the signal into multiple modal components with a certain frequency bandwidth from high frequency to low frequency. In this embodiment, the collected seismic signal is regarded as consisting of a series of seismic signals, and each seismic signal is defined as a seismic signal collected within a certain time T, which can be expressed as x(t). Each seismic signal x(t) has a wide frequency range, which includes effective signals and noise signals. The seismic signal x(t) can be regarded as the sum of K modal components with a certain bandwidth frequency. K is the preset number of decompositions. For example, if K=7 is set, the decomposed set of modal components with a limited bandwidth can be expressed as [u k (t)],[u k (t)]=[u1(t), u2(t),···u K (t)],u k (t) can be defined as an AM-FM signal and can be expressed as follows:

[0068]

[0069] And satisfy:

[0070]

[0071] Wherein, k represents the sequence number of each modal component in K modal components, k = 1, 2, ···, K; t represents the time t within the time T; A k (t) is the amplitude of the kth modal component at time t, and A k (t)≥0; is the phase of the kth modal component at time t, and u k (t) the instantaneous frequency ω k (t) can be expressed by the following formula:

[0072]

[0073] Each finite bandwidth modal component u of the decomposition can be obtained by the following derivation: k The expression of (t):

[0074] The frequency bandwidth of each modal component is estimated in the frequency domain. The mathematical expression is as follows:

[0075]

[0076] Among them, ω k (t) is the kth modal component u k (t) is the frequency center, and is the derivative with respect to time t, δ(t) is the unit impulse function, j is the complex unit, * is the convolution operator, Expressed as the square of the two-norm, The meaning of is to transform the modal component u through Hilbert transform k (t) becomes an analytical signal, making the real-valued signal u k (t) is transformed into a complex value.

[0077] Then, the penalty factor α and the Lagrange multiplication operator λ are introduced to transform the constrained variational model into an unconstrained variational model, and we get u k The expression of (t) is shown in formula (3):

[0078]

[0079] Among them, u k (t) represents the kth finite bandwidth mode of the seismic signal, u k n+1(t) represents the result of the n+1th iteration, arg min represents the value of the independent variable when the function takes the minimum value, ω k (t) is the corresponding finite bandwidth mode u k (t), α is the penalty factor, and λ(t) is the Lagrange multiplication operator.

[0080] At this point, we get the frequency domain solutions of K finite bandwidth modal components. By solving formula (3) using the Alternating Direction Method of Multipliers (ADMM), we can convert each finite bandwidth modal component into a time domain signal. The specific steps are as follows:

[0081] Take the inverse Fourier transform of formula (3) to obtain the modal component u k (t) and the center frequency ω of the modal component k (t) and the corresponding Lagrange multiplication operator expressions are as follows:

[0082]

[0083]

[0084]

[0085] in, for u k Fourier transform of (t), n+1 represents the n+1th iteration, ω is the frequency; The current remainder Wiener filtering; is the center of the modal frequency spectrum; τ represents the noise tolerance parameter, which needs to be determined according to the signal-to-noise ratio of the seismic signal.

[0086] Formula (3) can be solved by formulas (3-1), (3-2), and (3-3). The specific steps are as follows:

[0087] Set n = n + 1, k = k + 1, and cyclically update u according to formulas (3-1) and (3-2) k (t) and ω k (t), until k = K, the cycle ends;

[0088] Update λ(t) cyclically according to formula (3-3) and execute the above loop until formula (3-4) is satisfied:

[0089]

[0090] Wherein, ε is the preset accuracy, ε>0.

[0091] At this point, the frequency band decomposition of the seismic signal x(t) is completed, that is, K finite bandwidth modal components [u k (t)], the K finite bandwidth modal components [u k (t)] constitute the first set of eigenmode components.

[0092] S2. Extract independent source signals from the first group of intrinsic mode components using a preset independent component analysis method to obtain a separation matrix of the seismic signal, and obtain a first group of independent source components and a mixing matrix of the seismic signal based on the separation matrix.

[0093] The independent component analysis algorithm is a type of blind source separation method that searches for inherent independence and non-Gaussian factors in multidimensional data. Based on the fact that seismic signals are generated by different signal sources, this embodiment regards the seismic signal as the result of mixing a set S of n independent source signals with a mixing matrix A. Since seismic signals usually contain complex noise signals, it is not convenient to directly obtain independent source signals. In this embodiment, each seismic signal is first decomposed into K eigenmode component signals with a certain bandwidth frequency using step S1. k (t)] to reduce the residual noise in each mode and reduce the redundant modes in the seismic signal, and then the eigenmode component signal [u k (t)] to obtain independent source signal processing.

[0094] Based on the fact that seismic signals are generated by different signal sources, this embodiment converts the first group of eigenmode component signals [u k (t)] is regarded as the result of mixing the set S of n independent source component signals with the mixing matrix A, which can be expressed by formula (4-1):

[0095]

[0096] Where S represents n independent source signals [s1,s2,···,s n ], A represents a set of n matrix row vectors [a1, a2, ···, a n ] collection.

[0097] According to this characteristic, there must be a separation matrix W, W=[w1,w2,···,w n ], the seismic signal can be decomposed into a separated signal set y that is similar to the original independent source signal set S, which can be expressed by formula (4-2):

[0098] y=W*[u k (t)] (4-2)

[0099] Where y is n separated signals [y1,y2,...,yn ], and y i ≈s i ,i=1,2,···,n,W is n row vectors [w1,w2,···,w n ] collection.

[0100] By solving the first group of eigenmode component signals [u k (t)] can be used to obtain the independent source signals of the seismic signal. Before solving the separation matrix for the first group of eigenmode components, the first group of eigenmode components can also be preprocessed, including centering and whitening, to remove the correlation between the observed signals and simplify the extraction process of independent source components, making the process of solving the separation matrix simpler and the algorithm convergence faster.

[0101] In this embodiment, the preprocessed signal of the first group of eigenmode components is recorded as z, and the process of solving the independent source components is explained using the preprocessed signal z as an example.

[0102] First, it is necessary to determine the objective function for solving the separation matrix. The objective function can be obtained through the following derivation process:

[0103] The independent source component y′ of the preprocessed signal z can be expressed as y′=W T z, y′=[y′1,y′2,...,y′ n ], and y′≈y, where T represents the transpose of the vector. When solving the separation matrix W, negative entropy is selected as the measure of the non-Gaussianity of the preprocessed signal z, which is expressed as formula (4-3):

[0104]

[0105] Among them, J(·) represents the required negative entropy, z gauss is a Gaussian random vector with the same covariance matrix as the preprocessed signal z, H(·) represents the differential entropy of the random variable, and p(·) represents the probability of the random variable.

[0106] Using the fixed-point independent component analysis algorithm FastICA to solve, formula (4-3) is transformed into the following expression:

[0107] J(y′)∝[E{G(y′)}-E{G(v)}] 2 . (4-4)

[0108] Where v is a Gaussian variable with zero mean and unit variance, E(·) represents the expected value of the random variable, and G(·) is a non-quadratic function of some form. According to the Gaussian characteristics of seismic signals, the value of G(·) satisfies formula (4-5).

[0109]

[0110] Where x represents a random variable;

[0111] Combining formula (4-5) to derive formula (4-4), we can obtain the objective function for obtaining the separation matrix W, as shown in formula (4):

[0112]

[0113] Among them, w i and w i+1 denote the i-th and i+1-th row vectors of the separation matrix, respectively; g(·) denotes the derivative of the function G(·); and g'(·) denotes the derivative of g(·).

[0114] The specific process of deriving formula (4-4) includes:

[0115] In formula (4-4), in order to make the separated signal y′ satisfy the function J(W) and obtain the maximum value, combined with formula (4-5), we can know that E{G(y′) 2}=1, so the actual requirement is to take the maximum value of E{G(y′)}, and since y′=W T z, that is, to take E{G(w i T According to the Kuhn-Tucker condition, formula (4-4) can be solved by the following formula (4-6):

[0116] E{zg(w i T z)}-βw i =0 (4-6)

[0117] Among them, w i Represents a row vector of the separation matrix W, β is relative to w i The fixed value satisfies

[0118] Use Newton's method to solve w i , define an auxiliary function F(w i ), as shown in formula (4-7):

[0119]

[0120] Derivative of formula (4-7) is shown in formula (4-8):

[0121]

[0122] Where I is the unit matrix. Since the pre-processed signal z is the seismic signal after whitening and centering, E{zz T g'(w i T z)} can be approximately transformed into E{zz T}E{g'(w i T z)}I, that is, E{zz T g'(w i T z)}≈E{zz T}E{g'(w i T z)}I, so the iterative expression of Newton iteration can be obtained as:

[0123]

[0124] Known Multiply both sides of the equation (4-9) by After simplification, we get formula (4), which is the objective function for obtaining the separation matrix W.

[0125] When actually processing seismic signals, the objective function for obtaining the separation matrix has been derived. It is only necessary to input the preprocessed signal z into the function, and then, according to the set convergence conditions, the optimal solution of the objective function is iteratively obtained by the alternating direction multiplier algorithm to obtain all row vectors [w1, w2, ···, w n ], that is, the separation matrix W is obtained. The convergence conditions include: when || w i+1 -w i When ||<ε, the iteration terminates, where ε is the set accuracy threshold, usually a minimum value, and ε>0; in some cases, ||w is still not satisfied after a large number of iterations. i+1 -w i ||<ε, a preset number of iterations needs to be set at this time. When the preset number of iterations is reached, the obtained separation matrix is ​​close to the target separation matrix, and the iteration can be terminated at this time.

[0126] After obtaining the separation matrix W, the separation signals y′ can be solved according to the above formula (4-2) i , that is, [y′1,y′2,…,y′ n ], the [y′1,y′2,…,y′ n ] constitute the first group of independent source components; since the independent source signal set S is approximately the separated signal set y′, that is, y′≈S, the mixing matrix A can be obtained according to the above formula (4-1).

[0127] S3. Perform noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components.

[0128] After obtaining the first set of independent source components, denoising can be performed on each independent source component by threshold processing. The noise threshold can be used to filter out the noise signal in each independent source signal and suppress the amplitude of the noise signal, thereby obtaining the denoised seismic signal. The specific threshold processing process can be expressed by the following formulas (1) and (2):

[0129]

[0130]

[0131] Among them, y i ′ represents the i-th independent source component in the first group of independent source components, y″ i represents y′ i The new independent source component after threshold processing; t i is the threshold of the i-th independent source component; a is a given minimum value used to suppress the amplitude of the noise signal; σ i represents the mean square error of the i-th independent source component; median is the median; N refers to the component length, that is, the number of sampling points in each seismic signal; j represents the j-th sampling point among the n sampling points of each signal, y′ i (j) represents the signal value at the jth sampling point.

[0132] Formula (1) is to judge the noise of the signal value of each sampling point in each independent source component according to the threshold value when the threshold value is set. If the signal value is less than or equal to the threshold value, the signal of the sampling point is judged to be noise, and then the noise is amplitude suppressed to reduce its impact on the seismic signal; if the signal value is greater than the threshold value, it is judged to be non-noise, that is, a valid signal, and then the valid signal is amplitude contracted to make the overall display image of the seismic valid signal smoother. After processing by formula (1), a new set of independent source components is obtained, namely the second set of independent source components y″ i , y″ i =[y″1,y″2,···,y″ n ].

[0133] Formula (2) is a method for setting the noise threshold. The noise threshold is set manually, and different threshold selection methods can be used according to actual needs.

[0134] After threshold processing, the noise signals of each independent source component have been effectively removed and the effective signals have been retained. By reconstructing these independent source components, the denoised seismic signal can be obtained.

[0135] S4. Multiply the second group of independent source components by the mixing matrix to construct a second group of eigenmode components; sum the second group of eigenmode components to obtain a denoised seismic signal.

[0136] According to step S3, the second group of independent source components y″ is obtained i According to the above formula (4-1), the expression of the second group of eigenmode components corresponding to the second group of independent source components is as follows:

[0137]

[0138] The denoised seismic signal can be obtained by summing the second group of eigenmodal components.

[0139] This embodiment first decomposes the seismic signal into modes of different bandwidths using a preset variational modal decomposition method, and converts them into time-domain signals to obtain a first group of eigenmodal components, so as to reduce the residual noise in each mode and reduce redundant modes in the seismic signal; then, independent component analysis is performed on the first group of eigenmodal components to obtain independent source components of the seismic signal, so as to facilitate noise judgment on each independent source component; amplitude processing is performed on each independent source component based on the noise judgment result to reduce the amplitude of the noise signal; finally, each independent source component that has undergone amplitude processing is reconstructed to obtain a denoised seismic signal, thereby achieving effective denoising of the seismic signal.

[0140] Combined with Figure 2 to Figure 5 As shown, in this embodiment, in order to demonstrate the denoising effect, the result of denoising a preset synthetic seismic signal is taken as an example for explanation. The synthetic seismic signal can be regarded as composed of a series of synthetic seismic signals, and each of the synthetic seismic signals includes a noise signal and a preset effective seismic signal.

[0141] FIG2 is a seismic signal image of a synthetic seismic signal subjected to denoising processing using the seismic signal denoising method proposed in an embodiment of the present invention, wherein: Figure 2a It is a preset effective seismic signal image. Figure 2b Yes Figure 2a A synthetic seismic signal image after adding noise to the seismic signal in Figure 2c yes Figure 2b The de-noised seismic signal image obtained by processing the synthetic seismic signal in step S1 of the embodiment of the present invention is Figure 2d yes Figure 2b The synthetic seismic signal in is processed by the seismic signal denoising method proposed in an embodiment of the present invention to obtain a denoised seismic signal image.

[0142] Figure 3It is a first group of eigenmode components obtained by decomposing a preset synthetic seismic signal through the seismic signal denoising method proposed in an embodiment of the present invention.

[0143] Figure 4 yes Figure 3 The first group of eigenmode components in are reconstructed by extracting independent source signals and processing to obtain the second group of eigenmode components.

[0144] Figure 5 It is an effect diagram of the preset synthetic seismic signal after being processed by the seismic signal denoising method proposed in an embodiment of the present invention, wherein Figure (a) is a preset effective seismic signal image, Figure (b) is a synthetic seismic signal image after adding noise to the seismic signal in Figure (a), Figure (c) is a denoised seismic signal image obtained by processing the seismic signal in Figure (b) using the seismic signal denoising method proposed in an embodiment of the present invention, and Figure (d) is a noise signal image obtained by processing the signal in Figure (b) using the seismic signal denoising method proposed in an embodiment of the present invention.

[0145] According to Figure 2 to Figure 5 The seismic signal processing results shown in the figure show that the seismic signal denoising method proposed in the embodiment of the present invention first decomposes the seismic signal into different modal components from high frequency to low frequency, then extracts the independent source signal of the earthquake and denoises the independent source signal, so that the effective signal is not damaged while removing the noise, thereby achieving a good separation of the seismic signal and random noise. The seismic signal denoising method proposed in the embodiment of the present invention can achieve a significant noise reduction effect for the seismic signal, significantly improve the signal-to-noise ratio, and lay a good foundation for the subsequent processing and interpretation of the seismic signal.

[0146] Based on the same inventive concept, an embodiment of the present invention further provides a seismic signal denoising device, comprising:

[0147] A modal decomposition module is used to decompose the seismic signal into multiple finite bandwidth modes using a preset variational modal decomposition method, and convert it into a time domain signal to obtain a first set of eigenmodal components;

[0148] an independent source extraction module, configured to extract independent source signals from the first group of intrinsic mode components using a preset independent component analysis method to obtain a separation matrix of the seismic signal, and obtain a first group of independent source components and a mixing matrix of the seismic signal based on the separation matrix;

[0149] a noise suppression module, configured to perform noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components;

[0150] A signal reconstruction module is used to multiply the second group of independent source components by the mixing matrix to construct a second group of eigenmode components; and sum the second group of eigenmode components to obtain a denoised seismic signal.

[0151] In one embodiment, the noise suppression module includes a noise threshold calculation module and a signal amplitude suppression module, wherein:

[0152] The noise threshold calculation module is used to uniformly sample the time domain signal of the target independent source component at a preset time interval for each independent source component in the first group of independent source components, obtain a certain number of sampling point signal values, and determine the noise threshold of the target independent source component based on the certain number of sampling point signal values.

[0153] The specific noise threshold calculation method can be performed according to formula (2):

[0154]

[0155] The specific meaning of formula (2) can be seen from the implementation of the above-mentioned seismic signal denoising method.

[0156] The signal amplitude suppression module is configured to compare the signal values ​​of the certain number of sampling points with the noise threshold, determine whether the sampling point signal corresponding to each sampling point signal value is noise, perform amplitude processing on each sampling point signal based on the determination result, obtain a quantized signal value of each sampling point, and determine a new independent source component of the target independent source component based on the obtained quantized signal values ​​of each sampling point.

[0157] The specific amplitude suppression method can be implemented according to formula (1):

[0158]

[0159] The specific meaning of formula (1) can be seen from the implementation of the above-mentioned seismic signal denoising method.

[0160] In one embodiment, the seismic signal denoising device also includes a preprocessing module, which is used to perform centering and whitening on the first group of intrinsic mode components before the independent source extraction module processes the first group of intrinsic mode components, and input the first group of intrinsic mode components that have undergone the centering and whitening processing into the independent source extraction module.

[0161] Regarding the seismic signal denoising device in the above embodiment, the specific manner in which each module performs operations has been described in detail in the embodiment of the method and will not be elaborated on here.

[0162] According to an embodiment of the present invention, a non-transitory computer-readable storage medium is further provided, on which a computer program is stored. When the program is executed by a processor, the above-mentioned seismic signal denoising method can be implemented.

[0163] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Thus, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage and optical storage, etc.) containing computer-usable program code.

[0164] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowcharts and / or block diagrams, as well as combinations of processes and / or blocks in the flowcharts and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowcharts and / or block diagrams. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0165] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0166] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0167] Obviously, those skilled in the art may make various changes and modifications to the present invention without departing from the spirit and scope of the present invention. Thus, if such changes and modifications fall within the scope of the claims and their equivalents, the present invention is intended to include such changes and modifications.

Claims

1. A seismic signal denoising method, characterized in that: include: The seismic signal is decomposed into multiple finite bandwidth modes using a preset variational mode decomposition method and converted into a time domain signal to obtain the first set of eigenmode components; Extracting independent source signals from the first group of intrinsic mode components using a preset independent component analysis method to obtain a separation matrix of the seismic signal, and obtaining a first group of independent source components and a mixing matrix of the seismic signal based on the separation matrix; performing noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components; multiplying the second set of independent source components by the mixing matrix to construct a second set of eigenmode components; The second group of eigenmodal components is summed to obtain a denoised seismic signal.

2. The seismic signal denoising method according to claim 1, wherein: The performing noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components specifically includes: For each independent source component in the first set of independent source components: The time domain signal of the target independent source component is uniformly sampled at a preset time interval to obtain a certain number of sampling point signal values; Determining a noise threshold of the target independent source component according to the signal values ​​of the certain number of sampling points; Comparing the signal values ​​of the certain number of sampling points with the noise threshold respectively, determining whether the sampling point signal corresponding to each sampling point signal value is noise, and performing amplitude processing on the signal of each sampling point according to the determination result to obtain a quantized signal value of each sampling point; A new independent source component of the target independent source component is determined according to the obtained quantized signal values ​​of each sampling point.

3. The seismic signal denoising method according to claim 2, wherein: The step of comparing the signal values ​​of the certain number of sampling points with the noise threshold, determining whether the sampling point signal corresponding to each sampling point signal value is noise, and performing amplitude processing on the signal of each sampling point according to the determination result to obtain a quantized signal value of each sampling point includes: The signal values ​​of the certain number of sampling points are compared with the noise threshold value respectively by the following formula (1), and it is judged whether the sampling point signal corresponding to each sampling point signal value is noise. Then, the amplitude of each sampling point signal is processed according to the judgment result to obtain the quantized signal value of each sampling point: Among them, y′ i represents the i-th independent source component in the first group of independent source components, y″ i represents y′ i The new independent source component after noise judgment and amplitude processing, j represents the jth sampling point in the preset number of sampling points of the independent source component, y′ i (j) represents the signal value of the jth sampling point of the i-th independent source component, t i is the threshold of the i-th independent source component in the first group of independent source components, a is a given minimum value used to suppress the amplitude of the noise signal, and sign is the sign function.

4. The seismic signal denoising method according to claim 2, wherein: The determining the noise threshold of the target independent source component according to the certain number of sampling point signal values ​​includes: Substitute the signal values ​​of the certain number of sampling points into the following formula (2) to obtain the noise threshold of the target independent source component: Among them, t i is the threshold of the ith independent source component in the first group of independent source components, σ i represents the mean square error of the i-th independent source component, N refers to the component length, that is, the number of preset sampling points, median is the median function, abs is a function for finding the absolute value, y′ i represents the i-th independent source component in the first group of independent source components, y′ i (j) represents the signal value of the jth sampling point among the N sampling points of the i-th independent source component.

5. The seismic signal denoising method according to claim 1, wherein: The method of decomposing the seismic signal into multiple finite bandwidth modes by a preset variational modal decomposition method and converting them into time domain signals includes: Substituting the signal value of the seismic signal into the following formula (3) to obtain K finite bandwidth modes, where K is the preset number of finite bandwidth modes, and solving formula (3) by the alternating direction multiplier algorithm to obtain K eigenmodal components converted into time domain signals, the K eigenmodal components constituting the first group of eigenmodal components: Among them, t is the time variable of the seismic signal, u k (t) represents the kth finite bandwidth mode of the seismic signal, and k = 1, 2, ···, K, u k n+1 (t) represents the result of the n+1th iteration, arg min represents the value of the independent variable when the function takes the minimum value, is the derivative with respect to time t, δ(t) is the unit impulse function, j is the complex unit, * is the convolution operator, The meaning is to transform each eigenmode component u through Hilbert transform k (t) becomes an analytical signal, making the real-valued signal u k (t) is transformed into a complex value, represents the square of the second norm, ω k (t) is the corresponding finite bandwidth mode u k (t) is the frequency center, and α is the penalty factor, x(t) represents K finite bandwidth modes u k (t), and λ(t) is the Lagrange multiplication operator.

6. The seismic signal denoising method according to claim 1, wherein: The extracting of independent source signals from the first group of intrinsic mode components by a preset independent component analysis method includes: The first group of eigenmode components is substituted into the following formula (4), and according to the preset convergence condition, the optimal solution of formula (4) is iterated by the alternating direction multiplier algorithm to obtain all row vectors of the separation matrix. The separation matrix of the seismic signal is obtained by adding all the row vectors: Where z represents the set of the first group of eigenmode components, w represents a row vector of the separation matrix W, and w i and w i+1 denote the i-th and i+1-th row vectors of the separation matrix, respectively. E(·) denotes the expected value of the random variable. g(·) denotes the derivative of the function G(·). g'(·) denotes the derivative of g(·). The value of the function G(·) satisfies the following formula. T denotes the transpose of the vector. represents the square of the two-norm; Here, x represents a random variable.

7. The seismic signal denoising method according to claim 6, wherein: The preset convergence condition includes: when the following formula is satisfied or when the number of iterations reaches a preset number of iterations, the iteration is terminated; ||w i+1 -w i ||<e e>0 Among them, ε is the preset accuracy threshold.

8. The seismic signal denoising method according to claim 1, wherein: Before extracting independent source signals from the first group of intrinsic mode components using a preset independent component analysis method, the method further includes: Centralizing the first group of eigenmode components.

9. The seismic signal denoising method according to claim 8, wherein: Before extracting independent source signals from the first group of intrinsic mode components using a preset independent component analysis method, the method further includes: A whitening process is performed on the first group of eigenmode components.

10. A seismic signal denoising device, characterized in that: include: A modal decomposition module is used to decompose the seismic signal into multiple finite bandwidth modes using a preset variational modal decomposition method, and convert it into a time domain signal to obtain a first set of eigenmodal components; an independent source extraction module, configured to extract independent source signals from the first group of intrinsic mode components using a preset independent component analysis method to obtain a separation matrix of the seismic signal, and obtain a first group of independent source components and a mixing matrix of the seismic signal based on the separation matrix; a noise suppression module, configured to perform noise determination and amplitude processing on each independent source component in the first group of independent source components to obtain a second group of independent source components; a signal reconstruction module, configured to multiply the second group of independent source components by the mixing matrix to construct a second group of eigenmode components; The second group of eigenmodal components is summed to obtain a denoised seismic signal.

11. The seismic signal denoising device according to claim 10, wherein: The noise suppression module includes: a noise threshold calculation module, configured to uniformly sample the time domain signal of the target independent source component at a preset time interval for each independent source component in the first group of independent source components, obtain a certain number of sampling point signal values, and determine the noise threshold of the target independent source component based on the certain number of sampling point signal values; a signal amplitude suppression module, configured to compare the signal values ​​of the certain number of sampling points with the noise threshold, determine whether the sampling point signal corresponding to each sampling point signal value is noise, perform amplitude processing on the signal of each sampling point according to the determination result, obtain a quantized signal value of each sampling point, and determine a new independent source component of the target independent source component according to the obtained quantized signal values ​​of each sampling point.

12. The seismic signal denoising device according to claim 10 or 11, wherein: It also includes a preprocessing module, which is used to perform centering and whitening on the first group of intrinsic mode components before the independent source extraction module processes the first group of intrinsic mode components, and input the first group of intrinsic mode components that have undergone the centering and whitening processing into the independent source extraction module.

13. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the seismic signal denoising method according to any one of claims 1 to 9 is implemented.

Images