Abstract

This application implements a location procedure which allows imaging the distribution of seismic sources in both time and space. Using the predicted arrival times for a set of trial locations and origin times, the code calculates a brightness function by stacking a characteristic function observed at all the stations. The spatial and temporal distribution of sources is then identified by a systematic search throughout the model space and time for the maximum brightness. The main advantage of this procedure is that it exploits waveform information (both arrival times and relative amplitudes) without the need of pre-assembled phase-picking data. The method, therefore, is particularly well-suited for emergent signals and/or for real-time applications. The technique has been largely inspired by the previous works from Kao and Shan (2004), Grigoli et al. (2013, 2014), Langet et al., (2014), to which the reader is referred for a more comprehensive treatment of the methodology.

Method

Consider a seismic source positioned at x₀ acting at time t₀ and radiating a wavefield which is recorded by N receivers located at x. This signal is observed at the i-th recorder at the time t=t₀+T(x₀, x_i), where T is the source-to-receiver travel time.
Now consider a characteristic function CF(t) which results from a transformation of the original signal u(t) aimed at enhancing the energy of the incoming signal, in turn alleviating the waveform distortions due to path and source effects. A brightness function (Kao and Shan, 2004) is then defined as the stack of the time-shifted CFs:

(1)

which also corresponds to the standard beam-forming estimate for the case of arbitrarily-shaped wavefronts. From eq. (1), it is obvious that the B(t) function will attain a maximum when the time shifts x_i applied to the CFs from individual stations correspond to the actual travel times from the source. The procedure begins by parameterizing the volume which is supposed to contain the source(s) with a regular grid of M nodes located at X; travel times T(X, x) are calculated from each of these nodes to all the stations, and stored for subsequent use. Then, an exhaustive grid-search is conducted throughout all the nodes of the mesh: for each tentative source location, the CFs from all the stations are migrated (time-shifted) by the predicted travel times for that node, and then stacked to create a brightness function which is function of time and space (Fig. 1b):

(2)

Whenever the time t and grid node position X correspond respectively to the actual origin time and source location, the CFs will sum in phase thus yelding a maximum for B(t,X).
The procedure is eventually repeated using S-wave travel times, so that a final brightness function is obtained from multiplication of the two individual brightnesses (Grigoli et al., 2013; 2014):

(3)

Figure 1. Schematic diagram of the application. (a) Data processing: calculation of an appropriate characteristic function to maximize sensitivity to the presence of seismic events. (b) Migration: the processed waveforms are shifted back in time according to the predicted travel times from a potential hypocentral location, and then stacked. This procedure is repeated for a 3D grid of potential locations, and for both P- and S-waves travel times. (c) At each time step we store the maximum of the point stacks and the location it corresponds to, thereby allowing simple detection and event location.

Implementation

details

1. Read seismic data from the N stations. Data are read in SAC format (one file for each station). Station coordinates (lat, lon, elevation) are taken from SAC file's header.
2. Preprocessing. Detrend, demean, taper, band-pass filter within the [f1, f2] frequency band.
3. Decimation. Traces are down-sampled at 0.5 / f2 sampling interval (for instance, if f2 = 10Hz, the new sampling interval dt will be 0.05s, so that the Nyquist frequency corresponds to f2)
4. Alignment. Not necessarily all the traces start and end at the same time. Let's indicate with ts_i and te_i (i = 1, ... ,N)

1. the

1. starting

1. and

1. ending

1. times,

1. respectively,

1. of

1. the

1. different

1. seismograms.

1. The

1. application

1. seeks

1. the

1. smallest

1. starting

1. time

1. (tmin)

1. and

1. the

1. largest

1. ending

1. time

1. (tmax),

1. for

1. building

1. a

1. time

1. vector

1. between

1. of

1. time

1. samples

1. equallyspaced

1. by

1. dt spanning the tmin – tmax interval. Seismograms are then interpolated according to this time vector, so that all the traces start and end at the same time, and are exactly synchronized.
2. Calculation of the STA/LTA function S_i(t). The characteristic function used for the back-projection is given by the short-term-average/long-term-average

1. (STA/LTA)

1. ratio

1. of

1. the

1. energy

1. of

1. the

1. vertical-component

1. signals

1. e(t) = z(t)

1. ²

1. .

1. For

1. each

1. station,

1. the

1. STA/LTA

1. ratio

1. is

1. calculated

1. using

1. a

1. recursive

1. algorithm

1. described

1. by

1. the

1. following

1. equations (4):

1. • sta(i)

1. • =

1. • ks

1. • [c(i)

1. • ]

1. • +

1. • (1

1. • −

1. • ks)

1. • sta(i−1)

1. • ; lta(i)

1. • =

1. • kl

1. • [c(i)

1. • ]

1. • +

1. • (1−kl

1. • )

1. • lta(i−1)

1. • , where

1. • ns

1. • and

1. • nl

1. • are

1. • the

1. • number

1. • of

1. • samples

1. • of

1. • the

1. • short-

1. • and

1. • long-time

1. • windows,

1. • respectively,

1. • and

1. • the

1. • index

1. • i

1. • varies

1. • in

1. • the

1. • range

1. • between

1. • ns

1. • +

1. • nl

1. • and

1. • the

1. • last

1. • sample

1. • of

1. • the

1. • time

1. • series

1. • e(t)

1. • .

1. • Following

1. • Withers

1. • et

1. • al.

1. • (1999),

1. • the

1. • decaying

1. • constants

1. • Ks and Kl are set to 1/ns

1. •  and 1/nl

1. • ,

1. • respectively.

2. Calculation of predicted Travel Times. In the present implementation, the predicted travel times between the j-th trial source X_j and the i-th recorder located at x_i are calculated assuming a homogeneous medium, and hence: T_ij = || X_j - x_i || / v (5) where v is the user-supplied wave velocity and || indicates the L2-norm.
3. Calculation of the Brightness function B(x{~}j{~},t):

For

each

trial

source

j

For

each

receiver

i

Calculate

the

theoretical

travel

time

T

^P_ij for P-wave

Calculate

the

theoretical

travel

time

T

^S_ij for S-wave

as

T

^P_ij /1.73

Calculate

the

finite

time

sample

corresponding

to

Tij

as

tij = 1 + floor(

T_ij / dt

)

Time-shift

of

characteristic

function:

C_i (t +

T_ij)

end

Stack

C_i,

i=1:N

to

obtain

B(x

_j

, t)

end

For

each

time

sample

k

search

the

location

X^max at which B(t

_k

,

X

)

takes

a

maximum,

and

store

those

time-

dependent

coordinates

of

maximum

brightness

as

B'(t

_k

,

X

^max

)

.

end

Step by Step

1. Data preparation

• Access the website https://tcs.ah-epos.eu/, and log in or sign up for a new account (see the Quick Start Guide for IS-EPOS for more information on the log in and sign up procedure).
• Proceed to the episodes page (see the Quick Start Guide for IS-EPOS for more information on episodes). In this case, we are using the USCB: regional seismicity and ground motion associating underground coal mining
• From the menu at the left side, proceed to "Event Related Waveforms" (Figure 3).

Figure 3: Selection of 'Event related Waveforms'.

• Select an event of interest (here, event #9; Fig. 4); by clicking on the central button at the right, add a waveform related to it (seed data file) to your workspace. A prompt window will appear, asking for the directory where to place file in (Fig. 5).

Figure 4. Event waveform selection

Figure 5. Adding waveform to workspace

• Access your workspace (see the Quick Start Guide for IS-EPOS for more information on workspace) and proceed to the directory where you stored your seed data file. Click on it for visualization / checking. It is important that you take note of any bad / noisy channel, as this/these must be ignored in the subsequent location procedure (for instance, in the example of Figure 6, station MAR.MR , seventh from the top).

Figure 6. Seed data file visualization

2. Data conversion

Seismograms are contained in a seed volume which, for being used in the Waveform-based seismic event location, must be converted to a series of single-channel SAC data files.

WARNING:the format conversion also implies reading of station coordinates from the seed volume (extension .seed), and their writing into the SAC data files' header. If the original file is a miniseed volume (extension .mseed) rather than a full seed one, it will not contain stations metadata, and therefore the resulting SAC data files will miss that information and the location procedure will fail.

From the ACTIONS blue button at the top right, choose Use in Application → Seed Converter. At the program's prompt (Fig. 7), run the converter, selecting the desired ground motion component (usually, the Z one) and SAC format.

Figure 7. Performing conversion into SAC with Seed Converter application

At the end of the process, the SAC data files will be listed within the SeedConverter directory.

See also Seed converter User Guide for more information on the application parameters.

3. Running the application

Adding the application to workspace

• From the Applications page (see the Quick Start Guide for IS-EPOS for more information on applications), choose the Waveform-based seismic event location, which will lead you to the description page of the application (Figure 9)

Figure 8. Applications list

Figure 9. Description of the Waveform-based seismic event location application

• Click on Add to Workplace button at the upper right, and then proceed to your workspace by clicking on the above button, which leads you to the application input form (Figure 10). The application form could be later accessed also by choosing the WaveformBasedEventLocation from the tree to the left.

Figure 10. Application input form

Input Parameters

Sources

are

sought

over

a

volume

whose

corners

at

the

Earth’s surface

are

given

by

[latmin,

 lonmin

]

and

[latmax,

lonmax

].

The

volume

extends

from

zmin

to

zmax

kilometers

beneath the

sea

level.

 This volume

is

then

discretised

with

a

mesh

of

regularly-spaced

grid

nodes.

[nx,

ny,

nz

]

indicate

the

number

of

grid

nodes

to

use

along

the

longitude,

latitude

and

depth

directions,

respectively.

 The higher

the

number

of

nodes,

the

better

the

possible

source(s)

will

be

resolved.

At

the

cost,

however,

of

longer

computational

times.

It

is

strongly

recommended

not

using

number

of

grid

nodes

larger

than

50-70

along

the

three

directions.

To

have

an

idea

of

the

final

resolution,

just

consider

an

example

in

which

your

grid

extends

along

the

NS

direction

by

one

degree,

and

that

you

choose

ny=50.

Since

the

length

of

a

latitude

degree

is

on

the

order

of

110

km,

your

grid

nodes

along

the

NS

(y)

direction

will

be

spaced

by

about

110/5

=

2.2

km.

The

same

holds

for

the

EW

(longitude)

extension

of

the

grid,

with

the

difference

that

the

length

of

a

longitude

degree

depends

on

the

latitude.

The

grid

must

encompass

the

potential

source

volume,

it

is

necessary

an

a

priori

knowledge

about

the

possible

location

of

the

signals

under

examination.

If

a

catalog

is

available

for

the

selected

episode,

then

you

can

use

the

extreme

values

for

latitude

and

longitude

as

provided

in

the

informative

page

of

the

catalog

(

Figure 11).

If

a

catalog

is

not

available,

a

good

choice

is

to

use

a

volume

whose

horizontal

projection

includes

all

the

stations

which

are

going

to

be

used

in

the

calculation.

Image Added

Figure 11. Seismic catalog view

These are the length of the time windows for calculating the short- and long-term averages of the seismograms. There are no obvious choices for setting these values. In general, one wants the STA/LTA ratio to highlight the main phases (P, S), in turn suppressing both the noise preceding the signal of interest and the coda waves following both P and S-wave arrivals. A reasonable choice is to set both wl and ws shorter than the S-P differential time, with a wl / ws ratio on the order of 5. Figure 12 illustrates the effect of different choices of ws and wl.

Image Removed

Image Added

(a)	Image Modified (b)

Figure

12

. Effects

of

different

window

lengths

on

the

calculation

of

the

STA/LTA

function.

While

being

noisier,

shorter

windows

preserve

details

of

the

main

body-wave

arrivals.

Minimum and maximum frequency of the band-pass filter. Be careful that f2 must be lower than the Nyquist frequency Fn = 0.5/dt,

where

dt

is

the

sampling

interval

of

the

seismograms.

F1

must

be

higher

than

the

fundamental

frequency:

f1

>

1/T,

where

T

is

the

toal

time

length

of

the

seismogram.

A

option:

the

output

is

saved

as

an

ASCII

file

with

5

columns:

time,

x-location

of

grid

node,

ylocation

of

grid

node,

z-location

of

grid

node,

Brightness

as

a

function

of

[t, x, y, z

].

It

may

be

quite

a

large

file.

R

option:

the

output

is

saved

as

an

ASCII

file

with

4

columns:

time,

x,y,z

location

of

the

grid

node

taking

the

largest

brightness

at

that

particular

time.

Results

• An intermediate plot of the original seismograms and characteristic functions at all the stations:

• (x,y), (x,z) and (y,z) brightness sections in correspondence of the time of maximum brightness:

• A final plot of the time-dependent locations of maximum brightness:

Numeric:

• The complete data file containing the space-time dependent brightness B(t,X) or, alternatively, the time-dependent coordinates of maximum brightness B'(t,X^max)

References

Honn, K. and S.J. Shan, 2004. The Source-Scanning Algorithm: mapping the distribution of seismic sources in time and space. Geophys. J. Int. 157, 589–594 doi: 10.1111/j.1365246X.2004.02276.x

Grigoli, F., Cesca, S., Vassallo, M. & Dahm, T., 2013. Automated seismic event location by travel-time stacking: an application to mining induced seismicity, Seism. Res. Lett., 84, 666–677.

Grigoli, F., Cesca, S., Amoroso, O., Emolo, A., Zollo, A., & Dahm, T. (2013). Automated seismic event location by waveform coherence analysis. Geophys. J. Int., 196, 1742-1753.

Langet,N., A. Maggi, A. Michelini, and F. Brenguier, 2014. Continuous Kurtosis-Based Migration for Seismic Event Detection and Location, with Application to Piton de la Fournaise Volcano, La Réunion. Bull. Seism. Soc. Am., 104, 229–246. doi: 10.1785/0120130107.

Related Documents