Introduction¶
En géosciences, une grande partie des données que nous manipulons — ondes sismiques, marées océaniques, variations de pression atmosphérique — se présentent en séries temporelles. Ces signaux contiennent souvent de l’information à plusieurs fréquences, et une simple analyse dans le domaine temporel peut masquer des caractéristiques importantes des données.
Les spectrogrammes sont des outils utiles pour analyser et visualiser le contenu fréquentiel des séries temporelles. Un spectrogramme représente la densité spectrale d’un signal au fil de son évolution dans le temps. Il est particulièrement utile dans les applications géoscientifiques où les signaux sont non stationnaires, c’est-à-dire où leur contenu fréquentiel change au cours du temps. Exemples :
- Sismologie des tremblements de terre : les signaux de basse et de haute fréquence y sont pertinents à différents stades d’un événement.
- Sciences de l’atmosphère : les motifs périodiques comme les marées et les ondes y ont des composantes fréquentielles distinctes.
- Télédétection : différents processus, comme les fluctuations d’humidité du sol, peuvent y présenter des fréquences caractéristiques au cours du temps.
Pourquoi les spectrogrammes comptent
- Analyse temps-fréquence : les spectrogrammes montrent comment le contenu fréquentiel d’un signal évolue dans le temps, ce qui les rend adaptés à l’étude des données non stationnaires.
- Extraction de caractéristiques : les spectrogrammes mettent en évidence les événements transitoires et les motifs de longue durée, utiles pour détecter et caractériser des phénomènes géophysiques comme les séismes, les glissements de terrain ou les ondes atmosphériques.
- Analyse multi-échelle : de nombreux processus naturels opèrent à différentes échelles de temps. Les spectrogrammes permettent de visualiser et d’extraire des caractéristiques à travers ces échelles.
- Visualisation : ils offrent une visualisation compacte et interprétable de séries temporelles complexes.
Dans cette section, nous transformons les données en les projetant sur une base de fonctions. Les deux transformées les plus utilisées sont la transformée de Fourier et la transformée en ondelettes.
Avertissement. Filtrer des données, quelles qu’elles soient, demande du soin. Les artefacts de filtrage peuvent conduire à une interprétation complètement fausse. Pièges courants :
- Tous les filtres ne préservent pas la causalité : certains signaux peuvent apparaître avant les événements et être pris à tort pour des précurseurs.
- Filtrer par-dessus des lacunes de données introduit des artefacts à haute fréquence.
- Les effets de bord du filtrage des séries temporelles sont difficiles à atténuer.
La leçon couvre plusieurs niveaux et méthodes de transformation des données.
- Transformées de Fourier : 1D [niveau 1]
- Transformées de Fourier : 2D [niveau 3]
- Spectrogrammes [niveau 2]
- Transformées en ondelettes [niveau 3]
# Import modules for seismic data and feature extraction
import os
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.signal as signal
# seismic python toolbox
import obspy
import obspy.clients.fdsn.client as fdsn
from obspy import UTCDateTime
os.makedirs('data', exist_ok=True)Nous téléchargeons d’abord des données : des sismogrammes enregistrés dans le Puget Sound pour le séisme de Chignik (Alaska) de magnitude 8,2 du 29 juillet 2021. Nous interrogeons le centre de données FDSN avec le client IRIS. Notez que les services de données IRIS sont désormais opérés par EarthScope, mais le nom de client IRIS fonctionne toujours.
Le code de canal HHZ désigne un sismomètre large bande à gain élevé (H), échantillonné à 100 échantillons par seconde (H), composante verticale (Z). C’est un canal large bande : il enregistre le mouvement du sol sur une large bande de fréquences, pas à une fréquence unique.
Nous téléchargeons aussi les métadonnées de la station (la réponse instrumentale) et retirons la réponse avec remove_response(output="VEL"). Cela convertit les comptes bruts du numériseur en vitesse du sol en m/s, pour que les amplitudes aient des unités physiques.
# Download seismic data
network = 'UW'
station = 'RATT'
channel = 'HHZ' # broadband high-gain vertical channel, 100 samples per second
Tstart = UTCDateTime(2021, 7, 29, 6, 15)
Tend = Tstart + 7200
fdsn_client = fdsn.Client('IRIS') # client to query the EarthScope (formerly IRIS) DMC server
# download the station metadata, including the instrument response
inv = fdsn_client.get_stations(network=network, station=station, channel=channel,
starttime=Tstart - 7200, endtime=Tend, level='response')
# call to download the specific data: earthquake waveforms
Z = fdsn_client.get_waveforms(network=network, station=station, location='--', channel=channel,
starttime=Tstart, endtime=Tend)
# basic pre-processing: merge if there are gaps, detrend, taper,
# then remove the instrument response to convert digitizer counts to ground velocity (m/s).
Z.merge(); Z.detrend(type='linear'); Z[0].taper(max_percentage=0.05)
Z.remove_response(inventory=inv, output="VEL")
# call to download the specific data: noise waveforms (the two hours before the earthquake)
N = fdsn_client.get_waveforms(network=network, station=station, location='--', channel=channel,
starttime=Tstart - 7200, endtime=Tstart)
N.merge(); N.detrend(type='linear'); N[0].taper(max_percentage=0.05)
N.remove_response(inventory=inv, output="VEL")/Users/marinedenolle/Dropbox/CLASSES/ESS490/curriculum-book/.pixi/envs/default/lib/python3.12/site-packages/obspy/clients/fdsn/client.py:251: ObsPyDeprecationWarning: IRIS is now EarthScope, please consider changing the FDSN client short URL to 'EARTHSCOPE'.
warnings.warn(msg, ObsPyDeprecationWarning)
1 Trace(s) in Stream:
UW.RATT..HHZ | 2021-07-29T04:15:00.000000Z - 2021-07-29T06:14:59.990000Z | 100.0 Hz, 720000 samplesplt.plot(Z[0].data); plt.plot(N[0].data); plt.grid(True)
plt.ylabel('Velocity (m/s)'); plt.xlabel('Sample index')
plt.legend(['Earthquake', 'Noise'])
1. Transformées de Fourier [niveau 1]¶
Nous utilisons le module scipy.fft pour transformer les deux séries temporelles (séisme et bruit). L’ancien module scipy.fftpack est un module hérité et ne doit plus être utilisé.
La transformée de Fourier est une décomposition de la série temporelle sur une base orthonormée de fonctions cosinus et sinus. La transformée de Fourier d’une série temporelle (de même si la variable est l’espace ) est :
est la valeur de Fourier complexe à la fréquence . La transformée de Fourier détermine quelles fréquences dominent la série temporelle.
1.1 Nyquist¶
La transformée de Fourier que nous utilisons dans ce cours prend une série temporelle discrète de nombres réels. Si la série couvre secondes avec échantillons régulièrement espacés, l’intervalle d’échantillonnage est . La plus haute fréquence résoluble dans une série temporelle discrète, appelée fréquence de Nyquist, est limitée par :
En pratique, on ne peut pas contraindre des signaux qui varient plus vite que deux échantillons temporels. Ici s, donc Hz.
1.2 Incertitudes¶
- La transformée de Fourier discrète fournit une approximation de la TF. Plus la série temporelle est courte, moins la TF est exacte. Cela signifie que la TF sur des fenêtres temporelles courtes est moins exacte.
- La TF suppose (et exige) la périodicité de la série : la série temporelle finie/tronquée est censée se répéter dans le temps. Pour l’imposer, nous apodisons (taper) la série temporelle afin que le premier et le dernier point soient égaux (à zéro).
from scipy.fft import fft, ifft, fftfreq, next_fast_len
npts = Z[0].stats.npts
## FFT the signals
# pad up to a fast FFT length to speed up the FFT
Nfft = next_fast_len(int(Z[0].data.shape[0])) # this will be an even number
freqVec = fftfreq(Nfft, d=Z[0].stats.delta)[:Nfft//2]
Z.taper(max_percentage=0.05)
Zhat = fft(Z[0].data, n=Nfft)Consultez la documentation Obspy pour découvrir la fonction taper. Tracez les spectres d’amplitude et de phase.
fig, ax = plt.subplots(2, 1, figsize=(11, 8))
ax[0].plot(freqVec, np.abs(Zhat[:Nfft//2])/Nfft)
ax[0].grid(True)
ax[0].set_xscale('log'); ax[0].set_yscale('log')
ax[0].set_xlabel('Frequency (Hz)'); ax[0].set_ylabel('Amplitude (m/s)')
ax[0].set_title('Amplitude spectrum, earthquake')
ax[1].hist(np.angle(Zhat[:Nfft//2]))
ax[1].grid(True)
ax[1].set_xlabel('Phase (radians)'); ax[1].set_ylabel('Count')
ax[1].set_title('Distribution of the phase spectrum')
Vous noterez ci-dessus que les valeurs de phase sont distribuées aléatoirement entre et π. Nous pouvons le vérifier en affichant la distribution des spectres de phase et d’amplitude.
# your turn. Plot the histogram of the amplitude spectrum
plt.hist(np.log10(np.abs(Zhat[:Nfft//2])/Nfft), 100); plt.grid(True)
plt.xlabel('log10 amplitude'); plt.ylabel('Count')
plt.show()
Nous pouvons aussi analyser les caractéristiques spectrales de la série temporelle de bruit. Ci-dessous :
- calculez la transformée de Fourier
- tracez les spectres de phase et d’amplitude
- tracez la distribution des valeurs de phase et d’amplitude
# compute Fourier transform of the noise time series
npts1 = N[0].stats.npts
## FFT the signals
# pad up to a fast FFT length to speed up the FFT
Nfft1 = next_fast_len(int(N[0].data.shape[0])) # this will be an even number
freqVec1 = fftfreq(Nfft1, d=N[0].stats.delta)[:Nfft1//2]
# taper the data to enforce periodicity
N.taper(max_percentage=0.05)
# Fourier transform
Nhat = fft(N[0].data, n=Nfft1)# plot the phase and amplitude spectra
fig, ax = plt.subplots(2, 1, figsize=(11, 8))
ax[0].plot(freqVec1, np.abs(Nhat[:Nfft1//2])/Nfft1)
ax[0].grid(True)
ax[0].set_xscale('log'); ax[0].set_yscale('log')
ax[0].set_xlabel('Frequency (Hz)'); ax[0].set_ylabel('Amplitude (m/s)')
ax[0].set_title('Amplitude spectrum, noise')
ax[1].hist(np.angle(Nhat[:Nfft1//2]))
ax[1].grid(True)
ax[1].set_xlabel('Phase (radians)'); ax[1].set_ylabel('Count')
# Overlay the spectrum of the earthquake and the spectrum of the noise
fig, ax = plt.subplots(1, 1, figsize=(11, 8))
ax.plot(freqVec, np.abs(Zhat[:Nfft//2])/Nfft)
ax.plot(freqVec1, np.abs(Nhat[:Nfft1//2])/Nfft1)
ax.grid(True)
ax.set_xscale('log'); ax.set_yscale('log')
ax.set_xlabel('Frequency (Hz)'); ax.set_ylabel('Amplitude (m/s)')
ax.legend(['Earthquake', 'Noise'])
Superposez leurs densités de probabilité (PDF).
# your turn. Plot the histograms of the amplitude spectra
plt.hist(np.log10(np.abs(Zhat[:Nfft//2])/Nfft), 100)
plt.hist(np.log10(np.abs(Nhat[:Nfft1//2])/Nfft1), 100)
plt.grid(True)
plt.xlabel('log10 amplitude'); plt.ylabel('Count')
plt.legend(['Earthquake', 'Noise'])
plt.show()
Vous remarquez que leurs différences statistiques se situent dans les queues des distributions. Des métriques statistiques comme la mean ou la variance peuvent donc ne pas être discriminantes, mais le kurtosis pourrait l’être.
# print short float values
print(f"Skewness of earthquake {scipy.stats.skew(np.log10(np.abs(Zhat[:Nfft//2])))} and noise {scipy.stats.skew(np.log10(np.abs(Nhat[:Nfft1//2])))}")
print(f"Kurtosis of earthquake {scipy.stats.kurtosis(np.log10(np.abs(Zhat[:Nfft//2])))} and noise {scipy.stats.kurtosis(np.log10(np.abs(Nhat[:Nfft1//2])))}")
print(f"Mean of earthquake {np.mean(np.log10(np.abs(Zhat[:Nfft//2])))} and noise {np.mean(np.log10(np.abs(Nhat[:Nfft1//2])))}")
print(f"standard deviation of earthquake {np.std(np.log10(np.abs(Zhat[:Nfft//2])))} and noise {np.std(np.log10(np.abs(Nhat[:Nfft1//2])))}")Skewness of earthquake 2.4536951625773393 and noise 0.7840246929880705
Kurtosis of earthquake 11.193100408198315 and noise 2.7513010817145336
Mean of earthquake -4.772491127799235 and noise -4.872561055684669
standard deviation of earthquake 0.7077083875872446 and noise 0.4345104292391202
2. Transformées de Fourier 2D [niveau 3]¶
La transformée de Fourier 2D s’applique à une matrice 2D. Elle applique d’abord une transformée de Fourier 1D à chaque ligne de la matrice, puis une transformée de Fourier 1D à chaque colonne de la matrice intermédiaire.
Les transformées de Fourier 2D donnent les coefficients de Fourier qui dominent une image. On peut s’en servir pour filtrer les données. Une autre application est la compression des données : garder quelques coefficients au lieu de stocker l’image entière.
Nous nous exercerons sur un champ synthétique de type topographique. La vraie topographie a un spectre de nombres d’onde « rouge » : l’essentiel de la puissance se trouve aux grandes longueurs d’onde (les chaînes de montagnes), avec de moins en moins de puissance aux courtes longueurs d’onde (la rugosité de petite échelle). L’amplitude spectrale décroît à peu près comme une loi puissance du nombre d’onde , . Nous pouvons construire un tel terrain fractal directement dans le domaine de Fourier : tirer des phases aléatoires, mettre les amplitudes à l’échelle en et appliquer la transformée inverse.
from scipy.fft import fft2, ifft2, fftshift
from scipy.fft import fftfreq as fftfreq2d
# build a synthetic fractal terrain with a power-law (red) wavenumber spectrum
rng = np.random.default_rng(42)
nx, ny = 512, 512 # grid size
dx = 1.0 # grid spacing in km
beta = 2.0 # spectral decay exponent: amplitude ~ k^-beta
# wavenumber grids (cycles per km)
kx = fftfreq2d(nx, d=dx)
ky = fftfreq2d(ny, d=dx)
KX, KY = np.meshgrid(kx, ky, indexing='ij')
K = np.sqrt(KX**2 + KY**2)
K[0, 0] = np.inf # avoid division by zero at the zero wavenumber (mean)
# random phases on a white-noise field, amplitudes scaled by k^-beta
white = rng.standard_normal((nx, ny))
shaped = fft2(white) * K**(-beta)
elevation = np.real(ifft2(shaped))
# rescale to a plausible elevation range in meters
elevation = 2000 * (elevation - elevation.min()) / (elevation.max() - elevation.min())
# spatial coordinate vectors in km
xkm = np.arange(nx) * dx
ykm = np.arange(ny) * dx
plt.figure(figsize=(7, 6))
plt.contourf(xkm, ykm, elevation.T, 30, cmap='terrain')
plt.colorbar(label='Elevation (m)')
plt.xlabel('x (km)'); plt.ylabel('y (km)')
plt.title('Synthetic fractal terrain')
plt.axis('scaled')
plt.show()
Considérez elevation comme un jeu de données 2D. Nous pouvons effectuer une transformée 2D, qui donne un spectre dans les dimensions spatiales. Les axes de l’image transformée sont des nombres d’onde (cycles par km). Nous utilisons fftshift pour placer le nombre d’onde zéro au centre de l’image.
Zel = fft2(elevation)
kmax = 1 / (2 * dx) # Nyquist wavenumber in cycles per km
plt.figure(figsize=(7, 6))
plt.imshow(fftshift(np.log10(np.abs(Zel) / Zel.size)), cmap='RdYlBu',
extent=[-kmax, kmax, -kmax, kmax])
plt.colorbar(label='log10 amplitude')
plt.xlabel('$k_x$ (km$^{-1}$)'); plt.ylabel('$k_y$ (km$^{-1}$)')
plt.title('2D FT of elevation')
plt.show()
L’énergie est concentrée près du centre du graphique (bas nombres d’onde, grandes longueurs d’onde), comme voulu par construction. Nous allons maintenant compresser l’image en ne gardant que les plus grands coefficients de Fourier.
# Sort the Fourier coefficient amplitudes
Zsort = np.sort(np.abs(Zel).reshape(-1))
print(len(Zsort))
print(Zsort.shape)262144
(262144,)
fig, ax = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
for i, keep in enumerate((0.1, 0.05, 0.01)):
thresh = Zsort[int(np.floor((1 - keep) * len(Zsort)))]
ind = np.abs(Zel) > thresh
Atlow = Zel * ind # zero out the small coefficients
Alow = np.real(ifft2(Atlow))
ax[i].contourf(xkm, ykm, Alow.T, 30, cmap='terrain')
ax[i].set_title(f'keep {keep*100:.0f}% of coefficients')
ax[i].axis('scaled')
ax[i].set_xlabel('x (km)')
ax[0].set_ylabel('y (km)')
plt.show()
Comparons maintenant le jeu de données 2D original avec les données compressées par Fourier. Garder 1 % des coefficients correspond à un taux de compression de 100:1, et la reconstruction capture encore la structure de grande échelle du terrain. Le prix à payer est la perte de la rugosité de petite échelle, qui vit dans les coefficients à haut nombre d’onde écartés.
keep = 0.01
thresh = Zsort[int(np.floor((1 - keep) * len(Zsort)))]
ind = np.abs(Zel) > thresh
Atlow = Zel * ind # zero out the small coefficients
print("We are keeping %.2f%% of the Fourier coefficients, a compression ratio of %d:1" % (keep*100, int(1/keep)))
Alow = np.real(ifft2(Atlow))
fig, ax = plt.subplots(1, 2, figsize=(10, 5), sharex=True, sharey=True)
ax[0].contourf(xkm, ykm, elevation.T, 30, cmap='terrain'); ax[0].set_title('Original data')
ax[0].axis('scaled')
ax[1].contourf(xkm, ykm, Alow.T, 30, cmap='terrain'); ax[1].set_title('Compressed data (1% of coefficients)')
ax[1].axis('scaled')
# quantify the reconstruction error
rel_err = np.linalg.norm(elevation - Alow) / np.linalg.norm(elevation)
print(f"Relative reconstruction error: {rel_err:.3f}")We are keeping 1.00% of the Fourier coefficients, a compression ratio of 100:1
Relative reconstruction error: 0.006

3. Spectrogrammes [niveau 2]¶
Dans les problèmes dépendants du temps et multi-échelles, il peut être intéressant d’extraire des caractéristiques des données à partir de la transformée de Fourier à court terme (short-time Fourier transform, STFT).
La STFT est une transformée de Fourier appliquée à des fenêtres courtes (avec recouvrement) pour résoudre les fréquences à différents instants de la série.
from scipy.signal import stft
fs = Z[0].stats.sampling_rate
nperseg = 1000
z = np.asarray(Z[0].data)
f, t, Zxx = stft(z, fs=fs, nperseg=nperseg, noverlap=200)
fig, ax = plt.subplots(2, 1, figsize=(11, 8), sharex=True)
logZ = np.log10(np.abs(Zxx) + 1e-20)
ax[0].pcolormesh(t/3600, f, logZ, vmin=np.percentile(logZ, 50), vmax=np.percentile(logZ, 99.9), shading='gouraud')
ax[0].set_title('STFT Magnitude, earthquake')
ax[0].set_ylabel('Frequency [Hz]')
ax[0].set_yscale('log'); ax[0].set_ylim(0.1, 40)
n = np.asarray(N[0].data)
fn, tn, Nxx = stft(n, fs=fs, nperseg=nperseg, noverlap=200)
logN = np.log10(np.abs(Nxx) + 1e-20)
ax[1].pcolormesh(tn/3600, fn, logN, vmin=np.percentile(logN, 50), vmax=np.percentile(logN, 99.9), shading='gouraud')
ax[1].set_title('STFT Magnitude, noise')
ax[1].set_ylabel('Frequency [Hz]')
ax[1].set_xlabel('Time [Hours]'); ax[1].set_yscale('log'); ax[1].set_ylim(0.1, 40)(0.1, 40)
Le spectrogramme Zxx est une transformée des données originales. Il est courant d’utiliser les spectrogrammes comme entrées de réseaux de neurones, sous forme de tableaux 2D.
Zxx.shape(501, 902)4. Transformée en ondelettes continue [niveau 3]¶
Introduction¶
La transformée en ondelettes continue (CWT, Continuous Wavelet Transform) est un outil important de l’analyse des données géoscientifiques, en particulier pour l’analyse temps-fréquence des signaux non stationnaires. Comme le spectrogramme, la CWT éclaire la façon dont le contenu fréquentiel d’un signal varie au cours du temps. Mais la CWT offre une meilleure résolution aux différentes fréquences, ce qui la rend plus adaptée à l’analyse des signaux à changements fréquentiels transitoires ou localisés.
La transformée en ondelettes décompose un signal en versions dilatées et translatées d’une petite fonction oscillante appelée ondelette. Cela rend la CWT bien adaptée aux applications géoscientifiques, où de nombreux phénomènes, comme les séismes, les éruptions volcaniques et les régimes météorologiques, peuvent se manifester à différentes échelles et fréquences.

Figure : fonctions de base de Fourier et en ondelettes. Image tirée de cet article.
Il existe de nombreuses familles d’ondelettes canoniques. Les familles diffèrent typiquement par leur forme, leur compacité et leur régularité. On choisit typiquement une famille pour une série temporelle donnée. Les ondelettes ont une énergie finie et une moyenne nulle.

Figure : familles de fonctions de base en ondelettes. Image tirée de cet article
La transformée en ondelettes est :
où est l’ondelette mère dilatée d’un facteur et translatée/décalée de . Dans la transformée continue, et prennent des valeurs continues. La transformée en ondelettes discrète est la transformée en ondelettes effectuée sur un nombre fini d’échelles et de décalages.
La représentation temps-échelle d’une série temporelle est un scalogramme. Les échelles peuvent être converties en pseudo-fréquences : si est la fréquence centrale de l’ondelette, l’échelle et l’intervalle d’échantillonnage, alors la pseudo-fréquence est .
Pourquoi les transformées en ondelettes continues sont importantes¶
- Analyse multirésolution : la CWT capture à la fois les tendances de basse fréquence et de longue durée et les caractéristiques de haute fréquence et de courte durée. C’est particulièrement précieux en géosciences, où les processus se déroulent à différentes échelles temporelles et spatiales.
- Données non stationnaires : de nombreux signaux géoscientifiques sont non stationnaires, c’est-à-dire que leurs propriétés statistiques changent au cours du temps. La CWT révèle ces composantes fréquentielles variables dans le temps.
- Détection de caractéristiques locales : la CWT convient bien à l’identification d’événements localisés, comme les ondes sismiques, les glissements de terrain ou les perturbations atmosphériques, en analysant l’évolution du signal à la fois en temps et en fréquence.
- Meilleure résolution temps-fréquence : la transformée en ondelettes localise le temps et la fréquence plus précisément que des méthodes comme la transformée de Fourier ou le spectrogramme, surtout pour les événements de courte durée.
Usage en sismologie¶
Dans les données sismiques, les différents types d’ondes sismiques (ondes P, ondes S et ondes de surface) se produisent à différentes fréquences et durées. En appliquant la CWT, les sismologues peuvent détecter ces différentes phases d’ondes et leurs temps d’arrivée précis, ce qui compte pour la caractérisation des séismes.
Exemple Python : transformée en ondelettes continue de données sismiques¶
Nous utilisons le paquet PyWavelets (pywt). Les anciennes fonctions scipy.signal.cwt et scipy.signal.morlet2 ont été retirées de SciPy, donc pywt est désormais l’outil standard. Nous choisissons l’ondelette de Morlet complexe 'cmor1.5-1.0' (largeur de bande 1,5, fréquence centrale 1,0), un choix courant pour les données sismiques parce que sa forme complexe donne à la fois l’amplitude et la phase.
pywt.cwt travaille en échelles. Nous choisissons les fréquences que nous voulons résoudre (0,1 à 40 Hz, espacées en logarithme), puis les convertissons en échelles par la relation , ce que calcule pywt.frequency2scale. Pour alléger le calcul, nous appliquons la CWT aux 20 premières minutes de l’enregistrement du séisme, qui contiennent les ondes P, S et de surface.
import pywt
fs = Z[0].stats.sampling_rate
dt = 1 / fs
# analyze the first 20 minutes of the earthquake record
nsel = int(1200 * fs)
zsel = z[:nsel]
tsel = np.arange(nsel) / fs
# choose the frequencies to resolve, then convert them to scales
wavelet = 'cmor1.5-1.0'
freqs_target = np.logspace(np.log10(0.1), np.log10(40), 100)
scales = pywt.frequency2scale(wavelet, freqs_target * dt) # scales for the target frequencies
# compute the CWT; frequencies returned in Hz thanks to sampling_period
coefficients, frequencies = pywt.cwt(zsel, scales, wavelet, sampling_period=dt, method='fft')
print(coefficients.shape, frequencies.min(), frequencies.max())(100, 120000) 0.1 39.99999999999999
# plot the scalogram with a frequency axis in Hz
logC = np.log10(np.abs(coefficients) + 1e-20)
plt.figure(figsize=(11, 5))
plt.pcolormesh(tsel, frequencies, logC, cmap='viridis',
vmin=np.percentile(logC, 50), vmax=np.percentile(logC, 99.9), shading='auto')
plt.yscale('log')
plt.ylim(frequencies.min(), frequencies.max())
plt.colorbar(label='log10 |CWT coefficient|')
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.title('Scalogram of the Chignik earthquake at UW.RATT (cmor1.5-1.0)')
plt.show()
Le scalogramme montre l’onde P qui arrive en premier avec une énergie à haute fréquence, suivie de l’onde S et des ondes de surface longue période sous 0,1–0,5 Hz. Comparez avec le spectrogramme STFT ci-dessus : la transformée en ondelettes résout mieux les ondes de surface basse fréquence parce que sa fenêtre temporelle s’adapte à chaque fréquence.
Avantages de la CWT en géosciences :¶
- Détection d’événements localisés dans le temps : la CWT identifie les événements géophysiques de courte durée, comme l’arrivée des ondes sismiques pendant un séisme.
- Phénomènes multi-échelles : les processus naturels comme l’activité tectonique et les marées océaniques opèrent sur une large gamme d’échelles temporelles et spatiales. La CWT peut analyser les signaux de ces processus en utilisant une large gamme d’échelles.
- Détection de bords : la transformée en ondelettes est particulièrement bonne pour détecter les changements ou les discontinuités dans les signaux, comme le début abrupt des ondes sismiques ou les frontières entre différentes couches géophysiques.
Cas d’usage en géosciences :¶
- Alerte sismique précoce : en appliquant la CWT, les géoscientifiques peuvent détecter plus exactement la première arrivée des ondes sismiques, ce qui compte pour les systèmes d’alerte précoce.
- Détection de glissements de terrain : la CWT aide à identifier les signaux à haute fréquence révélateurs de glissements de terrain, ces événements se caractérisant souvent par de courtes bouffées d’énergie.
- Trémors volcaniques : la CWT peut servir à analyser les signaux de trémor volcanique, souvent complexes, qui présentent à la fois des caractéristiques de courte et de longue durée.
Conclusion¶
La transformée en ondelettes continue (CWT) est un outil utile pour l’analyse temps-fréquence en géosciences. Sa capacité à résoudre les signaux à travers différentes échelles et fréquences la rend bien adaptée à l’étude des processus non stationnaires. Les bibliothèques Python comme pywt, obspy et matplotlib permettent aux géoscientifiques d’appliquer la CWT à leurs données et d’en extraire une compréhension des phénomènes géophysiques complexes.
Les transformées temps-fréquence prennent du temps de calcul dans un flux de travail. Comparons le coût de la CWT et de la STFT sur le même segment de 20 minutes.
import time
tic = time.perf_counter()
pywt.cwt(zsel, scales, wavelet, sampling_period=dt, method='fft')
toc = time.perf_counter()
print(f"CWT (100 scales, 20 min of data): {toc - tic:.2f} s")
tic = time.perf_counter()
stft(zsel, fs=fs, nperseg=nperseg, noverlap=200)
toc = time.perf_counter()
print(f"STFT (same data): {toc - tic:.3f} s")CWT (100 scales, 20 min of data): 0.46 s
STFT (same data): 0.001 s
La STFT est bien moins coûteuse que la CWT. De ces transformées, nous pouvons extraire des caractéristiques statistiques similaires.