Rockfall localization¶
This notebook presents an algorithm for the localization of rockfalls by means of recorded seismic signals. You are very welcomed to play around and test the influence of different parameters on the rockfall location probabilities!
The here presented rockfalls occured at the flanks of Dolomieu crater on Piton de la Fournaise volcano, La Réunion. Seismic signals generated by the rockfalls were recorded at 4 seismic stations surrounding the crater: BON, BOR, DSO, and SNE. Have a look at the following map to see their positions as well as the trajectories of 2 exemplary rockfalls.
The recorded seismic signals generated by these 2 events are stored in the subdirectory data/signals/.
The localization method is based on the comparison between observed and simulated energy ratios of station pairs. Simulations of the seismic wave propagation were carried out using the Spectral Element Method (SEM). For this, a numerical domain was created with the topography of Dolomieu crater and realistic subsurface properties. In order to explore the rockfall locations which best fit the data, a grid of potential source positions was defined. This way, we hold information of the seismic energy generated by each source position for each component (east E, north N, vertical Z) of each seismic station (note that station DSO only contains the vertical component Z). The synthetic energy values are stored in data/simu/.
By analysing the observed seismic signals in a sliding time window, we can gain insight in the spatio-temporal evolution of the rockfall. How cool is this?! We will visualize the most probable rockfall positions as a function of time on a map using a 2D colormap.
This is a very brief introduction of the localization method. Nevertheless, you can start exploring and test for example the sensitivity of the localization on the width of the sliding time window, on the influence of site effects, or using only vertical components vs. using all 3 components.
Have fun!
# Adjust notebook style
from IPython.core.display import HTML
css_file = 'style.css'
HTML(open(css_file,'r').read())
%matplotlib inline
# Import essential libraries
import numpy as np
import os
import time
from obspy import read, Stream
from obspy.core import UTCDateTime
from matplotlib.dates import num2date
# Import functions from modules subdirectory
from modules.m_PLOTS import *
from modules.m_PickTW import cd, pickTW
from modules.m_PROB import calc_prob
from modules.m_SEdeconv import deconvolve_SE
# Reduce figure size
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 60.
Define Parameters¶
You can select the following values for each parameter:
| Parameter | Value | Function |
|---|---|---|
| event | '2016.348', '2017.022' | Choose between 2 rockfalls |
| chID | 'Z', 'ENZ' | Only vertical or 3 components |
| site_eff_deconv | True, False | Site effect correction of signals |
| lwin | float | Length (seconds) of sliding time window |
| swin | float | Step size (seconds) of sliding time window |
| freqband | [fmin, fmax] | Corner values (Hertz) of frequency band which is used for localization |
# You are welcomed to play around with the following parameters
event = '2016.348'
chID = 'ENZ'
site_eff_deconv = True
lwin = 4.
swin = 2.
# As of now you cannot play around with the frequency band
freqband = [13., 17.]
# You better DON'T play around with the grid dimensions
xdim = 121
ydim = 101
# Define closest station according to event
if event == '2016.348':
nearSTA = 'BOR'
elif event == '2017.022':
nearSTA = 'BON'
else:
print("Please select event '2016.348' or '2017.022'")
Load rockfall signals¶
The rockfall seismic signals are loaded using Obspy. All signals are stored as Traces in a Stream.
# Create Stream with all signals
st = Stream()
if chID == 'ENZ':
channels = [ 'Z','N','E']
elif chID == 'Z':
channels = [ 'Z']
path_event = 'data/signals/'+event+'/'
for ch in channels:
st += read(path_event+"/PF.*"+ch+".D.20*.VEL")
# Pre-filter the signal:
st.filter('bandpass', freqmin=1., freqmax=40.,
corners = 2, zerophase = 'true')
# Determine starttime t0
t0 = st[0].stats.starttime
eventDate = (str(t0.year)+'-'+str(t0.month).zfill(2)+'-'+
str(t0.day).zfill(2)+' '+str(t0.hour).zfill(2)+':'+
str(t0.minute).zfill(2)+':'+str(t0.second).zfill(2))
print(st)
Let's have a first look on the rockfall signals!
plot_traces([st],['raw'],eventDate)
Processing the signals¶
The signal at a station can be amplified due to the local soil structure. This can bias the localization when the stations are differently influenced. In order to aboid this bias, we have to remove these site effects from the recorded signals. This can be done by deconvolution in the time domain, which becomes a division in the frequency domain. The spectral site amplification functions were estimated using seismic signals of volcano-tectonic (VT) events.
# Site effect deconvolution
if site_eff_deconv:
st_raw = st.copy()
# We remove the site effects from each trace tr:
for tr in st:
tr = deconvolve_SE(tr)
# Plotting the site corrected signals
plot_traces([st_raw,st],['raw','corrected'],eventDate)
You can see that in particular station SNE was amplified. After the site effect correction we hope that the localization is less biased.
In the following we band-pass filter the signals in the frequency band in which observed and simulated energy ratios will be compared.
# Band-pass filtering the signals
freqmin = freqband[0]
freqmax = freqband[1]
freqID = 'BP'+str(int(freqmin))+'to'+str(int(freqmax))
st.filter('bandpass', freqmin=freqmin, freqmax=freqmax,
corners = 2, zerophase = 'true')
plot_traces([st],['filtered'],freqID+' - '+eventDate)
Pick start and end time of rockfall signal¶
Let's pick the time window of the signal which shall be used to localize the rockfall. The picking is performed interactively on the figure using the mouse curser. If the file with the picked start and end times for the current event already exists, this step will be ignored (just delete or rename the files pickedTW_*event*.txt if you want to do the picking yourself).
if os.path.isfile('pickedTW_'+event+'.txt'):
newPick = False
print('File with picked time window already exists.')
plot_trace(st.select(station=nearSTA,channel='*Z')[0],event)
else:
# Activate interactive plotting:
%matplotlib notebook
%matplotlib notebook
import matplotlib.pyplot as plt
newPick = True
pickTW(st.select(station=nearSTA,channel='*Z')[0])
In the following we load the previously picked start and end time.
# Only in case of newly picked time window
if newPick:
# Rename file with picked time window according to event:
os.rename('pickedTW.txt', 'pickedTW_'+event+'.txt')
# Deactivate interactive plotting:
%matplotlib inline
%matplotlib inline
import matplotlib.pyplot as plt
# Load time window
pickedTW = np.loadtxt('pickedTW_'+event+'.txt')
tstart = UTCDateTime(num2date(pickedTW[0]))
tend = UTCDateTime(num2date(pickedTW[1]))
# Define start of time window relative to signal starttime
offset = (tstart+swin-lwin/2) - t0
print('Starttime: ', tstart)
print('Endtime: ', tend)
Location optimization routine¶
Inter-station energy ratios simulated from each grid position are now compared to the observed energy ratios to find the most probable source locations. This way, the rockfall signal is analyzed over its whole duration using a sliding time window. Finally, at each grid position, the highest probability is stored together with the corresponding time. The result will be visualized subsequently.
# Pre-allocate arrays of probabilities and corresponding times
probs = np.zeros((ydim,xdim))
times = np.zeros((ydim,xdim))
# Mark start time
t1 = time.time()
print('Sliding time window (TW):')
# Analyze signal in sliding time window
for iTW, windowed_st in enumerate(st.slide(
window_length=lwin,
step=swin,
offset=offset,
include_partial_windows=True)):
# Condition to stop analysis at the end of the time window
if (windowed_st[0].stats.starttime+lwin/2) > tend:
break
print('TW ', iTW)
# Algorithm which compares observed energy ratios with
# synthetic energy ratios from each potential source position
probs_tw = calc_prob(windowed_st, channels, ydim, xdim, freqID)
# Find highest probability across the whole signal duration
# and store probability value and corresponding time
if iTW == 0:
t_rel = 0.
times[:,:] = t_rel
probs = probs_tw
else:
t_rel += swin
for ii in range(ydim):
for jj in range(xdim):
if probs[ii,jj] < probs_tw[ii,jj]:
times[ii,jj] = t_rel
probs[ii,jj] = probs_tw[ii,jj]
# Mark end time and print elapsed computation time
t2 = time.time()
print('Elapsed time: ', (t2-t1),'s')
Result Visualization¶
In the previous step we analyzed the rockfall signal and stored for each potential source position the highest probability and the corresponding time. In the following, these 2 values are visualized on a map using a 2D colormap which represents the probability in lightness and the time in color.
# Define array of size (2, nx, ny) which contains maps of
# times and probabilities.
loc_prob_time = np.array([ np.flipud(times),
np.flipud(probs)/10. ])
# Plot the spatio-temporal rockfall location probability.
plot_LOCA(loc_prob_time, t_rel, swin)
# Plot colorbar and trace to visualize the time window used for
# the localization.
plot_Tr_Cb(st.select(station=nearSTA,channel='*Z')[0],
pickedTW[0], pickedTW[1], t_rel, swin)
Conclusions¶
Congrats! You just located a rockfall :)
You can observe that the most probable rockfall locations are close to the crater rim in the the beginning of the signal. Towards later times, the rockfall moves towards the crater bottom.
Comparing the trajectory to the one drawn on the map in the beginning of this notebook, we can see that the rockfall is localized in the correct region. The uncertainty of the rockfall location is in the order of 100m and increases with increasing time. This is because the rockfall seismic source becomes more spatially distributed over time. The method however assumes a single source position. This is a limitation of the method.
The presented localizization metheod is still in developement. For example, in order to refine the localization we will exploit multiple frequency bands. Still, the method allows a first estimation of the rockfall location. The use of a sliding time window for the analysis allows monitoring of rockfall activity in real-time.