Forensic Seismology and Acoustics¶
Notebook for class on February 17th.
Stephen Arrowsmith ([email protected])
Data Exercise 1: Reading and analyzing DPRK data¶
This first class introduces you to analyzing seismic data with Python. We're going to dive in the deep end and try to tread water. First thing we'll do is import some 'modules'. These are different code packages that contain different things. We can import whole modules, parts of modules, and rename modules when we import them.
The module called 'forensic' is built for this class. The module called 'obspy' is an important one for analyzing seismic and acoustic data.
import forensic
from obspy import UTCDateTime
from pyproj import Geod
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import datetime, pickle
The next thing we'll do is define an event location and time. We're just defining the location and time of the DPRK nuclear test in September 2017, using the values from the Myers et al. (2018) paper that we read already. Defining times can be tricky, ObsPy (the main seismic data module we read in above) likes times in a type of time format called 'UTCDateTime'.
evla = 41.3; evlo = 129.08; evdp = 0.08
evt0 = UTCDateTime(datetime.datetime(2017,9,3,3,30,1))
evt0
Reading data¶
You can read data seismic directly into Python with an internet connection using the forensic module like below. We'll skip over the details of how to do this for now, because it's complex, but it might come in handy later. For now, we'll open the file called 'DPRK.p' that you downloaded already.
#st = forensic.getStreamFromIRIS(evt0-600, evt0+3000, '*', 'IU,II,IC', 'BHZ')
#pickle.dump(st, open('DPRK.p', 'wb'))
st = pickle.load(open('DPRK.p', 'rb'))
Exploring ObsPy¶
To explore data in ObsPy, you need to understand how the data are structured in an 'object' called an ObsPy Stream: https://docs.obspy.org/packages/obspy.core.html#waveform-data
print(st)
tr = st[0]
print(tr.stats)
tr.data
%pylab notebook
plt.plot(tr.data)
%pylab notebook
#fig, ax = plt.subplots()
fig = tr.plot(handle=True)
for tr in st:
print(tr.stats.station, tr.stats.network,
tr.stats.channel, tr.stats.sampling_rate,
round(forensic.get_dist(tr, evla, evlo)/111.1,2))
%pylab notebook
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
for tr in st:
ax.plot(tr.stats.sac.stlo, tr.stats.sac.stla, 'y^', transform=ccrs.Geodetic())
ax.plot(evlo, evla, 'r*', transform=ccrs.Geodetic())
ax.gridlines(draw_labels=True)
plt.show()
Exploring the DPRK data¶
We'll plot 'record sections' which are plots that show multiple seismograms, arranged according to the location and origin time of an event (in this case, the September 2017 DPRK test).
%pylab notebook
forensic.plot_recordsection_global(st, evlo, evla,
date2num(evt0.datetime),
scaleFactor=1.5)
plt.xlim([0., 160.])
plt.ylim([0., 50.])
%pylab notebook
st_f = st.copy()
st_f.taper(type='cosine', max_percentage=0.001)
st_f.filter('highpass', freq=1)
forensic.plot_recordsection_global(st_f, evlo, evla,
date2num(evt0.datetime),
scaleFactor=1.5)
plt.xlim([0., 160.])
plt.ylim([0., 50.])
%pylab notebook
st_f = st.copy()
st_f.taper(type='cosine', max_percentage=0.001)
st_f.filter('lowpass', freq=0.1)
forensic.plot_recordsection_global(st_f, evlo, evla,
date2num(evt0.datetime),
scaleFactor=1.5)
plt.xlim([0., 160.])
plt.ylim([0., 50.])
Estimating mb¶
Following the new IASPEI standards for magnitude estimation (Borman and Dewey, 2012), if $A$ is calculated in nm (Peak-to-Peak) and $T$ is calculated in s, the body wave magnitude can be calculated from:
$ mb = \log_{10} (A/T) + Q(\Delta, h) - 3.0$
Where:
$A$ = P-wave ground amplitude in nm calculated from the maximum trace-amplitude in the entire P-phase train (time spanned by P, pP, sP, and possibly PcP and their codas, and ending preferably before PP);
$T$ = period in seconds, T < 3 s; of the maximum P-wave trace amplitude.
$Q(\Delta, h)$ = attenuation function for PZ (P-waves recorded on vertical component seismographs) established by Gutenberg and Richter (1956a) in the tabulated or algorithmic form as used by the U.S. Geological Survey/National Earthquake Information Center (USGS/NEIC)
$\Delta$ = epicentral distance in degrees, 20° ≤ $\Delta$ ≤ 100°;
$h$ = focal depth in km;
and where both $T$ and the maximum trace amplitude are measured on output from a vertical-component instrument that is filtered so that the frequency response of the seismograph/filter system replicates that of a WWSSN short-period seismograph, with $A$ being determined by dividing the maximum trace amplitude by the magnification of the simulated WWSSN-SP response at period $T$.
We'll begin by plotting the attenuation function, $Q(\Delta, h)$:
%pylab notebook
delta, depth, mb_q = forensic.q_for_mb()
plt.pcolormesh(delta, depth, mb_q.transpose(), cmap=plt.get_cmap('bone_r'))
plt.xlabel('Distance $(^{\circ})$')
plt.ylabel('Depth (km)')
plt.colorbar()
plt.show()
Analyzing data from MAKZ¶
We'll start by analyzing data from MAKZ, a seismic station in Kazakstan that is part of the Global Seismographic Network (GSN) network. MAKZ is in the set of distance ranges in which it is suitable to measure a mb magnitude (33.8 degrees).
forensic.get_dist?
tr = st.select(station='MAKZ')[0].copy()
forensic.get_dist(tr, evla, evlo)/111.1
What value of $Q(\Delta, h)$ should we use in estimating an mb for MAKZ? We know the event is near the surface, so we can index into our matrix of mb_q values for the closest tabulated distance and a source at 0 km depth.
print(delta.shape)
print(mb_q.shape)
ix = np.argmin(np.abs(delta - forensic.get_dist(tr, evla, evlo)/111.1))
mb_q[ix,0]
Plotting the raw data¶
The amplitudes are in 'counts', which have no physical meaning but are simply related to the data logging system.
%pylab notebook
t = np.arange(0, tr.stats.delta*tr.stats.npts, tr.stats.delta)
plt.plot(t, tr.data, '-')
plt.xlabel('Time (s) after ' + str(tr.stats.starttime).split('.')[0])
plt.ylabel('Counts')
The metadata on the trace contains the response information, however. It tells us all the information we need to convert from counts to units of meters/second (because the sensor is measuring ground velocity).
print(tr.stats)
We're going to have to do some manipulation to the data in order to get it into units of displacement, and then simulate the response of the WWSSN SP seismometer, which is what the mb magnitude scale is calibrated for. Fortunately, ObsPy (using mathematical methods called deconvolution and convolution, which are covered in Digital Signal Processing) comes to the rescue!
We first remove the response of the instrument, and convert to displacement (this performs a numerical integration of the time series). We then filter the data using a highpass and lowpass filter in order to simulate the response of the WWSSN SP seismometer (what it would measure). The filtered seismogram has lost a lot of the low-frequency energy and accentuates the P waves.
We then convert the displacements from meters to nanometers and plot the seismogram. Now, we can measure the amplitude and period and apply the mb formula above.
%pylab notebook
# Converting the seismogram to physical displacement:
tr_disp = tr.copy()
tr_disp.remove_response(output='DISP')
# Simulating the response of a WWSSN short-period instrument with a high and low pass filter:
tr_disp_filt = tr_disp.copy()
tr_disp_filt.taper(type='cosine', max_percentage=0.001)
tr_disp_filt.filter('highpass', freq=1, corners=3)
tr_disp_filt.filter('lowpass', freq=2.8, corners=2)
# Plotting the data (in nanometers):
plt.plot(t, 1e9*tr_disp_filt.data, '-')
plt.ylabel('Displacement (nanometers)')
plt.xlabel('Time (s) after ' + str(tr.stats.starttime).split('.')[0])
Exercise¶
Use the recording above to calculate a mb for this station. Next, look at seismograms from a few other stations at suitable distance ranges and calculate magnitudes to estimate a mean and standard deviation of the magnitudes.
Analyzing data from MDJ¶
The New China Digital Seismograph network station MDJ is one of the closest open stations to the DPRK test site, and has been used extensively in analysis of those events.
tr = st.select(station='MDJ')[1].copy()
g = Geod(ellps='sphere')
_,_,dist = g.inv(evlo,evla,tr.stats.sac.stlo,tr.stats.sac.stla); dist = dist/1000
print(dist/111.1)
%pylab notebook
t = forensic.get_time_for_trace(tr)
plt.plot_date(t, tr.data, '-')
plt.xlim(date2num(evt0.datetime)+30/86400, date2num(evt0.datetime)+200/86400)
Exercise¶
Reproduce Figure 1 in Myers et al. (2018).
Hint: You can request data from only MDJ as follows:
st = forensic.getStreamFromIRIS(evt0-600, evt0+3000, 'MDJ', 'IC', 'BHZ')