See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/266633928

Automated Microseismic Event Location Using Finite Difference Traveltime

Calculation and Enhanced Waveform Stacking

Article · June 2013

DOI: 10.3997/2214-4609.20130398

CITATIONS READS

3 154

3 authors, including:

Junwei Huang Juan reyes-montes

Scotiabank University of Liverpool
32 PUBLICATIONS   103 CITATIONS    43 PUBLICATIONS   449 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Migration based microseismic event location View project

Passive traveltime tomography View project

All content following this page was uploaded by Junwei Huang on 04 April 2015.

The user has requested enhancement of the downloaded file.

Automated microseismic event location using finite difference traveltime calculation and
enhanced waveform stacking

Junwei Huang1*, Juan M. Reyes-Montes1and R. Paul Young2

1
Applied Seismology Consultants Ltd, 5 Claremont Buildings, Claremont Buildings, Shrewsbury, SY1 1RJ, UK
2
University of Toronto, 170 College Street, Toronto, Ontario M5S 3E3, Canada
*junwei@appliedseismology.com
Abstract

The Fast Sweeping Method (FSM) is a finite difference algorithm providing significant computational
efficiency and modelling capability in calculating the first arrivals of seismic P- or S-waves. Based on
the calculated travel-timetable, we stack waveform amplitude for every possible source location and
origin time followed by semblance weighting. With three-component (3-C) data, the stacked image
can be further improved by matching the waveform polarization with the modelled ray vector from
FSM and thus reduce the azimuth ambiguity in the location of microseismic (MS) events. We apply
the algorithm on synthetic surface and borehole data and show that the combination of FSM and
enhanced waveform stacking can efficiently locate MS events with high spatial resolution even when
realistic noise levels are added. Finally we test the algorithm on a set of field 3-C data from the
stimulation of a sandstone reservoir monitored using two borehole arrays. The method allowed the
location of twice the number of event compared to locations using each borehole individually, with
both sets defining a consistent fracture structure.

75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013

London, UK, 10-13 June 2013
Introduction

Passive microseismic (MS) monitoring is an established technology for the monitoring of oil and gas
reservoir treatment, engineered geothermal systems, geological storage and mining operations.
Monitoring is performed using surface or downhole geophone arrays. Downhole monitoring typically
deploys 10s of instruments close to the stimulated reservoir providing high quality signals. Surface
monitoring on the other hand involves ~100s of geophones with the given aperture and fold
requirements to compensate for low signal-to-noise ratio. Techniques of MS event location for both
monitoring techniques generally rely on automated or semi-automated picking of seismic arrivals.
However, automated arrival time picking becomes unreliable for low quality signals and unrealistic
for significant amount of real-time streaming data. This problem is most severe for long term surface
monitoring. The main challenge comes from the data handling and processing for typical size of data
sets from 100s of instruments. In this study, we present an automated MS event locating method
based on waveform stacking similar to the source-scanning algorithm (Kao and Shan, 2004). Our
method uses the Fast Sweeping Method (FSM) to calculate first arrivals through complex reservoir
models. The waveform is stacked following the calculated travel time of P-wave for each possible MS
location and initial time. When S-wave arrivals dominate, the S-wave travel-timetable can also be
used on its own or in combination with the P-wave travel-timetable.

The geometry of the monitoring array is usually constrained by the number of available monitoring
boreholes or the accessibility at the surface. Due to the limited aperture, it is likely that multiple
locations can match the travel time equally well and thus introduce ambiguity in its inversion for
location. This problem is more severe for monitoring projects using a single borehole array.
Conventional borehole monitoring utilizes polarization analysis to reduce location uncertainty (Pettitt
and Young, 2007), therefore we include the waveform polarization analysis in the stacking algorithm
and use the residual between the modelled and the measured polarization to suppress amplitude with
unmatched polarization and thus reduce the spatial ambiguity.

In this paper, we will first introduce the background of FSM and the enhanced stacking algorithm. We
then demonstrate the robustness and efficiency of this algorithm using synthetic surface and borehole
data. We finally test our algorithm on field data from the stimulation of a sandstone oil reservoir
monitored using geophone strings installed down two monitoring boreholes.

Fast Sweeping Method: a showcase

Conventionally, travel time is calculated by ray tracing, which can be problematic due to shadow
zones and multi-pathing problems in complex reservoirs and is computationally expensive when
many source–receiver pairs are involved. FSM is a finite difference method directly solving the
Eikonal equation (e.g., Chapman, 2004). This grid-based method presents numerous advantages over
conventional ray tracing. For example, it assigns a travel time to each grid point and thus avoids
shadow zone problems. In addition, unlike shooting or bending methods in ray tracing, it is
formulated as a forward problem, which is highly stable for first arrival travel time calculations.

Generally FSM advances wave fronts with iterative sweeping and updates travel time at each grid
point monotonically to ensure the causality (Zhao, 2004). The implementation of FSM is thus
significantly simplified and its computational efficiency is superior to other grid-based methods such
as the fast marching method (Sethian, 1996) even when velocity changes rapidly (Bak et al., 2010).
We demonstrate the modelling power of FSM by simulating a point source in a complex velocity
model as shown in Figure 1. The velocity of the model ranges from 4000 m/s to 0 m/s and the zero
velocity regions considered as barriers have no travel time defined. The model contains ~104 grid
points and the total computation time is less than 5 seconds on a single CPU core. Seismic rays can be
traced following the gradient of the travel time field and provide the modelled polarization as well as
travelling distance for both P- and S-wave to any geophone locations (Figure 1b and 1c).

75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013

London, UK, 10-13 June 2013
Figure 1 (a) The velocity model Figure 2 the three-dimensional view and the plane view
constructed from yield zone modelling of of the velocity model and the numerical monitor survey
a developing underground caving mine geometry for an example reservoir monitoring using
(Reyes-Montes et al., 2012), (b) the surface geophones (a) and borehole geophones (b). The
seismic rays traced based on the gradient colour represents the P-wave velocity, and the dots on
of travel time and (c) the contour of the the surface in (a) and the vertical blue cones in (b)
travel time rendered as isosurfaces. The denote the geophones. The red box marks the monitoring
full computation time is less than 5 volume and the red dot is the true location of the
seconds using a single CPU. modelled MS event.

Enhanced stacking of synthetic data

Figure 2 shows the numerical model of MS monitoring using surface geophones (Figure 2a) and
borehole geophones (Figure 2b). Both monitoring surveys are modified from field monitoring jobs. In
total 14 receiver lines each with 79 three-component geophones are deployed along the surface
covering a 4 km by 4 km surface area. The microseismic double-couple source is located at 3500 m
depth beneath the surface with 10 Hz primary frequency. Uncorrelated white noise was later added to
the synthetic seismogram (Figure 3b) to evaluate the robustness of the stacking algorithm against
realistic monitoring conditions (Figure 3d and 3e). The numerical borehole monitoring survey is
conducted on a layered 3D model based on field well logs. Unlike the surface monitoring project, it
focuses on a much smaller scale volume and the primary frequency of the MS event is ~150 Hz. The
modelled MS event is located at the centre of the model (red dot) with a double-couple source
mechanism and recorded by geophones from all boreholes (e.g., Figure 4a).

Although four monitoring boreholes are available in the simulation, only one borehole is used to show
the benefit of including polarization analysis in the stacking algorithm to mimic a typical monitoring
setup. Following the P-wave travel-timetable, we stacked 21 horizontal x-component traces from one
borehole (Figure 4b) and the residual between the modelled ray vector and waveform polarization is
used to suppress the amplitude with unmatched polarization. As shown in Figure 4b, among all points
satisfying the travel time table equally well, the points with a better match of the waveform
polarization stand out (marked by the three black dots: the maximum and minimum amplitude and the
middle point) and provide a good estimate of the true MS event location (red dot).

Enhanced stacking of field data

We applied the enhanced stacking algorithm to a set of field data recorded on two borehole arrays
(Figure 5) during the stimulation of a sandstone reservoir. The velocity model used in the forward

75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013

London, UK, 10-13 June 2013
modelling is based on borehole logs and calibrated using ball seat events. The axis of x, y, and z point
towards North, East and Downward directions, respectively. As shown in Figure 5a, the layered
model has a strike angle of 118 o and a dipping angle of 1.8 o, deepening towards the south-west. The
stacking algorithm located 45 events for one stage (compared to 21 events using each borehole
separately), with locations defining a fracture network developing towards the south-west. Plotting the
theoretical travel time on top of the waveform, we find it is easy to verify that the estimated location
is consistent with the observed S-wave arrivals (see Figure 5e).

Figure 3 (a) Synthetic seismogram of horizontal x-component without random noise. (b) Synthetic
seismogram of horizontal x-component with random noise. The noise amplitude is designed to be
slightly higher than S-wave amplitude thus both P- and S-wave arrivals are hardly visible. (c) Stacked
image using noise-free horizontal x-component as seen in (a) using P-wave travel time only. The three
black dots mark the grid points with maximum and minimum amplitudes and the middle point. The
true location, marked by a red dot, is overlapped by the middle point. (d) Stacked image using noise
contaminated horizontal x-component as seen in (b) using P-wave travel time only. The high noise
level significantly shifts the estimated location (black dots) away from the true location. (e) Stacked
image using noise contaminated horizontal x-component as seen in (a) using S-wave travel time only.
Due to the higher amplitude of the S-wave, the stacking is adequate to recover the source location.

Figure 4 (a) Synthetic seismogram of horizontal x-component from all monitoring wells. Synthetic
traces from well a (marked by dashed box) are stacked following P-wave travel time table. (b) The
plan view of the stacked image using 21 traces from a single borehole indicated by the blue dot at the
lower left corner. The estimated location range is denoted by the three black dots. It is important to
point out that all the grid points along the curve match the travel time equally well. The polarization
constrain from the ray vectors suppresses the waveform with unmatched polarization. The side view
from northeast direction marked by the arrow is shown in (c), where the blue cones mark the location
of geophones.

75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013

London, UK, 10-13 June 2013
 Conclusions

In this paper, we present an algorithm for the location of MS events based on the fast sweeping
method and waveform stacking. The fast sweeping method proved to be computational efficient in
calculating the travel times through complex velocity models. Seismic ray vectors can be calculated
using the travel time field and provide further constraint for the stacking algorithm using the modelled
polarizations. The algorithm was tested on synthetic surface and borehole data with modelled low
signal-to-noise conditions, providing accurate estimated locations. Finally, the test on a set of field
data from two monitoring boreholes showed that the stacking algorithm efficiently located more MS
events enabling a better statistical analysis of the fracturing process.

Figure 5 (a) The velocity model and

the located 45 events for one stage
from a field reservoir stimulation job.
The white triangle marks the location
of the initiation point. Axes x, y, and
z point towards North, East and
Downward directions, respectively.
The dashed line indicates the path of
the horizontal treatment well.
Located event are coloured scaled to
their occurrence time relative to
stage initiation. Two side views and
the plan view are shown in
(b),(c),and (d). (e) Sample waveform
of one event. The first 10 traces are
from well 1 and the rest from well 2.
The recorded S-wave from both
boreholes is much stronger than P-
wave. The stacking algorithm thus
uses S-wave travel time only. The
first red dashed line marks the
theoretical arrival time of S-wave
without time delay and the second
one is delayed by 113 ms, which is
the occurrence time of this event
relative to the zero time of the
triggered data.

References

Bak, S., McLaughlin, J. & Renzi, D., 2010. Some improvements for the fast sweeping method, SIAM
J. Sci. Comput., 32(5), 2853–2874
Chapman, C., 2004. Fundamentals of Seismic Wave Propagation. Cambridge University Press,
Cambridge, 172pp.
Kao, H. & Shan, S.-J., 2004, The Source-Scanning Algorithm: mapping the distribution of seismic
sources in time and space, Geophys. J. Int., Volume 157, Issue 2, pp. 589-594.
Pettitt, W. S., & R. P. Young, 2007, InSite seismic processor: User operations manual version 2.14:
Applied Seismology Consultants Ltd.
Reyes-Montes, J.M., Sainsbury, B., Andrews, J.R. & Young, R.P., 2012. Application of cave-scale
rock degradation models in the imaging of the seismogenic zone, Proc. 6th Int. Conf. and Ex. on Mass
Mining. Sudbury, Ontario, Canada.
Sethian, J.A., 1996. A fast marching level set method for monotonically advancing fronts, Proc. Natl.
Acad. Sci., 93, 1591–1595.
Zhao, H-K., 2004. A fast sweeping method for Eikonal equations, Math. Comput., 74, 603–627.

75th EAGE Conference & Exhibition incorporating SPE EUROPEC 2013

London, UK, 10-13 June 2013

View publication stats