MouseTrap module

MouseTrap is an open-source module for Python for detecting a special type of disturbance (so called mouse, ping, or fling step) in seismic records.

Copyright: Jiří Vackář.
Version: The most recent version is 0.2.1, released on 2015-04-28.
License: GNU Lesser General Public License, Version 3 (http://www.gnu.org/copyleft/lesser.html)

The code is described in paper Vackář et al., Seis. Res. Lett., 2015 (link, PDF)

Acknowledging

Please support the project by acknowledging the use of it. If you use MouseTrap for work resulting in an academic publication, we would be grateful if the following paper is cited:

Jiří Vackář, Jan Burjánek, and Jiří Zahradník.
Automated detection of long-period disturbances in seismic records; MouseTrap code,
Seismological Research Letters, 2015.

Download

  • MouseTrap.py module file
  • Example 0: creating a synthetic mouse
  • Example 1: mouse detection in single record, reading from file (needs record and poles-and-zeros files; the output should be this fit plot)
  • Example 1b: mouse detection in long record with no a priori assumption where the mouse should be (designed for detection in continuous records), reading from file (needs record and poles-and-zeros files; the output should be similar plot as before for each time window in which is record internaly divided)
  • Example 2: mouse detection in single record, reading via ArcLink
  • SwissMouse: massive detection of disturbances in a large dataset using SeisComP database (needs database structure: use SwissMouse.sql and module peak_values.py)

Function description

mouse.__init__(min_step = 2., min_step2 = 0.1, fit_time_before = 115, fit_time_after = 60)

Initialize class mouse. For using with default configuration, simply type:

my_mouse = mouse()
mouse.create(paz, length, dt, onset, typ=1, demean=True, integrate=True)

Calculates the instrument output for a given transfer function and a unit acceleration (variantly velocity or displacement) step on the input.

The transfer function is given by paz, which might come from ObsPy function sac.sacio.attach_paz or arclink.client.Client.getPAZ. The calculated time series has a length of length, which should be a power of 2, otherwise it is increased to next power of 2. Parameter dt specifies time interval between following two time samples. The step is placed onset seconds after the record start time. The input is velocity step for typ=2, displacement step for typ=3, and acceleration step for typ=1 or any other value (default). If the demean and integrate parameters are True, the functions mouse.demean() and mouse.integrate() are applied on the result, respectively. If the parameter integrate=True, the result, saved in the variable mouse.mouse, is displacement time series; integrate=False results in velocity time series.

mouse.demean()

Removes the mean value of the interval before mouse onset from entire synthetic mouse record.

Typically, it is not necessary to run it separately, it is called from mouse.create().

mouse.integrate()

Integrates the synthetic mouse record by rectangular rule.

Typically, it is not necessary to run it separately, it is called from mouse.create().

mouse.fit_mouse(st, t_min=0, t_max=0)

Fits a given three-component or one-component record by the synthetic mouse.

Synthetic mouse must be calculated before by mouse.create(). For three-component inversion, st is ObsPy stream object, which must contain three components with channels (st[i].stats.channel) named ??N, ??E, and ??Z (e.g. HHN, HHE, and HHZ). For one-component inversion, st is ObsPy trace object.

Parameters t_min and t_max affect the grid-search over time. If they are non-zero, they declares first and last tested time of mouse onset (in seconds from stream start-time). If t_min=0 and t_max=0, the first and last tested time is start-time and end-time of the record, respectively.

mouse.invert(N, E, Z, shift)

Analytically calculates best-fit mouse direction and amplitude based on partial derivatives (see Section Formulas behind).

It is called from mouse.fit_3D(), no need for using it alone.

mouse.invert_1D(D, shift)

Analytically calculates best-fit mouse amplitude based on partial derivatives in case of one-component inversion. It is simplified version of mouse.invert().

It is called from mouse.fit_3D(), no need for using it alone.

mouse.exist(Ts=120, likehood=False)

Distinguishes whether the mouse is present according the fit value and the difference between the mouse onset time and S-wave arrival time. The output is 2 if the mouse is present, 0 is the mouse is not present, and 1 if it is case in the middle of the way. Variantly, the output can be mouse likelihood value (see below).

Parameters:
  • Ts (float) – specifies the S-wave arrival time (second).
  • likehood (bool) – if set to True, returns mouse likehood value instead values described before.

Remark: Function uses some empirical constants chosen for one study. If you need to decide obaut mouse presence or absence, think about your own criterion based on amplitude and fit values and possibly time of the onset.

mouse.params(degrees=False)
Parameters:degrees (bool, optional) – return angles in degrees instead of radians

Returns tuple of the parameters and fit of the detected mouse: onset t_0 (second), amplitude A (\mathrm{m\,s^{-2}} for mouse type 1), azimuth \phi (radians, optionaly degrees), inclination \theta (radians, optionaly degrees), and fit (percent). If the angles were not inverted (in case of one-component inversion), returns None value for them, so the order of returned tuple does not change.

mouse.plot(st, outfile='', distance=5e3, ylabel='', xmin=100, xmax=300, legend=False, yaxes=True)

Plots the comparison between observed record (st) and synthetic mouse with fitted onset, amplitude and direction, so the function mouse.fit_3D() must be called first.

If outfile is specified, the figure is saved to this file, otherwise it is displayed on the screen. The parameter distance specifies the vertical distance between three traces in the figure (in the same units as the traces are). ylabel specifies the y-axis label, xmin and xmax range of x axis, legend turns on legend in the plot. yaxes turns on ticks and values on the y-axis and add y2-axis showing acceleration of the simulated step.

PrepareRecord(st, t, ratio_velocity=20, ratio_displacement=8, bitrate=10)

Analyses the signal-to-noise ratio of the record, removes the before-event mean value and integrates record into displacement.

It is shorthand for

if NoiseTest1_demean(st, t, ratio_velocity):
        return 'Record too noisy.'
ToDisplacement(st, bitrate)
if NoiseTest2(st, t, ratio_displacement):
        return 'Record too noisy after integration.'
NoiseTest1_demean(st, t, ratio=20)

Removes the before-event mean value from entire st stream. The before-event interval is specified by t parameter in seconds from start-time of the stream.

Then it tests simple signal-to-noise criterion. The ratio of before-event-maximum must be ratio times smaller than the record maximum. If the stream does not pass this test, the function returns True.

demean(st, t):

Calculates mean value from part of the records between starttime and the time t and removes this mean from entire stream st.

ToDisplacement(st, bitrate=10)

Integrates and downsamples the record st to given bitrate.

NoiseTest2(st, t, ratio=8)

Test simple signal-to-noise criterion. The ratio of before-event-maximum must be ratio times smaller than the record maximum. If the stream does not pass this test, the function returns True.

Formulas behind

Mathematical description of the disturbance

The synthetic mouse (forward modeling) in raw velocity m_v(t) is described by

m_v(t) = \int_{-\infty}^{t} s(t) \mathrm{d}t  \ast  h(t)     \; \rm{,}

where h(t) is an impulse response of the instrument to input velocity and s(t) is the input unit acceleration step, i.e. the Heaviside function (s(t)=1 for t\geq1 and s(t)=0 for t<0). The mouse in displacement is then

m_d(t) = \int_{-\infty}^{t} m_v(t) \mathrm{d}t

We use 4 parameters for description of a real mouse: time t_0 of the onset of the acceleration step, amplitude A of the acceleration step and two spatial angles – horizontal azimuth of acceleration pulse \phi and its inclination \theta from the horizontal plane. The east, north and vertical component of the displacement record are then

m_E(t) &= m_d(t-t_0) * A * \sin\phi * \cos\theta   \\
m_N(t) &= m_d(t-t_0) * A * \cos\phi * \cos\theta   \\
m_Z(t) &= m_d(t-t_0) * A * \sin\theta

Fitting the record by a synthetic disturbance

Fitting a record by the synthetic mouse formally means solving inverse problem with 4 parameters: t_0, A, \phi, and \theta. The inverse problem is solved by the least-squares method (LSQ) to minimize the L^2-norm difference between the observed record and synthetic mouse, assuming a value of t_0. Analytical expressions of the partial derivatives with respect to A, \phi, and \theta are used in the LSQ fitting. The proper value of t_0 is calculated by grid search. The fit is quantified by variance reduction (VR).

Let have a three component record with north-south, east-west, and vertical components s^N_i, s^E_i, and s^Z_i, respectively, where i index time samples. We want to minimize the difference between the record and a synthetic disturbance m_i by finding proper amplitude A, azimuth \phi, and inclination \theta of the disturbance. The difference in the L^2-norm should be minimal

\sum_i \left( s^N_i - A \cos \phi \cos \theta m_i \right)^2  +  \sum_i \left( s^E_i - A \sin \phi \cos \theta m_i \right)^2  +  \sum_i \left( s^Z_i - A \sin \theta m_i \right)^2  =  min     \; \mathrm{,}

so its partial derivatives should be zero. From \frac{\partial}{\partial \phi} = 0 we get

\frac{\sum_i m_i s^E_i}{\sum_i m_i s^N_i}  =  \frac{\sin \phi}{\cos \phi}  =  \tan \phi     \; \mathrm{.}

From \frac{\partial}{\partial \theta} = 0 we get

\frac{\sum_i m_i s^Z_i}{\sum_i m_i \left( s^N_i \cos \phi + s^E_i \cos \phi \right)}  =  \frac{\sin \theta}{\cos \theta}  =  \tan \theta     \; \mathrm{.}

From \frac{\partial}{\partial A} = 0 we get

\frac{\sum_i m_i \left( s^N_i \cos\phi\cos\theta + s^E_i \sin\phi\cos\theta + s^Z_i \sin\theta \right)}{\sum_i m_i^2}  =  A     \; \mathrm{.}

Examples

Basic usage

To start using MouseTrap, type:

from MouseTrap import *

To create synthetic mouse, we need poles and zeros for some instrument:

from obspy.core import AttribDict
paz = AttribDict({
        'sensitivity': 600000000.0,
        'poles': [(-0.1103+0.111j), (-0.1103-0.111j), (-86.3+0j), (-241+178j), (-241-178j), (-535+719j), (-535-719j)],
        'gain': 110400.0,
        'zeros': [0j, 0j, (-68.8+0j), (-323+0j), (-2530+0j)]
        })

Now, we can create synthetic a mouse:

m1 = mouse()
m1.create(paz, 4000, 0.1, 100, 1)

The m1 mouse has 4000 samples at sampling rate 0.1 second and the mouse onset is at time 100 second.

Finally, lets plot the result:

import matplotlib.pyplot as plt
t = np.arange(0, len(m1.mouse)/10., 0.1)
plt.plot(t,m1.mouse)
plt.show()

You can download source codes of some examples in Section Download.

SwissMouse

Example of use MouseTrap module for a massive detection of disturbances in a large dataset using SeisComP database.

For description, see flowchart of the algorithm and paper mentioned in section Acknowledging.

For running this example you need:
  • access to SeisComP database
  • access to ArcLink server (SED ArcLink server has free access, only set your e-mail address)
  • your own Postgre SQL database for saving results
  • specified structure prepared in your DB (use SwissMouse.sql)

Figures

_images/fit_VANNI.png

Example of a step-like disturbance in the integrated output (raw displacement) of the Nanometrics Trillium T40 seismometer. The disturbance is well fitted by the simulated instrument response to an acceleration step. There earthquake in the record has M_{Lh} = 1.7, its epicentral distance is 1.9 km. The acceleration step has amplitude 8.8\cdot10^{-7}\,\mathrm{m \cdot s^{-2}}.

_images/spectrum.png

Power spectrum of raw velocity (Z component in the previous figure) and of the disturbance. The disturbances cause significant spurious increase of the spectrum at low frequencies. It dominates up to frequency \sim0.12\,\mathrm{Hz}, that is much higher than the low-frequency corner of the instrument.

_images/mys_0_thin.png _images/mys_1_thin.png _images/mys_2_thin.png
_images/mys_3_thin.png

Modeling a disturbance. The acceleration step is equivalent to the velocity ramp. After applying the instrumental response we get a characteristic one-sided pulse in the output raw velocity, and a step-like disturbance in the raw displacement. Note that the step visible in the raw displacement is not the permanent displacement, but a permanent acceleration. It is because the integrated broad-band output is proportional to displacement only at frequencies above f_c (the low-frequency corner of the instrument), but it is proportional to acceleration at the low-frequency limit.

Indices and tables