# A Nonparametric Method for Automatic Denoising of Microseismic Data

Authors: Pingan Peng; Liguan Wang
Journal: Shock and Vibration
 (2018)

---

## Abstract

Noise suppression or signal-to-noise ratio (SNR) enhancement is often desired for better processing results from a microseismic dataset. In this paper, we proposed a nonparametric automatic denoising algorithm for microseismic data. The method consists of three major steps: (1) applying a two-step AIC algorithm to pick P-wave arrival; (2) subtracting the noise power spectrum from the signal power spectrum; (3) recovering the microseismic signal by inverse Fourier transform. The proposed method is tested on synthetic datasets with different signal types and SNRs, as well as field datasets. The results of the proposed method are compared against ensemble empirical mode decomposition (EEMD) and wavelet denoising methods, which shows the effectiveness of the method for denoising and improving the SNR of microseismic data.


---

## Body

## 1. Introduction

In the past few decades, microseismic monitoring has been widely used in a variety of applications, ranging from hydraulic fracturing to mining and geotechnical engineering [1, 2]. The signal is crucial for microseismic monitoring, and its quality will significantly affect the results of event detection, location, and source mechanism estimation. However, there are plenty of noise sources due to the complex situation of the project site, such as artificial activities, mechanical equipment, electrical interference, and equipment noise, which can severely affect the quality of the seismograms [3]. As a result, microseismic signals are often strongly contaminated by unwanted noise.Noise suppression is a key step in microseismic data processing, aiming at improving the signal-to-noise ratio (SNR) [4]. Hence, numerous methods have been presented for this in the past few decades that can be classified into four main categories: (1) infinite impulse response (IIR) filters, known as low-pass, high-pass, band-pass, and band-stop filter (e.g., [5]); (2) Fourier-based filters (e.g., [6–8]); (3) EMD-based (empirical mode decomposition) filters (e.g., [9–14]); (4) wavelet-based filters (e.g., [15–17]).Huang et al. [18] introduced the EMD algorithm to the signal processing community and it has been widely used since then. Bekara and Van der Baan [9] proposed a frequency-offset EMD technique to suppress the seismic random and coherent noise. In order to eliminate mode mixing present in the original EMD algorithm, ensemble empirical mode decomposition (EEMD) method has been recently proposed by Wu and Huang [19]. Han and van der Baan [10] developed a novel seismic and microseismic denoising method based on EEMD combined with adaptive thresholding. To et al. [20] compared Fourier-based and wavelet-based denoising techniques applied to geophysical data. Gaci [21] studied soft thresholding denoising techniques based on the discrete wavelet transform to enhance the first-arrival picking. In addition to Mousavi et al. [1], Mousavi and Langston [22] also used synchrosqueezed continuous wavelet transform (SS-CWT) for seismic denoising.The aforementioned methods have played an important role in enhancing the quality of microseismic signals, but there is still room for SNR improvement. On the other hand, many of these methods have a large number of parameters, which must be tuned for each data set, and the goal of our work is to reduce the number of adjustable parameters that make the process more automatic.In this paper, we introduce a nonparametric method for automatic denoising of microseismic data based on Akaike's Information Criterion (AIC) and Fourier transform, called PD (Pick & Denoise) method. In the following sections, a brief theoretical introduction to Akaike's Information Criterion and Fourier transform will be presented, as well as the details of our method. We analyze the performance of the proposed method applied to both synthetic and field microseismic data in terms of efficiency and signal preservation and contrast the results against standard EEMD and wavelet denoising methods.


## 2. Theoretical Background and Development

### 2.1. Akaike’s Information Criterion

AIC is one of the commonly used methods for P-wave arrival picking in both the seismic and microseismic field. It is assumed that a seismogram can be divided into locally stationary segments, each modeled as an autoregressive process, and that the intervals before and after the onset time are two different stationary processes [23]. The AIC function for seismogram of lengthN is defined as(1)AICk=k-Mlogσ1,max2+N-M-klogσ2,max2+C where k is thekth sampling point, k=1,2,3,…,N; M is the order of the autoregressive process; C is a constant; σ1,max2 and σ2,max2 indicate the variance of the seismogram in the two intervals not explained by the autoregressive, respectively. Since the calculation of (1) is very complex, Maeda [24] proposed a simpler method as follows:(2)AICk=klog⁡varx1,k+N-k-1log⁡varxk+1,Nwherevar represents the variance of the data, var=1/N-1∑k=1N|xk-1/N∑k=1Nxk|2. The P-wave arrival is the corresponding point where the AIC has a minimum value.AIC function will have a minimum value for any input. When SNR is relatively low, the global minimum cannot be guaranteed to indicate the P-wave arrival. Thus, the accuracy of Maeda’s AIC method depends heavily on the choice of the time window.


### 2.2. Fourier Transform and Its Inverse Transform

Fourier transform is an important tool to convert the signal from time domain to frequency domain, and its inverse transform can be used to restore the signal from frequency domain to time domain. For signalx(j), j=1,2,3,…,N, the Fourier transform and its inverse transform can be expressed as(3)Yk=∑j=1Nxje-i2πj-1k-1/N(4)xj=1N∑k=1NYkei2πj-1k-1/NBecause Discrete Fourier Transform (DFT) requires large computational cost, we generally reduce the amount of computation by decomposing the time series into odd and even subsequences by Fast Fourier transform (FFT).For a signalx(j) of length N=2M (M is a positive integer; zero-padding is used if the signal length is insufficient), we split theN-point data sequence into twoN/2-point data sequences x1(j) and x2(j), corresponding to the even-numbered and odd-numbered samples of x(j), respectively; that is,(5)x1j=x2j,j=1,2,…,N2x2j=x2j-1,j=1,2,…,N2Thusx1(j) and x2(j) are obtained by decimating x(j) by a factor of 2, and the above decimation-in-time decomposition can be performed M-1 times, until each DFT is length 2. A length 2 DFT requires no multiplies, and the computational cost is evidently reduced. Finally, the N-point DFT can be expressed as follows:(6)Yk=Y1k+WNkY2k,k=1,2,…,N2(7)Yk+N2=Y1k-WNkY2k,k=1,2,…,N2where WNk is the rotation factor, Y1(k) is the DFT of the odd sequence, and Y2(k) is the DFT of the even sequence.


### 2.3. The Proposed Method

A typical microseismic signal consists of three segments: the presignal noise, microseismic event, and background noise level, as shown in Figure1. Normally, the microseismic event starts from the P-wave arrival and ends at the S coda, containing the key information of rock fracture. If we can extract the noise-related information from the presignal noise and then subtract this part from the entire signal, we can restore the original signal better. According to this idea, a nonparametric automatic denoising method is proposed for microseismic data as well as earthquake data that includes presignal noise.Figure 1
Classical microseismic event waveform.The first step is to apply a modified AIC autopicker [24] to find the P-wave arrival. Then, the spectral component of the noise is subtracted from the spectrum of the noisy microseismic signal in the frequency domain. Finally, the denoised signal is restored by inverse Fourier transform. Our proposed method is automatic and does not need to tune any parameter. The implementation procedure for PD method is as follows:(1)
Lety=x(i) be the noisy microseismic signal, calculate the value of the characteristic function according to (8), and denote the sampling point corresponding to the maximum value of the characteristic function as R.(8)CFi=xi+4xi-xi-1(2)
For time interval[1,R], apply the AIC method and denote the primitive P-wave arrival as Q.(3)
For time interval[Q-(R-Q)/4,Q+(R-Q)/4], apply the AIC method again and denote the corrected P-wave arrival as T. Thus we get the presignal [1,T], denoted by n(i).(4)
Add Hamming window and apply Fast Fourier transform tox(i) and n(i), respectively, then we get Xw(i) and Nw(i).(5)
Calculate the average amplitudeXw(i) and phase angle φ(ω) of Xw(i); calculate the average amplitude Nw(i) of Nw(i). Since the microseismic event is assumed to be uncorrelated with the presignal and background noise, the clean signal can be estimated by subtracting a noise estimate from the noisy signal. Hence,(9)Dwi2=Xwi2-Nwi2(6)
Use the phase angleφ(ω) instead of the phase angle of the noise-free signal, then the denoising microseismic signal can be restored by inverse Fourier transform.(10)yi=ifftDwieiφω


## 2.1. Akaike’s Information Criterion

AIC is one of the commonly used methods for P-wave arrival picking in both the seismic and microseismic field. It is assumed that a seismogram can be divided into locally stationary segments, each modeled as an autoregressive process, and that the intervals before and after the onset time are two different stationary processes [23]. The AIC function for seismogram of lengthN is defined as(1)AICk=k-Mlogσ1,max2+N-M-klogσ2,max2+C where k is thekth sampling point, k=1,2,3,…,N; M is the order of the autoregressive process; C is a constant; σ1,max2 and σ2,max2 indicate the variance of the seismogram in the two intervals not explained by the autoregressive, respectively. Since the calculation of (1) is very complex, Maeda [24] proposed a simpler method as follows:(2)AICk=klog⁡varx1,k+N-k-1log⁡varxk+1,Nwherevar represents the variance of the data, var=1/N-1∑k=1N|xk-1/N∑k=1Nxk|2. The P-wave arrival is the corresponding point where the AIC has a minimum value.AIC function will have a minimum value for any input. When SNR is relatively low, the global minimum cannot be guaranteed to indicate the P-wave arrival. Thus, the accuracy of Maeda’s AIC method depends heavily on the choice of the time window.


## 2.2. Fourier Transform and Its Inverse Transform

Fourier transform is an important tool to convert the signal from time domain to frequency domain, and its inverse transform can be used to restore the signal from frequency domain to time domain. For signalx(j), j=1,2,3,…,N, the Fourier transform and its inverse transform can be expressed as(3)Yk=∑j=1Nxje-i2πj-1k-1/N(4)xj=1N∑k=1NYkei2πj-1k-1/NBecause Discrete Fourier Transform (DFT) requires large computational cost, we generally reduce the amount of computation by decomposing the time series into odd and even subsequences by Fast Fourier transform (FFT).For a signalx(j) of length N=2M (M is a positive integer; zero-padding is used if the signal length is insufficient), we split theN-point data sequence into twoN/2-point data sequences x1(j) and x2(j), corresponding to the even-numbered and odd-numbered samples of x(j), respectively; that is,(5)x1j=x2j,j=1,2,…,N2x2j=x2j-1,j=1,2,…,N2Thusx1(j) and x2(j) are obtained by decimating x(j) by a factor of 2, and the above decimation-in-time decomposition can be performed M-1 times, until each DFT is length 2. A length 2 DFT requires no multiplies, and the computational cost is evidently reduced. Finally, the N-point DFT can be expressed as follows:(6)Yk=Y1k+WNkY2k,k=1,2,…,N2(7)Yk+N2=Y1k-WNkY2k,k=1,2,…,N2where WNk is the rotation factor, Y1(k) is the DFT of the odd sequence, and Y2(k) is the DFT of the even sequence.


## 2.3. The Proposed Method

A typical microseismic signal consists of three segments: the presignal noise, microseismic event, and background noise level, as shown in Figure1. Normally, the microseismic event starts from the P-wave arrival and ends at the S coda, containing the key information of rock fracture. If we can extract the noise-related information from the presignal noise and then subtract this part from the entire signal, we can restore the original signal better. According to this idea, a nonparametric automatic denoising method is proposed for microseismic data as well as earthquake data that includes presignal noise.Figure 1
Classical microseismic event waveform.The first step is to apply a modified AIC autopicker [24] to find the P-wave arrival. Then, the spectral component of the noise is subtracted from the spectrum of the noisy microseismic signal in the frequency domain. Finally, the denoised signal is restored by inverse Fourier transform. Our proposed method is automatic and does not need to tune any parameter. The implementation procedure for PD method is as follows:(1)
Lety=x(i) be the noisy microseismic signal, calculate the value of the characteristic function according to (8), and denote the sampling point corresponding to the maximum value of the characteristic function as R.(8)CFi=xi+4xi-xi-1(2)
For time interval[1,R], apply the AIC method and denote the primitive P-wave arrival as Q.(3)
For time interval[Q-(R-Q)/4,Q+(R-Q)/4], apply the AIC method again and denote the corrected P-wave arrival as T. Thus we get the presignal [1,T], denoted by n(i).(4)
Add Hamming window and apply Fast Fourier transform tox(i) and n(i), respectively, then we get Xw(i) and Nw(i).(5)
Calculate the average amplitudeXw(i) and phase angle φ(ω) of Xw(i); calculate the average amplitude Nw(i) of Nw(i). Since the microseismic event is assumed to be uncorrelated with the presignal and background noise, the clean signal can be estimated by subtracting a noise estimate from the noisy signal. Hence,(9)Dwi2=Xwi2-Nwi2(6)
Use the phase angleφ(ω) instead of the phase angle of the noise-free signal, then the denoising microseismic signal can be restored by inverse Fourier transform.(10)yi=ifftDwieiφω


## 3. Tests on Synthetic Data

### 3.1. Results for Signals with Different Types

We first apply the algorithm to synthetic microseismic signals to test its performance. The procedure for creating synthetic signals is based on the work of Halldorsson and Papageorgiou [25]. A noisy sample was generated by adding white Gaussian noise to a synthetic signal. We have tested the algorithm for synthetic examples with different noise level and compared it with the EEMD denoising [26] and wavelet denoising method [27], which are widely used in noise reduction of microseismic signals.We use the mean absolute errorμ and the standard deviation σ as the evaluation criteria of the denoising performance. μ represents the extent of discrimination between the denoised signal and the real signal. The smaller the value is, the more similar the denoised signal is to the real signal. σ is used to quantify the amount of variation or dispersion. A low standard deviation indicates that the data points tend to be close to the mean value, while a high standard deviation indicates that the data points are spread out over a wider range of values, defined as follows:(11)μ=1N∑i=1Ndi-si(12)σ=1N-1∑i=1Ndi-si-λ2where λ=1/N∑i=1Ndi-si,d is the denoising signal, s is the real signal, and N is the number of signal points.Four types of signals that include presignal noise are synthesized: sinusoidal signal, double sinusoidal signal, blast vibration signal, and microseismic signal. The amplitude of the signal with additive white Gaussian noise is normalized to±1. After that, we denoise the signal by PD, EEMD, and wavelet denoising method. Results are shown in Figure 2. It can be seen that the presignal noise of the denoised signal by PD method is basically closer to zero, while the others still have some noise. PD method has a better performance in restoring the original signal compared with EEMD and wavelet-based algorithms.Figure 2
Comparison of signal denoising result with different types.
(a)
Denoising result for synthetic sinusoidal signal (b)
Denoising result for synthetic double sinusoidal signal (c)
Denoising result for synthetic blast signal (d)
Denoising result for synthetic microseismic signalTable1 summarizes the mean absolute error and standard deviation for three methods. It can be seen that the mean absolute error and standard deviation of PD method are less than that of EEMD and wavelet methods, which indicates that PD method has the ability to adapt to different noise types like white and colored ones.Table 1
Comparison of signal denoising indicator with different types.
Signal type PD method EEMD method Wavelet method μ σ μ σ μ σ Sinusoidal signal 0.0258 0.0416 0.0618 0.0775 0.0475 0.0592 Double sinusoidal signal 0.0179 0.0286 0.0330 0.0408 0.0332 0.0422 Blast signal 0.0083 0.0153 0.0316 0.0708 0.0274 0.0588 Microseismic signal 0.0075 0.0140 0.0278 0.0507 0.0184 0.0284


### 3.2. Results for Signals with Different SNRs

SNR is a crucial indicator of signal quality, defined as follows:(13)SNR=10log⁡∑i=1Nsi2∑i=1Ndi2The lower the SNR is, the more difficult to distinguish signal from noise is. Therefore, a signal with low SNR requires enhancement prior to the main process. In order to compare the performance of PD algorithm for low SNR signal, we add different level white Gaussian noise to the synthetic signals. 16 examples with different SNRs ranging from -6db to 10db are generated and denoised by PD, EEMD, and wavelet method, as shown in Figure3. The mean absolute error and standard deviation are summarized in Table 2.Table 2
Comparison of signal denoising indicator with different SNR.
SNR PD method EEMD method Wavelet method μ σ μ σ μ σ -6 0.0163 0.0282 0.0851 0.1073 0.0530 0.0738 -5 0.0149 0.0263 0.0760 0.0956 0.0423 0.0593 -4 0.0132 0.0245 0.0677 0.0852 0.0359 0.0500 -3 0.0128 0.0233 0.0410 0.0656 0.0306 0.0428 -2 0.0120 0.0215 0.0308 0.0585 0.0269 0.0375 -1 0.0099 0.0187 0.0355 0.0616 0.0233 0.0330 1 0.0084 0.0162 0.0383 0.0479 0.0203 0.0309 2 0.0080 0.0162 0.0245 0.0559 0.0186 0.0292 3 0.0066 0.0125 0.0278 0.0570 0.0185 0.0292 4 0.0054 0.0099 0.0256 0.0556 0.0159 0.0266 5 0.0048 0.0090 0.0248 0.0519 0.0149 0.0260 6 0.0044 0.0077 0.0215 0.0269 0.0147 0.0260 7 0.0037 0.0070 0.0236 0.0527 0.0144 0.0258 8 0.0036 0.0066 0.0171 0.0214 0.0130 0.0241 9 0.0030 0.0057 0.0222 0.0524 0.0127 0.0248 10 0.0028 0.0054 0.0212 0.0522 0.0124 0.0244Figure 3
Comparison of signal denoising result with different SNR by PD method.It can be seen that the mean absolute error and the standard deviation of PD, EEMD, and wavelet denoising increase with decrease of SNR (Table2). When SNR reaches -4 or lower, P- and S-wave have been almost covered by the noise; PD method can still suppress the noise effectively and enhance the P- and S-wave phase. When SNR is -6, the mean absolute error of PD method is 0.0132, which is relatively smaller than that of the EEMD algorithm (0.0677) and the wavelet method (0.0359). We can conclude that PD method can perform better in low SNR cases.


### 3.3. Results for Signals with Wrong Autopicking Arrivals

Our proposed method is fully automated and has not any parameter to tune, which is more suitable for automatic processing in microseismic monitoring. However, it is well known that no autopicker is 100% accurate especially in presence of noise, as well as our modified AIC picker used in PD method. That is to say, the picking error is always existing and can not be eliminated. In order to examine the influence of picking error caused by our modified AIC picker, we use the synthetic signal with the lowest SNR (-6db) generated in Section3.2 for experiments. The length of this seismogram is 3000 points, and the actual P-wave arrival is 1001 points. We assume that the picking P-wave arrival is ranging from the 601 to 1401 points and then compare the denoising results between the right P-wave arrival input and the wrong P-wave arrival input. The mean absolute percentage error is used to exhibit the deviations, defined as follows:(14)ε=1N∑i=1NDwi-DriDriwhere Dw and Dr are the denoised signals with wrong and right P-wave arrival input, respectively;N is the number of points of the whole signal.The results are summarized in Table3. We can observe that, as the picking error increases, the absolute percentage error increases. If the difference between the autopicking arrival and the actual arrival is less than 100 points, the absolute percentage error is less than 10%. As suggested by Leonard [28], for many applications, the time estimated by AIC picker can be used without review by a seismic analyst, with minor differences between analysts. Therefore, we can conclude that our proposed method is applicable for denoising of most microseismic data. However, it should also be noted that our method may not have a perfect denoising result for extreme seismograms with autopicking arrivals far away from the actual arrival.Table 3
Influence of picking error.
Picking error in samples -400 -200 -100 -50 0 50 100 200 400 ε 19.8% 18.2% 8.8% 6.5% 0 8.1% 9.6% 27.8% 46.8%


## 3.1. Results for Signals with Different Types

We first apply the algorithm to synthetic microseismic signals to test its performance. The procedure for creating synthetic signals is based on the work of Halldorsson and Papageorgiou [25]. A noisy sample was generated by adding white Gaussian noise to a synthetic signal. We have tested the algorithm for synthetic examples with different noise level and compared it with the EEMD denoising [26] and wavelet denoising method [27], which are widely used in noise reduction of microseismic signals.We use the mean absolute errorμ and the standard deviation σ as the evaluation criteria of the denoising performance. μ represents the extent of discrimination between the denoised signal and the real signal. The smaller the value is, the more similar the denoised signal is to the real signal. σ is used to quantify the amount of variation or dispersion. A low standard deviation indicates that the data points tend to be close to the mean value, while a high standard deviation indicates that the data points are spread out over a wider range of values, defined as follows:(11)μ=1N∑i=1Ndi-si(12)σ=1N-1∑i=1Ndi-si-λ2where λ=1/N∑i=1Ndi-si,d is the denoising signal, s is the real signal, and N is the number of signal points.Four types of signals that include presignal noise are synthesized: sinusoidal signal, double sinusoidal signal, blast vibration signal, and microseismic signal. The amplitude of the signal with additive white Gaussian noise is normalized to±1. After that, we denoise the signal by PD, EEMD, and wavelet denoising method. Results are shown in Figure 2. It can be seen that the presignal noise of the denoised signal by PD method is basically closer to zero, while the others still have some noise. PD method has a better performance in restoring the original signal compared with EEMD and wavelet-based algorithms.Figure 2
Comparison of signal denoising result with different types.
(a)
Denoising result for synthetic sinusoidal signal (b)
Denoising result for synthetic double sinusoidal signal (c)
Denoising result for synthetic blast signal (d)
Denoising result for synthetic microseismic signalTable1 summarizes the mean absolute error and standard deviation for three methods. It can be seen that the mean absolute error and standard deviation of PD method are less than that of EEMD and wavelet methods, which indicates that PD method has the ability to adapt to different noise types like white and colored ones.Table 1
Comparison of signal denoising indicator with different types.
Signal type PD method EEMD method Wavelet method μ σ μ σ μ σ Sinusoidal signal 0.0258 0.0416 0.0618 0.0775 0.0475 0.0592 Double sinusoidal signal 0.0179 0.0286 0.0330 0.0408 0.0332 0.0422 Blast signal 0.0083 0.0153 0.0316 0.0708 0.0274 0.0588 Microseismic signal 0.0075 0.0140 0.0278 0.0507 0.0184 0.0284


## 3.2. Results for Signals with Different SNRs

SNR is a crucial indicator of signal quality, defined as follows:(13)SNR=10log⁡∑i=1Nsi2∑i=1Ndi2The lower the SNR is, the more difficult to distinguish signal from noise is. Therefore, a signal with low SNR requires enhancement prior to the main process. In order to compare the performance of PD algorithm for low SNR signal, we add different level white Gaussian noise to the synthetic signals. 16 examples with different SNRs ranging from -6db to 10db are generated and denoised by PD, EEMD, and wavelet method, as shown in Figure3. The mean absolute error and standard deviation are summarized in Table 2.Table 2
Comparison of signal denoising indicator with different SNR.
SNR PD method EEMD method Wavelet method μ σ μ σ μ σ -6 0.0163 0.0282 0.0851 0.1073 0.0530 0.0738 -5 0.0149 0.0263 0.0760 0.0956 0.0423 0.0593 -4 0.0132 0.0245 0.0677 0.0852 0.0359 0.0500 -3 0.0128 0.0233 0.0410 0.0656 0.0306 0.0428 -2 0.0120 0.0215 0.0308 0.0585 0.0269 0.0375 -1 0.0099 0.0187 0.0355 0.0616 0.0233 0.0330 1 0.0084 0.0162 0.0383 0.0479 0.0203 0.0309 2 0.0080 0.0162 0.0245 0.0559 0.0186 0.0292 3 0.0066 0.0125 0.0278 0.0570 0.0185 0.0292 4 0.0054 0.0099 0.0256 0.0556 0.0159 0.0266 5 0.0048 0.0090 0.0248 0.0519 0.0149 0.0260 6 0.0044 0.0077 0.0215 0.0269 0.0147 0.0260 7 0.0037 0.0070 0.0236 0.0527 0.0144 0.0258 8 0.0036 0.0066 0.0171 0.0214 0.0130 0.0241 9 0.0030 0.0057 0.0222 0.0524 0.0127 0.0248 10 0.0028 0.0054 0.0212 0.0522 0.0124 0.0244Figure 3
Comparison of signal denoising result with different SNR by PD method.It can be seen that the mean absolute error and the standard deviation of PD, EEMD, and wavelet denoising increase with decrease of SNR (Table2). When SNR reaches -4 or lower, P- and S-wave have been almost covered by the noise; PD method can still suppress the noise effectively and enhance the P- and S-wave phase. When SNR is -6, the mean absolute error of PD method is 0.0132, which is relatively smaller than that of the EEMD algorithm (0.0677) and the wavelet method (0.0359). We can conclude that PD method can perform better in low SNR cases.


## 3.3. Results for Signals with Wrong Autopicking Arrivals

Our proposed method is fully automated and has not any parameter to tune, which is more suitable for automatic processing in microseismic monitoring. However, it is well known that no autopicker is 100% accurate especially in presence of noise, as well as our modified AIC picker used in PD method. That is to say, the picking error is always existing and can not be eliminated. In order to examine the influence of picking error caused by our modified AIC picker, we use the synthetic signal with the lowest SNR (-6db) generated in Section3.2 for experiments. The length of this seismogram is 3000 points, and the actual P-wave arrival is 1001 points. We assume that the picking P-wave arrival is ranging from the 601 to 1401 points and then compare the denoising results between the right P-wave arrival input and the wrong P-wave arrival input. The mean absolute percentage error is used to exhibit the deviations, defined as follows:(14)ε=1N∑i=1NDwi-DriDriwhere Dw and Dr are the denoised signals with wrong and right P-wave arrival input, respectively;N is the number of points of the whole signal.The results are summarized in Table3. We can observe that, as the picking error increases, the absolute percentage error increases. If the difference between the autopicking arrival and the actual arrival is less than 100 points, the absolute percentage error is less than 10%. As suggested by Leonard [28], for many applications, the time estimated by AIC picker can be used without review by a seismic analyst, with minor differences between analysts. Therefore, we can conclude that our proposed method is applicable for denoising of most microseismic data. However, it should also be noted that our method may not have a perfect denoising result for extreme seismograms with autopicking arrivals far away from the actual arrival.Table 3
Influence of picking error.
Picking error in samples -400 -200 -100 -50 0 50 100 200 400 ε 19.8% 18.2% 8.8% 6.5% 0 8.1% 9.6% 27.8% 46.8%


## 4. Field Data Example

In this section, we demonstrate the performance of the proposed method using a field dataset recorded at Dongguashan Copper Mine located in southern Anhui Province, China. Dongguashan is famous in China for its large and deep porphyry copper deposits. Although mining at Dongguashan currently takes place at depths ranging from 730m to 875m, its mining depth was extended to 1000m underground. Dongguashan was one of the first mines in China to deploy a microseismic monitoring system for rockburst management, in 2004 [2]. The system was abandoned since 2013. Recently, a newly microseismic monitoring system was installed in October 2017. Six accelerators with a sensitivity of 10V/g were embedded in -730m level and -790m level, as shown in Figure 4. The sampling frequency was set to 10KHz.Figure 4
Layout of six accelerators embedded in Dongguashan Copper Mine (a gray inverted cone is used to indicate the accelerator).Our dataset is composed of 348 highly noised seismograms associated with 58 microseismic events occurred during December 8~10, 2017. As the original signal cannot be restored, SNR is unable to calculate using (13). A new indicator is introduced to illustrate the quality of the signal, which is defined as follows:(15)PSNR=10log⁡∑i=pp+lxi2∑i=p-lpxi2where p is the exact P-wave arrival and l is the window length in sampling points.PSNR represents the discrimination between signal and background noise nearby P-wave arrival.PD, EEMD, and wavelet denoising methods were applied to the dataset. The averagePSNR after denoising with these methods was compared with the original one. Results are summarized in Table 4. The averagePSNR is improved from 5.8010 to 25.1369 after PD method, which is better than 14.8202 of EEMD method and 17.3983 of wavelet method.Table 4
Mean P-wave SNR.
Type Original signals After EEMD denoising After wavelet denoising After PD denoising Mean PSNR 5.8010 14.8202 17.3983 25.1369The improvement of SNR has a great influence on the subsequent microseismic data processing (e.g., arrival picking and event location). The accuracy of P-wave arrival is pivotal for event location and source parameters calculation. For high SNR signal, automatic algorithm can pick the first-arrival accurately. For low SNR signal, automatic picking may have a large deviation.The SNR of microseismic signal can significantly be improved by applying our proposed method, which will benefit the accuracy of automatic P- and S-wave picking. As shown in Figures5(a) and 5(b), P-wave arrival is automatically picked by STA/LTA method, resulting in 192 sampling points error (19.2 ms) compared with manual picking. As shown in Figures 5(c) and 5(d), the signal quality is improved obviously after PD algorithm denoising, and the P-wave arrival is picked again by STA/LTA method, resulting in 6 sampling points error (0.6 ms) compared with manual picking.Figure 5
Different P-wave autopicking result before and after PD denoising.
(a)
Comparison between actual P-arrival and P-arrival picked by STA/LTA on noisy microseismic signal (b)
Zoomed windows of (a) around P-arrival (c)
P-arrival picked by STA/LTA on microseismic signal denoised by PD method (d)
Zoomed windows of (c) around P-arrival


## 5. Conclusion

We have introduced a simple nonparametric denoising method based on two-step AIC and Fourier transform. The method consists of three major steps: (1) applying a two-step AIC algorithm to pick P-wave arrival; (2) subtracting the noise power spectrum from the signal power spectrum; (3) recovering the microseismic signal by inverse Fourier transform. Our proposed denoising method has been applied to both synthetic and real microseismic data. Results show that the new method outperforms EEMD and wavelet filtering methods and considerably improves the SNR of the waveforms.



---

*Source: 4367201-2018-07-16.xml*