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) 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.
5. 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 (STA/LTA) ratio of the energy of the vertical-component signals e(t) = z(t)². 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.
6. 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.
7. 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.

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.

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.

Back to top