- Email: [email protected]
Automatic microseismic event detection by band-limited phase-only correlation
Accepted Manuscript Automatic Microseismic Event Detection by Band-Limited Phase-Only Correlation Shaojiang Wu, Yibo Wang, Yi Zhan, Xu Chang PII: DOI: Reference:
S0031-9201(16)30214-X http://dx.doi.org/10.1016/j.pepi.2016.09.005 PEPI 5965
To appear in:
Physics of the Earth and Planetary Interiors
Received Date: Revised Date: Accepted Date:
2 October 2015 20 September 2016 21 September 2016
Please cite this article as: Wu, S., Wang, Y., Zhan, Y., Chang, X., Automatic Microseismic Event Detection by Band-Limited Phase-Only Correlation, Physics of the Earth and Planetary Interiors (2016), doi: http://dx.doi.org/ 10.1016/j.pepi.2016.09.005
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
Automatic Microseismic Event Detection by Band-Limited Phase-Only Correlation
2 3 4
Shaojiang Wu 1, 2, Yibo Wang1,*, Yi Zhan3, Xu Chang 1
5 6
1. Key Laboratory of Shale Gas and Geoengineering, Institute of Geology and
7
Geophysics, Chinese Academy of Sciences, Beijing 100029, China
8
2. University of Chinese Academy of Sciences, Beijing 100049, China
9
3. Geophysical Research Institute, Bureau of Oil Geophysical Prospecting, CNPC,
10
Zhuozhou, Hebei Province 072750, China
11
* Corresponding author. Tel. /fax: +86 1082998624.
12
E-mail address: [email protected] (Y. Wang).
13 14
15
ABSTRACT
16
Identification and detection of microseismic events is a significant issue in source
17
locations and source mechanism analysis. The number of the records is notably large,
18
especially in the case of some real-time monitoring, and while the majority of
19
microseismic events are highly weak and sparse, automatic algorithms are
20
indispensable. In this study, we introduce an effective method for the identification
21
and detection of microseismic events by judging whether the P-wave phase exists in a
22
local segment from a single three-component microseismic records. The new judging
23
algorithm consists primarily of the following key steps: 1) transform the waveform
24
time series into time-varying spectral representations using the S-transform; 2)
25
calculate the similarity of the frequency content in the time-frequency domain using
26
the phase-only correlation function; and 3) identify the P-phase by the combination
27
analysis between any two components. The proposed algorithm is compared to a
28
similar approach using the cross-correlation in the time domain between any two
29
components and later tested with synthetic microseismic datasets and real
30
field-recorded datasets. The results indicate that the proposed algorithm is able to
31
distinguish similar and dissimilar waveforms, even for low signal noise ratio and
32
emergent events, which is important for accurate and rapid selection of microseismic
33
events from a large number of records. This method can be applied to other
34
geophysical analyses based on the waveform data.
35 36
INTRODUCTION
37
Generally, the microseismic activity is induced by small magnitude hydraulic
38
fracturing (Maxwell and Urbancic, 2001) and acquired in low signal to noise (SNR)
39
environments. Meanwhile, the number of raw records is very large, especially from
40
real-time monitoring, while the majority of microseismic events seldom occur. Prior
41
to the localization and mechanism analysis of the source, the identification and
42
detection of events become an important challenge. Normally, the existing
43
event-detecting processes can be split into two steps: first, detection of the segment of
44
the records containing the seismic events, and then accurately identifying and picking
45
the target phase. The first stage of detection for a microseismic event is
46
time-consuming, as it requires examination of a large number of raw records; however,
47
all subsequent steps require less resources, as they are based on the previous
48
selections.
49 50
The correlation-type and energy analysis are two principal methods for event
51
detection. The energy analysis is a broadly applicable method due to few assumptions
52
about the data. The most-used approach is defined as the energy ratio of the
53
short-time-average to the long-time-average (STA/LTA) (Earle and Shearer 1994)
54
within moving time windows, which can detect the appearance of seismic events
55
when the ratio exceeds the given thresholds. However, this method is not sensitive to
56
the weak low-magnitude events that have loud background noise, similar to those
57
which occur during hydraulic fracturing. The correlation-type methods, such as
58
cross-correlation within the time domain, detect similar events to a referenced event,
59
known as the master event, by a calculated correlation coefficient to it (Gibbons and
60
Ringdal, 2006; Song et al. 2010). However, the referenced master template waveform
61
is difficult to make when the majority of seismic events occur from unknown and
62
dissimilar sources, and the calculated similarity coefficient is inaccurate without
63
sufficient data, meaning that a long data collection window is required for accuracy.
64 65
Both of the described methods are amplitude-biased, due to the strong dependence on
66
the waveforms and relative amplitudes (Schimmel and Paulssen, 1997). Therefore,
67
results may be ambiguous due to the waveform distortion or weak signals. Some
68
amplitude-normalized cross-correlation (Neidell & Taner 1971) and semblance-based
69
stacks (Taner and Koehler, 1969) are widely used. These methods can remove the
70
disturbing influence of energetic features by some additional pre-processing or
71
post-processing steps, such as normalization, but the waveform coherence may
72
deteriorate. Therefore, several instantaneous phase-based measurements are
73
introduced, which are amplitude-unbiased. The phase-weighted cross-correlations
74
(Schimmel, 1999) and phase-weighted stacks (Schimmel, 1997) are designed for
75
noise attenuation or signal extraction. They are weighted by the coherency of the
76
instantaneous phase, which implicitly contains the waveform information through the
77
Hilbert Transform.
78 79
Normally, when the seismic waveform is resolved into the time-varying frequency
80
components, the time-frequency spectrum can reveal some underlying periodicities
81
and express more abundant characteristics of the seismic waveform. In practice, the
82
useful seismic signals concentrate on some limited frequency band, which is
83
meaningful for noise attenuation and signal extraction. The phase-only correlation
84
(POC) is a simple and robust technique to evaluate the similarity of two images. The
85
technique has been successfully applied to electronics and communication
86
engineering, such as in fingerprint matching (Ito et al., 2004) and subpixel image
87
registration (Takita et al., 2004). The technique has also been expanded into the
88
geophysical exploration field, where POC is used to identify similar waveforms, such
89
as the events among the aftershocks of an earthquake (Moriya, 2011, Zhang et al,
90
2015).
91 92
In this study, we apply the POC for automatic event detection within a single
93
three-component microseismic record. The process workflow is similar to the
94
well-known polarization filtering approach (Jurkevics, 1988), except for identifying
95
similar waveforms with POC based on the transformed time-frequency spectrum. The
96
method first splits the raw records into some contiguous segments, and then it detects
97
the segment containing the event by judging whether the P-wave phase exists. For the
98
improved judging process of the proposed method, the following key steps are
99
necessary: decompose the seismic waveform into the time-varying frequency
100
components with the S-transform, which can retain the absolute phase of each
101
frequency and offer a high resolution, calculate the similarity of the frequency content
102
in the time-frequency domain using the POC function, and identify the P-phase by the
103
combination analysis between any two components. The proposed algorithm, together
104
with the polarization filtering approach using the cross-correlation in the time domain,
105
is tested with synthetic microseismic datasets and real field-recorded datasets.
106 107
METHODOLOGY
108 109
S-Transform
110
The S-transform enables a better resolution of low frequency components and time
111
resolution of high frequency signals by employing frequency-dependent windows,
112
similar to the wavelet transform. Meanwhile, it retains the absolute phase of each
113
frequency (Pinnegar and Mansinha, 2003) for multiresolution analysis and detection
114
of events in time series.
115 116
In normal signal processing, we need to determine the frequency spectrum, , of
117
the time series, , of any component of the raw records, by using the following
118
Fourier transform:
119
= . (1)
120 121
With the additional frequency dependence of the analyzing window, the S-transform
122
of the time series, , can be defined
123
, = − , , ≠ 0, (2)
124
, = 0 = → ! , (3)
!
125 126
where − , is a Gaussian function, which maintains the center at time, τ, of
127
the selected Gaussian window and a standard deviation proportional to ## for the
128
localized spectrum.
129
− , =
||
% √
'(! )'*! !+!
, , > 0, (4)
130 131
where k is the scaling factor, which controls the time-frequency resolution by the
132
number of oscillations in each Gaussian window.
133 134
Aimed to simplify the implementation by computers, equation (3) can also expressed
135
as a multiplication of spectrum:
136
, =
/
+
!1!
2! +! (! 23
/, ≠ 0, (5)
137 138
where / + is the shifted spectrum of the Fourier transform. The Gaussian
139
function in the above integrand equation has a frequency-dependent bandwidth. More
140
specifically, the phase-shifted time signal, , is bandpassed at the center frequency
141
with the deviation 4 =
%||
by the Gaussian window function.
142 143
The time-frequency representation of the S-transform is unique and invertible
144
(Schimmel and Gallart, 2005). The local time-frequency spectra, , , can be
145
transformed back to the time domain, since the Gaussian window for the S-transform
146
must satisfy the condition
147
− , = 1, (6)
148 149
and the time series, u789 t, is easily reconstructed as:
150
;<= = , . (7)
151 152
Band-limited POC
153
Next, we considered the transformed time-frequency spectra , and ,
154
of any two of the components of the raw records as two images, > , > and
155
> , > , respectively. Because the only similarity of the frequency content in the
156
time-frequency domain is needed, the two axes of time, , and frequency, , were
157
equal, familiar with the general 2-D pictures, and rewritten as > and > ,
158
respectively. Conveniently, the size of each of the images was defined as ? × ? ,
159
where the index ranges were > = −A , ⋯ , A and > = −A , ⋯ , A , and
160
therefore, ? = 2A + 1 and ? = 2A + 1 . Then, the 2-D discrete Fourier
161
transform of the images, > , > and > , > , were described as the following
162
equations (Ito et al., 2004):
163 F
F
% G
% G!
= L , , , MNI %I ,%! , (8)
F
F
% G
% G!
= L , , , MN! %I ,%! , (9)
164
D , , , = ∑GIIH FI ∑G!!H F! > , > JKII I JK!!
165
D , , , = ∑GIIH FI ∑G!!H F! > , > JKII I JK!!
166 M
!1 OI
, and JK! =
M
!1 O!
167
where , = −A , ⋯ , A , , = −A , ⋯ , A , JKI =
168
L , , , and L , , , are the amplitude components of the two images,
.
169
respectively, and P , , , and P , , , are the relative phase components.
170
Normally, the cross-spectrum, QR, , , , between D , , , and D , , , can be
171
defined as:
172
RRR% ,% S % ,% S QR, , , = |SI %I ,%! SRRR!%I ,%! | = MN%I ,%! , (10) I
I
!
!
I
!
173 174
where DR , , , represents the complex conjugate of D, , , , and P, , , is
175
the difference of the phases P , , , and P , , , . The POC correlation
176
function is denoted as the 2-D inverse discrete Fourier transform of the
177
cross-spectrum, QR , , , , which can be expressed as
178
T̅ > , > =
KI K!
F! %IGI % G I R ∑F JK! ! ! , (11) GI H FI ∑G! H F! Q , , , JKI
179 180
where T̅ > , > is the value of the POC function, which offers a range from 0 to 1.
181
The peak value of the function can indicate the similarity of the two images and the
182
position provides the translational displacement ∆> , ∆> between the two images:
183
∆> , ∆> = −> , −> . (12)
184 185
Additionally, the function gives a distinct sharp peak of value for two similar images,
186
and it decreases significantly when they are not similar, which demonstrates higher
187
discrimination capabilities than the traditional correlation. Meanwhile, as a phase
188
indication of two images, it works well against the images shifts and brightness
189
changes (Ito et al., 2004).
190
191
In practice, the significant information in the 2-D DFT of the given image is usually
192
clustered in a limited frequency band. Hence, some phase components in the higher or
193
lower frequency domain are meaningless or not reliable. The basic idea of the
194
band-limited POC is to improve the matching performance by eliminating
195
meaningless frequency components before the calculation of the cross-phase
196
spectrum, QR, , , . So, the band-limited POC function can be defined as
197
T̅W XYZ > , > = W
I W!
[! %I GI % G I R ∑[ JW! ! ! , (13) GI H [I ∑G! H [! Q , , , JWI
198 199
ehere k and k are the given inherent frequency bands with the ranges of
200
, = −\ , ⋯ , \ and , = −\ , ⋯ , \ , and 0 ≤ \ ≤ A and 0 ≤ \ ≤ A . In
201
the processing, the parameters \ and \ can be selected with an effective coverage
202
range of the meaningful phase components.
203 204
Combination analysis
205
Theoretically, the three-component seismogram can represent the motion in three
206
orthogonal directions of a ground detector, containing two in the horizontal plane and
207
one in the vertical. The observed seismogram is often represented as the convolution
208
of the source function, the transmission response of subsurface media, or the receiver
209
function. The source refers to a series of vibrations in the Earth due to an earthquake
210
or explosion. The subsurface media are assumed to be heterogeneous, where some
211
reflection and transmission takes place and make wave encounters discontinuous.
212
Some experiments (Morlet et al., 1982) have demonstrated that the above
213
phenomena of wave propagation can considered as some "scattering." The distortion
214
waveforms due to the scattering are frequency-dependent and thus isolated in certain
215
scales (Anant and Dowla, 1997). Therefore, the arrival waves are normally fairly
216
similar, with possible differences in the very high frequencies (very low scales)
217
(Anant and Dowla, 1997). Based on the frequency-dependent nature of the scattering,
218
the arrivals can be viewed in terms of the type of decomposition, and they are only
219
present in some scales, not consistently across all scales.
220 221
By definition, the P-phase is represented by the longitudinal or compressional motion,
222
while the S phase is represented by the transverse motion. The compressional body
223
wave is often linearly polarized. Therefore, it would be more accurate to detect P
224
arrivals based on the linear polarization and the frequency dependence. In our method,
225
the P-wave phase can be identified under the condition of simultaneous high
226
correlation values between any two components within the raw three-component
227
seismograms records.
228 229
EXAMPLES
230
As described above, the proposed workflow is much the same as the well-known
231
polarization filtering approach (Jurkevics, 1988), or other similar methods, except that
232
it identifies the similar waveforms using the POC based on the transformed
233
time-frequency spectrum, replacing the cross-correlation in the time domain or other
234
amplitude-based methods. As shown in Figure 1, the method first splits the raw
235
records into some contiguous segments, and then it detects the segment containing the
236
event by judging whether the P-wave phase exists. For improved judging processes,
237
the examples with similar and dissimilar waveforms were tested as the two
238
possibilities in each loop. The proposed algorithm-based band-limited POC and the
239
1-D cross-correlation in time domain of the raw records were compared mainly for
240
distinguishing capabilities between the similar and dissimilar waveforms. The
241
window length and the threshold were two important parameters in the proposed
242
method. In the proposed method, the lag of any two P-phase arrival times within the
243
three-component geophone was almost small. Meanwhile, the longer window may
244
contain the target data of the P-phase, together with other phases, such as the coda
245
wave, which may affect the accuracy of judgment. Therefore, in order to improve
246
the accuracy and robustness, the window length of the time series is recommended to
247
be between one-half and one waveform length. The window moves from left to right
248
with an overlapped length of window in each loop. For some batch processing, such
249
as real-time monitoring, the threshold can be obtained by pre-processing experiments,
250
which mainly depend on the performance of the hydraulic fracturing, and also by
251
measuring the background noise in the acquired environment.
252 253
(The position of Figure 1)
254 255
Synthetic microseismic datasets
256
A two-layer model was used for the synthetic microseismic datasets modeling. The
257
P-wave velocity of the model was set to 1500, 3000 m/s, the S-wave velocities to
258
1500, and 3000 m/s , and the density to 11, 2.5 kg/mh . Approximately 10-level
259
three-component geophones were placed in the vertical observation well, which were
260
spaced at intervals of 25 m within the array. The microseismic source was located
261
atx, y, z = 500, 300, 4000, with a dominant frequency of 400 Hz. The data were
262
recorded at 280 ms with a sampling rate of 4 samples per millisecond.
263 264
(The position of Figure 2)
265 266
The recorded three-component seismograms within an array are displayed in Figure
267
2(a). In practice, the microseismic events within real field records were weak and had
268
superimposing noise from background noise, such as vibrations from machinery. To
269
better simulate real data acquisition, Gaussian noise (approximately 50% of the
270
maximum signal amplitude) was added to the synthetic records. Figure 2(b)
271
demonstrates that the noise affected the microseismic records, which would present a
272
challenge to detect the events reliably using some automatic processes. The window
273
length of any selected segments was set to 32 ms.
274 275
(The position of Figure 3-5)
276 277
The first example demonstrates the segment from the black box in Figure 2(b), which
278
contains the major events from 60 to 92 ms. Figure 3(b)-(d) demonstrates the
279
time-varying spectral representation of the components x, y, and z, respectively. The
280
transformed spectral representation captures the features of seismic phase arrivals,
281
especially the predominant frequency distribution. Then, the frequency contents from
282
300Hz to 600Hz were selected as the significant information for the following steps,
283
which means that both the very low and high frequency noise were the most excluded.
284
Figure 4(a)-(c) demonstrates the POC coefficients of the coupled components x-y, x-z,
285
and y-z, respectively. The high values were mostly distributed at the center of the
286
images, which had a high correlation value with a frequency lag of 0 and a time lag of
287
0. The peak of amplitude was very sharp, which makes it easy to detect with given
288
threshold. Then, we calculated the similarity using the cross-correlation function in
289
time domain with the same window length as the records. Figure 5 demonstrates the
290
POC coefficient (red curve) at the frequency lag 0 and the cross-correlation
291
coefficient in time domain (blue curve) corresponding to Figure 4(a)-(c). Both
292
methods can indicate the high correlation values at the time lag 0. However, the
293
cross-correlation coefficient in the time domain demonstrated a lower resolution
294
compared to the POC with the short window length of the selected time series because
295
the calculated similarity coefficient of the cross-correlation in the time domain would
296
be more accurate with sufficient data, meaning a longer window size.
297
(The position of Figure 6-8)
298 299
The second example demonstrates the segment from 0 to 32 ms, from the red box in
300
Figure 2(b), which mainly contains random noise without appreciable signals. Figure
301
6(b)-(d) are the time-varying spectral representation of the components x, y, and z,
302
respectively. Different from the previous example, the transformed spectral
303
representations have a large difference in the frequency distribution. Then, the same
304
parameters were set for the calculation of the POC. Figure 7(a)-(c) demonstrates the
305
POC coefficients of the coupled components x-y, x-z, and y-z. The highest
306
correlations were not at the center of the image with the frequency lag 0 and the time
307
lag 0, and the values were reduced considerably. The values of the dissimilar images
308
were relatively small compared the peak value of the similar images, which was easy
309
to distinguish. In the next comparison with the cross-correlation coefficient in the
310
time domain, both methods were stable and demonstrated small coefficient values,
311
indicating the dissimilar waveforms of the pair of segments. However, the normalized
312
cross-correlation coefficient in the time domain distributes in a relatively wider range.
313 314
Real field-recorded datasets
315
Figure 9(a) demonstrates the continuous microseismic records for approximately 3000
316
ms with a sampling rate of 4 samples per millisecond. The 4th level of records
317
demonstrated in Figure 9(b) was used to show an actual process based the proposed
318
method. And the window length was also set to 32 ms . A large-magnitude
319
microseismic event (the black arrow B) was noticeable, which could be easily
320
selected by certain correlation-based methods. However, some potential small
321
magnitude events in the background noise are difficult to detect correctly.
322
323
(The position of Figure 9-10)
324 325
As a result of above processing algorithm, three segments with potential events were
326
detected in Figure 10(a) with red traces and magnified to more detail in Figure
327
10(b)-(d), separately. The segment from 2240 to 2272 ms in Figure 10(c) contains the
328
large-magnitude microseismic event mentioned before. But, the other two segments
329
containing potential events from 1400 to 1432 ms in Figure 10(b) and from 2930 to
330
2962 ms in Figure 10(d) were very weak, which still needs investigation. The
331
following example demonstrated the segment from 2240 to 2272 ms in Figure 10(c)
332
in more detail. Figure 11(a)-(c) were the time-varying spectral representations of
333
components x, y, and z, respectively. Figure 12(a)-(c) demonstrates the POC
334
coefficients of the coupled component x-y, x-z, and y-z. Figure 13 demonstrates the
335
normalized cross-correlation coefficient in the time domain (blue curve) and the POC
336
function coefficient (red curve) at the frequency lag 0. The cross-correlation
337
coefficient had multiple blunt peaks, which were close to the max value of the POC
338
coefficient. However, it is difficult to indicate the relative difference of the arrival
339
times of the compared waveforms in the selected segments with the constant
340
threshold.
341 342
(The position of Figure 11-13)
343 344
DISCUSSIONS
345
This study introduced an effective method for the identification and detection of the
346
microseismic events by judging whether the P-wave phase exists in a local segment
347
from single three-component microseismic records. The improved judging process
348
contained the following key steps: first, the waveform time series was correctly
349
transformed into time-varying spectral representations using the S-transform. Next,
350
the similarity of the frequency content in the time-frequency domain was calculated
351
using the POC function, and the P-phase was identified by the combination analysis
352
between any two components. The proposed algorithm was compared to similar
353
approaches using the cross-correlation in time domain between any two components,
354
and it was capable of distinguishing between the similar and dissimilar waveforms.
355 356
The proposed workflow is similar to the well-known polarization filtering approach,
357
except that it identifies the similar waveforms using the POC based on the
358
transformed time-frequency spectrum, replacing the cross-correlation in the time
359
domain, or other amplitude-based methods. For the proposed workflow, the lag of any
360
two P-phase arrival times within the three-component geophone is almost small.
361
Meanwhile, the longer window may contain the target data of the P-phase, which,
362
together with other phases, such as coda wave, may affect the accuracy of judgment.
363
So, the window length of the time series is recommended to be between one-half and
364
one waveform length for improved accuracy and robustness. However, the calculated
365
similarity coefficient by the cross-correlation coefficient in the time domain would be
366
more accurate with sufficient data, meaning a longer window size. Therefore, the
367
short window size or short data may result in lower resolution coefficients with the
368
cross-correlation in the time domain, where it is difficult to indicate the relative lag of
369
the arrival times of the compared waveforms in the selected segments. Therefore, the
370
POC-based workflow with high resolution lag cross-correlation coefficient is
371
a great choice for the proposed workflow, which is meaningful for accurate
372
identification and detection of the microseismic events. In addition, distinguishing
373
between the similar and dissimilar waveforms with high resolution was easily
374
achieved with some constant thresholds, which is ideal for automating a process.
375 376
However, processing with POC has a low probability of providing false indications,
377
mainly because the calculated coefficient in the time-frequency domain is sensitive
378
for some unreasonable frequency content. In consideration of the computation
379
efficiency and stability, it is suggested that the practical process should recombine and
380
combine the bright side of the traditional correlation-type and energy analysis
381
methods and the POC method, especially for rapid field analysis and hydro-fracture
382
monitoring.
383 384
CONCLUSIONS
385
This work presents an effective method for the detection of the microseismic events
386
by judging if the P-wave phase exists in a local segment using the POC function. The
387
results of the synthetic data and the real hydraulic-fracture case show that the
388
proposed method is sufficiently stable to distinguish the similar and dissimilar
389
waveforms, even for low SNR, which means that it can reliably handle the two
390
possibilities during the judging algorithm. It is important to identify and detect the
391
microseismic events within a large number of records accurately and rapidly, and this
392
method can be applied to other geophysical analyses based on the waveform data.
393 394
ACKNOWLEDGMENTS
395
This research was funded by the National Natural Science Foundation of China (Grant
396
no. 41230317) and the National Basic Research Program of China (Grant no.
397
2015CB258500).
398
399
REFERENCES
400
Anant K S, Dowla F U. Wavelet Transform Methods for Phase Identification in
401
Three-Component Seismograms [J]. Bulletin of the Seismological Society of America,
402
1997, 87(6):1598-1612.
403 404
Earle P S, Shearer P M. Characterization of global seismograms using an
405
automatic-picking algorithm [J]. Bulletin of the Seismological Society of America,
406
1994, 84(2): 366-376.
407 408
Gibbons S J, Ringdal F. The detection of low magnitude seismic events using
409
array-based waveform correlation [J]. Geophysical Journal International, 2006, 165(1):
410
149-166.
411 412
Ito K, Nakajima H, Kobayashi K, et al. A Fingerprint Matching Algorithm Using
413
Phase-Only Correlation [J]. Ieice Transactions, 2004: 682-691.
414 415
Jurkevics A. POLARIZATION ANALYSIS OF THREE-COMPONENT ARRAY
416
DATA[J]. Bulletin of the Seismological Society of America, 1988, 78(5): 1725-1743
417 418
Maxwell S C, Urbancic T I. The role of passive microseismic monitoring in the
419
instrumented oil field [J]. Geophysics, 2012, 20(6): 636-639.
420
Morlet J, Arens G, Fourgeau E, et al. Wave propagation and sampling theory—Part I:
421
Complex signal and scattering in multilayered media [J]. Geophysics, 2012, 47,
422
203-221.
423 424
Moriya H. Phase-only correlation of time-varying spectral representations of
425
microseismic data for identification of similar seismic events [J]. Geophysics, 2011,
426
76(6).
427 428
Neidell N S, Taner M T. Semblance and other coherency measures for multichannel
429
data [J]. Geophysics, 1971, 36(3): 482-497.
430 431
Pinnegar C R, Mansinha L. The S-transform with windows of arbitrary and varying
432
shape [J]. Geophysics, 2003, 68(1):381.
433 434
Schimmel M, Gallart J. The inverse S-transform in filters with time-frequency
435
localization [J]. IEEE Transactions on Signal Processing, 2005, 53(11): 4417-4422.
436 437
Song F, Kuleli H S, Toksoz M N, et al. An improved method for
438
hydrofracture-induced microseismic event detection and phase picking [J].
439
Geophysics, 2010, 75, A47–A52.
440 441
Stockwell R G, Mansinha L, Lowe R P, et al. Localization of the complex spectrum:
442
the S transform [J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.
443 444
Takita K, Aoki T, Sasaki Y, et al. High-accuracy subpixel image registration based on
445
phase-only correlation [J]. Ieice Transactions, 2003: 1925-1934.
446 447
Zhang L, Wang Y, Chang X, et al. Wigner distribution based phase-only correlation[J].
448
Journal of Applied Geophysics, 2015: 277-282.
449 450
Figure 1. The workflow of the proposed algorithm of identification and detection of
451
the microseismic events
452
453 454
(a)
455 456
(b)
457
Figure 2. (a) Synthetic microseismograms produced by 3D elastic wave equation
458
modeling. Every triplet of x (red), y (blue), and z (green) components is plotted for
459
the 3-component seismograms of each geophone. (b) Gaussian noise was added with
460
a given SNR of approximately 3. The segments in the red and black boxes are used
461
for the following tests of identification capability.
462
463 464
(a)
465 466
(b)
467 468
(c)
469
(d)
470
Figure 3. (a) The selected segment contains the major events with three components
471
from the black box in Figure 2(b); (b)-(d) time-varying spectral representation of
472
components x, y, and z, respectively.
473
474 475
(a)
476 477
(b)
478 479
(c)
480
Figure 4. (a) - (c) The POC function of the couple components x-y, x-z, and y-z,
481
respectively. The sharp peaks of amplitude are significant, which are easy to be
482
detected with some given threshold.
483 484
Figure 5. The comparison of the POC function coefficient (red curve) at the frequency
485
lag 0 and the normalized cross-correlation coefficient in time domain (blue curve)
486
corresponding to Figure 4(a)-(c). Both methods can indicate the max correlation
487
coefficients at the same point. However, the cross-correlation coefficient in the time
488
domain is less resolved.
489
490 491
(a)
492 493
(b)
494 495
(c)
496 497
(d)
498
Figure 6. (a) The selected segment contains the major events with three components,
499
from the red box in Figure 2(b); (b)-(d) time-varying spectral representation of
500
components x, y, and z, respectively.
501
502 503
(a)
504 505
(b)
506 507
(c)
508
Figure 7. (a) - (c) The phase-only correlation function of the couple components x-y,
509
x-z, and y-z, respectively. The high values show a scattered distribution, which is
510
treated as a weak correlation.
511 512
Figure 8. The comparison of the phase-only correlation function coefficient (red curve)
513
at the frequency lag 0 and the normalized cross-correlation coefficient in the time
514
domain (blue curve) corresponding to the Figure 7(a)-(c). Both methods show stable
515
results with small values in a tight distribution.
516
517 518
(a)
519 520
(b)
521
Figure 9. (a) The 7-level seismograms of the downhole geophone array with some
522
triplets of x (red), y (blue), and z (green) components. (b) The 4th level of the records
523
in Figure 9 is selected for the following testing. The red box demonstrates a split
524
segment, which moves from left to right with an overlapped half length.
525
A
526 527
(a)
528 529
(b)
530 531
(c)
B
C
532 533
(d)
534
Figure 10. (a) Three detected segments with potential events extracted in Figure 10(a)
535
with red traces; (b)-(d) the extracted and magnified traces in more details. The
536
segments containing the large-magnitude microseismic event (arrow B) and some
537
potential weak events (arrow A and C) are detected, though they will be identified by
538
some additional information.
539
540 541
(a)
542 543
(b)
544 545
(c)
546
Figure 11. The analysis of the selected segment from 2240 to 2272 ms from the 4th
547
level in Figure 10 (b); (a)-(c) time-varying spectral representation of components x, y,
548
and z, respectively.
549
550 551
(a)
552 553
(b)
554 555
(c)
556
Figure 12. (a) - (c) The phase-only correlation function of the couple components x-y,
557
x-z, and y-z, respectively, which indicate the maximum correlation coefficients with
558
the relative peak values.
559 560
Figure 13. The comparison of the phase-only correlation function coefficient (red
561
curve) at the frequency lag 0 and the normalized cross-correlation coefficient in the
562
time domain (blue curve) corresponding to Figure 11(a)-(c). The cross-correlation
563
coefficient in the time domain has multiple blunt peaks, which make it difficult to
564
indicate the relative difference of the arrival times of the compared waveforms in the
565
selected segments.
S0031-9201(16)30214-X http://dx.doi.org/10.1016/j.pepi.2016.09.005 PEPI 5965
To appear in:
Physics of the Earth and Planetary Interiors
Received Date: Revised Date: Accepted Date:
2 October 2015 20 September 2016 21 September 2016
Please cite this article as: Wu, S., Wang, Y., Zhan, Y., Chang, X., Automatic Microseismic Event Detection by Band-Limited Phase-Only Correlation, Physics of the Earth and Planetary Interiors (2016), doi: http://dx.doi.org/ 10.1016/j.pepi.2016.09.005
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
Automatic Microseismic Event Detection by Band-Limited Phase-Only Correlation
2 3 4
Shaojiang Wu 1, 2, Yibo Wang1,*, Yi Zhan3, Xu Chang 1
5 6
1. Key Laboratory of Shale Gas and Geoengineering, Institute of Geology and
7
Geophysics, Chinese Academy of Sciences, Beijing 100029, China
8
2. University of Chinese Academy of Sciences, Beijing 100049, China
9
3. Geophysical Research Institute, Bureau of Oil Geophysical Prospecting, CNPC,
10
Zhuozhou, Hebei Province 072750, China
11
* Corresponding author. Tel. /fax: +86 1082998624.
12
E-mail address: [email protected] (Y. Wang).
13 14
15
ABSTRACT
16
Identification and detection of microseismic events is a significant issue in source
17
locations and source mechanism analysis. The number of the records is notably large,
18
especially in the case of some real-time monitoring, and while the majority of
19
microseismic events are highly weak and sparse, automatic algorithms are
20
indispensable. In this study, we introduce an effective method for the identification
21
and detection of microseismic events by judging whether the P-wave phase exists in a
22
local segment from a single three-component microseismic records. The new judging
23
algorithm consists primarily of the following key steps: 1) transform the waveform
24
time series into time-varying spectral representations using the S-transform; 2)
25
calculate the similarity of the frequency content in the time-frequency domain using
26
the phase-only correlation function; and 3) identify the P-phase by the combination
27
analysis between any two components. The proposed algorithm is compared to a
28
similar approach using the cross-correlation in the time domain between any two
29
components and later tested with synthetic microseismic datasets and real
30
field-recorded datasets. The results indicate that the proposed algorithm is able to
31
distinguish similar and dissimilar waveforms, even for low signal noise ratio and
32
emergent events, which is important for accurate and rapid selection of microseismic
33
events from a large number of records. This method can be applied to other
34
geophysical analyses based on the waveform data.
35 36
INTRODUCTION
37
Generally, the microseismic activity is induced by small magnitude hydraulic
38
fracturing (Maxwell and Urbancic, 2001) and acquired in low signal to noise (SNR)
39
environments. Meanwhile, the number of raw records is very large, especially from
40
real-time monitoring, while the majority of microseismic events seldom occur. Prior
41
to the localization and mechanism analysis of the source, the identification and
42
detection of events become an important challenge. Normally, the existing
43
event-detecting processes can be split into two steps: first, detection of the segment of
44
the records containing the seismic events, and then accurately identifying and picking
45
the target phase. The first stage of detection for a microseismic event is
46
time-consuming, as it requires examination of a large number of raw records; however,
47
all subsequent steps require less resources, as they are based on the previous
48
selections.
49 50
The correlation-type and energy analysis are two principal methods for event
51
detection. The energy analysis is a broadly applicable method due to few assumptions
52
about the data. The most-used approach is defined as the energy ratio of the
53
short-time-average to the long-time-average (STA/LTA) (Earle and Shearer 1994)
54
within moving time windows, which can detect the appearance of seismic events
55
when the ratio exceeds the given thresholds. However, this method is not sensitive to
56
the weak low-magnitude events that have loud background noise, similar to those
57
which occur during hydraulic fracturing. The correlation-type methods, such as
58
cross-correlation within the time domain, detect similar events to a referenced event,
59
known as the master event, by a calculated correlation coefficient to it (Gibbons and
60
Ringdal, 2006; Song et al. 2010). However, the referenced master template waveform
61
is difficult to make when the majority of seismic events occur from unknown and
62
dissimilar sources, and the calculated similarity coefficient is inaccurate without
63
sufficient data, meaning that a long data collection window is required for accuracy.
64 65
Both of the described methods are amplitude-biased, due to the strong dependence on
66
the waveforms and relative amplitudes (Schimmel and Paulssen, 1997). Therefore,
67
results may be ambiguous due to the waveform distortion or weak signals. Some
68
amplitude-normalized cross-correlation (Neidell & Taner 1971) and semblance-based
69
stacks (Taner and Koehler, 1969) are widely used. These methods can remove the
70
disturbing influence of energetic features by some additional pre-processing or
71
post-processing steps, such as normalization, but the waveform coherence may
72
deteriorate. Therefore, several instantaneous phase-based measurements are
73
introduced, which are amplitude-unbiased. The phase-weighted cross-correlations
74
(Schimmel, 1999) and phase-weighted stacks (Schimmel, 1997) are designed for
75
noise attenuation or signal extraction. They are weighted by the coherency of the
76
instantaneous phase, which implicitly contains the waveform information through the
77
Hilbert Transform.
78 79
Normally, when the seismic waveform is resolved into the time-varying frequency
80
components, the time-frequency spectrum can reveal some underlying periodicities
81
and express more abundant characteristics of the seismic waveform. In practice, the
82
useful seismic signals concentrate on some limited frequency band, which is
83
meaningful for noise attenuation and signal extraction. The phase-only correlation
84
(POC) is a simple and robust technique to evaluate the similarity of two images. The
85
technique has been successfully applied to electronics and communication
86
engineering, such as in fingerprint matching (Ito et al., 2004) and subpixel image
87
registration (Takita et al., 2004). The technique has also been expanded into the
88
geophysical exploration field, where POC is used to identify similar waveforms, such
89
as the events among the aftershocks of an earthquake (Moriya, 2011, Zhang et al,
90
2015).
91 92
In this study, we apply the POC for automatic event detection within a single
93
three-component microseismic record. The process workflow is similar to the
94
well-known polarization filtering approach (Jurkevics, 1988), except for identifying
95
similar waveforms with POC based on the transformed time-frequency spectrum. The
96
method first splits the raw records into some contiguous segments, and then it detects
97
the segment containing the event by judging whether the P-wave phase exists. For the
98
improved judging process of the proposed method, the following key steps are
99
necessary: decompose the seismic waveform into the time-varying frequency
100
components with the S-transform, which can retain the absolute phase of each
101
frequency and offer a high resolution, calculate the similarity of the frequency content
102
in the time-frequency domain using the POC function, and identify the P-phase by the
103
combination analysis between any two components. The proposed algorithm, together
104
with the polarization filtering approach using the cross-correlation in the time domain,
105
is tested with synthetic microseismic datasets and real field-recorded datasets.
106 107
METHODOLOGY
108 109
S-Transform
110
The S-transform enables a better resolution of low frequency components and time
111
resolution of high frequency signals by employing frequency-dependent windows,
112
similar to the wavelet transform. Meanwhile, it retains the absolute phase of each
113
frequency (Pinnegar and Mansinha, 2003) for multiresolution analysis and detection
114
of events in time series.
115 116
In normal signal processing, we need to determine the frequency spectrum, , of
117
the time series, , of any component of the raw records, by using the following
118
Fourier transform:
119
= . (1)
120 121
With the additional frequency dependence of the analyzing window, the S-transform
122
of the time series, , can be defined
123
, = − , , ≠ 0, (2)
124
, = 0 = → ! , (3)
!
125 126
where − , is a Gaussian function, which maintains the center at time, τ, of
127
the selected Gaussian window and a standard deviation proportional to ## for the
128
localized spectrum.
129
− , =
||
% √
'(! )'*! !+!
, , > 0, (4)
130 131
where k is the scaling factor, which controls the time-frequency resolution by the
132
number of oscillations in each Gaussian window.
133 134
Aimed to simplify the implementation by computers, equation (3) can also expressed
135
as a multiplication of spectrum:
136
, =
/
+
!1!
2! +! (! 23
/, ≠ 0, (5)
137 138
where / + is the shifted spectrum of the Fourier transform. The Gaussian
139
function in the above integrand equation has a frequency-dependent bandwidth. More
140
specifically, the phase-shifted time signal, , is bandpassed at the center frequency
141
with the deviation 4 =
%||
by the Gaussian window function.
142 143
The time-frequency representation of the S-transform is unique and invertible
144
(Schimmel and Gallart, 2005). The local time-frequency spectra, , , can be
145
transformed back to the time domain, since the Gaussian window for the S-transform
146
must satisfy the condition
147
− , = 1, (6)
148 149
and the time series, u789 t, is easily reconstructed as:
150
;<= = , . (7)
151 152
Band-limited POC
153
Next, we considered the transformed time-frequency spectra , and ,
154
of any two of the components of the raw records as two images, > , > and
155
> , > , respectively. Because the only similarity of the frequency content in the
156
time-frequency domain is needed, the two axes of time, , and frequency, , were
157
equal, familiar with the general 2-D pictures, and rewritten as > and > ,
158
respectively. Conveniently, the size of each of the images was defined as ? × ? ,
159
where the index ranges were > = −A , ⋯ , A and > = −A , ⋯ , A , and
160
therefore, ? = 2A + 1 and ? = 2A + 1 . Then, the 2-D discrete Fourier
161
transform of the images, > , > and > , > , were described as the following
162
equations (Ito et al., 2004):
163 F
F
% G
% G!
= L , , , MNI %I ,%! , (8)
F
F
% G
% G!
= L , , , MN! %I ,%! , (9)
164
D , , , = ∑GIIH FI ∑G!!H F! > , > JKII I JK!!
165
D , , , = ∑GIIH FI ∑G!!H F! > , > JKII I JK!!
166 M
!1 OI
, and JK! =
M
!1 O!
167
where , = −A , ⋯ , A , , = −A , ⋯ , A , JKI =
168
L , , , and L , , , are the amplitude components of the two images,
.
169
respectively, and P , , , and P , , , are the relative phase components.
170
Normally, the cross-spectrum, QR, , , , between D , , , and D , , , can be
171
defined as:
172
RRR% ,% S % ,% S QR, , , = |SI %I ,%! SRRR!%I ,%! | = MN%I ,%! , (10) I
I
!
!
I
!
173 174
where DR , , , represents the complex conjugate of D, , , , and P, , , is
175
the difference of the phases P , , , and P , , , . The POC correlation
176
function is denoted as the 2-D inverse discrete Fourier transform of the
177
cross-spectrum, QR , , , , which can be expressed as
178
T̅ > , > =
KI K!
F! %IGI % G I R ∑F JK! ! ! , (11) GI H FI ∑G! H F! Q , , , JKI
179 180
where T̅ > , > is the value of the POC function, which offers a range from 0 to 1.
181
The peak value of the function can indicate the similarity of the two images and the
182
position provides the translational displacement ∆> , ∆> between the two images:
183
∆> , ∆> = −> , −> . (12)
184 185
Additionally, the function gives a distinct sharp peak of value for two similar images,
186
and it decreases significantly when they are not similar, which demonstrates higher
187
discrimination capabilities than the traditional correlation. Meanwhile, as a phase
188
indication of two images, it works well against the images shifts and brightness
189
changes (Ito et al., 2004).
190
191
In practice, the significant information in the 2-D DFT of the given image is usually
192
clustered in a limited frequency band. Hence, some phase components in the higher or
193
lower frequency domain are meaningless or not reliable. The basic idea of the
194
band-limited POC is to improve the matching performance by eliminating
195
meaningless frequency components before the calculation of the cross-phase
196
spectrum, QR, , , . So, the band-limited POC function can be defined as
197
T̅W XYZ > , > = W
I W!
[! %I GI % G I R ∑[ JW! ! ! , (13) GI H [I ∑G! H [! Q , , , JWI
198 199
ehere k and k are the given inherent frequency bands with the ranges of
200
, = −\ , ⋯ , \ and , = −\ , ⋯ , \ , and 0 ≤ \ ≤ A and 0 ≤ \ ≤ A . In
201
the processing, the parameters \ and \ can be selected with an effective coverage
202
range of the meaningful phase components.
203 204
Combination analysis
205
Theoretically, the three-component seismogram can represent the motion in three
206
orthogonal directions of a ground detector, containing two in the horizontal plane and
207
one in the vertical. The observed seismogram is often represented as the convolution
208
of the source function, the transmission response of subsurface media, or the receiver
209
function. The source refers to a series of vibrations in the Earth due to an earthquake
210
or explosion. The subsurface media are assumed to be heterogeneous, where some
211
reflection and transmission takes place and make wave encounters discontinuous.
212
Some experiments (Morlet et al., 1982) have demonstrated that the above
213
phenomena of wave propagation can considered as some "scattering." The distortion
214
waveforms due to the scattering are frequency-dependent and thus isolated in certain
215
scales (Anant and Dowla, 1997). Therefore, the arrival waves are normally fairly
216
similar, with possible differences in the very high frequencies (very low scales)
217
(Anant and Dowla, 1997). Based on the frequency-dependent nature of the scattering,
218
the arrivals can be viewed in terms of the type of decomposition, and they are only
219
present in some scales, not consistently across all scales.
220 221
By definition, the P-phase is represented by the longitudinal or compressional motion,
222
while the S phase is represented by the transverse motion. The compressional body
223
wave is often linearly polarized. Therefore, it would be more accurate to detect P
224
arrivals based on the linear polarization and the frequency dependence. In our method,
225
the P-wave phase can be identified under the condition of simultaneous high
226
correlation values between any two components within the raw three-component
227
seismograms records.
228 229
EXAMPLES
230
As described above, the proposed workflow is much the same as the well-known
231
polarization filtering approach (Jurkevics, 1988), or other similar methods, except that
232
it identifies the similar waveforms using the POC based on the transformed
233
time-frequency spectrum, replacing the cross-correlation in the time domain or other
234
amplitude-based methods. As shown in Figure 1, the method first splits the raw
235
records into some contiguous segments, and then it detects the segment containing the
236
event by judging whether the P-wave phase exists. For improved judging processes,
237
the examples with similar and dissimilar waveforms were tested as the two
238
possibilities in each loop. The proposed algorithm-based band-limited POC and the
239
1-D cross-correlation in time domain of the raw records were compared mainly for
240
distinguishing capabilities between the similar and dissimilar waveforms. The
241
window length and the threshold were two important parameters in the proposed
242
method. In the proposed method, the lag of any two P-phase arrival times within the
243
three-component geophone was almost small. Meanwhile, the longer window may
244
contain the target data of the P-phase, together with other phases, such as the coda
245
wave, which may affect the accuracy of judgment. Therefore, in order to improve
246
the accuracy and robustness, the window length of the time series is recommended to
247
be between one-half and one waveform length. The window moves from left to right
248
with an overlapped length of window in each loop. For some batch processing, such
249
as real-time monitoring, the threshold can be obtained by pre-processing experiments,
250
which mainly depend on the performance of the hydraulic fracturing, and also by
251
measuring the background noise in the acquired environment.
252 253
(The position of Figure 1)
254 255
Synthetic microseismic datasets
256
A two-layer model was used for the synthetic microseismic datasets modeling. The
257
P-wave velocity of the model was set to 1500, 3000 m/s, the S-wave velocities to
258
1500, and 3000 m/s , and the density to 11, 2.5 kg/mh . Approximately 10-level
259
three-component geophones were placed in the vertical observation well, which were
260
spaced at intervals of 25 m within the array. The microseismic source was located
261
atx, y, z = 500, 300, 4000, with a dominant frequency of 400 Hz. The data were
262
recorded at 280 ms with a sampling rate of 4 samples per millisecond.
263 264
(The position of Figure 2)
265 266
The recorded three-component seismograms within an array are displayed in Figure
267
2(a). In practice, the microseismic events within real field records were weak and had
268
superimposing noise from background noise, such as vibrations from machinery. To
269
better simulate real data acquisition, Gaussian noise (approximately 50% of the
270
maximum signal amplitude) was added to the synthetic records. Figure 2(b)
271
demonstrates that the noise affected the microseismic records, which would present a
272
challenge to detect the events reliably using some automatic processes. The window
273
length of any selected segments was set to 32 ms.
274 275
(The position of Figure 3-5)
276 277
The first example demonstrates the segment from the black box in Figure 2(b), which
278
contains the major events from 60 to 92 ms. Figure 3(b)-(d) demonstrates the
279
time-varying spectral representation of the components x, y, and z, respectively. The
280
transformed spectral representation captures the features of seismic phase arrivals,
281
especially the predominant frequency distribution. Then, the frequency contents from
282
300Hz to 600Hz were selected as the significant information for the following steps,
283
which means that both the very low and high frequency noise were the most excluded.
284
Figure 4(a)-(c) demonstrates the POC coefficients of the coupled components x-y, x-z,
285
and y-z, respectively. The high values were mostly distributed at the center of the
286
images, which had a high correlation value with a frequency lag of 0 and a time lag of
287
0. The peak of amplitude was very sharp, which makes it easy to detect with given
288
threshold. Then, we calculated the similarity using the cross-correlation function in
289
time domain with the same window length as the records. Figure 5 demonstrates the
290
POC coefficient (red curve) at the frequency lag 0 and the cross-correlation
291
coefficient in time domain (blue curve) corresponding to Figure 4(a)-(c). Both
292
methods can indicate the high correlation values at the time lag 0. However, the
293
cross-correlation coefficient in the time domain demonstrated a lower resolution
294
compared to the POC with the short window length of the selected time series because
295
the calculated similarity coefficient of the cross-correlation in the time domain would
296
be more accurate with sufficient data, meaning a longer window size.
297
(The position of Figure 6-8)
298 299
The second example demonstrates the segment from 0 to 32 ms, from the red box in
300
Figure 2(b), which mainly contains random noise without appreciable signals. Figure
301
6(b)-(d) are the time-varying spectral representation of the components x, y, and z,
302
respectively. Different from the previous example, the transformed spectral
303
representations have a large difference in the frequency distribution. Then, the same
304
parameters were set for the calculation of the POC. Figure 7(a)-(c) demonstrates the
305
POC coefficients of the coupled components x-y, x-z, and y-z. The highest
306
correlations were not at the center of the image with the frequency lag 0 and the time
307
lag 0, and the values were reduced considerably. The values of the dissimilar images
308
were relatively small compared the peak value of the similar images, which was easy
309
to distinguish. In the next comparison with the cross-correlation coefficient in the
310
time domain, both methods were stable and demonstrated small coefficient values,
311
indicating the dissimilar waveforms of the pair of segments. However, the normalized
312
cross-correlation coefficient in the time domain distributes in a relatively wider range.
313 314
Real field-recorded datasets
315
Figure 9(a) demonstrates the continuous microseismic records for approximately 3000
316
ms with a sampling rate of 4 samples per millisecond. The 4th level of records
317
demonstrated in Figure 9(b) was used to show an actual process based the proposed
318
method. And the window length was also set to 32 ms . A large-magnitude
319
microseismic event (the black arrow B) was noticeable, which could be easily
320
selected by certain correlation-based methods. However, some potential small
321
magnitude events in the background noise are difficult to detect correctly.
322
323
(The position of Figure 9-10)
324 325
As a result of above processing algorithm, three segments with potential events were
326
detected in Figure 10(a) with red traces and magnified to more detail in Figure
327
10(b)-(d), separately. The segment from 2240 to 2272 ms in Figure 10(c) contains the
328
large-magnitude microseismic event mentioned before. But, the other two segments
329
containing potential events from 1400 to 1432 ms in Figure 10(b) and from 2930 to
330
2962 ms in Figure 10(d) were very weak, which still needs investigation. The
331
following example demonstrated the segment from 2240 to 2272 ms in Figure 10(c)
332
in more detail. Figure 11(a)-(c) were the time-varying spectral representations of
333
components x, y, and z, respectively. Figure 12(a)-(c) demonstrates the POC
334
coefficients of the coupled component x-y, x-z, and y-z. Figure 13 demonstrates the
335
normalized cross-correlation coefficient in the time domain (blue curve) and the POC
336
function coefficient (red curve) at the frequency lag 0. The cross-correlation
337
coefficient had multiple blunt peaks, which were close to the max value of the POC
338
coefficient. However, it is difficult to indicate the relative difference of the arrival
339
times of the compared waveforms in the selected segments with the constant
340
threshold.
341 342
(The position of Figure 11-13)
343 344
DISCUSSIONS
345
This study introduced an effective method for the identification and detection of the
346
microseismic events by judging whether the P-wave phase exists in a local segment
347
from single three-component microseismic records. The improved judging process
348
contained the following key steps: first, the waveform time series was correctly
349
transformed into time-varying spectral representations using the S-transform. Next,
350
the similarity of the frequency content in the time-frequency domain was calculated
351
using the POC function, and the P-phase was identified by the combination analysis
352
between any two components. The proposed algorithm was compared to similar
353
approaches using the cross-correlation in time domain between any two components,
354
and it was capable of distinguishing between the similar and dissimilar waveforms.
355 356
The proposed workflow is similar to the well-known polarization filtering approach,
357
except that it identifies the similar waveforms using the POC based on the
358
transformed time-frequency spectrum, replacing the cross-correlation in the time
359
domain, or other amplitude-based methods. For the proposed workflow, the lag of any
360
two P-phase arrival times within the three-component geophone is almost small.
361
Meanwhile, the longer window may contain the target data of the P-phase, which,
362
together with other phases, such as coda wave, may affect the accuracy of judgment.
363
So, the window length of the time series is recommended to be between one-half and
364
one waveform length for improved accuracy and robustness. However, the calculated
365
similarity coefficient by the cross-correlation coefficient in the time domain would be
366
more accurate with sufficient data, meaning a longer window size. Therefore, the
367
short window size or short data may result in lower resolution coefficients with the
368
cross-correlation in the time domain, where it is difficult to indicate the relative lag of
369
the arrival times of the compared waveforms in the selected segments. Therefore, the
370
POC-based workflow with high resolution lag cross-correlation coefficient is
371
a great choice for the proposed workflow, which is meaningful for accurate
372
identification and detection of the microseismic events. In addition, distinguishing
373
between the similar and dissimilar waveforms with high resolution was easily
374
achieved with some constant thresholds, which is ideal for automating a process.
375 376
However, processing with POC has a low probability of providing false indications,
377
mainly because the calculated coefficient in the time-frequency domain is sensitive
378
for some unreasonable frequency content. In consideration of the computation
379
efficiency and stability, it is suggested that the practical process should recombine and
380
combine the bright side of the traditional correlation-type and energy analysis
381
methods and the POC method, especially for rapid field analysis and hydro-fracture
382
monitoring.
383 384
CONCLUSIONS
385
This work presents an effective method for the detection of the microseismic events
386
by judging if the P-wave phase exists in a local segment using the POC function. The
387
results of the synthetic data and the real hydraulic-fracture case show that the
388
proposed method is sufficiently stable to distinguish the similar and dissimilar
389
waveforms, even for low SNR, which means that it can reliably handle the two
390
possibilities during the judging algorithm. It is important to identify and detect the
391
microseismic events within a large number of records accurately and rapidly, and this
392
method can be applied to other geophysical analyses based on the waveform data.
393 394
ACKNOWLEDGMENTS
395
This research was funded by the National Natural Science Foundation of China (Grant
396
no. 41230317) and the National Basic Research Program of China (Grant no.
397
2015CB258500).
398
399
REFERENCES
400
Anant K S, Dowla F U. Wavelet Transform Methods for Phase Identification in
401
Three-Component Seismograms [J]. Bulletin of the Seismological Society of America,
402
1997, 87(6):1598-1612.
403 404
Earle P S, Shearer P M. Characterization of global seismograms using an
405
automatic-picking algorithm [J]. Bulletin of the Seismological Society of America,
406
1994, 84(2): 366-376.
407 408
Gibbons S J, Ringdal F. The detection of low magnitude seismic events using
409
array-based waveform correlation [J]. Geophysical Journal International, 2006, 165(1):
410
149-166.
411 412
Ito K, Nakajima H, Kobayashi K, et al. A Fingerprint Matching Algorithm Using
413
Phase-Only Correlation [J]. Ieice Transactions, 2004: 682-691.
414 415
Jurkevics A. POLARIZATION ANALYSIS OF THREE-COMPONENT ARRAY
416
DATA[J]. Bulletin of the Seismological Society of America, 1988, 78(5): 1725-1743
417 418
Maxwell S C, Urbancic T I. The role of passive microseismic monitoring in the
419
instrumented oil field [J]. Geophysics, 2012, 20(6): 636-639.
420
Morlet J, Arens G, Fourgeau E, et al. Wave propagation and sampling theory—Part I:
421
Complex signal and scattering in multilayered media [J]. Geophysics, 2012, 47,
422
203-221.
423 424
Moriya H. Phase-only correlation of time-varying spectral representations of
425
microseismic data for identification of similar seismic events [J]. Geophysics, 2011,
426
76(6).
427 428
Neidell N S, Taner M T. Semblance and other coherency measures for multichannel
429
data [J]. Geophysics, 1971, 36(3): 482-497.
430 431
Pinnegar C R, Mansinha L. The S-transform with windows of arbitrary and varying
432
shape [J]. Geophysics, 2003, 68(1):381.
433 434
Schimmel M, Gallart J. The inverse S-transform in filters with time-frequency
435
localization [J]. IEEE Transactions on Signal Processing, 2005, 53(11): 4417-4422.
436 437
Song F, Kuleli H S, Toksoz M N, et al. An improved method for
438
hydrofracture-induced microseismic event detection and phase picking [J].
439
Geophysics, 2010, 75, A47–A52.
440 441
Stockwell R G, Mansinha L, Lowe R P, et al. Localization of the complex spectrum:
442
the S transform [J]. IEEE Transactions on Signal Processing, 1996, 44(4): 998-1001.
443 444
Takita K, Aoki T, Sasaki Y, et al. High-accuracy subpixel image registration based on
445
phase-only correlation [J]. Ieice Transactions, 2003: 1925-1934.
446 447
Zhang L, Wang Y, Chang X, et al. Wigner distribution based phase-only correlation[J].
448
Journal of Applied Geophysics, 2015: 277-282.
449 450
Figure 1. The workflow of the proposed algorithm of identification and detection of
451
the microseismic events
452
453 454
(a)
455 456
(b)
457
Figure 2. (a) Synthetic microseismograms produced by 3D elastic wave equation
458
modeling. Every triplet of x (red), y (blue), and z (green) components is plotted for
459
the 3-component seismograms of each geophone. (b) Gaussian noise was added with
460
a given SNR of approximately 3. The segments in the red and black boxes are used
461
for the following tests of identification capability.
462
463 464
(a)
465 466
(b)
467 468
(c)
469
(d)
470
Figure 3. (a) The selected segment contains the major events with three components
471
from the black box in Figure 2(b); (b)-(d) time-varying spectral representation of
472
components x, y, and z, respectively.
473
474 475
(a)
476 477
(b)
478 479
(c)
480
Figure 4. (a) - (c) The POC function of the couple components x-y, x-z, and y-z,
481
respectively. The sharp peaks of amplitude are significant, which are easy to be
482
detected with some given threshold.
483 484
Figure 5. The comparison of the POC function coefficient (red curve) at the frequency
485
lag 0 and the normalized cross-correlation coefficient in time domain (blue curve)
486
corresponding to Figure 4(a)-(c). Both methods can indicate the max correlation
487
coefficients at the same point. However, the cross-correlation coefficient in the time
488
domain is less resolved.
489
490 491
(a)
492 493
(b)
494 495
(c)
496 497
(d)
498
Figure 6. (a) The selected segment contains the major events with three components,
499
from the red box in Figure 2(b); (b)-(d) time-varying spectral representation of
500
components x, y, and z, respectively.
501
502 503
(a)
504 505
(b)
506 507
(c)
508
Figure 7. (a) - (c) The phase-only correlation function of the couple components x-y,
509
x-z, and y-z, respectively. The high values show a scattered distribution, which is
510
treated as a weak correlation.
511 512
Figure 8. The comparison of the phase-only correlation function coefficient (red curve)
513
at the frequency lag 0 and the normalized cross-correlation coefficient in the time
514
domain (blue curve) corresponding to the Figure 7(a)-(c). Both methods show stable
515
results with small values in a tight distribution.
516
517 518
(a)
519 520
(b)
521
Figure 9. (a) The 7-level seismograms of the downhole geophone array with some
522
triplets of x (red), y (blue), and z (green) components. (b) The 4th level of the records
523
in Figure 9 is selected for the following testing. The red box demonstrates a split
524
segment, which moves from left to right with an overlapped half length.
525
A
526 527
(a)
528 529
(b)
530 531
(c)
B
C
532 533
(d)
534
Figure 10. (a) Three detected segments with potential events extracted in Figure 10(a)
535
with red traces; (b)-(d) the extracted and magnified traces in more details. The
536
segments containing the large-magnitude microseismic event (arrow B) and some
537
potential weak events (arrow A and C) are detected, though they will be identified by
538
some additional information.
539
540 541
(a)
542 543
(b)
544 545
(c)
546
Figure 11. The analysis of the selected segment from 2240 to 2272 ms from the 4th
547
level in Figure 10 (b); (a)-(c) time-varying spectral representation of components x, y,
548
and z, respectively.
549
550 551
(a)
552 553
(b)
554 555
(c)
556
Figure 12. (a) - (c) The phase-only correlation function of the couple components x-y,
557
x-z, and y-z, respectively, which indicate the maximum correlation coefficients with
558
the relative peak values.
559 560
Figure 13. The comparison of the phase-only correlation function coefficient (red
561
curve) at the frequency lag 0 and the normalized cross-correlation coefficient in the
562
time domain (blue curve) corresponding to Figure 11(a)-(c). The cross-correlation
563
coefficient in the time domain has multiple blunt peaks, which make it difficult to
564
indicate the relative difference of the arrival times of the compared waveforms in the
565
selected segments.