Improvements in seismic event locations in a deep western 

U.S. coal mine using tomographic velocity models and 

an evolutionary search algorithm

LURKA Adam1, 

1

Central Mining Institute, Katowice 40-166, Poland 

SWANSON Peter2 

2

NIOSH/Spokane Research Laboratory, Spokane WA, USA 

Abstract: Methods of improving seismic event locations were investigated as part of a research study aimed at reducing ground 

control safety hazards. Seismic event waveforms collected with a 23-station three-dimensional sensor array during longwall coal 

mining provide the data set used in the analyses. A spatially variable seismic velocity model is constructed using seismic event 

sources in a passive tomographic method. The resulting three-dimensional velocity model is used to relocate seismic event positions. An evolutionary optimization algorithm is implemented and used in both the velocity model development and in seeking 

improved event location solutions. Results obtained using the different velocity models are compared. The combination of the to

mographic velocity model development and evolutionary search algorithm provides improvement to the event locations. 

­

1 Introduction 

Developing strategies to deal with ground control 

hazards requires an accurate picture of the time-dependent, spatially-distributed, deformation processes 

attending rock mass failures. To the extent that these 

deformation processes generate detectable seismic 

emission, location of these events allows these processes to be delineated. Accurate event locations are 

therefore essential to developing an accurate characterization of critical deformation processes. 

The location accuracy of seismic events depends 

on three main factors: spatial distribution of seismic 

stations, accuracy of P-wave first-arrival picks and 

accuracy of the velocity model used in the location 

algorithm[1–5]. When the geometry of seismic stations 

is fixed the greatest influence on the location of seismic events is the accuracy of the first arrival picks 

and accuracy of the velocity model. 

P-waves are most often used in event locations because they are easy to identify on seismograms, the 

picking process is readily automated in software, and 

the pick errors are generally small. The precision of 

the first-arrival pick time is limited by the time resolution of the seismic recording equipment. Accuracy 

of the pick is sensitive to the rate at which the first 

arriving energy emerges from the background noise. 

The latter is a function of both source and propagation-path properties. In this regard, possibilities for 

improvement of pick times are rather limited. Therefore the last factor that can be improved is the velocity model used in the location procedure. This task is 

especially important when we locate seismic events 

in mines where the seismic network has a size on the 

order of a few kilometers and the required location 

accuracy is a few tens of meters. In seismically active 

mine areas one observes relatively large variability of 

P-wave velocity with time thus necessitating periodic 

updates of the velocity model used in the location 

process. This article presents an application of velocity images obtained from passive tomography calculations to the problem of seismic event location in a 

deep western U.S. coal mine. 

2 Method 

Passive travel-time tomography utilizes seismic 

events as sources in a simultaneous inversion for 

event locations and velocity structure. It has been 

used extensively in studies of earth structure and is 

being increasingly applied to mining problems[6–9]. 

The passive tomography problem can be solved in 

many different ways[6–8,10–11]. Every specific tomography algorithm depends on factors like: the type of 
input data, measure of data misfit, minimization 

method, ray tracing technique, etc. 

The present approach to solving passive tomography calculations has been performed using a misfit 

function based on the L1 norm. The L1 norm was selected over the least-squares norm L2 as it is much 

less sensitive to the effects of large errors in outlier 

data[10]. L1 is given by: 

m 

ˆ ˆ n

L v obs theor 

1 (, h ) = ¦¦ t ij − tij (1) i=1 j =1 

where vˆ is the vector of velocity model parameters; 

hˆ vector of seismic sources parameters (coordinates 

and origin times); n number of seismic events; m 

number of seismic stations; t

obs 

ij observed P-wave 

travel time; t

theor 

ij P-wave travel time described by 

model. 

The number of variables of the function (1) is a 

few hundred or more for typical passive tomography 

problems. To find the minimum of this function we 

use evolutionary algorithms that are well adapted for 

optimization problems with a large number of variables. Evolutionary algorithms use the principles of 

biological progression and try to mimic organic evolution that nature constructed during millions of 

year[12]. One can assert, in a simplification, that nature acts randomly to some extent but also aims at 

self perfection. Evolutionary algorithms mimic some 

of these features in optimization problems. 

Evolutionary algorithms have at least three advantages that are especially useful for passive tomography calculations. First, convergence approaches a 

global minimum. Second, it is easy to apply additional constraints to variables of the optimized function. For example, we can apply additional constraints to velocity values from borehole velocity logs 

that increase the stability of the algorithm. And third, 

evolutionary algorithms are adapted to operate on 

variables that have different physical units, i.e. they 

are not scaled. In the present application there are 

variables connected with seismic velocity, seismic 

event coordinates and event origin times. 

A short and very simplified scheme of evolutionary 

algorithms can be described as follows[12]: 

Step 0: Initialization 

There should be storage in a digital computer for 

two points of an n-dimensional Euclidean space 

Step 1: Variation 

Iteration g. Starting from point E(g)

 (position signed 

as a vector x (g)

E ) a second point N(g)

 (position signed 

as a vector x (g)

N,i ) is generated by using random perturbations. This point differs slightly from E(g)

. 

Step 2: Filtration 

The two points are associated with different values 

of the objective function F(x). Only one of them 

serves as a starting point in the next iteration g+1. 

Random perturbations are made by using an 

n-dimensional Gaussian probability density function 

with expected values equal to zero and standard deviations equal to σi (i=1, 2, ..., n): 

§ 2 n · 1 1 § z ·

wz ( 1 , ..., zn ) = exp  ̈ −  ̈ i  ̧  ̧

(2ʌ ) n n /2 ∏  ̈ ¦2 σ  ̧

σ © i =1 © i 1 i 1 i =1 

(2) 

Graphical illustration of the above steps is shown 

in Fig. 1. 

( ) g 

E x

( ) g 

Nx

( 2) g 

E x +

Fig. 1 Schematic illustration of the evolutionary 

algorithm[12] 

Passive tomography calculations, based on evolutionary algorithms and used in this paper, are less 

time consuming if one assumes that seismic rays are 

straight lines. This assumption would be a quite good 

approximation if the velocity contrasts were small 

and the area of study not large[8]. 

3 Data 

The seismic network consisted of 23 stations 

equipped with 4.5-Hz geophones. Nine stations were 

installed on the surface and 14 were installed underground (Table 1). Underground stations were largely 

located at seam level at a depth of 0.8 km. Both surface and underground seismic arrays operated as 

separate and independent networks and covered an 

area measuring 2.3 by 1.5 km. A GPS timing system 

was used to synchronize both networks, enabling 

them to operate as one 3-dimensional array. Additional details about the data acquisition and processing system can be found in Swanson[13]. 

Calculations were performed for seismic events 

that occurred in the vicinity of an active longwall 

measuring 250 meters wide by 1500 meters long. 

About 5000 seismic events were recorded and located 

using conventional techniques during the study interval. Out of this data set 64 seismograms from wellrecorded events with magnitudes M>0.7 were selected for the passive tomography analyses. P-wave 

first arrivals were manually picked on every seismogram and then the events were located using the layered velocity model given in Table 2. Calculated 
event coordinates, origin times and velocity values 

(Table 2) were used as starting parameters in the passive tomography calculations. 

Table 1 Surface and underground seismic stations

 coordinates

 Station’s name X (feet) Y (feet) Z (feet) 

Surface 

P1G1

P1G2

P1G3

P2G1

P2G2

P2G3

P3G1

P3G2

P3G3

 81320 

 82191 

 82941 

 82357 

 83326 

 83070 

 81621 

 82285 

 83154 

17799 

17459 

17031 

15768 

15940 

14613 

18792 

19220 

18563 

7698

8019

7822

8099

7595

7859

7685

7891

7991

Underground 

M39L

M33L

M26L

M21L

M13L

H21L

B10L

B18U 

B26L

B36L

B44L

B62L

B71L

TS1L

 86329 

 85691 

 84996 

 84518 

 83525 

 84030 

 83350 

82536 

 81466 

 80504 

 79796 

 79808 

 80647 

 83802 

15024 

14400 

13705 

13205 

12627 

16474 

13829 

14605 

15471 

16519 

17236 

18769 

19700 

12898 

5860

5937

6010

6062

6164

5613

5992 

5984 

5725

5588

5489

5242

5120

6125

Table 2 Layered velocity model used in starting location 

procedure 

Elevation (feet) Velocity (feet/second) 

7500~9000 8000 

6000~7500 11000 

4500~6000 14500 

Passive tomography calculations were performed 

in the area of the active longwall using both surface 

and underground seismic stations. Calculations were 

carried out in a 3-layer framework where every layer 

was subdivided into 8x8 cells. Additional subdivision 

was applied in the deepest layer where the seismicity 

was concentrated. The size of this finer scale velocity 

grid was 12x12 cells. The range of coordinates and 

the number of cells are shown in Tables 3 and 4. 

Table 3 Geometry of velocity grid 

X (feet) Y (feet) Z (feet) 

Begining of the grid 79000 12000 4500 

End of the grid 87000 20000 9000 

Numer of cells 8 8 3 

Table 4 Geometry of additional velocity grid in the 

deepest layer 

X (feet) Y (feet) Z (feet) 

Begining of the grid 81000 17000 4500 

End of the grid 83000 19000 6000 

Numer of cells 12 12 1 

4 Results 

Fig. 2a shows the seismic ray distribution in the 

area of interest and Fig. 2b shows the distribution in 

the region of finer velocity cell subdivision. Fig. 3 

illustrates the time residua between measured P-wave 

first arrivals and theoretical first arrivals for the layered velocity model shown in Table 2 and for the tomographic velocity model (Tables 3 and 4). Fig. 3 

allows an assessment of whether the travel time data 

fit the tomographic velocity model more accurately 

than the simple layered velocity model. Comparing 

both histograms of time residua one can see significant improvement between the theoretical and observed first arrivals for the passive tomography velocity model. 

Figs. 4a, 4b and 4c show calculated velocity images in the three layers. One can observe horizontal 

velocity variation in each velocity layer. In the deepest layer (Fig. 4a) P-wave velocity varies from 14000 

feet/sec to 14800 feet/sec. In the middle layer (Fig. 4b) 

velocity varies from 10000 feet/sec to 12000 feet/sec 

and in the shallowest layer velocity values vary from 

7300 feet/sec to 8700 feet/sec. A slightly different 

velocity image is shown in Fig. 4d for the dense velocity grid in the deepest layer. Horizontal velocity 

variations are rather minimal and lie between 14460 

feet/sec and 14510 feet/sec in the larger part of the 

region. 

The calculated velocity model presented in Fig. 4 

has been used to relocate seismic tremors with magnitude M0.7. Fig. 5 shows, in the planes XY and XZ, 

seismic tremors located by using both the velocity 

model from Table 3 and the velocity model determined in the tomographic inversion. The spatial differences in the event locations reach values exceeding 

200 feet (about 60 meters). Therefore, it is clearly 

seen that the tomographic velocity model fits the arrival-time data better and the relocation using the tomographic velocity model differs substantially from 

the locations using the layered constant velocity 

model. 

However, the relocation of seismic tremors did 

encounter some difficulties due to the existence of 

many minima in the solution procedure. In such cases 

the location process was repeated many times using 

different starting points and the position with the 

smallest value of the misfit function was assigned to 

be the event location. The occurrence of many minima in the misfit function was caused probably by 

applying tomographic velocity model that is strongly 

heterogeneous. 
(a) Region including all seismic stations (b) Area with finer velocity grid 

Fig. 2 Ray path coverage for seismic events with magnitude M0.7 

 

 

1XPEHURIWLPHUHVLGXD

LQVSHFLILHGLQWHUYDO

1XPEHURIWLPHUHVLGXD

LQVSHFLILHGLQWHUYDO

(a) 3-layer constant-velocity model (b) Tomographic velocity model 

Fig. 3 Time residua for the 3-layer constant-velocity model and for tomographic velocity model 

 

< IHHW

< IHHW



































































 



(a) The deepest layer: from Z=4500 to Z=6000 feet (b) The middle layer: from Z=6000 to Z=7500 feet 

 

 

< IHHW









 















 





  



 



























  









< IHHW

(c) The surface layer: from Z=7500 to Z=9000 feet (d) Area with increased of grid cells density 

Fig. 4 Map showing velocity model (feet/sec) obtained from passive tomography calculations 
 

 

 

< IHHW

= IHHW

(a) Horizontal plane XY (b) Vertical plane Z 

Fig. 5 Horizontal plane XY and vertical plane Z 

(Cross symbols show positions of seismic events calculated by using 3 layer constant velocity model; 

circle symbols show positions of seismic events calculated by using tomographic velocity model) 

5 Conclusions 

Methods to improve seismic tremor location accuracy in a deep western U.S. longwall coal mine by 

using tomographic velocity images have been considered. Tomographic velocity model and event location calculations have been carried out using arrival-time data from digital seismograms of mining-related seismic activity recorded by a local seismic network. Velocity images have been obtained 

using a passive tomography method. The calculated 

tomographic velocity model has been used to relocate 

seismic events. Passive tomography calculations have 

been carried out by using evolutionary algorithms 

that are well-suited to seeking the global minimum of 

the multidimensional functions. 

Differences between locations of seismic tremors 

obtained by using constant layered and tomographic 

velocity models were on the order of a few tens of 

meters. It follows that the tomographic P-wave velocity model has a significant impact on the seismic 

tremor locations. As the three-dimensional velocity 

model fits the theoretical P-wave travel times better 

than the constant velocity layered model one can say 

that this location procedure improves the event location accuracy. Additionally, in seismically active 

mines where one observes changes of velocity in time 

due to mining, periodic calculations of the velocity 

structure should be conducted using such passive tomographic methods. 

During relocations of seismic events using the tomographic velocity model one can observe unfavorable phenomenon associated with a relatively large 

number of local minima in the misfit function. It is 

likely caused by using complicated heterogeneous 

velocity model in the location algorithm. 

This paper was partially supported by the advanced 

NATO fellowship program. 

References 

[1] Rabinowitz N, Steinberg D M. Optimal configuration of 

seismographic network: a statistical approach. Bull 

Seism Soc Am, 1990. 

[2] Swanson P, Estey L, Boler F, Billington S. Mining-induced microseismic event location errors: accuracy 

and precision of analog and digital systems. Pure and 

Appl Geophys, 1992(139): 375–404. 

[3] Mendecki A J. Seismic Monitoring in Mines. Cambridge: 

Chapman & Hall, 1997. 

[4] DubiĔski J, DrzĊĨla B. Location of seismic events in 

mines. In: Biblioteka Szkoly Eksploatacji Podziemnej. 

Kraków, 1995. (In Polish) 

[5] Gibowicz J G, Kijko A. An Introduction to Mining 

Seismology. London: Academic Press, 1994. 

[6] Iyer H M, Hirahara K. Seismic tomography. Theory and 

Practice. London: Chapman & Hall, 1993. 

[7] Maxwell S C, Young R P. A comparison between controlled source and passive source seismic velocity images. Bull Seism Soc Am, 1993, 83(6): 1813–1834. 

[8] Maxwell S C, Young R P. Seismic velocity inversion 

from microseismic data. Seismic Monitoring in Mines. 

London: Chapman & Hall, 1997. 

[9] Lurka A. Seismic hazard assessment in the Bielszowice 

coal mine using passive tomography. Seismogenic Process Monitoring. Tucson: A. A. Balkema Publisher, 2002. 

[10] Tarantola A. Inverse Problem Theory. Amsterdam: Elsevier, 1986. 

[11] Shalev E, Lees M. Cubic B-splines tomography at loma 

prieta. Bull Seism Soc Am, 1998, 88(1): 256–269. 

[12] Schwefel H P. Evolution and Optimum Seeking. London: 

John Wiley & Sons, 1995. 

[13] Swanson P. Development of an automated pc-networkbased seismic monitoring system in dynamic rock mass 

response to mining. In: Proceedings of the 5th International Symposium on Rockbursts and Seismicity in Mines. 

Johannesburg, 2001: 11–17.