Example on how to use zizou to detect anomalous signals¶
We will request data from station MAVZ at Mt. Ruapehu which recorded a volcano-seismic signal on 7 March around 10 am UTC. We will train an autoencoder on data from the preceding three days and then check whether the autoencoder loss will increase during the time of the volcano-seismic signal.
In [1]:
Copied!
from obspy import UTCDateTime
from zizou import SSAM, RSAM, DSAR, SpectralFeatures, FDSNWaveforms, DataSource, AutoEncoder,EnergyExplainedByRSAM
from zizou.visualise import multi_feature_plot_mpl
from tonik import Storage
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from obspy import UTCDateTime
from zizou import SSAM, RSAM, DSAR, SpectralFeatures, FDSNWaveforms, DataSource, AutoEncoder,EnergyExplainedByRSAM
from zizou.visualise import multi_feature_plot_mpl
from tonik import Storage
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
In [2]:
Copied!
fw = FDSNWaveforms(url="GeoNet")
ds = DataSource(clients=[fw])
starttime_test = UTCDateTime("2024-03-07T00:00:00")
endtime_test = starttime_test + 86400
starttime_training = UTCDateTime("2024-03-04T00:00:00")
endtime_training = UTCDateTime("2024-03-06T23:59:59")
rawdata = next(ds.get_waveforms("NZ", "MAVZ", "10", "HHZ", starttime_training, endtime_test, cache=True))
fw = FDSNWaveforms(url="GeoNet")
ds = DataSource(clients=[fw])
starttime_test = UTCDateTime("2024-03-07T00:00:00")
endtime_test = starttime_test + 86400
starttime_training = UTCDateTime("2024-03-04T00:00:00")
endtime_training = UTCDateTime("2024-03-06T23:59:59")
rawdata = next(ds.get_waveforms("NZ", "MAVZ", "10", "HHZ", starttime_training, endtime_test, cache=True))
Compute Features¶
Next we will compute a few typical seismic features and store them using the tonik package:
In [3]:
Copied!
r = RSAM(interval=600)
spec = SSAM(interval=60, resample_int=(None, '10min'))
sf = SpectralFeatures(interval=600)
re = EnergyExplainedByRSAM(interval=600)
dsar = DSAR(interval=600)
s = Storage('anomaly_detection_test', rootdir='/tmp',
starttime=starttime_training.datetime,
endtime=endtime_training.datetime)
st = s.get_substore('NZ', 'MAVZ', '10', 'HHZ')
for feature in [r, spec, sf, re, dsar]:
xds = feature.compute(rawdata)
st.save(xds, mode='w')
r = RSAM(interval=600)
spec = SSAM(interval=60, resample_int=(None, '10min'))
sf = SpectralFeatures(interval=600)
re = EnergyExplainedByRSAM(interval=600)
dsar = DSAR(interval=600)
s = Storage('anomaly_detection_test', rootdir='/tmp',
starttime=starttime_training.datetime,
endtime=endtime_training.datetime)
st = s.get_substore('NZ', 'MAVZ', '10', 'HHZ')
for feature in [r, spec, sf, re, dsar]:
xds = feature.compute(rawdata)
st.save(xds, mode='w')
2024-10-24 22:01:44:zizou.rsam:3982179:96 INFO Computing RSAM for NZ.MAVZ.10.HHZ between 2024-03-04T00:00:00.003130 and 2024-03-07T23:59:59.993130. 2024-10-24 22:01:51:zizou.ssam:3982179:288 INFO Computing spectrograms for NZ.MAVZ.10.HHZ between 2024-03-04T00:00:00.003130 and 2024-03-07T23:59:59.993130. 2024-10-24 22:02:01:zizou.spectral_features:3982179:64 INFO Computing spectral features for NZ.MAVZ.10.HHZ between 2024-03-04T00:00:00.003130 and 2024-03-07T23:59:59.993130. 2024-10-24 22:02:07:zizou.rsam:3982179:208 INFO Computing Energy explained by RSAM for NZ.MAVZ.10.HHZ between 2024-03-04T00:00:00.003130 and 2024-03-07T23:59:59.993130. 2024-10-24 22:02:14:zizou.dsar:3982179:91 INFO Computing DSAR for NZ.MAVZ.10.HHZ between 2024-03-04T00:00:00.003130 and 2024-03-07T23:59:59.993130.
Once features have been computed we can request them using tonik's Storage class:
In [4]:
Copied!
sono = st('sonogram')
rsam = st('rsam')
_ = multi_feature_plot_mpl(sono, rsam)
sono = st('sonogram')
rsam = st('rsam')
_ = multi_feature_plot_mpl(sono, rsam)
Now we can run the autoencoder for anomaly detection. Note, that it depends on the scikit-learn and pytorch packages. These are not intalled with the package so you will have to install them manually.
In [5]:
Copied!
ac = AutoEncoder(st, epochs=100)
ac.fit()
ac = AutoEncoder(st, epochs=100)
ac.fit()
2024-10-24 22:02:36:zizou.autoencoder:3982179:209 INFO Using cpu device
2024-10-24 22:02:36:zizou.anomaly_base:3982179:123 WARNING Number of samples is not consistent.Expected 1271952, got 1271953for feature: sonogram
2024-10-24 22:02:36:zizou.anomaly_base:3982179:123 WARNING Number of samples is not consistent.Expected 1271952, got 1271953for feature: sonogram
2024-10-24 22:02:36:zizou.autoencoder:3982179:217 INFO TiedAutoEncoder(
(weights): ParameterList(
(0): Parameter containing: [torch.float64 of size 2000x14]
(1): Parameter containing: [torch.float64 of size 500x2000]
(2): Parameter containing: [torch.float64 of size 200x500]
(3): Parameter containing: [torch.float64 of size 6x200]
)
(enc_biases): ParameterList(
(0): Parameter containing: [torch.float64 of size 2000]
(1): Parameter containing: [torch.float64 of size 500]
(2): Parameter containing: [torch.float64 of size 200]
(3): Parameter containing: [torch.float64 of size 6]
)
(dec_biases): ParameterList(
(0): Parameter containing: [torch.float64 of size 200]
(1): Parameter containing: [torch.float64 of size 500]
(2): Parameter containing: [torch.float64 of size 2000]
(3): Parameter containing: [torch.float64 of size 14]
)
(loss): MSELoss()
)
2024-10-24 22:02:36:root:3982179:270 INFO Loaded pretrained model.
/workspaces/zizou/src/zizou/autoencoder.py:268: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature. torch.load(os.path.join(store.path, self.files["modelfile"]))
Once we have trained the autoencoder on 'normal' data we will apply it to the test data containing the volcano-seismic signal.
In [9]:
Copied!
ac.transform(starttime_test.datetime, endtime_test.datetime)
ac.transform(starttime_test.datetime, endtime_test.datetime)
2024-10-24 22:03:44:root:3982179:270 INFO Loaded pretrained model. 2024-10-24 22:03:44:zizou.autoencoder:3982179:307 INFO Starting transform process.
/workspaces/zizou/src/zizou/autoencoder.py:268: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature. torch.load(os.path.join(store.path, self.files["modelfile"]))
2024-10-24 22:03:44:zizou.anomaly_base:3982179:123 WARNING Number of samples is not consistent.Expected 1271952, got 1271953for feature: sonogram 2024-10-24 22:03:44:zizou.anomaly_base:3982179:123 WARNING Number of samples is not consistent.Expected 1271952, got 1271953for feature: sonogram 2024-10-24 22:03:44:zizou.autoencoder:3982179:313 INFO Transforming data. 2024-10-24 22:03:45:zizou.autoencoder:3982179:324 INFO Calculating loss between input and decoded output. 2024-10-24 22:03:45:zizou.autoencoder:3982179:329 INFO Clustering embeddings into 5 cluster. 2024-10-24 22:03:45:root:3982179:117 INFO Loaded cluster model. 2024-10-24 22:03:45:root:3982179:96 WARNING Could not set all meta info for autoencoder_loss: 'station' 2024-10-24 22:03:45:root:3982179:96 WARNING Could not set all meta info for autoencoder_cluster: 'station' 2024-10-24 22:03:45:root:3982179:96 WARNING Could not set all meta info for autoencoder_embedding: 'station'
In [10]:
Copied!
sono = st('sonogram')
ac = st('autoencoder_loss')
_ = multi_feature_plot_mpl(sono, ac)
sono = st('sonogram')
ac = st('autoencoder_loss')
_ = multi_feature_plot_mpl(sono, ac)
In [15]:
Copied!
rawdata.trim(UTCDateTime(2024,3,7,9,0), UTCDateTime(2024,3,7,11,0,0))
rawdata.detrend('demean')
rawdata.taper(max_percentage=0.05)
rawdata.filter('bandpass', freqmin=1, freqmax=6)
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(rawdata.times(type='matplotlib'), rawdata.data, color='deepskyblue')
ax1 = ax.twinx()
acloss = ac.loc[{'datetime': slice('2024-03-07T9:00:00', '2024-03-07T11:00:00')}]
ax1.plot(mdates.date2num(acloss.datetime), acloss.data, color='black')
ax.tick_params(axis='y', labelcolor='deepskyblue')
date_format = mdates.DateFormatter('%H:%M')
ax.xaxis.set_major_formatter(date_format)
ax.set_xlabel('Time')
ax.set_ylabel('Seismic amplitude', color='deepskyblue')
ax1.set_ylabel('Autoencoder loss', color='black')
rawdata.trim(UTCDateTime(2024,3,7,9,0), UTCDateTime(2024,3,7,11,0,0))
rawdata.detrend('demean')
rawdata.taper(max_percentage=0.05)
rawdata.filter('bandpass', freqmin=1, freqmax=6)
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(rawdata.times(type='matplotlib'), rawdata.data, color='deepskyblue')
ax1 = ax.twinx()
acloss = ac.loc[{'datetime': slice('2024-03-07T9:00:00', '2024-03-07T11:00:00')}]
ax1.plot(mdates.date2num(acloss.datetime), acloss.data, color='black')
ax.tick_params(axis='y', labelcolor='deepskyblue')
date_format = mdates.DateFormatter('%H:%M')
ax.xaxis.set_major_formatter(date_format)
ax.set_xlabel('Time')
ax.set_ylabel('Seismic amplitude', color='deepskyblue')
ax1.set_ylabel('Autoencoder loss', color='black')
Out[15]:
Text(0, 0.5, 'Autoencoder loss')
