Introduction

Seismology provides an attractive tool to investigate physical processes in deforming systems. The seismic signals from active glaciers, for example, could enable monitoring of mechanisms including basal sliding [6, 13, 24], fracturing [9], melt water drainage [4], and iceberg calving [14, 15]. The detection of seismic events from the recorded continuous seismic waveform data is a vital first step in any analysis. Event catalogues thus constructed are needed for local seismicity studies, comparisons between locations, or detection of change over extended time periods. However, the automated detection of seismic events is complicated in environmental and geotechnical seismology by the diverse populations of signal generation mechanisms. Those generated by active glaciers, for example, can be expected to span several orders of magnitude in duration and amplitude [16]. Volcanoes [10], landslides [19] and mining activity [25] similarly produce a broad range of seismic signals. As a result, the majority of cryoseismology studies to date use manual identification of events [6, 17, 18]. Manual techniques are not readily scalable nor exactly reproducible: in particular, they are not a first choice for monitoring applications nor for the data-driven detection of change.

Established event detection algorithms have largely been developed for earthquake seismology [1, 2, 3, 7]. These algorithms, including STA/LTA (short-term average/long-term average) [1] and template matching [3], are applied in real-time to seismic data to detect earthquakes and produce event catalogues [2, 7, 8, 23]. The core classes in the obspy software package are therefore designed assuming event metadata is available online alongside the waveform data [22]; additional functions are included outside the core classes to select events directly from waveform data. These standard algorithms are generally only applicable to non-earthquake signals by employing an experimental approach to parameter selection [5, 16]. The seismic attributes library provides software tools to download seismic waveform data from online repositories and detect environmental and geotechnical seismic events using a choice of algorithms. Algorithm options include the classic, recursive and delayed STA/LTA algorithms [22], and the newly developed multi-STA/LTA algorithm [11]. The multi-STA/LTA can simultaneously extract both short and long duration events of very different signal-to-noise levels and potentially enables real-time monitoring of cryoseismic events.

Event catalogues constructed using the seismic attributes library would be expected to comprise diverse signals generated by a range of mechanisms as noted above. The various signals typically need to be separated into related clusters based on the characteristics of their waveforms in order to study the events further [10, 19]. For example, Provost et al. [19] consider 71 attributes based on waveform data in their study of landslide seismicity. These are broadly split into four categories: (1) waveform attributes (e.g. duration, energy, kurtosis); (2) spectral and spectrogram attributes (e.g. discrete Fourier transform); (3) network attributes (e.g. station with maximum amplitude); and (4) polarity attributes (e.g. azimuth, inclination). The seismic attributes library includes functions to calculate a number of standard signal properties: duration, ratio between ascending and descending time, energy in the autocorrelation function, energy in the frequency filtered spectrum, and the direction of wave propagation. These are provided in three bundles of attribute functions describing the waveform, spectrum and polarity of the signal (Table 1); we do not include network attributes as these are unordered, discrete variables (and largely application dependent). User-defined functions can be added to derive customised characteristics within our software architecture. The correlation between attributes can be investigated using a plotting function to inform the removal of redundant variables, if appropriate, for the subsequent clustering application to a given set of events. Data-driven techniques more generally, such as machine learning algorithms, may be readily applied to the calculated attributes to inform the current and future state of the glacier (i.e. identify signals in the lead-up to a large event [20]).

Implementation and architecture

The analysis of waveform data from environmental and geotechnical seismic deployments requires multiple distinct steps. Following the principle that one should write programs that do one thing well, and write programs that work together [21], the seismic waveforms library is correspondingly split into three primary functions to streamline the workflow. The name and purpose of these functions are as follows:

• get_waveforms(); this function downloads waveform data from an online repository or alternatively reads data stored locally.
• get_events(); this function uses an event triggering algorithm to produce an event catalogue based on the seismic signals in one or more components of one or more seismic stations.
• get_attributes(); this function produces a pandas DataFrame of the attributes for each event using functions included in the library, or user-defined attribute functions.

The library includes the numerous attribute functions that are called by get_attributes(). Additional functions are also included to analyse the output of the primary functions (e.g. plotting), together with several private functions that support the primary workflow (documented in the code itself). The implementation of the three primary functions is outlined below.

get_events

As noted in the introduction, several algorithms exist to trigger events from single component waveform data. The get_events code first separates the signals from multiple seismometers and combines their (usually) three components into a single waveform. The user can optionally specify whether the component waveforms are combined as the sum-of-squares of their amplitude (i.e. Euclidean norm) to give the wave energy, or the absolute value of the wave amplitude. Further, to ensure that taking the absolute value of the amplitude does not affect the results (e.g. doubling frequency) we tested a computationally intensive option to fit the (time-varying) principle component of the direction of wave oscillation to obtain a signed amplitude; this gave identical event detections to the absolute amplitude and so is not included in our published version of the code due to the significantly longer computation time.

We use an adapted version of the obspy coincidence_trigger function to create a reference catalogue of events based on the STA/LTA characteristic function at each seismometer for their combined component waveforms. Our version of the coincidence_trigger function includes adjustments in the algorithm to better align with the outputs needed in our workflow and improvements to computational performance that are possible for our narrower use of this function (see Figure 1 for schematic of algorithm). This function can use the standard STA/LTA algorithms available within obspy [22], and optionally the recently developed multi-STA/LTA algorithm [11]. In all cases, location-specific algorithm parameters are likely to be more successful for environmental seismology applications. The start and stop times of the detected event from a single seismometer can be extended to consider small gaps in between triggers (up to length thr_event_join), creating a single, longer-duration record. The reference event catalogue is created by finding times when n (i.e. thr_coincidence_sum) of the N seismometers in the array have temporally coincident records. The physical dimensions of the seismic array are calculated at this step to estimate an upper bound on the delay in arrival time of the signal at different seismometers. This delay is added to the duration of the single seismometer records to ensure the reference events are not artificially shortened. The reference event catalogue, which includes a reference start time and duration for each event, is output as a pandas DataFrame to provide a convenient format to write to file and interrogate.

Figure 1 

Top: Diagrammatic representation of algorithm used to identify events from a single seismometer comprising one or more channels (traces show their Euclidean norm). The characteristic function for an STA/LTA algorithm is used to ‘trigger’ events (shown in red) and small gaps between these events (shorter than thr_event_join) are ignored (orange). No events are present at other times (shown in green). Bottom: Representation of algorithm used to identify events in the ‘reference event’ and ‘trace’ catalogues for seismic arrays (with multiple seismometers). The example shows an indicative array with N = 3 seismometers and thr_coincidence_sum = 2 simultaneous detections. Events identified for the single seismometers traces are shown in red (as given in the top panel). The duration of these events is extended at each end by half the delay in arrival time of a wave between the most distant seismometers in the array (shown in deep red). The reference event is identified as the times when thr_coincidence_sum = 2 seismometers have a detection (shown in red or deep red); small gaps with fewer seismometers are joined over (shown in orange). The events at each seismometer are simply the events identified from the single seismometer traces (red but not deep red), but with any times between events that occur during the reference event also included (shown in grey).

Our adapted version of the coincidence_trigger function provides additional utility for other researchers by outputting a secondary catalogue of trace metadata, for each seismometer, for all identified reference events. The triggered records at each seismometer that occur during (at least) part of the time period covered by a given reference event are joined together (if not already a single record). This record may extend beyond the time period of the reference event, be contained entirely within it, or have no detection at all in the case of weaker events. This trace (metadata) catalogue, which includes the trigger start time and duration for all seismometers, is similarly output as a pandas DataFrame. The trace metadata can be used to extract waveform data for identified events in the catalogue in a format that will be familiar to seismologists, for example, as shown in Figure 2 for a low-frequency event detected using four seismometers on the Whillans Ice Stream in Antarctica from 13:35:37 on 16 December 2010.

Figure 2 

Waveforms of an event detected using four seismometers (BB01, BB03, BB04 and BB06) comprising part of a seismic array on the Whillans Ice Stream in Antarctica from 13:35:37 on 16 December 2010. The top, middle and bottom panels show the vertical, north and east component of the signal respectively. The waveforms for each seismometer start 30 seconds prior to the start time in the trace catalogue and terminate 60 seconds after the stop time. The code to reproduce this plot is provided as a worked example.

get_attributes

The get_attributes function is written following further principles [21], that one should write flexible and open programs. This function extracts the waveform data for a given seismometer for the duration of a given event listed in the reference event, or trace (metadata) catalogue. The waveforms for each component are stored as separate Trace objects in a single stream. This Stream object is passed to one or multiple attribute functions to derive the characteristics of the given event as measured by a given seismometer. The values of the requested attributes for the array of seismometers are output as a pandas DataFrame. The correlation between the spectral attributes included in the library are shown in Figure 3 for events detected at Ilulissat, Greenland (DK.ILULI) on 1 January 2018 using the recursive STA/LTA algorithm.

Figure 3 

Corner (or pair) plot illustrating the correlation between spectral attributes for the events detected at Ilulissat, Greenland on 1 January 2018 using the recursive STA/LTA algorithm. The attribute names are defined to match those used by Provost et al. [19]; see Table 1 for a description of each attribute. In this example, attribute 26 considers a frequency range close to the sampling rate of the recorded signal, therefore, it has null value for all events. The code to reproduce this plot is provided as a worked example in the detailed documentation provided with the software.

Optionally, custom attribute functions may be added to calculate any chosen characteristic of a waveform. Custom functions must take a Stream object containing one or more components as an input, and output the attribute name and value. Public functions to combine the component waveforms into the wave energy or absolute value of the wave amplitude are included in the seismic attributes library to aid the user. The attribute functions (as many as are required) are passed to the get_attributes function as optional parameters, with return values stored as above.

Quality control

We have tested the seismic attributes library to find any software bugs and to ensure outputs are reproducible by other researchers running the code or using different platforms [12]. We first created test cases based on glacier seismic signals from a field campaign seismic array on the Whillans Ice Stream in West Antarctica [24]. Plausible uses of the code were tested by selecting some or all of the seismometers in the array, some or all of the three components of the signal recorded at each seismometer (i.e. vertical, north and east), and considering different durations of requested time (with and without small data gaps). The optional function inputs were also tested in this manner. In application, e.g. [11], the event catalogue produced using the seismic attributes workflow compared favourably to events that were visually classified using the same set of seismometers [18]. Finally, the attribute functions included in the library were tested on mock waveform data to ensure they produce the theoretically expected results.

Testing was primarily carried out using a Python 3.8 installation on a macOS 10/11 system. The key features of the seismic attributes library were tested on other platforms for compatibility, including both Python 2.7 and 3.7/8 for the most up to date version of obspy, numpy, pandas and seaborn available for that installation. The code was further tested on both macOS 10/11 and Microsoft Windows operating systems, especially to verify functionality of the different file systems. The reference event and trace (metadata) catalogues were compared between these platforms to assess the consistency of code functionality.

Operating system

macOS 10/11, GNU/Linux and Microsoft Windows.

Programming language

Python 2 or 3.

Additional system requirements

Memory and disk space will limit the volume of data that can be processed in a contiguous chunk; the code will work on any system that can support obspy. Internet access is required to download new waveform data from online repositories.

Dependencies

The minimal dependency for use is Python 2/3 with obspy, numpy, pandas and seaborn packages installed.

Software location

Archive

Name: An ObsPy library for event detection and seismic attribute calculation: preparing waveforms for automated analysis

Persistent identifier: https://doi.org/10.5281/zenodo.5496794

Licence: GPL version 3

Publisher: Zenodo

Date published: 10/09/2021

Code repository

Name: seismic_attributes

Persistent identifier: https://github.com/rossjturner/seismic_attributes

Licence: GPL version 3

Date published: 10/09/2021

Language

Python 2 or 3.