Arrival Time Measurement Techniques in MsPASS#

Gary L. Pavlis#

Overview#

Arrival time measurements are the oldest and most fundamental data of seismology. They are the fundamental data used to produce all current generation earthquake catalogs and were used in one way or the other in everything we know about earth structure from seismic waves. There are two fundamentally different problems a research computing framework like MsPASS needs to address to support the near universal dependence of seismology research on arrival time data:

  1. Data management tools to efficiently handle all types of timing data.

  2. Tools to measure timing relationships from waveform data. Note there are two classes of timing measurements that we need to handle: (a) absolute time tags like traditional phase picks, and (b) relative timing relationships between parts of the same waveform or between two or more different waveforms.

For item 1, the fact that timing data of one kind or another are a universal concept of seismology means that every scientist has a set of both specialized and generic tools they use for handling their data. For that reason the role MsPASS can play for most people is as a tool for low-level data handling used in combination with one or more other tools they use as a component of their research. For that reason much of the issue with item 1 is an import/export problem. One of the most important features of the MongoDB database in MsPASS is the complete flexibility it provides in what it stores. It provides a generic way to store any format of timing data we have yet encountered. Every source we know of stores such data either as a set of text files or as tables in a relational database. Import/export of such data is discussed in the Importing Tabular Data section of the MsPASS documentation.

This section is focused mainly on item 2. That is, it describes tools in MsPASS for creating arrival time measurements from waveform data. We emphasize, on the other hand, that there are multiple commercial and open-source systems for processing and data management of traditional earthquake monitoring networks. We view the problem of earthquake catalog generation as solved and best handled by one of the systems designed specifically to do that job well. MsPASS is a system for more specialized research problems. The tools currently available are:

  1. A multichannel correlation algorithm building on earlier work by Pavlis and Vernon (2011). That algorithm can be used for estimating travel time residuals for teleseismic phases. It also computes a robust stack of the inputs after they are time aligned by cross-correlation. It differs from the original application as it is designed for automated processing while the original application was used with an integrated graphical user interface.

  2. ObsPy has several trigger and picking methods described here. Below we discuss how they can be used with data managed by MsPASS.

  3. MsPASS provides a wrapper around SeisBench’s PhaseNet model for annotating P-wave arrival probabilities. It is described in MsPASS PhaseNet Picking below.

The rest of this section has subsections on each of the items above.

Teleseismic Arrival Time Measurement#

Fundamentals#

There are some fundamental properties of teleseismic body-wave data that any algorithm for working with such signals must account for:

  1. Teleseismic body-wave signals vary in space at spatial scales far larger than signals from local and regional events. The reasons are a fundamental property of the earth and a topic outside the scope of this manual. The key fact, however, is that coherent processing of signals from local and regional events can only be done over the scale of one or two wavelengths. In contrast, many studies demonstrate coherent processing of teleseismic P and S waves are possible over distances of thousands of kilometers, which correspond to hundreds of wavelengths.

  2. The combination of source spectrum dependence on magnitude and the strong frequency dependence of the microseisms create variations in data spectra that can create large errors in measurements if not handled properly. One perspective is that a large fraction of traditional signal processing algorithms inherited from oil and gas processing have explicit or implicit assumptions that noise is “white”. Modern broadband data noise is anything but white. It is very “colored” by microseisms. Similarly, traditional algorithms from oil and gas processing have an implicit assumption the source is constant for all data. That is rarely true with earthquake data (the only exception is repeating earthquakes). As a result a fundamental thing teleseismic processing must handle is that the optimal bandwidth for signal processing is always dependent on earthquake magnitude and the noise spectrum of the seismic station being analyzed.

A corollary of item 2 is that any automated processing algorithm needs to be “robust”. In this case that means it will automatically handle a mix of data of variable quality and extreme outliers. In my experience automated arrival time estimation of any type can be treated as two different problems best handled by different processing algorithms:

  1. Automated discarding of the lost causes. Data can be a lost cause for a long list of reasons from the completely dead channel or always noisy channel to a one-off problem created by some local noise source that overwhelms the signal you want to see during a time window of interest.

  2. Robust handling of data with signals of variable quality. For the largest earthquakes this is a minor concern. The issue is fundamental, however, for the smallest events that may have usable signals only on a fraction of the receivers being processed. Unfortunately, due to the magnitude frequency relationship there are always far more of the marginal signals to handle than the larger, easier ones.

The approach in MsPASS is shaped by my experience in developing and using the dbxcor program described in Pavlis and Vernon (2011). In dbxcor we addressed both of these issues with the graphical user interface. The key idea was to iteratively improve a stack of the array data by repeating the steps of (a) compute a robust stack, (b) sort by some quality metric, (c) kill the junk, and repeat until satisfied. MsPASS provides tools for accomplishing that same process in an automated workflow. Examples showing how the tools can be linked together to accomplish that are found in mspass tutorials. I would emphasize, however, the tools there are not the last word on this problem and creative solutions are needed to improve performance and efficiency. In particular, using machine learning to auto-edit data is an obvious and likely valuable way to accomplish auto-editing.

MsPASS Multichannel Correlation Method#

Data Preparation#

An assumption of the multichannel correlation algorithm in MsPASS is that the data have been preprocessed through some variation of the following steps:

  1. The working dataset is a version of a common-source gather. By that I mean the data set is indexed in a way that all the waveforms linked to a particular seismic source can be assembled into a set of ensembles. For parallel processing that means the data set is an RDD/bag of ensembles. The members of the ensembles are assumed to span a time range around the seismic phase that is to be analyzed.

  2. All waveforms are normalized so source and receiver coordinates can be loaded from the database during the initial read operation or within the workflow when needed.

  3. Although not required, I have found that in practice all waveforms in the ensemble should normally have some basic low-level processing. For high quality data like USArray data that can be as simple as demeaning and scaling the data by a constant to compensate for gain variations. For more heterogeneous data a more sophisticated response correction may be necessary to assure the data are all normalized to a common response.

  4. The data are required to be resampled to a common sample rate to run the main processor called align_and_stack. The generic MsPASS function to accomplish this task is called resample.

  5. Time tags need to be defined in the Metadata container of each waveform to define at least an initial estimate of the arrival time of the phase of interest. These can be previously measured arrivals that are to be refined or model-based estimates computed from source coordinates, receiver coordinates, and an earth model. The MsPASS tutorial notebooks contain many examples of how to do this using obspy’s tau-p travel time calculator.

  6. The data should be shifted from UTC times to what in the docstrings we call the “arrival time reference” frame. That is, we expect the data have been shifted with the function ator with the shift time as the initial arrival time estimate set previously. That means a plot using the time_axis method for each atomic datum will have 0 as the initial arrival time estimate.

  7. The working ensembles are TimeSeriesEnsemble objects containing a component of data appropriate for the phase being analyzed. For P data that can be simple vertical components while for S it always demands the data were processed to Seismogram objects and oriented to radial or transverse components. For most modern data I would recommend all data be assembled to begin with as SeismogramEnsemble objects, rotated to LQT or with the free_surface_transformation operator, and then the appropriate component extracted using the ExtractComponent function.

  8. The data should be passed through the algorithm called broadband_snr_QC. That function computes Metadata attributes that are required for two purposes described below: (a) a single value that can be used as the “best” signal used as a seed for the multichannel algorithm, and (b) a pair of attributes that can be used to filter the data to an optimal frequency band for each ensemble. Plug in replacements are possible but would require careful looks at the python functions that utilize those attributes in the MCXcorStacking module.

Algorithm Background#

The top-level function for processing teleseismic body waves is align_and_stack. As the name suggests, it does two distinctly different tasks:

  1. Aligns a set of input waveforms by cross-correlation so all the waveforms match as closely as possible when plotted on a common, relative time base. (see e.g. Figure 1 of Pavlis and Vernon (2011))

  2. Produce a robust stack of the time-aligned data. After the first alignment stage the robust stack is used as the correlation reference. The algorithm name contains “MCXcor” which is shorthand for “Multichannel X(cross-)correlation” to contrast it with pairwise cross-correlation algorithms following the older work of VandeCarr and Crosson (1990).

Any user who needs to use the align_and_stack function of MsPASS should plan to have a copy of the Pavlis and Vernon (2011) paper for reference. The primary theory behind align_and_stack is documented in that paper. There are, however, some large differences between the MsPASS implementation and the publication. Pavlis and Vernon’s paper(2011) describes an analyst tool for measuring teleseismic body wave phase arrival times using multichannel cross-correlation they called dbxcor. dbxcor has an integrated graphical user interface that allows the user to set some key processing parameters interactively and graphically edit the data to discard bad and unacceptably noisy data. align_and_stack, in contrast, was designed as a purely automated tool that could be applied to large data sets and produce quality results without any human intervention. The next section documents the current algorithms used to automate what was done interactively in dbxcor. Note this topic is a research area that could be improved with alternative algorithms. For example, it is an obvious candidate for machine learning.

MsPASS Automation#

Adapting the algorithm of dbxcor to work as an automated tool required developing algorithms to define three key parameters dbxcor required the user to set via the graphical user interface:

  1. The algorithm used for cross-correlation ultimately uses the array stack to correlate with each signal in an input TimeSeriesEnsemble. However, to start the iterative sequence the input ensemble has to be time aligned to first order or the stack will be grossly distorted. dbxcor solved that issue by requiring the user to select the signal in the ensemble that would be used for initial correlations. The MsPASS implementation uses a python function in the module mspasspy.algorithms.MCXcorStacking called extract_initial_beam_estimate. That function is designed to use the output of the MsPASS function broadband_snr_QC, which creates a “subdocument” (aka python dictionary) with a specified key containing a suite of waveform signal-to-noise estimate attributes. extract_initial_beam_estimate uses the datum with the largest value of specified attribute as the initial stack estimate.

  2. The algorithm requires a definition of what I call the correlation time window. As noted above, the first stage of the multichannel algorithm is to align the data by cross-correlation with the current estimate of the stack (aka beam - a jargon term in array processing). For automated processing it is never a good idea to use a fixed window that is the duration of the input signals for two reasons: (a) variations in source properties (size, location, near source structure, and source complexity) drastically vary the optimal duration for correlation, and (b) interference from secondary phases (e.g. P or pP or even S for teleseismic P) is a nearly universal problem unless the initial data time span was carefully trimmed previously. The MsPASS solution to setting this time window uses a coda duration estimation algorithm appropriate only for teleseismic P wave data. That function is used within the higher level function called MCXcorPrepP discussed below. The algorithm used there automatically avoids interference from pP or PP phases and sets the end time of the correlation window based on a well-established coda decay algorithm. Specifically, the envelope of each signal is computed and the algorithm searches backward in time until the envelope exceeds an amplitude threshold based on a specified signal-to-noise ratio. In MCXcorPrepP the correlation window is derived from an average (multiple estimates of center are supported) coda duration from all ensemble members.

  3. The algorithm requires a different time window that in our original paper we called the robust time window. The robust time window is the time window used to define the robust stack I discuss below. The idea of the robust stack is to create a stack that automatically discards outliers and focuses on the better data while reducing the impact of low signal-to-noise data. Experience with dbxcor has shown that for reliable results the robust time window needs to be smaller than the time window defined for the correlation time window. A rule of thumb I always suggested for interactive processing is to define the robust window as the first 2 or 3 cycles of the dominant frequency of the phase being analyzed. For P wave data MCXcorPrepP standardizes this rule of thumb in the private helper _compute_default_robust_window. That algorithm is more-or-less a generalization of the “rule of thumb” I noted for dbxcor. It sets the robust window as a specified number of cycles of an input frequency. I recommend that normally be set as the low frequency band edge computed by broadband_snr_QC. An alternative appropriate for most data is a fixed time window spanning a time range from a small negative number large enough to exceed the maximum residual from earth model time estimates (typically from -2 to -1 s) to a few seconds (2 or 3 times the duration of the center frequency of the traditional short-period band). A fixed window performs well on smaller events that have significant signal only in the short period band but will work badly on large events. As a result an alternative is to process the data with different robust window lengths for different magnitude ranges.

Teleseismic P-wave automation#

For teleseismic P wave data the above automated algorithms are encapsulated in a single function called MCXcorPrepP. It provides a top-level interface for automatically preparing the filtered ensemble and initial beam required by align_and_stack. It returns those two objects as (ensemble, beam). The parameters this function defines are discussed in the section above. See the docstring for the function for guidance on use and examples below and in the mspass tutorial repository.

Robust stacking#

align_and_stack has several tunable parameters that were fixed constants in the original dbxcor program. To understand the context it is helpful to show the full function signature for align_and_stack:

def align_and_stack(
  ensemble,
  beam,
  *args,
  correlation_window=None,
  correlation_window_keys=["correlation_window_start", "correlation_window_end"],
  window_beam=True,
  robust_stack_window=None,
  robust_stack_window_keys=["robust_window_start", "robust_window_end"],
  robust_stack_method="dbxcor",
  use_median_initial_stack=True,
  output_stack_window=None,
  robust_weight_key="robust_stack_weight",
  time_shift_key="arrival_time_correction",
  time_shift_limit=2.0,
  abort_irregular_sampling=False,
  convergence=0.01,
  residual_norm_floor=0.1,
  demean_residuals=True,
  handles_ensembles=True,
  checks_arg0_type=True,
  handles_dead_data=True,
  **kwargs,
  ) -> tuple:

Note

The current implementation accepts use_median_initial_stack and abort_irregular_sampling, but does not forward either value to the underlying stack or sampling-validation calls. Consequently, the robust stack currently always starts from a median stack, and irregularly sampled members follow the default behavior of being logged and killed. Do not rely on setting either argument to change those behaviors until the implementation is updated.

The parameters of note are:

  1. The robust stack currently starts from the median stack of the ensemble aligned to the beam computed in the previous iteration. This is a stable starting estimate. As noted above, use_median_initial_stack=False is accepted by the API but does not currently select the previous beam as the starting estimate.

  2. residual_norm_floor implements a concept not recognized when we developed the original dbxcor program. It relates to a subtle feature of the robust weighting scheme we discussed in the original paper but we didn’t realize then how the concept this parameter implements would impact the results. In dbxcor it was fixed constant. To understand its use this is the weight formula used in align_and_stack that is enabled when robust_stack_method is set to the (default) of “dbxcor”:

\[w_i = \frac{1}{\| \mathbf{r}_i \|} \frac{\mid \mathbf{b} \cdot \mathbf{d}_i \mid} { \| \mathbf{d}_i \|}\]

where \(\mathbf{b}\) is unit-normalized and \(\mathbf{r}_i = \mathbf{d}_i - (\mathbf{b} \cdot \mathbf{d}_i )\mathbf{b}\). As we noted in the original paper there can be an issue for very consistent data if \(\| \mathbf{r}_i \|\) gets too small. A floor on that value is required, for example, to avoid a divide-by-zero error with simulation data that have no variance. There is a more subtle issue that we now know has a major impact on the result. To understand why it is helpful to note something we didn’t recognize when the dbxcor paper was published. That is, the second term in the weight formula, \(\frac{\mid \mathbf{b} \cdot \mathbf{d}_i \mid}{ \| \mathbf{d}_i \|}\), should be understood as the peak value of the cross-correlation function between the \(i^{th}\) datum with the beam (stack). Hence, the dbxcor weight function is the product of a cross-correlation weight and the inverse of the residual norm term. As we noted in the original paper the residual term makes the weighting more aggressively downweight any datum that differs significantly from the stack. That is desirable and a reason this algorithm can handle wildly variable quality data. The dark side to recognize, however, is that it makes the result strongly history dependent. The current median-started behavior produces a stack that is close to the median stack. How “close” is controlled by the setting of residual_norm_floor. When residual_norm_floor is 1.0 the residual weighting term is disabled. As the floor becomes smaller, the residual term more aggressively concentrates weight on signals closest to the evolving, median-started stack. The default 0.1 is appropriate for quality data. For ensembles with a large fraction of marginal signals a smaller value may be appropriate. See dbxcor_weights for the implemented floor definition and weighting details.

Examples#

This example shows the skeleton of a serial job reading from data that were previously preprocessed to bundle data into Seismogram objects. It is a “skeleton” as it shows the typical steps to process data with the multichannel correlation algorithm, but is far from complete and untested. The idea is you can use this as a starting point to work with your data set:

def prep_ensemble(e):
  """
  Applies required preprocessing to an input SeismogramEnsemble.
  This function is specialized to this workflow with
  function call arguments fixed.   A more generic function
  would use kwargs to change some arguments.
  """
  qcnw = TimeWindow(-120.0,-5.0)
  qcsw = TimeWindow(-5.0,60.0)
  for i in range(len(e.member)):
    # d is shorthand used for readability of this example
    d = e.member[i]
    if d.dead():
      continue
    d = rotate_to_standard(d)
    # assume default handling of input slowness using Metadata
    # attributes ux and uy set previously
    d = free_surface_transformation(d,vp0=5.0,vs0=3.5)
    d = broadband_snr_QC(d,
      noise_window=qcnw,
      signal_window=qcsw,
      use_measured_arrival_time=True,
      measured_arrival_time_key="Ptime",
    )
    # align the relative time origin to the measured P arrival
    d = ator(d, "Ptime")
    e.member[i]=d
  return e

######################### MAIN ############
# assume db is Database handle object
site_matcher=ObjectIdMatcher(db,
  collection='site',
  attributes_to_load=['lat','lon','elev','_id'])
source_matcher=ObjectIdMatcher(db,
  collection='source',
  attributes_to_load=['lat','lon','depth','time','_id'])

srcids = db.wf_Seismogram.distinct('source_id')
nw = TimeWindow(-100.0,-5.0)
# code above would define a query and run find to generate cursor
for sid in srcids:
  cursor = db.wf_Seismogram.find({'source_id' : sid})
  ens = db.read_data(cursor,
    collection='wf_Seismogram',
    normalize=[site_matcher,source_matcher],
    )
  ens = prep_ensemble(ens)
  ens = ExtractComponent(ens,2)
  ens, beam = MCXcorPrepP(ens,
    nw,
    station_collection='site',
  )
  if beam.dead():
    continue
  ens, beam = align_and_stack(ens, beam)
  db.save_data(ens,collection='wf_TimeSeries',data_tag='MCXcorProcessed')
  db.save_data(beam,collection='wf_TimeSeries',data_tag='MCXcorProcessed')

Local Earthquake Arrival Time Measurement#

Fundamentals#

Processing local earthquake data has similarities to processing teleseismic data, but there are some fundamental differences that require very different handling:

  1. Local earthquake signals can only be processed by coherent signal processing methods like that used in align_and_stack for dense arrays of stations. For most of the Earth that means a group of instruments with an aperture of the order of 1 km or less.

  2. Common receiver gathers of sources located in a small volume of the order of a few square km can sometimes be processed with coherent processing like that used in align_and_stack. In practice that tends to work only for events of a similar size and focal mechanism.

An under-appreciated, in my opinion, fundamental property of the Earth is that local earthquake signals are so fundamentally different from teleseismic signals. The reason is not actually known but the prevailing model is that waves with frequencies over around 1 Hz are more strongly scattered than the lower-frequency signals that define teleseismic body wave phases. There is indirect evidence that the crust and near-surface are the main source of the stronger scattering but that may be an observational gap due to the fact that earthquakes are not observed below depths of about 700 km. In any case, from a data processing perspective, coherent wavefield processing of local earthquake data is rarely feasible and other approaches have proven more successful. The methods currently available can be grouped into two broad categories:

  1. Traditional analyst-based picking. In the dark ages of seismology that means times measured from paper records. In the digital data era that has always meant manual picks made from a computer screen with a graphical user interface. That approach has been the bread-and-butter of seismic network operations worldwide for decades.

  2. Automated waveform processing methods using some form of detector and associator. By detector I mean an algorithm like the tried and true sta/lta detector that has been the standard since the earliest days of digital seismic network operations. A detector flags a signal transient based on looking at one and only one channel of data. An associator is an algorithm that collects a group of detector outputs and makes a decision on which detections are possible seismic events that should be processed to estimate a possible location. Many modern associators have an integrated, fast preliminary earthquake location algorithm. The net output is a preliminary earthquake location and the set of detections that are consistent with that location. In almost all cases a large fraction of detections are discarded as spurious by a good associator.

The detector/associator paradigm is central to local and regional seismic network operations. As stated multiple times in this User Manual, MsPASS was not designed to address the seismic-network operation problem. Multiple, robust software systems exist to address this problem both in the private and open-source worlds. MsPASS views this as a solved problem best handled by other tools if done on a large scale. On the other hand, I know from experience a typical research problem may need to do some form of manual or automated picking of local earthquake data. The remainder of this section addresses how MsPASS can be used in combination with some other packages for addressing that issue.

ObsPy Picking#

ObsPy does not provide a complete graphical analyst workstation or a general earthquake associator/location system. It does provide single-station triggers and pickers plus a network coincidence trigger. These are best understood by reading the Trigger/Picker Tutorial. Like MsPASS, ObsPy is a research toolkit rather than a complete seismic-network operations system, so association and location normally require additional software.

Monitoring Network Software#

A large fraction of research problems in seismology may need to use one of the specialized packages used for seismic network operations. At present the ones I know of are:

  1. In the US scientists at academic institutions can obtain a license for the commercial package called antelope provided they are not the operators of an operational seismic network. Antelope has a full suite of network processing capabilities. For research applications a particularly valuable feature is its capability to operate the real-time system on previously recorded data. Antelope uses a nonstandard relational database called Datascope that uses ascii files to hold the relational database tables. Because I am a long term user of Antelope I developed a specialized python class for handling these tables called DatascopeDatabase. It can be used to import and internally manage antelope database tables. How to do that is discussed in the Importing Tabular Data section of this manual. Using Antelope for processing is a larger topic addressed in Antelope’s documentation.

  2. Earthworm is an open-source earthquake monitoring system. Often as the processing core of AQMS, it is widely used by regional and volcano monitoring networks in the U.S. Earthworm builds on several applications originally developed by the U.S. Geological Survey that have been around for decades. A typical example is Hypoinverse. Earthworm releases and development are maintained with support from ISTI. Earthworm is organized as communicating modules and does not provide a single MsPASS-compatible database interface. Integration therefore depends on the archive or export modules used by a particular installation. MiniSEED waveform data can be indexed normally by MsPASS; other Earthworm products may require custom import/export adapters.

  3. seiscomp is more-or-less the European equivalent of Earthworm. It has similar functionality for real-time data acquisition, detection, event association, and location. It also has a collection of graphical user interface tools useful for network operations. I have zero experience with seiscomp but the same statement I made about Earthworm applies: this is a place someone from the Seiscomp community could help MsPASS development by producing import/export functions from that system.

MsPASS PhaseNet Picking#

MsPASS provides annotate_arrival_time as a thin wrapper around a SeisBench waveform model. It operates on a TimeSeries, applies a supplied model (or loads SeisBench’s pretrained PhaseNet stead model), and writes the P-wave probabilities that meet the requested threshold to the datum’s p_wave_picks Metadata value. That value is a dictionary mapping UTC timestamps to probabilities. An optional absolute-time TimeWindow can restrict the annotation interval. The function mutates the input datum and does not return a separate picks object.

SeisBench is an optional dependency. Install MsPASS with the seisbench extra when this feature is needed, for example pip install 'mspass[seisbench]'. The wrapper annotates candidate P-wave times; it is not an event associator or earthquake-location workflow.