Page History
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 x0 acting at time t0 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=t0+T(x0, xi), 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 xi 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
- 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.
- Preprocessing. Detrend, demean, taper, band-pass filter within the [f1, f2] frequency band.
- 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)
- Alignment. Not necessarily all the traces start and end at the same time. Let's indicate with tsi and tei (i = 1, ... ,N)
- the
- starting
- and
- ending
- times,
- respectively,
- of
- the
- different
- seismograms.
- The
- application
- seeks
- the
- smallest
- starting
- time
- (tmin)
- and
- the
- largest
- ending
- time
- (tmax),
- for
- building
- a
- time
- vector
- between
- of
- time
- samples
- equallyspaced
- by
- 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.
- Calculation of the STA/LTA function Si(t). The characteristic function used for the back-projection is given by the short-term-average/long-term-average
- (STA/LTA)
- ratio
- of
- the
- energy
- of
- the
- vertical-component
- signals
- e(t) = z(t)
- 2
- .
- For
- each
- station,
- the
- STA/LTA
- ratio
- is
- calculated
- using
- a
- recursive
- algorithm
- described
- by
- the
- following
- equations (4):
- sta(i)
- =
- ks
- [c(i)
- ]
- +
- (1
- −
- ks)
- sta(i−1)
- ; lta(i)
- =
- kl
- [c(i)
- ]
- +
- (1−kl
- )
- lta(i−1)
- , where
- ns
- and
- nl
- are
- the
- number
- of
- samples
- of
- the
- short-
- and
- long-time
- windows,
- respectively,
- and
- the
- index
- i
- varies
- in
- the
- range
- between
- ns
- +
- nl
- and
- the
- last
- sample
- of
- the
- time
- series
- e(t)
- .
- Following
- Withers
- et
- al.
- (1999),
- the
- decaying
- constants
- Ks and Kl are set to 1/ns
- and 1/nl
- ,
- respectively.
- Calculation of predicted Travel Times. In the present implementation, the predicted travel times between the j-th trial source Xj and the i-th recorder located at xi are calculated assuming a homogeneous medium, and hence: Tij = || Xj - xi || / v (5) where v is the user-supplied wave velocity and || indicates the L2-norm.
- 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~_ forPij for P-wave
Calculate
the
theoretical
travel
time
_T
{^}S{^}{~}ij~_ forSij for S-wave
as
_T
{^}P{^}{~}ij~_ /1.73 Calculate the finite time sample corresponding to _Tij_ asPij /1.73
Calculate the finite time sample corresponding to Tij as tij = 1 + floor(
{_}TijTij / dt
{_})
Time-shift
of
characteristic
function:
CiCi (t +
TijTij)
end
Stack
CiCi,
i=1:N
to
obtain
B(x
{~}j
{~}, t)
end
For
each
time
sample
_k
_search
the
location
*_X{*}{^}max^_ at which _Xmax at which B(t
{~}k
{~},
{*}X
{*})
_takes
a
maximum,
and
store
those
time-
dependent
coordinates
of
maximum
brightness
as
_B'(t
{~}k
{~},{*}X{^}max{^}{*}){_}., Xmax).
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.Data Preparation
The User can obtain required SAC files via data conversion as well as by uploading them into his/her workspace from his/her computer.
Acquiring SEED file of interest:
- go to 'Data Search' (red field 1, Figure 1), select an episode of interest (red field 2, Figure 1) and select 'Seed Signal' or/and 'Accelerogram Seed Signal' from 'Data Type' (red field 3, Figure 1),
Figure 1. Event waveform selection.
- select an event of interest (here, event #7, Figure 1) by clicking on the 'ADD TO WORKSPACE' button at the right. A prompt window will appear, asking for the directory where to place file in (Figure 2),
Figure 2. Adding waveform to the workspace.
- access your workspace (for more information see My 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 (e.g. station MAR.MR, last signal).
Figure 3. Seed data file visualization.
Anchor | ||||
---|---|---|---|---|
|
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 Seed Signal view click on the 'ACTIONS' blue button at the top right, choose 'USE IN APPLICATION' tab and than 'SEED to SAC converter'. After successfully running the application the SAC data files will be listed within the 'SeedToSac' directory (Figure 4). See also SEED to SAC converter for more information on the application parameters.
Figure 4. Performing conversion into SAC with the 'SEED to SAC converter' application.
Running the application
Adding the application to workspace:
- from the 'Applications' page (for more information see Applications), choose the 'Waveform-based seismic event location', which will lead you to the description page of the application (Figure 5),
Figure 5. Description of the 'Waveform-based seismic event location' application.
- click on 'Add to Workplace' button at the upper right. A prompt window will appear, asking for the directory where to place the application in, then the User is redirected to his/her workspace and to the application input form (Figure 6). The application form could be later accessed also by choosing the 'WaveformBasedEventLocation'
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
106.
ApplicationThe application input form.
Input Parameters
Vp | P-wave velocity in km/s of homogeneous model. Travel times are calculated for an homogenous, isotropic Earth's model. For source located in the upper 10 km of crust, a value of Vp in between 4 and 6 is generally appropriate. Shear-wave velocities Vs are hard-coded as Vs=Vp/sqrt (3). |
Latmin | Minimum latitude of grid within which sources are sought (°, positive N) |
Lonmin | Minimum longitude of grid within which sources are sought (°, positive E) |
Latmax | Maximum latitude of grid within which sources are sought |
Lonmax | Maximum longitude of grid within which sources are sought |
zmin | min depth of grid (km below sea level, positive downward) |
zmax | max depth of grid (km below sea level, positive downward) |
nx | N. of grid nodes along EW direction |
ny | N. of grid nodes along NS direction |
nz | N. of grid nodes along Z direction |
Sources
are
sought
over
a
volume
whose
corners
at
the
Earth'sEarth’s surface
are
given
by
\[latmin,
lonmin
\]
and
\[latmax,
lonmax
\].
The
volume
extends
from
zmin
to
zmax
kilometers
beneathbeneath the
sea
level.
ThisThis 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.
TheThe 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
(Fig. 11(Figure 7).
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
._ !worddav37841a8f216c24695ab2b9f51b089d20.png|height=362,width=643! Figure 11 *ws* short window length for sta/lta (s) *wl* long window length for sta/lta (s) _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._
(a)
(b)
Wiki Markup |
---|
_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._
*f1* low corner frequency for filter (Hz) *f2* high corner frequency for filter (Hz)
_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._
*Otype* Output type:
*A{*}=save the entire B(x,y,z,t) grid
*R{*}=save only (x,y,z) maxima at each time step
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
Graphic files:- 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,Xmax)
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
.
Figure 7. Seismic catalog view.
ws | short window length for sta/lta (s) |
wl | long window length for sta/lta (s) |
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 8 illustrates the effect of different choices of ws and wl.
Figure 8. 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.
f1 | low corner frequency for filter (Hz) |
f2 | high corner frequency for filter (Hz) |
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.
Otype | Output type |
A = save the entire B(x,y,z,t) grid R = save only (x,y,z) maxima at each time step |
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
Graphic files:- 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,Xmax)