Estimation of the Relative Arrival Time of Microseismic Events Based on Phase-Only Correlation

The arrival time of a microseismic event is an important piece of information for microseismic monitoring. The accuracy and efficiency of arrival time identification is affected by many factors, such as the low signal-to-noise ratio (SNR) of the records, the vast amount of real-time monitoring records, and the abnormal situations of monitoring equipment. In order to eliminate the interference of these factors, we propose a method based on phase-only correlation (POC) to estimate the relative arrival times of microseismic events. The proposed method includes three main steps: (1) The SNR of the records is improved via time-frequency transform, which is used to obtain the time-frequency representation of each trace of a microseismic event. (2) The POC functions of all pairs of time-frequency representations are calculated. The peak value of the POC function indicates the similarity of the traces, and the peak position in the time lag axis indicates the relative arrival times between the traces. (3) Using the peak values as weighting coefficients of the linear equations, consistency processing is used to exclude any abnormal situations and obtain the optimal relative arrival times. We used synthetic data and field data to validate the proposed method. Comparing with Akaike information criterion (AIC) and cross-correlation, the proposed method is more robust at estimating the relative arrival time and excluding the influence of abnormal situations.


Introduction
Microseismic monitoring is widely used to predict rockbursts in coal mining, minimize the risks of CO 2 leakage, characterize geothermal energy reservoir, and evaluate hydraulic fracturing [1][2][3][4][5][6][7][8].The arrival time of microseismic event is an important piece of information of microseismic monitoring for phase identification, source location, source mechanism analysis, and microseismic interpretation [9][10][11][12][13].The arrival time is often identified by a change from background noise [14].Multiple attributes of the record are used to identify the arrival time, including the amplitude, energy, frequency content, and polarization.The ratio of the short-term average over long-time average (STA/LTA) is the most commonly used algorithm for identifying arrival times [15].This method is sensitive to background noise and the length of the window [16].Some modified algorithms have been proposed, such as modified energy ratio, recursive STA/LTA, and characteristic function preprocessing [17][18][19][20][21].The Akaike Information Criterion (AIC) assumes that a time series can be broken down into two pseudo-stationary segments.A pseudo-stationary segment has unchanging mean, variance, and autocorrelation [22,23].The properties of higher-order statistic parameters, such Energies 2018, 11, 2527 2 of 16 as kurtosis or skewness, have also been used to detect the arrival time [24].Many time-frequency approaches for picking arrival times have been proposed, such as continuous wavelet transform, synchrosqueezed continuous wavelet transform, empirical mode decomposition, and Cohen's class distribution [25][26][27][28].Coherence-based methods have been used to estimate the relative arrival time among traces.Artificial neural network, entropy function, and image processing method were used to pick arrivals automatically [29][30][31].With the development of machine learning, some unsupervised and supervised algorithms were used to detect the arrival times of microseismic events [32,33].Since a single method may fail in accurately identifying arrival times, several hybrid methods have been proposed.Zhang et al. designed a method that combines AIC and wavelet transform [34].Saragiotis et al. discussed a hybrid method that combines the advantages of the coherence method and time-frequency analysis [35].Akram and Eaton proposed an STA/LTA-AIC hybrid approach [36].Li et al. used empirical mode decomposition to improve the signal-to-noise ratio of record and applied AIC to pick the arrival times [37].
It is well known that a low signal-to-noise ratio of record remains a challenge for accurate arrival time identification.Meanwhile, real-time microseismic monitoring has been widely applied.
Due to the high sample rate, real-time monitoring may generate a vast number of continuous records.In addition, abnormal situations of monitoring equipment may arise in the field (for example, abnormal geophone coupling, damaged transmission lines, or geophone faults).These factors render arrival time identification more complex.Accuracy and efficiency become two important indicators of the arrival time identification method.Therefore, we propose a new method based on phase-only correlation (POC) to estimate the relative arrival times of downhole microseismic data.
It is hard to know the accurate origin time of microseismic events before source location [38].This means that the travel times of microseismic events cannot be directly converted from the arrival times.Thus, the relative arrival times are of more practical value for monitoring processing.Due to similar source locations and focal mechanisms, a cluster of seismic events show similar waveforms on the record [39].Correspondingly, the traces of one event recorded by a geophone array also reveal waveform similarity.Coherence-based methods can be used to pick the relative arrival times with the information of multi-trace records [40][41][42].The proposed method applies POC to estimate the relative arrival time.Previously, POC was an image-matching method for fingerprint identification [43].POC evaluates the similarity level by calculating the peak value of the cross-phase spectrum of two two-dimensional (2D) images.Later, POC was applied to identify the similarity of seismic events, diagnose complex diseases, and detect microseismic events [44][45][46].The two main properties of the POC are as following: (1) the POC function is not influenced by the intensity difference of data because POC only uses the phase information of data; and (2) the peak position of POC will shift when one image is shifted, and the peak value is invariant [44].So, we can estimate relative arrival time among the traces by detecting the position of the peak value.Meanwhile, the peak values of the POC are used as weighting coefficients to improve the consistency of estimated relative arrival times and exclude the interferences of abnormal situations.
We used synthetic data and a high signal-to-noise ratio (SNR) microseismic event to validate the proposed method.Firstly, we used synthetic data to validate the proposed method.Secondly, we used a high SNR microseismic event to generate three testing datasets by adding different noise levels to the microseismic event and replacing one trace of the event with background noise to simulate abnormal situations.Lastly, AIC, cross-correlation, and POC were used to calculate the relative arrival times of the two testing datasets.Moriya used moving time window spectral estimation (that is, short time Fourier transform (STFT)) for POC [44].Zhang replaced the STFT with Wigner-Ville distribution (WVD) to overcome the low resolution of STFT [47].We used these two time-frequency transforms to implement the POC function and named them POC (STFT) and POC (WVD), respectively.The list of symbols used in this article is given in Table 1.

Symbol Definition
f (n, w) the time-frequency representation of one trace the phase component of the phase-only correlation (POC) function between f (n, w) and g(n, w) r f g the peak value of the POC function r f g (n, ω) ∆t ij the relative arrival time between the ith and jth traces

Time-Frequency Transform
Given the two traces of a microseismic event, R 1 (n) and R 2 (n), suppose that the waveforms of two traces are identical expect for the time shifting, and the arrival time difference of this two traces is Using Wigner-Ville distribution (WVD) [48], the time-frequency representations of the two time series are obtained.Suppose that the time-frequency representations of R 1 and R 2 are f (n, ω) and g(n, ω), respectively.
where n is time, ω is angular frequency, k is time shifting, R * denotes the complex conjugate of R.
The POC function r f g (n, ω) is obtained by using the inverse 2D discrete Fourier transform of RFG (k 1 , k 2 ).
Energies 2018, 11, 2527 4 of 16 where n denotes the time lag and ω denotes the frequency lag.The range of r f g (n, ω) is from 0 to 1.When f (n, ω) and g(n, ω) are similar, the POC function r f g (n, ω) exhibits a distinct sharp peak.The position of the peak value (the maximum value of the POC function) is located at the center point of r f g (n, ω).When f (n, ω) and g(n, ω) are not similar, the peak drops obviously [43].Therefore, the peak value of the POC can be used to indicate the similarity between f (n, ω) and g(n, ω).
According to the arrival time difference in the time domain between the two traces, the time shifting of the event in the time-frequency domain is n delay (that is, g(n, ω) = f (n + n delay , ω)).By using the time-shift property of DFT, the POC function between f (n, ω) and g(n, ω) can be given by: As previously described, the only difference between R 1 (n) and R 2 (n) is the time shifting n delay , and the peak value of the POC function r f f (n, ω) is located in the center point.So, Equation ( 6) implies that the position shifting between the peak value of r f g (n, ω) and the peak value of r f f (n, ω) in the time lag axis can indicate the time shifting between R 1 (n) and R 2 (n).Therefore, the relative arrival times can be estimated by detecting the peak location of the POC.It is obvious that the position shifting of the peak value in the frequency lag axis can indicate the main frequency shifting between R 1 (n) and R 2 (n).To suppress the interferences in the time-frequency domain, we performed a low-pass filtering for the cross-phase spectrum with a 2D Hamming window.

Consistency Processing
Let t i and t j be the arrival time of the ith and jth traces, respectively.Then, the arrival time difference ∆t ij can be given by: Due to noise interference or abnormal situations of monitoring equipment, the relative arrival times are not always consistent (for example, ∆t 12 − ∆t 23 = ∆t 13 ).Meanwhile, in the normal situation, there are high similarities among the traces of a microseismic event record.The peak value of POC is high.If the waveform of one trace is abnormal (this implies the abnormal situations of monitoring equipment), the similarities will be low between this trace and other normal traces.Accordingly, the peak value of the POC will be low.So, in order to improve the consistency of the estimated values of relative arrival time and exclude the influence of abnormal situations, we used the over-determined linear equations to obtain the optimal result [39].
The peak values of the POC functions are added into Equation ( 7) (shown as Equation ( 8)).
where r ij denotes the peak value of the POC between the ith and jth traces and ∆t ij denotes the arrival time difference obtained using POC.
Energies 2018, 11, 2527 5 of 16 Although Equation ( 8) is identical to Equation ( 7) for a single pair of traces, we can obtain the optimized result of arrival times by using all pair of traces to construct linear equations (shown as Equation ( 9)).
where t i denotes the optimized relative arrival time.The mean of the optimized relative arrival times must be equal to zero (that is, t 1 + t 2 + t 3 + t 4 = 0).Then, t 1 ∼ t 4 are actually the relative arrival times, not the arrival times in the corresponding traces.
As shown in Figure 1, the workflow of the proposed method consists of three main steps: (1) calculating the time-frequency representations of traces; (2) using the POC function of all pairs of traces to obtain peak values (or similarities between traces) and peak positions (or arrival time differences among traces); and (3) using weighting linear equations to obtain the optimal relative arrival times.

Synthetic Data Analysis
To validate the proposed method, we synthesized a record with four traces and a sampling of 0.5 ms (Figure 2a).Each trace included an exponential decaying sinusoidal signal, the frequency of which was 300 Hz.The length of a single trace was 150 ms (or 300 sampling points).The amplitude ratio of the four traces was set as 4:3:2:1, and the lag times of adjacent traces were set as 15 ms (or 30 sampling points).Then, we generated noisy data by adding Gaussian white noise to the previous record.The signal-to-noise ratio of the noisy data was 0 dB (Figure 2b).The signal-to-noise ratio is defined as where σ signal denotes the standard deviation of the signal and σ noise denotes the standard deviation of the noise.The first trace (or Trace 1) and the fourth trace (or Trace 4) of the noisy data were selected to illustrate the process of the proposed method.Figure 3 shows the time-frequency representations of Trace 1 and Trace 4 obtained using the Wigner-Ville Distribution.It was difficult to detect the sinusoidal signal of Trace 4 in the time domain.However, we were able to observe the signal in the time-frequency domain due to the time-frequency sparseness of this signal.So, using time-frequency transform can improve the signal-to-noise ratio of a record.This result is consistent with the theoretical value.Meanwhile, it is clear that the peak value of the POC function in Figure 4 is less than 1.This result is attributed to the effect of the low-pass filtering for the cross-phase spectrum.
As shown in Figure 4, the range of the time lag axis is from −75 to 75, which is half the length of a single trace.If the arrival time difference between two traces were bigger than half the length of a single trace, the peak location would not directly indicate the arrival time difference, and we would have to recalculate to obtain the accurate arrival time difference based on the peak location, the periodicity property of Fourier transform, and other previously determined information.The analyzed data should be long enough to avoid this situation.For microseismic monitoring, the length can be estimated by using the prior information of the monitoring survey and the downhole velocity structure.By calculating all of the pairs of traces, we can obtain all arrival time differences among all traces (Table 2).The estimation of the relative arrival time of the noisy data was obtained after the consistency processing.The estimated relative arrival times of the four traces are −22.6 ms, −7.6 ms, 7.5 ms, and 22.8 ms, respectively.Comparing the theoretical values with the estimated values, the relative arrival time of the first trace was forced to zero.Then, the estimated relative arrival times were 0 ms, 15 ms, 30.1 ms, and 45.4 ms, respectively.If the relative arrival time of one trace is positive, it indicates that the arrival time of this trace is lagging behind that of the first trace.Otherwise, the arrival time of this trace is ahead of that of the first trace.As shown in Table 3, the estimated values of relative arrival times are highly consistent with the theoretical values.The maximum error is 0.4 By calculating all of the pairs of traces, we can obtain all arrival time differences among all traces (Table 2).The estimation of the relative arrival time of the noisy data was obtained after the consistency processing.The estimated relative arrival times of the four traces are −22.6 ms, −7.6 ms, 7.5 ms, and 22.8 ms, respectively.Comparing the theoretical values with the estimated values, the relative arrival time of the first trace was forced to zero.Then, the estimated relative arrival times were 0 ms, 15 ms, 30.1 ms, and 45.4 ms, respectively.If the relative arrival time of one trace is positive, it indicates that the arrival time of this trace is lagging behind that of the first trace.Otherwise, the arrival time of this trace is ahead of that of the first trace.As shown in Table 3, the estimated values of relative arrival times are highly consistent with the theoretical values.The maximum error is 0.4 ms, which is smaller than one sampling interval (0.5 ms).This result indicates that POC can be used to estimate the relative arrival time.

Field Data Analysis
The studied microseismic event was acquired during a multi-stage hydraulic fracture treatment of a shale gas reservoir.A 12-level geophone array was used to record with a sampling of 0.25 ms in the monitoring well.The distance between the fracture stage and the geophone array was about 1500 m.The length of the record time was 750 ms (or 3000 sampling points).According to the waveforms of the event (Figure 5), it is obvious that the arrival time of the P-wave was at approximately 375 ms and the arrival time of the S-wave was at approximately 700 ms.The P-wave record of this event was intercepted from 325.25 ms to 575 ms (or 1000 sampling points) (Figure 5b).Four methods (AIC, cross-correlation, POC (STFT), and POC (WVD)) were used to estimate the relative arrival times of the P-wave record.
The formula of AIC is defined as: where x(i) is the time series of one trace, n is the length of the trace, and var() denotes the variance of a segment of x.
The estimated values of the first trace were forced to zero.Due to the high SNR of the microseismic event record (Figure 5), the four sets of calculation results had a high degree of consistency (Table 4).These results were used as adjusting value to eliminate the arrival time difference among the traces.As shown in the four subgraphs of Figure 6, it is obvious that the P-wave waveforms of the traces were flattened.Because the accurate arrival times of the field data were not known, the average values of the four sets of results were viewed as the theoretical relative arrival times (Table 4).To simulate abnormal situations of monitoring equipment, we used background noise to replace the 10th trace of the P-wave record, then added three different level Gaussian white noises to the previous P-wave record (Figure 7).Table 5 shows the peak value of the POC (WVD) for the first noisy data (Figure 7a).Each peak value corresponds to one pair of traces.The peak values related to the 10th trace are obviously smaller than others.This result is consistent with the fact that there is not a P-wave signal in the 10th trace.So, the peak value of the POC function was used to detect abnormal situations in the record. Figure 8 shows the POC function of the first trace and the 10th trace, and Figure 9 shows the POC function of the first trace and the 12th trace.Comparing Figure 8 with Figure 9, there is a distinct peak in Figure 9 that does not exist in Figure 8.Therefore, the peak value of the POC function is an important constraint to estimate the relative arrival time of microseismic event records.Figure 8 shows the POC function of the first trace and the 10th trace, and Figure 9 shows the POC function of the first trace and the 12th trace.Comparing Figure 8 with Figure 9, there is a distinct peak in Figure 9 that does not exist in Figure 8.Therefore, the peak value of the POC function is an important constraint to estimate the relative arrival time of microseismic event records.The estimated values of relative arrival times were obtained using the four methods (Table 6).Then, we used the estimated values to eliminate the arrival time differences of the first noisy dataset.As shown in Figure 10, the estimated values of AIC and cross-correlation cannot eliminate the arrival time differences among the traces.The 11th trace in Figure 10c is not consistent with other traces and the flattening effect of POC (WVD) is better than that of POC (STFT).Meanwhile, the estimated values of the 10th trace are random values for all four methods (Table 6).The estimated values of relative arrival times were obtained using the four methods (Table 6).Then, we used the estimated values to eliminate the arrival time differences of the first noisy dataset.As shown in Figure 10, the estimated values of AIC and cross-correlation cannot eliminate the arrival time differences among the traces.The 11th trace in Figure 10c is not consistent with other traces and the flattening effect of POC (WVD) is better than that of POC (STFT).Meanwhile, the estimated values of the 10th trace are random values for all four methods (Table 6).For the second noisy dataset (Figure 7b), the estimated values of the relative arrival times were obtained using the four methods (Table 7) and then used to eliminate the arrival time differences of the second dataset (Figure 11).Comparing Figure 10 with Figure 11, it can be observed that the POC (WVD) is more stable in estimating the relative arrival times when the abnormal situations exist in the monitoring record.For the second noisy dataset (Figure 7b), the estimated values of the relative arrival times were obtained using the four methods (Table 7) and then used to eliminate the arrival time differences of the second dataset (Figure 11).Comparing Figure 10 with Figure 11, it can be observed that the POC (WVD) is more stable in estimating the relative arrival times when the abnormal situations exist in the monitoring record.For the third noisy dataset (Figure 7c), the estimated value of the relative arrival times were obtained using the four methods (Table 8) and then used to eliminate the arrival time differences of the third dataset (Figure 12).Compared with other three results, the result of POC (WVD) is more robust in Figure 12.For quantitatively evaluate the robustness of the four methods, we calculated the differences between the estimated arrival times obtained using the four in addition to the theoretical relative arrival times, and calculated the standard deviations of these differences.
where σ denotes the standard deviation, N denotes the number of traces, t estimate denotes the estimate value of the relative arrival time, and t true denotes the true value of the relative arrival time.
Table 9 shows the standard deviations of the results (excluding the 10th trace) based on the four methods.For the interferences of noise and the code wave of the microseismic event, there are obvious errors for the estimated arrival times based on AIC.For the cross-correlation, it is hard to handle abnormal situation of the monitoring record.Due to the advantage of using the peak values of the POC functions as weighting coefficients to improve the consistency, POC (STFT), and POC (WVD) were found to have better stability in estimating the relative arrival times.By using a time window, STFT can obtain the time-frequency representation of the data.The resolution of STFT is fixed due to fixed sliding window length.Wigner-Ville distribution has the ideal time-frequency resolution [47].For the difference of the time-frequency resolution, the estimation of the relative arrival time based on WVD is better than that based on STFT.

Discussion and Conclusions
In order to handle the low SNR record, real-time monitoring, and abnormal situations of monitoring equipment, we proposed a new method based on phase-only correlation to estimate the relative arrival times of microseismic events.We initially used time-frequency transform to suppress noise.Furthermore, the position of the peak value of the POC function of the time-frequency spectra was used to indicate the relative arrival times between traces.Since the peak value of the POC function can indicates the similarity between traces, we finally used the peak values as weighting coefficients to improve the consistency of the estimated relative arrival times and eliminate the effect of abnormal situations of the monitoring equipment.Beacause the peak value and the position of the peak value of the POC function are automatically obtained, the whole processing does not require manual intervention.Therefore, the proposed method can meet the requirements of real-time monitoring.
According to the synthetic data and field data employed in this study, the proposed method is valid.In particular, it was found that POC (WVD) is the most stable method when abnormal situations arise in the monitoring equipment.The advantages of the proposed method include the following.
(1) Due to the time-frequency sparseness of the microseismic event, the time-frequency representation of the event record can be used to suppress random noise.The time-frequency resolution of WVD is better than that of STFT.So, POC (WVD) is more stable than POC (STFT).(2) Excluding the effect of the amplitude components, the POC function merely calculates the difference between the phase components.It is useful to process two signals with a significant energy difference.(3) The peak values of the POC can be used as weighting coefficients of linear equations to improve the consistency of the relative arrival time among the traces.

Figure 1 .
Figure 1.The flowchart of the proposed method for estimating relative arrival times.The procedure of this method starts with time-frequency transform, followed by phase-only correlation (POC) and consistency processing.

Figure 2 .
Figure 2. The waveforms of the synthetic data.

Figure 3 .
Figure 3.The time-frequency representations of two traces: (a) the result of time-frequency transform of Trace 1; (b) the result of time-frequency transform of Trace 4.

Figure 4
Figure 4 shows the POC between Trace 1 and Trace 4.There is a distinct peak in the POC function.The peak value is located at −45 of the time lag axis.The minus sign indicates that the arrival time of Trace 1 is ahead of that of Trace 4. The location of the peak value directly indicates the arrival time difference between Trace 1 and Trace 4. This result is consistent with the theoretical value.Meanwhile, it is clear that the peak value of the POC function in Figure 4 is less than 1.This result is attributed to the effect of the low-pass filtering for the cross-phase spectrum.As shown in Figure4, the range of the time lag axis is from −75 to 75, which is half the length of a single trace.If the arrival time difference between two traces were bigger than half the length of a single

Figure 4 .
Figure 4.The POC of the first trace and the fourth trace.The peak value of the POC is located in −45 of the time lag axis.

Figure 4 .
Figure 4.The POC of the first trace and the fourth trace.The peak value of the POC is located in −45 of the time lag axis.

Figure 5 .
Figure 5.The waveforms of the microseismic event.(a) the P-wave and S-wave waveforms; (b) the P-wave record, which is a segment of the microseismic event from 325.25 ms to 575 ms.

6 .
Using the calculation results of relative arrival times as the adjusting value to eliminate the arrival time difference among traces.(a-d) Akaike Information Criterion (AIC), cross-correlation, phase-only correlation (short time Fourier transform) (POC (STFT)), and phase-only correlation (Wigner-Ville distribution) (POC (WVD)) were used to estimate the relative arrival times of the P-wave record.

Figure 7 .
Figure 7.The 10th trace of the P-wave record is replaced by background noise, and the P-wave record is aminated by different levels of Gaussian white noises.(a) The signal-to-noise ratio of the first noisy data is 5 dB; (b) the signal-to-noise ratio of the second noisy data is 0 dB; (c) the signal-to-noise ratio of the third noisy data is −2 dB.

Figure 8 .
Figure 8.The POC function between the first trace and the 10th trace.

Figure 8 .
Figure 8.The POC function between the first trace and the 10th trace.

Figure 9 .
Figure 9.The POC function between the first trace and the 12th trace.

Table 6 .
The estimated values of the relative arrival times of the first noisy dataset (signal-to-noise ratio (SNR) = 5 dB).

Figure 10 .
Figure 10.Using estimation values of the relative arrival times to eliminate the arrival time difference among the traces of the first noisy dataset.(a) The result based on AIC; (b) the result based on cross-correlation; (c) the result based on POC (STFT); (d) the result based on POC (WVD).

Figure 9 .
Figure 9.The POC function between the first trace and the 12th trace.

Figure 10 .
Figure 10.Using estimation values of the relative arrival times to eliminate the arrival time difference among the traces of the first noisy dataset.(a) The result based on AIC; (b) the result based on cross-correlation; (c) the result based on POC (STFT); (d) the result based on POC (WVD).

Figure 11 .
Figure 11.Using estimation values of the relative arrival times to eliminate the arrival time difference among the traces of the second noisy data.(a) The result based on AIC; (b) the result based on cross-correlation; (c) the result based on POC (STFT); (d) the result based on POC (WVD).

Figure 12 .
Figure 12.Using estimation values of the relative arrival times to eliminate the arrival time difference among the traces of the third noisy data.(a) The result based on AIC; (b) the result based on cross-correlation; (c) the result based on POC (STFT); and (d) the result based on POC (WVD).

Table 1 .
Symbols present in this article.

Table 2 .
The arrival time differences of all the pairs of traces.

Table 2 .
The arrival time differences of all the pairs of traces.

Table 3 .
The comparison of the theoretical values and the estimated values of the relative arrival times.

Table 4 .
The relative arrival times of the microseismic event using the four methods.

Table 7 .
The estimated values of the relative arrival times of the second noisy dataset (SNR = 0 dB).

Table 6 .
The estimated values of the relative arrival times of the first noisy dataset (signal-to-noise ratio (SNR) = 5 dB).

Table 7 .
The estimated values of the relative arrival times of the second noisy dataset (SNR =

Table 8 .
The estimated values of the relative arrival times of the third noisy dataset (SNR = −2 dB).

Table 9 .
The standard deviations of the relative arrival time of the microseismic event under different signal-to-noise ratios, while one abnormal trace exists in the record.