Seismic event relocation techniques

Intro discussion

Seismic event relocation overview

In the study of mining induced seismicity, the accuracy of event location and magnitudes has improved through the technological advancement of hardware and processing algorithms. Seismic events are typically located using P and S wave arrival times at different receivers paired to a velocity model. One of the most influential developments in micro seismic data processing is event relocation. Research in the field of crustal seismology has developed a wide array of techniques to increase the accuracy of event location. These techniques fall into the broad category of “event relocation” or “relative location methods”. Adding event relocation to a micro seismic data processing flow aids interpretation by concentrating previously diffuse event clouds by reducing positional uncertainty. This increased confidence in event location can add context to the geometry of the source mechanism; adding value to the seismic data set.

Initial (absolute) event location accuracy is generally hampered by the following

• Errors in P and S wave arrival time picks either from auto-picking or human error
• Limited understanding or over-simplification of regional velocity model
• Limited spatial coverage of geophones in mining applications

The techniques outlined below are deemed “relocation” techniques because they are done after an initial (absolute) location solution has been calculated. These types of process can also be considered “post-processing” or “reprocessing” in the seismic data processing workflow. The case examples presented are related specifically to mining induced seismicity, though the development and use of the techniques transcend the field, with applications in hydraulic fracture monitoring and crustal seismology.

When is relocation completed?

In the seismic data processing workflow, event relocation can be performed immediately after an initial (absolute) solution for an event’s location has been calculated. Often a relocation algorithm will be applied to heritage data as a relatively inexpensive way to add value to pre-existing datasets. Due to the advance in computational power in seismic monitoring solutions relocation algorithms can be completed “on the fly” as part of the standard event location processing flow to reduce uncertainty in calculated event locations.

How are is it done?

Seismic event relocation algorithms describe a style of process, with many variations on a handful of themes. The most popular and arguably most effective theme being the Double-difference technique which differs from other relocation algorithms due to its assumptions and simultaneous relocation of large numbers of events over large distances. Other themes include “Master Event Relocation”, where events are moved relative to single “master” event which is computationally straight forward, but propagates location errors from the initial placement of the master event through to the relocated “relative events”. Alternatively, there are “Simultaneous Relocation” approaches. First implemented by Got et al, they determined cross-correlation time delays for all possible event pairs in question and combined them into a system of linear equation which was solved via least squares approximation and converted to positional corrections. For the remainder of this article reference is made to the three aforementioned classifications; Double-difference, Master Event, and Simultaneous Event relocation algorithms to summarize the current state of event relocation techniques and enable direct comparison between them.

Double-difference (DD) method

The double-difference approach, first introduced by Waldhauser in 2000 has been applied successfully as an event relocation algorithm option in mining induced seismicity applications.

Double-difference (DD) workflow

The following procedure for executing the double difference algorithm is adapted from (Waldhauser, Ellsworth 2000).

First, both P and S-wave differential travel times are derived from cross-spectral (cross correlation) methods with travel-time differences formed from catalog data and minimize residual differences (or double differences) for pairs of earthquakes by adjusting the vector difference between their hypocenters. From this calculated “double difference” we are able to determine interevent distances between correlated events that form a single multiplet (group of similar events) to the accuracy of the cross-correlation data while simultaneously determining the relative locations of other multiplets and uncorrelated events to the accuracy of the absolute travel-time data, without the use of station corrections.

The arrival time, T, for an earthquake, i, to a seismic station, k, is expressed using ray theory as a path integral along the ray,

$T^{i}=\tau ^{i}+\int _{i}^{k}uds\qquad (1)$

where τ is the origin time of event i, u is the slowness field, and ds is an element of path length. Due to the nonlinear relationship between travel time and event location, a truncated Taylor series expansion (Geiger, 1910) is generally used to linearize equation (1). The resulting problem then is one in which the travel-time residuals, r, for an event i are linearly related to perturbations, Δm, to the four current hypocentral parameters for each observation k:

${\frac {\partial t_{k}^{i}}{\partial \mathbf {m} }}\Delta \mathbf {m} ^{i}=r_{k}^{i}\qquad (2)$

where rik = (tobs - tcal)ik, tobs and tcal are the observed and theoretical travel time, respectively, and Δmi = (Δxi, Δyi, Δzi, Δ τ i). Equation (2) is appropriate for use with measured arrival times. However, cross-correlation methods measure travel-time differences between events, (tki - tkj)obs, and as a consequence, equation (2) can not be used directly. Frechet (1985) obtained an equation for the relative hypocentral parameters between two events i and j by taking the difference between equation (2) for a pair of events,

${\frac {\partial t_{k}^{ij}}{\partial \mathbf {m} }}\Delta \mathbf {m} ^{i}j=dr_{k}^{ij}\qquad (3)$

where Δmij_(Δdxij, Δdyij, Δdzij, Δdsij) is the change in the relative hypo central parameters between the two events, and the partial derivatives of t with respect to m are the components of the slowness vector of the ray connecting the source and receiver measured at the source (e.g., Aki and Richards, 1980). Note that in equation (3) the source is actually the centroid of the two hypocenters, assuming a constant slowness vector for the two events. drikj in equation (3)m L. Ellsworth is the residual between observed and calculated differential travel time between the two events defined as

$dr_{k}^{ij}=(t_{k}^{i}-t_{k}^{j})^{obs}-(t_{k}^{i}-t_{k}^{j})^{cal}\qquad (4)$

We define equation (4) as the double-difference. Note that equation (4) may use either phases with measured arrival times where the observables are absolute travel times, t, or cross-correlation relative travel-time differences. The assumption of a constant slowness vector is valid for events that are sufficiently close together, but breaks down in the case where the events are farther apart. A generally valid equation for the change in hypo central distance between two events i and j is obtained by taking the difference between equation (2) and using the appropriate slowness vector and origin time term for each event

${\frac {\partial t_{k}^{i}}{\partial \mathbf {m} }}\Delta \mathbf {m} ^{i}-{\frac {\partial t_{k}^{i}}{\partial \mathbf {m} }}\Delta \mathbf {m} ^{j}=dr_{k}^{ij}\qquad (5)$

or written out in full

${\frac {\partial t_{k}^{i}}{\partial {x^{i}}}}\Delta {x}^{i}+{\frac {\partial t_{k}^{i}}{\partial {y}}}\Delta {y}^{i}+{\frac {\partial t_{k}^{i}}{\partial {z^{i}}}}\Delta {z}^{i}+\Delta \tau ^{i}-{\frac {\partial t_{k}^{j}}{\partial {x^{i}}}}\Delta {j}^{i}-{\frac {\partial t_{k}^{j}}{\partial {y}}}\Delta {y}^{j}-{\frac {\partial t_{k}^{j}}{\partial {z^{j}}}}\Delta {z}^{i}-\Delta \tau ^{j}=dr_{k}^{ij}\qquad (6)$

The partial derivatives of the travel times, t, for events i and j, with respect to their locations (x, y, z) and origin times (s), respectively, are calculated for the current hypocenters and the location of the station where the kth phase was recorded. Δx, Δy, Δz, and Δs are the changes required in the hypo central parameters to make the model better fit the data. We combine equation (6) from all hypocentral pairs for a station, and for all stations to form a system of linear equations of the form

$\mathbf {WGm} =\mathbf {Wd} \qquad (7)$

where G defines a matrix of size M _ 4N (M, number of double-difference observations; N, number of events) containing the partial derivatives, d is the data vector containing the double-differences (4), m is a vector of length 4N, [Δx, Δy, Δz, ΔT]T, containing the changes in hypocentral parameters we wish to determine, and W is a diagonal matrix to weight each equation. We may constrain the mean shift of all earthquakes during relocation to zero by extending (7) by four equations so that

$\sum _{i=0}^{N}i\qquad (8)$

for each coordinate direction and origin time, respectively. Note that this is a crude way to apply a constraint, but appropriate for a solution constructed by conjugate gradients (see Lawson and Hanson [1974] for more exact solutions of constrained least squares). As shown later, the double difference algorithm is also sensitive to errors in the absolute location of a cluster. Thus, equation (8) is usually down weighted during inversion to allow the cluster centroid to move slightly and correct for possible errors in the initial absolute locations.