Iterative interferometry-based method for picking microseismic events

TL;DR

This paper introduces an interferometry-based iterative algorithm for accurately picking microseismic first breaks in noisy data, improving signal detection with minimal iterations.

Contribution

The study proposes a novel seismic interferometry and time-frequency analysis method for microseismic first-break picking, reducing noise influence and requiring only two iterations for convergence.

Findings

Significantly improves first-break detection in noisy microseismic data

Converges within two iterations, avoiding noise correlation issues

Effective on both synthetic and real datasets

Abstract

Continuous microseismic monitoring of hydraulic fracturing is commonly used in many engineering, environmental, mining, and petroleum applications. Microseismic signals recorded at the surface, suffer from excessive noise that complicates first-break picking and subsequent data processing and analysis. This study presents a new first-break picking algorithm that employs concepts from seismic interferometry and time-frequency (TF) analysis. The algorithm first uses a TF plot to manually pick a reference first-break and then iterates the steps of cross-correlation, alignment, and stacking to enhance the signal-to-noise ratio of the relative first breaks. The reference first-break is subsequently used to calculate final first breaks from the relative ones. Testing on synthetic and real data sets at high levels of additive noise shows that the algorithm enhances the first-break picking considerably. Furthermore, results show that only two iterations are needed to converge to the true first breaks. Indeed, iterating more can have detrimental effects on the algorithm due to increasing correlation of random noise.

Equations

$$Y(p,q)=\sum _{j=0}^{L_w-1}{y}_{j+qh}{w}_j\exp \left(\frac{-i2\pi jp}{L_w}\right),p=0,1,\cdots ,L_w-1,q=0,1,2,\cdots ,\frac{N-L_w}{h}$$

$${\Phi }_{lm}[\tau ]=x_l\otimes x_m=\sum _n{x}_l[n]x_m[n+\tau ],$$

$${\tau }_{rm}^i=\{\begin{array}{ll}\underset{\tau }{\mathrm{argmax}}{\Phi }_{rm}^i[\tau ]& r<m\\ 0& r=m\\ -\underset{\tau }{\mathrm{argmax}}{\Phi }_{rm}^i[\tau ]& r>m\end{array}$$

$$\text{ISSE}(i)=\sum _{m=1}^M{\left(n_m^i-n_m^{i-1}\right)}^2=\sum _{m=1}^M{\left({\tau }_{rm}^i-{\tau }_{rm}^{i-1}\right)}^2$$

$${\Phi }_s^i[\tau ]=\frac{1}{Q}\sum _{l<m}{\Phi }_{lm}^i[\tau -{\tau }_{lm}^i],\text{mbox}where{\tau }_{lm}^i=\underset{\tau }{\mathrm{argmax}}{\Phi }_{lm}^i[\tau ]$$

$${\hat{\Phi }}_{lm}^{i+1}[\tau ]$$

$${\Phi }_{lm}^{i+1}[\tau ]$$

$$T[\tau]=\left\{\begin{array}[]{ll}1&\mbox{ for }-N_{T}\leq\tau\leq N_{T}\\ 0&\mbox{ otherwise}\end{array}\right.$$

$$x_m[n]=x_m[n]p_m[n],m=1,\dots ,M$$

$$p_{m}[n]=\left\{\begin{array}[]{ll}1&\mbox{ for }1\leq n\leq S_{m,arrival}\\ 0&\mbox{ otherwise}\end{array}\right.$$

$$W(t)=A{t}^n{e}^{-\alpha t}\cos (2\pi ft+\varphi ),$$

$$\text{SNR}=10{\log }_{10}\left(\frac{{\sigma }_s^2}{{\sigma }_n^2}\right),$$

$$\text{PSNR}=10{\log }_{10}\left(\frac{P^2}{{\sigma }^2}\right),$$

$$\text{ISSE}(i)=\sum _{m=1}^M{\left(n_m^i-n_m^{i-1}\right)}^2$$

$$\text{TSSE}(i)=\sum _{m=1}^M{\left(n_m^i-n_{\text{true}}\right)}^2$$

Full text

Key words— cross correlation, first-break, time-frequency representation, P- and S-arrivals.

1 Introduction

Renewed interest in microseismic monitoring of hydraulic fracturing for the development of unconventional reservoirs has led the petroleum industry to invest in microsesismic technology. The technology involves installing downhole and surface sensors to detect and localize microseismic events in order to image fracturing associated with fluid injection or extraction. The microseismic fracture image is used to understand the extent of fracture network. However, microseismic events are known to have very low signal-to-noise ratio (SNR) making it difficult to detect and accurately locate microseismic events. Many researchers have proposed seismic interferometry to increase the SNR of seismic data (Snieder,, 2004; Wapenaar et al.,, 2008; Schuster,, 2009; Wapenaar et al., 2010a, ; Wapenaar et al., 2010b, ). Recent application of seismic interferometry in first-break enhancement through cross-correlation, summation, and convolution showed promising results when applied on seismic signals generated by an active source, e.g., (Mallinson et al.,, 2011; Alshuhail et al.,, 2012; Bharadwaj et al.,, 2013; Al-Hagan et al.,, 2014).

Xiao et al., (2009) used seismic interferometry to locate microseismic events by cross-correlating the direct P and S arrivals from repeated sources. They tested their approach using an elastic model and found that the repeated sources enhance arrivals after stacking. Song et al., (2010) proposed an array-based waveform correlation method to enhance the detectability of microseismic events. They used a transformed spectrogram to identify the arrivals and found an improvement over an array-stacked short-time average/long-time average (STA/LTA) approach. There are various versions of STA/LTA method, which are based on energy, amplitude or entropy functions (Wong et al.,, 2009; Xiantai et al.,, 2011). There are other similar methods, e.g., multi-window method (Chen and Stewart,, 2005) and modified energy ratio (MER) (Wong et al.,, 2009). Mousa et al., (2011) discussed a method based on the digital image segmentation. Al-Shuhail et al., (2013) proposed a workflow to enhance microseismic events and reported that the last step in their workflow (i.e., convolving the enhanced and raw records) seems responsible for leaking noise from the raw records into the enhanced data. This problem is eliminated by Iqbal et al., (2016) using singular-value decomposition in place of convolution.

Among the existing methods, STA/LTA method is the most widely used in earthquake seismology. The STA/LTA method processes the signals in two moving windows (long and short). The ratio of the average energies in the windows are calculated. In case of noise-free seismograms, the maximum value of numerical derivative of the ratio is close to the first-arrival time. In the MER method, the concept of STA/LTA is modified in such a way that two windows are of equal length and move adjacent to each other. The peak of the cubic power of energies ratio in two windows is close to a first-break. The other technique to detect the first breaks includes interferometry approach where cross-correlations of all the distinct pair of traces are carried out. Then, the cross-correlations are aligned to a specific instant and summed up. The resulting stacked cross-correlation is considered as a filter to denoise the noisy traces by convolving the stacked cross-correlation with the traces.

In this study, we propose a first-break picking approach based on seismic interferometry. The proposed approach applies the basic steps of seismic interferometry (cross-correlation and summation) iteratively in order to minimize the sum of squared errors (SSE) between successive iterations. We use both synthetic and real microseismic data to demonstrate that this new approach significantly enhances first breaks, thus making it easy to pick event arrival times. Consequently, the enhanced signals with accurate picking of first breaks will likely improve microseismic event localization over the reservoir.

2 Proposed Method

The proposed method has two phases: reliable/accurate first-break estimation on a reference trace and estimation of the relative time delays from the reference trace to all other traces with an enhanced cross correlation function (CCF). Assume there are $M$ raw microseismic traces, $x_{m}[n]$ for $m=1,\ldots,M$ of length $L_{T}$ and sampling interval $\Delta t$. With $M$ traces, the number of unique CCFs is $Q=(M/2)(M-1)$.

The success of the proposed method depends on the correct picking of a reference first-break on at least one trace; therefore, we have given special attention to this step. Consequently, we use the time-frequency contents of the traces simultaneously in the time-frequency domain, which will make the process of first-break manual picking more reliable. Time-frequency representation is not a new concept and many high resolution time-frequency decompositions have been developed. The short-time Fourier transform (STFT) and the continuous wavelet transform are well-known time-frequency transforms that can recover the signal contents if they do not overlap in the time-frequency domain (Diallo et al.,, 2005; Kulesh et al.,, 2007; Roueff et al.,, 2004). There are other transforms such as the empirical mode decomposition (Huang et al.,, 1998), synchrosqueezing transform (Daubechies et al.,, 2011), matching pursuit (Mallat and Zhang,, 1993) and basis pursuit (Bonar and Sacchi,, 2010; Vera Rodriguez et al.,, 2012; Zhang and Castagna,, 2011), to name a few. The time-frequency representation has been used for travel-time picking in the past (Saragiotis et al.,, 2013; Zhang and Zhang,, 2015; Herrera et al.,, 2015). In this work, we have tested various time-frequency transforms but find spectrogram method to be most effective. Hence, we use the spectrogram method proposed by Auger and Flandrin, (1995) to pick a reference trace $x_{r}[n]$ whose first-break $n_{r}$ is the most clear in the time-frequency domain.

For Spectrogram, Fast Fourier Transform (FFT) is applied on subset of data points $(N)$. This subset of data is selected using a window $\mathbf{w}=[w_{0},w_{1},\cdots,w_{L_{w}}]$. First, FFT is computed for data points of length $L_{w}$, then the window in moved $h$ data points and again FFT is calculated. Here, $h$ is called as step size. This procedure is repeated until the window covers the last $L_{w}$ data points. Mathematically,

[TABLE]

The spectrogram is calculated as$|Y(l,n)|^{2}$. We have used the Hamming window of length $(L_{w})$ quarter of the trace length and the overlap is $L_{w}-1$, i.e., $h=1$ (these parameters are optimized using simulations).

The reference time index $n_{r}$ will be utilized to calculate the absolute timings of the enhanced first breaks after the second step which involves multiple cross correlations. The units of $n_{r}$ are number of samples, which can be converted to time using the sampling interval, i.e., $t_{r}=(n_{r}-1)\Delta t$. Another step in the method is to accurately find the relative time delays of all other traces to the reference trace. Then, the first-break for an arbitrary trace is the sum of $t_{r}$ and the relative delay to the reference trace. See Algorithm 1 for a detailed description of the proposed scheme.

Notation: For two distinct traces, the CCF is defined as

[TABLE]

where $\otimes$ denotes the convolution operator and $\tau$ is the lag index, and $-L_{T}\leq\tau\leq L_{T}$, ($L_{T}$ denotes the length of the traces. i.e., total number of samples).

The algorithm use $\displaystyle\tau_{lm}=\operatornamewithlimits{argmax}_{\tau}\{\Phi_{lm}[\tau]\}$ as the initial picks of the relative delays. These delays enable CCF stacking (step 3 of the Algorithm 1) which is a key step in this proposed scheme, where we iteratively perform cross-correlation of the CCFs to mitigate the negative impact of the high noise on the first-break picking. However, in the sequel we will also see that over iterating this scheme will ultimately harm the results.

It is worth mentioning here that in case of both P- and S-arrivals, the algorithm is applied sequentially to find both arrivals. Since S-waves have relatively high amplitudes than P-waves, application of the algorithm in the first place leads us to estimate S-arrivals. After getting the S-arrivals ($S_{m,arrival}$), the P-arrivals are found by again applying the algorithm on the data set with S-waves removed. This is done by taking the part of data up to the $S_{m,arrival}$, i.e,

[TABLE]

where

[TABLE]

Accordingly, the reference trace will be used to pick S-arrival reference in the first place and P-arrival reference in the second place.

3 Results

This section outlines the performance of the proposed cross-correlation based method for the cases of extremely noisy synthetic and real data records.

3.1 Test on synthetic microseismic data

First, we apply our work flow on a synthetic microseismic record generated using the following parameters:

• •

Receivers are located in a vertical borehole at $x_{m}=0$, $y_{m}=0$, $z_{m}=-1000-(m-1)\Delta z$ for $m=1,\ldots,M$, where $M=14$ and $\Delta z=50$ m as depicted in Figure 1.

• •

Constant medium velocity of $2000$ m/s.

• •

Source wavelet is a normalized minimum-phase Berlage wavelet (Aldridge,, 1990)

[TABLE]

with $A=1$, $n=0.001$, $\alpha=15$, $f=5$ Hz, and $\phi=-\pi/2$. This wavelet is frequently used to represent a seismic source in the literature (e.g., Poiata et al., (2016)).

• •

Sampling interval is $\Delta t=1$ ms with 1001 samples per trace.

• •

Source coordinates are $x_{s}=500$ m, $y_{s}=250$ m, and $z_{s}=-1800$ m.

Figure 2a shows the traces before adding noise in order to illustrate the correct first breaks. Figure 2b shows the traces after adding noise and trace normalization. The additive white Gaussian noise has zero mean and variance $\sigma_{n}^{2}$, and the noise less signal has variance $\sigma_{s}^{2}$. The signal-to-noise ratio

[TABLE]

is set to $-$12 dB for the synthetic data and the corresponding peak-signal-to-noise ratio

[TABLE]

is 1 dB, where $P$ is the peak value of the known wavelet. When the wavelet duration is short, the PSNR is more intuitive than SNR, because its value is not affected by the signal length. At this PSNR level, automatic first-break picking from the noisy records is challenging. Figure 3 depicts the spectrogram or time-frequency transform (TFT) of trace 1 of the synthetic data. The time-frequency processing uses a Hamming window whose length is $L_{T}/4$, where $L_{T}$ is the total length of a trace. The time-frequency plot of trace 1 in Figure 3 shows the clear case for manually picking the arrival time of the microseismic event due to the confinement of the seismic energy to a dominant region (shown in detail by the rectangle and its zoomed version). We pick the first-break of the microseismic event on trace 1 and save it as $n_{r}=289$ samples. In most cases, manual picking of the arrival time of the microseismic event can be done easily with the TFT of the first trace. For cases with even lower PSNR, when the visibility of the signal is not clear in time-frequency domain, more than one trace can be transformed to the time-frequency domain to make a robust decision.

Next, the iteration is started by performing all pairwise cross-correlations. Raw cross-correlations are not shown here since they do not give any useful information to the reader; instead aligned CCFs are shown in Figure 4a (and to make the picture clear only 15 out of 91 CCFs are shown). The alignment of CCFs is obtained by shifting each cross-correlation such that its maximum value occurs at lag $\tau=0$. In succeeding iterations, modified CCFs are produced in the update of Step 4, so Figure 4b, 4c, and 4d also depict the aligned modified cross-correlations for each iteration. It is clear that the aligned modified CCFs after iteration 1 are less noisy when compared to the aligned CCFs after iteration 0. One notable feature that can also be seen in these figures is that, as the number of iterations increases, the modified CCFs exhibit an oscillating behavior. The reason is that as we iterate, the Fourier transform of the modified CCFs is essentially the Fourier transform of the signal raised to higher powers. As a result the dominant frequency amplitudes will be enhanced multiple times, and finally yields modified cross-correlations that are very close to a cosine function. Therefore, the stacked cross-correlation after a number of iterations does not help to locate an accurate first-break, because of its oscillatory behavior.

After alignment the (modified) cross-correlations are stacked with the result being shown in Figure 5. This figure validates the points highlighted for Figure 4. The spike in the stacked cross-correlation after iteration 0 is due to the noise and misalignment of the CCFs. These problems are reduced significantly and removed automatically after iteration 1 because of increasing signal contribution and decreasing noise contribution with iteration number.

The raw traces after posting the first breaks from each iteration are shown in Figure 6. Since we have to iterate through step $1$ and step $2$ of the proposed algorithm, therefore, we use the iterative sum of squared errors (ISSE) as a stopping criteria. ISSE is calculated as the sum of the difference of first breaks obtained for two consecutive iterations (current and previous)

[TABLE]

We also compute the true sum of squares error (TSSE) as the sum of the differences of first breaks and the true first breaks

[TABLE]

Plots of the ISSE and TSSE versus iteration number are depicted in Figure 7. For plotting ISSE/TSSE, 20 Monte Carlo simulations are performed. The behavior of TSSE confirms that we can base our stopping criteria on the ISSE. It can be seen that TSSE first decreases and then increases after iteration 2 and ISSE follows the same pattern. The increase in TSSE is due to the enhanced dominant frequency amplitudes mentioned earlier for Figures 4 and 5.

The performance of the algorithm for different PSNRs and iterations is shown in Figure 8. A 5Hz minimum-phase Ricker wavelet is used as the source waveform and 500 Monte Carlo simulations were performed. Figure 8 shows that the TSSE increases after first iteration for all PSNRs. Figure 8 also compares the performance of the proposed technique with and without truncation window $T[\tau]$ (represented by solid lines and dotted line, respectively). With truncation, the cross-correlations in step 4 (see Table 1) of the algorithm are taken as zero after $N_{T}$ (or $-N_{T}$ on the negative lag) samples. The value of $N_{T}$ = 350 is used in the simulation, which is quarter of the total signal duration (1000 samples). On the other hand, without truncation cross-correlations are left unchanged. The performance improvement due to truncation window $T[\tau]$ is evident from the figure.

We also apply our proposed algorithm to zero-phase Ricker wavelet with dominant frequency of 30Hz. Figure 9 shows the algorithm converges in first iteration.

In order to show the effectiveness of the proposed method, we also present the result of the STA/LTA method (with zero-phase Ricker wavelet of 30Hz) in Figure 10. It is apparent from the figure that the STA/LTA method is not applicable for low SNR signal, since many false peaks are comparable or even surpass the local peaks on the correct arrivals.

3.2 Test on real data

The algorithm was originally tested on the real data set from the local industry in the Middle East. However, due to proprietary issues, we are unable to include those results. In order to demonstrate the effectiveness of the proposed method, we now apply our work-flow on a field data set (Figure 11a). The data used in this study comes from the High Resolution Seismic Network (HRSN) operated by Berkeley Seismological Laboratory, University of California, Berkeley, and the Northern California Seismic Network (NCSN) operated by the U.S. Geological Survey, Menlo Park, and is distributed by the Northern California Earthquake Data Center (NCEDC). The sampling frequency is 250 Hz. We repeat a single waveform 13 times with time shifts that simulate multiple traces recorded at distances far from the source. The microseismic event originally has a good SNR and did not need first-break enhancement; but for testing we add zero-mean white Gaussian noise. The PSNR is set to 8 dB in order to make the first-break challenging for automatic picking, as shown in Figure 11b.

In Figure 12a, we examine the TFT computed from the real data. It appears that trace 1 gives a good indication of the arrival time of P- and S-wave of the microseismic event in the time-frequency representation as depicted in Figure 12a. Therefore, trace 1 was used as the reference trace to manually pick the arrival time of the P-wave (or S-wave) microseismic event and saved as $n_{r}=355$ (or $=896$) samples. It should be noted here that first, we pick the S-arrival reference and apply the aforementioned algorithm to get the S-arrivals of all the traces. Next, P-arrival reference is picked and algorithm is again applied after removing the S-waves. This gives the P-arrivals of all the traces. In Figure 12b we show cumulative frequency amplitude in each time sample, which serves as a good reference for the time picking in the spectrogram.

Figures 13 and 14 show the CCFs after alignment and stacking, respectively, during the iterative algorithm of our proposed approach, respectively. The figures for CCFs and stacked CCFs are plotted for S-arrivals only, however, for P-arrivals same procedure should be followed. The raw traces after posting the S-arrivals from each iteration are depicted in Figure 16, as well as plots of the ISSE and TSSE versus iteration number in Figure 15 (resulting from 20 Monte Carlo simulations). Next, to find P-arrivals, S-waves are muted according to Eqs. (3) and (4) and the P-waves only data set is shown in figure 17. After application of the algorithm again, P-arrivals are shown in Figure 18, which confirms the algorithm’s performance for P-arrivals as well. It can be seen from these figures that results of testing the proposed algorithm on real data agree well with tests on synthetic.

4 Discussion and conclusion

We have developed a new approach for picking first breaks of noisy microseismic records using an iterative approach based on seismic interferometry. Unlike existing interferometry-based approaches our approach does not involve convolution of the stacked CCF with raw traces, which ensures that the raw data does not mix with the enhanced stacked record. Using synthetic data as well as real data, we demonstrate a significant improvement in first break picking. Consequently, the enhanced first breaks will improve microseismic event localization over the reservoir. Results from synthetic data show that the SSE between the true and estimated first breaks does not change considerably after the first iteration, which indicates the fast convergence rate of the proposed method. Therefore, we recommend iterating only once and using more iterations only if the SSE changes slowly. The reason for this observation is that in the first CCF between the raw traces the signal could be highly variable from trace to trace due to effects of noise and attenuation. In contrast, in the first iteration, we cross-correlate the aligned stacked CCF which essentially consists of an average auto-correlation of the signals with the cross-correlations of trace pairs which consists of an approximate auto-correlation of the signal in this trace pair; hence, they should exhibit more similarity than raw traces. We believe that the approach proposed in this study is important and useful for the petroleum industry as well as academic and research institutions interested in continuous reservoir monitoring using microseismic data. Thus, we recommend using our newly developed workflow as a first-break picking tool during the processing of microseismic data recorded over the reservoir.

5 Acknowledgements

We thank King Fahd University of Petroleum & Minerals (KFUPM) for supporting this study through the Center for Energy and Geo Processing (CeGP) at KFUPM under grant number GTEC1311.