4.3 Réseaux de neurones convolutifs - Machine Learning en géosciences

Les réseaux de neurones convolutifs (CNN) sont les chevaux de trait de l’apprentissage profond sur données structurées : images, spectrogrammes et séries temporelles. Ce carnet couvre :

Les convolutions et les noyaux, appliqués à la main à une image.
Les briques d’un CNN : convolution, activation, pooling (sous-échantillonnage par agrégation), couches entièrement connectées.
LeNet-5 sur les chiffres MNIST, écrit en PyTorch — l’échauffement.
Un CNN 2-D qui régresse des tendances de réchauffement à partir d’un champ climatique synthétique maillé, face à un modèle de référence par moindres carrés.
Un détecteur de séismes par CNN 1-D entraîné sur des sismogrammes synthétiques, son plancher de détection en fonction du rapport signal/bruit, et le détecteur classique STA/LTA mesuré sur les mêmes traces.
Une épreuve de réalité : le détecteur entraîné sur synthétiques, évalué sur de vraies formes d’onde miniPNW étiquetées.
Comment traduire en code fonctionnel le tableau d’architecture d’une publication.

🖥️ Diapositives du cours — Séance 22 (ven. 20 nov.)

import os

import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torchinfo
from scipy.signal import convolve2d
from skimage import data

# Device-agnostic setup: CUDA GPU, Apple silicon (MPS), or CPU.
device = torch.device("cuda" if torch.cuda.is_available()
                      else "mps" if torch.backends.mps.is_available()
                      else "cpu")
print(f"Using device: {device}")

Using device: mps

0. À revoir¶

À ce stade, vous devriez comprendre :

Le perceptron
Le rôle des fonctions d’activation
Pourquoi et comment implémenter une descente de gradient simple
Les fonctions de perte
La constitution des lots (batching)
Les perceptrons multicouches

1. Convolutions¶

Tout comme les perceptrons et les perceptrons multicouches, chaque élément d’un CNN reçoit des entrées, effectue un produit scalaire et fait passer ce produit par une fonction d’activation. La différence est qu’un CNN applique le même petit jeu de poids, appelé noyau, partout dans l’entrée. Ces réseaux fonctionnent bien sur des données structurées 1-D, 2-D et 3-D comme les séries temporelles et les images.

Comme leur nom l’indique, les CNN reposent sur des convolutions :

$(f*g)(x)= \int\limits^{+\infty}_{\tau=-\infty}f(\tau)g(x-\tau)d\tau$

où $g$ est l’entrée, $f$ le noyau de convolution et τ une variable muette.

Pour des fonctions discrètes, on somme simplement :

$(f*g)(x)= \sum\limits^{+\infty}_{\tau=-\infty} f(\tau)g(x-\tau)$

La convolution s’étend à plusieurs dimensions.

Par exemple, considérez le filtre moyenneur (ici, pour un noyau 3 x 3) :

\begin{bmatrix} \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\ \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\ \frac{1}{9} & \frac{1}{9} & \frac{1}{9} \\ \end{bmatrix}

(1)

# Load an example image from scikit-image
image = data.coins()

# Kernel size
filterDimension = 3

# Define a mean filter
mean_filter = np.ones((filterDimension, filterDimension)) / (filterDimension * filterDimension)

# Apply the mean filter to the image using convolution
filtered_image = convolve2d(image, mean_filter, mode='same', boundary='symm')

fig, ax = plt.subplots(1, 2, figsize=(8, 3))
ax[0].imshow(image, cmap='gray')
ax[0].set_title('Original image')
ax[1].imshow(filtered_image, cmap='gray')
ax[1].set_title('Mean-filtered image')
for a in ax:
    a.axis('off')
fig.tight_layout()

Le filtre moyenneur floute. D’autres noyaux détectent la structure. Les filtres de Sobel approximent le gradient de l’image et répondent aux contours :

# Sobel filters
sobel_filter_horizontal = np.array([
    [-1, -2, -1],
    [0, 0, 0],
    [1, 2, 1]
])

sobel_filter_vertical = np.array([
    [-1, 0, 1],
    [-2, 0, 2],
    [-1, 0, 1]
])

# Apply the Sobel filters to the image
result_horizontal = convolve2d(image, sobel_filter_horizontal)
result_vertical = convolve2d(image, sobel_filter_vertical)

# Combine the horizontal and vertical edges to get the overall edges
edges = np.sqrt(result_horizontal**2 + result_vertical**2)

# Thresholding to accentuate the edges
threshold = 10
edges[edges < threshold] = 0

fig, axes = plt.subplots(1, 4, figsize=(13, 3.2))
axes[0].imshow(image, cmap='gray')
axes[0].set_title('Original image')
axes[1].imshow(result_horizontal, cmap='gray')
axes[1].set_title('Horizontal Sobel')
axes[2].imshow(result_vertical, cmap='gray')
axes[2].set_title('Vertical Sobel')
axes[3].imshow(edges, cmap='gray')
axes[3].set_title('Combined edges (thresholded)')
for a in axes:
    a.axis('off')
fig.tight_layout()

L’intérêt d’un CNN est justement de ne pas utiliser des noyaux conçus à la main comme ceux-ci. Il apprend les poids du noyau à partir des données, par descente de gradient. Les premières couches finissent souvent par apprendre d’elles-mêmes des détecteurs de contours.

Exercice 1. Dans la cellule du filtre moyenneur, passez filterDimension de 3 à 15. Qu’arrive-t-il aux pièces de monnaie, et pourquoi ?

Solution

Un filtre moyenneur 15 x 15 moyenne sur un voisinage beaucoup plus grand : l’image devient fortement floutée et les petits détails (la texture à l’intérieur des pièces, les espaces entre elles) disparaissent. La taille du noyau fixe l’échelle spatiale des motifs auxquels une convolution répond.

2. Assembler un CNN¶

L’architecture d’un CNN comprend classiquement :

Une couche convolutive
Une couche d’activation
Une couche de pooling
Une couche entièrement connectée

2.1 La couche convolutive¶

Nous considérons ici un CNN qui prend des images en entrée et renvoie une sortie.

Les aspects de cette couche :

Entrées : dans le cas de convolutions 2D, trois dimensions : hauteur, largeur, profondeur.

Filtres : les noyaux de convolution, chacun avec un biais et un poids pour chaque élément du noyau. La sortie d’un filtre est une carte de caractéristiques (feature map). En général, le nombre de filtres est supérieur à la profondeur d’entrée. La taille du noyau est souvent carrée, mais ce n’est pas obligatoire.

Pas (stride) : le pas de déplacement du filtre.

Remplissage (padding) : des valeurs ajoutées aux bords de l’entrée.

Convolution avec fenêtre de noyau (fig. 7.2.1 de Dive into Deep Learning).

Dans l’exemple ci-dessus, la valeur de sortie en haut à gauche vaut 0×0 + 1×1 + 3×2 + 4×3 = 19. On répète le procédé jusqu’à remplir tous les éléments de la carte de sortie, et on recommence pour chaque filtre.

_{Un bon cours sur les CNN est le pense-bête du cours CS230 de Stanford.}

Comment interpréteriez-vous le code suivant ?

# pytorch version
torch.nn.Conv2d(in_channels=3, out_channels=64, kernel_size=6)

Conv2d(3, 64, kernel_size=(6, 6), stride=(1, 1))

Cette couche prend une entrée à 3 canaux (par exemple une image RGB), applique 64 filtres de taille 6 x 6 (chacun couvrant les 3 canaux d’entrée) et produit 64 cartes de caractéristiques.

2.2 La couche d’activation¶

La couche convolutive est suivie d’une couche d’activation. Celle-ci applique une fonction d’activation, comme ReLU (unité linéaire rectifiée), aux sorties de la convolution.

# Define the ReLU function
def relu(x):
    return np.maximum(0, x)

x = np.linspace(-5, 5, 100)
y_relu = relu(x)

fig, ax = plt.subplots(figsize=(5, 3.5))
ax.plot(x, y_relu, label='ReLU')
ax.set_xlabel('x')
ax.set_ylabel('ReLU(x)')
ax.axhline(0, color='black', linewidth=0.5)
ax.axvline(0, color='black', linewidth=0.5)
ax.grid(color='gray', linestyle='--', linewidth=0.5)
ax.legend()
plt.show()

2.3 La couche de pooling¶

La couche de pooling sous-échantillonne les cartes de caractéristiques. En combinant des valeurs, le pooling réduit la taille du modèle et le rend moins sensible aux petits décalages de l’entrée.

Les couches de max pooling prennent la valeur maximale dans un voisinage donné (fixé par la taille de pooling).

2.4 La couche entièrement connectée¶

Enfin, nous incluons une couche entièrement connectée (chaque entrée connectée à chaque sortie). Pour la classification, les CNN se terminent typiquement par un softmax qui transforme les sorties de la dernière couche entièrement connectée en probabilités de classes (des valeurs entre 0 et 1 dont la somme vaut 1).

2.5 Quelques remarques¶

Dans la couche convolutive, les neurones ne sont pas connectés à toutes les parties des données d’entrée.

Une couche dense apprend des motifs globaux. Une couche de convolution apprend des motifs locaux et est équivariante par translation : décalez l’entrée et la carte de caractéristiques se décale avec elle, si bien qu’un motif appris dans une partie de l’image (ou de la série temporelle) est détecté partout ailleurs. L’invariance par translation approchée d’un CNN complet — la même prédiction où que se trouve le motif — vient des couches de pooling, qui jettent la position à mesure qu’elles sous-échantillonnent ; le détecteur 1-D de la section 5 tient sa tolérance au décalage de son pooling moyen global. Les CNN apprennent aussi des motifs hiérarchiques : une première couche apprend un motif local, une deuxième couche combine les caractéristiques locales en caractéristiques d’échelle plus large.

2.6 Aparté : vocabulaire de la segmentation d’images¶

La classification attribue une étiquette à l’image entière. La segmentation attribue des étiquettes au niveau du pixel, et se décline en trois variantes. L’illustration classique est une scène de rue montrée de trois façons : la segmentation sémantique colore chaque pixel selon sa classe (route, ciel, personne, voiture) sans séparer les individus ; la segmentation d’instances délimite séparément chaque objet individuel (cette personne-ci, cette voiture-là) ; la segmentation panoptique combine les deux, en étiquetant chaque pixel par classe tout en gardant distinctes les instances d’objets.

3. LeNet-5 sur MNIST¶

MNIST est le « hello world » de l’apprentissage profond : petit, propre, et sans rien de commun avec un jeu de données géoscientifique. Nous l’utilisons comme échauffement, pour mettre au point la mécanique d’entraînement sur un problème où rien d’autre ne peut mal tourner ; la charge utile géoscientifique commence à la section 4. Nous construisons maintenant l’architecture LeNet-5 (LeCun et al., 1998), l’un des premiers CNN à succès, et l’entraînons à classifier les chiffres manuscrits de MNIST. Le réseau est une pile séquentielle de 2 couches convolutives et 3 couches entièrement connectées. Une représentation graphique courante :

3.1 Charger les données¶

Nous utilisons torchvision pour télécharger MNIST. L’ensemble d’entraînement complet compte 60 000 images ; pour garder ce carnet rapide, nous entraînons sur un sous-ensemble de 6 000 images et évaluons sur 1 500 images de test. Sur votre propre machine, vous pouvez augmenter ces nombres (et le nombre d’époques) pour une meilleure exactitude.

from torch.utils.data import DataLoader, Subset, TensorDataset
from torchvision import datasets
from torchvision.transforms import Compose, Normalize, ToTensor

data_root = os.path.expanduser("~/.cache/mlgeo-mnist")
transform = Compose([ToTensor(), Normalize([0.5], [0.5])])

train_full = datasets.MNIST(root=data_root, train=True, download=True, transform=transform)
test_full = datasets.MNIST(root=data_root, train=False, download=True, transform=transform)

# Subsets for speed
train_set = Subset(train_full, range(6000))
test_set = Subset(test_full, range(1500))

loaded_train = DataLoader(train_set, batch_size=64, shuffle=True)
loaded_test = DataLoader(test_set, batch_size=256)

X, y = next(iter(loaded_train))
print("One batch of images:", X.shape, "labels:", y.shape)

One batch of images: torch.Size([64, 1, 28, 28]) labels: torch.Size([64])

# Display 3 images from one batch
images, labels = next(iter(loaded_train))
fig, axes = plt.subplots(1, 3, figsize=(8, 3))
for i in range(3):
    axes[i].imshow(images[i].numpy().squeeze(), cmap='gray')
    axes[i].set_title(f"Label: {labels[i].item()}")
    axes[i].axis('off')
plt.show()

3.2 Le modèle¶

Nous écrivons LeNet-5 comme une pile torch.nn.Sequential. Le module Reshape garantit que l’entrée a la forme (batch, channels, height, width), ce qu’attend Conv2d.

class Reshape(torch.nn.Module):
    def forward(self, x):
        return x.view(-1, 1, 28, 28)

model_lenet = torch.nn.Sequential(
    Reshape(),
    nn.Conv2d(1, 6, kernel_size=5, padding=2), nn.Sigmoid(),
    nn.AvgPool2d(kernel_size=2, stride=2),
    nn.Conv2d(6, 16, kernel_size=5), nn.Sigmoid(),
    nn.AvgPool2d(kernel_size=2, stride=2),
    nn.Flatten(),
    nn.Linear(16 * 5 * 5, 120), nn.Sigmoid(),
    nn.Linear(120, 84), nn.Sigmoid(),
    nn.Linear(84, 10),
)

# Walk a dummy input through the network and watch the shape change
Xd = torch.rand(size=(1, 1, 28, 28), dtype=torch.float32)
print('Initial input shape: \t', Xd.shape)
for layer in model_lenet:
    Xd = layer(Xd)
    print(layer.__class__.__name__, 'output shape: \t', Xd.shape)

Initial input shape: 	 torch.Size([1, 1, 28, 28])
Reshape output shape: 	 torch.Size([1, 1, 28, 28])
Conv2d output shape: 	 torch.Size([1, 6, 28, 28])
Sigmoid output shape: 	 torch.Size([1, 6, 28, 28])
AvgPool2d output shape: 	 torch.Size([1, 6, 14, 14])
Conv2d output shape: 	 torch.Size([1, 16, 10, 10])
Sigmoid output shape: 	 torch.Size([1, 16, 10, 10])
AvgPool2d output shape: 	 torch.Size([1, 16, 5, 5])
Flatten output shape: 	 torch.Size([1, 400])
Linear output shape: 	 torch.Size([1, 120])
Sigmoid output shape: 	 torch.Size([1, 120])
Linear output shape: 	 torch.Size([1, 84])
Sigmoid output shape: 	 torch.Size([1, 84])
Linear output shape: 	 torch.Size([1, 10])

torchinfo.summary(model_lenet, input_size=(1, 1, 28, 28))

==========================================================================================
Layer (type:depth-idx)                   Output Shape              Param #
==========================================================================================
Sequential                               [1, 10]                   --
├─Reshape: 1-1                           [1, 1, 28, 28]            --
├─Conv2d: 1-2                            [1, 6, 28, 28]            156
├─Sigmoid: 1-3                           [1, 6, 28, 28]            --
├─AvgPool2d: 1-4                         [1, 6, 14, 14]            --
├─Conv2d: 1-5                            [1, 16, 10, 10]           2,416
├─Sigmoid: 1-6                           [1, 16, 10, 10]           --
├─AvgPool2d: 1-7                         [1, 16, 5, 5]             --
├─Flatten: 1-8                           [1, 400]                  --
├─Linear: 1-9                            [1, 120]                  48,120
├─Sigmoid: 1-10                          [1, 120]                  --
├─Linear: 1-11                           [1, 84]                   10,164
├─Sigmoid: 1-12                          [1, 84]                   --
├─Linear: 1-13                           [1, 10]                   850
==========================================================================================
Total params: 61,706
Trainable params: 61,706
Non-trainable params: 0
Total mult-adds (Units.MEGABYTES): 0.42
==========================================================================================
Input size (MB): 0.00
Forward/backward pass size (MB): 0.05
Params size (MB): 0.25
Estimated Total Size (MB): 0.30
==========================================================================================

3.3 Configuration de l’entraînement¶

Nous devons choisir quelques paramètres d’entraînement :

La métrique d’erreur : l’exactitude
La fonction de perte : l’entropie croisée pour la classification multi-classes
La taille de lot (64 ici)
Le nombre d’époques (2 ici, pour la vitesse)
L’optimiseur : Adam, une variante de la descente de gradient stochastique qui utilise des estimations des premier et second moments du gradient pour adapter le taux d’apprentissage de chaque poids

La fonction d’entraînement ci-dessous est générique : elle accepte n’importe quel modèle et n’importe quelle paire de chargeurs de données, et nous la réutiliserons à la section 5. Elle déplace chaque lot vers device, exécute la passe avant (forward pass), rétropropage la perte et met à jour les poids. Après chaque époque, elle évalue la perte et l’exactitude sur le chargeur de validation (gradients désactivés).

La fonction conserve aussi un checkpoint (point de sauvegarde) : une copie des poids de l’époque à la meilleure exactitude de validation, restaurée à la fin. L’entraînement stochastique ne s’améliore pas de façon monotone : la dernière époque n’est donc pas toujours la meilleure. En production, vous écririez le checkpoint sur disque avec torch.save({'state_dict': model.state_dict()}, 'checkpoint.pt') ; ici nous le gardons en mémoire.

def train_model(model, train_loader, val_loader, n_epochs=2, learning_rate=1e-3, device=device):
    """Train a classifier, record per-epoch metrics, and restore the best-validation checkpoint."""
    model.to(device)
    criterion = nn.CrossEntropyLoss()
    optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)
    history = {"train_loss": [], "val_loss": [], "val_acc": []}
    best_acc, best_state = -1.0, None

    for epoch in range(n_epochs):
        # --- training pass ---
        model.train()
        running_loss = 0.0
        for inputs, labels in train_loader:
            inputs = inputs.float().to(device)
            labels = labels.long().to(device)
            optimizer.zero_grad()
            loss = criterion(model(inputs), labels)
            loss.backward()
            optimizer.step()
            running_loss += loss.item()
        history["train_loss"].append(running_loss / len(train_loader))

        # --- validation pass (no gradients) ---
        model.eval()
        val_loss, correct, total = 0.0, 0, 0
        with torch.no_grad():
            for inputs, labels in val_loader:
                inputs = inputs.float().to(device)
                labels = labels.long().to(device)
                outputs = model(inputs)
                val_loss += criterion(outputs, labels).item()
                correct += (outputs.argmax(dim=1) == labels).sum().item()
                total += labels.size(0)
        history["val_loss"].append(val_loss / len(val_loader))
        history["val_acc"].append(100 * correct / total)

        print(f"[Epoch {epoch + 1}] train loss: {history['train_loss'][-1]:.3f} - "
              f"val loss: {history['val_loss'][-1]:.3f} - val accuracy: {history['val_acc'][-1]:.1f}%")

        # checkpoint: remember the weights of the best epoch so far
        if history["val_acc"][-1] > best_acc:
            best_acc = history["val_acc"][-1]
            best_state = {k: v.detach().clone() for k, v in model.state_dict().items()}

    model.load_state_dict(best_state)
    print(f"Restored checkpoint with best validation accuracy: {best_acc:.1f}%")
    return history

torch.manual_seed(42)
history_lenet = train_model(model_lenet, loaded_train, loaded_test, n_epochs=2, learning_rate=0.005)

[Epoch 1] train loss: 2.272 - val loss: 1.948 - val accuracy: 37.1%

[Epoch 2] train loss: 0.927 - val loss: 0.544 - val accuracy: 84.1%
Restored checkpoint with best validation accuracy: 84.1%

epochs = np.arange(1, len(history_lenet["train_loss"]) + 1)
fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
ax[0].plot(epochs, history_lenet["train_loss"], marker='o', label='train loss')
ax[0].plot(epochs, history_lenet["val_loss"], marker='s', label='validation loss')
ax[0].set_xlabel('Epoch')
ax[0].set_ylabel('Cross-entropy loss')
ax[0].set_xticks(epochs)
ax[0].legend()
ax[1].plot(epochs, history_lenet["val_acc"], marker='o', color='tab:green', label='validation accuracy')
ax[1].set_xlabel('Epoch')
ax[1].set_ylabel('Accuracy (%)')
ax[1].set_xticks(epochs)
ax[1].legend()
fig.suptitle('LeNet-5 on MNIST (6,000-image subset)')
fig.tight_layout()

Deux époques sur un sous-ensemble de 6 000 images font décoller le LeNet à sigmoïdes, mais le laissent loin de son plafond. Avec les 60 000 images complètes et ~10 époques, cette architecture atteint environ 99 % d’exactitude ; essayez sur votre propre machine. Remplacer les sigmoïdes par ReLU accélère aussi considérablement l’entraînement.

L’échauffement est terminé. La recette — des tenseurs en entrée, des blocs de convolution, une tête, une boucle d’entraînement — passe maintenant à des données qui ressemblent aux nôtres.

4. Un CNN 2-D sur un champ géoscientifique maillé¶

Les champs maillés — température de réanalyse, luminances satellitaires, sorties de modèles — sont les données 2-D natives des géosciences, et elles ne se comportent en rien comme des images de chiffres. Nous utilisons le générateur du cours mlgeo_synth.climate_field : des anomalies mensuelles de température sur une grille latitude-longitude de 40 x 80, avec un cycle saisonnier antisymétrique entre hémisphères, un mode zonal lent, un bruit météorologique spatialement corrélé et une tendance linéaire de réchauffement amplifiée aux hautes latitudes nord. Le générateur étant le nôtre, la vraie tendance locale de chaque pixel est connue exactement.

La tâche. Régresser la tendance locale de réchauffement, en °C par décennie, à partir de vignettes (patches) de 8 x 8 pixels. Chaque vignette entre dans le réseau comme une image à 30 canaux : 30 cartes d’anomalies moyennes annuelles, un canal par année. Les canaux-comme-temps sont l’astuce standard pour donner un champ évoluant dans le temps à un CNN 2-D.

Le découpage. Toutes les vignettes d’un même champ partagent la valeur de tendance globale de ce champ : des vignettes du même champ ne doivent donc jamais chevaucher la frontière entraînement/test — la doctrine du découpage par groupes du chapitre 3.8, appliquée. Nous générons 24 champs avec des tendances tirées uniformément entre 0,05 et 0,5 °C/décennie, et découpons par champ : 18 pour l’entraînement, 6 pour le test.

La référence d’abord. Moyennez chaque vignette spatialement et ajustez une droite sur ses 30 moyennes annuelles — un seul appel à np.polyfit. Si le CNN ne fait pas mieux que cela, nous voulons le savoir.

import mlgeo_synth

n_fields, n_years, patch = 24, 30, 8
rng_f = np.random.default_rng(0)
true_trends = rng_f.uniform(0.05, 0.5, n_fields)  # one global trend per field, deg C / decade

X_patches, y_trend, field_id = [], [], []
for fi, trend in enumerate(true_trends):
    field, truth = mlgeo_synth.climate_field(n_lat=40, n_lon=80, n_months=12 * n_years,
                                             trend_c_per_decade=float(trend), seed=100 + fi)
    annual = field.reshape(n_years, 12, 40, 80).mean(axis=1)  # (30, 40, 80) annual means
    # the generator amplifies the trend poleward in the north:
    # local trend = global trend x (1 + 1.5 * max(latitude in radians, 0))
    amp = 1.0 + 1.5 * np.clip(np.deg2rad(truth["lat"]), 0, None)
    local_trend = trend * amp
    for r in range(0, 40, patch):
        for c in range(0, 80, patch):
            X_patches.append(annual[:, r:r + patch, c:c + patch])
            y_trend.append(local_trend[r:r + patch].mean())
            field_id.append(fi)

X_patches = np.array(X_patches, dtype=np.float32)
y_trend = np.array(y_trend, dtype=np.float32)
field_id = np.array(field_id)
print("patches:", X_patches.shape,
      f"| true local trends {y_trend.min():.2f} to {y_trend.max():.2f} C/decade")

train_mask = field_id < 18  # grouped split: fields 0-17 train, 18-23 test
Xc_train, yc_train = X_patches[train_mask], y_trend[train_mask]
Xc_test, yc_test = X_patches[~train_mask], y_trend[~train_mask]
print(f"train: {len(Xc_train)} patches from 18 fields | test: {len(Xc_test)} patches from 6 fields")

# Classical baseline: least-squares slope of the patch-mean annual series
years = np.arange(n_years)

def lsq_trend(patches):
    """Least-squares warming trend (deg C / decade) of each patch's spatial-mean series."""
    series = patches.mean(axis=(2, 3))          # (n, 30)
    slopes = np.polyfit(years, series.T, 1)[0]  # deg C / year
    return slopes * 10.0

mae_lsq = np.abs(lsq_trend(Xc_test) - yc_test).mean()
print(f"least-squares baseline, test MAE: {mae_lsq:.3f} C/decade")

patches: (1200, 30, 8, 8) | true local trends 0.05 to 1.36 C/decade
train: 900 patches from 18 fields | test: 300 patches from 6 fields
least-squares baseline, test MAE: 0.083 C/decade

torch.manual_seed(0)
model_clim = nn.Sequential(
    nn.Conv2d(n_years, 32, kernel_size=3, padding=1), nn.ReLU(),
    nn.Conv2d(32, 32, kernel_size=3, padding=1), nn.ReLU(),
    nn.AdaptiveAvgPool2d(1), nn.Flatten(),
    nn.Linear(32, 1),
)
print(f"parameters: {sum(p.numel() for p in model_clim.parameters()):,}")

clim_loader = DataLoader(TensorDataset(torch.from_numpy(Xc_train), torch.from_numpy(yc_train)),
                         batch_size=64, shuffle=True)
model_clim.to(device)
opt = torch.optim.Adam(model_clim.parameters(), lr=1e-3)
mse = nn.MSELoss()
for epoch in range(30):
    model_clim.train()
    running = 0.0
    for xb, yb in clim_loader:
        xb, yb = xb.to(device), yb.to(device)
        opt.zero_grad()
        loss = mse(model_clim(xb).squeeze(1), yb)
        loss.backward()
        opt.step()
        running += loss.item() * len(xb)
    if (epoch + 1) % 10 == 0:
        print(f"epoch {epoch + 1:2d}: train MSE {running / len(Xc_train):.4f} (C/decade)^2")

parameters: 17,953

epoch 10: train MSE 0.0008 (C/decade)^2

epoch 20: train MSE 0.0001 (C/decade)^2

epoch 30: train MSE 0.0001 (C/decade)^2

model_clim.eval()
with torch.no_grad():
    pred_cnn = model_clim(torch.from_numpy(Xc_test).to(device)).squeeze(1).cpu().numpy()
pred_lsq = lsq_trend(Xc_test)
mae_cnn = np.abs(pred_cnn - yc_test).mean()

# A fairer classical opponent: the generator's slow interannual mode (a 5-year cycle,
# identical in every field) aliases into a straight-line fit. Adding that known mode
# as one extra regressor is all it takes.
zon_annual = np.sin(2 * np.pi * np.arange(12 * n_years) / 60.0).reshape(n_years, 12).mean(axis=1)
G = np.column_stack([np.ones(n_years), years, zon_annual])

def lsq_trend_modeaware(patches):
    """Trend (deg C / decade) with the known interannual mode as a co-regressor."""
    series = patches.mean(axis=(2, 3))
    coef, *_ = np.linalg.lstsq(G, series.T, rcond=None)
    return coef[1] * 10.0

pred_lsq2 = lsq_trend_modeaware(Xc_test)
mae_lsq2 = np.abs(pred_lsq2 - yc_test).mean()

lims = [0, 1.05 * max(yc_test.max(), pred_cnn.max(), pred_lsq.max())]
fig, ax = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True)
for a, pred, name, mae in [(ax[0], pred_lsq, 'Least squares', mae_lsq),
                           (ax[1], pred_lsq2, 'Least squares + known mode', mae_lsq2),
                           (ax[2], pred_cnn, '2-D CNN', mae_cnn)]:
    a.scatter(yc_test, pred, s=12, alpha=0.6)
    a.plot(lims, lims, 'k--', lw=1, label='1:1')
    a.set_xlabel('True local trend (C/decade)')
    a.set_title(f'{name}\ntest MAE {mae:.3f} C/decade')
    a.legend(loc='upper left')
    a.grid(alpha=0.3)
ax[0].set_ylabel('Estimated trend (C/decade)')
fig.tight_layout()
print(f"test MAE (C/decade)  naive least squares: {mae_lsq:.3f}   "
      f"+ known mode: {mae_lsq2:.3f}   2-D CNN: {mae_cnn:.3f}")

test MAE (C/decade)  naive least squares: 0.083   + known mode: 0.009   2-D CNN: 0.008

Les deux estimateurs, appris et classique, retrouvent la structure de tendance amplifiée vers les pôles, mais les marges méritent une lecture attentive. L’ajustement naïf d’une droite aboutit à une erreur absolue moyenne (MAE) de 0,083 °C/décennie tandis que le CNN atteint 0,008 — un facteur dix. Avant de créditer le réseau de magie, demandez-vous ce qu’il a appris. Le générateur superpose à chaque champ un mode zonal lent de 5 ans, et sur un enregistrement de 30 ans ce mode ne se moyenne pas à zéro dans un ajustement linéaire : il se replie dans la pente. Le panneau central donne à l’estimateur classique ce morceau de physique — le mode connu comme simple régresseur supplémentaire — et sa MAE tombe à 0,009 °C/décennie, statistiquement indiscernable du CNN.

D’où la conclusion honnête : la victoire du CNN sur les moindres carrés naïfs est en réalité une victoire sur un modèle de référence sous-spécifié. Ce que le réseau a extrait de 18 champs d’entraînement, c’est le mode interannuel de confusion — ce qu’il n’a pu faire que parce que ce mode est identique dans chaque champ synthétique. Les modèles profonds gagnent leur pain quand les facteurs de confusion sont inconnus ou impossibles à écrire ; quand vous pouvez les écrire, l’ajustement classique avec les bons régresseurs égale le CNN pour une fraction du coût. Notez aussi ce que le découpage par groupes nous a apporté : si des vignettes d’un même champ s’étaient retrouvées des deux côtés de la frontière, le CNN aurait pu mémoriser le bruit propre au champ et la comparaison n’aurait eu aucun sens. Le détecteur de séismes ci-dessous reçoit exactement le même traitement face à sa propre référence classique.

5. Un détecteur de séismes par CNN 1-D¶

Les convolutions ne se limitent pas aux images. Les sismologues utilisent des CNN 1-D pour balayer des sismogrammes continus et classifier de courtes fenêtres en « séisme » ou « bruit ». ConvNetQuake (Perol et al., 2018) en fut un exemple précoce ; les pointeurs automatiques modernes comme EQTransformer et PhaseNet suivent la même idée à plus grande échelle.

Les jeux de données de formes d’onde réelles sont volumineux et leurs étiquettes imparfaites. Nous utilisons ici le générateur de sismogrammes synthétiques du cours, mlgeo_synth. Chaque trace d’événement contient une arrivée P suivie d’une arrivée S plus forte, sur fond de bruit coloré ; chaque trace de bruit ne contient que du bruit coloré. Comme nous contrôlons le générateur, chaque trace vient avec des métadonnées exactes : temps d’arrivée et rapport signal/bruit (SNR, défini comme l’amplitude crête du signal divisée par l’écart-type du bruit).

5.1 Générer et inspecter les données¶

import mlgeo_synth

fs = 100.0  # sampling rate (Hz)
X_seis, y_seis, metas = mlgeo_synth.seismogram_dataset(n_events=800, n_noise=800, fs=fs,
                                                       duration_s=30.0, seed=0)
print("X:", X_seis.shape, "- y:", y_seis.shape)
print("Class balance: %d events, %d noise" % ((y_seis == 1).sum(), (y_seis == 0).sum()))
print("Example event metadata:", {k: round(float(v), 2) for k, v in metas[0].items()})

X: (1600, 3000) - y: (1600,)
Class balance: 800 events, 800 noise
Example event metadata: {'t_p': 9.21, 't_s': 12.21, 'peak_amplitude': 0.5, 'snr': 0.58}

t = np.arange(X_seis.shape[1]) / fs
fig, axes = plt.subplots(4, 1, figsize=(9, 7), sharex=True)

# three events with their P and S arrival times
event_idx = [0, 1, 2]
for ax, i in zip(axes[:3], event_idx):
    ax.plot(t, X_seis[i], color='tab:gray', linewidth=0.8)
    ax.axvline(metas[i]['t_p'], color='tab:blue', linestyle='--', label='P arrival')
    ax.axvline(metas[i]['t_s'], color='tab:red', linestyle='--', label='S arrival')
    ax.set_ylabel('Amplitude')
    ax.set_title(f"Event, SNR = {metas[i]['snr']:.1f}", fontsize=10, loc='left')
axes[0].legend(loc='upper right')

# one noise-only window
noise_idx = np.where(y_seis == 0)[0][0]
axes[3].plot(t, X_seis[noise_idx], color='tab:gray', linewidth=0.8)
axes[3].set_title('Noise window', fontsize=10, loc='left')
axes[3].set_ylabel('Amplitude')
axes[3].set_xlabel('Time (s)')
fig.tight_layout()

Les événements à faible SNR sont difficiles à voir à l’œil. C’est exactement le régime où nous voulons savoir comment se comporte un détecteur.

Avant l’entraînement, nous standardisons chaque trace (retirer sa moyenne, diviser par son écart-type). C’est important : sans cela, le réseau pourrait classifier sur la seule amplitude absolue au lieu de la forme de l’onde.

from sklearn.model_selection import train_test_split

def standardize(X):
    """Per-trace standardization: zero mean, unit standard deviation."""
    X = X - X.mean(axis=1, keepdims=True)
    return X / (X.std(axis=1, keepdims=True) + 1e-10)

Xn = standardize(X_seis).astype(np.float32)

# 70% train, 15% validation, 15% test
X_train, X_tmp, y_train, y_tmp = train_test_split(Xn, y_seis, test_size=0.3,
                                                  random_state=42, stratify=y_seis)
X_val, X_test, y_val, y_test = train_test_split(X_tmp, y_tmp, test_size=0.5,
                                                random_state=42, stratify=y_tmp)

def make_loader(X, y, batch_size, shuffle):
    # unsqueeze(1) adds the channel dimension: (batch, 1, n_samples)
    ds = TensorDataset(torch.from_numpy(X).unsqueeze(1), torch.from_numpy(y).long())
    return DataLoader(ds, batch_size=batch_size, shuffle=shuffle)

train_loader = make_loader(X_train, y_train, 64, shuffle=True)
val_loader = make_loader(X_val, y_val, 256, shuffle=False)
test_loader = make_loader(X_test, y_test, 256, shuffle=False)
print(f"train {len(X_train)}, val {len(X_val)}, test {len(X_test)}")

train 1120, val 240, test 240

5.2 Le modèle¶

L’analogue 1-D du CNN d’images : trois blocs Conv1d → ReLU → MaxPool1d, puis un AdaptiveAvgPool1d qui moyenne chaque carte de caractéristiques en une seule valeur, et une tête linéaire à 2 sorties (bruit, événement). Le pooling moyen global maintient un petit nombre de paramètres et rend le modèle indépendant de la longueur d’entrée.

torch.manual_seed(1)
model_seis = nn.Sequential(
    nn.Conv1d(1, 8, kernel_size=7, padding=3), nn.ReLU(), nn.MaxPool1d(4),
    nn.Conv1d(8, 16, kernel_size=7, padding=3), nn.ReLU(), nn.MaxPool1d(4),
    nn.Conv1d(16, 32, kernel_size=7, padding=3), nn.ReLU(), nn.MaxPool1d(4),
    nn.AdaptiveAvgPool1d(1),
    nn.Flatten(),
    nn.Linear(32, 2),
)
torchinfo.summary(model_seis, input_size=(1, 1, 3000))

==========================================================================================
Layer (type:depth-idx)                   Output Shape              Param #
==========================================================================================
Sequential                               [1, 2]                    --
├─Conv1d: 1-1                            [1, 8, 3000]              64
├─ReLU: 1-2                              [1, 8, 3000]              --
├─MaxPool1d: 1-3                         [1, 8, 750]               --
├─Conv1d: 1-4                            [1, 16, 750]              912
├─ReLU: 1-5                              [1, 16, 750]              --
├─MaxPool1d: 1-6                         [1, 16, 187]              --
├─Conv1d: 1-7                            [1, 32, 187]              3,616
├─ReLU: 1-8                              [1, 32, 187]              --
├─MaxPool1d: 1-9                         [1, 32, 46]               --
├─AdaptiveAvgPool1d: 1-10                [1, 32, 1]                --
├─Flatten: 1-11                          [1, 32]                   --
├─Linear: 1-12                           [1, 2]                    66
==========================================================================================
Total params: 4,658
Trainable params: 4,658
Non-trainable params: 0
Total mult-adds (Units.MEGABYTES): 1.55
==========================================================================================
Input size (MB): 0.01
Forward/backward pass size (MB): 0.34
Params size (MB): 0.02
Estimated Total Size (MB): 0.37
==========================================================================================

Moins de 5 000 paramètres, contre ~62 000 pour LeNet. Les petits modèles s’entraînent vite et surapprennent difficilement sur 1 120 traces d’entraînement.

5.3 Entraîner et évaluer¶

Nous réutilisons train_model de la section 3. Douze époques s’exécutent en quelques secondes ; augmentez le nombre d’époques sur votre propre machine pour gagner encore quelques pour cent.

torch.manual_seed(1)
history_seis = train_model(model_seis, train_loader, val_loader, n_epochs=12, learning_rate=2e-3)

[Epoch 1] train loss: 0.693 - val loss: 0.684 - val accuracy: 50.0%

[Epoch 2] train loss: 0.672 - val loss: 0.637 - val accuracy: 75.4%
[Epoch 3] train loss: 0.600 - val loss: 0.521 - val accuracy: 77.9%

[Epoch 4] train loss: 0.505 - val loss: 0.543 - val accuracy: 71.7%
[Epoch 5] train loss: 0.496 - val loss: 0.430 - val accuracy: 85.8%

[Epoch 6] train loss: 0.468 - val loss: 0.414 - val accuracy: 84.2%
[Epoch 7] train loss: 0.429 - val loss: 0.406 - val accuracy: 87.5%

[Epoch 8] train loss: 0.411 - val loss: 0.377 - val accuracy: 85.8%
[Epoch 9] train loss: 0.419 - val loss: 0.384 - val accuracy: 82.1%

[Epoch 10] train loss: 0.388 - val loss: 0.363 - val accuracy: 87.9%
[Epoch 11] train loss: 0.392 - val loss: 0.390 - val accuracy: 82.1%

[Epoch 12] train loss: 0.376 - val loss: 0.432 - val accuracy: 72.1%
Restored checkpoint with best validation accuracy: 87.9%

epochs = np.arange(1, len(history_seis["train_loss"]) + 1)
fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
ax[0].plot(epochs, history_seis["train_loss"], marker='o', label='train loss')
ax[0].plot(epochs, history_seis["val_loss"], marker='s', label='validation loss')
ax[0].set_xlabel('Epoch')
ax[0].set_ylabel('Cross-entropy loss')
ax[0].legend()
ax[1].plot(epochs, history_seis["val_acc"], marker='o', color='tab:green', label='validation accuracy')
ax[1].set_xlabel('Epoch')
ax[1].set_ylabel('Accuracy (%)')
ax[1].legend()
fig.suptitle('1-D CNN detector on synthetic seismograms')
fig.tight_layout()

def predict(model, X, device=device):
    """Predicted class (0 noise, 1 event) for an array of standardized traces."""
    model.eval()
    with torch.no_grad():
        xb = torch.from_numpy(X.astype(np.float32)).unsqueeze(1).to(device)
        return model(xb).argmax(dim=1).cpu().numpy()

test_pred = predict(model_seis, X_test)
test_acc = (test_pred == y_test).mean()
print(f"Test accuracy: {100 * test_acc:.1f}%")

Test accuracy: 87.9%

Une exactitude dans les 80 % : bien au-dessus du hasard, mais loin de la perfection. Pourquoi pas 100 % ? Le jeu de données inclut délibérément des événements de SNR bien inférieur à 1, quasi indétectables. D’où la vraie question : où le détecteur échoue-t-il ?

5.4 La récompense : le plancher de détection, avec un adversaire classique¶

Avec de vrais sismogrammes, on ne connaît jamais le vrai SNR d’un événement manqué : le comportement à faible SNR d’un détecteur est donc difficile à caractériser. Les synthétiques nous donnent le bouton de réglage que les données réelles n’offrent jamais : nous pouvons générer des événements à n’importe quel SNR et mesurer exactement où la détection s’effondre.

Le CNN n’a pas la scène pour lui seul. Le déclencheur STA/LTA (moyenne à court terme sur moyenne à long terme ; Allen, 1978) — le détecteur classique qui tourne depuis des décennies dans les chaînes de traitement des observatoires, et dont nous avons mesuré le plancher de détection sur des ondelettes de Ricker au chapitre 2.10 — tourne sur les mêmes traces, avec la même définition du SNR (amplitude crête du signal sur écart-type du bruit). Son seuil est fixé de la manière honnête, à partir du bruit seul : le 99e centile du pic de STA/LTA sur les 600 fenêtres de bruit pur, ce qui fixe par construction son taux de fausses alarmes à 1 %.

Le balayage ci-dessous génère 60 événements à chacune de 10 valeurs de SNR d’environ 0,3 à 20 (espacées en logarithme), en variant la magnitude, la distance et la graine de bruit, plus 60 fenêtres de bruit pur appariées par valeur de SNR. Le réseau entraîné tourne en mode inférence ; aucun réentraînement. Chaque probabilité de détection représente 60 épreuves de Bernoulli : les deux courbes portent donc des barres d’erreur binomiales (Wilson).

from obspy.signal.trigger import classic_sta_lta

snrs = np.logspace(-0.5, 1.3, 10)  # ~0.32 to ~20
n_per = 60
rng = np.random.default_rng(42)

nsta, nlta = int(1.0 * fs), int(5.0 * fs)  # 1 s / 5 s: arrivals sit 6-28 s into the window

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

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

detect_prob = np.zeros(len(snrs))   # CNN: fraction of events flagged as events
false_alarm = np.zeros(len(snrs))   # CNN: fraction of noise windows flagged as events
stalta_ev_peaks, stalta_no_peaks = [], []

for i, snr in enumerate(snrs):
    # events at this SNR, with varied magnitude, distance, and noise realization
    traces = []
    for k in range(n_per):
        mag = rng.uniform(1, 4)
        dist = rng.uniform(5, 80)
        _, trace, _ = mlgeo_synth.synthetic_seismogram(magnitude=mag, distance_km=dist,
                                                       snr=float(snr), seed=int(10000 * i + k))
        traces.append(trace)
    traces = np.array(traces)

    # matched noise-only windows
    X_no, _, _ = mlgeo_synth.seismogram_dataset(n_events=0, n_noise=n_per, seed=7000 + i)

    pred_ev = predict(model_seis, standardize(traces))
    pred_no = predict(model_seis, standardize(X_no))
    detect_prob[i] = pred_ev.mean()
    false_alarm[i] = pred_no.mean()

    # STA/LTA peaks on the very same traces (the ratio is amplitude-invariant)
    stalta_ev_peaks.append(np.array([peak_stalta(tr) for tr in traces]))
    stalta_no_peaks.append(np.array([peak_stalta(tr) for tr in X_no]))

# STA/LTA threshold from noise alone: 99th percentile of 600 noise windows -> 1% false alarms
stalta_thresh = np.quantile(np.concatenate(stalta_no_peaks), 0.99)
stalta_prob = np.array([(pk >= stalta_thresh).mean() for pk in stalta_ev_peaks])
stalta_fa = np.array([(pk >= stalta_thresh).mean() for pk in stalta_no_peaks])

print(f"STA/LTA threshold (99th percentile of noise-only peaks): {stalta_thresh:.2f}")
print("   SNR   CNN detect   CNN false alarm   STA/LTA detect")
for snr, d, f, sd in zip(snrs, detect_prob, false_alarm, stalta_prob):
    print(f"{snr:6.2f}   {d:10.2f}   {f:15.2f}   {sd:14.2f}")

STA/LTA threshold (99th percentile of noise-only peaks): 3.11
   SNR   CNN detect   CNN false alarm   STA/LTA detect
  0.32         0.02              0.00             0.02
  0.50         0.02              0.05             0.03
  0.79         0.07              0.00             0.00
  1.26         0.23              0.00             0.03
  2.00         0.92              0.00             0.07
  3.16         1.00              0.02             0.05
  5.01         1.00              0.00             0.57
  7.94         1.00              0.00             1.00
 12.59         1.00              0.00             1.00
 19.95         1.00              0.00             1.00

def binom_err(p, n):
    """Asymmetric Wilson error bars for an array of proportions."""
    lo, hi = np.array([wilson_interval(int(round(pi * n)), n) for pi in p]).T
    return p - lo, hi - p

cnn_lo, cnn_hi = binom_err(detect_prob, n_per)
sta_lo, sta_hi = binom_err(stalta_prob, n_per)

fig, ax = plt.subplots(figsize=(7, 4.2))
ax.errorbar(snrs, detect_prob, yerr=[cnn_lo, cnn_hi], marker='o', capsize=3,
            label='CNN detection probability')
ax.errorbar(snrs, stalta_prob, yerr=[sta_lo, sta_hi], marker='s', capsize=3,
            color='tab:orange', label='STA/LTA detection probability (1% false alarms)')
ax.plot(snrs, false_alarm, color='tab:red', ls='--', lw=1, marker='^', ms=4,
        label='CNN false alarm rate')
ax.plot(snrs, stalta_fa, color='tab:brown', ls=':', lw=1, marker='v', ms=4,
        label='STA/LTA false alarm rate')
ax.set_xscale('log')
ax.set_xlabel('SNR (peak signal amplitude / noise standard deviation)')
ax.set_ylabel('Fraction of windows')
ax.set_ylim(-0.05, 1.05)
ax.set_title('Detection floor: 1-D CNN vs STA/LTA, same traces')
ax.legend(loc='center left', fontsize=9)
ax.grid(alpha=0.3, which='both')
fig.tight_layout()

Lisez d’abord les deux courbes à leurs points de fonctionnement : les taux de fausses alarmes en tireté sont comparables (1 % fixé pour STA/LTA, 0–5 % mesuré pour le CNN), les courbes de détection peuvent donc être comparées directement. Trois régimes. En dessous de SNR ≈ 0,8, les deux détecteurs se rejoignent à leur plancher de fausses alarmes — aucun algorithme, appris ou classique, ne sauve un signal enfoui aussi profond dans un bruit de même couleur. Au-dessus de SNR ≈ 8, ils se rejoignent à nouveau, à 1,0 : là, le déclencheur classique est le bon choix d’ingénierie, parce qu’il détecte tout sans rien coûter à entraîner, régler ou maintenir, et qu’il ne peut pas dériver silencieusement quand la distribution des données change. Le CNN ne mérite sa complexité que dans la bande intermédiaire : à SNR 2, il détecte 92 % des événements quand STA/LTA, au même budget de fausses alarmes, en attrape 7 %, et son croisement à 50 % se situe près de SNR 1,5 contre un croisement près de 5 pour le déclencheur — le même plancher STA/LTA que l’expérience aux ondelettes de Ricker avait trouvé au chapitre 2.10 avec un signal différent et du bruit réel apparié en spectre. Environ une demi-décade de SNR : voilà tout le territoire que gagne le détecteur appris. Que ce territoire compte ou non est une question scientifique, pas architecturale ; en sismologie, il se trouve que c’est là que vivent la plupart des petits séismes, et c’est pourquoi les détecteurs appris ont supplanté les déclencheurs en énergie dans la construction des catalogues.

Ce type de mesure contrôlée est l’argument principal en faveur des bancs d’essai synthétiques : avec des données réelles, on ne peut rapporter la performance que sur les événements qu’on a réussi à cataloguer, et le catalogue lui-même est biaisé au détriment des événements à faible SNR.

Pourquoi pas l’appariement de gabarits (template matching) ? Il existe une méthode classique qui bat ces deux détecteurs à faible SNR : le template matching, qui corrèle une forme d’onde d’événement connue avec les données continues et déclenche sur le coefficient de corrélation. Pour une forme d’onde connue dans un bruit stationnaire, c’est le filtre adapté — statistiquement optimal — et en pratique il détecte des événements répétitifs environ un ordre de grandeur sous le plancher STA/LTA (Gibbons & Ringdal, 2006) ; c’est ainsi que les séismes basse fréquence ont été extraits du trémor (Shelly et al., 2007). Le piège est dans le nom : il lui faut un gabarit, il ne trouve donc que des événements qui répètent une forme d’onde déjà en votre possession — séquences de répliques, essaims, séismes répétitifs, sismicité induite. Il répond à une question plus étroite (« cette source a-t-elle rompu de nouveau ? ») que les détecteurs ci-dessus (« un événement, quel qu’il soit, est-il présent ? »), et l’évaluer équitablement exigerait un jeu de données de sources répétitives, ce qui explique qu’il reste hors du périmètre travaillé ici.

Exercice 2. À partir des résultats du balayage, estimez le SNR auquel la probabilité de détection du CNN franchit 50 %, et faites de même pour STA/LTA. Puis réentraînez le modèle avec seulement 2 époques et relancez le balayage. Le plancher de détection du CNN bouge-t-il ? Celui de STA/LTA ?

Solution

Lisez les croisements sur le tableau imprimé ou sur la figure (vous pouvez interpoler avec np.interp(0.5, detect_prob, snrs), les courbes étant à peu près monotones) ; le CNN croise à un SNR plus bas que le déclencheur classique. Avec seulement 2 époques d’entraînement, le réseau est sous-entraîné : sa transition se décale vers des SNR plus élevés et son plateau supérieur peut rester sous 1,0, si bien que son plancher de détection remonte. La courbe STA/LTA, elle, ne bouge pas du tout — elle n’a aucun paramètre à sous-entraîner, et c’est exactement ce qui en fait un modèle de référence stable. Entraîner plus longtemps, ajouter des données ou augmenter avec des exemples à faible SNR abaisse le plancher du CNN, jusqu’à la limite fixée par le bruit lui-même.

6. Épreuve de réalité : le détecteur évalué sur de vraies formes d’onde miniPNW¶

Tous les chiffres jusqu’ici ont été mesurés sur le générateur même qui a produit les données d’entraînement. La question inconfortable : que se passe-t-il sur de vrais sismogrammes ?

Le chapitre 2.11 a téléchargé miniPNW, un sous-ensemble étiqueté du jeu de référence « prêt pour l’IA » du Pacifique Nord-Ouest (Ni et al., 2023) : de vraies formes d’onde avec pointés d’analystes et types de sources. Nous évaluons sur quelques centaines de fenêtres le CNN entraîné sur données synthétiques de vrais séismes, tel quel — pas de réentraînement, pas d’affinage, pas d’ajustement de seuil. Le prétraitement reproduit exactement celui de l’entraînement et rien de plus : une fenêtre de 30 s de la composante verticale, découpée pour que le pointé P tombe à 7 s du début (dans la plage utilisée par les synthétiques), puis la même standardisation par trace. Les fenêtres de bruit proviennent des 30 premières secondes de chaque trace, qui se terminent bien avant l’arrivée P. Quel que soit le score, nous le rapportons.

import urllib.request
import h5py
import pandas as pd

# Same files and loader pattern as notebook 2.11, which caches them in its data/ folder
pnw_dir = os.path.join('..', 'Chapter2-DataManipulation', 'data')
os.makedirs(pnw_dir, exist_ok=True)
metadata_path = os.path.join(pnw_dir, 'miniPNW_metadata.csv')
waveform_path = os.path.join(pnw_dir, 'miniPNW_waveforms.hdf5')
base_url = 'https://dasway.ess.washington.edu/shared/niyiyu/PNW-ML'

have_pnw = True
try:
    if not os.path.exists(metadata_path):
        urllib.request.urlretrieve(f'{base_url}/miniPNW_metadata.csv', metadata_path)
    if not os.path.exists(waveform_path):
        print('Downloading miniPNW waveforms (about 670 MB, one-time; shared with notebook 2.11)...')
        urllib.request.urlretrieve(f'{base_url}/miniPNW_waveforms.hdf5', waveform_path)
except Exception as err:
    have_pnw = False
    print('miniPNW cache is missing and the download failed, so the reality check below is skipped.\n'
          'Run notebook 2.11 first (it downloads and caches the files), then rerun this section.\n'
          f'Reason: {err!r}')
have_pnw = have_pnw and os.path.exists(waveform_path)
print('miniPNW available:', have_pnw)

miniPNW available: True

if have_pnw:
    meta = pd.read_csv(metadata_path)
    eq = meta[(meta['source_type'] == 'earthquake') & meta['trace_P_arrival_sample'].notna()].copy()
    n_win, pre = 3000, 700  # 30 s at 100 Hz; P pick 7 s into the window
    p_samp = eq['trace_P_arrival_sample'].astype(int)
    # keep traces where the event window and a non-overlapping pre-P noise window both fit
    eq = eq[(p_samp - pre >= n_win) & (p_samp - pre + n_win <= 15001)].head(300)

    def read_z(f, trace_name):
        """Vertical component of one miniPNW trace (same reader as notebook 2.11)."""
        bucket, narray = trace_name.split('$')
        x, _, z = (int(v) for v in narray.split(',:'))
        return f['/data/' + bucket][x, 2, :z]  # channel order N, E, Z

    X_real_ev = np.zeros((len(eq), n_win), dtype=np.float64)
    X_real_no = np.zeros((len(eq), n_win), dtype=np.float64)
    snr_db = np.zeros(len(eq))
    with h5py.File(waveform_path, 'r') as f:
        for i, (_, row) in enumerate(eq.iterrows()):
            tr = read_z(f, row['trace_name']).astype(np.float64)
            pk = int(row['trace_P_arrival_sample'])
            X_real_ev[i] = tr[pk - pre : pk - pre + n_win]
            X_real_no[i] = tr[:n_win]  # the trace starts 50 s before the pick: pre-event noise
            snr_db[i] = np.mean([float(v) for v in str(row['trace_snr_db']).split('|')])

    # drop dead windows (flat traces cannot be standardized)
    alive = (X_real_ev.std(axis=1) > 0) & (X_real_no.std(axis=1) > 0)
    X_real_ev, X_real_no, snr_db = X_real_ev[alive], X_real_no[alive], snr_db[alive]
    print(f"real earthquake windows: {len(X_real_ev)}, matched pre-event noise windows: {len(X_real_no)}")

real earthquake windows: 300, matched pre-event noise windows: 300

if have_pnw:
    pred_real_ev = predict(model_seis, standardize(X_real_ev))
    pred_real_no = predict(model_seis, standardize(X_real_no))
    real_det = pred_real_ev.mean()
    real_fa = pred_real_no.mean()
    real_acc = 0.5 * (real_det + (1 - real_fa))

    print(f"synthetic test accuracy (Section 5.3):        {100 * test_acc:.1f}%")
    print(f"real miniPNW detection rate (earthquakes):    {100 * real_det:.1f}%")
    print(f"real miniPNW false alarm rate (pre-P noise):  {100 * real_fa:.1f}%")
    print(f"real miniPNW balanced accuracy:               {100 * real_acc:.1f}%")
    print(f"synthetic-to-real gap (balanced accuracy):    "
          f"{100 * (test_acc - real_acc):.1f} percentage points")

    # detection rate vs catalog SNR (quartile bins), with binomial error bars
    edges = np.quantile(snr_db, [0, 0.25, 0.5, 0.75, 1.0])
    centers, det_b, err_b = [], [], []
    for a, b in zip(edges[:-1], edges[1:]):
        m = (snr_db >= a) & (snr_db <= b)
        k, n = int(pred_real_ev[m].sum()), int(m.sum())
        lo, hi = wilson_interval(k, n)
        centers.append(0.5 * (a + b))
        det_b.append(k / n)
        err_b.append((k / n - lo, hi - k / n))
    err_b = np.array(err_b).T
    for ctr, d, (a, b) in zip(centers, det_b, zip(edges[:-1], edges[1:])):
        print(f"catalog SNR {a:5.1f} to {b:5.1f} dB: detection rate {d:.2f}")

    fig, ax = plt.subplots(figsize=(6, 3.8))
    ax.errorbar(centers, det_b, yerr=err_b, marker='o', capsize=3, color='tab:blue',
                label='real earthquakes (quartile bins)')
    ax.axhline(real_fa, color='tab:red', ls='--', lw=1,
               label=f'false alarm rate on real noise ({100 * real_fa:.0f}%)')
    ax.set_xlabel('Catalog SNR (dB, channel average)')
    ax.set_ylabel('Detection rate')
    ax.set_ylim(-0.05, 1.05)
    ax.set_title('Synthetic-trained CNN on real miniPNW waveforms')
    ax.legend(loc='center right', fontsize=9)
    ax.grid(alpha=0.3)
    fig.tight_layout()

synthetic test accuracy (Section 5.3):        87.9%
real miniPNW detection rate (earthquakes):    86.3%
real miniPNW false alarm rate (pre-P noise):  92.7%
real miniPNW balanced accuracy:               46.8%
synthetic-to-real gap (balanced accuracy):    41.1 percentage points
catalog SNR -12.1 to   1.7 dB: detection rate 0.99
catalog SNR   1.7 to   4.8 dB: detection rate 0.85
catalog SNR   4.8 to  11.2 dB: detection rate 0.80
catalog SNR  11.2 to  59.1 dB: detection rate 0.81

L’écart est le résultat : 87,9 % d’exactitude équilibrée sur les traces de test synthétiques, 46,8 % sur les fenêtres réelles miniPNW — le niveau du hasard, une chute de 41 points de pourcentage. Le mode de défaillance est spécifique et mérite d’être nommé. Les fenêtres de vrais séismes sont signalées à 86 %, ce qui paraît honorable isolément, mais le détecteur signale aussi comme séismes 93 % des fenêtres de vrai bruit pré-événement. La courbe par classes de SNR dit la même chose plus crûment : la détection est la plus haute (99 %) dans le quartile de SNR le plus bas, où la fenêtre n’est presque que du bruit, et la courbe entière est indiscernable du taux de fausses alarmes. Le réseau ne répond pas du tout aux séismes. Il a appris « événement synthétique contre le seul modèle de bruit synthétique de son entraînement », et le vrai bruit du Pacifique Nord-Ouest — non stationnaire, dominé par le microséisme océanique — ne ressemble ni à l’un ni à l’autre : presque tout déclenche.

Nous rapportons ce chiffre tel que mesuré : pas de bande passante choisie pour le corriger, pas de réentraînement, pas de seuil déplacé après avoir vu les données de test. Faire disparaître l’écart par réglage effacerait la leçon, qui est la règle d’admissibilité du chapitre 2.10 rendue concrète : un modèle validé uniquement sur données synthétiques n’a pas été validé, et voilà à quoi ressemble cette phrase une fois transformée en mesure. Combler l’écart exige des données réelles côté entraînement — entraîner ou affiner sur des formes d’onde miniPNW étiquetées (elles attendent dans le cache), ou informer par la physique le bruit du générateur avec le bruit apparié en spectre de 2.10. Le carnet 4.6 mène l’expérience complémentaire : quelle part d’un encodeur pré-entraîné sur synthétiques survit au contact des mêmes données réelles.

7. Lire et recoder des réseaux publiés¶

Supposons qu’un article de recherche décrive l’architecture du CNN que ses auteurs ont utilisé pour leur analyse, mais ne fournisse aucun code. Comment le reproduiriez-vous ?

Considérez cet article :

Rouet-Leduc, B., Hulbert, C., McBrearty, I. W., Johnson, P. A. (2020). Probing slow earthquakes with deep learning. Geophysical Research Letters, 47, e2019GL085870. Rouet‐Leduc et al. (2020)

(La figure 1 de Rouet-Leduc et al. (2020) montre l’architecture du réseau ; voir l’article à Rouet‐Leduc et al. (2020) — la figure n’est pas reproduite ici pour des raisons de licence.)

Schéma du CNN et de son architecture (figure 1 de Rouet-Leduc et al., 2020).

En lisant l’article et son matériel supplémentaire, nous pouvons lister les couches :

Entrée : image de spectrogramme de 129 x 95 x 1 pixels
Conv2D : noyau de taille 16 x 16, profondeur 32 (nombre de canaux), produisant des cartes de caractéristiques de taille 114 x 80 ; activation ReLU (trouvée dans le matériel supplémentaire)
Max pooling de taille 2
Dropout de 5 % (trouvé dans le matériel supplémentaire)
Conv2D : noyau de taille 8 x 8, profondeur 64
Max pooling de taille 2
Dropout de 5 % (trouvé dans le matériel supplémentaire)
Conv2D : noyau de taille 4 x 4, profondeur 128
Couche entièrement connectée (dense) aplatissant en 36 608 valeurs (trouvé dans le matériel supplémentaire)
Couche entièrement connectée (dense) de 10 neurones
Couche entièrement connectée (dense) de 1 neurone, à activation sigmoïde, donnant la probabilité que la fenêtre contienne du trémor

Nous traduisons maintenant cette liste, ligne à ligne, en une pile torch.nn.Sequential :

model_tremor = torch.nn.Sequential(
    nn.Conv2d(in_channels=1, out_channels=32, kernel_size=16),
    nn.ReLU(),
    nn.MaxPool2d(kernel_size=2),
    nn.Dropout(0.05),
    nn.Conv2d(in_channels=32, out_channels=64, kernel_size=8),
    nn.ReLU(),
    nn.MaxPool2d(kernel_size=2),
    nn.Dropout(0.05),
    nn.Conv2d(in_channels=64, out_channels=128, kernel_size=4),
    nn.Flatten(),
    nn.Linear(36608, 10),
    nn.Sigmoid(),
    nn.Linear(10, 1),
    nn.Sigmoid(),
)

Pour vérifier la traduction, nous exécutons une passe avant sur un tenseur factice à la forme d’entrée de l’article, (batch, channels, height, width) = (1, 1, 129, 95), et inspectons la forme après chaque couche :

Xd = torch.randn(1, 1, 129, 95)
print('Input shape: \t\t', Xd.shape)
for layer in model_tremor:
    Xd = layer(Xd)
    print(layer.__class__.__name__, 'output shape: \t', Xd.shape)

Input shape: 		 torch.Size([1, 1, 129, 95])
Conv2d output shape: 	 torch.Size([1, 32, 114, 80])
ReLU output shape: 	 torch.Size([1, 32, 114, 80])
MaxPool2d output shape: 	 torch.Size([1, 32, 57, 40])
Dropout output shape: 	 torch.Size([1, 32, 57, 40])
Conv2d output shape: 	 torch.Size([1, 64, 50, 33])
ReLU output shape: 	 torch.Size([1, 64, 50, 33])
MaxPool2d output shape: 	 torch.Size([1, 64, 25, 16])
Dropout output shape: 	 torch.Size([1, 64, 25, 16])
Conv2d output shape: 	 torch.Size([1, 128, 22, 13])
Flatten output shape: 	 torch.Size([1, 36608])
Linear output shape: 	 torch.Size([1, 10])
Sigmoid output shape: 	 torch.Size([1, 10])
Linear output shape: 	 torch.Size([1, 1])
Sigmoid output shape: 	 torch.Size([1, 1])

torchinfo.summary(model_tremor, input_size=(1, 1, 129, 95))

==========================================================================================
Layer (type:depth-idx)                   Output Shape              Param #
==========================================================================================
Sequential                               [1, 1]                    --
├─Conv2d: 1-1                            [1, 32, 114, 80]          8,224
├─ReLU: 1-2                              [1, 32, 114, 80]          --
├─MaxPool2d: 1-3                         [1, 32, 57, 40]           --
├─Dropout: 1-4                           [1, 32, 57, 40]           --
├─Conv2d: 1-5                            [1, 64, 50, 33]           131,136
├─ReLU: 1-6                              [1, 64, 50, 33]           --
├─MaxPool2d: 1-7                         [1, 64, 25, 16]           --
├─Dropout: 1-8                           [1, 64, 25, 16]           --
├─Conv2d: 1-9                            [1, 128, 22, 13]          131,200
├─Flatten: 1-10                          [1, 36608]                --
├─Linear: 1-11                           [1, 10]                   366,090
├─Sigmoid: 1-12                          [1, 10]                   --
├─Linear: 1-13                           [1, 1]                    11
├─Sigmoid: 1-14                          [1, 1]                    --
==========================================================================================
Total params: 636,661
Trainable params: 636,661
Non-trainable params: 0
Total mult-adds (Units.MEGABYTES): 329.27
==========================================================================================
Input size (MB): 0.05
Forward/backward pass size (MB): 3.47
Params size (MB): 2.55
Estimated Total Size (MB): 6.07
==========================================================================================

Deux vérifications par rapport à l’article : la première convolution produit des cartes de caractéristiques de 114 x 80, exactement comme indiqué, et le vecteur aplati qui entre dans les couches denses compte 128 x 22 x 13 = 36 608 valeurs, ce qui correspond aux « 36 608 neurones » du matériel supplémentaire. Quand ces nombres concordent, le recodage est cohérent avec le tableau publié.

Notez que nous n’entraînons pas ce réseau. L’exercice consiste ici à traduire la description d’une architecture publiée en code fonctionnel et à la confronter aux formes et aux nombres de paramètres rapportés par l’article. L’entraîner exigerait le jeu de données de l’article (des années de données sismiques et GPS continues de Cascadia) et des moyens de calcul considérables — ce cours n’a ni l’un ni l’autre, et aucun des deux n’est nécessaire pour apprendre le geste.

8. Régler les réseaux CNN¶

Il y a beaucoup d’hyperparamètres et de choix de modèle à faire :

Entraînement : taux d’apprentissage, optimiseur, taille de lot, fonction de perte, régularisation
Architecture : nombre de couches, nombre de canaux (profondeur) par couche, tailles de noyaux, fonctions d’activation, normalisation par lots, dropout

Nous revenons au réglage systématique des hyperparamètres, y compris la recherche automatisée avec Optuna, au carnet 4.5 (entraînement de modèles).

Résumé¶

Une convolution fait glisser un petit noyau sur l’entrée ; des noyaux conçus à la main floutent ou détectent les contours, et un CNN apprend ses noyaux à partir des données.
Un bloc de CNN, c’est convolution, activation, pooling ; une tête de classification, une ou plusieurs couches entièrement connectées.
LeNet-5 en PyTorch tient en une douzaine de lignes ; un détecteur de séismes par CNN 1-D est encore plus petit.
Sur un champ climatique maillé, un petit CNN 2-D bat d’un facteur dix l’estimation naïve de tendance par moindres carrés — mais seulement parce qu’il a appris le mode interannuel de confusion du générateur ; donnez ce mode comme régresseur à l’ajustement classique et les deux s’égalent à ~0,01 °C/décennie.
Mesuré sur les mêmes traces à des taux de fausses alarmes comparables, STA/LTA franchit 50 % de détection près de SNR 5 et le CNN près de 1,5 ; la bande intermédiaire est tout l’argumentaire du détecteur appris.
Évalué tel quel sur de vraies formes d’onde miniPNW, le détecteur entraîné sur synthétiques tombe de 87,9 % à 46,8 % d’exactitude équilibrée — le hasard — parce que le vrai bruit riche en microséisme le déclenche en permanence. L’écart synthétique-réel est la mesure ; la validation sur synthétiques seuls n’est pas une validation.
Un tableau d’architecture publié plus des vérifications de formes avec torchinfo.summary suffisent à recoder un réseau dont vous n’avez pas le code source.

La suite : les réseaux récurrents pour données séquentielles au carnet 4.4.

References¶

Gibbons, S. J., & Ringdal, F. (2006). The detection of low magnitude seismic events using array-based waveform correlation. Geophysical Journal International, 165(1), 149–166. 10.1111/j.1365-246x.2006.02865.x
Shelly, D. R., Beroza, G. C., & Ide, S. (2007). Non-volcanic tremor and low-frequency earthquake swarms. Nature, 446(7133), 305–307. 10.1038/nature05666
Ni, Y., Hutko, A., Skene, F., Denolle, M., Malone, S., Bodin, P., Hartog, R., & Wright, A. (2023). Curated Pacific Northwest AI-ready Seismic Dataset. Seismica, 2(1). 10.26443/seismica.v2i1.368
Rouet‐Leduc, B., Hulbert, C., McBrearty, I. W., & Johnson, P. A. (2020). Probing Slow Earthquakes With Deep Learning. Geophysical Research Letters, 47(4). 10.1029/2019gl085870