2.10 Bruit synthétique - Machine Learning en géosciences

Génération de bruit synthétique et de données synthétiques pour les géosciences

1. Introduction¶

Le bruit des données géoscientifiques vient de l’environnement, de l’instrument et de l’activité humaine. Si nous savons simuler le bruit, nous pouvons construire des jeux de données contrôlés à vérité terrain connue. Cela permet de tester des filtres, d’évaluer des détecteurs par rapport à une référence et d’augmenter les ensembles d’entraînement pour l’apprentissage automatique.

Cette leçon traite du bruit synthétique dans deux contextes :

Séries temporelles 1D : enregistrements sismiques, pression atmosphérique, température de surface de la mer.
Champs géospatiaux 2D : topographie, humidité du sol, cartes de température.

Nous terminons par un exemple sismologique réaliste et par une discussion sur les cas où les données synthétiques sont recevables en recherche.

Types de bruit¶

Bruit blanc : densité spectrale de puissance constante à toutes les fréquences. Typique du bruit instrumental.
Bruit rose : densité spectrale de puissance proportionnelle à $1/f$ . Courant dans les systèmes naturels comme les enregistrements climatiques.
Bruit gaussien : bruit dont la distribution d’amplitude est normale. Le bruit blanc est souvent tiré d’une gaussienne.
Bruit spatialement corrélé : bruit qui varie de façon régulière dans l’espace, comme les artefacts atmosphériques ou de capteur dans les cartes.

Dans tout ce carnet (notebook), nous utilisons le générateur aléatoire moderne de numpy, initialisé par une graine pour la reproductibilité. La graine compte : un jeu de données synthétique qu’on ne sait pas régénérer n’est pas reproductible.

🖥️ Diapositives du cours — Séance 09 (lun. 19 oct.)

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)

2. Modèles de bruit en 1D et en 2D¶

2.1 Bruit blanc et bruit rose (1D)¶

Le bruit blanc est une suite de tirages gaussiens indépendants. Le bruit rose se construit dans le domaine fréquentiel : nous imposons un spectre d’amplitude proportionnel à $1/\sqrt{f}$ (de sorte que la puissance varie en $1/f$ ), nous tirons des phases aléatoires et nous revenons dans le domaine temporel. Nous utilisons np.fft.irfft, qui ne prend le spectre qu’aux fréquences positives ou nulles et renvoie par construction un signal réel.

n = 1000  # number of time points
time = np.linspace(0, 100, n)

# White noise: independent Gaussian draws, zero mean
white_noise = rng.standard_normal(n)


def generate_pink_noise(size, rng):
    """Pink (1/f) noise: power spectral density proportional to 1/f.

    Build the one-sided spectrum with amplitude 1/sqrt(f) and random
    phases, then invert with irfft to get a real time series.
    """
    freqs = np.fft.rfftfreq(size)
    amplitude = np.zeros_like(freqs)
    amplitude[1:] = 1.0 / np.sqrt(freqs[1:])  # skip f=0 to avoid division by zero
    phases = rng.uniform(0, 2 * np.pi, len(freqs))
    spectrum = amplitude * np.exp(1j * phases)
    pink = np.fft.irfft(spectrum, n=size)
    return pink / np.std(pink)


pink_noise = generate_pink_noise(n, rng)

plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(time, white_noise, label='White Noise')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Synthetic White Noise')
plt.grid()

plt.subplot(2, 1, 2)
plt.plot(time, pink_noise, label='Pink Noise', color='orange')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Synthetic Pink Noise (1/f Noise)')
plt.grid()

plt.tight_layout()
plt.show()

Le bruit blanc fluctue rapidement, sans mémoire d’un échantillon au suivant. Le bruit rose vagabonde : les basses fréquences portent l’essentiel de la puissance, donc la série dérive sur de longues échelles de temps.

Applications. Le bruit blanc simule le bruit propre des instruments dans les enregistrements sismiques ou atmosphériques. Le bruit rose imite la variabilité naturelle, dominée par les basses fréquences dans de nombreux processus géophysiques.

2.2 Bruit spatialement corrélé (2D)¶

Le bruit géospatial est rarement indépendant d’un pixel à l’autre. Les perturbations atmosphériques et la dérive des capteurs varient de façon régulière dans l’espace. Une manière simple de le simuler : engendrer un champ 2D de bruit blanc, puis le lisser avec un filtre gaussien. La largeur du filtre fixe la longueur de corrélation.

from scipy.ndimage import gaussian_filter

ngrid = 100  # size of the grid
x = np.linspace(0, 10, ngrid)
y = np.linspace(0, 10, ngrid)
X2d, Y2d = np.meshgrid(x, y)

# 2D white noise
white_noise_2d = rng.standard_normal((ngrid, ngrid))

# Spatially correlated noise: smooth the white noise
spatially_correlated_noise = gaussian_filter(white_noise_2d, sigma=3)

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(white_noise_2d, extent=[0, 10, 0, 10], cmap='viridis')
plt.colorbar(label='Amplitude')
plt.title('2D White Noise')

plt.subplot(1, 2, 2)
plt.imshow(spatially_correlated_noise, extent=[0, 10, 0, 10], cmap='viridis')
plt.colorbar(label='Amplitude')
plt.title('2D Spatially Correlated Noise')

plt.tight_layout()
plt.show()

Applications. Le bruit spatialement corrélé simule la distorsion atmosphérique dans l’imagerie satellitaire, les artefacts d’interpolation dans les produits maillés et les erreurs de capteur régulières dans les cartes géophysiques.

3. Événements transitoires synthétiques¶

Nous construisons maintenant une série temporelle qui contient un événement transitoire plus du bruit. C’est l’ingrédient de base d’un problème de classification binaire : une fenêtre contient-elle du signal ou seulement du bruit ?

3.1 Le signal de l’événement : l’ondelette de Ricker¶

L’ondelette de Ricker, aussi appelée ondelette « chapeau mexicain », est un modèle standard de source sismique impulsive. C’est l’opposé de la dérivée seconde d’une gaussienne. scipy.signal.ricker a été retiré des versions récentes de scipy, nous la définissons donc nous-mêmes.

import scipy.signal as sig


def ricker(t, a):
    """Ricker (Mexican-hat) wavelet.

    Parameters
    ----------
    t : array
        Time axis in seconds, centered on zero.
    a : float
        Width parameter in seconds. The amplitude spectrum peaks near
        sqrt(2) / (2 * pi * a) Hz.

    Returns
    -------
    array
        The wavelet 2/(sqrt(3a) * pi**0.25) * (1 - (t/a)**2) * exp(-t**2/(2a**2)).
    """
    prefactor = 2.0 / (np.sqrt(3.0 * a) * np.pi**0.25)
    return prefactor * (1.0 - (t / a) ** 2) * np.exp(-(t**2) / (2.0 * a**2))

fs = 100.0   # sampling rate in Hz
twin = 50.0  # window length in seconds
t = np.linspace(0, twin, int(twin * fs))

a = 0.5  # width parameter in seconds
t_wavelet = np.arange(-2, 2, 1 / fs)  # 4 seconds of support
sa = ricker(t_wavelet, a)

# place the wavelet in the middle of the time series
pad = len(t) // 2 - len(sa) // 2
s = np.concatenate((np.zeros(pad), sa, np.zeros(len(t) - pad - len(sa))))

plt.figure()
plt.plot(t, s)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.xlim([0, 50])
plt.title('Signal: Ricker wavelet')
plt.show()

L’ondelette de Ricker est régulière et à bande limitée. Traçons son spectre d’amplitude de Fourier en valeur absolue.

from scipy.fft import fft, fftfreq, next_fast_len

# pad to a fast FFT length
Nfft = next_fast_len(len(s))  # this will be an even number

freqVec = fftfreq(Nfft, d=1 / fs)[: Nfft // 2]
Zhat = fft(s, n=Nfft)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].plot(freqVec, np.abs(Zhat[: Nfft // 2]))
ax[0].grid()
ax[0].set_title('FFT of the signal')
ax[0].set_xlabel('Frequency in Hz')
ax[1].plot(freqVec, np.abs(Zhat[: Nfft // 2]))
ax[1].set_xscale('log')
ax[1].set_xlabel('Frequency in Hz')
ax[1].set_title('FFT of the signal in log')
ax[1].grid()

À quoi ressemble la distribution des données de l’événement ?

plt.hist(sa, bins=10)
plt.xlabel('Amplitude')
plt.ylabel('Counts')
plt.title('Histogram of the wavelet')
plt.show()

3.2 Bruit gaussien et rapport signal sur bruit¶

Nous avons créé un signal pur. Créons maintenant une série temporelle de bruit à lui ajouter.

Une précaution : np.random.uniform(0, 1) tire des valeurs entre 0 et 1, dont la moyenne vaut 0,5. Ajouter cela à un signal injecte une composante continue — un pic parasite à la fréquence nulle dans le spectre. Nous utilisons plutôt rng.standard_normal, de moyenne nulle par construction.

Une seule définition, partout. Nous définissons le rapport signal sur bruit (SNR) comme l’amplitude absolue maximale du signal divisée par l’écart-type du bruit. C’est la définition intégrée à mlgeo_synth et employée par les expériences de détection du chapitre 4.3. Garder une définition unique compte : une affirmation comme « le détecteur échoue en dessous d’un SNR de 1 » n’a aucun sens si deux chapitres mesurent le SNR différemment. Dans le code, nous normalisons le signal à une amplitude crête unitaire et le bruit à un écart-type unitaire, de sorte que s + noise / SNR ait exactement le SNR annoncé.

noise = rng.standard_normal(len(s))
plt.figure()
plt.plot(t, noise)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.title('Gaussian noise (zero mean)')
plt.show()

Comparez le spectre d’amplitude de Fourier du bruit à celui du signal.

nhat = fft(noise, n=Nfft)
plt.plot(freqVec, np.abs(nhat[: Nfft // 2]))
plt.plot(freqVec, np.abs(Zhat[: Nfft // 2]))
plt.xscale('log')
plt.legend(['noise', 'signal'])
plt.xlabel('Frequency in Hz')
plt.ylim([0, 100])
plt.grid()

Ils paraissent très différents dans le domaine spectral. Le bruit étale son énergie sur toutes les fréquences ; le signal concentre la sienne dans une bande étroite.

Ajoutons maintenant le bruit au signal, mis à l’échelle par un rapport signal sur bruit (SNR). Nous définissons ici le SNR comme le rapport de l’amplitude absolue maximale du signal à celle du bruit. Normalisons d’abord les deux, puis divisons le bruit par le SNR.

SNR = 100  # signal-to-noise ratio

s /= np.max(np.abs(s))  # unit peak amplitude
noise /= noise.std()    # unit standard deviation, so s + noise / SNR has the stated SNR

noisy_signal = s + noise / SNR
plt.plot(t, noisy_signal)
plt.grid()
plt.title(f'Noisy signal, SNR = {SNR}')
plt.xlabel('Time in s')

3.3 Du bruit d’une couleur choisie : randomiser la phase¶

Le bruit peut avoir un contenu fréquentiel différent, ou couleur. Une recette générale construit une série temporelle de bruit à partir d’un spectre d’amplitude de Fourier choisi :

Choisir l’amplitude à chaque fréquence positive ou nulle (plate pour un bruit blanc, décroissance en $1/f$ pour un bruit coloré, ou le spectre de données réelles).
Tirer une phase aléatoire à chaque fréquence, uniforme entre 0 et $2\pi$ .
Inverser avec np.fft.irfft.

Une série temporelle à valeurs réelles exige un spectre à symétrie hermitienne : la valeur en $-f$ doit être le conjugué complexe de la valeur en $+f$ . Coder cette symétrie à la main est source d’erreurs. irfft ne prend que les fréquences positives ou nulles et impose la symétrie par construction : c’est donc ce que nous utilisons. Les termes de fréquence nulle et de Nyquist doivent être réels, nous mettons donc leurs phases à zéro.

n_rfreq = Nfft // 2 + 1  # number of non-negative frequency bins
phases = rng.uniform(0, 2 * np.pi, n_rfreq)
phases[0] = 0.0    # DC bin must be real
phases[-1] = 0.0   # Nyquist bin must be real (Nfft is even)

amplitude = np.ones(n_rfreq)  # flat spectrum: white noise
amplitude[0] = 0.0            # no DC component

spectrum = amplitude * np.exp(1j * phases)
noise = np.fft.irfft(spectrum, n=Nfft)[: len(s)]

noise /= noise.std()  # unit standard deviation
plt.plot(t, noise)
plt.title('Noise with random phase and white spectrum')

Ajoutez le nouveau bruit au signal (l’ondelette de Ricker), cette fois avec un SNR bien plus faible, et tracez dans les domaines temporel et fréquentiel.

SNR = 1
news = s + noise / SNR

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].plot(t, news)
ax[0].set_title(f'Noisy signal, SNR = {SNR}')
ax[0].set_xlabel('Time in s')
ax[0].grid()
ax[1].plot(freqVec, np.abs(fft(news, n=Nfft)[: Nfft // 2]))
ax[1].set_xscale('log')
ax[1].grid()
ax[1].set_xlabel('Frequency in Hz')
ax[1].set_title('FFT of the noisy signal in log');

Comparez les distributions de données du signal pur, du bruit et de leur somme.

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].hist(noise, bins=20, label='noise')
ax[0].hist(s, bins=20, label='signal')
ax[0].legend()
ax[0].grid()
ax[0].set_title('Signal and noise')
ax[1].hist(news, bins=20)
ax[1].grid()
ax[1].set_title('Noisy signal');

3.4 Moments statistiques¶

Calculez les moments statistiques du signal propre et du bruit. Quels moments discriminent le signal du bruit, et quelle est leur sensibilité au niveau de bruit ?

from scipy.stats import moment

print('The first moment of the signal is:', moment(s, 1))
print('The second moment of the signal is:', moment(s, 2))
print('The third moment of the signal is:', moment(s, 3))
print('The fourth moment of the signal is:', moment(s, 4))

print('The first moment of the noise is:', moment(noise, 1))
print('The second moment of the noise is:', moment(noise, 2))
print('The third moment of the noise is:', moment(noise, 3))
print('The fourth moment of the noise is:', moment(noise, 4))

The first moment of the signal is: 0.0
The second moment of the signal is: 0.013293333067181288
The third moment of the signal is: 0.006430940666682904
The fourth moment of the signal is: 0.007489819567378666
The first moment of the noise is: 0.0
The second moment of the noise is: 1.0
The third moment of the noise is: 0.005114776505715557
The fourth moment of the noise is: 3.134334832790316

Le premier moment centré est nul par définition. Le deuxième moment (la variance) dépend de la normalisation : c’est donc un discriminant faible. Le quatrième moment est le plus intéressant : une série temporelle presque partout nulle, avec un unique transitoire bref, a des queues lourdes par rapport à une gaussienne, donc son kurtosis est grand. Le kurtosis est une caractéristique (feature) de premier plan pour détecter des événements impulsifs dans le bruit, et nous le retrouverons dans la leçon d’ingénierie des caractéristiques (2.11).

4. Bruit synthétique informé par la physique, à partir de données réelles¶

Les bruits blanc et coloré sont des idéalisations. Le vrai bruit sismique a une structure : pics microsismiques, bruit anthropique en journée, réponse instrumentale. Pour rendre le bruit synthétique réaliste, nous pouvons emprunter le spectre d’amplitude d’un bruit réel et n’en randomiser que la phase.

Nous téléchargeons deux heures de bruit sismique enregistré à la station UW.RATT avant un séisme, avec le client FDSN d’obspy.

import obspy
import obspy.clients.fdsn.client as fdsn
from obspy import UTCDateTime

# Download seismic noise data
network = 'UW'
station = 'RATT'
channel = 'HHZ'  # high-gain broadband, vertical component, 100 Hz
Tstart = UTCDateTime(2021, 7, 29, 6, 15)
fdsn_client = fdsn.Client('IRIS')  # client to query the IRIS DMC server
# download the two hours of noise before the event
N = fdsn_client.get_waveforms(network=network, station=station, location='--',
                              channel=channel, starttime=Tstart - 7200,
                              endtime=Tstart, attach_response=True)
N.merge()
N.detrend(type='linear')
N[0].taper(max_percentage=0.05)
print(N)

/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)
/var/folders/js/lzmy975n0l5bjbmr9db291m00000gn/T/ipykernel_25811/286277248.py:8: ObsPyDeprecationWarning: attach_response is deprecated and will be removed in a future release. Use remove_response() instead.
  N = fdsn_client.get_waveforms(network=network, station=station, location='--',

1 Trace(s) in Stream:
UW.RATT..HHZ | 2021-07-29T04:15:00.000000Z - 2021-07-29T06:14:59.990000Z | 100.0 Hz, 720000 samples

4.1 Bruit à spectre apparié¶

Nous engendrons maintenant un bruit synthétique dont le spectre d’amplitude correspond à celui du bruit réel de RATT. Le paquet du cours mlgeo_synth fournit spectrum_matched_noise, qui conserve le spectre d’amplitude d’un enregistrement de référence et tire des phases aléatoires. Il s’appuie sur rfft/irfft, si bien que la symétrie hermitienne nécessaire à une sortie à valeurs réelles est garantie par construction. Exécutez help(mlgeo_synth.spectrum_matched_noise) pour en voir la signature.

import mlgeo_synth

real_noise = N[0].data.astype(float)
real_noise -= real_noise.mean()
fs_real = N[0].stats.sampling_rate

synthetic_noise = mlgeo_synth.spectrum_matched_noise(real_noise, seed=42)
print(real_noise.shape, synthetic_noise.shape)

(720000,) (720000,)

Comparez les deux dans le domaine temporel et dans le domaine fréquentiel. Les séries temporelles doivent paraître différentes — les phases diffèrent, donc les formes d’onde ne coïncident nulle part. Les spectres d’amplitude doivent se superposer exactement, puisque le spectre synthétique est une copie du spectre réel.

t_real = np.arange(len(real_noise)) / fs_real

fig, ax = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
ax[0].plot(t_real, real_noise, lw=0.5)
ax[0].set_title('Real noise (UW.RATT.HHZ)')
ax[0].set_ylabel('Counts')
ax[0].grid()
ax[1].plot(t_real, synthetic_noise, lw=0.5, color='orange')
ax[1].set_title('Spectrum-matched synthetic noise (random phase)')
ax[1].set_xlabel('Time in s')
ax[1].set_ylabel('Counts')
ax[1].grid()
plt.tight_layout()
plt.show()

Nfft_real = next_fast_len(len(real_noise))
freq_real = np.fft.rfftfreq(Nfft_real, d=1 / fs_real)

spec_real = np.abs(np.fft.rfft(real_noise, n=Nfft_real))
spec_synth = np.abs(np.fft.rfft(synthetic_noise, n=Nfft_real))

plt.figure(figsize=(10, 5))
plt.loglog(freq_real[1:], spec_real[1:], lw=2, label='real noise')
plt.loglog(freq_real[1:], spec_synth[1:], lw=0.7, label='synthetic noise')
plt.xlabel('Frequency in Hz')
plt.ylabel('Fourier amplitude')
plt.title('Amplitude spectra: real vs spectrum-matched synthetic noise')
plt.legend()
plt.grid()
plt.show()

Les deux spectres se superposent : l’information d’amplitude est identique. L’information de phase, elle, ne l’est pas. Nous pouvons le vérifier directement.

phase_real = np.angle(np.fft.rfft(real_noise, n=Nfft_real))
phase_synth = np.angle(np.fft.rfft(synthetic_noise, n=Nfft_real))

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].plot(freq_real[1:200], phase_real[1:200], '.', label='real')
ax[0].plot(freq_real[1:200], phase_synth[1:200], '.', label='synthetic')
ax[0].set_xlabel('Frequency in Hz')
ax[0].set_ylabel('Phase [rad]')
ax[0].set_title('Phase spectra (first 200 bins)')
ax[0].legend()
ax[0].grid()
ax[1].hist(phase_synth, bins=30)
ax[1].set_xlabel('Phase [rad]')
ax[1].set_ylabel('Counts')
ax[1].set_title('Synthetic phases: uniform')
ax[1].grid()
plt.tight_layout()
plt.show()

4.2 Augmenter un événement avec du bruit réaliste à divers SNR¶

Nous combinons maintenant les deux idées : un événement synthétique (l’ondelette de Ricker) plus un bruit réaliste à spectre apparié, sur une gamme de niveaux de SNR. C’est une stratégie standard d’augmentation de données : un seul gabarit d’événement propre fournit de nombreux échantillons d’entraînement à niveaux de bruit contrôlés.

spectrum_matched_noise accepte une longueur de sortie n, ce qui nous permet de tirer une réalisation de bruit de la même longueur que la fenêtre de signal.

matched_noise = mlgeo_synth.spectrum_matched_noise(real_noise, n=len(s), seed=7)
matched_noise /= matched_noise.std()  # unit standard deviation, per our SNR definition

snr_levels = [np.power(10, i / 5) for i in range(10)]
fig, ax = plt.subplots(10, 1, figsize=(10, 20), sharex=True)
for i, snr in enumerate(snr_levels):
    ax[i].plot(t, s + matched_noise / snr)
    ax[i].set_title(f'SNR = {snr:.2f}')
    ax[i].grid()
ax[-1].set_xlabel('Time in s')
plt.tight_layout()
plt.show()

Aux faibles SNR, l’événement disparaît dans le bruit. Quelque part le long de cette échelle, tout détecteur — humain ou algorithmique — commence à échouer. Trouver ce point de rupture est exactement ce que nous faisons ensuite, avec le plus ancien détecteur de la boîte à outils sismologique.

4.3 Le générateur de données synthétiques du cours¶

Le paquet mlgeo_synth employé ci-dessus est le générateur de données synthétiques du cours. Outre spectrum_matched_noise, il sait engendrer des sismogrammes synthétiques complets avec une onde P, une onde S, une coda et du bruit coloré, ainsi que des jeux de données étiquetés pour les exercices de classification. Nous nous appuierons dessus dans les chapitres suivants chaque fois qu’il nous faudra des données à vérité terrain connue.

t_syn, trace, meta = mlgeo_synth.synthetic_seismogram(
    duration_s=30.0, fs=100.0, magnitude=3.0, distance_km=30.0, snr=5.0, seed=3)

plt.figure(figsize=(10, 4))
plt.plot(t_syn, trace, lw=0.7)
plt.xlabel('Time in s')
plt.ylabel('Amplitude')
plt.title('mlgeo_synth.synthetic_seismogram: M3.0 at 30 km, SNR 5')
plt.grid()
plt.show()
print(meta)

{'t_p': 10.0, 't_s': 13.571428571428571, 'peak_amplitude': np.float64(0.4841325757612952), 'snr': 5.0}

4.4 Exemple traité : le plancher de détection de STA/LTA¶

Avant qu’un réseau de neurones ne touche à de telles traces (le chapitre 4.3 en entraîne un), nous leur devons un modèle de référence classique. Le détecteur STA/LTA (moyenne à court terme sur moyenne à long terme ; Allen, 1978) fait glisser deux fenêtres sur la trace : une fenêtre courte (ici 1 s) suit le niveau instantané du signal, une fenêtre longue (10 s) suit le fond, et leur rapport s’envole à l’arrivée d’un transitoire. Il tourne dans les chaînes de déclenchement des observatoires depuis des décennies, il n’a rien à entraîner, et tout détecteur appris doit le battre pour justifier sa complexité.

Un détecteur a besoin d’un seuil, et le seuil doit provenir du bruit seul. Nous engendrons 200 fenêtres de bruit pur, relevons le pic de STA/LTA de chacune et plaçons le seuil au 99e centile de cette distribution : un taux de fausses alarmes de 1 % par construction, avant même d’avoir regardé un signal. Nous exécutons ensuite le détecteur sur s + noise / SNR le long de l’échelle de SNR np.logspace(-1, 2, 20), avec 60 réalisations de bruit fraîches par SNR, et comptons la fraction d’essais où le pic de STA/LTA franchit le seuil. Chaque point représente 60 épreuves de Bernoulli, nous lui attachons donc un intervalle binomial (de Wilson).

from obspy.signal.trigger import classic_sta_lta

nsta, nlta = int(1 * fs), int(10 * fs)  # 1 s short window, 10 s long window


def peak_stalta(trace, nsta=nsta, nlta=nlta):
    """Peak STA/LTA of a trace, ignoring the first LTA window (not yet filled)."""
    cft = classic_sta_lta(trace, nsta, nlta)
    return cft[nlta:].max()


# Threshold choice: peak STA/LTA of 200 windows that contain no signal at all
n_noise_trials = 200
noise_peaks = np.array([
    peak_stalta(mlgeo_synth.spectrum_matched_noise(real_noise, n=len(s), seed=50_000 + k))
    for k in range(n_noise_trials)])
threshold = np.quantile(noise_peaks, 0.99)
print(f"noise-only peak STA/LTA: median {np.median(noise_peaks):.2f}, "
      f"99th percentile {threshold:.2f}")

noise-only peak STA/LTA: median 4.34, 99th percentile 6.09

def wilson_interval(k, n, z=1.0):
    """Wilson score interval for a binomial proportion (z=1: roughly 68% coverage)."""
    p = k / n
    denom = 1 + z**2 / n
    center = (p + z**2 / (2 * n)) / denom
    half = z * np.sqrt(p * (1 - p) / n + z**2 / (4 * n**2)) / denom
    return center - half, center + half


snrs_ex = np.logspace(-1, 2, 20)
n_trials = 60
detect_prob_stalta = np.zeros(len(snrs_ex))
err_lo, err_hi = np.zeros(len(snrs_ex)), np.zeros(len(snrs_ex))
for i, snr in enumerate(snrs_ex):
    hits = 0
    for k in range(n_trials):
        noise_k = mlgeo_synth.spectrum_matched_noise(real_noise, n=len(s), seed=1000 * i + k)
        noise_k /= noise_k.std()  # unit noise std; s has unit peak, so the SNR is exact
        hits += peak_stalta(s + noise_k / snr) >= threshold
    detect_prob_stalta[i] = hits / n_trials
    lo, hi = wilson_interval(hits, n_trials)
    err_lo[i], err_hi[i] = detect_prob_stalta[i] - lo, hi - detect_prob_stalta[i]

snr50 = np.interp(0.5, detect_prob_stalta, snrs_ex)
for snr, dp in zip(snrs_ex, detect_prob_stalta):
    print(f"SNR {snr:6.2f}: detection probability {dp:.2f}")
print(f"detection probability crosses 50% near SNR {snr50:.2f}")

fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].hist(noise_peaks, bins=30, color='tab:gray')
ax[0].axvline(threshold, color='tab:red', ls='--',
              label=f'threshold = {threshold:.2f} (99th pct)')
ax[0].set_xlabel('Peak STA/LTA of noise-only windows')
ax[0].set_ylabel('Count')
ax[0].set_title('Threshold choice from noise alone')
ax[0].legend()
ax[1].errorbar(snrs_ex, detect_prob_stalta, yerr=[err_lo, err_hi], marker='o',
               capsize=3, label=f'STA/LTA ({n_trials} trials per SNR)')
ax[1].axvline(snr50, color='tab:red', ls=':', label=f'50% crossing, SNR = {snr50:.2f}')
ax[1].set_xscale('log')
ax[1].set_xlabel('SNR (peak signal amplitude / noise standard deviation)')
ax[1].set_ylabel('Detection probability')
ax[1].set_title('STA/LTA detection floor')
ax[1].grid(alpha=0.3, which='both')
ax[1].legend(loc='lower right')
fig.tight_layout()

SNR   0.10: detection probability 0.02
SNR   0.14: detection probability 0.02
SNR   0.21: detection probability 0.02
SNR   0.30: detection probability 0.02
SNR   0.43: detection probability 0.05
SNR   0.62: detection probability 0.00
SNR   0.89: detection probability 0.02
SNR   1.27: detection probability 0.03
SNR   1.83: detection probability 0.00
SNR   2.64: detection probability 0.07
SNR   3.79: detection probability 0.22
SNR   5.46: detection probability 0.62
SNR   7.85: detection probability 0.85
SNR  11.29: detection probability 0.98
SNR  16.24: detection probability 1.00
SNR  23.36: detection probability 1.00
SNR  33.60: detection probability 1.00
SNR  48.33: detection probability 1.00
SNR  69.52: detection probability 1.00
SNR 100.00: detection probability 1.00
detection probability crosses 50% near SNR 4.97

Le panneau de gauche est la logique du seuil : le bruit pur produit une distribution de valeurs de pic de STA/LTA (médiane autour de 4,3 pour ce bruit riche en microséismes), et placer le seuil à son 99e centile, 6,1, fixe le taux de fausses alarmes à 1 % avant que le moindre signal n’entre dans l’expérience. Le panneau de droite est le plancher de détection : la courbe reste au niveau de fausses alarmes visé de 1 à 2 % jusqu’à un SNR d’environ 2,6, franchit les 50 % vers un SNR de 5 et sature à 1 au-dessus d’un SNR d’environ 16 — une transition large de moins d’une décade. À retenir : avec son taux de fausses alarmes fixé à 1 %, ce détecteur classique a besoin d’une amplitude crête d’environ cinq écarts-types de bruit avant de se déclencher de façon fiable sur cet événement. Ce nombre est le modèle de référence que tout détecteur appris doit battre ; au chapitre 4.3, nous superposons cette courbe au plancher de détection d’un CNN entraîné, mesuré de la même façon.

5. Exercices¶

Exercice 1 : bruit à spectre apparié depuis une autre fenêtre de référence¶

Téléchargez une autre fenêtre de bruit de deux heures à UW.RATT (24 heures plus tôt, par exemple), engendrez-en un bruit à spectre apparié et comparez avec la fenêtre utilisée ci-dessus :

Comparez la variance et le kurtosis des deux fenêtres réelles et de leurs deux versions synthétiques.
Superposez les quatre spectres d’amplitude. Les pics microsismiques se déplacent-ils d’un jour à l’autre ?
Commentez : une seule fenêtre de bruit synthétique suffit-elle à représenter le bruit de cette station ?

# Exercise 1 scaffold
# Tstart2 = Tstart - 86400  # 24 hours earlier
# N2 = fdsn_client.get_waveforms(network=network, station=station, location='--',
#                                channel=channel, starttime=Tstart2 - 7200,
#                                endtime=Tstart2, attach_response=True)
# N2.merge(); N2.detrend(type='linear'); N2[0].taper(max_percentage=0.05)
# real_noise2 = N2[0].data.astype(float) - N2[0].data.mean()
# synthetic_noise2 = mlgeo_synth.spectrum_matched_noise(real_noise2, seed=42)
# from scipy.stats import kurtosis
# ... compare np.var and kurtosis of the four series, then overlay their spectra

Auto-vérification

Les deux fenêtres réelles comptent 720 000 échantillons (deux heures à 100 Hz), et leurs versions synthétiques aussi.
Chaque fenêtre synthétique reproduit fidèlement la variance de sa propre fenêtre de référence ; les deux jours peuvent différer entre eux d’un facteur notable.
scipy.stats.kurtosis renvoie l’excès de kurtosis, donc les fenêtres synthétiques (à phase aléatoire, donc gaussiennes) se situent près de 0 ; les fenêtres réelles peuvent s’en écarter si des transitoires sont présents.
Sur les spectres superposés, chaque courbe synthétique se superpose à sa propre courbe réelle. D’un jour à l’autre, les pics microsismiques (grosso modo 0,1–1 Hz) restent aux mêmes fréquences mais changent d’amplitude.
Partie 3 : les variations de niveau et de forme spectrale d’un jour à l’autre font qu’une seule fenêtre sous-représente le bruit de la station ; tirer de plusieurs fenêtres est plus sûr.

Exercice 2 : déplacer le plancher de détection de STA/LTA¶

Le plancher de détection de la section 4.4 a été mesuré avec un couple de fenêtres particulier : 1 s pour la courte, 10 s pour la longue. Ces longueurs sont les seuls réglages du détecteur, et les observatoires les choisissent en fonction des signaux qu’ils traquent.

Répétez le balayage de la section 4.4 pour trois couples de fenêtres : (0,5 s ; 5 s), (1 s ; 10 s) et (2 s ; 20 s). Pour chaque couple, recalculez d’abord le seuil à partir des fenêtres de bruit pur — la distribution du bruit change avec les fenêtres — puis mesurez la courbe de détection et son franchissement des 50 %.
Tracez les trois courbes de détection sur un même axe. Quel couple de fenêtres détecte l’événement de Ricker au SNR le plus bas, et quel rapport cela a-t-il avec la durée d’environ 1 s de l’événement ?
Qu’est-ce qui déraille si vous gardez le seuil de 1 s / 10 s tout en changeant les fenêtres ?

# Exercise 2 scaffold
# for sta_s, lta_s in [(0.5, 5.0), (1.0, 10.0), (2.0, 20.0)]:
#     n_s, n_l = int(sta_s * fs), int(lta_s * fs)
#     # 1) recompute the threshold: peak_stalta of 200 noise-only windows with
#     #    these window lengths, take the 99th percentile
#     # 2) rerun the SNR sweep of Section 4.4 with peak_stalta(..., n_s, n_l)
#     # 3) record the 50% crossing with np.interp(0.5, detect_prob, snrs_ex)
# plot the three curves on one semilogx axis and compare the crossings

Auto-vérification

Chaque couple de fenêtres a besoin de son propre seuil : la distribution des pics sur bruit pur se décale quand les fenêtres changent (une fenêtre courte plus longue moyenne davantage d’échantillons de bruit, donc ses pics sont plus bas). Attendez-vous à trois 99e centiles différents.
Les trois courbes de détection montent du plancher de fausses alarmes de 1 % jusqu’à 1, et chaque franchissement des 50 % se situe à un facteur deux près de la valeur de la section 4.4 (SNR proche de 5).
La fenêtre courte la plus proche de la durée de l’événement (environ 1 s pour cette ondelette de Ricker) donne le franchissement le plus bas : une STA de 0,5 à 2 s intègre l’essentiel de l’énergie de l’événement, tandis qu’une STA bien plus longue la dilue dans le bruit environnant.
Réutiliser le seuil de 1 s / 10 s avec d’autres fenêtres casse la calibration à 1 % de fausses alarmes : la comparaison entre courbes devient alors déloyale, car chaque détecteur opère à un taux de fausses alarmes différent. Une probabilité de détection ne signifie quelque chose qu’à un taux de fausses alarmes annoncé.

6. Quand les données synthétiques sont-elles recevables ?¶

Les données synthétiques sont un outil, pas un substitut à l’observation. Servez-vous de cette liste de contrôle avant d’introduire des données synthétiques dans un projet.

Recevable :

Développement de méthode. Prototyper un détecteur, un filtre ou un modèle sur des données à vérité terrain connue avant de toucher aux données réelles.
Évaluation par rapport à une référence. Mesurer l’exactitude, le rappel ou les points de rupture (comme dans l’exercice 2) suppose de connaître la vraie réponse. Les données synthétiques la fournissent.
Augmentation à provenance déclarée. Ajouter des réalisations de bruit ou des événements synthétiques à un ensemble d’entraînement est une pratique courante, à condition que l’article ou le rapport indique quels échantillons sont synthétiques et comment ils ont été engendrés.
Ensembles de test cachés. Les enseignants et les responsables de campagnes d’évaluation utilisent des données synthétiques pour bâtir des ensembles de test que les modèles ne peuvent pas avoir mémorisés.

Non recevable :

Remplacer la validation réelle dans une affirmation scientifique. Un modèle validé uniquement sur données synthétiques n’a pas été validé. Une conclusion scientifique sur la Terre doit reposer sur des observations réelles.

Déclarez toujours. Chaque fois que des données synthétiques entrent dans une analyse, dites-le : quels échantillons, quel générateur, quelle graine. Le code du générateur et les graines ont leur place dans le dépôt, avec le reste du projet. Dans ce cours, utilisez le paquet mlgeo_synth comme générateur, pour que la provenance ne soit qu’à un import de distance, et suivez les consignes des instructions du projet MLGeo.