Noise suppression of microseismic data based on a fast singular value decomposition algorithm

Journal of Applied Geophysics 170 (2019) 103831

Contents lists available at ScienceDirect

Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo

Noise suppression of microseismic data based on a fast singular value decomposition algorithm Hui Lv a,b,⁎ a b

College of Civil Engineering and Architecture, Zhejiang University, 38 Zheda Road, Hangzhou, Zhejiang Province 310027, China School of Civil Engineering and Architecture, Nanchang Hangkong University, 696 South Fenghe Avenue, Nanchang, Jiangxi Province 330063, China

a r t i c l e

i n f o

Article history: Received 23 March 2019 Received in revised form 30 July 2019 Accepted 30 July 2019 Available online 09 August 2019

a b s t r a c t Raw data from microseismic monitoring is usually contaminated with ambient noise, which severely decreases the quality of the data and further affects the accuracy of arrival-based or amplitude-based source imaging performance., and thusThus, noise suppression of microseismic data plays an indispensable role in event detection during hydraulic fracturing. We investigatedeveloped a new way to effectively suppress noise from microseismic data via cross-correlation of traces from multiple microseismic channels. The cross-correlation between a reference trace and all available traces is computed to estimated the optimal time shifts between the waveforms in different channels. When the time shifts are estimated, the microseismic data can be rearranged to enhance the spatial correlation of the microseismic events. Next, a singular value decomposition (SVD) step is applied to extract the eigen-components or eigen-images of the noisy traces and to reject the noise in the data. When neighbor waveforms are too close, a windowed processing strategy is required to obtain well-flattened gather for the SVD processing. Considering that the SVD calculation is computationally expensive, we propose a fast SVD decomposition algorithm by derivation. The proposed method is easy to implement and can obtain robust performance in denoising multi-channel microseismic dataset. We validate the performance on both synthetic and real microseismic datasets. The proposed method can also be combined with active-source seismic data to efficiently suppress strong random background noise. © 2019 Published by Elsevier B.V.

1. Introduction Microseismic monitoring is useful in characterizing the fracturing process when developing unconventional reservoirs (Rigo et al., 1996; Friedrich et al., 1998; Maxwell et al., 2002; Fisher et al., 2004; Mayerhofer et al., 2006; Warpinski et al., 2009; Šílený et al., 2009; Reshetnikov et al., 2010; Flewelling et al., 2013; Hajati et al., 2015; Hansen and Schmandt, 2015; Mousavi and Langston, 2016; Grigoli et al., 2016). The fracture distribution can be inverted from observed microseismic data using the wave theory. The source localization is one of most important tasks in the process of microseismic monitoring, which has attracted attention from both academia and industry (House, 1987; Rutledge et al., 1998; Flewelling et al., 2013). Theoretical study of the source localization problem is significant for research development of seismology and rock mechanics (Dong and Li, 2013; Dong et al., 2014, 2016a, 2016b; Rutledge et al., 2004; Shapiro et al., 2006). However, the small amplitude of most microseismic signals make the effective ⁎ Corresponding author at: School of Civil Engineering and Architecture, Nanchang Hangkong University, 696 South Fenghe Avenue, Nanchang, Jiangxi Province 330063, China. E-mail address: [email protected].

https://doi.org/10.1016/j.jappgeo.2019.103831 0926-9851/© 2019 Published by Elsevier B.V.

usage of the observed data challenging as in many cases the signal-tonoise ratio (SNR) of the data is extremely low, which causes difficulties for many conventional seismic processing techniques (Gan et al., 2016c; Siahsar et al., 2017; Zu et al., 2019). For example, a low SNR may result in inaccurately picked arrivals and negatively affect the microseismic imaging performance. To fully take advantage of the benefits from microseismic monitoring technique, it is crucial to enhance the microseismic waveforms in many the microseismic monitoring applications (Pearson, 1981; Phillips et al., 1998; Fehler et al., 1998; Bolton and Masters, 2001; Chouet et al., 2003; Reyes-Montes et al., 2005; Capilla, 2006; Fischer et al., 2008; Šílený et al., 2009; Wu and Huang, 2009; Baig and Urbancic, 2010; Maxwell et al., 2010; Han and van der Baan, 2015; Chen, 2016; Li et al., 2016a; Huang et al., 2017b; Chen et al., 2017; Zhang et al., 2017; Chen, 2018; Huang et al., 2018). Single-channel processing strategy is widely used in the microseismic data processing community as most microseismic data are sparsely distributed in space. Decomposition and reconstruction is a commonly used method to denoise microseismic datasets. The principle of this type of methods is to decompose a noisy microseismic trace into several components with different oscillating frequencies and then reconstruct the microseismic trace by carefully select the dominant components based on some criteria. Among these methods, the empirical mode

2

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

Table 1 Comparison of computing time with different model sizes between the SVD and fast SVD calculations. Model sizes

SVD (s)

FSVD (s)

Speedup

1000×8 5000×8 5000×40 10,000×40

0.036 0.487 1.832 8.123

0.00163 0.00287 0.0231 0.069

22.086 169.686 79.307 117.725

decomposition (EMD) (Huang et al., 1998; Wu and Huang, 2009; Chen, 2016; Chen et al., 2017; Zhang et al., 2017) is usually used to decompose microseismic datasets. The benefit of the EMD method is that the decomposition is fully automatic, meaning that no predefined parameters need to be set before the decomposition. This adaptivity makes EMD very convenient to use for massive data processing. More advanced decomposition is the ensemble empirical mode decomposition (EEMD) method, which is an improved variant of EMD. The EEMD method solves the mode-mixing problem of the traditional EMD method by introducing well-controlled Gaussian white noise. Han and van der Baan (2015) combined the EEMD method and block thresholding strategy to denoise microseismic data and obtain satisfying results. However, EEMD has its inherent drawbacks, e.g., the reconstructed signal deviate from the clean signal seriously. EEMD was then improved into the complete EEMD method and was recently replaced by the variational mode decomposition (VMD) method (Dragomiretskiy and Zosso, 2014; Liu et al., 2016a, 2017, 2018), which enjoys the benefits of EMD method but with better control on the behavior. A recently popular EMD-like method is the morphological decomposition method, developed by Li

et al. (2016a) and Li et al. (2016b). The morphological decomposition method is based on the theory of mathematical morphology and decompose an input 1D signal into an ensemble of components according to the morphological scale. Li et al. (2016a) applied the morphological decomposition method to suppress the strong low-frequency noise by carefully designing a filter bank according the target frequency of the noise component. The morphological decomposition method is more computationally efficient than the EMD method but requires some empirical knowledge in designing the scale operator. Huang et al. (2017b) improved the morphological decomposition method by developing a more adaptively filter bank designing strategy. They also applied the morphological decomposition method to active-source seismic data processing. Huang et al. (2018) combined the local orthogonalization (Chen and Fomel, 2015) method with the morphological decomposition method to denoise microseismic dataset. This hybrid method is referred to as the orthogonalized morphological decomposition method, where the local orthogonalization method is used to increased the orthogonality among different morphological components. With the better orthogonality, the denoised waveform signals are easier to be picked by an automatic arrival picking algorithm and is beneficial for more accurate source location imaging (Huang et al., 2017a). More recently, Chen (2018) introduced another EMD-like but more mathematical solid decomposition method called least-squares complex decomposition (LSCD). The LSCD method decompose a microseismic trace by assuming the frequency spectrum of different components are smoothly variable and use the smoothness constraint to regularize the behaviors of the resulted components. In addition to the decomposition and reconstruction method, many other denoising algorithms that are well studied in the active-source seismic literature are also potential to achieve higher

Fig. 1. Synthetic data. (a) Noise-free microseismic data. (b) Noisy microseismic data (with SNR as −4.58 dB).

30 25 20

Shifts

Channel NO

15

10

15 10

5

5 0

50

100

150

200

Time (ms)

250

300

350

0

5

10

Channel NO

Fig. 2. Synthetic example. (a) The data after rearrangement. (b) Time shifts calculated from cross-correlation.

15

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

3

Fig. 3. Synthetic example. (a) Result from the new method. (b) Result from the empirical mode decomposition method. (c) Filtered noise by the new method. (d) Filtered noise by the empirical mode decomposition method. The SNRs of the proposed and empirical mode decomposition methods are 12.23 dB and −0.95 dB, respectively.

SNR for microseismic datasets (Chen and Ma, 2014; Gan et al., 2016a, 2016b; Xue et al., 2017; Chen and Song, 2018; Zhou and Han, 2018; Zhou et al., 2018). Most microseismic denoising algorithms cannot deal with extremely noisy data. In this paper, we propose an effective method for removing the strong noise in a multi-channel microseismic data. We use crosscorrelation between neighbor channels to find the optimal relative shift between the two neighbor traces. With the cross-correlation time shifts among all the traces, we can rearrange the multi-channel microseismic gather to enhance the spatial coherency. When the multichannel microseismic traces are aligned well, we proposedpropose to apply the singular value decomposition (SVD) method to obtain the rank-one approximation of the aligned microseismic waveforms.

However, the SVD operation is computationally intensive. In order to reduce the overburden of SVD calculation, we develop an efficient SVD algorithm because of the special shape of the microseismic data matrix. The SVD method is powerful in attenuating random noise, even when there is strong noise energy. When a multiplenumber of microseismic events exist, we first segment the microseismic gather into small time-space windows, and apply the cross-correlation and rearrangement in local segmented windows. We apply the proposed method in both synthetic and real microseismic datasets and obtain very appealing performance, especially when compared with traditional alternative methods. It is also worth noting that the proposed method can be combined with active-source seismic data to suppress strong random background noise.

15

Channel NO

Channel NO

15

10

5

5

0

10

50

100

150

200

Time (ms)

250

300

350

0

50

100

150

200

250

300

Time (ms)

Fig. 4. Synthetic example with erratic incoherent noise. (a) Noisy data. (b) Denoised data using the proposed method.

350

4

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

15

Channel NO

Channel NO

15

10

5

5

0

10

50

100

150

200

250

300

0

350

50

100

Time (ms)

150

200

250

300

350

Time (ms)

Fig. 5. Synthetic example. (a) Clean data with two events. (b) Noisy data with two events (with SNR as −1.29 dB).

2. Theory 2.1. Enhancing the energy correlation via crosscorrelation The cross-correlation between two 1D signals s and g are formulated as: cðτÞ ¼

∞ X

sðnÞ  g ðn þ τ Þ;

two signals to reach the maximal correlation. In this way, we can recursively shift the neighbor traces to make their waveform energy most correlative. We can use cross-correlation to estimate the relative time shifts between two similar traces. The time shifts can be used to rearrange the two corresponding traces to have better energy correlation. Getting the time shifts can be expressed as:

ð1Þ

n¼−∞

τn ¼ arg max

where c(τ) denotes the cross-correlation sequence. τ is the time lag. n is the time index of signal s and g. The cross-correlation sequence obtains its maximum when the two signals are mostly correlative. The time lag corresponding to the maximum indicates the optimal shift between the

τ

Channel NO

Channel NO

10

5

5

50

100

150

200

250

300

0

350

50

100

Time (ms)

150

200

250

300

350

250

300

350

Time (ms)

15

15

Channel NO

Channel NO

ð2Þ

15

10

10

5

0

r ðmÞdðm þ τ; nÞ;

m ¼ minð1−τ;1Þ

where τn is the time shift. More specifically, it represents the shift of the nth trace. Mw is the number of samples in the data. r denotes the

15

0

maxðM w ;M w −τ Þ X

10

5

50

100

150

200

Time (ms)

250

300

350

0

50

100

150

200

Time (ms)

Fig. 6. Synthetic example with two events. (a) Result from the new method. (b) Filtered noise by the new method. The SNR of the filtered data in (a) is 15.39 dB. (c) Result from the EMD method. (d) Filtered noise by the EMD method. The SNR of the filtered data in (d) is 2.33 dB.

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

15

Channel NO

Channel NO

15

10

10

5

5

0

50

100

150

200

250

300

0

350

50

100

Time (ms)

200

250

300

350

250

300

350

15

Channel NO

Channel NO

150

Time (ms)

15

10

5

0

5

10

5

50

100

150

200

250

300

350

0

Time (ms)

50

100

150

200

Time (ms)

Fig. 7. Synthetic example with two very close events. (a) Result from the new method. (b) Filtered noise by the new method. The SNR of the filtered data in (a) is 15.39 dB. (c) Result from the EMD method. (d) Filtered noise by the EMD method. The SNR of the filtered data in (d) is 2.33 dB.

Channel

8

6

4

2

500

1000

1500 2000 Time (ms)

Fig. 8. Real data example.

2500

3000

6

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

reference trace. m denotes the index. d(m + τ, n) denotes the microseismic data at coordinate (m + τ, n). When τn N 0, it means we shift the data downward. When τn b 0, it means we shift the data upward. Considering that there might be a multiple of micro-seismic events recorded in the data, the optimal time shift obtained by the crosscorrelation method may not be optimal for arranging all the recorded microseismic waveform events. In this way, an local temporal window is created along the time axis to segment the whole seismic data into different time periods and by this way we can maximize the spatial correlation of the microseismic signals for each event.

  P ¼ ½p1 ; p2 ; ⋯; pr ; V ¼ diag v1 ; v; ⋯; vr ; Q ¼ ½q1 ; q2 ; ⋯; qr :

r X

vk pk qTk ;

ð5Þ

k¼1

When the multi-channel microseismic are shifted to obtain its maximum correlation among all traces, we apply the lowrank approximation method to further enhance the energy of the waveforms. We can use the matrix S to represent the 2D multi-channel microseismic data. The size of the matrix is M × N. M is the number of time samples and N is the number of channels in the microseismic data. We apply the SVD to the data as:

S ¼ PVQ T :

ð3Þ

Here, P is the left singular matrix. It can be computed by the eigenvectors of SST. Q is the right singular matrix and can be computed by STS. V is the singular value matrix. The singular values are stored in the singular value matrix with a decreasing order. The singular values are calculated via the eigenvalues of the covariance matrix SST. We can compute the positive square roots of the eigenvalues to obtain the

where pkqTk is the kth rank-one approximation matrix of the S. We can remove the random noise in seismic data in order to enhance the seismic reflections by only selecting the first several eigenimages (Xue et al., 2016): The random noise in the microseismic data can be easily filtered by selecting K eigenimages from the multi-channel data (Xue et al., 2016): ^¼ S

vk pk qTk :

ð6Þ

The SVD method is very powerful in suppressing random noise. As long as the microseismic waveforms are aligned well spatially, the SVD method mentioned above can help extract the principal components from the data matrix easily. The SVD method, however, is depending on the spatial coherence of the microseismic waveforms. If the waveforms are not aligned well, it may be difficult to obtain a satisfactory separation between signal and noise. It is also the reason why

8

6

6

Channel

8

4

K X k¼1

2

4

2 500

1000

1500 2000 Time (ms)

2500

3000

8

8

6

6

Channel

Channel

ð4Þ

r here is called the rank. The rank of S is assumed to be r. pk is sometimes referred to as the propagation vectors. qk is usually referred to as the eigen-wavelets. The singular values vk are sorted such that v1 ≥ v2 ≥ ⋯ ≥ vr. Eq. 3 can be formulated as: S¼

2.2. Noise attenuation using SVD

Channel

singular values. Let us denote P, V, and Q as:

4

2

500

1000

1500 2000 Time (ms)

2500

3000

500

1000

1500 2000 Time (ms)

2500

3000

4

2

500

1000

1500 2000 Time (ms)

2500

3000

Fig. 9. Denoising comparison. (a) Denoised microseismic data by the new method. (b) Denoised microseismic data by the smoothing method. (c) Denoised microseismic data by the SSA method. (d) Denoised microseismic data by the empirical mode decomposition method.

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

we need a cross-correlation and rearrangement process to precondition the microseismic data. From the perspective of band-pass filtering, the SVD filter is similar to a low-pass wavenumber filter, which passes only the low-wavenumber and spatially correlative signal energy and rejects the high-wavenumber random noise. Next, we will test this method on several examples. 2.3. Fast SVD calculation The SVD method introduced above, however, is computationally expensive. The main computational cost comes from the calculation of eigenvalues and corresponding eigen vectors from the covariance matrix SST. Note that the matrix SST differs from the matrix STS in that the size of SST is much smaller than STS. Suppose the size of S is M × N and M = 300N, the size of SST is 300N × 300N but the size of STS is N × N. Considering that in the microseismic datasets, the number of channels is usually much smaller than the number of time samples. We propose a fast way to calculate the SVD. Since Q is orthonormal, thus QTQ = I. Then SQ ¼ PVQ T Q ¼ PV:

ð7Þ

We can also derive that  T ST S ¼ PVQ T PVQ T ;

ð8Þ

Let σi be the ith eigenvalue of matrix STS, then from eq. 9, we know that σ i ¼ v2i :

ð10Þ

Eq. (10) is important in reduce the computational cost of the traditional SVD since we can use this relation to calculate vi instead of computing it directly from matrix SST. The ith eigenvector can then be obtained by SQ i ¼ vi Pi :

3. Examples One synthetic and one real microseismic datasets will be used in this part to demonstrate the validity of the new denoising method. To better measure the performance, we use the signal-to-noise ratio (SNR) that is widely used in the literature (Chen et al., 2014; Liu et al., 2016b; Zu et al., 2016): SNR ¼ 20 log10

T

S S ¼ Q ΣQ ;

ð11Þ

All the computation of SST are thus saved and we obtain a large boost of the computational efficiency. Note that the computational saving is only meaningful when the shape of the matrix S is ‘narrow and tall’, which is exactly a property of our microseismic data. A comparison of computing efficiency is shown in Table 1.

which can be simplified as T

7

ð9Þ

∥X∥ F ; ∥X−Y∥ F

ð12Þ

where X denotes the exact solution in the denoising problem and Y denotes the microseismic data after filtering.

8

8

6

6

Channel

Channel

where Σ ¼ V2 .

4

2

2 500

1000

1500 2000 Time (ms)

2500

3000

8

8

6

6

Channel

Channel

4

4

500

1000

1500 2000 Time (ms)

2500

3000

500

1000

1500 2000 Time (ms)

2500

3000

4

2

2 500

1000

1500 2000 Time (ms)

2500

3000

Fig. 10. Denoising comparison. (a) Filtered noise by the new method. (b) Filtered noise by the smoothing method. (c) Filtered noise by the SSA method. (d) Filtered noise by the empirical mode decomposition method.

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

1 2 NO

Fig. 1 presents the first example, where 1a plots the clean data that is generated from the finite-difference modeling and 1b plots the noisy data that is constructed by putting some random noise into the clean data. The SNR using the definition in eq. 12 is −4.58 dB. From the SNR value, it is intuitive that the noisy energy is much stronger than the signal. Due to the very strong noise, most useful waveforms are nearly buried in the noise. The noise level in this synthetic example is very close to the noise level in a real microseismic data. The First, we show the rearranged microseismic section using the crosscorrelation method in Fig. 2a. It can be seen that after the rearrangement the microseismic waveforms are aligned well along the axis of “Channel NO”. In Fig. 2b shows the time-shifts calculated from the cross-correlation. The red line plots the time shifts (in point) of each microseismic trace with respect to the first trace in order to make the microseismic data well aligned. The time shifts increases from 0 point to about 28 points and then drops to 25 points for the 17th microseismic trace. We then apply the singular value decomposition (SVD) operation onto the well arranged microseismic gather and then use the time shifts to rearrange back to the original shape. The filtered microseismic data is shown in Fig. 3a, where the microseismic waveforms become very clean. Most random noise has been removed using the SVD filtering. Fig. 3c displays the removed strong random noise. It is difficult to find any waveform signals in the removed noise, indicating the fact that the proposed method is very successful in rejecting the noise without damaging the useful waveform signals. As a comparison, we also test the performance of an empirical mode decomposition method on the same dataset (Han and van der Baan, 2015). In the empirical mode decomposition method, we first decompose each trace into several intrinsic mode functions (IMF), and then we remove the most oscillating IMF which corresponds to the high-frequency noise to output a denoised signal. The result from the empirical mode decomposition method is shown in Fig. 3b and its removed noise is shown in Fig. 3d. It is obvious that the empirical mode decomposition method is not effective in removing the noise. There is still a lot of residual nosie in the denoised microseismic data and the useful waveforms are distorted very seriously. The removeremoved noise is of low-amplitude and indicates a insufficient removal of the noise component. The calculated SNRs of the proposed and the empirical mode decomposition methods are 12.23 dB and −0.95 dB, respectively. The SNR of the proposed method turns out to be much higher than the empirical mode decomposition method. To consider the influences from the erratic incoherent noise which are the real crucial noise to the denoising performance, we conduct another test based on this example to demonstrate the performance in the case of the such incoherent and erratic noise. The test is presented in Fig. 4, where there are two channels corrupted by the very strong erratic noise. In this case, the proposed method still obtain an appealing result.

3 4

500

1000

1500 2000 Time (ms)

2500

3000

Fig. 12. One-trace comparison (the 3rd trace in Fig. 9a) for the real data example. NO 1 denotes the filtered noise by the new method. NO 2 shows the filtered noise by the smoothing method. NO 3 shows the filtered noise by the SSA method. NO4 shows the filtered data by the empirical mode decomposition method.

In the first synthetic example, there is only one microseismic event. It is a little more complicated when there are a multiple of microseismic events. For example, a synthetic microseismic data with two events is shown in Fig. 5a. The two events correspond to microseismic source ignitions. By adding some noise, we have the noisy data as shown in Fig. 5b. We compute the SNR of the noisy data with two events as −1.29 dB. In this case, where several microseismic events are in the same gather, we need to apply a windowing strategy. We segment the whole microseismic gather into small windows, each window with 200 ms in length. In each segmented window, we apply the crosscorrelation and rearragement, and then the SVD filtering. The resulted denoised data using the proposed windowing strategy is shown in

8

Channel

8

6

4

2

500

1000

1500 2000 Time (ms)

2500

3000

8

Channel

1

NO

2 3

6

4

4

2

5

500 500

1000

1500 2000 Time (ms)

2500

3000

Fig. 11. One-trace comparison (the 3rd trace in Fig. 9a) for the real data example. NO 1 shows the noisy data. NO 2 shows the filtered data by the new method. NO 3 shows the filtered data by the smoothing method. NO 4 shows the filtered data by the SSA method. NO5 shows the filtered data by the empirical mode decomposition method.

1000

1500 2000 Time (ms)

2500

3000

Fig. 13. Arrival picking comparison. (a) Automatically picked arrival using the STA/LTA ratio method from the raw data. (b) Automatically picked arrival using the STA/LTA ratio method from the denoised data.

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

Next, we apply the new approach to a low-quality microseismic data that is obtained from a microseismic monitoring project. The real microseismic data is plotted in Fig. 8. There are 9 channels in this microseismic dataset. There are 3000 time points in the dataset. The random noise in this example is extremely strong and the useful signals cannot be possibly detected without being arranged in a multichannel gather. Most of time, a microseismic dataset with a low SNR like this is almost useless. We try to apply the proposed method onto this dataset and to show if the new method can be robust even in the case of extremely strong noise. The denoised data is shown in Fig. 9a, where most random noise has been rejected and the microseismic arrivals become very clear. Fig. 9b shows the denoised data using the smooth filter. Although most random nosie has also been removed using the smoothing filter, most signals are no longer existing. Fig. 9c shows the denoised data using the singular spectrum analysis (SSA) filter. It seems that the SSA filter fails in removing enough noise and the noise level in the denoised data is still very strong. Fig. 9d shows the denoised data using the empirical mode decomposition method, where the useful waveform signals become much clearer but a lot of random noise is still remaining. We plot the removed noise gathers corresponding to the four methods in Fig. 10. From Fig. 10, it is salient that all methods except for the smoothing filter does not cause obvious damages to the useful signals. However, the

0

0

0.5

0.5

1

1

1.5

1.5

Time (s)

Time (s)

Fig. 6a. Correspondingly, Fig. 6b displays the noise that is rejected using the proposed method. In this example, the SNR has increased to 15.39 dB (from originally −1.29 dB). The results from the EMD method are shown in Fig. 6c and d. In the presence of two very close microseismic events, we not only need to segment the data in time but also in space. The space window makes the close events easier to separate in large offset. Fig. 7 shows an example of such case. Fig. 7a shows the clean microseismic dataset containing two close events. Fig. 7b shows the noisy data with strong random noise. The two microseismic events are too close to be separated solely by a time window. In this case, we use a temporal window size of 100 samples and space window size of 5 traces. The resulted denoised data based on the proposed fast algorithm is plotted in Fig. 7c. Although it contains more residual noise as compared with the last example due to the weaker constraint from the spatial coherency, it is still much cleaner than the noisy data. Moreover, when observing the removed noise plotted in Fig. 7d, we found that the proposed method does not damage obvious signal energy. These two examples demonstrate that the proposed method can also be effective when a multiple of microseismic events are existing in the same dataset. However, it is indeed more challenging when the two events are close to each other. Although there are some leaked signals in the noise section, they are very week compared with other signal components.

2

2

2.5

2.5

3

3

3.5

3.5

4

4 20

40

60

80

100

20

Trace

40

60

80

100

80

100

Trace

0

0

0.5

0.5

1

1

1.5

1.5

Time (s)

Time (s)

9

2

2

2.5

2.5

3

3

3.5

3.5

4

4 20

40

60

Trace

80

100

20

40

60

Trace

Fig. 14. Active-source seismic data example. (a) The raw seismic data. (b) Filtered data using the f − x prediction based method. (c) Filtered data using the SSA method. (d) Filtered data using the proposed method.

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

1.6

1.6

1.8

1.8

Time (s)

Time (s)

10

2

2

2.2

2.2

2.4

2.4

40

50

60

70

80

40

50

1.6

1.6

1.8

1.8

2

2.2

2.4

2.4 50

60

70

80

70

80

2

2.2

40

60

Trace

Time (s)

Time (s)

Trace

70

80

40

Trace

50

60

Trace

Fig. 15. Zoomed comparison of the active-source seismic data example. (a) The raw seismic data. (b) Filtered data using the f − x prediction based method. (c) Filtered data using the SSA method. (d) Filtered data using the proposed method.

proposed removemethod removes the most noise thanks to the denoising power of the SVD filtering. The SSA filter removes the least amount of noise, which is followed by the empirical mode decomposition method. For a better view in comparison, we plot a trace-by-trace comparison in Fig. 11. The traces in Fig. 11 are extracted from the 3rd trace in Figs. 8 and 10. The first trace in Fig. 11 corresponds to the noisy data. The second to fifth traces correspond to the proposed method, smoothing filtering method, SSA filtering method, and empirical mode decomposition method, respectively. From this comparison, we further confirm that the proposed method obtain a very satisfactory performance. The removed noise are are plotted in Fig. 12 in a trace-by-trace form. All the noise traces are plotted in the same amplitude scale, thus, we conclude again that the SSA filter fails in removing enough noise while the empirical mode decomposition method removes the secondly least amount of noise. The smoothing filtering method removes the similar level of noise as the new approach. However, the damages of useful signals are more serious for the smoothing filtering method, e.g., damages of the arrival around 1250 ms. Note that, although it is possible to observe the arrivals when watching from the multi-channel gather, but it is extremely difficult to pick the arrival from the single-trace dataset, as seen in Fig. 11 (the second trace). Even though we can observe the arrivals from the multi-channel gather, it is not possible for us to manually pick the arrivals by eyes when thousands of microseismic traces are to be processed. To demonstrate the advantage of the proposed method in automatic arrival picking. We apply the STA/LTA ratio method to the raw microseismic data and denoised data, respectively. The comparison is plotted in Fig. 13. The red circles in Fig. 13 denote the automatically picked arrivals. It is obvious that the raw data is extremely noisy such that the picked arrivals using STA/LTA method are totally incorrect. However, the picked arrivals using the same picking algorithm from the denoised data are very accurate.

It is worth mentioning that the proposed method can also be combined with active-source seismic data to suppress the random background noise. Fig. 14a shows a seismic profile from an active seismic dataset. The active-source data is borrowed from Fomel and Backus (2003). This dataset is a multi-component seismic data. The data presented here denotes the PP component. There are 1024 time points and 101 spatial traces are selected for demonstrating the denoising performance. Fig. 14b, c, and d correspond to the denoised data using the f − x predictive filtering method, the SSA method, and the proposed method, respectively. To make the comparison more apparent, we zoom parts from the these results and show them in Fig. 15. It is quite obvious that although all three methods obtain improvement, the proposed method obtain the cleanest result. It is very clear that the proposed method obtains the best performance in removing the strong background noise. 4. Conclusions Raw microseismic data can be extremely noisy in oil field environment. Because of the strong random ambient noise, the effective signals triggered by the hydraulic perforation are difficult to detect, which will further affect the microseismic event localization. We have presented a novel denoising method for improving the quality of the microseismic data based on cross-correlation and singular value decomposition (SVD). The cross-correlation between each trace and the reference is conducted to estimated the relative shifts between neighbor microseismic channels. The time shifts calculated from the cross-correlation step can be used to optimally rearrange the microseismic data so that the spatial correlation among the channels reaches the maximum. The SVD is applied in the rearranged microseismic data to extract the principal components, i.e., the effective microseismic signals, and thus suppress the random noise. The new method is validated by one synthetic

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

and, one real microseismic, and one real active-source datasets. Results show that the performance of the proposed method is superior to some state-of-the-art methods. Acknowledgements The research is supported by Natural Science Foundation of Jiangxi Province (Grant No. 20181BBG78008). References Baig, A., Urbancic, T., 2010. Microseismic moment tensors: a path to understanding frac growth. Lead. Edge 29, 320–324. Bolton, H., Masters, G., 2001. Travel times of p and s from the global digital seismic networks: implications for the relative variation of p and s velocity in the mantle. J. Geophys. Res. 106, 13527–13540. Capilla, C., 2006. Application of the haar wavelet transform to detect microseismic signal arrivals. J. Appl. Geophys. 59, 36–46. Chen, Y., 2016. Dip-separated structural filtering using seislet thresholding and adaptive empirical mode decomposition based dip filter. Geophys. J. Int. 206 (1), 457–469 no. Chen, Y., 2018. Non-stationary least-squares complex decomposition for microseismic noise attenuation. Geophys. J. Int. 213 (3), 1572–1585. Chen, Y., Fomel, S., 2015. Random noise attenuation using local signal-and-noise orthogonalization. Geophysics 80 (6), WD1–WD9. Chen, Y., Ma, J., 2014. Random noise attenuation by f-x empirical mode decomposition predictive filtering. Geophysics 79, V81–V91. Chen, W., Song, H., 2018. Automatic noise attenuation based on clustering and empirical wavelet transform. J. Appl. Geophys. 159, 649–665. Chen, Y., Fomel, S., Hu, J., 2014. Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization. Geophysics 79 (5), V179–V189. Chen, W., Xie, J., Zu, S., Gan, S., Chen, Y., 2017. Multiple reflections noise attenuation using adaptive randomized-order empirical mode decomposition. IEEE Geosci. Remote Sens. Lett. 14 (1), 18–22. Chouet, B., Dawson, P., Ohminato, T., Martini, M., Saccorotti, G., Giudicepietro, F., De Luca, G., Milana, G., Scarpa, R., 2003. Source mechanisms of explosions at stromboli volcano, Italy, determined from moment-tensor inversions of very-long-period data. J. Geophys. Res. Solid Earth 108. Dong, L., Li, X., 2013. A microseismic/acoustic emission source location method using arrival times of ps waves for unknown velocity system. Int. J. Distrib. Sensor Netw. 9, 307489. Dong, L., Li, X., Xie, G., 2014. Nonlinear Methodologies for Identifying Seismic Event and Nuclear Explosion Using Random Forest, Support Vector Machine, and Naive Bayes Classification: Presented at the Abstract and Applied Analysis. Hindawi Publishing Corporation. Dong, L., Wesseloo, J., Potvin, Y., Li, X., 2016a. Discrimination of mine seismic events and blasts using the fisher classifier, naive bayesian classifier and logistic regression. Rock Mech. Rock. Eng. 49, 183–211. Dong, L.-J., Wesseloo, J., Potvin, Y., Li, X.-B., 2016b. Discriminant models of blasts and seismic events in mine seismology. Int. J. Rock Mech. Min. Sci. 282–291. Dragomiretskiy, K., Zosso, D., 2014. Variational mode decomposition. IEEE Trans. Signal Process. 62, 531–544. Fehler, M., House, L., Phillips, W.S., Potter, R., 1998. A method to allow temporal variation of velocity in travel-time tomography using microearthquakes induced during hydraulic fracturing. Tectonophysics 289, 189–201. Fischer, T., Hainzl, S., Eisner, L., Shapiro, S., Calvez, J. Le, 2008. Microseismic signatures of hydraulic fracture growth in sediment formations: Observations and modeling. J. Geophys. Res. 113. Fisher, M., Heinze, J., Harris, C., Davidson, B., Wright, C., Dunn, K., et al., 2004. Optimizing Horizontal Completion Techniques in the Barnett Shale Using Microseismic Fracture Mapping: Presented at the SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Flewelling, S.A., Tymchak, M.P., Warpinski, N., 2013. Hydraulic fracture height limits and fault interactions in tight oil and gas formations. Geophys. Res. Lett. 40, 3602–3606. Fomel, S., Backus, M.M., 2003. Multicomponent seismic data registration by least squares. 73th Annual International Meeting. vol. 22, pp. 781–784 SEG, Expanded Abstracts. Friedrich, A., Krueger, F., Klinge, K., 1998. Ocean-generated microseismic noise located with the gräfenberg array. J. Seismol. 2, 47–64. Gan, S., Wang, S., Chen, Y., Chen, X., Huang, W., Chen, H., 2016a. Compressive sensing for seismic data reconstruction via fast projection onto convex sets based on seislet transform. J. Appl. Geophys. 130, 194–208. Gan, S., Wang, S., Chen, Y., Chen, X., Xiang, K., 2016b. Separation of simultaneous sources using a structural-oriented median filter in the flattened dimension. Comput. Geosci. 86, 46–54. Gan, S., Wang, S., Chen, Y., Qu, S., Zu, S., 2016c. Velocity analysis of simultaneous-source data using high-resolution semblance-coping with the strong noise. Geophys. J. Int. 204, 768–779. Grigoli, F., Cesca, S., Krieger, L., Kriegerowski, M., Gammaldi, S., Horalek, J., Priolo, E., Dahm, T., 2016. Automated microseismic event location using master-event waveform stacking. Sci. Rep. 6, 1–13. Hajati, T., Langenbruch, C., Shapiro, S., 2015. A statistical model for seismic hazard assessment of hydraulic-fracturing-induced seismicity. Geophys. Res. Lett. 42. Han, J., van der Baan, M., 2015. Microseismic and seismic denoising via ensemble empirical mode decomposition and adaptive thresholding. Geophysics 80 (6), KS69–KS80.

11

Hansen, S.M., Schmandt, B., 2015. Automated detection and location of microseismicity at Mount St. Helens with a large-n geophone array. Geophys. Res. Lett. 42, 7390–7397. House, L., 1987. Locating microearthquakes induced by hydraulic fracturing in crystalline rock. Geophys. Res. Lett. 14, 919–921. Huang, N.E., Shen, Z., Long, S.R., Wu, M.C., Shih, H.H., Zheng, Q., Yen, N.-C., Tung, C.C., Liu, H.H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. Ser. A 454, 903–995. Huang, W., Wang, R., Li, H., Chen, Y., 2017a. Unveiling the signals from extremely noisy microseismic data for high-resolution hydraulic fracturing monitoring. Sci. Rep. 7, 11996. Huang, W., Wang, R., Zu, S., Chen, Y., 2017b. Low-frequency noise attenuation in seismic and microseismic data using mathematical morphological filtering. Geophys. J. Int. 211, 1318–1340. Huang, W., Wang, R., Chen, Y., 2018. Regularized non-stationary morphological reconstruction algorithm for weak signal detection in micro-seismic monitoring: Methodology. Geophys. J. Int. 213 (2), 1189–1211. Li, H., Wang, R., Cao, S., Chen, Y., Huang, W., 2016a. A method for low-frequency noise suppression based on mathematical morphology in microseismic monitoring. Geophysics 81 (3), V159–V167. Li, H., Wang, R., Cao, S., Chen, Y., Tian, N., Chen, X., 2016b. Weak signal detection using multiscale morphology in microseismic monitoring. J. Appl. Geophys. 133, 39–49. Liu, W., Cao, S., Chen, Y., 2016a. Applications of variational mode decomposition in seismic time-frequency analysis. Geophysics 81 (5), V365–V378. Liu, W., Cao, S., Chen, Y., Zu, S., 2016b. An effective approach to attenuate random noise based on compressive sensing and curvelet transform. J. Geophys. Eng. 13 (2), 135–145. Liu, W., Cao, S., Wang, Z., Kong, X., Chen, Y., 2017. Spectral decomposition for hydrocarbon detection based on vmd and teager-kaiser energy. IEEE Geosci. Remote Sens. Lett. 14 (4), 539–543. Liu, W., Cao, S., Jin, Z., Wang, Z., Chen, Y., 2018. A novel hydrocarbon detection approach via high-resolution frequency-dependent avo inversion based on variational mode decomposition. IEEE Trans. Geosci. Remote Sens. 56 (4), 2007–2024. Maxwell, S.C., Urbancic, T., Steinsberger, N., Zinno, R., et al., 2002. Microseismic Imaging of Hydraulic Fracture Complexity in the Barnett Shale: Presented at the SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Maxwell, S.C., Rutledge, J., Jones, R., Fehler, M., 2010. Petroleum reservoir characterization using downhole microseismic monitoring. Geophysics 75, 75A129–75A137. Mayerhofer, M.J., Lolon, E.P., Youngblood, J.E., Heinze, J.R., et al., 2006. Integration of Microseismic-Fracture-Mapping Results with Numerical Fracture Network Production Modeling in the Barnett Shale: Presented at the SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Mousavi, S.M., Langston, C.A., 2016. Hybrid seismic denoising using higher-order statistics and improved wavelet block thresholding. Bull. Seismol. Soc. Am. 106 (4), 1380–1393. Pearson, C., 1981. The relationship between microseismicity and high pore pressures during hydraulic stimulation experiments in low permeability granitic rocks. J. Geophys. Res. 86, 7855–7864. Phillips, W., Fairbanks, T., Rutledge, J., Anderson, D., 1998. Induced microearthquake patterns and oil-producing fracture systems in the Austin chalk. Tectonophysics 289, 153–169. Reshetnikov, A., Buske, S., Shapiro, S., 2010. Seismic imaging using microseismic events: results from the San Andreas fault system at safod. J. Geophys. Res. 115. Reyes-Montes, J., Rietbrock, A., Collins, D., Young, R., 2005. Relative location of excavation induced microseismicity at the underground research laboratory (aecl, Canada) using surveyed reference events. Geophys. Res. Lett. 32. Rigo, A., Lyon-Caen, H., Armijo, R., Deschamps, A., Hatzfeld, D., Makropoulos, K., Papadimitriou, P., Kassaras, I., 1996. A microseismic study in the western part of the gulf of Corinth (Greece): implications for large-scale normal faulting mechanisms. Geophys. J. Int. 126, 663–688. Rutledge, J.T., Phillips, W.S., Schuessler, B.K., 1998. Reservoir characterization using oilproduction-induced microseismicity, Clinton county, Kentucky. Tectonophysics 289, 129–152. Rutledge, J.T., Phillips, W.S., Mayerhofer, M., 2004. Faulting induced by forced fluid injection and fluid flow forced by faulting: an interpretation of hydraulic-fracture microseismicity, Carthage cotton valley gas field, Texas. Bull. Seismol. Soc. Am. 94, 1817–1830. Shapiro, S., Dinske, C., Rothert, E., 2006. Hydraulic-fracturing controlled dynamics of microseismic clouds. Geophys. Res. Lett. 33. Siahsar, M.A.N., Abolghasemi, V., Chen, Y., 2017. Simultaneous denoising and interpolation of 2D seismic data using data-driven non-negative dictionary learning. Signal Process. 141, 309–321. Šílený, J., Hill, D.P., Eisner, L., Cornet, F.H., 2009. Non–double-couple mechanisms of microearthquakes induced by hydraulic fracturing. J. Geophys. Res. 114. Warpinski, N., et al., 2009. Microseismic monitoring: inside and out. J. Pet. Technol. 61, 80–85. Wu, Z., Huang, N.E., 2009. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Adv. Adapt. Data Anal. 1 (1), 1–41. Xue, Y., Chang, F., Zhang, D., Chen, Y., 2016. Simultaneous sources separation via an iterative rank-increasing method. IEEE Geosci. Remote Sens. Lett. 13 (12), 1915–1919. Xue, Y., Man, M., Zu, S., Chang, F., Chen, Y., 2017. Amplitude-preserving iterative deblending of simultaneous source seismic data using high-order radon transform. J. Appl. Geophys. 139, 79–90. Zhang, P., Dai, Y., Wang, R., Tan, Y., 2017. A quantitative evaluation method based on emd for determining the accuracy of time-varying seismic wavelet extraction. J. Seism. Explor. 26, 267–292.

12

H. Lv / Journal of Applied Geophysics 170 (2019) 103831

Zhou, Y., Han, C., 2018. Seismic data restoration based on the grassmannian rank-one update subspace estimation method. J. Appl. Geophys. 159, 731–741. Zhou, Y., Han, C., Chi, Y., 2018. Deblending of simultaneous-source data using iterative seislet frame thresholding based on a robust slope estimation. J. Appl. Geophys. 138, 17–37.

Zu, S., Zhou, H., Chen, Y., Qu, S., Zou, X., Chen, H., Liu, R., 2016. A periodically varying code for improving deblending of simultaneous sources in marine acquisition. Geophysics 81 (3), V213–V225. Zu, S., Zhou, H., Ru, R., Jiang, M., Chen, Y., 2019. Dictionary learning based on dip patch selection training for random noise attenuation. Geophysics 84 (3), V169–V183 no.