Microseismic Signal Denoising via Empirical Mode Decomposition, Compressed Sensing, and Soft-thresholding

Microseismic signal denoising is of great significance for P wave, S wave first arrival picking, source localization, and focal mechanism inversion. Therefore, an Empirical Mode Decomposition (EMD), Compressed Sensing (CS), and Soft-thresholding (ST) combined EMD_CS_ST denoising method is proposed. First, through EMD decomposition of the noise signal, a series of Intrinsic Mode Functions (IMF) from high frequency to low frequency are obtained. By calculating the correlation coefficient between each IMF and the original signal, the boundary component between the signal and the noise was identified, and the boundary component and its previous components were sparsely processed in the discrete wavelet transform domain to obtain the original sparse coefficient θ. Second, θ is filtered by ST to get the reconstruction coefficient θnew after denoising. Then, CS was used to recover and reconstruct θnew to get the denoised IMFnew component and then recombined with the remaining IMF components to get the signal after denoising. In the simulation experiment, the denoising process of EMD_CS_ST method is introduced in detail, and the denoising ability of EMD_CS_ST, DWT, EEMD, and VMD_DWT under 10 different noise levels is discussed. The signal-to-noise ratio, signal standard deviation, correlation coefficient, waveform diagram, and spectrogram before and after denoising are compared and analyzed. Moreover, the signals obtained from the underground cavern of the Shuangjiangkou hydropower station were denoised by the EMD_CS_ST method, and the signals before and after denoising were analyzed by time-frequency spectrum. These results show that the proposed method has better denoising ability.


Introduction
Microseismic (MS) monitoring technology, as an advanced dynamic disaster monitoring method, has been widely used in deep-buried tunnels [1], coal mines [2], underground chambers [3], rock slopes [4] and other engineering fields. Sensors are used to receive vibration signals induced by mass rock rupture. However, monitoring systems often work in high-noise environments. The signals picked up by each channel contain a large amount of external noise. These external noises are used for P and S wave pick-up, source location, source mechanism interpretation, and so on, which has caused significant interference for subsequent research work. Therefore, the signal recorded by the vibrator must be filtered to separate the MS effective signal from the noise. experiments, the application process of the algorithm is introduced in detail, the advantages and disadvantages of different noise reduction methods are compared, and the effect of noise reduction at different noise levels is discussed. The fourth part briefly introduces the source of the actual MS signal and mainly analyzes the denoising of the actual MS signal. Finally, the conclusions of this article are explained in detail.

Empirical Mode Decomposition
The EMD method assumes that any complex signal consists of some simple Intrinsic Mode Functions (IMF) and that each IMF is independent of each other [26].
The IMF component must meet the following two conditions: (1) the number of extreme values and the number of zero crossings are the same or differ by at most one; (2) the upper and lower envelopes have local symmetry about the time axis. Any signal x(t) can be decomposed and can be expressed as follows: After decomposition, n IMF components are arranged in order of frequency from high to low, and the last one is the residual component r(t). Because the MS signal has the characteristics of low frequency and high energy, the noise in the MS signal is generally distributed in the higher frequency component, and the MS effective signal dominates the lower frequency component. Therefore, after the decomposition is completed, the components with more noise in the high-frequency region are discarded, and the residual components can be reconstructed to achieve denoising. The EMD method based on local features of signals has the advantages of a multi-resolution wavelet transform. It overcomes the difficulties in choosing a wavelet basis and determining the decomposition scale in the wavelet transform. It is suitable for signal analysis, non-linear and non-stationary. However, there may be some effective MS signals in the higher frequency components. If the high-frequency components are directly removed, signal distortion will inevitably occur [27,28].
According to the advantages of EMD decomposition, some scholars denoise the noise-containing components by determining the boundary component imf k of the noise and the effective signal, and then reconstruct the residual components to obtain the denoised signal. The method of determining the boundary is as follows: Calculate the correlation coefficient Cofe between each IMF component and the original signal according to the following formula, denoted as C = {cofe 1 , cofe 2 , . . . , cofe k , . . . , cofe n }. Find the component that corresponds to the first minimum value in the coefficient C, and its next component is the boundary component between signal and noise, which is recorded as imf k .
where x(t) is the original signal, imf i (t) is the i-th IMF component, N is the number of sampling points, and mean is the mean function.

Compressed Sensing
CS theory has been widely used in information theory, image optimization, and other fields. It has accurate reconstruction ability for compressible target signals and has low requirements on sampling frequency. The basic idea of compressed sensing is the following: If the sparsity of an unknown signal is K or is transformed into K sparse by a sparse basis, then based on K sparse, combined with non-linear transformation, the original signal can be accurately reconstructed [21].
Let the expression of the noisy MS signal be: where x is a noisy MS signal, X is a valid signal, and e is random noise. The basic ideas of CS are as follows: In the first step, a set of orthogonal basis matrices Ψ(N × N) is found to thin the MS signal x and obtain the sparse vector X 0 : In the second step, select a group of measurement matrix bases Φ (M × N, K < M < N), which are not related to the orthogonal basis, to project MS signal in low dimension, from N-dimension signal to M-dimension signal. Measured value y: Combining Equations (4) and (5) gives: where θ is the coefficient to be reconstructed.
In the third step, in the case where the measured value y and the matrix Φ are known, the MS signal x is recovered from the measured value y with a high probability by an optimization method.
In the sparse transform domain, the effective signal X has sufficient sparsity, and the noise e does not have sparsity, so the noise cannot be restored in the process of signal reconstruction. This is the denoising principle of compressed sensing [22].
If the matrix Φ meets the Restricted Isometry Property (RIP) criterion, K sparse can obtain an optimal solution from the M measured values y. Since the dimension N of the coefficient to be reconstructed θ is much larger than the dimension M of the measured value y, the Equation (6) is an underdetermined system of formulas with an infinite set of solutions. Taking advantage of the sparsity of theta and using L1 norm regularization formula to find the optimal solution, it is a typical optimization algorithm to solve the underdetermined formula [29]: Equation (7) is the theoretical basis of the compressed sensing algorithm, where || || p denotes the L p norm; the first term in the braces is the constraint measurement term; the second term is the sparse constraint term, and the parameter α is the balance parameter. This paper chooses a sparse Bayesian learning algorithm that is more suitable for solving Equation (7) than the traditional L 1 norm regularization algorithm [30].
At present, commonly used measurement matrices include Gaussian random measurement matrices, random Bernoulli matrices, partially orthogonal matrices, Topliz and cyclic matrices and sparse random matrices, and other matrices. The Gaussian random measurement matrix is used as a universal compressed sensing matrix [31]. For this reason, it is used as the measurement matrix in this paper. The core idea is to construct the M × N measurement matrix Φ and make each element in the measurement matrix Φ obey the N(0,1/M) Gaussian distribution and conform to the RIP criterion. There is no correlation with most orthogonal bases [32]. The classical sparse transform methods include discrete cosine transform [33], discrete Fourier transform [34], discrete wavelet transform [35], and other transforms. In this paper, to organically combine with the following soft-thresholding denoising methods, a discrete wavelet transform is selected.

Soft-Thresholding
Wavelet threshold denoising is to distribute the energy of MS signals mainly in larger wavelet coefficients, and noise primarily distributed in smaller wavelet coefficients. Based on this distribution characteristic, the denoising can be achieved by converting the time domain signal into the wavelet domain and filtering the noise by setting an appropriate threshold value [36].
The key of the wavelet threshold is to select the appropriate threshold T to perform threshold quantization on the wavelet coefficients w i,j of the noisy signal, to be as close as possible to the wavelet coefficient w of the original signal. Processed wavelet coefficients are represented as w 0 , and wavelet reconstruction is performed to obtain a denoising signal.
Currently widely used is the data-driven threshold proposed by Donoho [36], the expression is as follows: where Y i is the i-th IMF component, median is the median function, and T i is the threshold of the i-th reconstruction coefficient θ i . Threshold denoising has a hard threshold and a soft thresholding. The hard thresholding function is better than the soft thresholding function in the sense of mean square error, but the signal will produce additional shocks, generate jumping points, and do not have the smoothness of the original signal, which greatly affects the initial pick-up of P wave and S wave of MS signal and the location of MS source. The overall continuity of the wavelet coefficients obtained from the soft threshold estimation is excellent, so that the estimated signal will not generate additional shocks, which is conducive to the later research work. Therefore, soft-thresholding denoising is used in this paper. The expression is: where sgn is a symbol-function, and θ i,new is a reconstruction coefficient after filtering the i-th IMF component.

The Shortcomings of CS Denoising
To illustrate the problem of poor denoising performance when CS is facing a signal with a low signal-to-noise ratio. This paper verifies by digital simulation experiments. The Picker wavelet of the simulated seismic signal is used in the simulation experiment. The Picker wavelet can be expressed as: where f (t) is the amplitude, and f p is the peak frequency. In the simulation analysis, f p = 30 Hz, the sampling number is 1000, the sampling frequency is 1000 Hz, and random noise is added. The noise belongs to Gaussian distribution, and the variance is 0.09. Their waveforms and spectrums are shown in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 16 wavelet coefficient w of the original signal. Processed wavelet coefficients are represented as w0, and wavelet reconstruction is performed to obtain a denoising signal. Currently widely used is the data-driven threshold proposed by Donoho [36], the expression is as follows: where Yi is the i-th IMF component, median is the median function, and Ti is the threshold of the i-th reconstruction coefficient θi. Threshold denoising has a hard threshold and a soft thresholding. The hard thresholding function is better than the soft thresholding function in the sense of mean square error, but the signal will produce additional shocks, generate jumping points, and do not have the smoothness of the original signal, which greatly affects the initial pick-up of P wave and S wave of MS signal and the location of MS source. The overall continuity of the wavelet coefficients obtained from the soft threshold estimation is excellent, so that the estimated signal will not generate additional shocks, which is conducive to the later research work. Therefore, soft-thresholding denoising is used in this paper. The expression is: where sgn is a symbol-function, and θi,new is a reconstruction coefficient after filtering the i-th IMF component.

The Shortcomings of CS Denoising
To illustrate the problem of poor denoising performance when CS is facing a signal with a low signal-to-noise ratio. This paper verifies by digital simulation experiments. The Picker wavelet of the simulated seismic signal is used in the simulation experiment. The Picker wavelet can be expressed as: where f(t) is the amplitude, and fp is the peak frequency. In the simulation analysis, fp = 30 Hz, the sampling number is 1000, the sampling frequency is 1000 Hz, and random noise is added. The noise belongs to Gaussian distribution, and the variance is 0.09. Their waveforms and spectrums are shown in Figure 1.
Compressed sensing reconstruction of the noisy Ricker wavelet is shown in Figure 2. It can be known from Figure 2b that the original signal does not have sufficient sparsity in the sparse domain. After passing the CS algorithm, it can be seen from Figure 2c,d that the noise is removed, but the effect is not obvious, and the NF phenomenon (red circle) appears.
If soft-thresholding filtering is performed on the reconstruction coefficient θ and the sparse coefficient corresponding to the noise is removed, the obtained result is as shown in Figure 3. Although significant denoising was achieved, some NF phenomena (red circles) still appeared.  Compressed sensing reconstruction of the noisy Ricker wavelet is shown in Figure 2. It can be known from Figure 2b that the original signal does not have sufficient sparsity in the sparse domain. After passing the CS algorithm, it can be seen from Figure 2c,d that the noise is removed, but the effect is not obvious, and the NF phenomenon (red circle) appears.

Improvement of CS Denoising
In view of the problems in Section 2.4.1, the MS signal is decomposed by EMD, and combined with the advantage of CS theory with certain denoising and soft-thresholding denoising, a method of denoising based on EMD_CS_ST is proposed. The specific steps of the denoising method are shown below.
Step 1: The EMD method is used to decompose the MS signal x(t) into a series of high-to-low frequency components, containing n IMF components and a residual component, denoted as Step 2: Calculated the correlation coefficient of n IMF components by Equation (2), denoted as C = {cofe1, cofe2, ..., cofek, ..., cofen}, and found the first local minimum value of C, and the next IMF component corresponding to it is the boundary between noise and signal, recorded as imfk; Step 3: Performed discrete wavelet transform sparse processing on the imfk component and the components before it, and calculated the coefficient θ to be reconstructed of k components, which are denoted as θ = {θ1, θ2, ..., θk}; Step 4: Performed soft-thresholding filtering on the θ to obtain k reconstruction coefficients θnew, which are denoted by θnew = {θ1,new, θ2,new, ..., θk,new }; Step 5: Reconstructed the reconstruction coefficient θnew by CS to obtain k IMFnew, denoted as IMFnew = {imf1,new, imf2,new, ..., imfk,new }; Step 6: Recombined the k IMFnew components with the Remaining Components (RCs) to obtain the denoised signal, which is the effective MS signal x(t)new after denoising.
The EMD_CS_ST algorithm flow chart is shown in Figure 4. If soft-thresholding filtering is performed on the reconstruction coefficient θ and the sparse coefficient corresponding to the noise is removed, the obtained result is as shown in Figure 3. Although significant denoising was achieved, some NF phenomena (red circles) still appeared.

Improvement of CS Denoising
In view of the problems in Section 2.4.1, the MS signal is decomposed by EMD, and combined with the advantage of CS theory with certain denoising and soft-thresholding denoising, a method of denoising based on EMD_CS_ST is proposed. The specific steps of the denoising method are shown below.
Step 1: The EMD method is used to decompose the MS signal x(t) into a series of high-to-low frequency components, containing n IMF components and a residual component, denoted as I = {imf1, imf2, ..., imfk, ..., imfn, r(t)}; Step 2: Calculated the correlation coefficient of n IMF components by Equation (2), denoted as C = {cofe1, cofe2, ..., cofek, ..., cofen}, and found the first local minimum value of C, and the next IMF component corresponding to it is the boundary between noise and signal, recorded as imfk; Step 3: Performed discrete wavelet transform sparse processing on the imfk component and the components before it, and calculated the coefficient θ to be reconstructed of k components, which are denoted as θ = {θ1, θ2, ..., θk}; Step 4: Performed soft-thresholding filtering on the θ to obtain k reconstruction coefficients θnew, which are denoted by θnew = {θ1,new, θ2,new, ..., θk,new }; Step 5: Reconstructed the reconstruction coefficient θnew by CS to obtain k IMFnew, denoted as IMFnew = {imf1,new, imf2,new, ..., imfk,new }; Step 6: Recombined the k IMFnew components with the Remaining Components (RCs) to obtain the denoised signal, which is the effective MS signal x(t)new after denoising.
The EMD_CS_ST algorithm flow chart is shown in Figure 4.

Improvement of CS Denoising
In view of the problems in Section 2.4.1, the MS signal is decomposed by EMD, and combined with the advantage of CS theory with certain denoising and soft-thresholding denoising, a method of denoising based on EMD_CS_ST is proposed. The specific steps of the denoising method are shown below.
Step 1: The EMD method is used to decompose the MS signal x(t) into a series of high-to-low frequency components, containing n IMF components and a residual component, denoted as I = {imf 1 , imf 2 , . . . , imf k , . . . , imf n , r(t)}; Step 2: Calculated the correlation coefficient of n IMF components by Equation (2), denoted as C = {cofe 1 , cofe 2 , . . . , cofe k , . . . , cofe n }, and found the first local minimum value of C, and the next IMF component corresponding to it is the boundary between noise and signal, recorded as imf k ; Step 3: Performed discrete wavelet transform sparse processing on the imf k component and the components before it, and calculated the coefficient θ to be reconstructed of k components, which are denoted as θ = {θ 1 , θ 2 , . . . , θ k }; Step 4: Performed soft-thresholding filtering on the θ to obtain k reconstruction coefficients θ new , which are denoted by θ new = {θ 1,new , θ 2,new , . . . , θ k,new }; Step 5: Reconstructed the reconstruction coefficient θ new by CS to obtain k IMF new , denoted as IMF new = {imf 1,new , imf 2,new , . . . , imf k,new }; Step 6: Recombined the k IMF new components with the Remaining Components (RCs) to obtain the denoised signal, which is the effective MS signal x(t) new after denoising.
The EMD_CS_ST algorithm flow chart is shown in Figure 4.

EMD_CS_ST Algorithm Application
Using the digital simulation analysis in Section 2.4.1, the EMD adaptive decomposition of the noisy Ricker wavelet is performed to obtain 7 IMF components. The waveforms of each component and its spectrum are shown in Figure 5.

EMD_CS_ST Algorithm Application
Using the digital simulation analysis in Section 2.4.1, the EMD adaptive decomposition of the noisy Ricker wavelet is performed to obtain 7 IMF components. The waveforms of each component and its spectrum are shown in Figure 5.
The calculated cross-correlation coefficient is shown in Figure 6 with bar and dot charts. According to the demarcation point discrimination method of noise and signal, it can be known that imf 6 is a boundary. Then, the first six IMF components are respectively subjected to sparse transformation to obtain coefficients to be reconstructed θ = {θ 1 , θ 2 , . . . , θ 6 }.
The threshold T is calculated according to the algorithm of Section 2.3, and the calculation result is shown in Figure 7, and soft-thresholding filtering is performed on the reconstruction coefficient θ. Finally, the denoised IMF new component is reconstructed by discrete wavelet inverse transform, and the remaining IMF component is combined with IMF new to obtain the denoised MS signal. The calculated cross-correlation coefficient is shown in Figure 6 with bar and dot charts. According to the demarcation point discrimination method of noise and signal, it can be known that imf6 is a boundary. Then, the first six IMF components are respectively subjected to sparse transformation to obtain coefficients to be reconstructed θ = {θ1, θ2, ..., θ6}. The threshold T is calculated according to the algorithm of Section 2.3, and the calculation result is shown in Figure 7, and soft-thresholding filtering is performed on the reconstruction The calculated cross-correlation coefficient is shown in Figure 6 with bar and dot charts. According to the demarcation point discrimination method of noise and signal, it can be known that imf6 is a boundary. Then, the first six IMF components are respectively subjected to sparse transformation to obtain coefficients to be reconstructed θ = {θ1, θ2, ..., θ6}. The threshold T is calculated according to the algorithm of Section 2.3, and the calculation result is shown in Figure 7, and soft-thresholding filtering is performed on the reconstruction  To further verify the effectiveness and feasibility of the proposed method, this paper compares with traditional discrete wavelet transform (DWT) denoising, EEMD denoising, and VMD_DWT [37] denoising methods. The simulation experiments of these denoising methods are carried out on the MATLAB platform. For the DWT denoising method of MS signals, the perform function is wdenoise, sym is used as the wavelet base, the number of filters is 8, and the number of decomposition layers is 5. The threshold rule is soft-threshold; estimate noise using only finest-scale wavelet coefficients. According to the signal-to-noise ratio (SNR), the signal standard deviation (SD), the correlation coefficient (CC) of the pure signal, and the denoising signal, these three indicators quantitatively illustrate the denoising effect of this paper. The waveforms and their spectra after denoising by these four methods are shown in Figure 8.
The SNR is defined as: x t x t (11) where x(t) is a pure signal, x(t)new is a denoised signal, and N is the number of sampling points. If the signal has a high SNR after denoising, the denoising effect is better. The SD is defined as: If the SD is smaller, it means that the smaller the fluctuation between the denoising signal and the pure signal, the closer to the pure signal.
The CC between the pure signal and the denoising signal is defined as: x t x t (13) where Cov is the covariance function, Var is the variance function. If the CC value is closer to 1, the denoising signal is similar to the clean signal. Apply the Equations (11)- (13) to calculate the SNR before and after denoising of the noisy Picker wavelet, and the SD and CC after denoising. These results are shown in Table 1. To further verify the effectiveness and feasibility of the proposed method, this paper compares with traditional discrete wavelet transform (DWT) denoising, EEMD denoising, and VMD_DWT [37] denoising methods. The simulation experiments of these denoising methods are carried out on the MATLAB platform. For the DWT denoising method of MS signals, the perform function is wdenoise, sym is used as the wavelet base, the number of filters is 8, and the number of decomposition layers is 5. The threshold rule is soft-threshold; estimate noise using only finest-scale wavelet coefficients. According to the signal-to-noise ratio (SNR), the signal standard deviation (SD), the correlation coefficient (CC) of the pure signal, and the denoising signal, these three indicators quantitatively illustrate the denoising effect of this paper. The waveforms and their spectra after denoising by these four methods are shown in Figure 8.   It can be seen from Table 1 and Figure 8 that the EMD_CS_ST has a better denoising effect than the other three methods. Analysis from Table 1 shows that the SNR of the noisy Ricker wavelet is 4.651 dB. After applying this method, the SNR is increased to 27.170 dB, the signal standard deviation is the smallest, and the correlation coefficient is the largest. Transient non-stationary characteristics of the MS signal are fully preserved compared to the other three methods.

Comparison of Denoising Effects at Different Noise Levels
To test the performance of the four methods at different noise levels, we added 10 sets of random noise at different levels to the pure signal, and the noise variance changed from 0.02 to 0.20. Signals with different input noise levels, denoised signals, and SNR are shown in Figure 9. Figure 9f shows the SNR at different noise levels and the SNR after denoising. The red line is the proposed method. It can be seen directly from Figure 9f that the red line corresponding to this method is at the top of all other lines, indicating that the denoising effect is the best. When faced with high SNR signals, the proposed method is equivalent to DWT denoising, but when faced with low SNR The SNR is defined as: where x(t) is a pure signal, x(t) new is a denoised signal, and N is the number of sampling points. If the signal has a high SNR after denoising, the denoising effect is better.
The SD is defined as: If the SD is smaller, it means that the smaller the fluctuation between the denoising signal and the pure signal, the closer to the pure signal.
The CC between the pure signal and the denoising signal is defined as: where Cov is the covariance function, Var is the variance function. If the CC value is closer to 1, the denoising signal is similar to the clean signal. Apply the Equations (11)- (13) to calculate the SNR before and after denoising of the noisy Picker wavelet, and the SD and CC after denoising. These results are shown in Table 1. It can be seen from Table 1 and Figure 8 that the EMD_CS_ST has a better denoising effect than the other three methods. Analysis from Table 1 shows that the SNR of the noisy Ricker wavelet is 4.651 dB. After applying this method, the SNR is increased to 27.170 dB, the signal standard deviation is the smallest, and the correlation coefficient is the largest. Transient non-stationary characteristics of the MS signal are fully preserved compared to the other three methods.

Comparison of Denoising Effects at Different Noise Levels
To test the performance of the four methods at different noise levels, we added 10 sets of random noise at different levels to the pure signal, and the noise variance changed from 0.02 to 0.20. Signals with different input noise levels, denoised signals, and SNR are shown in Figure 9. Figure 9f shows the SNR at different noise levels and the SNR after denoising. The red line is the proposed method. It can be seen directly from Figure 9f that the red line corresponding to this method is at the top of all other lines, indicating that the denoising effect is the best. When faced with high SNR signals, the proposed method is equivalent to DWT denoising, but when faced with low SNR signals, the proposed method seems to perform better than other methods. This phenomenon indicates that the proposed method is more effective for signals with higher noise levels. Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 16

Description of the MS Monitoring System and Actual MS Signal
The actual MS signal comes from the safety construction monitoring system of the underground chamber of Shuangjiangkou Hydropower Station in Southwest China. In view of the severe problem of rock mass deformation and rockburst disasters in the underground main powerhouse, the seismic monitoring system of the Canadian Engineering Seismology Group (ESG) was adopted [38][39][40]. The seismic monitoring system mainly includes Paladin digital signal acquisition system, Hyperion digital signal processing system, acceleration sensor, cable fiber, and so on, as shown in Figure 10a. In the underground main powerhouse, two groups of acceleration sensors were installed; each group has six. They were installed in the upper and middle drainage gallery, as shown in Figure 10b. The monitoring region mainly covers the main powerhouse, involving the main transformer chamber, the surrounding rock of part of the pressure pipeline. was adopted [38][39][40]. The seismic monitoring system mainly includes Paladin digital signal acquisition system, Hyperion digital signal processing system, acceleration sensor, cable fiber, and so on, as shown in Figure 10a. In the underground main powerhouse, two groups of acceleration sensors were installed; each group has six. They were installed in the upper and middle drainage gallery, as shown in Figure 10b. The monitoring region mainly covers the main powerhouse, involving the main transformer chamber, the surrounding rock of part of the pressure pipeline.  The acceleration sensors with a frequency response of 20 Hz to 5000 Hz were arranged in a spatial network distribution. The sensitivity of the acceleration sensor is 30 V/g, and the Paladin digital signal acquisition system has a sampling frequency of 10 kHz. According to the Short Time Average and Longtime Average algorithms (STA/LTA), the threshold of the trigger channel is set to 4. If multiple sensors detect the seismic waves generated by a rock fracture and the corresponding value exceeds the threshold, the sensors record them immediately. These recorded electrical signals will be sent to the digital acquisition system for processing and transmitted to WaveVis software via fiber optic cables to obtain a valid MS signal [41]. Often, the signal received by the system is not only the signal generated by rock fracture but also some blasting, mechanical vibration, current sound, and other noise signals. These mixed signals need to be denoised. After the acquisition, these electrical signals are transmitted to the processing system via fiber optic cables.

Field Data Examples
The ESG system has been operating in the main powerhouse of the Shuangjiangkou hydropower station for two years, during which thousands of MS events were obtained. Three groups of effective MS signals were randomly selected for EMD_CS_ST de-noising processing and analysis. To facilitate calculation and comparison, the effective MS signals are intercepted 1000 ms from each time-window of 1500 ms. Firstly, three groups of signals are denoised by the method of EMD_CS_ST, and then using the MATLAB platform to perform FFT spectrum analysis, it can more intuitively judge the denoising effect of MS signals. The waveform and time spectrum before and after denoising are shown in Figure 11.
will be sent to the digital acquisition system for processing and transmitted to WaveVis software via fiber optic cables to obtain a valid MS signal [41]. Often, the signal received by the system is not only the signal generated by rock fracture but also some blasting, mechanical vibration, current sound, and other noise signals. These mixed signals need to be denoised. After the acquisition, these electrical signals are transmitted to the processing system via fiber optic cables.

Field Data Examples
The ESG system has been operating in the main powerhouse of the Shuangjiangkou hydropower station for two years, during which thousands of MS events were obtained. Three groups of effective MS signals were randomly selected for EMD_CS_ST de-noising processing and analysis. To facilitate calculation and comparison, the effective MS signals are intercepted 1000 ms from each time-window of 1500 ms. Firstly, three groups of signals are denoised by the method of EMD_CS_ST, and then using the MATLAB platform to perform FFT spectrum analysis, it can more intuitively judge the denoising effect of MS signals. The waveform and time spectrum before and after denoising are shown in Figure 11. It can be seen from Figure 11 that the three groups of original MS signals all contain a large amount of random noise. Even if the main spectrum distribution of these signals can be seen from the time spectrum, it is not obvious due to noise interference. After the denoising processing of this method, the waveform and time spectrum of the MS signal are clearer, the dominant spectrum is more visible, and there is almost no NF phenomenon. The EMD_CS_ST method was adopted to It can be seen from Figure 11 that the three groups of original MS signals all contain a large amount of random noise. Even if the main spectrum distribution of these signals can be seen from the time spectrum, it is not obvious due to noise interference. After the denoising processing of this method, the waveform and time spectrum of the MS signal are clearer, the dominant spectrum is more visible, and there is almost no NF phenomenon. The EMD_CS_ST method was adopted to denoise MS signal, and the signal denoised has the advantages of high SNR and stable waveform. The effectiveness and practicability of the proposed method are proved qualitatively.

Conclusions
Considering that the MS signal contains a large amount of random noise, an EMD_CS_ST denoising method was proposed, and the advantages and disadvantages of the EMD_CS_ST, DWT, EEMD, and VMD_DWT are compared and analyzed through a numerical simulation experiment. The denoising analysis of MS signals obtained from the underground cavern of Shuangjiangkou hydropower station was carried out. The following conclusions are drawn.
The numerical simulation experiment demonstrates the NF phenomenon of the compressed sensing denoising method in the process of signal reconstruction with low SNR. Although the theory of CS has a certain ability to reduce noise in the denoising of MS signals, the effect is weak.
The CS denoising method is not suitable for direct application to the denoising processing of MS signals. The EMD_CS_ST method can better compensate for the deficiency of CS in MS signal denoising and greatly weaken the NF phenomenon caused by CS.
Simulation experiments were conducted to compare the effects of four kinds of MS signal denoising methods such as EMD_CS_ST, DWT, EEMD, and VMD_DWT. The FFT time spectrum of MS signals before and after denoising were compared and analyzed. The feasibility of the EMD_CS_ST method is demonstrated qualitatively.
By introducing the SNR, the SD, and the CC, the rationality and effectiveness of the EMD_CS_ST method were proved quantitatively. The actual three MS signals were denoised, and the effect was remarkable.
The CS technology itself has matured, and we can apply EMD_CS_ST to data collection, storage, and noise reduction in the later stage to achieve integrated data processing of the three, which can greatly reduce the space occupied by data.