- © 2014 by the Seismological Society of America
Online Material: Script and input files illustrating the TDD package.
Many types of seismic analysis require ground‐motion data in physical units of displacement, velocity, or acceleration. Because seismometers do not respond uniformly to all frequencies, scientists must deconvolve the instrument response from recorded data. However, popular deconvolution methods violate causality and introduce acausal artifacts to deconvolved traces, which are unacceptable in some applications. In this paper, we discuss a new method using recursive filters that we consider superior for many purposes.
One simple and common method of instrument response deconvolution is spectral division by the analog instrument response (1)in which S is the instrument’s sensitivity, K is the normalization constant, and P and Z are the poles and zeros. This method should be avoided because it fails to consider effects of digitization on the instrument response. Failing to consider these effects results in inaccurate phase and amplitude spectra, and can introduce acausal precursory oscillations to the signal.
Another common method is that of the TRANSFER command in the widely used software package SAC (University of California, 2012). This command uses a zero‐phase high‐pass or band‐pass filter (through the freqlimits option) to ensure the stability of the result. This filter is acausal and, therefore, problematic as well.
Acausal artifacts can cause serious analytic problems. One obvious case is picking wave‐arrival times and polarities. If causality is not preserved in deconvolution, phases may appear to arrive earlier in deconvolved data than in raw data, causing picks to be inaccurate and biased. Additionally, some applications analyze small signals that are partially overwritten by later‐arriving, higher‐amplitude (even clipped) phases; we show an example of this in the Deconvolution of Partially Clipped Data section. Acausal precursory oscillations of the stronger signal can contaminate the weaker signal, reducing the quality of the results.
One new alternative involves calculating an approximate discrete instrument response with the bilinear transform and spectrally dividing the response from the data. Discretization effects make the spectrum of the discrete signal diverge from that of an analog signal at higher frequencies; however, this effect can be mitigated by oversampling before deconvolving (in practice, nearest‐neighbor interpolation works well) followed by decimation. Oversampling is necessary when the upper limit of the frequency band of interest exceeds about one fifth of the Nyquist frequency (Haney et al., 2012).
This method is a major improvement over ordinary spectral division because it accurately deconvolves instrument responses while maintaining causality. However, frequency‐domain methods are inefficient and cumbersome in real‐time applications in which data must be deconvolved as they arrive. Automated, real‐time analysis of incoming data is becoming increasingly feasible and important; in some applications like Earthquake Early Warning systems, lives may depend on rapid and accurate data analysis.
We propose recursive filters as an alternative to deconvolution by spectral division. These are appealing because they are necessarily causal, run quickly with little memory required, and can be performed as data arrive in real‐time applications like Earthquake Early Warning.
For each sample, recursive filters define the current output value yn as a linear combination of the current input value xn and past input and output values: (2)
Because yn is calculated without consideration of later terms in x, the filter is causal. The sequences a (autoregressive coefficients) and b (moving‐average coefficients) control the response of the filter. Suitable selection of a small number of filter coefficients can produce filters with arbitrary, causal responses (Karl, 1989).
Recursive filters have previously been used for real‐time analysis of earthquake data. They have proven useful for generating Wood–Anderson seismograms (Kanamori et al., 1999), and removing instrument responses from ultra‐long‐period W phases (Kanamori and Rivera, 2008). The approach used in these studies was to approximate seismometer responses using a gain, damping coefficient, and natural frequency, and to construct a second‐order recursive filter using these coefficients. This approach is suitable for some simple short‐period seismometers such as the Wood–Anderson. It can work for broadband instruments as well when the frequency band of interest is low enough that it excludes high‐frequency corners and other structures, but it cannot accurately process the full response spectrum.
We have generalized this approach to calculate recursive filters corresponding to instrument responses over any frequency band, and have created the R package “TDD” (Anderson, 2013) to facilitate time‐domain deconvolution. TDD includes documented functions for performing time‐domain deconvolution and calculating filter coefficients for arbitrary instruments, along with precalculated recursive filters for common seismometers (see electronic supplement for demonstration). R provides a powerful open‐source environment for data processing, and we encourage its use in seismology with this and other relevant packages such as RSEIS (Lees, 2012). TDD uses the command filter from the signal package (Signal developers, 2013) to perform recursive filtering; similar functions are available in Python (http://www.python.org/, last accessed November 2013) and MATLAB (http://www.mathworks.com/products/matlab/, last accessed November 2013) (lfilter and filter in their signal toolboxes), so users of those languages will not find the routines difficult to translate.
In this paper, we describe the methods involved in implementing recursive filters, and apply them to relevant data. In addition to the applications presented in this paper, we believe that these methods would also be suitable for real‐time seismic‐data processing.
We discuss two methods of calculating recursive filter coefficients. The first method uses a finite‐difference approximation to derivatives to convert the differential equation describing the instrument response into a recursive filter, and then uses a nonlinear inversion to adjust the coefficients to optimize fit to the analog spectrum. The second is an adaptation of the bilinear transform method of Haney et al. (2012). For each seismometer and sample rate considered here, we calculated discrete responses using each method.
To determine the best method for each seismometer‐sample rate pair, we found the highest frequency for which the discrete and analog spectra match within 1%. By this measure, the finite‐difference approximation method outperformed the bilinear transform method in most of the responses studied here (72 out of 84). Each precalculated discrete instrument response given in the TDD package was calculated using the best method for that seismometer and sample rate.
Converting Sensor Poles and Zeros to Filter Coefficients
In the first discretization method, we calculate the coefficients of the recursive filter that corresponds to the seismometer’s response. Instrument responses are normally given in terms of poles (P), zeros (Z), normalization constant (K), and sensitivity (S), such that the seismometer response Y to an input ground‐motion velocity X is defined as (3)with Y(2πif0)=SX(2πif0) for some normalization frequency f0. This can be rewritten as the differential equation: (4)in which the coefficients β and α are calculated as the polynomial expansion of the numerator and denominator of equation (3).
We use a finite‐difference approximation for the ith derivative: (5)in which is the binomial coefficient. So, equation (4) can be rewritten as (6)
Then, autoregressive (a) and moving‐average (b) coefficients of a recursive filter in the form of equation (2) can be calculated: (7)(8)
However, some additional adjustment is required to obtain an accurate filter (described in the following section) because of effects of discretization.
Improving Filter Fit
When Δt is not infinitesimal, the response of the filter generated by this method does not exactly match that of the instrument. Small adjustments to input values can reduce this inaccuracy. Using a Markov Chain Monte Carlo (MCMC) routine (Aster et al., 2012), we tried adjusting each of the following: filter coefficients (a and b) of equation (2), differential equation coefficients (α and β) of equation (4), and poles and zeros (P and Z) of equation (3). We determined that varying the poles and zeros found an accurate solution most quickly and reliably.
The MCMC model consists of the poles and zeros of the instrument response. For paired poles or zeros of the form a±bi, a and b are varied independently so that the product of the paired parameters remains real. In each iteration, a model parameter is selected at random, and a random perturbation is added to it. The perturbation is drawn from a normal distribution centered at zero with standard deviation equal to some factor times the starting value of that parameter; we found that setting that factor between 0.1 and 1 generally works well, although trying multiple factors is occasionally necessary to find a good fit.
This MCMC implementation minimizes the misfit between the analog amplitude spectrum of the instrument response and the amplitude spectrum of the filter’s response. (Explicit consideration of the phase spectrum here is unnecessary because the phase spectrum of these responses is strictly dependent on the amplitude spectrum; therefore, once the amplitude spectra match, so do the phase spectra.) In each iteration, the perturbed poles and zeros are converted into differential equation coefficients, which are in turn converted to recursive filter coefficients. The filter is then applied to an impulse, and the sum‐of‐squares misfit is calculated between the amplitude spectrum of the filtered impulse and the true instrument amplitude spectrum. Each perturbation to the model is accepted or rejected with a probability dependent on its effect on this misfit. Because of the effects of discretization, matching frequencies higher than about half the Nyquist frequency is sometimes impossible; therefore, we ignore misfit at frequencies above that.
For each sensor, we matched frequencies down to 0.1 times its low‐corner frequency. This means that in each iteration of the MCMC, we must calculate the FFT of a filtered impulse whose duration is the inverse of that low frequency. As a result, calculating recursive filters can be a time‐consuming process, especially at high sample rates for instruments with very low corner frequencies. Consequently, we precalculated recursive filters for the following common seismometers: Streckeisen STS‐1 360s and 20s, and STS‐2 (generations 1–3); Nanometrics Trillium 40, 120, and 240 (generations 1 and 2); and Guralp CMG‐3T, CMG‐3ESP, and CMG‐40T 30 and 1 s. Coefficients for recursive filters are particular to a sample interval; for each instrument, we calculated filter coefficients corresponding to 1, 0.1, 0.05, 0.025, 0.02, and 0.01 s. We did not calculate coefficients for customizable short‐period sensors; however, coefficients for short‐period sensors can be calculated fairly quickly given the damping constant, natural frequency, and sensitivity.
Calculating Filter Coefficients with the Bilinear Transform
In the second discretization method, we use the bilinear transform to convert the analog poles and zeros of the Laplace transform from equation (3) to K digital poles and zeros of the Z‐transform. This is essentially a time‐domain implementation of the method of Haney et al. (2012). It is noteworthy that, unlike the finite‐difference approximation, the number of digital poles and zeros returned by the bilinear transform is generally different from the number of poles and zeros in the analog instrument response.
The discrete response can then be written: (9)
Filter coefficients may then be calculated by multiplying both sides by the denominator of the fraction, expanding the polynomials, and performing the inverse Z‐transform.
Implementation of Recursive Filter
Once filter coefficients are calculated, they may be applied to data. They can be used directly for convolution (for obtaining Wood–Anderson responses, for example). However, the autoregressive and moving‐average coefficients must be swapped for deconvolution. Additionally, high‐pass filtering is required to make the response stable. Our deconvolution routine implements a customizable Butterworth band‐pass filter by concatenating its digital poles and zeros to those of the deconvolution filter. Deconvolution and band‐pass filtering are both implemented using a recursive filter, so the output of this process remains causal. For most seismometers, both the finite difference and bilinear approximations become unreliable at high frequencies (Fig. 1). Oversampling can improve deconvolution accuracy in cases for which the frequency band of interest extends above the maximum frequency that can be reliably deconvolved, but it also increases runtimes.
Efficiency of Recursive Filters
Time‐domain deconvolution is the fastest means of removing instrument responses from long time series. Recursive filters run by iterating through the terms in the input data sequence and calculating an output term for each. Every iteration requires the same number of operations, so the filter theoretically runs in linear time . Only a small number of terms (equal to the total number of poles and zeros) need to be considered in each iteration, so the memory required is only slightly more than that needed to store the input and output traces. Additionally, seismic data are typically real, so no expensive complex number operations are required.
In contrast, traditional frequency‐domain methods, such as the one used in SAC, are of slightly higher complexity (Cooley and Tukey, 1965) and require complex operations, making them more expensive in terms of time and memory for longer data traces. Additionally, in frequency‐domain methods, doubling the length of time series by zero‐padding is necessary to prevent wrap‐around effects; this increases the expense as well.
We compared runtimes of deconvolution performed in the frequency and time domains. We calculated the median runtime of 100 trials of each method for several time‐series lengths; results for the CMG‐3T are shown in Figure 2. To ensure that these tests measured the speed of the deconvolution method and nothing else, we precalculated instrument responses and timed only the deconvolution, using minimal R code, without subsequent filtering. For time series greater than about 104 samples, the recursive filter was faster than spectral division.
APPLICATIONS AND DISCUSSION
Demonstration on Data from Distinct Collocated Sensors
We compare data from collocated sensors of different types to demonstrate the effectiveness of these recursive filters (Fig. 3). These data are from a Guralp CMG‐3T and CMG‐40T‐1s, which were deployed in the same vault and digitized at 100 samples per second and 24‐bit resolution by a Reftek RT‐130 data logger. Instrument responses were deconvolved from recorded data using precalculated recursive filters and were high‐pass filtered above 0.1 Hz. These instruments have very different low‐corner periods (120 s for the 3T, 1 s for the 40T), and raw data from the two sensors reflect this difference. However, the two traces are almost indistinguishable after deconvolution.
The 10‐s cutoff period of the high‐pass filter is well within the passband of the CMG‐3T, but a factor of 10 below the low corner of the CMG‐40T‐1s. Despite a 40‐dB drop in sensitivity at this low frequency compared to its passband, its deconvolved records agree very well with those from the CMG‐3T. This excellent match between data from such different seismometers demonstrates the accuracy of these recursive filters.
Deconvolution of Partially Clipped Data
One exciting application of causal instrument corrections is the opportunity to use data that have been partially clipped. Seismic signals commonly arrive in distinct phases that arrive close to each other in time, but sometimes with different amplitudes; later‐arriving phases, like S waves, airwaves, and surface waves, are often of higher amplitude than earlier phases. Sometimes, these later phases clip whereas earlier data are completely intact. In these cases, causal deconvolution can be used to convert data to physical units of ground motion; unlike acausal deconvolution, artifacts of deconvolving off‐scale data will not affect any part of the signal occurring before the first instance of clipping.
It is likely that a wealth of such data exists that has not been fully exploited because of potential artifacts from acausal deconvolution. We show one example from Tungurahua volcano, Ecuador, whose explosions produce intense infrasonic signals with amplitudes up to 180 Pa as far as 3.5 km from the vent (Ruiz et al., 2006). In many explosions, these powerful airwaves overprint seismic signals, and some stations use short‐period sensors and recording systems that can clip upon their arrival. We show one such explosion in which a long‐period signal arriving before the airwave is visible in deconvolved traces, but not in raw data (Fig. 4); these long‐period signals can be used to understand seismic‐source processes.
Recursive filters can be used for time‐domain instrument corrections that preserve causality. Common methods of deconvolution, including that of the popular SAC package, introduce acausal artifacts, which are unacceptable in many seismic‐analysis methods. For long time series, recursive filters are faster than spectral division. Additionally, they are much better suited for real time use than frequency‐domain methods.
One important application of causal deconvolution is extracting useful information from seismograms that have been partially clipped by late‐arriving phases. Acausal deconvolution of instrument responses causes artifacts related to clipping to appear in the nonclipped part of the signal; therefore, instrument corrections of partially clipped data have traditionally been avoided. This is not a problem for causal methods of deconvolution, in which artifacts related to off‐scale recordings cannot occur prior to the onset of clipping. As a result, seismologists should consider using causal deconvolution methods to correct partially clipped traces up to the first clipped arrival.
We have released an R package (TDD) for time‐domain deconvolution and generation of instrument responses. This package includes precalculated discrete instrument responses for 14 common seismometers, calculated at six common sample rates. TDD integrates well with the general‐purpose seismology package RSEIS to provide a powerful, open‐source environment in which many routine seismic‐data processing tasks, including deconvolution, can be accomplished.
We thank H. Ortiz and the Instituto Geofisico–Escuela Politecnica Nacional (Ecuador) for providing seismic data and deployment information for the RETU station, and R. Aster, S. Ingate, and PASSCAL for providing data from the collocated CMG‐3T and CMG‐40T‐1s. We gratefully acknowledge a thorough and helpful review by M. Haney. This work was funded by NSF Grant EAR‐0838562.