A signal processing method for adaptively suppressing abnormal interference of microseismic signals

Abstract

The present invention belongs to the field of microseismic signal data processing and is a signal processing method for adaptively removing abnormal interference. It combines the particle swarm optimization algorithm, the variational mode decomposition algorithm, the single-frequency interference judgment and suppression algorithm, and the dictionary learning algorithm to achieve the purpose of adaptively removing abnormal interference. The method first uses the particle swarm optimization algorithm to determine the global optimal parameter K, then decomposes the signal into K intrinsic mode functions (imfs) through the variational mode decomposition algorithm, then accurately finds the imf components with single-frequency interference through the single-frequency interference judgment algorithm, and then specifically removes the single-frequency interference. Finally, the dictionary learning algorithm is used to further suppress random interference, so as to obtain the microseismic signal after noise reduction.

Description

Technical Field

[0001] The present invention belongs to the field of microseismic signal data processing, and is a signal processing method that uses the particle swarm optimization algorithm, variational mode decomposition algorithm, single-frequency interference suppression algorithm, and dictionary learning algorithm to adaptively remove abnormal interference. Background Art

[0002] As a highly sensitive geophysical monitoring means, microseismic monitoring technology plays an important role in the fields of oil and gas field development, mine safety, geothermal resource exploration, and geological disaster warning. It realizes the dynamic perception and stability assessment of underground structures by capturing the weak vibration signals generated by underground rock mass fractures or fluid activities. However, the microseismic signals collected in the actual engineering environment are often contaminated by various complex interference sources, including mechanical vibration, electromagnetic noise, environmental background noise, and instrument self-noise. These interference signals are highly aliased with the effective microseismic signals in the time-frequency domain. Especially under low signal-to-noise ratio conditions, traditional filtering methods are difficult to effectively distinguish and suppress abnormal interference, resulting in missed detection or misjudgment of microseismic events, directly affecting the reliability of subsequent inversion and interpretation. Traditional band-stop filters are suitable for suppressing fixed-frequency interference. The Fourier transform cannot dynamically capture the interference characteristics that change with time in non-stationary signals. And common LMS or RLS algorithms may show insufficient adaptability and high computational overhead in complex environments. Therefore, it is particularly important to find a signal processing method that can adaptively suppress abnormal interference in microseismic signals. The variational mode decomposition method can perform fine modal decomposition on signals, thus effectively separating and suppressing single-frequency interference. The dictionary learning method reconstructs the core features of signals through constructing sparse representations, further enhancing the noise reduction effect. This method not only effectively suppresses single-frequency interference but also removes random noise, which helps to retain the effective components of microseismic signals. Summary of the Invention

[0003] The present invention aims to provide a method that uses the particle swarm optimization algorithm (PSO), variational mode decomposition algorithm (VMD), single-frequency interference suppression algorithm, and dictionary learning algorithm (K-SVD) to adaptively suppress single-frequency interference and remove random noise to improve the signal-to-noise ratio of microseismic signals. First, the optimal number of modes K of the variational mode decomposition (VMD) of microseismic signals is obtained through the particle swarm algorithm, so as to obtain K intrinsic mode functions (imf). According to the ratio of the area of each imf component to the maximum value threshold, single-frequency interference is found by setting a threshold, and then substitution processing is performed on the frequency domain of the single-frequency interference. Subsequently, the dictionary learning method is used to denoise random noise, thereby improving the signal-to-noise ratio of microseismic signals. The present invention includes the following steps:

[0004] S1. Import the microseismic signal x(t), use the particle swarm optimization algorithm to optimize the parameters before variational mode decomposition (VMD), and obtain the global optimal mode parameter K, where t is time;

[0005] S11. Determine the sampling frequency fs, sampling time interval dt, frequency resolution df, the number of particles SearchAgents_no participating in the optimization process, the maximum number of iterations Max_iteration, specify the parameter dimension dim to be optimized, the minimum value lb of the search space, the maximum value ub of the search space, and select the fitness function fitness;

[0006]

[0007] In the formula, K * is the mode parameter to be optimized, u i represents the i-th imf component, i ∈ {1, …, K *}, σ(u i ) is the standard deviation of the i-th imf component, U is the reconstructed signal of all mode components after variational mode decomposition (VMD), SampEn(u i , 2, 0.2σ(u i ), 1) is the sample entropy of the i-th mode component, Cov(x(t), U) is the covariance between the input signal x(t) and the reconstructed signal U, σ x(t) , σ U are the standard deviations of x(t) and U, ρ is the correlation coefficient between the input signal and the reconstructed signal, log 10 (d) is the dimension factor, and d is the metric parameter for measuring the signal dimension;

[0008] S12. Perform parameter optimization through the particle swarm optimization algorithm (PSO) to determine the global optimal mode parameter K, and use the fitness function fitness to measure the optimization effect of the current parameter K;

[0009] Or

[0010] In each iteration, the particle swarm algorithm continuously optimizes K by searching the positions of the particles where the velocity

[0011]

[0012] Among them, the velocity is given by:

[0013] When the change rate of the fitness function value is less than 1e-5 or the maximum number of iterations is reached, terminate to obtain the global optimal mode parameter K;

[0014] In the formula, argmin[·] / argmax[·] represents finding the parameter K that minimizes or maximizes the objective function. represents the position of the i-th particle in the a-th generation. is the velocity of the i-th particle in the a-th generation, w is the inertia weight, which controls the search range of the particle and affects the convergence speed and global search ability. c1 and c2 are learning factors, and r1 and r2 are random numbers (uniformly distributed between 0 and 1). is the modal parameter value of the i-th particle in the a-th generation, p i is the historical optimal solution of the i-th particle, K g represents the global optimal position. represents the updated position of the i-th particle in the next iteration. represents the updated velocity of the i-th particle in the next iteration, and a is the index of the iteration number.

[0015] S2. Substitute the global optimal modal parameter K into the variational mode decomposition (VMD) to obtain K intrinsic mode functions (imfs). The specific formula is as follows:

[0016]

[0017] where, u k is the set of intrinsic mode functions (imfs), k represents the summation index variable, u k (t) represents the k-th intrinsic mode function, ω k is its central frequency, j is the imaginary unit, δ(t) is the Dirac function, * is the convolution operator, represents the partial derivative operator with respect to time t. represents the L2 norm. represents finding the set of values that can minimize the objective function among all possible combinations of {u k} and {ω k}.

[0018] S3. Perform Fourier transform (FT) on each imf component to obtain the amplitude spectrum diagram of each component, then normalize its amplitude spectrum, and judge whether there is over-decomposition based on the normalized area.

[0019] AMimf i = FT(imf i )

[0020]

[0021] If Aimf i is greater than E totalIf it is 0.5 times of, it is considered that there is an over-decomposition situation, and steps S1 and S2 can be repeated, where E total is expressed as follows:

[0022]

[0023] In the formula, imf i represents the i-th intrinsic mode function component, and AMimf i represents the amplitude spectrum obtained after performing Fourier transform (FT) on the i-th imf i component, FT[·] represents Fourier transform, and Nimf i represents the frequency-domain normalization result of the i-th imf component, and Aimf i represents the area after frequency-domain normalization of the i-th imf component, and E total is the total energy of the entire signal x(t), f1 is the starting point of the number of points in the frequency domain, and f n is the end point of the number of points in the frequency domain, x(t) is the original time-domain signal, and max[·] represents finding the maximum value;

[0024] S4. Calculate the ratio ratio of the maximum value to the adjacent value of the amplitude spectrum in each imf frequency domain;

[0025] V max = max(A data )

[0026] I max = argmax(A data )

[0027]

[0028] ratio = max(L ratio , R ratio )

[0029] In the formula, A data represents all values in the amplitude spectrum of the frequency domain, V max is the maximum value in the amplitude spectrum corresponding to the imf frequency domain, I max is the position corresponding to the maximum value in the amplitude spectrum of the frequency domain, L ratio is the ratio of the maximum value in the amplitude spectrum of the imf frequency domain to the left adjacent value, R ratio is the ratio of the maximum value in the amplitude spectrum of the imf frequency domain to the right adjacent value, ratio is the maximum value in R ratio and L ratio , and argmax[·] represents finding the frequency point corresponding to the amplitude maximum value;

[0030] S5. Calculate the ratio between ratio and the corresponding Aimf iThe ratio area_ratio is used to determine whether there is single - frequency abnormal interference in the IMF component. The specific judgment formula is:

[0031]

[0032] If area_ratio is greater than 0.5, it is determined that there is single - frequency interference in the IMF component, and then a substitution operation is performed on the amplitude spectrum in the frequency domain of this IMF component. If there is no single - frequency interference in the IMF component, no processing is required for its frequency domain;

[0033] S6. Remove the anomaly from the IMF component with single - frequency abnormal interference. Calculate the average value of 10 point positions at adjacent positions when the amplitude is maximum, specifically as follows:

[0034]

[0035] In the formula, L is the signal length, L avg is the average value of the five adjacent points on the left of the maximum value in the amplitude spectrum of the frequency domain, R avg is the average value of the five adjacent points on the right of the maximum value in the amplitude spectrum of the frequency domain, A data [j] represents the j - th value in the amplitude spectrum of the frequency domain, A data [I max is the amplitude value after removing the single - frequency interference point. Then the updated amplitude spectrum of the IMF after removing the single - frequency abnormal interference can be expressed as AMimf n =AMimf[1,2,…,A data [I max ,…,L]. When there are multiple single - frequency abnormal interferences, repeat step S6;

[0036] S7. Superimpose the IMF frequency - domain signals after single - frequency anomaly suppression, obtain the updated frequency - domain signal, and perform the inverse Fourier transform;

[0037] y(t)=ifft[F n (f)]

[0038] In the formula, y(t) is the signal after removing the single - frequency abnormal interference, F n (f) is the updated frequency - domain signal, and ifft[·] represents the inverse Fourier transform;

[0039] S8. Perform sparse dictionary learning on the signal y(t) after removing single - frequency interference in step S7 to remove random noise, specifically:

[0040] Divide the signal y(t) into NSample overlapping blocks and construct the matrix Y:

[0041] Y=[y1(t),y2(t),…,yNsample (t) ∈ R n×NSample

[0042] where y i (t) = y[1+(i - 1):1+(i - 1)+n - 1] represents the i-th block, and the dictionary D and the sparse coefficient matrix X are solved through the optimization problem:

[0043] The constraint conditions are:

[0044] Reconstruct the denoised block using the trained dictionary and sparse coefficients:

[0045] Y d = DX ∈ R N×Sample

[0046] Overlap and add the blocks and normalize them to obtain the final signal y out (t):

[0047]

[0048] In the formula, y(t) is the input noise signal, y i (t) is the i-th block of the signal y, n is the length of the signal block y i (t), m is the number of dictionary atoms, E is the maximum allowable error of sparse representation, NSample is the number of overlapping blocks into which the signal is divided, Y is the matrix after block division, D is the trained dictionary, X is the sparse matrix, X i is the sparse coefficient vector corresponding to the i-th signal block y i , Y d is the denoised block matrix, Y d [·] represents the i-th matrix after denoising, Weight[·] is the normalization weight, y out (t) is the denoised signal, Start i represents the starting superposition position index of the i-th block in the final signal y i , R is the set of real numbers, where i represents any number.

[0049] An adaptive signal processing method for suppressing abnormal interference of microseismic signals according to the present invention has the following advantages:

[0050] (1) In the optimization parameter search step described above, the globally optimal number of modes K determined using the particle swarm optimization algorithm (PSO) will not cause incomplete decomposition or over-decomposition of the signal;

[0051] (2) In the step of determining whether there is over-decomposition or incomplete decomposition of the IMF, it is determined by setting a threshold value threshold. The specific value of the threshold value threshold is 0.5 times the total energy of the signal. This step can avoid the situations of insufficient mode, mode mixing, and pseudo-mode in the subsequent processed IMF components;

[0052] (3) In the step of removing abnormal interference, first, the single-frequency interference is adaptively found through the judgment criterion, and then the substitution processing is used to suppress the single-frequency interference. Compared with directly removing the single-frequency interference, this method retains the continuity of the original signal energy spectrum to a greater extent. The processed frequency domain is transformed into the time-domain signal through the inverse Fourier transform more naturally and continuously. After the microseismic signal with the single-frequency interference suppressed is processed by the sparse dictionary learning algorithm, the random noise interference in the microseismic signal can be removed. The three cooperate to adaptively achieve the hierarchical removal of complex noise. Description of the Drawings

[0053] To more clearly illustrate the technical solutions and their embodiments in the present invention, the drawings required for the technical description and embodiments are briefly introduced below. The drawings are used to provide a further understanding of the embodiments of the present invention and constitute a part of the specification, and are used together with the subsequent specific implementation manners to explain the embodiments of the present invention, but do not constitute a limitation to the embodiments of the present invention.

[0054] Figure 1 It is the processing flowchart of the signal processing method based on adaptively suppressing abnormal interference;

[0055] Figure 2 a is the time-domain diagram of the noise-free signal;

[0056] Figure 2 b is the time-frequency spectrum diagram of the noise-free signal;

[0057] Figure 2 c is the frequency-domain diagram of the noise-free signal;

[0058] Figure 3 a is the time-domain diagram of the noisy signal;

[0059] Figure 3 b is the time-frequency spectrum diagram of the noisy signal;

[0060] Figure 3 c is the frequency-domain diagram of the noisy signal;

[0061] Figure 4 a is the time-domain diagram of the signal processed by the method in this article;

[0062] Figure 4 b is the time-frequency spectrum diagram of the signal processed by the method in this article;

[0063] Figure 4 c is the frequency domain diagram of the signal processed by the method of this article. Specific embodiments

[0064] In order to more clearly express the technical advantages of the present invention, taking the synthetic microseismic signal as an example, the embodiments of the present invention will be further described in detail with reference to the accompanying drawings. The specific embodiments of the present invention are as follows:

[0065] S1. Import the microseismic signal x(t), and use the particle swarm optimization algorithm to perform parameter optimization before variational mode decomposition (VMD) to determine the global optimal mode parameter K, where t is time;

[0066] S11. Determine the sampling frequency fs, sampling time interval dt, frequency resolution df, the number of particles SearchAgents_no participating in the optimization process, the maximum number of iterations Max_iteration, specify the parameter dimension dim to be optimized, the minimum value lb of the search space, the maximum value ub of the search space, and select the fitness function fitness;

[0067] S12. Perform parameter optimization through the particle swarm optimization algorithm (PSO) to determine the global optimal mode parameter K, and use the fitness function fitness to measure the optimization effect of the current parameter K:

[0068] S2. Substitute the global optimal mode parameter K into the variational mode decomposition (VMD) to obtain K intrinsic mode functions (imf);

[0069] S3. Perform Fourier transform (FT) on each imf component to obtain the amplitude spectrum diagram of each component, then normalize its amplitude spectrum, and judge whether there is over-decomposition based on the normalized area;

[0070] S4. Calculate the ratio ratio of the maximum value to the adjacent value of the amplitude spectrum in the frequency domain of each imf;

[0071] S5. Calculate the ratio area_ratio of ratio to the corresponding Aimf i to judge whether there is single-frequency abnormal interference in the imf component. If area_ratio is greater than 0.5, it is judged that the imf component has single-frequency interference, and then a substitution operation is performed on the amplitude spectrum in the frequency domain of the imf component. If there is no single-frequency interference in the imf component, no processing is required for its frequency domain;

[0072] S6. Remove the abnormality of the imf component with single-frequency abnormal interference, and calculate the average value of the positions of 10 points at adjacent positions when the amplitude is maximum. Then the updated amplitude spectrum of the imf after removing the single-frequency abnormal interference can be expressed as AMimf n = AMimf[1,2,…,Adata [I max , …, L], when there are multiple single - frequency abnormal interferences, just repeat step S6;

[0073] S7. Superimpose the imf frequency - domain signals after single - frequency abnormal suppression, obtain the updated frequency - domain signals, and perform inverse Fourier transform;

[0074] S8. Perform sparse dictionary learning on the signal y(t) after removing single - frequency interference in step S7 to remove random noise and obtain the final signal y out (t)( Figure 4 a).

[0075] Illustration of the implementation example of the present invention:

[0076] Figure 1 It is the processing flow chart of the signal processing method based on adaptive suppression of abnormal interference in microseismic signals;

[0077] Figure 2 It is the synthesized microseismic signal and its time - frequency spectrum and amplitude spectrum diagram, Figure 2 a is the time - domain diagram of the noise - free signal. As can be seen from the figure, the time range is from 0 to 3 seconds. The amplitude value fluctuates between - 0.1 and 0.1, showing stable periodic characteristics; Figure 2 b is the time - frequency spectrum diagram of the noise - free signal. As can be seen from the figure, the time axis is from 0 to 3 seconds and the frequency range is from 0 to 150 Hz. The energy is concentrated in the frequency band below 50 Hz, and the time - frequency structure is clear without spurious frequency components; Figure 2 c is the amplitude spectrum diagram of the noise - free signal. As can be seen from the figure, the frequency range is from 0 to 150 Hz, the amplitude peaks are concentrated in the low - frequency region (0 - 50 Hz), the amplitude distribution is uniform, and the energy in the high - frequency band (> 100 Hz) approaches zero;

[0078] Figure 3 It is the noisy microseismic signal and its time - frequency spectrum and amplitude spectrum diagram, Figure 3 a is the time - domain diagram of the noisy signal. As can be seen from the figure, high - frequency random noise is superimposed on the waveform, resulting in partial masking of the original periodic characteristics; Figure 3 b is the time - frequency spectrum diagram of the noisy signal. As can be seen from the figure, strong noise interference appears at 100 HZ, the time - frequency energy distribution is diffuse, and the clarity of the low - frequency main component (0 - 50 Hz) is reduced; Figure 3 c is the frequency - domain diagram of the noisy signal. As can be seen from the figure, there is a single - frequency interference at 100 HZ in the frequency domain. Compared with the noise - free signal ( Figure 2 c), the amplitude of the low - frequency main peak decreases, and the overall signal - to - noise ratio drops significantly;

[0079] Figure 4 It is the signal processed by the technology of the present invention and its time - frequency spectrum and amplitude spectrum diagram, Figure 4a is the time-domain diagram of the signal processed by the method in this paper. It can be seen from the figure that the high-frequency noise components are effectively suppressed, and the waveform is close to that of the noise-free signal( Figure 2 a)'s periodic characteristics; Figure 4 b is the time-frequency diagram of the processed signal. By comparing it with the time-frequency diagram of the noise-added signal( Figure 3 b), it can be seen that both the single-frequency interference and the abnormal random noise in the processed time-frequency diagram are significantly suppressed, and the energy distribution of the effective signal is more focused. It has a good similarity with the time-frequency diagram of the original signal( Figure 2 b), demonstrating the advantages of the technology of this invention; Figure 4 c is the frequency-domain diagram of the processed signal. Compared with Figure 2 c, the key characteristics of the frequency domain of the noise-free signal are basically restored, showing the effectiveness of the method in this paper in restoring the frequency characteristics of the signal.

[0080] The above embodiments are only used to illustrate the present invention. The implementation steps of the method and the like can all be changed. Any equivalent transformation and improvement based on the technical solution of the present invention should not be excluded from the protection scope of the present invention.

Claims

1. An adaptive signal processing method for suppressing abnormal interference of microseismic signals, characterized in that The following specific steps are adopted: S1. Import the microseismic signal x(t), and use the particle swarm optimization algorithm to optimize the parameters before variational mode decomposition to obtain the global optimal mode parameter K, where t is time; S2. Substitute the global optimal mode parameter K into the variational mode decomposition to obtain K intrinsic mode functions. The specific formula is as follows: where, u k is a set of intrinsic mode functions, k represents the summation index variable, u k (t) represents the k-th intrinsic mode function, ω k is its central frequency, j is the imaginary unit, δ(t) is the Dirac function, * is the convolution operator, represents the partial derivative operator with respect to time t, represents the L2 norm, represents finding the combination of values of all possible {u k} and {ω k} that minimizes the objective function; S3. Perform Fourier transform on each intrinsic mode function component to obtain the amplitude spectrum diagram of each component, then normalize its amplitude spectrum, and judge whether there is over-decomposition based on the normalized area; AMimf i = FT(imf i ) If Aimf i is greater than total 0.5 times, it is considered that there is an over-decomposition situation, and steps S1 and S2 can be repeated, where E total is expressed as follows: where, imf i represents the i-th intrinsic mode function component, AMimf i represents the amplitude spectrum obtained after performing Fourier transform on the i-th intrinsic mode function component, FT[·] represents Fourier transform, Nimf i represents the frequency-domain normalization result of the i-th intrinsic mode function component, Aimf i represents the area after frequency-domain normalization of the i-th intrinsic mode function component, E total is the total energy of the entire signal x(t), f1 is the starting point of the number of points in the frequency domain, f n is the ending point of the number of points in the frequency domain, x(t) is the original time-domain signal, max[·] represents obtaining the maximum value; S4. Calculate the ratio ratio of the maximum value to the adjacent value of the amplitude spectrum in the frequency domain of each intrinsic mode function; V max = max(A data ) I max = argmax(A data ) ratio = max(L ratio , R ratio ) Where, A data represents all values in the amplitude spectrum in the frequency domain, V max is the maximum value in the amplitude spectrum corresponding to the frequency domain of the intrinsic mode function, I max is the position corresponding to the maximum value in the amplitude spectrum in the frequency domain, L ratio is the ratio of the maximum value in the amplitude spectrum in the frequency domain of the intrinsic mode function to the left neighboring value, R ratio is the ratio of the maximum value in the amplitude spectrum in the frequency domain of the intrinsic mode function to the right neighboring value, ratio is the maximum value of R ratio and L ratio ; argmax[·] represents finding the frequency point corresponding to the maximum amplitude; S5. Calculate the ratio of ratio to the corresponding Aimf i and the ratio area_ratio, and determine whether there is single-frequency abnormal interference in the intrinsic mode function components. The specific judgment formula is as follows: If area_ratio is greater than 0.5, it is judged that there is single-frequency interference in this intrinsic mode function component, and thus a value substitution operation is performed on the amplitude spectrum in the frequency domain of this intrinsic mode function component. If area_ratio is less than 0.5, it is judged that there is no single-frequency interference in this intrinsic mode function component, and no processing is required for its frequency domain; S6. Remove the anomalies from the intrinsic mode function components with single-frequency abnormal interference, and calculate the average value of the positions of 10 points at adjacent positions when the amplitude is maximum, specifically as follows: where L is the signal length, L avg is the average value of the five adjacent points to the left of the maximum value in the amplitude spectrum in the frequency domain, R avg is the average value of the five adjacent points to the right of the maximum value in the amplitude spectrum in the frequency domain, A data [j] represents the j-th value in the amplitude spectrum in the frequency domain, A data [I max is the amplitude value after removing the single-frequency interference points. Then, the updated amplitude spectrum of the intrinsic mode function after removing the single-frequency abnormal interference can be expressed as AMimf n = AMimf[1, 2, …, A data [I max , …, L]. When there are multiple single-frequency abnormal interferences, just repeat step S6; S7. Superimpose the frequency domain signals of the intrinsic mode functions after single-frequency anomaly suppression, obtain the updated frequency domain signal, and perform inverse Fourier transform; y(t) = ifft[F n (f)] where \(y(t)\) is the signal after removing the single - frequency abnormal interference, and \(F n (f)\) is the updated signal in the frequency domain, and ifft[\(\cdot\)] represents the inverse Fourier transform; S8. Perform sparse dictionary learning on the signal y(t) with single-frequency interference removed in step S7 to remove random noise, specifically: Divide the signal y(t) into NSample overlapping blocks and construct the matrix Y: Y = [y1(t), y2(t), …, y Nsample (t)] ∈ R n×NSample where y i (t) = y[1+(i-1):1+(i-1)+n-1] represents the i-th block, and the dictionary D and the sparse coefficient matrix X are solved by an optimization problem: The constraint condition is: Reconstruct the denoised block with the trained dictionary and sparse coefficients: Y d = DX ∈ R N×Sample Sum the blocks with overlap and normalize to obtain the final signal y out (t): where y(t) is the input noise signal, y i (t) is the i-th block of the signal y, n is the length of the signal block y i (t), m is the number of dictionary atoms, E is the maximum allowable error of the sparse representation, NSample is the number of overlapping blocks into which the signal is segmented, Y is the matrix after segmentation, D is the trained dictionary, X is the sparse matrix, X i is the sparse coefficient vector corresponding to the i-th signal block y i , Y d is the block matrix after denoising, Y d [·] represents the i-th matrix after denoising, Weight[·] is the normalized weight, y out (t) is the denoised signal, Start i represents the position index at which the i-th block starts to be superimposed in the final signal y i , R is the set of real numbers, where i represents an arbitrary number and j represents the position index of the signal.

Images