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.
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 length N is defined as where is the kth sampling point, ; is the order of the autoregressive process; is a constant; and 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:where var represents the variance of the data, . 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 signal , , the Fourier transform and its inverse transform can be expressed as
Because 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 signal of length (M is a positive integer; zero-padding is used if the signal length is insufficient), we split the N-point data sequence into two N/2-point data sequences and , corresponding to the even-numbered and odd-numbered samples of , respectively; that is,
Thus and are obtained by decimating by a factor of 2, and the above decimation-in-time decomposition can be performed -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:where is the rotation factor, is the DFT of the odd sequence, and 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 Figure 1. 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.
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)Let 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 .(2)For time interval , apply the AIC method and denote the primitive P-wave arrival as .(3)For time interval , apply the AIC method again and denote the corrected P-wave arrival as . Thus we get the presignal , denoted by .(4)Add Hamming window and apply Fast Fourier transform to and , respectively, then we get and .(5)Calculate the average amplitude and phase angle of ; calculate the average amplitude of . 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,(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.
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:where , d is the denoising signal, is the real signal, and 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 . 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.
(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 signal
Table 1 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.
3.2. Results for Signals with Different SNRs
SNR is a crucial indicator of signal quality, defined as follows:
The 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 Figure 3. The mean absolute error and standard deviation are summarized in Table 2.
It can be seen that the mean absolute error and the standard deviation of PD, EEMD, and wavelet denoising increase with decrease of SNR (Table 2). 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 Section 3.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:where and 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 Table 3. 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.
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.
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:where is the exact P-wave arrival and 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 average PSNR after denoising with these methods was compared with the original one. Results are summarized in Table 4. The average PSNR 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.
The 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 Figures 5(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.
(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.
Data Availability
The data used to support the findings of this study were provided by Dongguashan Copper Mine under license and so cannot be made freely available. Access to these data will be considered by the author upon request, with permission of Dongguashan Copper Mine.
Conflicts of Interest
On behalf of all authors, the corresponding author states that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
This work was supported by the National Key R&D Program of China (no. 2017YFC0602905) and the Fundamental Research Funds for the Central Universities of Central South University (no. 2015zzts073).