Traitement de signaux pour la détection de séismes et leur multilatération

Dalibard Louis 15 mai 2026

Table des matières

1 Séismes  2
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Ondes P et S . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Différence de vitesse de propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Théorie  3
2.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 DSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.1 Acquisition des données . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.2 Extraction des temps d’arrivée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2.3 Convolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2.4 Propriétés algébriques de la convolution . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.5 Moyennage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.6 Transformée de Fourier discrète . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.7 Radix-2 decimation-in-time (DIT) - Factorisation Cooley-Tukey . . . . . . . . . . . . . 12
2.2.8 IFFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.9 Passage du domaine temporel au domaine fréquentiel . . . . . . . . . . . . . . . . . . 14
2.2.10 Corrélation croisée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Modélisation de la propagation des ondes sismiques . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.1 Calcul du temps de propagation selon iasp91 . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.2 Tabulation et interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 Multilatération . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 Fonction d’erreur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.2 Méthode de Nelder-Mead . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Magnitude sismique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3 Toutes les pièces mises ensembles  19
3.1 Fonctionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4 Précision et erreur  23
4.1 Cas réduit avec propagation rectiligne à longues distances et 4 detecteurs . . . . . . . . . . . 23
4.2 Méthode de Monte-Carlo pour les cas réels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5 Bonus  23
6 Bibliographie  23

1 Séismes

1.1 Introduction

Un séisme ou tremblement de terre est une secousse du sol résultant de la libération brusque d’énergie accumulée par les contraintes exercées sur les roches. Cette libération d’énergie se fait par rupture le long d’une faille, généralement préexistante. Plus rares sont les séismes dus à l’activité volcanique ou d’origine artificielle (explosions par exemple). Le lieu de la rupture des roches en profondeur se nomme le foyer ; la projection du foyer à la surface est l’épicentre du séisme. Le mouvement des roches près du foyer engendre des vibrations élastiques qui se propagent, sous la forme de paquets d’ondes sismiques, autour et au travers du globe terrestre.

1.2 Ondes P et S

Les ondes sismiques, ou ondes élastiques, sont des mouvements vibratoires qui se propagent à travers un milieu matériel et peuvent le modifier irréversiblement si leur amplitude est suffisante. Elles sont engendrées par un événement initial, généralement un séisme. L’impulsion de départ déplace des atomes, qui en poussent d’autres à leur tour avant de reprendre leur place, ces déplacements oscillatoires se propageant ensuite de proche en proche. Un séisme émet des ondes dans toutes les directions de l’espace.

La physique des milieux élastiques est régie par l’équation de Navier, qui implique l’existence de deux grands types d’ondes, mises en évidence expérimentalement : les ondes de volume qui traversent l’intérieur de la Terre et les ondes de surface qui se propagent dans une couche d’épaisseur limitée en suivant la surface terrestre. Les ondes de volume se subdivisent par ailleurs en ondes longitudinales (ondes P) et transversales (ondes S).

Figure 1 – Ondes P  Figure 2 – Ondes S

1.3 Différence de vitesse de propagation

Les ondes P ou ondes primaires, appelées aussi ondes de compression (en anglais Pressure waves) ou ondes longitudinales. Le déplacement du sol qui accompagne leur passage se fait par des dilatations et des compressions successives. Ces déplacements du sol sont parallèles à la direction de propagation de l’onde. Elles se propagent dans tous les milieux et sont les plus rapides (autour de 6 km/s près de la surface), parcourant le chemin le plus court, même par noyau terrestre, et sont donc les premières enregistrées sur les sismogrammes. Elles sont responsables du grondement sourd que l’on peut entendre au début d’un tremblement de terre.

Les ondes S ou ondes secondaires, appelées aussi ondes de cisaillement (Shear waves) ou ondes transversales.

À leur passage, les mouvements du sol s’effectuent perpendiculairement au sens de propagation de l’onde. Ces ondes ne se propagent pas dans les milieux liquides, elles sont en particulier arrêtées par le noyau externe de la Terre. Leur vitesse est d’environ 4 km/s. Elles apparaissent en deuxième sur les sismogrammes.

Figure 3 – 6 km/s (ondes P) vs 4 km/s (ondes S)

On peut donc exploiter les temps d’arrivée de ces différentes ondes pour localiser un séisme.

2 Théorie

2.1 Principe

Différentes étapes :

1. Acquisition de données en temps réel (SeedLink)

2. Reconnaissance d’un séisme et mesure automatique des temps

3. Calcul de la position et de la magnitude (a) Modélisation de la propagation des ondes sismiques

(b) Méthode numérique d’optimisation de fonction à plusieurs variables pour la multilatération (c) Calcul de la magnitude

2.2 DSP

2.2.1 Acquisition des données

On utilise le protocole Seedlink pour récupérer les données des sismographes en temps réel.

Name  Host
AusPass  auspass.edu.au
BGR  eida.bgr.de
CISMID  www.cismid.uni.edu.pe
ENS  ephesite.ens.fr
GEOFON, GFZ  geofon.gfz-potsdam.de
GEONET  link.geonet.org.nz
Geoscience Australia  seis-pub.ga.gov.au
GSRAS ( ?)  89.22.182.133
Helsinki  finseis.seismo.helsinki.fi
Haiti  ayiti.unice.fr
ICGC  ws.icgc.cat
IDA Project  rtserve.ida.ucsd.edu
IFZ  data.ifz.ru
IPGP  rtserver.ipgp.fr
IRIS DMC  rtserve.iris.washington.edu
IRIS Jamaseis  jamaseis.iris.edu
ISNET - UNINA  185.15.171.86
LMU  erde.geophysik.uni-muenchen.de
NIGGG  195.96.231.100

NRCAN  earthquakescanada.nrcan.gc.ca
OBSEBRE  obsebre.es
OGS  nam.ogs.it
Oklahoma University  rtserve.ou.edu
ORFEUS  eida.orfeus-eu.org
PLSN (IGF Poland)  hudson.igf.edu.pl
Red Sìsmica de Puerto Rico  161.35.236.45
Red Sìsmica Baru  helis.redsismicabaru.com
RESIF  rtserve.resif.fr
SANET  147.213.113.73
RSIS  rsis1.on.br
SCSN-USC (South Carolina Seismic Network)  eeyore.seis.sc.edu :6382
Seisme IRD  rtserve.ird.nc
Staneo  vibrato.staneo.fr
SNAC NOA  snac.gein.noa.gr
TexNet  rtserve.beg.utexas.edu
Thai Meteorological Department  119.46.126.38
UFRN (Universidade Federal do Rio Grande do Norte) sislink.geofisica.ufrn.br
Unical Universita Della Calabria  www.sismocal.org
UNITS Università degli studi di Trieste  rtweb.units.it
UNIV-AG Université des Antilles  seedsrv0.ovmp.martinique.univ-ag.fr
Universidade de Évora  clv-cge.uevora.pt
Universidad de Colima  148.213.24.15
UPR  worm.uprm.edu
USGS  cwbpub.cr.usgs.gov
USP-IAG  seisrequest.iag.usp.br

2.2.2 Extraction des temps d’arrivée

Figure 4 – Exemple d’un enregistrement de sismographe

On a un signal comme ça, on peut facilement voir l’onde P et l’onde S. Mais comment faire pour extraire ces temps d’arrivée automatiquement ? Notamment quand on aura des centaines de mesures venant de partout dans le monde ? Pour cela on va faire un détour par des maths.

2.2.3 Convolution Si on a deux fonctions, quelles opérations peut-on faire pour avoir une nouvelle fonction ? La plus simple semble être les additionner.

−6 −5 −4 −3 −2 −1 123456−1.5 −1 −0.5 0.5 1 f(x) = sin x g(x) = cos 2x h(x) = f(x) + g(x) x y  \begin{matrix}f(x) = sinx \\ g(x) = cos2x \\ h(x) = f(x) + g(x) \end{matrix}

Une autre façon serait de les multiplier.

−6 −5 −4 −3 −2 −1 123456−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 f(x) = sin x g(x) = cos 2x h(x) = f(x) · g(x) x y  \begin{matrix}f(x) = sinx \\ g(x) = cos2x \\ h(x) = f(x) · g(x) \end{matrix}

Mais il existe aussi une autre opération, la convolution.

−6 −5 −4 −3 −2 −1 123456−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 f(x) = sin x g(x) = cos 2x h(x) = (f ∗ g)(x) x y  \begin{matrix}f(x) = sinx \\ g(x) = cos2x \\ h(x) = (f ∗ g)(x) = 0 \end{matrix}

Qui a de nombreuses applications dans le traitement d’images, la théorie de la probabilité, dans la solution d’équations différencielles (température), et comme nous le verrons (aussi la multiplication de polynômes). On définit formellement la convolution comme un produit sur l’espace des fonctions intégrables. Soit f et g deux fonctions intégrables sur ℝ. Pour tout t ∈ ℝ,

(f ∗ g)(t) = ∫_{−∞}^{+∞} f(x) · g(t − x)dx

Dans le cas discret on utilisera une somme. Soit f et g deux fonctions de ℤ dans ℂ. Pour tout n ∈ ℤ,

(f ∗ g)[n] = ∑_{m=−∞}^{+∞} f[m] · g[n − m]

et dans le cas de fonction périodiques, il suffit de faire l’intégrale sur une période. Pour tout t ∈ ℝ, (f ∗ g)(t) = ∫_{0}^{T} f(x) · g(t − x)dx

(on peut également généraliser cela à des fonctions discretès périodiques)

2.2.4 Propriétés algébriques de la convolution

— Commutatif

On remarquera que si

(f ∗ g)(t) = ∫_{−∞}^{+∞} f(x) · g(t − x)dx

Et on fait le changement de variable u = t − x On a

(f ∗ g)(t) = ∫_{+∞}^{−∞} f(t − u) · g(u) − du = ∫_{−∞}^{+∞} f(t − u) · g(u)du = (g ∗ f)(t)

— Associatif

(f ∗ g) ∗ h = f ∗ (g ∗ h)

C’est une conséquence du théorème de Fubini.

— Distributif

f ∗ (g + h) = f ∗ g + f ∗ h

Par linéarité de l’intégrale. L’espace des fonctions intégrables muni de ∗ forme un demi-groupe commutatif (car pas d’élement neutre). Pourquoi une telle construction ? Déjà à quoi ça ressemble ? Si on voulait compter le nombre de combinaisons de dés qui donnent un certain total, on les disposerait en liste, on metterait une d’elle à l’envers et on regarderait le nombre de dés qui s’alignent à chaque fois.

Somme donnant 2 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 3 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 4 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 5 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 6 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 7 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 8 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 9 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 10 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 11 :

1 2 3 4 5 6 6 5 4 3 2 1 Somme donnant 12 :

1 2 3 4 5 6 6 5 4 3 2 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 1 2 3 4 5 6 5 4 3 2 1 Somme Fréquence

Distribution des sommes

On obtient les fréquences d’apparition des différentes combinaisons. Nous pouvons également remarquer que si à la place de rajouter 1 à chaque fois, on ponderait les probabilités et faisait le produit des probabilités qui se trouvent face à face, on aurait également le résultat pour le cas non équiprobable.

On vient de faire la convolution entre cette fonction et elle même.

−1 1 2 3 4 5 6 7 8 1 x f(x)

Fonction nulle puis constante égale à 1 puis nulle

Figure 5 – numpy confirme ce résultat

On remarque que si on le fait dans le cas discret, on a des valeurs en plus sur les bords que l’on peut ignorer. Faire la convolution revient donc à retourner une des listes des coefficients et le résultat est une fonction qui à un décalage t entre les deux listes associe une somme des produits deux à deux des coefficients qui se retrouvent en face l’un de l’autre par ce décalage.

−6 −5 −4 −3 −2 −1 123456−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 x y  \begin{matrix}cos(2x) \\ sin(x) \end{matrix}

−6 −5 −4 −3 −2 −1 123456−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 x y  \begin{matrix}cos(2x) \\ sin(−x) \end{matrix}

−6 −5 −4 −3 −2 −1 123456−0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 x y  cos(2x) sin(−x)

On voit que c’est nul dans ce cas. 2.2.5 Moyennage Si on applique la fonction \frac{f}{6} à une fonction g quelconque on observe une sorte de moyennage qui se fait entre 6 valeurs proches.

0246810 12 14 16 18 20 22 02468x y Moyennage par convolution \begin{matrix}g \\ g ∗ \frac{f}{6} \end{matrix}

Ce processus est fait en deux dimensions sur des images avec differentes fonctions pour avoir des effets comme la détection de bord ou le floutage. 2.2.6 Transformée de Fourier discrète Si on doit faire la convolution de deux listes a et b telles que a = [a_{0}, a_{1}, ..., a_{n−1} ] et b = [b_{0}, b_{1}, ..., b_{k−1} ],

On a (a ∗ b)[p] = ∑_{\begin{matrix}i∈ℤ \\ 0≤p−i≤k−1 \\ 0≤i≤n−1 \end{matrix}} a_{i}b_{p−i}

On reconnait la formule du produit des polynômes P (X) = ∑ _{i=0}^{n−1} a_{i}X^{i} et Q(X) = ∑ _{j=0}^{k−1} a_{j}X^{j} (a ∗ b)[p] est le coefficient du terme de degré p dans le produit

P Q(X) = ∑ _{j=0}^{p} ( ∑_{\begin{matrix}i∈ℤ \\ 0≤p−i≤k−1 \\ 0≤i≤n−1 \end{matrix}} a_{i}b_{j−i} )X^{j}

Bon, de premier abord, ça n’a pas l’air de beaucoup nous aider, le produit de polynômes est O(n^{2}). Mais pour un polynôme, du fait de leur rigidité, on peut juste évaluer leur produit en un nombre fini de points pour retrouver les coefficients. Et donc l’idée si on avait un nombre discret de coefficients serait de les traiter comme les coefficients d’un polynôme, on les évalue un certain nombre de fois et on multiplie ensemble les évaluations, puis on retrouve les coefficients en résolvant le système de Vandermonde correspondant. Le problème c’est que cela semble prendre plus de temps, un pivot de Gauss est en O(n^{3}) Mais si on choisit d’évaluer nos valeurs en les racines de l’unité, on a un système plus simple à résoudre. On va pouvoir faire les opérations en O(n · log(n)) Les coefficients que l’on obtient quand on évalue en les ω_{n}^{k} = e^{− \frac{2kiπ}{n}} valeurs sont appellés la transformation de Fourier discrète TFD (un équivalent discret de la transformation de Fourier continue). Il existe aussi d’autres algorithmes.

2.2.7 Radix-2 decimation-in-time (DIT) - Factorisation Cooley-Tukey Evaluer notre polynôme P (x) = ∑ _{i=0}^{n−1} a_{i}x^{i} en les ω_{n}^{k} revient à faire la multiplication matricielle suivante :

R = [ f_{0} ; f_{1} ; ⋮ ; f_{n−1} ] = [ P (ω_{n}^{0}) ; P (ω_{n}^{1}) ; ⋮ ; P (ω_{n}^{n−1} )] = [ \begin{matrix}1 & 1 & 1 & ⋯ & 1 \\ 1 & ω_{n}^{1} & ω_{n}^{2} & ⋯ & ω_{n}^{n−1} \\ 1 & ω_{n}^{2} & ω_{n}^{4} & ⋯ & ω_{n}^{2(n−1)} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 1 & ω_{n}^{n−1} & ω_{n}^{2(n−1)} & ⋯ & ω_{n}^{(n−1)(n−1)} \end{matrix} ] [ a_{0} ; a_{1} ; ⋮ ; a_{n−1} ]

On reconnaît la matrice de Vandermonde transposée associée à l’évaluation en les ω_{n}^{k} On se restreint au cas où n = 2^{p} (p ∈ ℕ) pour simplifier l’explication même si en théorie ça ne change pas grand chose vu qu’il y a des algorithmes généraux possibles et on peut rajouter des 0 à la fin pour avoir une liste de taille 2^{p}. On note la matrice de Vandermonde transposée, qui permet de calculer le DFT pour tous nos coefficients,

F_{2^{p}} = [ \begin{matrix}1 & 1 & 1 & ⋯ & 1 \\ 1 & ω_{n}^{1} & ω_{n}^{2} & ⋯ & ω_{n}^{n−1} \\ 1 & ω_{n}^{2} & ω_{n}^{4} & ⋯ & ω_{n}^{2(n−1)} \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ 1 & ω_{n}^{n−1} & ω_{n}^{2(n−1)} & ⋯ & ω_{n}^{(n−1)(n−1)} \end{matrix} ]

On peut permuter les lignes, à la fin, ce qu’on a c’est que, du fait des propriétés des racines de l’unité,