Three-Component Microseismic Data Denoising Based on Re-Constrain Variational Mode Decomposition

: Microseismic monitoring is an important technology used to evaluate hydraulic fracturing, and denoising is a crucial processing step. Analyses of the characteristics of acquired three-component microseismic data have indicated that the vertical component has a higher signal-to-noise ratio (SNR) than the two horizontal components. Therefore, we propose a new denoising method for three-component microseismic data using re-constrain variational mode decomposition (VMD). In this method, it is assumed that there is a linear relationship between the modes with the same center frequency among the VMD results of the three-component data. Then, the decomposition result of the vertical component is used as a constraint to the whole denoising effect of the three-component data. On the basis of VMD, we add a constraint condition to form the re-constrain VMD, and deduce the corresponding solution process. According to the synthesis data analysis, the proposed method can not only improve the SNR level of three-component records, it also improves the accuracy of polarization analysis. The proposed method also achieved a satisfactory effect for ﬁeld data.


Introduction
Microseismic monitoring is a technology used to monitor the subsurface fracturing of rocks due to human activities [1][2][3][4]. This is the first step in recognizing microseismic event signals for microseismic monitoring. However, owing to the characteristics of microseismic event signals (e.g., low energy, high frequency, and short duration) they are easily disrupted by various factors, which increases the difficulty of the subsequent processing and interpretation. Therefore, it is necessary to denoise the monitoring record.
The denoising process for microseismic data is mainly divided into two basic strategies according to the type of data: multi-channel denoising and single-channel denoising [5]. For the multi-channel strategy, microseismic data are denoised by the spatial distribution information of the geophone. Gan et al. [6] used a median filter based on the seismic profile structure to suppress blending noise. Chao et al. [7] used a method based on the correlation between the three components of microseismic to identify the 3D shearlet transform coefficient and effectively remove the noise. Bai et al. [8] established a model based on the least squares method to decompose seismic data using the spatial and temporal relationships among multi-channel seismic data. The single-channel strategy is mainly used for denoising via time-frequency transformation and signal decomposition. Chakrabort et al. [9] used wavelet time-frequency analysis to suppress the seismic noise. Mousavi et al. [10] combined the synchrosqueezed continuous wavelet transform and a detection function to remove noise from the signal.
Noise reduction methods based on decomposition and reconstruction are widely used in single-channel noise reduction [11], such as empirical mode decomposition (EMD) [12,13] which is used to process the real and imaginary parts of the frequency domain of the signal in the time window to reduce random and coherent noise. Chen et al. [14] applied an AR mode to process the results of EMD in the frequency domain to reduce random noise. However, EMD has some drawbacks in practical applications, such as mode mixing [15]. Han et al. [16] used ensemble empirical mode decomposition (EEMD), which was proposed by Wu et al. [17] to overcome the problem of mode mixing and to process the real and imaginary parts of the frequency domain of the signal in the time window to reduce the noise. Chen et al. [18] applied wavelet threshold filtering to process the high-frequency intrinsic mode function (IMF) of the EEMD, and then reconstructed the high frequency and low frequency to reduce the noise. Dong et al. [19] combined the complex curvelet transform and complementary ensemble empirical mode decomposition (CEEMD) which was proposed by Yeh et al. [20] to improve the effect of decomposition. Zuo et al. [21] selected the IMFs of the CEEMD using the self-correlation coefficient, then used the wavelet packet threshold to process the IMFs, and finally reconstructed the IMFs to suppress the noise. Peng et al. [22] selected the EEMD mode by the variance contribution rates (VCRs), which removed the modes with VCR < 0.01 and others retained. Then, each retained mode was constructed as a Hankel matrix, and it was processed by principal component analysis (PCA) to complete denoising.
Compared with EMD, EEMD, and CEEMD, VMD can further decrease the redundant modes and solve the mode-mixing problem, and has a more powerful anti-noise performance [23]. Since then, many studies have applied VMD to reduce noise. Liu et al. [24] used VMD to perform time-frequency analysis of seismic data to suppress noise and highlight the geologic characteristics. Li et al. [25] selected IMFs to reconstruct a signal based on detrended fluctuation analysis (DFA). Huang et al. [26] performed a correlation analysis on the IMFs of the VMD to suppress the noise. Zhou et al. [27] combined VMD and odd spectrum analysis to remove the remaining low-frequency noise of the VMD. Li et al. [28] applied time-frequency peak filtering for the IMFs of VMD and reconstructed the signal to reduce noise. Liu et al. [29] used VMD to suppress ground rolling waves to achieve denoising. Zhang et al. [30] used the Akaike information criterion to judge the results of the VMD and selected the threshold to suppress the noise. In the above application research, most scholars have focused on processing IMFs to achieve a better denoising effect while using VMD. Li et al. [31] applied DFA to determine whether the VMD modes of the water inrush signals were random noise or a valid signal, and eliminated the noise in the mode using DFA.
Analyses of the characteristics of the acquired microseismic data found that the Pwave of the vertical component has a higher SNR than the horizontal components, which is not conducive to polarization analysis [32,33]. Therefore, a noise reduction method for the horizontal component is needed. When the three-component microseismic data are processed by VMD, the horizontal component is more easily affected by noise of similar frequency than the vertical component. This is because the level of random noise of the horizontal component is higher than that of the vertical component. Rodriguez et al. used a redundant dictionary to denoise the three components of microseismic data based on the sparse coefficient distribution of the three components [34]. Therefore, to improve the SNR of the horizontal component, it is assumed that each component consists of multiple IMFs, and there is a linear relationship between the modes of the same center frequency to limit the P-wave of the horizontal component by the vertical component. In this study, we propose a new method that adds a new constraint to the original formulas of VMD.
The proposed method simultaneously processes three-component microseismic data based on the VMD. However, the formulas for processing are different between the vertical component and the horizontal components. Based on the linear relationship among the three components of microseismic data, and the non-linear relationship of the random noise in the three-component data, the vertical component is processed by the original formula of VMD, and the horizontal components are processed by the re-constrained formula of VMD.
The aim of this processing is to make the signal amplitude distribution of the horizontal components consistent with that of the vertical component as much as possible.
In the verification phase of the proposed method, a synthetic three-component data is designed in which the SNR horizontal component is relatively low, and that of the vertical component is relatively high. After denoising, the polarization characteristics of the data were analyzed to validate the effectiveness of the proposed method. In the practical phase, field microseismic data were processed based on the proposed method. The corresponding results show that the proposed method can effectively improve the SNR of three-component data as well as the precision of the polarization information.

Variational Mode Decomposition
VMD is a signal decomposition method in which the IMF and corresponding center frequencies are obtained through iteration. The method first assumes that the original signal is decomposed into k mode components. It then takes the minimum Gaussian smoothness of each mode component as the goal. The sum of the modes is equal to the decomposed signal as the constraint for optimization. Finally, the decomposition results are obtained. The VMD expression is defined as follows: where K is the number of decomposition modes, u k is the kth IMF, ω k is the center frequency of the kth mode, f is the original signal, * is the sign of convolution, and δ(t) is the Dirac delta function which indicates a unit impulse symbol.
The Lagrange multiplier method is used to solve the problem. The Lagrange multiplication operator is introduced to transform the constrained variational problem into an unconstrained variational problem. The updated formula is derived as follows: where α is the penalty factor, f(ω) is the Fourier transform of the original signal, ω k is the center frequency of the kth IMF of the horizontal component, u k is the Fourier transform of the kth IMF, and λ(ω) is the Fourier transform of the Lagrangian multipliers.

Re-Constraint of Variational Mode Decomposition
In the ideal state, the relationship between the three components of the microseismic data can be regarded as a linear relationship: where H 1 and H 2 are horizontal components, V is the vertical component, and a and b are coefficients that represent the linear relationships between the components. Since the center frequencies of the P-wave and S-wave of the same microseismic event are different, this study assumed a linear relationship for the modes of the same (or close to) center frequencies in the VMD decomposition results of different components. Based on this relationship, we re-constrained the VMD and obtained a new formula that adds a new constraint condition: where v k is the kth IMF after the decomposition of the vertical component, u k is the kth IMF after the decomposition of the horizontal component, and a k is a constant that expresses the linear relationship between the kth IMFs. The third line of Formula (4) is the re-constrained condition, indicating that there is a linear relationship between the modes.
By using the Lagrange multiplier method, we get the corresponding updated formulas: where α is the penalty factor, f(ω) is the Fourier transform of the original signal, ω k is the center frequency of the kth IMF of the horizontal component, and u n+1 k is the Fourier transform of the kth IMF. λ(ω) and µ(ω) are the Fourier transforms of the Lagrangian multipliers.
Under the noise-free condition, the VMD results for three-component data share the same frequency distributions. However, owing to the influence of various factors, the frequency distributions of the VMD results will be different in field data. Because the characteristic of the SNR of the vertical component is better than that of the horizontal component, we updated the horizontal component IMFs by using the center frequency of the vertical component IMFs in the iteration processing. The premise of using this method to process three-component microseismic data is that the SNR of the vertical component is higher than that of the horizontal components. Thus, we obtained the new updated formulas: where α is the penalty factor, f(ω) is the Fourier transform of the original signal, ω zk is the center frequency of the kth IMF of the vertical component, and u n+1 k is the Fourier transform of the kth IMF. λ(ω) and µ(ω) are the Fourier transforms of the Lagrangian multipliers.
As shown in Figure 1, the processing flow of the re-constrain VMD is as follows: 1.
Load the initial three-component microseismic data.

2.
Expend the data using a mirror image against the boundary effect, and perform the Fourier transform for the original data. The related parameters of the subsequent processing are initialized.

3.
Update the vertical component variable using VMD updated formulas as Formula (2).

4.
Update the parameters of the horizontal components using updated re-constrain VMD formula as in Formula (6).

5.
Judge whether the end condition of iteration is satisfied. The end condition is determined by the maximum number of iterations. If the condition is not satisfied, repeat Step 3.

6.
Perform the inverse Fourier transform for the previous output, and remove the mirror data. Then, the final result is obtained.

Synthetic Data
We used an attenuated simple harmonic to verify the effectiveness of this method. Figure 2 shows pure attenuated harmonic three-component data with frequencies of 30 Hz and 70 Hz. Figure 3 shows waves contaminated with random noise. The SNR of the horizontal component is set to −10 dB, and the SNR of the vertical component is set to 5 dB. In this study, we used the SNR, which is defined as: where x(t) is the pure synthetic signal, n(t) is random noise, and N is the number of sampling points. Figure 4 shows the VMD results of three-component data. Since the parameter k of VMD only impacts the frequency representation performance, it has no noticeable impact on the reconstruction performance [24]. In addition, the proposed method only focuses on the frequency representation of the vertical component, which is viewed as a re-constraint condition. Therefore, the parameter k for each component was set to the same value in this study. Figure 5 shows the decomposition results of the three-component synthetic data obtained by the re-constrain VMD. As shown in Figure 5, the redundant amplitude of the nonlinear relationship information was removed. The remaining modes were also suppressed because they are all random signals and have no direct correlation.
To prevent the re-constrained modes of two horizontal components (of which the SNR is low) from being an illusion of noise fitting (in other words, noise is regarded as an effective component of the signal to participate in the reconstruction processing), we removed the pure synthetic signal and performed the re-constrain VMD for these two components. Figure 6 shows the corresponding results. For the re-constrained modes in the two horizontal components, there was a weak correlation with the vertical component after re-constrained decomposition, indicating that the re-constrained modes of the two horizontal components were also affected by random noise, and the degree of influence was quite limited.   The first and second modes were reconstructed as a result of simple noise reduction. In Figure 7, it can be seen that compared with the VMD results, the results of re-constrain VMD could obtain better polarization characteristics and remove more background noise.
To better evaluate its denoising effect, we analyzed the polarization characteristics of the 30 Hz signal in the three-component data. The length of the selected time window was an entire period, and the 230th sampling point was manually selected as the starting point of this window, as shown in Figure 8.      In the synthetic data in Figure 8, the two red lines indicate the range of the selected windows for polar analysis. It is evident that the polarization directions were changed for the synthetic data contaminated with random noise. After denoising based on VMD and re-constrain VMD, the accuracy of the calculated polarization directions was greatly improved. In addition, the results of hodograms based on re-constrain VMD were more compact than those based on VMD.

Field Microseismic Data
The field three-component microseismic event data were acquired in the hydraulic fracturing of a shale gas reservoir. The original waveforms of this microseismic event are shown in the first row of Figure 9. These microseismic data include the P-wave and S-wave of the microseismic event. The range of the P-wave is approximately from the 1000th sampling point to the 1450th sampling point, and the range of S-wave is approximately from the 1450th sampling point to the 2000th sampling point. Because of the formation mechanism of microseismic events, there is a partial overlap between the P-wave and S-wave. The SNR of the vertical component is better than that of the two horizontal components. The noise reduction results obtained using VMD and re-constrain VMD are shown. Comparing the waveforms between the original data and denoising results, the random noise was effectively suppressed by both methods. As shown in Figure 9, these two methods could suppress random noise to a certain extent. According to the original waveform of the horizontal component H1, the SNR of the P-wave is relatively low. Comparing the H1 components of the two groups' denoising results, it can be seen that the P-wave waveforms of microseismic data were not well recovered using VMD. Correspondingly, the waveforms of the P-wave were more obvious after using re-constrain VMD. It can also be seen that the background noise was better suppressed after using re-constrain VMD, but this effect was controlled by the background noise intensity of the vertical component.
Because the P-wave is commonly used to analyze the polarization direction of microseismic events, we selected from the 1000th point to the 1070th point as the time window range, which is indicated by the red lines in Figure 9. The data in the time window were used to analyze the polarization and to judge the denoising effect. After denoising using VMD and re-constrain VMD, the polarization characteristics of the denoising results are shown in Figures 10 and 11.  In this study, the enhancement of polarization information by denoising was evaluated by determining whether the polarization characteristic of the curve of the hodogram was clear, compact, and self-consistent. It was largely straightforward to determine whether the polarization characteristic of the curve was clear, except for the H1-H2 hodogram in Figure 10. Regarding compactness, all hodograms in Figure 11 are better than those in Figure 10. Regarding self-consistency, according to the positive or negative correlation among the three components, the hodograms in Figure 11 are better than those in Figure 10.

Conclusions
The proposed method mainly focuses on solving the challenge of low-SNR horizontal components of microseismic data. Among the different components, there is no correlation for random noise, and there is a linear relationship for the signal. For most field threecomponent microseismic data, the SNR of the vertical component is better than that of the two horizontal components. Therefore, by adding a new constraint to the VMD and modifying the new updated formulas, we propose the re-constrain VMD. The new constraint is used to represent the characteristics of linear relationships for signal and random noise among the three components of the data.
According to the denoising results for the synthetic data and field microseismic data, we found that the proposed method could suppress random noise to a certain extent, and the effect was better than that of VMD. The subsequent polarization analysis also showed that the proposed method could obtain better polarization characteristics.