AC

A

T
B

____ Figure 8.—Seismic event locations using intersections of relative-arrival-time isochrons. A, B, and C are re-
‘celvers; AB, AC, and BC are receiver pairs. A, Well-constrained location with small equidimensional (schematic)
SVD error ellipsoid; B, poorly constrained location with elongate SVD error ellipsoid axis parallel to direction of

least constraint.

are to be avoided. These provide little constraint per-
pendicular to the axis and/or plane of receivers because
the relative arrival-time isochron surfaces largely intersect
at very shallow angles for all sources placed off axis or off
plane. Slight measurement errors can result in very large
displacements of the solution along the long axis of the
intersection volume (or SVD ellipsoid; fig. 8B).”

These array geometries lead to difficulties in three-dimensional loca-
tion solutions even without errors. For example, two-dimensional (pla-
nar) arrays result in a mirroring of the location solution about the sym-
metry plane and linear arrays lead to an infinite number of solutions
along a co-axial circle. Another solution problem arises when stations
are distributed on the surface of a right-circular cone with the source at

its apex (11).

The effects of linear or planar receiver geometries can
be present in an otherwise three-dimensional array when
only a subset of the receivers have usable arrivals for the
location process. This is especially true for array ge-
ometries that are initially biased toward linear or planar
receiver distributions, Only the largest events, represent-
ing a rather small fraction of all detected events, provide
usable signals on all receivers, and thus actually use the
full three-dimensional extent of an array.

Since many mining geometries provide access points
for receiver installation on a plane (e.g. coal and vein
deposits), it is common for underground arrays to possess
groups of stations that approximate planar distributions.


In practice, the magnitude of the gradient vector for
either gradient basis can be limited to some specified max-
imum size. Thus, if |dij| > max, one finds a new dx
such that |di’ | = max and diy’ is parallel to dy. This
helps stabilize convergence in cases of marginal converg-
ence by not overshooting the local minimum. In cases of
divergence (as when the source is far outside of the receiv-
er array), limiting the gradient vector can allow one to get
an estimate of the direction of the source relative to the

35

initial starting point by tracking the path of all of the ap-
proximate solutions x,’ (away from the array). For stope
arrays with characteristic dimensions of about 100 m, a
good maximum to the magnitude of the gradient vector is
about 3 m. (Note that this is not the most efficient
method of reaching the minimum; see Press, (13), Sec-
tion 14.4 for improved strategies such as the Levenberg-
Marquardt method.)


has, however, recently been used to image kilometer-scale
velocity heterogeneities in Czechoslovakian coal fields (30).

In a hybrid approach, active and passive sources may
be combined to compensate for some of the disadvantages
of each method. Passive microseismicity sources can pro-
vide ray coverage near active working areas (.g., stope),
where it is extremely difficult to maintain electric cable-
transducer systems, and active sources can provide cover-
age wherever and whenever needed.

A three- or four-component velocity model (fig. 25) has
telative simplicity as a virtue; however, the effectiveness of

31

this approach must first be assessed by characterizing the
three-dimensional velocity structure in detail using some
other method, This approach may be attractive for routine
monitoring of microseismicity, particularly if it can be
demonstrated that there is little variation in velocity within
a given model element relative to the velocity differences
between elements,

An additional consideration in all of these approaches
is the need to continually update the velocity model as
mining progresses.

CONCLUSIONS

Various sources of error associated with the micro-
seismic event location process have been examined. Er-
rors in measurements of accelerometer position, travel-
time pick and average velocity, including the neglect of the
effect of mine openings, were estimated. These errors
were used with various numerical location techniques to
investigate the accuracy and precision of calculated event
locations. One numerical location technique has been
found that minimizes spatial location discrepancies for
synthetic events embedded in a spatially uniform velocity
medium with representative random and systematic errors.
This technique uses a direct solution basis function similar
to that of Godson (12). In general, iterative techniques
were found to minimize event mislocation when random
measurement errors are present in accelerometer-position
and arrival-time data. Event location errors produced by
systematic errors in the value of isotropic velocity were
found to be minimized when using direct solution meth-
ods, The best location solutions, in the presence of a
systematic velocity error, are not necessarily those pro-
ducing a minimum travel-time residual, as is most often
assumed,

Through the analysis of a series of calibration test
blasts, it has been shown that (1) the location precision
of blasts and microseismic events at the test site can be

+1 m and is currently limited by the precision of the first
arrival time picks; (2) the location accuracy, for both the
analog and digital systems, can be +10 m and is currently
limited by inadequate modeling of the velocity structure;
and (3) spatial variations in the apparent velocity exceed
50% due largely to velocity reductions associated with
fractured rock surrounding mine openings. To improve
location accuracy, the constant isotropic velocity model and
straight raypaths must be abandoned for a spatially
variable velocity structure with curved raypaths.

While the event location precision for the digital sys-
tem is approximately a factor of 5 better than for the
analog system, the difference in accuracy is much less
pronounced.

The observed accuracy is sufficient for routine
microseismic monitoring in the mine environment, but is
insufficient for detailed analyses of the deformation mech-
anics in these media.

With the increase in event location accuracy that is
possible with a more realistic (heterogeneous) velocity
model, significant improvements can be made in delin-
eating the position, orientation, and areal extent of dis-
continuities actively participating in mining-induced defor-
mation using microseismic event locations.

REFERENCES

1. Swanson, P. L., and F.M. Boler. Application of AE/MS Wave-
form Analysis to Hazard Detection, Evaluation, and Control in Under-
ground Mining. Paper in Progress in Acoustic Emission IV (Proc. 9th
Intl, Acoustic Emission Symp., Kobe, Japan, Nov, 14-17, 1988). The
Japanese Society for Non-Destructive Inspection, 1988, pp. 303-310.

2. Fairhurst, C. (ed.). Rock Bursts and Seismicity in Mines (Proc.
2nd Intl. Symp. Rockbursts and Seismicity in Mines, Minneapolis, MN,
June 8-10, 1988). Balkema, 1990, 439 pp.

3. Billington, S, F. M. Boles, P. L. Swanson, and L. H. Estey.
P-Wave Polarity Patterns from Mining-Induced Microseismicity in @
Hard-Rock Mine. Paper in Proceedings of the 31st U.S. Symposium on
Rock Mechanics, ed. by W. A. Hustrulid and G. A. Johnson (Golden,
CO, June 18-20, 1990). Balkema, 1990, pp. 931-938.

4. Goodman, R. E., and G, Shi, Block Theory and Its Application
to Rock Engineering. Prentice-Hall, 1985, 338 pp.

5. Coughlin, J. P., and C. D. Sines, Field Trials of a Portable Mi-
croseismic Processor Recorder, BuMines IC 9022, 1985, 8 pp.

6. Steblay, B. J., T. Swendseid, and B. Brady. Innovative Micro-
seismic Rock Burst Monitoring System. Paper in Rock Bursts and Seis-
micity in Mines (Proc. 2nd Intl. Symp. Rockbursts and Seismicity i
Mines, Minneapolis, MN, June 8-10, 1988). Balkema, 1990, pp. 259-262.

7, Steblay, B. J., B. T. Brady, and E. E. Hollop. A Networked
Minewide Microseismic Rock Burst Monitoring System. Paper in Min-
ing for the Puture-Trends and Expectations (Pres. at 14th World Min.
Congr., Beijing, China, May 14-18, 1990). Pergamon, 1990, pp. 861-867.



The location based on the spatial-gradient basis
function (appendix A) is now computed. By way of
illustration, the starting point for the iterative process is
selected to be the revised SW-GMB solution from above
(ta = 0.039074738 s, x, = 3412,905762 m, y, = 2798,638184
m, and 2, = -362,666046 m), though almost any starting
point in or around the receiver array could have been
selected as a starting point. For the first iteration, the
linear system Ax = b is set up and solved, by the QR
algorithm:

0.189479 0.060893 0.025504 0.778861
-.063410 120467 =.111098 849044
-190580 ~.049664 030849 105471
043424 .213069 027282 501761
A= -101470 ~100726 128624 b = -.812930
176394 081186 014134 485710
195818 +.036953 015256 +.150945
-038918 +173794 012140 453275
001221 119656 +.160850 968294
127313 169171 018160 059605
with

+1.857610

Xx = + 756399

+ 598912

The solution vector x indicates the change in the estimate
of spatial coordinates of the source. The first correction
to the origin-time estimate is obtained independently via
equation (A-9) and is found to be -0.000032 s. After mak-
ing the corrections to the source origin time and coordi-
nates, 10 new equations are formed using the spatial-
gradient basis function and solved iteratively until the
corrections to the source coordinates are sufficiently small
(approaching zero). The final iteration of the spatial-
gradient basis function for this example gave a solution of

* US.QPO: e-611oreenees

t, = 0.039026 s, x) = 3410.91 m, y, = 2797.77 m, and z =
-363.41 m. This solution yields the following travel-time
residuals:

station P(calo)—-P(o-c)
m 0.044603 0.000477
13 .044520.000540
14.1 041377000243
6 043074 — .000606
7 045454 = -.000194
ic} 050914 000226
m4 042932 = .000228
fal) 059755 -.000975
m2 047473 -,000153
5.058209 -,000999

Thus for the spatial-gradient basis function, the effective
mms of the travel-time residuals is reduced to 0.000553 s
normalized by 1/m, or 0.000714 s normalized by 1/(m-4).

The final solution for the full-gradient basis func-
tion (appendix A) yields a solution of t, = 0.038555 s,
% = 3414.46 m, y, = 2800.86 m, and z) = -359.70 m, giv-

ing the following travel-time residuals:
station P(calo) Plo-c)
2 0.043688 0.001412
3 042910 .002150
14.1 041500 .000120
6 043108 .000571
7 045322 -,000062
8 050976 .000164
19.1 041750 .001410
0 058944 -.000164
n2 047039 — -.000281
5 057521 -.000221

Thus for the full-gradient basis function, the effective rms
of the travel-time residuals is 0.000955 s normalized by
1/m, or 0.001233 s normalized by 1/(m-4).

INT.BU.OF MINES,PGH.,PA 29479


In the following sections, the magnitude of each uncer-
tainty term on the right-hand-side of equation 6 is esti-
mated from field measurements for both the analog and
digital systems. To examine quantitatively the influence of
these measurement errors on three-dimensional locations,
several common source location techniques are first
examined, and then one method is selected which mini-
mizes bias introduced by numerical methods.

NUMERICAL SOURCE LOCATION TECHNIQUES
Basis Functions

The spatial location (X, Yo, %) and origin time (t,) of
microseismic events is desired, For sources occurring in
a constant velocity medium, equation 1 relates the travel
time to the distance between the source and the ith re-
ceiver. For m arrivals, there are m nonlinear equations
that must be solved simultaneously (m linear equations for
the one-dimensional case). Since this is difficult to do in
a general way, the usual approach is to find some system
of linear approximations which can be solved. The term
“basis function" or basis will be used to describe the equa-
tion actually used in the system of equations to be solved.

For applications typical of the mining environment in
which a spatially constant velocity field is assumed, at least
two methods have been described (9, 11-12) where the
source location can be found directly by solving a system
of linear equations, which is obtained after some algebraic
manipulation of equation 1. The direct equation develo-
ped by Salamon and Wiebols (11; which was originally in-
troduced to provide a starting solution for an iterative
technique) and Godson (12) will be referred to as the
SW-GBM basis; that developed by Blake (9) will be re-
ferred to as the BLD basis. Generalized forms of the ba-
sis functions for these two direct methods are given in the
appendixes and will be discussed in more detail in a sepa-
rate publication, Use of these direct solution methods
requires special care.

The generalized form of the SW-GBM basis yields
m(m-1)/2 different linear equations (equation A-1), since
each depends on information from two different receivers
(arrival times of the signal and receiver coordinates).
From these m(m-1)/2 equations (A-1), many different
subsets of m-1 to as many as m(m-1)/2 linearly independ-
ent equations can be selected depending upon whether or
not weighting is used. The generalized form of the BLD
basis yields m(m-1)(m-2)/6 different linear equations
(equation A-6) since each depends on information from
three different receivers and three different equations can
be obtained from each set of three receivers; of these
m(m-1)(m-2)/6 equations (A-6), many different subsets
of (m-1)(m-2)/2 to as many as m(m-1)(m-2)/6 linearly

independent equations can be selected, depending upon
whether or not weighting is used. Use of the SW-GBM
and BLD bases requires at least five arrivals at different
receivers to obtain a meaningful three-dimensional loca-
tion, At least six arrivals are needed to determine some
measure of confidence in the location accuracy. For even
modest numbers of receivers, these direct solution meth-
ods can result in (1) many linear equations to be solved
simultaneously, or (2) some selection of a particular subset
of equations out of many different possibilities,

Solving large systems of equations (number of equations
>> m) is obviously undesirable. However, if a particular
subset of the equations of one of the direct basis functions
is used, there is a risk of introducing a bias into the
solution, For example, Blake (9) considered a particular
subset of m-2 equations with the BLD basis, where all m-2
equations of their proposed subset depend on the arrivals
at the first and second hit receivers (as well as a third
arrival at the other receivers). Errors in the information
from the first two stations (i.¢., errors in the arrival times
and/or receiver coordinates) overinfluence the source solu-
tion resulting in unnecessary bias. In a similar fashion,
Salamon and Wiebols (17) considered a particular subset
of m-1 equations with the SW-GBM basis, where all m-1
equations of their proposed subset depend on arrivals from
the first hit receiver; this, like Blake’s use of the BLD
basis, results in a solution overinfluenced by the first hit
station. Godson’s (12) implementation of the SW-GBM
basis is considerably better, since all arrivals but two are
represented twice, but there is still a slight underrepresen-
tation of information from these two stations. In appen-
dix B, some of the different subsets of the direct basis
functions are discussed further.

In general seismological applications, where the velocity
field is often known to be heterogeneous (spatially vari-
able) and in some cases anisotropic (directionally variable)
and there is considerable incentive to accurately model this
velocity field, iterative solution techniques are widely used.
The most common approach is to apply some form of the
gradient technique (13). A Taylor series expansion can be
used to find a linear approximation which will be called
the gradient basis. The gradient basis is used to form a
system of m simultaneous linear equations for m different
arrivals at different receivers, where the unknowns are now
corrections to some approximate solution which will min-
imize some function, e.g., equation 1. These corrections,
icc. the gradient of the function being minimized, are then
added to the approximate solution to find a new approxi-
mate solution, The process is repeated until the correction
vector (the gradient) is essentially zero or some other cri-
terion is met. A gradient basis of this type is given in the
appendix equation A-8 and determines the source location
and origin time; this will be referred to as the full-gradient


channel for each event. Hardware-determined arrival
times are thus suited to near-real-time calculation of event
location, Another advantage is that computer mass stor-
age can be conserved. A variety of electronic circuits can
be used in hardware timing. The errors in picks can be
attributed to both the inherent errors in the electronics
and the inability of the hardware to readily handle wide
variations in signal shapes and signal-to-noise ratios. A
simple voltage threshold exceedance is used for the analog
system. For a signal whose onset approximates a step
function, this method will have an inherent error equal to
the time resolution, including channel synchronization,
designed into the electronics. This resolution is 100 us for
the analog system (23). For signals which are emergent,
i.e. which have a relatively long rise time, the time to the
threshold exceedance is correspondingly longer, and the
pick will be in error by an amount proportional to the rise
time. As high frequencies are selectively attenuated with
propagation distance (¢.g., fig. 5), rise times increase; thus
picking errors systematically increase with wave propaga-
tion distance. Systematic timing errors also occur at in-
dividual receiver stations when they experience low-pass
filtering due to poor coupling into the host rock and thus
produce emergent signals,

Another common source of error results when the volt-
age threshold is not exceeded upon receipt of a low-
amplitude compressional wave arrival, yet is subsequently
exceeded by the shear wave or other later-arriving phase;
this is potentially a systematic source of arrival-time error.
The magnitude of this error is proportional to the distance
to the event and the difference in velocities of the two
phases.

Manual Timing

If a permanent record of the signals exists, either on
chart paper, in a digitizing oscilloscope memory, or in
computer mass storage, manual arrival-time picking from
visual examination of the record is possible. (Both manual
and software timing offer the advantage of being able to
reprocess the data if permanent records are kept.) The
errors inherent in this method will correspond to how the
signal is stored, the skill and consistency of the person
doing the picking, and the quality and shape of the signal.
Further discussion of manual picking will be restricted to
the case of computer mass storage of digitized signals,
which can then be displayed at any desired magnification.
Assuming no other errors in the data acquisition and
display electronics, the digitization interval will be the
smallest inherent error in the pick. For the digital system,
this interval was 20 ys for the blast tests and 10 ys for the
microseismic monitoring. The person reading the records

19

will generally be able to judge from relative arrival times
which phases are reasonable to pick, and is unlikely to
inadvertently mix picks of different phases. There will be
some variability in the absolute picks from different per-
sons (pickers), but the internal consistency can be very
high. Similar statements apply to the assignment of quality
descriptors to a given pick. The picker can readily say
whether a low-quality pick is due to a low signal-to-noise
ratio, or to an emergent signal, or to some other cause.

Software Timing

Software picks for arrival-time determination are similar
to hardware picks in that there are many picking algo-
rithms (24). A distinct advantage is the ability to cus-
tomize, at any desired time, a particular picking algorithm
to suit a wide variety of (possibly time-variable) signal con-
ditions, Assignment of quality is also possible, and can
overcome some of the difficulty in picking first arrivals
using an algorithm that may not be optimized for certain
signal types.

Example Comparison

Two methods of obtaining arrival-time picks were com-
pared for ~40 microseismic events which had acceptable
location quality. The analog system used a simple voltage
threshold exceedance hardware pick which was subse-
quently screened in software as described earlier. The
digital system used interactive manual picks from visual
examination of digitized waveforms for the same set of
signals used for the hardware picks.

The analog and digital systems were not tied to the
same clock, It was necessary, therefore, to arbitrarily
adjust the timebase from one set of arrival-time picks to
the other. The hardware pick method sets the first hit sta-
tion (accelerometer with the first exceedance of threshold
voltage) to time zero and determines all the other arrival
times relative to the first hit accelerometer. The digital
waveform system assigns a time to each event from the
computer clock.

To obtain a comparison between these two methods, it
was assumed that there was no difference between the two
pick methods at least one accelerometer (not necessarily
the first one hit). The accelerometer with no discrepancy
was chosen such that (1) the arrival was impulsive to avoid
an obvious late pick associated with an emergent signal,
(2) tying to that station did not result in a hardware pick
which arrived ahead of the manual pick by more than
10 ps at any other station, and (3) the differences between
the two sets of picks were otherwise minimized.


36

APPENDIX B.—SUBSETS OF DIRECT BASIS FUNCTIONS

Each of the above direct basis functions can be used to
form a’system with more than m equations. Whether or
not the full system or a subset of the system of equations
can be used depends on the weighting of each basis
function.

SW-GBM BASIS SUBSETS

Each nonzero SW-GBM basis function involves infor-
mation from two different receivers, j and k. (If the
indices j and k are interchanged in equation A-1, the same
equation is formed.) Therefore, the set of all possible
SW-GBM equations can be referenced by specifying all
possible pairs of j and k. As an example, let m = 5:

2345

1

12
1 3
1

1

4
5
34

YEE

5
465

It is easy to see that there are m(m-1)/2 different
pairings since an even permutation of the indices Kl results
in the same equation (A-1) and an odd permutation re-
sults in only a sign change of the equation A-1. The
m(m-1)/2 equations, however, are not necessarily inde-
pendent. It is also fairly easy to see that a subset of as few
as m-1 independent equations can be selected as long as
all of the information is weighted equally. For example,

12 3 4 5 1 2
12 12
1 3 2
1
1

4
5 4
represent two different subsets of m-1 independent equa-
tions which can be used in linear combination to form all
the other (m-1)(m/2 - 1) equations (assuming equal
weighting). Suppose that the ordering of the receiver
indices in the above example indicates the order in which
signals arrived at the receivers (positive arrival-time
ordering). Then the first m-1 subset above represents that
suggested by Salamon and Wiebols (11), and the second
m-1 subset above represents that suggested by Godson
(12). Notice that the first subset drastically overem-
phasizes the information from receiver 1; imagine the
results if the receiver coordinates and/or the arrival-time

pick of receiver 1 are significantly in error! Note that
these two cases are not the only choices of m-1 independ-
ent equations; for example,

123 4 5 #1
2

1
1

wun]
@

3
3.4
4 6 45
represent only two more of the many different possible
subsets of m-1 independent equations that can be selected.

A subset like that suggested by Godson results in a set
of equations that are information poor for two receivers.
This can be altered by always including an equation for
that receiver pair, making m equations:

1
1

3445

vv
oo

4
45
1 5

For such subsets, one equation can always be formed by
a linear combination of the other m-1 equations, but
the representation of information from each receiver is
balanced.

So far, this discussion has assumed that all the in-
formation for the m receivers has equal weight, ie., all
arrival times and all station coordinates have the same
relative uncertainties—or at least there are no criteria by
which to assign different weights. If, however, different
weights, Wy are given to each basis function A-1 (i.e, the
weighted basis function is wy (A-1)), then, in general,
there could be as many as m(m-1)/2 independent SW-
GBM basis functions. Here w, represents some weighting -
function for the information for receiver j combined with
the information for receiver k.

By using different subsets of the SW-GBM basis func-
tion, it is easy to see that there is a staggering number of
different possible source location solutions that could be
obtained. If all the basis functions are weighted inde-
pendently, resulting in m(m-1)/2 independent basis equa-
tions, then there are

m(m=1)/2
>>

[mo], @.)

re4

n}]l_oat  , n(a-t) (n-2)...(a-r+1)
where ("] FC; |

possible different solutions.


. Influence of mine openings on seismic wave travel times . .

). Influence of 4% systematic velocity error on calculated location of 1,000 defin

.. Vertical cross section through four-component velocity mo

. Performance comparison between numerical source location methods, hypothetical idealized array .

. Difference between hardware and manual picks of first arrival times at nine stations .
.. Average calculated location-discrepancy vector magnitudes

ILLUSTRATIONS—Continued

Histogram showing difference between hardware and manual travel-time picks .

Comparison between 40 event locations calculated using hardware and manual picks of ar
with all other parameters identical 60.0.0... 0... cccsssseeee eee cceeeeeeneeeeeeeeeneues

Comparison between 1,000 defined event locations and locations calculated with hardware-based first-
arrival time pick errors added to exact travel times ... an

Influence of mine-opening-related travel-time delays on calculated event locations .
Apparent seismic velocities calculated for straight raypaths in calibration blast tests . .

). Measured travel times and straight-line distances between all source and receiver sites for calibration

blast tests

VOLUME oe eee e eee e nent eee e eee et EE eee ee EEE EEE eee EEE EE EEE EE

. Influence of 13% systematic velocity error on calculated location of 1,000 defined events in control

VOLUME cece eee eect ete e eee e eee EEE ENE E EEE EEO O Oe EEA E ENED EEE E EGE e EEE
Composite plot for analog system showing location discrepancies calculated using measured errors in
accelerometer position, arrival-time picks, and assumed isotropic velocity ..........sseeeeeeeeees

. Composite plot for digital system showing location discrepancies calculated using measured errors in

accelerometer position, arrival-time picks, and assumed isotropic velocity...........0..cseeeeues

. Calibration blast locations from digital system illustrating magnitude, shape, and orientation of 95%

confidence SVD error ellipsoids .

TABLES

Performance comparison between numerical source location methods, 4300-115 array

. Difference between surveyed and tape-and-map-measured accelerometer coordinates before and after

transformation

Page

8 RF RKPS BB

8

BB

BBS

UNIT OF MEASURE ABBREVIATIONS USED IN THIS REPORT

cm centimeter m meter

dB decibel m/s meter per second

ft foot us microsecond

ft/s foot per second ms millisecond

kHz kilohertz mV/g__ millivolt per acceleration
km Kilometer 8 second

km/s kilometer per second



If the significant sources of uncertainty in microseismic
event locations have been correctly identified and char-
acterized, then the discrepancy vectors in figure 23 should
be consistent with the observed discrepancies between
surveyed and calculated calibration blast locations. Fig-
ure 24 shows the calibration blast locations for the digital
system and axes of the SVD spatial uncertainty ellipsoids
of 95% confidence. A comparison between figures 23 and
24 reveals that the calculated and observed discrepancy
vectors are not equal. Near the calibration blasts, the
composite discrepancy vectors estimated for the digital
system (<1 m; fig, 23) are much smaller than the observed
discrepancies (~10 m; fig. 24). This discrepancy is pre-
sumably attributable to the fact that the heterogencous
seismic velocity structure at the experiment site does not
adequately conform to the assumed constant velocity
model.

Other observations from figure 24 are noted as follows:
(1) The computed 95% confidence error ellipsoids do not
encompass the surveyed position of the blast site. As
previously noted, the ~10 m offset of the surveyed position
from the calculated locations indicates a location accuracy

A va

\

a
0 10
a
Scale, m
B
A eT

Figure 24.—Calibration blast locations (skewed triad symbols)
from digital system illustrating magnitude, shape, and orientation
‘of 95% confidence SVD error ellipsoids. Surveyed blast position
Is represented by filled triangle. A, Plan view; B, vertical section
looking northwest.

29

of approximately +10 m. (2) The precision of the loca-
tions (+1 m) is limited by the precision of the travel-time
picks and is illustrated by the +1 m clustering of four of
the calculated blast positions. These four locations were
determined using arrivals from the same set of 17 re-
ceivers with only slight variations in travel-time picks.
(3) A smaller subset of 10 arrivals was available for the
location of the fifth blast, which is offset from the other
four. As the relative arrival times were the same as for
the other four locations, this offset indicates the significant
influence of receiver array geometry. (4) The error ellip-
soids of all of the blast locations have the same approxi-
mate magnitude, aspect ratio, and spatial orientation.

The fact that none of the spatial uncertainty error
ellipsoids encompass the surveyed blast position is consist-
ent with the contention that the isotropic velocity model
does not adequately conform to the actual velocity struc-
ture, Recall that the SVD error ellipsoid provides confi-
dence limits on the location uncertainty if the errors are
normally distributed. Since the deviation of a hetero-
geneous velocity structure from a uniform velocity model
does not in general result in normally distributed velocity
errors (or resultant measured travel times), the calculated
SVD error ellipsoids do not, as observed in figure 24,
adequately describe the constraint on the source locations
in this situation. Use of a more realistic (heterogeneous)
velocity structure is required so that the error ellipsoids
more closely reflect the constraint on the solution due to
the actual measurement errors.

A location discrepancy of 10 m is not expected for all
controlled-source tests. The various factors providing con-
straint on the source locations (¢.g,, array geometry, num-
ber of reporting receivers, degree of local deviation of ve-
locity structure from assumed model, etc.) vary over all
combinations so that a single number cannot effectively
characterize source location accuracy for all sources, Con-
straint on the source location solution will be better than
10 m in some situations and much worse in others.

In summary, the failure of the isotropic velocity model
to adequately describe wave propagation at this site ap-
pears to represent the largest systematic source of uncer-
tainty in event locations for both the analog and digital
systems. To significantly increase the accuracy of event
locations in any part of the mine covered by the array, a
more accurate representation of the velocity structure is
required.

While considerable effort toward an improved velocity
model for research purposes is warranted, continued use
of constant isotropic velocity models in routine microseis-
mic monitoring is appropriate. The value of this health,
safety, and productivity tool is negligibly affected by the
observed location uncertainty, at least when an appropriate
average velocity is used. An event location accuracy of
+10 m in a mine where working areas are separated by
hundreds of meters and are individually monitored, readily


16

Subsets of these receiver stations also tend to take on lin-
ear distributions within this plane (e.g,, fig. 4) due to re-
stricted access (i.e., drifts, raises, crosscuts, etc.), As the
number of receivers reporting usable first arrivals is re-
duced, due (¢.g., to low amplitudes of small magnitude
events), the probability of forming a predominantly linear
or planar array increases.*

Figure 9 illustrates how the constraint on the solution
can be influenced by the geometrical distribution of a sub-
set of reporting stations. The overall receiver array is the
same as shown in figure 4, Figures 9A and 9B show two
different sets of five receivers which are used in graphical
locations of the same event. Similar measurement errors
were added to the exact data from each receiver. From
the angles of intersection of isochrons near the source
location, it is immediately recognizable that there is much
more uniform constraint for the receiver distribution in
figure 9A than in 9B.

A comparison of the calculated event locations and
their associated SVD error ellipsoids is shown in the figure
inserts (lower left-hand corners), As expected from the
isochron intersections, the SVD error ellipsoid in figure 94
is quite equidimensional compared to figure 9B. In both
cases, the calculated event location is displaced away from
the exact location along the direction of the long ellipsoid
axis, The magnitude of the mislocation observed with the
linear array subset of figure 9B (1.8 m) is approximately
2.5 times that observed with the array subset that sur-
rounds the source in figure 9A (0.7 m). This illustrates the
strong influence of source and receiver-array geometry on
the sensitivity of a location position to measurement
errors.

While the array subset in figure 94 appears to yield
adequate constraint in the horizontal plane, there is no
constraint for either array subset in the vertical direction.
Thus the mislocations represent an absolute lower bound,
The reporting receivers need only weakly approximate a
linear or planar distribution to form poorly constrained
solutions.

"This aspect of the effect of array geometry is partly simulated in the
error analyses in previous and subsequent sections by randomly ex-
cluding 15% of the stations in each simulated-source location calculation
for both the analog and digital systems. On average, 7.5 stations were
used with the digital system in locating the microseismie events shown
in figure 7, Occasionally, only five stations were used (less than 5% of
this particular data set). Since first-arrival picking with hardware results
in significantly fewer usable signals than with manual waveform picking,
the number of reporting stations for the analog system has been overes-
timated. As a result, the location discrepancies reported for the analog
system have been underestimated.

0

W——_ Scale, m
Scale, m

B
Oo 20

0 5 uN
—— Scale, m
Scale, m

Figure 9.—Influence of receiver array geometry on event
location errors. A synthetic source location (open circle In
insets) Is calculated (center of SVD ellipsoid axes) using small
amounts of error added to exact travel-time data for two different
arrangements of five (circled) receivers. Relative-arrival-time
Isochrons are shown for comparison. A, Uniform receiver cover-
age yields uniform Isochron intersections and equidimensional
SVD error ellipsoid; B, quasi-linear array ylelds large elongate
error ellipsoid and larger misiocation for equivalent measurement
errors.



Table 4.—Difference between hardware and manual picks of first arrival times at nine stations, milliseconds

Event 1 3 4 5 6 7 8 9 12
0.010 - T 0.030 0.030 0.130 ad -
+040 T 0.010 +170 ~ -070 = 0.020 -
+270 - T +090 - +250 +200 ~

.000 0.010 +010 =110 - ~020 .000 -

+120 010 T +120 - +110 A -

- - T +130 +150 - -

~.090 - T +060 +500 010 -

+140 - +030 +240 - ~050 - -

+370 - 000 +290 - +050 +020 -

+110 020 T +290 ~.060 +210 000 -

+320 +010 T +100 +420 +060 +230 -

- +030 010 T =.100 040 ~.060 -

T - +030 ~.060 +240 +010 +060 -

T +74 +030 ~.030 +140 = 110 490 -

T - 000 010 +290 070 +090 -

T +350 +050 050 +270 000 =130 -

T 010 +040 000 +160 +140 +070 -

T - - 020 +130 020 ~040 -

T - 000 ~.030 +190 000 ~840 -

- - - +010 +210 T +030 -

- - - 010 +170 =020 +050 -

~090 - - T +320 ~050 +070 ”

+170 - +060 +100 ~580 +410 180 -

010 - - T +060 -.030 =010 -

.000 - T 010 +550 +050 ~.060 -

+010 - T 010 460 +030 060 -

+130 - T 000 +130 050 +150 -

- - -.030 T +210 ~.020 ~.280 -

+010 - T ~.020 +310 100 +060 -

010 - T “170 +250 000 +260 -

T - 040 +040 050 +020 +010 ~.260 -

060 - T +040 +150 +030 +040 040 -

.000 - 1.970 T 290 +130 +080 010 -

~.060 - T 000 - +030 100 +040 -

020 - - T +190, +500 +540 -.020 -
~.080 - T 000 +030 - +050 +020 0.290
- - - 040 - - 1,070 +050 T
020 - - +010 310 - =060 T +200
.000 - 010 T 130 - +570 - +240

010 - -.040 T 230 = +050 ~.360 -
Mean . 050 164 097 031 213 075 183 075 243
Std dev .. 077 208 287 .039 132 116 233 +103 045

- Arrival time not available for one or both systems,
T Time bases tied at this station.

An example of the effect of the difference between
the hardware-determined arrival-time picks of the ana-
Jog system and the arrival-time picks determined by vis-
ual inspection of waveforms is seen in figure 14. Here, a
subset of the 52 events, which were well-located by both
the analog and digital systems, were relocated using the
same values of velocity and accelerometer coordinates and
using identical numerical location techniques. Arrival-time

measurements from both the prescreened analog picks and
the manual picks were culled to include only data from
accelerometers in common, The differences between the
locations shown in pairs in figure 14 are solely due to
the difference between hardware arrival-time picks and
unweighted manual picks of arrival time. The magnitude

of the average three-dimensional discrepancy vector is
3.7m.


The differences in simulated microseismic event loca-
tions, caused by position errors of the same magnitude as
the differences in the two sets of accelerometer position
measurements, are illustrated in figure 10. This figure
shows 1,000 defined event locations selected at random in
a volume that is typically monitored by an array such as
found in figure 4. Exact travel-time data were calculated
from the distance between the defined event locations and
the surveyed station coordinates using a constant velocity,
In figure 10, event locations calculated using these travel-
time data and the tape-and-map-measured accelerometer
positions are connected with lines to the exact locations:

The interconnecting lines graphically illustrate the
spatial variation in location discrepancy; the flow-line trend
is generally consistent with the trend of major axes of the
calculated SVD error ellipsoids. Of course, the discrep-
ancies become greatest when the event is located outside
of the array volume. In most routine monitoring situa-
tions, due to the restrictions on where receivers can
conveniently be located, many events occur below most, or
all, of the receivers in a typical array (ie., in the active
workings of the stope) where constraint in the vertical di-
rection is typically at a minimum? The average magnitude
of the three-dimensional discrepancy vector connecting the
two locations in the particular volume of figure 10 is 2.3 m.

This number is very dependent upon the particular
accelerometer array geometry and the particular control
volume selected and is thus identified only to assist in
making relative comparisons of discrepancies calculated
using the same array and control volume. The same sur-
veyed accelerometer array and random event volume is
used in the next several sections to allow comparison of
the three-dimensional location errors resulting from errors
in travel time and velocity. For these comparisons, the
surveyed station coordinates are considered to be exact.

ARRIVAL TIME DETERMINATION

Arrival time determination methods can be divided into
hardware timing (picking), manual timing, and software
timing. Each method has advantages and disadvantages.
A related topic is the assignment of some sort of quality
descriptor to the pick which is used to weight the arrivals
in the numerical solution for the location. Weighting is
possible, in theory, for all of the above methods.

The several very long vertically oriented discrepancy vectors in the
vertical section may indicate a mirroring of the solution, or convergence
to one of several possible solution sites. In general, such behavior was
observed when either increasing the measurement errors or decreasing
the number of stations used in the location,

4 ns
ht ifn ae

Figure 10.—Defined event locations selected at random in
control volume centered on receiver array (from fig. 4) connected
with straight lines (flow lines) to event locations calculated with
accelerometer position measurement errors. Flow lines indicate
directions of least constraint on event location solutions. A, Plan
view; B, vertical section looking northwest.

Hardware Timing

Using hardware electronics to examine transient voltage
signals at each accelerometer offers the ability to deter-
mine arrival times without tedious examination of every


30

allows the identification of the specific working area that
is associated with microseismic activity (and those that are
not). In individual stopes of dimensions 60 by 100 m, it
can still be determined approximately where in the stope
this activity is taking place. Furthermore, the deter-
mination of rates of microseismic activity for daily mine
stability assessments, are not at all affected.

The improvements in event location accuracy which are
attainable by reducing random measurement errors, €.g.,
surveying the accelerometer coordinates and manually
picking first-arrival times, are small in comparison to the
improvements which are possible if the spatial variation in
seismic velocity is characterized and used in the event
location process. The increased effort required by the
precision measurement procedures is not warranted for
routine mine monitoring purposes until more accurate
seismic velocity models are developed. For the develop-
ment of rock burst prediction and control strategies based
on an understanding of the mechanics of rock burst proc-
esses, increases in event location accuracy are essential.
For example, since large stress gradients are sustained
over distances of less than 20 m, especially near mine
openings, discontinuity structures cannot be located con-
fidently in even the appropriate stress regime with a
location accuracy of +10 m. Thus, mine stability analyses
based on the proposed uses of microseismicity (e.g.,
discontinuity delineation, determination of slip directions,
etc.) require an increase in event location accuracy.

There are several possible approaches for reducing
event location uncertainty via characterization of the in-
homogeneous velocity structure, Seismic wave travel times
can be measured in individual monitoring areas using ac-
tive seismic sources (27-29), and/or relative arrival times
can be measured from passive seismic sources (30) and
used in curved raypath tomographic reconstructions of the
three-dimensional velocity structure. Alternatively, a sim-
plified model structure using constant values of velocity
(fig. 25) representing intact rock, air in mine openings,
sand in backfilled areas and a fractured halo region sur-
rounding mine openings, may provide acceptable improve-
ments in event location accuracy.

There are advantages and disadvantages to each ap-
proach, Active seismic tomography allows the controlled
placement of seismic sources at known positions and
known times. For three-dimensional velocity structure,
however, source-receiver placement must also be three-
dimensional, which is difficult and expensive since long
holes through hard rock generally must be drilled. High
resolution in areas of strong velocity gradient is also
difficult to achieve unless many measurements are ob-
tained at closely spaced source and receiver positions.

The use of passive sources for tomographic velocity
reconstruction, such as microseismicity sources, has the
advantage that these seismic sources are generally plentiful

vy

Figure 25.—Vertical cross section through four-component
velocity model: V, Intact rock, V, alr, V; sand, V, fractured rock.

_ Geometry Is typical of mine openings near stoping areas.

and usually accompany the mining process. Also, mines
which are presently monitoring microseismicity do not
have to acquire additional hardware for data acquisition.
A distinct disadvantage in this particular application is that
many events are required to be located in the specific
volume of rock to be characterized. Thus the large ve-
locity gradients near sensors mounted to ribs in previously
excavated mine openings will not be imaged effectively us-
ing these sources since most high-angle raypath intersec-
tions, which largely determine resolution, are roughly
limited to the seismic source region (i.e., stope). Further-
more, one does not have (ready) control over when these
sources are activated, Simultaneous determination of
event hypocenter and three-dimensional velocity structure
is also a much more difficult computing task. This method


Table 2.-Performance comparison between numerical source location methods, 4300-115 array

Random errors’ ‘Systematic errors”
Method (0,000020-8 seed 0.10-m seed 1.00-m seed 5.2-km/s velocity

,__
,m_<Jdt]>,8_<Jdx[>,m__</dt[>,8<|dx]>,m__<jdt[>,8_
,m
o 0.000011 0.146400 0.000114 1.500710 0.000414 1
807230 000010 142696 = 000111 1.462420 ©» 000420 1.843113
681949 000013 000131 1.600344. .000401 1.001503
(687508 000012 000129 1.582105 .000409 913786
Al cccceeeeees 000074 x 900011146412 .000114 1.500718 = 000414 1.436337
Fis. 000071 807285 = 000010, 142737 000111, 1.462448 000420, 1.849138
(m-2)BLD 000075 627874 022649 7413.42168 3005678 729.49322 31018368 °95,147606
[000013 164976] [.000131 1.674416) [ 1.837818)
m, stag’ , 900111 994049 1000022234148 000220 2.359569 .000410 1.987207
(m-2)stag*. . 000133 1.222574 000020220254 000200 2.271821 = .000415 850042
Spatial gradient 000052 “377889 000007096421 000076 969821 000480 2.166870
Full gradient .. . 000060 467821 000011 130575 000116 ___1.319788_.000495 3.201569
I m/s vel
0.10-m seed, 1,00-m seed
<idt|>,8__
,m_<Jdt[>,s_<|dx]>,m
++ 0.000479 1.866042 0.000512 2.112272
+ 000482 2.238079 000518 2.391626
000459 (1.477663 0004911
000463 «1.423559 000498
. 000479 1.866079 000512
Fis. 000482 2.298074 000518
(-2)BLp 000489 264748 000528
m, stag’ 000400 2.599844 000559
(m-2)stag’ ... 000517 1.818175 000706
Spatial gradient 000534 2.675963 000539
Full gradient . + __.000558_3.579921 000559

 Mean of absolute differences between synthetic and calculated source origin times.

 Mean of absolute differences between synthetic and calculated source locations.
BLD Blake-Leighton-Duval generalized basis function (see text),
SW-GBM —_Salamon-Welbols-Godson-Bridges-McKavanagh basis function (see text).
FiS Full independent subset (one possiblity, see text).

1,900020-s seedargest Possible error for stations approximately within 20 m is 0.020 ms, within 20 to 40 m Is 0.040 ms, within 40 to

100 m is 0.080 ms, and with

In 100 to 200 m Is 0,160 ms. 0,000080-s seed—4 times the time uncertainty of 0.000020-s

seed. 0.10-m seed

largest spatial error of any Cartesian receiver coordinate Is 0.10 m. 1.00-m seed—argest spatial error of any Cartesian receiver co-

ordinate Is 1.00 m.

24.8. and 5.2-km/s velocity—used for average P velocity rather than 5.0 km/s. 5.64-km/s velocity—6.6388 km/s (18,500 ft/s) used for

average P velocity rather than 5.0 km/s.

31 location was mislocated by a very large amount, Bracksted data are averages based on 999 locations (large mislocation removed),
‘Staggered ordering; that is, triplets 123, 234, 345, .. (m-2)(m-1)m (see appendix B).

Summary and Recommendations

When the errors are randomly distributed (such as is
the case for most receiver coordinates and for arrival
times), a good approach is -

1, use one of the direct bases (either SW-GBM or
BLD) to form a subset of linear equations (¢.g., using one
of the subsets suggested in appendix B), solve using a good
QR or SVD code, and obtain a first approximate solution
for the location;

2. use a gradient basis with this approximate solu-
tion to form new system of m linear equations, solve for

gradient using QR or SVD code, form new approximate
solution, and iterate until gradient is essentially zero or
other criterion is met (note: if a minimization of travel-
time residuals is desired, use the spatial-gradient basis);

3. if error estimates are desired (assuming that all
errors are normally distributed), perform one (or two)
more iterations with SVD code and compute confidence
regions of the solution.

When the errors are mainly systematic, as in the case
of modest deviations of the trial velocity from synthetic
known velocity, the best approach may be to use the SW-
GBM basis with the staggered subset suggested originally


Figure 4.—Plan view of uniaxial (filled squares) and triaxial (open squares) receiver array. Accelerometers 3, 6, 10, and 15 are,
respectively, +19, -9, +93, and -53 m above and below the 4300 level.

uniform rock structure, or absence of fracturing and other
sources of heterogeneity, abrupt first arrivals, etc.) may not
apply at all to microseismic events occurring in a different
part of the stope covered by the same array, particularly
towards the edge of the array. For this reason, controlled
source tests were conducted in the stope toward the bot-
tom edge of the array, approximately 5 m below the main
level.

Calibration blasts were made by detonating one-quarter
to one-half stick of dynamite placed untamped in a

2-m-long drill hole. With the analog system, the blast
waves were detected by the 13 rib-mounted uniaxial
accelerometers (fig. 4). Waveforms from both the uniaxial
and triaxial accelerometers were recorded with the digital
system. An electrical contactor placed in the dynamite
with the blasting cap was used as a trigger signal for the
digital system. The time history of the trigger signal was
recorded and used to determine blast origin time. An ex-
ample of the blast waveforms and the trigger signal are
shown in figure 5.


Figure 20.—Influence of 4% systematic velocity error (5,200 vs.
5,000 m/s) on calculated location of 1,000 defined events in con-
trol volume. A, Plan view; B, vertical section looking northwest.

limit of location uncertainty associated with a systematic
velocity error in a truly constant velocity medium.

A velocity of 5,640 m/s was previously determined for
use with the analog system in similar calibration blast sur-
veys in a different part of the Galena Mine. Observed ap-
parent velocities ranged from 4,200 to 6,700 m/s. The me-
dian value of 5,640 m/s (dashed line in figure 19), which
is a representative value for high-velocity undisturbed

27

Figure 21.—Influence of 13% systematic velocity error (5,640
vs. 5,000 m/s) on calculated location of 1,000 defined events
In control volume. A, Plan view; B, vertical section looking
northwest.

paths, has subsequently been used as the velocity for the
mine-wide-network monitoring experiment. Figure 21 il-
lustrates the location discrepancies which result when
5,640 m/s is used in an area where 5,000 m/s may be a
better average velocity (i.e., a 13% systematic discrepancy
in the velocity), The average magnitude of the location
discrepancy vectors is 7.0 m. This value is reduced
to 2.1 m when a 4% error is assumed for an average ve-
locity of 5,000 m/s and events are located using the BLD
method in the analog system.


Since the first signal detected at the receiver travels
along a minimum-time travel path (Fermat’s principle), the
fastest path, in the short wavelength limit, is through the
solid rock medium around the mine opening. This travel
time was estimated for various seismic source positions
and compared with the travel time for a straight raypath
without the mine opening. Figure 163 illustrates the per-
centage difference in travel time between these two paths
as a function of position. The largest discrepancies
(delays) occur for seismic sources positioned immediately
adjacent to mine openings. Note the reduced error as-
sociated with the more rounded opening. In practice,
mine openings often fall between these two extreme
geometries.

Since microseismic events associated with mining
activity often occur within a very few mine-opening-
dimensions of a mine opening, the travel times measured
by source-receiver geometries as shown above are sys-
tematically biased to longer times than those associated
with the straight line propagation path assumed in the
event location process. For events occurring on the same
level as the mine openings, where most receiver stations
are typically located (e.g., fig. 4), there can be a large
proportion of the total number of the received signals that
experience interaction with mine openings. The same is
true for events occurring near the working stope opening.
As the propagation path becomes very long in proportion
to the smallest mine-opening dimensions, however, the
travel-time delay due to interaction with the mine openings
diminishes.

Three additional sources of systematic bias exist for
taypaths interacting with mine openings, most of which do
not diminish rapidly with increasing distance between the
source and the opening. First, there is a reduction in
amplitude of received signals due to (1) a factor of 5 re-
duction in amplitude response in the shadow zone region
of the mine opening (source placed on right-hand side of
figure 164) (25) for typical wavelengths associated with
observed microseismic events, and (2) the amplitude re-
sponse of the rib-mounted accelerometers is negligible at
angles of incidence close to 90°. Second, there is an in-
crease in the seismic wave rise time of approximately 40%
for the shadow zone position compared with an accelerom-
eter mounted on the wave-incident side of the mine open-
ing (25). Third, a blasting-related low-velocity damage
zone surrounds mine openings. Additional shadowing is
produced when raypaths encounter this low-velocity region
and suffer a further reduction in wave amplitudes. The
reduced amplitudes and slower rise time can provide a
systematic late arrival-time pick. The low-velocity zone
also yields longer travel times than would be found for
paths far removed from mine openings due either to travel
through the slower velocity structure or through a faster,
but longer path. The net result of all the above influences

25

is to further increase the travel time over the straight-
taypath travel time through a constant velocity medium.

On the basis of the two-dimensional calculations pre-
sented in figure 16B, one might expect that an event loca-
tion, calculated assuming an isotropic velocity, would be
displaced away from the mine opening. The actual picture
is much more complicated as illustrated in the following
example, Filled triangles in figure 17 show 12 arbitrarily
selected event locations near mine openings with vertical
coordinates positioned at the mine opening center. Travel
times to each receiver were estimated assuming a (short
wavelength) raypath which travels around the mine open-
ing along the fastest path; no other travel-time bias asso-
ciated with mine openings was included. The correspond-
ing calculated event locations are shown as plus symbols.
Some events are displaced away from, and some toward,
mine openings; most discrepancy vectors (linking actual
and calculated positions) have a similar apparent trend. In
this particular example, which was chosen to maximize the
influence of mine openings, the average three-dimensional
discrepancy vector is on the order of 1.2 m. The relative
positions of the event location, the mine openings, and the
nearest receivers with delayed signals determine the direc-
tion of mislocation.

Seismic Velocity Measurements

Seismic wave velocities were calculated using the
calibration-blast test data. The apparent velocity was ob-
tained from the ratio of wave propagation path lengths,
measured along a straight line between the blast site and
each accelerometer, to the measured P-wave travel times.

——
NS

Figure 17.Influence of mine-opening-related travel-time
delays on calculated event locations. Defined event locations
(filled triangles); event locations calculated with estimated
travel-time delays (+); nearby receivers (filled squares).



INTRODUCTION

Deformation in fractured and faulted hard rock is ac-
commodated to some degree by normal and shear dis-
placements along discontinuity interfaces. A deterministic
approach to characterizing the mechanics of discontinuous
geologic structures requires (1) specification of the lo-
cations of major fractures and faults, (2) characterization
of their interface (constitutive) properties, (3) specifica-
tion of the discontinuity-interface and other boundary con-
ditions, and (4) a numerical method for solving specific
problems. As part of its effort to gain insight into how
mechanical instabilities develop in underground mines and
result in rock bursts, the U.S. Bureau of Mines is devel-
oping such methods for quantitative characterization of
discontinuous rock masses (1)? The ability to predict the
mechanical response of discontinuous rock masses to
stress, whether originating in gravity, tectonic, thermal,
residual, or applied-induced forces (e.g, mining), is of
fundamental importance to scientists and engineers in-
volved in a wide variety of disciplines.

This report addresses one aspect of the mechanics char-
acterization problem, namely that of locating geological
discontinuities in areas that are not directly accessible for
visual examination. Visual observations made from within
mine openings can provide much direct information on the
locations of significant discontinuity structures. However,
preexisting fractures and faults often terminate abruptly
and generally cannot be extrapolated beyond mine open-
ings with a high degree of confidence. This is particularly
true in the complex geology of the Coeur d’Alene mining
district of northern Idaho where ubiquitous discontinuity
structures are observed over a wide range in size scales, In
addition, mere visual detection of a fracture or fault does
not necessarily provide any information on the influence of
the discontinuity on the response of the mine rock to
mining-induced stress.

One method, which has the potential to overcome some
of these problems, is the use of three-dimensional plots of
microseismic event locations in mapping active movements
along fractures and faults. An example of a significant
fault structure that appears to be delineated by the spatial
distribution of microseismicity is shown in figure 1. To
adequately resolve similar and smaller discontinuities using
this method, it is necessary to locate microseismic events
accurately and with high precision. The magnitude of the
typical location error sets a lower limit to the dimensions
of structures that can be resolved using spatial distribu-
tions of microseismic activity. For mine-stability analyses
related to rock bursting, the uncertainty in event location
coordinates should be less than the smallest mine and

Italic numbers in parentheses refer to items in the list of references
preceding the appendixes at the end of this report,

geological structures playing important roles in the rock
burst process, Accuracy and precision of better than 1m
would be desirable for analyses of small microseismic
events in Coeur d’Alene mines,

A variety of receiver array geometries are typically used
to study mining-induced seismicity. The arrays range from
surface-deployed low-frequency seismic networks to moni-
tor large (local Richter magnitude, M,, ~5.5) rock bursts
occurring throughout a mining district to underground
high-frequency arrays with interstation spacing of tens of
meters to monitor microseismic events (down to M, ~-6)
in individual stopes (2). A single array of several tens of
receivers, which is deployed to cover an entire mine, is an
intermediate-scale network. Few analyses of errors in seis-
mic event locations from any of these networks have been
published. In certain applications, such as the Bureau’s
rock burst research, this issue is critical.

The motivation for this study stems from observations
(3) made during the first underground test of a computer-
automated-measurement-and-control (CAMAC) digital
microseismic data acquisition system (J). In this test, at-
tempts were made to compare the locations of 256 micro-
seismic events with known geologic structures. The test
site geologic structures are typical of Coeur d’Alene mines,
As an example from the test area, figure 2 shows two
faults (among others) intersecting near mine openings,
possibly forming key block structures (4). Accelerometer

x
x

x LEGEND

mz

Raise
Fault

mmm Vein

0 5
Lo

‘ Seale, m

Figure 1.~Spatial distribution of microseismic events (X)
delineating & significant and actively adjusting fault (plan view).
Microselsmic events occurred in response to nearby mining
activity.


ACCURACY AND PRECISION OF MICROSEISMIC EVENT
LOCATIONS IN ROCK BURST RESEARCH STUDIES

By P. L. Swanson,' L. H. Estey,’ F. M. Boler,’ and S. Billington!

ABSTRACT

Stability analyses of fractured and faulted rock masses require delineation of the position, extent, and
orientation of geologic discontinuities. The size of the smallest active discontinuity that may potentially
be resolved using the spatial distribution of microseismic event locations is limited by the accuracy and
precision of the location methods. ‘ At a hard-rock mine in the Coeur D’Alene mining district of
northern Idaho, two data sets consisting of calibration blast signals from a known source site and origin
time and microseismic event signals were recorded using a stope-wide accelerometer array. These
seismic signals are used to quantify various sources of error in event location, Five factors influencing
source location errors are examined in this U.S. Bureau of Mines study: (1) biases of the numerical
source location techniques, (2) receiver array geometries, (3) uncertainties in receiver positions,
(4) errors in picking arrival times, and (5) uncertainties in seismic velocity structure, including the effect
of mine openings. In addition, synthetic data (accelerometer positions, travel-time picks and wave
velocity) are used to determine the effect of known systematic and random errors on source location
calculations. It is shown that the commonly accepted association of minimum travel-time residuals with
the best location solution does not necessarily hold true when there is a systematic error in seismic
velocity. Recommendations are made for increasing the accuracy and precision of locations of
microseismic events detected under similar field conditions.

‘\Geophysicist, Denver Research Center, U.S,-Bureau of Mines, Denver, CO.


14

by Godson (12), i.c., m-1 pairs ordered by arrival times.
Apparently, by ordering the arrivals and selecting time-
ordered pairs, the terms of equation A-1 involving the
velocity in the SW-GBM basis are kept at a minimum and,
in this way, impact the total solution the least. Moreover,
this subset of the SW-GBM basis also does reasonably
well when the errors are random. (Note: Even though the
m-2 BLD subset of staggered triplets with arrival-time
ordering gives better results with only systematic velocity
errors, this subset does not perform very well with random
errors.)

Thus, if the dominant source of error is random, the
iterative gradient bases give the best solution, the spatial
gradient basis being slightly better than the full gradient
basis. If the dominant source of error is a systematic un-
certainty of the velocity, a special subset of the SW-GBM
basis functions may give the best solution; this has been
selected as the preferred location technique.

For the rest of this report, sources are located using
the SW-GBM basis with ordered arrivals (m-1 equations)
to get a quantitative assessment of the influence of the
measurement errors on event locations.

INFLUENCE OF RECEIVER ARRAY GEOMETRY
ON EVENT LOCATION ERRORS

The importance of installing a receiver array that pro-
vides thorough coverage of the volume to be monitored is
usually appreciated by those monitoring mining-induced
seismicity, but difficult to achieve in practice. Receiver
positions are typically constrained by available access,
Poor array geometries, and/or poor coverage by adequate
geometries, can greatly amplify the effect of small meas-
urement errors (e.g., Ge (10)) and lead to poorly con-
strained locations.

A graphical event location technique that illustrates the
way in which array geometry affects the sensitivity of event
locations to measurement errors has been used by the au-
thors in analyzing the degree of constraint placed on event
locations by receiver array geometries (27). In an iso-
tropic velocity medium, three-dimensional surfaces of con-
stant relative arrival time, or isochrons, are constructed
for each pair of receivers reporting first arrivals (e.g,
curve AB for receiver pair AB in figure 84). This three-
dimensional hyperboloid surface, which is symmetric about
the line connecting the receivers, is a solution to an equa-
tion formed by subtracting P-wave arrival times for the
receiver pair in question using equation 1, It represents all

possible positions of a source that are consistent with the
measured difference in P-wave arrival time for the receiver
pair. The infinite number of possible solutions for the
source location on this surface is reduced by considering
additional relative-arrival-time isochron surfaces from
other receiver pairs. The mathematical solution for the
event location can be graphically interpreted as the com-
mon point of intersection of all of the isochron surfaces.
In some situations, multiple solutions are possible.

When errors in travel time are added to the exact
relative arrival-time data, the isochron surfaces do not
intersect at a point, but instead cluster throughout a vol-
ume‘ The size, shape, and orientation of this volume de-
pends not only on the magnitude of the errors but also on
the angles of isochron surface intersection; the geometry
of the reporting receivers determines these angles. If the
errors are normally distributed, the isochron intersection
volume is characterized by the SVD error ellipsoid’ The
lengths of its axes are inversely proportional to the amount
of constraint on the solution in that direction and thus pro-
vide a visual indication of the influence of array geometry
on possible event location errors. If the SVD ellipsoid
is equidimensional, the center of the intersection volume
(fig. 84) does not significantly change with small adjust-
ments in isochron position due to small measurement
errors.’ However, if the SVD error ellipsoid has one
dominant axis (fig. 8B), then small measurement errors
(eg. receiver pair AC) can produce large displacements
of the intersection point.

Maximum constraint on the solution can be achieved by
minimizing measurement errors and selecting an optimum
receiver geometry such that the number of isochrons that
intersect at small angles is minimized. This is achieved
by completely surrounding the events to be located with
a three-dimensional distribution of receivers. Linear and
planar array geometries for three-dimensional locations

‘Errors in travel-time pick, velocity along the raypath to the receiver,
and the station coordinates all result in uncertainty in the position of a
relative-arrival-time isochron.

various numerical location techniques evaluated in the previous
section were concerned with determining which kind of isochron surfaces
were to be used in the location process. For a given set of arrival times,
each basis function produces slightly different volumes of intersecting
solution surfaces. The different numerical solution algorithms (QR,
SVD, ete.) define a best location point within these volumes.
errors may completely change the position, shape, and ori-
entation of the error ellipsoid and may also result in nonintersecting
isochron surfaces.


Analog
system

Tape-and-map

tae First-arrival times
measured positions

from hardware

circuitry BLD method
Acceleroneters O— PU Event
———— location
Digital First-arrival times SW-GBM method
5 waveform from visual examination
posits acquisition of waveforms

Figure 3.—Main features of analog and digital microseismic monitoring systems. CPU, central processing unit; BLD, Blake; and SW-
GBM, Salamon-Wiebols-Godson-Bridges-McKavanagh source location.

In the digital system, data acquisition is based on
CAMAC? (8). Transient accelerometer waveforms from
a stope-wide array are captured on up to 32 12-bit analog-
to-digital conversion (ADC) channels that are digitized at
rates as high as 100 kilosamples per second and trans-
ferred via GPIB interface to workstation disk storage.
Compressional-wave first-arrival times are determined
using a graphics workstation which allows the display
of the stored waveforms at any desired magnification.
Accelerometer positions are determined using electronic
distance-measuring underground surveying techniques
where possible. Event locations are calculated using
the Salamon-Wiebols-Godson-Bridges-McKavanagh
(SW-GBM) method, described later in the report.

ARRAY DESCRIPTION

The set of high-frequency accelerometers from which
the arrays for both systems were drawn was distributed in
a volume approximately 100 by 100 by 70 m on a side, cen-
tered approximately 1.4 km below Earth’s surface at the

SCAMAC is an international digital interface and modular instru-
mentation standard (IEEE-583) that readily allows the combination of
high-speed transient recorders for acoustic emission, medium-speed
transient recorders for microseismic-seismic emission, multiple channels
of arbitrary voltage measurement for rock mechanics instrumentation
and any other off-the-shelf or custom-designed instrumentation to be
interfaced to generic computers through one or more interfaces.

4300 level of the mine (fig. 4). Thirteen of the receivers
(340 mV/g, frequency response up to 5 kHz) were mount-
ed to competent rock surfaces in mine openings with
epoxy (rib mounted). An additional 15 accelerometers
were mounted triaxially at five locations. Four triaxial
borehole probes were installed near the bottom of four
NX-size diamond drill holes, 4 m in length. One triaxial
unit was rib mounted. Vertical accelerometer coverage
was less than ideal, as is typical in a routine monitoring
array; all but four accelerometers were located within
+4 m of the 4300 level. Amplifiers close to the acceler-
ometers were used to increase the signal level by a mini-
mum of 40 dB, The amplified signals were brought via
cable to an underground instrumentation room 50 m from
the center of the array where both the digital and analog
data acquisition hardware were located.

CALIBRATION BLAST TESTS

Seismic source location accuracy is often judged by how
well the position of a controlled active seismic source, such
as a high explosive, can be located. In such tests, it is im-
portant to place the sources close to the zone of expected
microseismicity, because location accuracy is known to be
a function of the array geometry and position within a
given array (10). A measure of accuracy obtained from
controlled source location tests performed under ideal con-
ditions (i.¢., source placed in the center of the array,


880409.043331

arena

2000 2150 2300 2450 2600

2750 2900 3050 3200 3350 35.00 ms

Figure 11.—Example of comparison between manual (first vertical bar on accelerometer waveform) and hardware (second vertical
bar on waveform) picks of first arrival times shows typical observed discrepancies.

Figure 11 shows examples of the picks for the two
methods. It is clear that the hardware picker is unable to
do a good job on the pick at station 6. Station 6 is
thought to be bonded to a semi-detached block of rock
which provided systematically emergent signals. Station 6
happens to be the first hit station for most of the events of
the side-by-side monitoring test: the recorded arrival
times at the rest of the stations in the analog system are

systematically offset by the error in the first hit station.
This leads to significant location errors when the first hit
station data is used to form all of the equations for the
BLD basis.

Examples of the hardware picking a later phase, a fairly
common occurrence, are evident in figure 12. This event
was subsequently rejected through the software screening
in the analog system.


RI| 9395 REPORT OF INVESTIGATIONS /1992

Accuracy and Precision of Microseismic Event
Locations in Rock Burst Research Studies

By P. L. Swanson, L. H. Estey, F. M. Boler, and S. Billington

UNITED STATES DEPARTMENT OF THE INTERIOR

BUREAU OF MINES



«
< «
x”
—)
0 10
a
Scale, m
B °
ra

Figure 6.—Comparison between surveyed blast position (filied
triangle) and calculated locations using analog system (open cir-
cles) and digital system (+ symbols). Diamonds show nearby ac-
celerometer positions. A, Plan view; B, vertical section in the
plane of the vein, looking northwest.

a particular working area becomes active and approximate-
ly where in the area the activity is located. Also, ob-
servations of the rates of activity in a given area are not
influenced by these location uncertainties. What cannot be
determined, in this particular part of the stope, is a precise
correlation between microseismic events and a previously
identified structural feature, such as a fault, smaller than
10 to 20 m in size. Also, in this part of the array, pre-
viously unidentified microseismically active features cannot
be spatially defined to an accuracy of better than +10 m.

MICROSEISMIC MONITORING

The performances of the analog and digital systems
were evaluated using microseismic data recorded during a
week in which both systems were operating side-by-side.
Usually, accelerometer signals were derived from 12 of the
uniaxial stations described previously. In the analog sys-
tem, hardware-determined arrival times are first screened
in an attempt to eliminate extraneous picks and picks
associated with late-arriving shear waves and/or arrivals
from multiple events. Following the calculation of the

Figure 7—Comparison of 52 microselsmic event locations
determined using analog and digital systems. Calculated
locations for each event are connected by lines. A, Plan view; B,
vertical section through vein looking northwest.

event location, additional processing ensures that the
arrival-time sequence expected on the basis of the calcu-
lated location, is consistent with the observed sequence of
signal detection; otherwise the event is rejected. During
the week of monitoring, the analog system was triggered
by 2,070 events, of which 772 were well-located by the cri-
teria just described.

The digital system, using different triggering logic and
sensitivity, recorded 266 events. Arrival times at each ac-
celerometer were picked by visual inspection and assigned
a quality weighting factor ranging from 0 to 4, Acceptable
quality locations for 242 events were obtained. Acceptable
quality, in this usage, means that each of the several nu-
merical solution techniques tested (described later) gener-
ated locations with relatively low root-mean-square (rms)
travel-time residual values. The overlap between the ac-
ceptable quality events of the analog and digital systems
consisted of 52 events.

A comparison of the locations of these 52 events is
shown in plan and vertical section views in figure 7. Most
of the 52 event locations for the two systems are different


890107

100445

trig

Bih
Bix
rl2
ril

ria
9
8
7

333

rs
5
r4
3
r2
rl

3

rata NAV AAS
fat fat diem
PANTY

70.00 75.00 80,00 ms

Figure 5.—Example of calibration-blast waveforms and trigger signal (top).

COMPARISON OF SEISMIC SOURCE LOCATIONS

CALIBRATION BLAST TESTS

A series of test blasts was made by repeated detonation
in a single drill hole. Figures 64 and 68 illustrate the
calculated blast locations in plan and section views using
both the analog (circles) and digital (+) data in compari-
son to the surveyed blast position (filled triangle).

Two conclusions pertinent to this particular array and
site follow from the data presented in figure 6.

1. With either the analog or digital systems, the blast
location precision is approximately +1 m with an accuracy
of +10 m, Average mislocations for analog and digital
systems are 10 and 8 m, respectively.

2. The discrepancy between the actual and calculated
blast locations is sufficiently small that it does not detract
from the demonstrated utility of routine microseismic
monitoring using a stope-level array with the present anal-
ysis methods. A determination can readily be made when


38

all m(m-1)(m-2)/6 equations being independent if all have
different weights. Thus for m = 4 and equal weighting,
there are (m-1)(m-2)/2 = three independent equations
from which to form the m(m-1)(m-2)/6 = four equations
of the full system. For m = 5 and equal weighting, there
are six independent equations; for example,

346
3

yew
*

oo

represents one possible selection.

By using different subsets of the BLD basis function,
the number of different possible source locations is even
larger than with the SW-GBM basis. If all the BLD basis
functions are weighted independently, resulting in
m(m-1)(m-2)/6 independent basis equations, then there
are in general

m(m=1)(m-2)/6

[nejen-},

re3

possible different solutions.


)

Figure 18 shows in plan view the apparent velocities meas-
ured along unobstructed raypaths. (For the other stations,
some part of a hypothetical straight raypath would have
traversed a drift or a mined-out stope.) The values shown
are average values from five tests. The repeatability of the
travel-time measurements in the calibration blast tests was
typically 1.5%.

The apparent P-wave velocity varies by as much as 50%
along the different unobstructed raypaths. The lowest ve-
locity measurement (3,260 m/s or 10,900 ft/s) occurs over
the shortest path, which is entirely near the active stope.
The highest velocity (4,950 m/s) occurs over one of the
longest paths which travels through a greater proportion of
undisturbed material far removed from mine openings.
Additional velocity measurements made over different
paths between the stopes of figure 4 using hammer sources
ranged from 4,100 to 5,700 m/s. The lowest values of
velocity were observed over paths most affected by mining-
induced fracturing near mine openings.

The travel path from source to receiver will not be a
straight line, but a curved path when there are such large
spatial variations in seismic velocity (>15%, 26). The lo-
cation calculations assume simple straight raypaths, Thus
the constant seismic velocity model is only an approxima-
tion to the actual medium." In accepting the use of a
constant velocity model, two questions must be posed:
(1) which value of velocity should be used; and (2) how
does this simplification of the velocity structure affect the

MUThe relative amounts of spatial variation in velocity due to geologic
iti ining-induced fracturing are now being investi-
sional seismic tomography experiments at this site.

10

Scale, m

accuracy of the location? This latter question is deferred
to the discussion section,

An average velocity was determined by plotting the
measured travel time against the straight-line distance be-
tween the calibration blast site and each receiver (fig. 19).
For raypaths interacting with mine openings, a small path-
dependent adjustment has been applied to the straight-ray
Propagation distance. The slope of a straight-line fit to the
data, through the origin, is taken to represent the average
velocity. At the 95% confidence level, this average value
of velocity is 5,000 +200 m/s. It must be kept in mind
that since the velocity varies as a function of position
(fig. 18) there is no single correct value of average velocity.
An appropriate average value for a source located at the
calibration blast site is 5,000 m/s. Significant location er-
rors will result for any one average value of velocity since
raypaths will in general traverse mine rock with seismic
velocities varying over the observed range (up to 50%).

Location discrepancies for the analog and digital sys-
tems due to systematic uncertainties in the average value
of velocity are now examined. Figure 20 shows the control
volume with 1,000 defined event locations and event loca-
tions calculated using the SW-GBM method with a 4%
systematic error in velocity (equivalent to the 200 m/s un-
certainty determined from the data in figure 19). The av-
erage magnitude of the resulting three-dimensional dis-
crepancy vectors is 1.2 m, This figure represents the lower

i 1
50 100
DISTANCE, m

150

Figure 18.—Apparent selamic velocities (kllometer per second)
calculated for straight raypaths in calibration blast tests. Filled
trlangle Indicates calibration blast location. Filled squares are
receiver locations.

Figure 19.—Measured travel times and straight-line distances
between all source and receiver sites for calibration blast tests.
Dashed line represents 5,640 m/s (13% higher than 5,000 m/s),
a velocity that may be more appropriate for unfractured rock far
removed from mine openings.



Since the first signal detected at the receiver travels
along a minimum-time travel path (Fermat’s principle), the
fastest path, in the short wavelength limit, is through the
solid rock medium around the mine opening. This travel
time was estimated for various seismic source positions
and compared with the travel time for a straight raypath
without the mine opening. Figure 163 illustrates the per-
centage difference in travel time between these two paths
as a function of position. The largest discrepancies
(delays) occur for seismic sources positioned immediately
adjacent to mine openings. Note the reduced error as-
sociated with the more rounded opening. In practice,
mine openings often fall between these two extreme
geometries.

Since microseismic events associated with mining
activity often occur within a very few mine-opening-
dimensions of a mine opening, the travel times measured
by source-receiver geometries as shown above are sys-
tematically biased to longer times than those associated
with the straight line propagation path assumed in the
event location process. For events occurring on the same
level as the mine openings, where most receiver stations
are typically located (e.g., fig. 4), there can be a large
proportion of the total number of the received signals that
experience interaction with mine openings. The same is
true for events occurring near the working stope opening.
As the propagation path becomes very long in proportion
to the smallest mine-opening dimensions, however, the
travel-time delay due to interaction with the mine openings
diminishes.

Three additional sources of systematic bias exist for
taypaths interacting with mine openings, most of which do
not diminish rapidly with increasing distance between the
source and the opening. First, there is a reduction in
amplitude of received signals due to (1) a factor of 5 re-
duction in amplitude response in the shadow zone region
of the mine opening (source placed on right-hand side of
figure 164) (25) for typical wavelengths associated with
observed microseismic events, and (2) the amplitude re-
sponse of the rib-mounted accelerometers is negligible at
angles of incidence close to 90°. Second, there is an in-
crease in the seismic wave rise time of approximately 40%
for the shadow zone position compared with an accelerom-
eter mounted on the wave-incident side of the mine open-
ing (25). Third, a blasting-related low-velocity damage
zone surrounds mine openings. Additional shadowing is
produced when raypaths encounter this low-velocity region
and suffer a further reduction in wave amplitudes. The
reduced amplitudes and slower rise time can provide a
systematic late arrival-time pick. The low-velocity zone
also yields longer travel times than would be found for
paths far removed from mine openings due either to travel
through the slower velocity structure or through a faster,
but longer path. The net result of all the above influences

25

is to further increase the travel time over the straight-
taypath travel time through a constant velocity medium.

On the basis of the two-dimensional calculations pre-
sented in figure 16B, one might expect that an event loca-
tion, calculated assuming an isotropic velocity, would be
displaced away from the mine opening. The actual picture
is much more complicated as illustrated in the following
example, Filled triangles in figure 17 show 12 arbitrarily
selected event locations near mine openings with vertical
coordinates positioned at the mine opening center. Travel
times to each receiver were estimated assuming a (short
wavelength) raypath which travels around the mine open-
ing along the fastest path; no other travel-time bias asso-
ciated with mine openings was included. The correspond-
ing calculated event locations are shown as plus symbols.
Some events are displaced away from, and some toward,
mine openings; most discrepancy vectors (linking actual
and calculated positions) have a similar apparent trend. In
this particular example, which was chosen to maximize the
influence of mine openings, the average three-dimensional
discrepancy vector is on the order of 1.2 m. The relative
positions of the event location, the mine openings, and the
nearest receivers with delayed signals determine the direc-
tion of mislocation.

Seismic Velocity Measurements

Seismic wave velocities were calculated using the
calibration-blast test data. The apparent velocity was ob-
tained from the ratio of wave propagation path lengths,
measured along a straight line between the blast site and
each accelerometer, to the measured P-wave travel times.

——
NS

Figure 17.Influence of mine-opening-related travel-time
delays on calculated event locations. Defined event locations
(filled triangles); event locations calculated with estimated
travel-time delays (+); nearby receivers (filled squares).



0 5

—
Scale, m

Figure 2.—Block structures formed by Intersection of faults
with other discontinuities and mine openings (plan view).

waveforms from microseismic activity were recorded in
this area and analyzed in terms of first-motion polarity
patterns for estimating rupture orientations and slip di-
tections (3). The consistency of first-motion patterns
indicated that slip was taking place along preferred planes
and that the preferred slip directions changed with time
and nearby blasting. Unfortunately, the boundaries of
these and other similar block structures, which must be
determined to evaluate the local mine stability, could not
be adequately resolved because of systematic and random

The report is organized in the following fashion. First,
the data acquisition and analysis methods for two different
microseismic monitoring systems deployed at the field site
are described. Second, estimates of the source location
accuracy are provided for each system by comparing calcu-
lated and measured calibration blast locations, The prob-
Jem of location accuracy is further illustrated in a compari-
son of microseismic event locations calculated using data
and analysis methods for each system. In the Error Analy-
sis section, several common source location techniques are
first examined and then one method is selected which min-
imizes bias introduced by numerical methods. The meth-
od which minimizes the bias is determined by comparing
defined synthetic event locations with locations calculated
using controlled amounts of error in travel time, receiver
position, and seismic velocity. The influence of receiver
array geometry on event locations through amplification of
measurement errors is then investigated. Next, the actual
or estimated errors in accelerometer position, seismic
travel time and velocity associated with the two micro-
seismic monitoring systems are used to calculate average
source mislocations for each source of error using random
distributions of 1,000 synthetic source sites, The com-
posite accuracy and precision of source locations are then
estimated for each system using the combined sources of
error, Recommendations are made for increasing the lo-
cation accuracy of microseismic events for similar arrays.
This work is directed at improving safety in the mining

microseismic event location errors. A primary goal of this industry, a primary goal of the Bureau,
report is to illustrate how these location errors arise.
ACKNOWLEDGMENTS

The authors thank ASARCO Inc. for their long-
standing cooperation with the Bureau in the field ex-
periments at the Galena Mine, Wallace, ID. The con-
siderable expertise and assistance brought to these field

experiments by C. D. Sines was essential for their suc-
cessful completion. J. Coughlin provided the authors with
many details of the hardware and software in the analog
microseismic monitoring system.

EXPERIMENTAL PROCEDURES

DATA ACQUISITION AND ANALYSIS METHODS

Two stope-level location systems are investigated.
The first is an analog system that is in use at the Galena
Mine in Wallace, ID, to characterize space-time statis-
tics of microseismicity. After field trials as a single-stope
monitoring system (5), it is now being expanded to an
industrial-scale computer networked mine-wide system
(6-7). The second system is a research-grade digital
system (8) deployed in a single stope and used in char-
acterizing deformation mechanics in fractured rock-
burst-prone ground.

Figure 3 summarizes the general features of each
location system. In the analog system, accelerometer
positions are estimated from available mine maps. Up to
32 signals from a stope-level accelerometer array are
amplified and fed into a hardware unit that determines the
relative arrival times of the first detected phases (5). This
information is then sent via an RS-232 line to a computer
in the instrumentation room for event location calculation
and display. The Blake (9) source location method (BLD)
is used with the assumption of a uniform isotropic seismic
wave velocity structure,


ACCELEROMETER POSITION MEASUREMENT

To describe the positions of receivers and microseismic
sources, a frame of reference is needed. Mines often es-
tablish survey control points on a system of east and north
and elevation coordinates, and produce drift outline maps
as workings proceed. Drift outline maps show the plan
view relative to east and north grid lines at one horizontal
level or in vertical section. The accuracies of such maps
vary tremendously depending upon the intended applica-
tion. In the study area, the mine survey and drift maps
were established approximately 30 years ago. As this area
has been subjected to continual deformation associated
with vertical stope mining under high horizontal stress
(sandfilled stope closures as great as 0.5 m (22)), including
tib sloughage and rock burst deformation, the relative po-
sitions of certain survey control points are expected to
have changed significantly since the original survey.

Existing (30-year old) mine maps provide the frame of
reference for accelerometer positions in the analog system.
The usual method is to find a recognizable feature on the
map that is close to the accelerometer in question; the dis-
tance to the accelerometer from this feature is then meas-
ured and plotted on the map. The vertical coordinate is
typically referenced to the track or floor level.

Several sources of error are associated with this
method, First, the accelerometer cannot be located to
better than the accuracy of the original mine map, even
with perfect identification of a position on the map.
Second, the map representation of a rib (a single line) is
only an approximation to a rough three-dimensional
surface. Third, changes in mine opening geometry occur
with time due to rib sloughage and other mining-induced
deformation such as rock bursting. Rock bursting also
produces ground deformation in excess of obvious surficial
modifications to mine openings. Overall, the uncertainty
for the accelerometer position determinations used by the
analog system is conservatively estimated at 1 m.

Control points of the 30-year old mine survey provide
the frame of reference for accelerometer position de-
termination in the digital system. Two surviving survey
control points (most of the survey spads have not survived
30 years of mining activity) were used to establish a
baseline to which all calculated accelerometer positions
were referenced. The coordinates of one point also serve
as the tie-in to the mine coordinate system. Accelerom-
eter positions are established by measuring angles and dis-
tances using a theodolite equipped with an electronic
(infrared) distance-measuring (EDM) device. Occasion-
ally, in difficult-to-survey or obstructed areas, short dis-
tances to accelerometers are estimated using tape, com-
pass, and inclinometer readings.

The overall uncertainty in accelerometer positions using
these techniques is estimated to be +5 cm. This location

17

uncertainty is referenced back to the particular mine sur-
vey spads used to establish the reference baseline and not
the original mine maps themselves. Thus, systematic dis-
crepancies may exist between the locations of accelerom-
eters plotted on these maps and the locations of mine
openings on the map due to systematic errors (1) in the
original mine map, (2) from subsequent deformation in the
surveyed area including displacement of survey control
points, and (3) in the accelerometer survey itself.

Table 3 shows the differences between accelerometer
position coordinates determined using the analog and dig-
ital systems. Differences in individual station coordinates
ranged from 0.02 to 1.9 m (0.07 to 6.21 ft). The average
three-dimensional discrepancy vector has a magnitude of
approximately 1 m, Systematic differences between the
two sets of measurements were minimized by finding a
best fit coordinate transformation (rotation and transla-
tion) between the two sets of accelerometer coordinates,
One is left with an estimate of the quasi-random differ-
ences between the two methods shown as residuals in the
bottom half of table 3. These differences are repre-
sentative of both the random measurement errors asso-
ciated with the two methods of determining accelerometer
position and deformation-related displacements occurring
since establishing the two reference frames. The individ-
ual station coordinate residuals range from 0.02 to 1.43 m
(0.08 to 4.7 ft), The average best fit three-dimensional dis-
crepancy vector has a magnitude of 0.79 m (2.6 ft).

Table 3.—Difference between surveyed and tape-
and-map-measured accelerometer coordinates
before and after transformation, feet

Station ax ay az lar]
BEFORE
1.82 0.09 -1,71 2.50
«72 1.42 “71 1.74
+34 4 08 64
272 196 1.10 3.83
1.43 189 0 -1.73 2.93
5.58 «78 07 5.63
1.70 621 -1.03 6.52
12 21.21 83 1.47

NAp NAp  NAp 13.42
AFTER

2.16 1.34 3.91
55-411 1.28

1.2200 +147 221
+56 1.29 2.81
~87 «30 4.70
3.84 62 4.35
71.20 6 1.97

2.61



33

APPENDIX A.—BASIS FUNCTIONS FOR SOURCE LOCATION

Here the generalized basis functions for determining
the source location are given along with an outline of their
derivation. In general, it is necessary to obtain both the
location (x, yo, and z,) and the origin time (t,) of the
source assuming a homogeneous and isotropic (spatially
constant) medium which transmits (either P or S) seismic
signals at a velocity v, Receiver stations having arrivals of
the signal are indexed 1,....m.

SW-GBM BASIS

To derive the SW-GBM basis function, the travel-time
equation 1! is squared:

5 ~~ eh 0x8 +8 - V5 - 29m
~ iyo ~ 2329 + 20% ty

Note that here, unlike the outlines given by Salamon and
Wiebols (12)? or Godson (12), the arrivals at receivers are
not necessarily ordered in any particular way and the un-
known time is cast explicitly in terms of the source origin
time, tp, rather than the travel time to receiver j = 1, t,.
Nonlinear terms in the above equation are removed by
subtracting equations of any two different receivers, say
the jth and kth. A system of linear equations can be ob-
tained using the resulting basis function.

Sot vet a B- VG - B= mq - x)
+ 2yol ~ Yi) + 2zolq - %) ~ Aov"(t; - ty. (Art)
BLD BASIS

The BLD basis function is obtained similarly in that
nonlinear equations of form equation 1 are found for two
different receivers, say, j and k, which are then squared
and differenced. Then, the origin time of the source (t,)
is replaced with the spatial coordinates of another re-
ceiver 1. Generalizing the development of Blake (9), the
starting equation for receiver j is, as stated, essentially
equation 1;

V(t ~ th + ~ te) = [05 - %)* + ( - Yo)?

+G- ny}? , (A-2)

4Bquation numbers without an A- prefix refer to equations in the
main text.

*ttalic numbers in parentheses refer to items in the list of references
preceding the appendix.

where t, - tg (i.¢., zero) has been added to the travel-time
difference from the source to receiver j, and of course tg
is the arrival time of the signal at receiver 1, The
substitutions

dig TM v(t, ~ te) = 10% - %)* + - Yo)? + - IM,
dig = 10% ~ %) + (% - Yo)? + - 2)",
dy = v(y - t), (43)

are then used where d,, 2 0 is the distance from receiver
1 to the source, d, 2 0 is the distance from receiver j to
the source, and |d,| is the distance the signal has traveled
between the time of arrival at receiver j and the time of
arrival at receiver 1 (note that d, can be either positive or
negative or zero). After substituting A-2 into A-3 and
squaring,

af +43 + 2ajdy = 43

The same steps are taken for receiver k and receiver 1,
yielding

(a4)

2 + 2 + dud = 4,2. (A-5)

Equation A-4 is then multiplied by d,, and equation A-5
is multiplied by dj, The difference of (A-4) * dy and
(A-5) * dy is taken which, after some simplification, gives
a basis function for a different system of linear equations
with the unknowns X, Yy 2:

2 2 2 2 2 2 2,
(Gj - 3) - yp ap exp typ + dy |

Gd -8-s- ded of exdg =
Adal - %) - du - x)
+ Adu - y) ~ dun - yD¥o
+ Adu - %) - diy - yo
which is a generalized version of Blake’s equation C-11,
After the spatial coordinates x), yp, and z have been found,
an estimate of the origin time t,, can easily be obtained as

follows. Estimates of ty are calculated for each receiver
i=1.,m

(A-6)

to = 5 1G 90? # = WP + @ PM, (AD

which are then averaged to obtain a value for ty.


FULL-GRADIENT BASIS

The basic premise for the gradient basis is that it has a
function of n unknowns (here n = 4, i.e. x, Yo, %, and ty)
with m constraints (here m = number of receivers with ar-
rival time information). If m 2 n, a solution is attempted.
The starting equation is again equation 1 which will be
called a function f() of the unknowns:

£(to Xo Yoo 20) = [0% - 0)? + (i - Yo)?
+@~ aI? ~ vt - t) = 0.
A Taylor series expansion can be used to obtain a linear
approximation to f() in terms of an approximate solution
(to's 0" Yo’ 0")
0 = f(to, Xo Yor 20) = fh" %9” Yo" » 20’)
~ try 2f ae ex) 2 ge
+ (to 0) ae to") + Go 10) a, M0")

r) aw, ~ at) 8f wy
* Go - W) Sr Go’) + 20) Far &')
where the partial derivatives are simply
of

— Vv

35
Fe 7 OIG + O19? + G-29"T?
Fe 7-1) Ue? + Gi-¥9? + age?

Fe 7 OD Gad? + Gir? + G29"?

(One could also take the derivative with respect to
the velocity v—in which case there would be n = 5
unknowns—and hence attempt to find corrections to the
starting velocity.) Thus a system of m linear equations is
constructed from which an attempt is made to find the dif-
ference or gradient vector, & between. the improved solu-
tion, X, and the approximate solution, x,’

v(t, ~ 0) - di] = divat + (x ~ x) ax
* Wo ~ yay + % -2)Az, (48)

where At, = ty- ty’ , Ax = %-%" , Ay = Yo-Yo" » Azy
=h-h!.

When di, is found, it is then added to % yielding the im-
proved solution %. This improved solution %’ then be-
comes a new X,’ from which a new gradient vector is
found and so on iteratively until the gradient vector is
essentially zero. At this point X/ and X, are essentially
equal, a local minimum of equation 1 has been found, and
the microseismic source is at to, X Yor Zo

SPATIAL-GRADIENT BASIS
The derivation of the spatial-gradient basis function fol-

lows that given in Lienert (16), Here equation 1 is first
rewritten as a function of the travel-time residuals

aT;
t - ty - Ti(%y,z) = Aty + Ax x
ay oT gy Th
+ byy — + Az—t
Yo oy * oz"

where T, (x, y, z) are the calculated travel times from the
source to the ith receiver. Since the source location and
the source origin time are interrelated through equation 1,
the full-gradient vector can be decomposed into two sep-
arate gradients: one is three-dimensional—the spatial
gradient—and increments the estimate of the source loca-
tion; the other is one-dimensional and increments the es-
timate of the source origin time using the solution of each
spatial gradient. The one-dimensional gradient is

oT;
ty + Aty = <tj> - <Tj> - Ax< >
x

4 <a, A «oh, (A-9)
Say? OSH?
where the brackets, < >, represent weighted means in gen-
eral. This equation is used to get an update of the estim-
ate of the origin time, i.e., ty + At, for each increment of
the spatial gradient, ic, (Ax, Ay, Az). The basis
function for the spatial gradient is

t; <t> aT; oT] a
BTS ie a7]

Tm a),
—-<

ay ay ‘Yo
am eT),
—-<

az az 0
where n = 3 since the unknowns involve only the correc-
tions to the spatial coordinates of the source.

(A-10)


T T T T T

100 F 4
e
e
N 10> 4
e
e
e

tr e 0@

i 1 L L 1

0 0.2 0.4 0.6 0.8 1.0 1.2
HARDWARE-P1CK ERROR, ms
Figure 13.—Histogram showing difference between hardware

and manual travel-time picks, Straight line indicates exponential
fit to errors less than 0.8 ms.

To allow comparison of travel-time pick errors with
discrepancies resulting from accelerometer position meas-
urement errors, figure 15 shows another defined set of
1,000 events in the control volume. Event locations cal-
culated using exact travel-time data with controlled ran-
dom errors added are shown with connecting lines. The
errors were randomly selected subject to the constraint of
the observed exponential distribution of travel-time-pick
errors. The average three-dimensional discrepancy vector
for the events in this control volume has a magnitude of
28 m. Location calculations using error values appropri-
ate to manual first-arrival picking resulted in average
three-dimensional discrepancy vector magnitudes of 0.7 m.

\

Figure 14.—Comparison between 40 event locations calculated
using hardware and manual picks of arrival times with all other
parameters Identical, Lines connect each set of solutions. A,
Pian view; B, vertical section looking northwest.

SEISMIC WAVE VELOCITY
influence of Mine Openings

Mine openings represent a distinct source of heteroge-
neity in the seismic velocity structure. The seismic velocity
in air (330 m/s) is so low that most of the seismic energy
travels around the periphery of the opening. As a prac-
tical means of assessing the influence of openings on mi-
croseismic event locations, the measured travel times were
investigated (as opposed to the much more involved task
of calculating accurate raypaths through heterogeneous



CONTENTS

Abstract ....
Introduction .
Acknowledgments .
Experimental procedures ..
Data acquisition and analysis met
Array description ......
Calibration blast tests . . .
Comparison of seismic source locations . . .
Calibration blast tests .....
Microseismic monitoring .
Error analysis
Numerical source location t
Basis functions ..........
Solution of linear systems
Numerical evaluation ..,
Summary and recommendations ..
Influence of receiver array geometry on event locati
Accelerometer position measurement ....
Arrival time determination .
Hardware timing ..
Manual timing ..
Software timing .
Example comparison .
Seismic wave velocity ...
Influence of mine openings . .
Seismic velocity measurements ...
Discussion ....
Conclusions .
References
Appendix A.—Basis functions for source location
Appendix B.—Subsets of direct basis functions . .
Appendix C—Examples of source location process

ILLUSTRATIONS

Spatial distribution of microseismic events delineating a significant and actively adjusting fault . .
Block structures formed by intersection of faults with other discontinuities and mine openings
Main features of analog and digital microseismic monitoring systems . .
Plan view of uniaxial and triaxial receiver array .
Example of calibration-blast waveforms and trigger signal . vee
Comparison between surveyed blast position and calculated locations using analog system and digital
SYSLEM eee eee eee ee eee e Er NEED EEE EEE E EEE EEE EEE EEE EEE EEE EE EE EEE
Comparison of 52 microseismic event locations determined using analog and digital systems . .
Seismic event locations using intersections of relative-arrival-time isochrons .
Influence of receiver array geometry on event location errors .
Defined event locations selected at random in control volume centered on receiver array connected with
straight lines to event locations calculated with accelerometer position measurement errors .........
11. Example of comparison between manual and hardware picks of first arrival times shows typical observed
discrepancies... 6... 6. cee e eee e cece e tree eee e reteset eee ee eee tere ete eee ee EE
12. Example of hardware picks missing first phase and picking later phase ..........ssscseeeeeeeees

AwSSh

~
Beem


39

APPENDIX C.—EXAMPLES OF SOURCE LOCATION PROCESS

Explicit numerical examples of the source location
process are given for (1) the SW-GBM basis function with
the m-1 time-ordered GBM subset, (2) the spatial-gradient
basis function, and (3) the full-gradient function. Selected
for this numerical example is one of the calibration blasts
discussed in the text which had the following set of ob-
served arrival time readings, P(obs), for each corres-
ponding station whose coordinates are given in meters:

station P(obs)(s)_x (m) y (m) 2 (m)
id 0.04508 3438.534 2793.316 -364.462
8 04506 3421.74 2816882 -346.750

44 04162 3400.516 2803.324 363.934
i] 04368 3417.924 2778723 -364.260
m5 .05730 3464.037 2856.791 -418.682
7 104526 3384.110 2815.730 -363.047
8 05114 3354.194 2816.320 -963.709
19.4 104316 3416.935 2816420 -362,691
no 05878 3414.929 2746.971  -272.683
n2 04732 340.619 2767.520 -363,983

The arrival time readings have a precision of +0.00002 s
and the coordinates have a precision of +0.005 m. The
P-wave velocity is assumed to be a constant 5.020 km/s.

The set of linear equations formed by the SW-GBM
basis function (appendix A) with the m-1 time-ordered
GBM subset (appendix B) is given below. Here, all
weighting for individual stations coordinates and arrival-
time readings are equal (ic., unity). The solution x has
been found using the QR (hence, “least-squares") algo-
rithm described in the text which employees partial
pivoting. The units used for this example are meters for
distance, milliseconds for time, and meters per millisecond
(equals kilometers per second) for velocity. Since there
are m = 10 arrival-time readings, 9 equations are formed.
See below. The source is located (using this set of equa-
tions) to be at x = 3,412.9 m, y, = 2,798.6 m, and
% = -362.7 m. The origin time is directly computed to be
0.03777 seconds (relative to the same time origin as the
arrival-time readings).

The SW-GBM solution does not, in general, yield a
very good direct estimate of the origin time of the source.
In this example with the origin time at 0.03777 s, the
travel-time residuals are computed to be

station —P(calo) —P(0-<)
12 0,042909 0.002081
3 042909 .002151
14.1 040424 001196
5 041876 .001804
7 044444 000816
8 049989 .001151
19.1 041405 .001755
110 058446 .000334
2 046078 001242
15 056811 .000489

where P(o-c) is the difference between the observed and
calculated arrival times, (P(calc)), for the P-wave, i.e., the
P-wave travel-time residuals. The effective rms of the
travel-time residuals is 0.001435 s normalized by 1/m, or
0.001853 s normalized by 1/(m-4), where 4 is the effective
number of unknowns in the problem (origin time and
three spatial coordinates). However, the origin time can
be computed independently from the estimate of the
source coordinates and the observed arrival times at each
receiver; for this example, the mean of the values for the
origin time (A-7) yields a revised estimate of the origin
time of the source to be 0.039074738 s, with the revised
travel-time residuals to be

station Peale) Poe)
m2 0.044301 0.000779
8 044211 000849
44 041725 0.000105
6 043178 000502
7 045746 +000486
6 4.051291 ~000151
19.1 042707 000453,
10059748 -,000968
n2 .047380 ~000060
ns 058113. +000813

The effective rms of the travel-time residuals is now
reduced to 0.000605 s normalized by 1/m, or 0.000781 s
normalized by 1/(m-4).

-77.617088 32,837891 26.191895 2.485962 181337.859375,
-26.208632 1.978027 -75,394043. 3.138000 +204157.828125
9.552963 7.699707 76.318359 35.020020 (224319.671875
1.008039 33,520020 ~47.132324 = -35.424011 ~4680.454102 7.772820
A= — -9.072160 -108.847656 44828125 2.830017 bb = -247032.125000 x=  3412,905762
-103.825714 — 113.017578 -96.419922 1.872009 112363.328125, 2798.638184
-192.530853  -172.849609 97.600098 0.547974 ~324426.406250 -362.666046
-310.469116 — 219.686035 80.941895 — -109.945984 1004708,000000
-74,592972 — -98.216309 _-219.640137 _291.997986 +1058486.750000


10

basis. Another gradient basis is cast directly only in terms
of the source location; the source origin time is in-
cremented separately. This basis will be referred to as the
spatial-gradient basis. These types of iterative basis
functions are often used in broad application seismological
location codes; e.g., a type of full-gradient basis is used in
HYPOINVERSE (14) whereas a type of spatial-gradient
basis is used in HYPO71 (15) and HYPOCENTER (16).
According to Professor M. Salamon, Colorado School of
Mines (personal communication), the spatial-gradient basis
function solved with an L2-norm minimization is equi-
valent to the iterative solution method used in Salamon
and Wiebols (11). Solutions using the spatial-gradient
basis function explicitly minimize the sum of the squares
of the travel-time residuals with respect to location
coordinates and origin time. The constant velocity
medium is a special case yielding a simple form for the
full-gradient basis equation A-8 and the spatial-gradient
basis equation A-9; more complex velocity structures are
also amenable to these bases, The gradient bases nat-
urally yield m independent equations, even when weighting
is employed,

Solution of Linear Systems

Any of the basis functions previously described can be
used, or some variation of them, to form a linear system
of equations which can be denoted by the matrix equation
Ax = b. A is the so-called design matrix, b is the data
parameter vector, and x is the unknown parameter vector.
Moreover, weighting, which can reflect the measurement
uncertainties of the coordinates of each receiver and/or
arrival times at each receiver, can easily be incorporated.

There are, of course, many different ways to obtain a
solution to the system Ax = b. Monte Carlo methods
search randomly and other brute force methods search
systematically (¢.g., G-block, (9)) for a solution, Other
numerical algorithms allow one to obtain more rigorous
estimates of the solution without random or systematic
searches.

If the system Ax = b is overdetermined, ic., the num-
ber of equations is greater than the number of unknowns
(which is generally the case for microseismic location
problems), numerical algorithms such as Choleski’s meth-
od, QR-decomposition, singular value decomposition
(SVD), and linear simplex methods can'be used to obtain
a solution, An algorithm such as Choleski’s method is
undesirable since it requires one to set up a new linear
problem that has the square of the condition number of
the original problem; with this new linear problem (the
so-called normal equations), twice as much floating point
precision is required to find the solution to the same

numerical precision as is possible with the original linear
problem (13, 17-18).

On the other hand, an algorithm such as SVD is desir-
able since with it one can get an error estimate for the so-
lution if the measurement errors are normally distributed.
This error estimate can either be in the form of a single
number representing the uncertainty of each unknown, or
in the form of an error ellipsoid which shows how the
error is distributed in solution space (13). This error
estimate is proportional to |Ax - b] (Euclidean L2 norm)
and can also be scaled to represent different confidence
limits of the uncertainty.

For any set of algorithms, there may be no unique solu-
tion x to Ax = b, especially when the system is overdeter-
mined. Worse yet, different code—even for the same gen-
eral type of algorithm (different SVD codes, say)—may
give different solutions since there are good ways and bad
ways of implementing a numerical algorithm.

Given all of the different possible basis functions (i.c.,
equations A-1, A-6, A-8, A-9, and others) and all of the
different algorithms for solving the system of equations re-
sulting from these bases (and different codes for the same
algorithm), what is the best numerical strategy for locating
a source? As a result of this investigation, it was dis-
covered that the best strategy seems to depend upon the
dominant type of error in the parameters of equation 1.

Numerical Evaluation

The performance of the different combinations of basis
functions and linear systems solution algorithms was in-
vestigated systematically in the following way. Two dif-
ferent receiver arrays were examined. One was based on
an actual array that had been set up to monitor micro-
seismic activity around a specific stope near the 4300 level
of the Galena Mine; this was a typical stope array largely
concentrated on one mine level with a few receivers lo-
cated off level, The other was a hypothetical array around
an adjacent stope defined using realistic constraints found
in the research site (i.¢., locations of accessible mine open-
ings, minimum number of borehole sites, etc.) in an at-
tempt to maximize location resolution in all spatial di-
rections inside the array.

A set of 1,000 possible event locations were chosen at
random inside a cubic volume 90 m on a side roughly cen-
tered on each receiver array and roughly the same size as
the array. Two sets of arrival times were calculated, One
set was the perfectly known arrival times for these event
locations calculated for all receivers assuming a perfect-
ly known constant velocity medium and perfectly known
event location and receiver coordinates. Another set of
arrival times was produced by adding random errors (using


2

Table 1.—Performance comparison between numerical source location methods, hypothetical Idealized array

Random errors’
Method 0.000020-8 seed (0.000080-s seed 0.10-m seed 1.00-m seed
<)at]>,8__<]ax]>,m_<|dt]>,s__<Jdx]>,m_<Jat[>,s8_<|dx[>,m_<]dt}>,8 <|dxj>,m

0.000048 0.321587 0,000207 1.286343 0.000013 0.100008 0.000118 0.913274
000047 310034 000203 1.242515 000012 poy 000112 871034
000051 348318 000217 1.399101 000012 093599 888562
000054 .873019 000227 1.484210 000015 191051
000055, «379864 000230 1.517062 000014 120781
-000054 (874245 000228 1.491588 000015 111640
.000057 393101 000238 1.567563 000017 137835

. 000048 821577 000207 000013 +100009
. 000047 310039 000203 095408
. 000051 348302 000217

FIS, unordered

000012, 093589

(m-2)BLD .. 200 000051 1328467 ‘000210 7000883 °4.573452
[000018 135734]

(m-2)BLD, unordered ...  .000071 490502000289 1000016 126929
My SAG) «scree tse 90008 -839688 C0084 000024 200586
m, ag yordered 1000064 456368 000020 “161289
(rn-2) 000070 ‘08479 .000285 0000171139787
Im2\stag’, unordered’ :. “Oooo “aeess = “b00286 1000019 1156417
Spatial racient 1900033 282034 += .000133 000009 072413
Full gradient . ++ 1000042 ___-295770___.000168__1:188689__000012__090023

Systematic errors”
4.8-km/s velocity 5.2-km/s velocity 5.64-km/s velocity
<|dt{>,s__<|dxj/>,m__<|dt}>,s_<|dx/>,m__
,s_ 
.m

SW-GBM:
All O.goo4so a.44s4ss — c.pgases —.4ste89 0.001158 jsgrare
(m-i)SW 000410 000375 001093
(rt}8W, unordered 000408 ‘Soioes 2002472
000465 001263
+ 000450 001210 1.85%
ts .000485 1263488149
s+ 000451 001216 1.884122
000430 001156 1.597029
000410 001093 2.689776
000408 2.
3001353 31001970 °8.247983
[000442 [001197 3.768638)
000476 2.902851
000483 001315 5.008366
000485 001313 2.187363,
000475 1203254514
2)stag’, 000484 001309 2.264886
Spatial gradient 000216 5.673956
Full gradient . 900244 000651_6.577181

0,000020-s seed, 5.2-km/s velocity
0.10-m seed 1.00-m seed
<ldtl>,s__
,m__<Jdt]>,s_</dx}>,m

SW-GBI
All

0.000421 0.655878 0.000373 1.254788
000408 924318 A
000408 -784399 000393
458182 000394
766041 000406
000448 1457550 000394
000439 .781610 .000383

000421 655882 .000372
000408 924340 000368

(mei)SW ooo
(-1)SW, unordered
mee

000408 -784390 000393

000452 1.235036 000466

000459 1.072379 -000479

000479 1.632003 -000523

000484 918189 000467

th -2)stagy . 000461 549581 000454
-2)stag’, unordered": . 000481 1 S588 000452
Spatial radient . 000234 000231

Full gradient . 1000257 2 197501 000237.

‘See notes at end of table 2



Library of Congress Cataloging in Publication Data:

Accuracy and precision of microseismic event locations in Rock Burst research
studies / by P.L, Swanson .. . [et al.]

Pp. cm. — (Report of investigations; 9395)

Includes bibliographical references (p. 31).

Supt. of Docs. no.: I 28.23:9395.

1. Rock bursts, 2. Microseisms, I. Swanson, P, L. (Peter L:) Il. Series: Report

of investigations (United States. Bureau of Mines); 9395.

TN23.U43 [ITN317] 622s—dc20 [622’.28]  91-23343 CIP



32

8. Boler, F. M., and P, L. Swanson. Computer-Automated
‘Measurement- and Control-Based Workstation for Microseismic and
Acoustic Emission Research. BuMines IC 9262, 1990, 9 pp.

9. Blake, W., F. Leighton, and W. I. Duvall. Microseismic Tech-
niques for Monitoring the Behavior of Rock Structures, BuMines B 665,
1974, 65 pp.

10. Ge, M., and H. R. Hardy, Jr. A Statistical Method of Evalua-
tion of AB/MS Source Location Accuracy and Transducer Array Geom-
etry. Paper in Proceedings of 30th U.S. Symposium on Rock Mechanics
(WV Univ., Morgantown, WV, June 19-22, 1989). Balkema, 1989,
Pp. 663-670.

11. Salamon, M. D. G., and G. A. Wiebols. Digital Location of
Seismic Events by an Underground Network of Seismometers Using the
Arrival Times of Compressional Waves. Rock Mech., v. 6, 1974,
pp. 141-166.

12. Godson, R.A. M. C. Bridges, and B, M. McKavanagh. A
32-Channel Rock Noise Location System. Paper in Second Conference
on Acoustic Emission/Microseismic Activity in Geologie Structures and
Materials, ed. by H. R. Hardy, Jr., and F. W. Leighton (Proc. Conf. PA
State Univ., University Park, PA, Nov. 13-15, 1978).  Clausthal-
Zellerfeld, FRP, Trans. Tech, Publ. D-3392, 1980, pp. 117-161.

13. Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T.
Vetterling. Numerical Recipes: The Art of Scientific Computing,
Cambridge Univ. Press, 1986, 818 pp.

14. Klein, R. W. Hypocenter Location Program HYPOINVERSE
Part I: Users Guide to Versions 1, 2, 3, and 4. U.S. Geol. Surv. OFR
78-694, 1978, 113 pp.

15. Lee, W. H. K,, and J. C. Lahr, HYPO71: A Computer Program
for Determining Hypocenter, Magnitude, and First Motion Pattern of
Local Earthquakes. U.S, Geol. Surv. OFR 75-311, 1972, 113 pp.

16. Lienert, B. R,, E. Berg, and L. N. Frazer. HYPOCENTER: An
Earthquake Location Method Using Centered, Scaled and Adaptively
Damped Least Squares. Bull. Seism, Soc. Am, v. 76, 1986, pp. 771-783.

17. Dahlquist, G., and A. Bjorck. Numerical Methods. Prentice-Hall,
1974, 573 pp.

18, Niewiadomski, J. Application of Singular Value Decomposition
Method for Location of Seismic Events in Mines. Pure and Appl. Geo-
phys., v. 129, 1989, pp. 551-570.

19. Riefenberg, J. A Simplex Method-Based Algorithm for Source
Location of Microscismic Events Associated With Underground Mining.
Paper in Proceedings of 30th U.S. Symposium on Rock Mechanics (WV
Univ., Morgantown, WV, June 19-22, 1989). Balkema, 1989, pp. 655-662.

20. Prugger, A. F,, and D. J. Gendzwill. Micro-Earthquake Location:
A Non-Linear Approach That Makes Use of a Simplex Stepping Proce-
dure. Bull. Seism, Soc. Am, v. 78, 1988, pp. 799-815.

21, Swanson, P. L. Stress Corrosion Cracking in Westerly Granite:
‘An Examination by the Double Torsion Technique. M.S. Thesis, Univ.
CO, Boulder, CO, 1980, 142 pp.

22, Corson, D. R, and W. R. Wayment. Load-Displacement Meas-
urement in a Backfilled Stope of a Deep Vein Mine. BuMines RI 7038,
1967, 51 pp.

23. Science Applications Intemational Corporation. Operations Man-
ual: Microseismic Processor Recorder Mode! DP-138, 1988, p. 6.

24, Lee, W. H. K,, and S. W. Stewart, Principles and Applications of
Microearthquake Networks. Academic, 1981, 293 pp.

25. Blair, D. P. Dynamic Modelling of In-Hole Mounts for Seismic
Detectors. Geophys. J. R. Astron. Soc, v. 69, 1982, pp. 803-818.

26, Dines, K. A., and R. J. Lytle. Computerized Geophysical Tomog-
raphy. Proc. IEEE, v. 67, 1979, pp. 1065-1073.

27, Mason, I. M. Algebraic Reconstruction of a Two-Dimensional
Velocity Inhomogeneity in the High Hazles Seam of Thoresby Colliery.
Geophys., v. 46, 1981, pp. 298-308.

28. Kormendi, A., T. Bodoky, L. Hermann, L. Dianski, and T. Kal-
man. Seismic Measurements for Safety in Mines. Geophys. Prospect.,
v. 34, 1986, pp. 1022-1037.

29. Young, R. P., D. A. Hutchins, M. Bostock, and P, Mottahed.
Geotomographic Imaging of Collapse/Solution Structures in Potash
Mines Using Seismic Methods. Paper in Proceedings of 28h US.
Symposium on Rock Mechanics (Univ. AZ, Tucson, AZ, June 29-July 1,
1987), Balkema, 1987, pp. 19-27.

30. Jech, J. Seismic Tomography in the Ostrava-Karvina Mining
Region, Pure and Appl. Geophys., v. 129, 1989, pp. 597-608.



2

880413.025912

wo

3100 3200 3300 3400 3500 3600 37.00 3800 3000 40.00

4100 ms

Figure 12.—Example of hardware picks missing first phase (P-wave) and picking later phase.

For the 40 events with acceptable-quality locations,
table 4 shows the average difference between the manual
and the software-screened hardware picks.'° (Hardware-
pick preprocessing eliminated two picks which differed

1An early prototype version of the first-arrival picking hardware pro-
vided four microsecond time resolution and was used in this comparison.
‘Travel-time pick resolution in the standard analog system (100 ys)
necessarily yields larger errors,

from the manual picks by as little as 10 us.) The
differences for the 216 arrivals which were analyzed ranged
from 0 to 1,170 ys with a distribution that approximated
an exponential function with a 142-ys time constant
(fig. 13). The average difference is 126 us, which repre-
sents a modest fraction of the average travel time to the
first hit accelerometer.


Report of Investigations 9395

Accuracy and Precision of Microseismic Event
Locations in Rock Burst Research Studies

By P. L. Swanson, L. H. Estey, F. M. Boler, and S. Billington

UNITED STATES DEPARTMENT OF THE INTERIOR
Manuel Lujan, Jr., Secretary

BUREAU OF MINES
T S Ary, Director


DISCUSSION

The relative contributions of measurement errors
in accelerometer positions, travel-time picks, and aver-
age velocity to errors in microseismic event locations
can be assessed for the analog system by comparing fig-
ures 10, 15, and 21, Table 5 summarizes the average
magnitudes of the discrepancy vectors of the 1,000 random
events in the control volume for gach source of error, A
composite figure incorporating these three sources of error
is shown in figure 22. ‘A. sanilar figure for the digital
system is shown in figure 23 and summarized in table 5.

‘he overall spatial discrepancies for the digital system are
approximately five times smaller than for the analog sys-
tem. When the lower average velocity value of 5.0 km/s

is assumed with a 4% error for the analog system, the
average composite discrepancy vector itude is
reduced to a little less than three times the value for the
digital system.

Table 5.—Average calculated location-diecrepancy
vector magnitudes, metere

Digital

Figure 22.—Composite plot for analog system showing loca-
tion discrepancies calculated using measured errors in accel-
‘erometer position, arrival-time picks, and assumed Isotropic
velocity. A, Plan view; B, vertical section looking northwest.

Figure 23.—Composite plot for digital system showing location
discrepancies calculated using measured errors in accelerometer
position, arrival-time picks, and assumed Isotropic velocity. A,
Pian view; B, vertical section looking northwest.


This subset problem leads to a vast number of potential
source location solutions. For m = 4 and equal weighting,
a three-dimensional location cannot be attempted since
only three equations can be formed and there are four un-
knowns. For m = 4 and unequal weighting, six independ-
ent equations can be formed from which 15 subsets can be
formed using four equations, six subsets using five
equations, and the total set of six equations; thus 22 dif-
ferent source locations would be possible. For m = 5 and
equal weighting, a unique three-dimensional source loca-
tion is obtained. For m = 5 and unequal weighting and
using all possible different subsets of 4, 5, 6, 7, 8, 9, and 10
equations, there are, in general, 848 different possible
source locations!

BLD BASIS SUBSETS

Each nonzero BLD basis function involves information
from three different receivers, say, j, k, and 1. Therefore,
the set of all possible BLD equations can be referenced by
specifying all possible triplets of j, k, and 1. Again, let
m=5:

yeni
ole

5

oo
=

5
45
364

3

yeN

5
46
3465

It is easy to see that there are at least m(m-1)(m-2)/6
different triplets. Since it is obvious from equation A-6
that interchanging the indices j and k results in only a sign
change of the equation, it may be thought at first that
there could be as many as m(m-1)(m-2)/2 different equa-
tions depending on the choice of Il. Fortunately, this turns
out not to be the case: the indices j, k, and | are all inter-
changeable; an even permutation of the indices jkl results
in the same equation A-6, an odd permutation results in
only a sign change of the equation A-6. So in the above
example of m = 5, the 10 different combinations for j, k,
and | give all the possible different equations A-6 to within
a sign change, regardless of which index is selected for j,
k, orl.

The subset of m-2 equations suggested by Blake (9)
could be indicated by the following for m = 5:

3°45

3

12
12
12 4
12

8,

37

where, for this subset of Blake’s, the ordering of the re-
ceiver indices indicates the order in which signals arrived
at the receivers (again, positive arrival-time ordering).
This type of subset of the BLD basis function, like the
Salamon and Wiebol subset of the SW-GBM basis func-
tion, over-emphasizes information for some receivers.
Here, the information of the first and second hit receivers
is overemphasized in the final solution.

A better subset of m-2 equations might be something
like

123 45
123
234
345

with or without arrival-time ordering. A balanced subset
of m equations might be something like

nd
ooo

4

45
1 45
with or without arrival-time ordering.

LINEAR DEPENDENCE OF DIRECT
BASIS FUNCTIONS

The question arises as to whether or not any of the
above subsets of the BLD basis function represents a basis
of independent equations from which the rest of the
m(m-1)(m-2)/6 equations can be formed, The answer is
no, even for the case of equal information weighting.

For equally weighted SW-GBM basis functions, two
equations of different jk indices can be selected to form a
third SW-GBM equation as long as one index is repeated.
For example, equation (A-1),, + (A-1), = (A-1),; with W,,
= Wy = Wy, noting that (A-1), = -(A-1),. Thus, for
equal weighting, either the Salamon and Wiebols m-1 sub-
set or the Godson m-1 subset gives a basis subset from
which the other equations can be formed.

For equally weighted BLD basis functions, three equa-
tions of different jkl indices can be selected to form a
fourth BLD equation as long as one index is repeated in
all three equations and the other indices are repeated in
two of the three equations. For example, equation (A-6),.,
+ AG)au + (A-6)iay = (A-6)aa With Was = Wisy = Wise =
Way, noting that (A-6),,, = -(A-6),. Thus, for the BLD
basis, there are at least (m-1)(m-2)/2 independent equa-
tions (more if weighting is unequal), with a maximum of


Figure 15—Comparison between 1,000 defined event lo-
cations and locations calculated with hardware-based first-arrival
time pick errors added to exact travel times, Lines connect
defined locations with calculated locations. A, Plan view; B,

vertical section view looking northwest.

velocity structures with extreme acoustic impedance mis-.
matches) and contrasted with the travel times that would
be expected without mine openings. Consider the sche-
matic in figure 164 showing a receiver attached to a mine
opening and seismic waves emanating from a source within
the rock on the opposite side of the opening. For suf-
ficiently long wavelengths (greater than 10 times the
mine-opening circumference), the presence of the mine

1%

: Lo.

20%

10%

3%
1% .

ee

Figure 16.—Influence of mine openings on selsmic wave travel
times. A, Approximate short-wavelength source (S) to receiver
(F) raypaths near mine openings; B, percentage difference in
travel times between virtual straight raypath and path around
mine opening as function of selsmic source location.

opening does not affect the wave propagation (25). In the
present study, however, the observed wavelengths of first
arrivals from microseismic events (1 to 5 m) are generally
on the order of the smallest mine-opening dimensions
(1 to 3 m), and hence their influence on first-arrival times
must be estimated.


from each other by less than 4 m in the horizontal plane.
The average component of discrepancy shown in the verti-
cal section is significantly larger and is vertically oriented.
The average magnitude of the three-dimensional discrep-
ancy vector between the locations is 7.1 m. This figure
provides one estimate of the minimum uncertainty in the
location of microseismic events using these techniques.

The previous comparisons between calculated locations
of blasts and microseismic events using the analog and
digital systems amply illustrate the need for improvement
in location accuracy. The remainder of this report dis-
cusses some of the factors which contribute to errors in
the calculated locations of calibration blasts and micro-
seismic events,

ERROR ANALYSIS

The following four steps are required for locating
seismic sources:

1, Determination of the spatial coordinates (x, y, 2.) of
all transducers in the receiving array.

2. Measurement of the arrival times, t,, at each station,
i, of the compressional wave (or other phase) generated by
the transient event.

3. Establishment of an appropriate model of seismic
wave velocity, v, structure in the area covered by the array.

4, Solving for the location (x, yp, Z) and origin time,
t,, using an appropriate technique—numerical or otherwise.

The nonlinear travel-time equation of the seismic signal
from the source to the ith receiver is

O = [(- %9)” + (= Yo)” + = m)*1M? - VEG; - ty). (1)

The velocity field is assumed to be both isotropic and
homogeneous, i.e., the signals (either compressional, P, or
shear, S, waves) propagate along straight raypaths in all
directions at a single value of velocity, v.

It is instructive to quantitatively examine the errors
associated with the microseismic source location problem.
From equation 1, one of the spatial coordinates of the
source can be expressed as

= % tL) O-¥0)?- @- 2) @)

By restricting the problem to the one-dimensional case, it
can readily be shown how uncertainties in the various
known quantities result in errors in the location solution.
The one-dimensional analog is

%H = % £ v(tj-t). @)

It is customary to obtain an error estimate, 6f, of the
function f = f(a, b, c, .. .), which in this case is x above,
by

Sea +

|*|-

at a |
Hoo | | 0 (2%)

+ | O(6b’) | +... (4)

where the terms on the righthand side of the equation are
just the magnitude of the terms of the Taylor series expan-
sion of f. As long as the uncertainties of a, b, c,... are
small in comparison to a, b, ¢, the higher order terms can
be neglected.

For the one-dimensional location problem,equation 3,
the error estimate is

J6x1 + 18v(t,- GI + Iv" 65 - tI 6)

Alternatively, the fractional uncertainty in location can
be expressed in terms of the source distance, d, from the
ith receiver,

| oe 5 -

[% Co) ; “at ®

where all the higher order terms are zero since the one-
dimensional equation is linear in the quantities x, t, t,, and
v. Equations 5 and 6 show how uncertainties in measure-
ments of the accelerometer position, 5x, seismic velocity,
6y, and travel time, 6(t; - t,), lead directly to uncertainties
in source location, 6%. This system of one-dimensional
linear travel-time equations leads to a coordinate solu-
tion with an average uncertainty that is weighted by the
measurement uncertainties for each parameter in each
equation.

For the two- and three-dimensional cases, the higher
order terms of equation 4 generally cannot be ignored,
particularly in the presence of significant uncertainties. As
measurement errors become large, the dependence of the
error in source coordinates (including t) on the mea-
surement uncertainties becomes increasingly nonlinear.

[6x] =

= bv

v



a pseudo-random number generator) to the perfect arrival
times in such a way as to simulate the effects of attenua-
tion and multipathing on the signals received at the re-
ceivers, €.g., arrivals at nearby stations were given smaller
random error than arrivals at distant stations. To help
simulate the effect of changes in receiver subarray ge-
ometry, on average 15% of all calculated arrivals were
randomly excluded.

Attempts were made to locate the 1,000 synthetic events
with each set of arrival times. Attempts included perfectly
and imperfectly known values of the velocity and receiver
coordinates. Solutions using all four linear basis func-
tions discussed above were examined with a varying num-
ber of m receivers. For the BLD basis, the full set of
m(m-1)(m-2)/6 equations and various subsets of equations
(those suggested by Blake and others subsets discussed
in appendix B) were tried, all with and without arriv-
al ordering. For the SW-GBM basis, the full set of
m(m-1)/2 equations and various subsets of equations were
tried, again, all with and without arrival ordering. For the
gradient bases, different starting solutions (¢.g., the results
of the direct solutions, location of first hit station,
geometric center of the array, etc.) were tried.

Several different numerical algorithms for solving the
linear systems of equations were tried. The code for the
OR algorithm given in Dahlquist and Bjorck (17) was re-
written by the second author in the programming language
C; the code for the SVD algorithm was adapted from the
Fortran code given in Press (13); a simplex code was sup-
plied by Riefenberg (19). Of these different codes of
numerical algorithms, the QR and SVD codes gave essen-
tially the same results (to floating-point accuracy) and both
gave more consistent results than the linear simplex code
used, since the simplex code was not robust. In principle,
however, a robust linear simplex code could give better
results than least squares algorithms such as QR or SVD
since the resulting solution of a simplex algorithm should
be less sensitive to outlying data points (19-20).

All of the resulting event locations were examined
graphically in three dimensions on an interactive graphics
workstation. To present the results in a simple way here,
two quantities are tabulated in tables 1 and 2 which apply
for both the QR- and SVD-based solution. The first is the
average of the magnitude of the difference between the
true origin time and the calculated origin time for the
1,000 hypothetical events, The second is the average of
the magnitude of the spatial discrepancy vector between
the true event location and the calculated location. Ta-
ble 1 gives thé average discrepancy for different amounts
of error added to original perfectly known arrival times,
receiver coordinates, and/or spatially constant velocity.
The results given in tables 1 and 2, and other results
summarized below, are based on close to 1,000 different

sue

trials of 1,000 synthetic event locations (close to 1 million
individual synthetic event locations).

On average, the direct solutions using a SW-GBM basis
give smaller discrepancies than the direct solutions using
a BLD basis for any type of error. For random errors in
arrival time picks or receiver coordinates, the lowest dis-
crepancies of the direct bases tended to be returned by the
full independent subsets, followed next by the complete
sets of equations. A full independent subset contains only
enough independent equations necessary to construct the
complete set of equations. Iteration with the gradient
bases reduces the discrepancies further; the same average
discrepancies for the gradient bases were obtained regard-
less of the starting solutions that were tried. Also, in-
creasing the arrival time errors by a factor of 4 increases
both the temporal and spatial discrepancies of the loca-
tions by a factor of close to 4, indicating that the arrival
time uncertainties are still small enough to be in the linear
realm of equation 4. The arrival time errors were in-
creased by yet another factor of 4 with still a near-linear
increase in the total average temporal and spatial dis-
crepancies of the locations. Similarly, increasing the
spatial uncertainties of the receivers by a factor of 10 in-
teases both the temporal and spatial discrepancies of the
locations by almost a factor of 10, indicating the larger of
the receiver uncertainties are just beginning to affect the
solutions in a nonlinear manner.

For random errors in arrival times or receiver coordi-
nates, the gradient bases gave smaller discrepancies than
the direct bases.

For systematic errors in the velocity, the results were
somewhat different. The smallest location discrepancies
were from particular subsets of the direct bases; these
subsets of m-1 SW-GBM and m-2 BLD functions used
arrival-time ordering and staggered combinations of in-
dices, The largest spatial discrepancies were returned by
the gradient bases, though these stilll yielded very low tem-
poral discrepancies. In fact, although not shown directly
in table 1, nearly all solutions based on the spatial- and
full-gradient bases tended to give the lowest travel-time re-
siduals, This is a surprising (significant?) result since a
minimum in the travel-time residuals is most commonly
accepted by seismologists as the solution providing the best
location. As indirectly suggested by the temporal dis-
crepancies in table 1, when there is a systematic error in
the velocity structure the best location in spatial coordi-
nates is not necessarily provided by the solution with mini-
mum travel-time residuals,

Weighting, if based on reasonable criteria, seems to re-
duce the location discrepancies, even when only a single
parameter is weighted. An explicit numerical example of
the event location process is given for each basis function
in appendix C.