In this electronic supplement, we describe in detail the kurtosis-based time-reversal (migration) method for earthquake location used in the main paper and include tables of parameters used to synthesize waveforms for tests and estimated locations and figures illustrating schematics of the method, synthetic waveforms, and kurtosis characteristics.
Time-reversal (TR) theory is a powerful tool that permits us to compute the sensitivity of a physical objective function or model with respect to its parameters. Propagating waveforms backward in time is an example of a TR calculation in geophysics. Backpropagating recorded information in time has long been used in seismic migration algorithms in geophysics (Baysal et al., 1983; Loewenthal and Mufti, 1983; McMechan et al., 1985; Tarantola, 1988) and in recent years has been used in adjoint imaging methods (e.g., Larmat et al., 2006; Larmat et al., 2009). In seismology, a TR experiment will often consist of transforming a network of seismic stations into an active array. In other words, a station that is generally a receiver and a recorder of the displacement field becomes an emitter of time-reversed waves, and the displacement field is propagated back as an excitation force (Larmat et al., 2006).
TR applications are mathematically equivalent to adjoint operations (e.g., Menke, 1989; Tarantola, 2005). Let us assume that data and model parameters are linearly related to each other by the expression di = Gij mj, in which d and m indicate the data and model parameters respectively, G is the kernel that is projected onto the model parameters in order to predict the data, and we use Einstein’s summation convention. Following Claerbout (2001; see Data and Resources), the adjoint operation of the above equation is given by m′j = Gji di. For a wavefield, we can re-express the adjoint operation as a TR operation (Kawakatsu and Montagner, 2008). For an observed waveform at a point, the result of the TR operation is given by a convolution operation:
Jn(t) = Gni(x′, t; x, t0) * di(t - t0),
in which J denotes the time-reversed seismogram or energy at time t, t0 denotes an arbitrary reference time, G denotes the Green’s function (GF) between the points x’ and x, and * denotes the convolution operator.
In a TR application for earthquake location, due to source-receiver reciprocity and the time-invariance of the elastic wave equation, the time-reversed wavefield migrates toward the source in reverse time along the same propagation path used by the forwardpropagating wave travelling from the source to the receiver. The waves interfere during backpropagation and focus at the source location at the correct origin time (Larmat et al., 2006; O’Brien et al., 2011).
The quality of a TR application depends strongly on the quality of the GFs used, and the latter depends both on our ability to simulate waveforms and our knowledge of the Earth’s structure. At long periods, the GFs are weakly dependent on 3D Earth structure, and TR event location works well at the global scale (e.g., Larmat et al., 2006). When we transpose the problem to a local context, as in this study, we have to work at much higher frequencies and are faced with the dual problem of simulating waveforms at such high frequencies and determining the local Earth structure sufficiently well so that the GFs will be consistent with the observed data.
We can greatly simplify the TR procedure for earthquake location by processing both observed and GFs in order to enhance that part of the signal that is most sensitive to changes in origin time and event location (i.e., the P-wave arrival time). The ideal processed seismogram should therefore contain strong peaks at the P-wave arrival times and be of low amplitude elsewhere. The GF for such data could then be approximated as similarly shaped peaks arriving at the times predicted by the Earth structure of the region of interest.
Our time-reversal (TR) method is a three-step procedure inspired by the work of Withers et al. (1999), Maggi and Michelini (2010), and Langet et al. (2014), and by the Waveloc algorithm. It is implemented using Python 2.7 with ObsPy libraries (Beyreuther et al., 2010).
The first step simplifies the observed waveforms by computing a characteristic function that produces high peaks at the first P-wave arrival times. We have followed Maggi and Michelini (2010) and Langet et al. (2014) in choosing the kurtosis characteristic function for this step. The kurtosis is the fourth statistical moment of the data and indicates how non-Gaussian a distribution is. A window within which the signal amplitudes are distributed according to a Gaussian distribution has a kurtosis of zero. In contrast, outliers due to non-Gaussian signals (such as seismic events, and especially P arrivals) produce high kurtosis values, because their amplitude distributions have sharp peaks and long, symmetric tails. Applying kurtosis within a sliding window on a raw seismogram will give a characteristic function that has low values while the sliding window only contains noise. It will peak sharply on the P arrival and then decay as the sliding window moves along the seismogram. Kurtosis can be an effective statistical tool for identifying P-wave arrivals and in recent years has been increasingly used to improve or replace the short-term average/long-term average (STA/LTA) characteristic function in automated arrival picking algorithms (e.g., Hadjileontiadis and Panas, 2002; Baillard et al., 2014).
Calculating an STA/LTA characteristic function is much faster than calculating an equivalent kurtosis characteristic function. In order to improve the speed of our TR algorithm, we first determine approximate P-wave arrival times by the STA-LTA technique, then window the data around these arrival times, remove mean and trend of the signal, merge all such windows together to form a single signal, apply a band-pass filter, and finally calculate the kurtosis characteristic function.
The second step of our TR method was inspired by Withers et al. (1999); it convolves these characteristic functions with Green’s functions calculated through a regional Earth model and then stacks them appropriately on a 3D grid of test hypocenters. The GFs should be calculated for a set of grid points in advance. If we were to use full-waveform GFs, this step would require a large storage volume and high computational costs. Because we are using kurtosis characteristic functions to simplify the data, we can approximate the GFs using a travel-time table of P-wave first arrivals for a given velocity model: we create a zero-value signal of the appropriate length and add an impulse of the appropriate width at the computed arrival time.
Performing an independent correlation for each grid point is time consuming and may include redundant calculations. In the absence of lateral variation in the velocity model, the GF for any station can be assumed identical for all grid points that are located at equal distances from that station. We use a 1D velocity model for the Rigan region, so we can take advantage of the reduction in the number of GFs to speed up the location algorithm. To do so, we modified the central migration routine of the Langet et al. (2014) Waveloc implementation to make it more similar to that of Withers et al. (1999).
Our migration routine is as follows. We store the 1D GFs row-wise in a Green’s functions matrix and store the processed data column-wise in a data matrix (see Fig. S1). Simple multiplications of the GF and data matrices, repeated on time-based sliding windows of the same length as the GFs, result in a matrix containing time-dependent correlation coefficients for each distance and each station. We then simply rearrange and stack these coefficients in geographical coordinates, summing the contributions of all stations for each test hypocenter.
The third step of our TR method was taken directly from Maggi and Michelini (2010) and Langet et al. (2014) and involves determining the local space–time maxima in the 4D stacked volume produced in the previous step in order to detect and approximately locate the seismic events. Grid points corresponding to local maxima in space and time are considered as the most probable hypocenters. We extract the maximum stacked correlation coefficient and its coordinates at each time step and save them as four independent time series (one for the maximum stack value and three for the geographical coordinates). We smooth the maximum stack value time series by a low-pass filter, then declare an event detection for each local maximum that exceeds a given threshold. We define left and right uncertainty bounds around this origin time by taking the times at which the stack descends down to 95% of its value at the local maximum. We then calculate the hypocentral coordinates and their uncertainties respectively as the mean and standard deviation of the coordinates within these left and right uncertainty bounds. We also record the absolute local maximum values and their corresponding coordinates to determine if they are more accurate than the averaged ones.
We tested our TR location method on four synthetic events recorded by six stations in the same configuration as the temporary stations deployed in the Rigan region. The synthetic seismograms are generated using the wavenumber-integration approach provided by Computer Programs in Seismology (see Data and Resources) and are based on real earthquake locations in the Rigan region. Each event has its own focal mechanism to test the effect of first-arrival clarification by the kurtosis characteristic function. The parameters of all events used in the synthetic test are listed in Table S1. All parameters except the origin times are based on real aftershocks of the first moderate earthquake in the Rigan region, as estimated by the Iranian Seismological Center (IRSC). The duration of each waveform is about 80 s, to ensure that all wavetrains are present, and its sampling rate is 100 samples/s. We use the CRUST2.0 velocity model (Bassin et al., 2000) for the Rigan region for the simulations, and we calculate the travel-time grids for migration through the same velocity model using routines from the NonLinLoc suite (see Data and Resources). In Figure S2a, the third and fourth events are barely visible because their magnitudes are much lower than those of the first two events.
We used two different band-pass filters for the synthetic tests: 1–10 and 1–4 Hz. We also varied the window length for the kurtosis between 1 and 3 s (this window length refers to the sliding window used to compute the kurtosis characteristic function and should not be confused with the length of the data segments kept around each preliminary P-wave arrival in the processing stage). We tested all combinations of filters and kurtosis windows to find the combination that results in the best impulse shape of the P arrival. As expected, using a 1–10 Hz band-pass filter for P phase identification in all tests is more robust than using a 1–4 Hz band-pass filter, because the impulsivity of the original arrival is more pronounced at higher frequencies. The shape of the kurtosis characteristic function depends on window length: wider windows produce larger kurtosis values; however, they also delay the arrival time by a few hundredths of a second. We have chosen to use a 3 s kurtosis window for our real-data application. Kurtosis characteristic functions for the synthetic events are shown in Figure S2b.
In order to enhance the quality of location, reduce the effect of unwanted kurtosis peaks such as S-wave arrivals, and eliminate artifacts due to the merging of data segments, we follow Langet et al. (2014) in adding some postprocessing to the kurtosis characteristic functions. We take the positive gradient of the kurtosis time series to reduce the arrival time delay induced by the length of the sliding window. We then convolve the processed characteristic function with a Gaussian source pulse of width 0.1 s at the positions where the derivative value exceeds a threshold (set to 500 in this study from the minimum value obtained from JZMN station). The resulting characteristic functions, which contain fewer unwanted peaks, are more symmetrical around the P-wave arrival time and produce smoother and more robust stacks.
We tested three low-pass smoothing filters for application to the maximum stack value time-series (with corner frequencies of 1, 2, or 3 Hz). The choice of the smoothing filters affects the location uncertainties. These uncertainties seem to be better for 1 Hz low-pass filter than for 2 and 3 Hz low-pass ones. We decided to use a 1 Hz low-pass filter for the real data.
We also tested several spatial sampling intervals from 2 to 5 km for potential source grids and the number of stations used for TR analysis. Unsurprisingly, the time and space uncertainties in the event locations are increased when larger sampling intervals are used. The choice of sampling interval is essentially a question of computational resources. For the real data application, we chose a 2 km sampling interval. The number of stations controls not only the location uncertainties but also the detectability of events (fewer stations implies fewer time series to stack and smaller overall stack values). The smaller events (EV3 and EV4 in Fig. S2 and Table S1) are no longer detected if we descend below six recording stations.
In summary, our synthetic tests enabled us to determine that best results should be obtained using a 1–10 Hz band-pass filter, a 3 s window length for kurtosis processing, a 1 Hz smoothing filter, and a 2 km spatial sampling interval for the potential hypocenter grid. Table S2 contains the locations of the synthetic events obtained using these parameters. Table S3 contains location information for events located using all six temporary stations of the Rigan network with an azimuthal gap less than 270° and a horizontal uncertainty of less than 20 km.
Table S1. Parameters used to synthesize waveforms for the synthetic tests.
Table S2. Locations retrieved by the synthetic location test using a kurtosis window of 3 s, a band-pass filter of 1–10 Hz, and a grid spacing of 2 km.
Table S3. NonLinLoc locations of the 46 well-constrained aftershocks of the 20 December 2010 Rigan earthquake.
Figure S1. A schematic illustrating the multiplication of Green’s functions and processed data for estimation of equidistant matrix coefficients (after Withers et al., 1999). (Asterisk denotes convolution.)
Figure S2. Synthetic data. (a) Synthetic waveforms for the four events listed in Table S1. (b) The kurtosis characteristic functions obtained after band-pass filtering the synthetic waveforms between 1 and 10 Hz and using a 3 s sliding window for the kurtosis calculation. Blue and green arrows indicate P- and S-wave phases highlighted by the kurtosis analysis. Red arrows show some artificial phase enhancement due to disturbances caused by merging event waveforms.
The dataset used for the analysis was collected by the Iranian Seismological Center (IRSC). The stations were equipped with broad-band Trillium 40 s sensors and Taurus digitizers configured with a sampling frequency of 100 Hz. The Waveloc code cited in this paper is open-source, is released under the CeCILL license and can be downloaded from http://amaggi.github.io/waveloc (last accessed December 2014). Waveloc is written in Python and was developed using the Enthought Python distribution under an academic license. The NonLinLoc package of Lomax et al. (2000) can be downloaded from http://alomax.free.fr/nlloc (last accessed December 2014). The synthetic seismograms used for the tests described in the electronic supplement were synthesized using the Computer Programs in Seismology package that can be downloaded from http://www.eas.slu.edu/eqc/eqccps.html (last accessed December 2014). The Claerbout software, Basic Earth Imaging v.2.4 is accessible http://sepwww.stanford.edu/sep/prof/ (last accessed June 2014).
Bassin, C., G. Laske, and G. Masters (2000). The current limits of resolution for surface wave tomography in North America , Eos Trans. AGU 81, no. 48, Fall Meet. Suppl., Abstract S12A–03.
Baillard, C., W. C. Crawford, V. Ballu, C. Hibert, and A. Mangeney (2014). An automatic kurtosis-based P- and S-phase picker designed for local seismic networks, Bull. Seismol. Soc. Am. 104, 394–409.
Baysal, E., D. Kosloff, and J. W. C. Sherwood (1983). Reverse time migration, Geophysics 48, no. 11, 1514–1524.
Beyreuther, M., R. Barsch, L. Krischer, T. Megies, Y. Behr, and J. Wassermann (2010). ObsPy: A Python toolbox for seismology, Seismol. Res. Lett. 81, no. 3, 530–533.
Hadjileontiadis, L. J., and S. M. Panas (2002). PAI-S/K: A robust automatic seismic P phase arrival identification scheme, IEEE Trans. Geosci. Rem. Sens. 40, no. 6, 1395–1404.
Kawakatsu, H., and J. P. Montagner (2008). Time-reversal seismic-source imaging and moment-tensor inversion, Geophys. J. Int. 175, 686–688.
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. Seismol. Soc. Am. 104, no. 1, 229–246.
Larmat, C., J. P. Montagner, M. Fink, Y. Capdeville, A. Tourin, and E. Clévédé (2006). Time-reversal imaging of seismic sources and application to the great Sumatra earthquake, Geophys. Res. Lett. 33, no. 19, L19312, doi: 10.1029/2006GL026336.
Larmat, C. S., R. A. Guyer, and P. A. Johnson (2009). Tremor source location using time reversal: Selecting the appropriate imaging field, Geophys. Res. Lett. 36, L22304, doi: 10.1029/2009GL040099.
Loewenthal, D., and I. R. Mufti (1983). Reverse time migration in spatial frequency domain, Geophysics 48, 627–635.
Maggi, A., and A. Michelini (2010). Waveloc: An algorithm for the detection and location of seismic sources within large, continuous waveform data volumes: The case of the L’Aquila earthquake sequence, Geophys. Res. Abstr. 12, Abstract 2010–12021.
McMechan, G., J. Luetgert, and W. Mooney (1985). Imaging of earthquake source in Long-Valley Caldera, California, 1983, Bull. Seismol. Soc. Am. 75, no. 4, 1005–1020.
Menke, W. (1989). Geophysical Data Analysis: Discrete Inverse Theory, Academic Press, New York, New York, 289 pp., ISBN: 0 12 490921-3.
O’Brien, G. S., I. Lokmer, L. De Barros, C. J. Bean, G. Saccorotti, J. P. Metaxian, and D. Patane (2011). Time reverse location of seismic long-period events recorded on Mt. Etna. Geophys. J. Int 184, 452–462.
Tarantola, A. (1988). Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation, Pure Appl. Geophys. 128, no. 1/2, 365–399.
Tarantola, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania.
Withers, M., R. Aster, and C. Young (1999). An automated local and regional seismic event detection and location system using waveform correlation, Bull. Seismol. Soc. Am. 89, no. 3, 657–669.
[ Back ]
