\documentclass[a4paper, 11pt]{article} \usepackage[a4paper, total={18cm, 24cm}]{geometry} \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage[french]{babel} \usepackage{lmodern} \usepackage{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsthm} \usepackage{graphicx} \usepackage{color} \usepackage{xcolor} \usepackage{url} \usepackage{theorem} \usepackage{textcomp} \usepackage{listings} \usepackage{hyperref} %\usepackage{glossaries} \usepackage{parskip} \usepackage{amsmath, amsthm, amssymb} \usepackage{stmaryrd} \usepackage{graphicx} \usepackage{subfig} \usepackage{color} \usepackage{longtable} \usepackage{pgfplots} \usepackage{nicematrix} \usepackage[table]{xcolor} \usepackage{tikz} \usepackage{tikz-3dplot} \usetikzlibrary{arrows.meta, 3d} \usepackage{assets/texpackages/annotate-equations} \newenvironment{graytext}{\color{gray}}{\ignorespacesafterend} \graphicspath{ {./assets/} } % Define custom colors \definecolor{mmi1}{HTML}{FFFFFF} \definecolor{mmi2}{HTML}{BFCCFF} \definecolor{mmi3}{HTML}{A0E6FF} \definecolor{mmi4}{HTML}{80FFFF} \definecolor{mmi5}{HTML}{7AFF93} \definecolor{mmi6}{HTML}{FFFF00} \definecolor{mmi7}{HTML}{FFC800} \definecolor{mmi8}{HTML}{FF9100} \definecolor{mmi9}{HTML}{FF0000} \definecolor{mmi10}{HTML}{C80000} \definecolor{mmi11}{HTML}{A40000} \definecolor{mmi12}{HTML}{800000} \begin{document} \title{Traitement de signaux pour la détection de séismes et leur multilatération} \author{Dalibard Louis} \date{\today} \maketitle \tableofcontents \newpage \section{Séismes} \subsection{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 s'exerçant sur les roches. Cette libération d'énergie se produit 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 où se produit la rupture des roches en profondeur est appelé 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 à travers le globe terrestre.\\ \includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{eq-ed-fault-labeled.png} \subsection{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). \noindent\begin{minipage}{.5\textwidth} \centering \includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_compression_impulsion_1d_30_petit.png} \captionof{figure}{Ondes P} \end{minipage}% \begin{minipage}{.5\textwidth} \centering \includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_cisaillement_impulsion_1d_30_petit.png} \captionof{figure}{Ondes S} \end{minipage} \subsection{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. \begin{center} \includegraphics[width=7cm, trim={0 12.5mm 0cm 2cm}, clip]{2025-03-26T23:29:53,572838335+01:00.png}\\ \captionof{figure}{6 km/s (ondes P) vs 4 km/s (ondes S)} \end{center} On peut donc exploiter les temps d'arrivée de ces différentes ondes pour localiser un séisme. \section{Théorie} \subsection{Principe} Différentes étapes: \begin{enumerate} \item Acquisition de données en temps réel (SeedLink) \item Reconnaissance d'un séisme et mesure automatique des temps \item Calcul de la position et de la magnitude \begin{enumerate} \item Modélisation de la propagation des ondes sismiques \item Méthode numérique d'optimisation de fonction à plusieurs variables pour la multilatération \item Calcul de la magnitude \end{enumerate} \end{enumerate} \subsection{DSP} \subsubsection{Acquisition des données} On utilise le protocole Seedlink pour récupérer les données des sismographes en temps réel. \renewcommand{\arraystretch}{1.2} \begin{longtable}{|l|l|} \hline \textbf{Name} & \textbf{Host} \\ \hline 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 \\ \hline \end{longtable} \subsubsection{Extraction des temps d'arrivée} \begin{center} \includegraphics[width=15cm, trim={0 0 0 0}, clip]{R0ED0.EHZ-1743142835749.9944-cropped.png}\\ \captionof{figure}{Exemple d'un enregistrement de sismographe} \end{center} 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. \subsubsection{Convolution} Si on a deux fonctions, quelles opérations peut-on faire pour avoir une nouvelle fonction? La plus simple semble être les additionner. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=300, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] % Define the functions \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$}; \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$}; \addplot[green, thick] {sin(deg(x)) + cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) + g(x)$}; % Add legend \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) + g(x)$} \end{axis} \end{tikzpicture} \end{center} Une autre façon serait de les multiplier. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=300, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] % Define the functions \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$}; \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$}; \addplot[green, thick] {sin(deg(x)) * cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) \cdot g(x)$}; % Add legend \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) \cdot g(x)$} \end{axis} \end{tikzpicture} \end{center} Mais il existe aussi une autre opération, la convolution. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=300, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] % Define the functions \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$}; \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$}; \addplot[green, thick] {0} node[right] {\footnotesize $h(x) = (f*g)(x)$}; % Add legend \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = (f*g)(x) = 0$} \end{axis} \end{tikzpicture} \end{center} 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 $\mathbb{R}$. Pour tout $t \in \mathbb{R}$, \[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \] Dans le cas discret on utilisera une somme. Soit $f$ et $g$ deux fonctions de $\mathbb{Z}$ dans $\mathbb{C}$. Pour tout $n \in \mathbb{Z}$, \[ (f*g)[n]=\sum_{m=-\infty}^{+\infty} f[m] \cdot 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 \in \mathbb{R}$, \[ (f*g)(t)=\int_{0}^{T} f(x) \cdot g(t-x) dx \] (on peut également généraliser cela à des fonctions discretès périodiques) \subsubsection{Propriétés algébriques de la convolution} \begin{itemize} \item Commutatif On remarquera que si \[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \] Et on fait le changement de variable $u=t-x$ On a \[ (f*g)(t)=\int_{+\infty}^{-\infty} f(t-u) \cdot g(u) -du=\int_{-\infty}^{+\infty} f(t-u) \cdot g(u) du=(g*f)(t)\] \item Associatif \[ (f*g)*h= f*(g*h) \] C'est une conséquence du théorème de Fubini. \item Distributif \[ f*(g+h)= f*g+f*h \] Par linéarité de l'intégrale. \end{itemize} 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. \begin{center} \begin{tikzpicture}[scale=0.5] % Define Dice Faces \def\diceA{1,2,3,4,5,6} % Top Die (Static) \def\diceB{6,5,4,3,2,1} % Bottom Die (Sliding) % Define step spacing \def\stepSpace{3} % Loop through each step \foreach \step in {-5,-4,-3,-2,-1,0,1,2,3,4,5} { \pgfmathsetmacro{\sommedonnant}{int(7 + \step)} % Label step number \node[anchor=east] at (-0.5, -\step*3) {Somme donnant \sommedonnant:}; % Top Row - Fixed Dice \foreach \x [count=\i] in \diceA { \node[draw, minimum size=0.5cm] at (\i, -\step*3) {\x}; } % Bottom Row - Sliding Dice \foreach \x [count=\i] in \diceB { \pgfmathsetmacro{\sumval}{int(\i + \step)} \ifnum \sumval > 0 \ifnum \sumval < 7 \node[draw=red, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x}; \else \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x}; \fi \else \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x}; \fi } } \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ xlabel={Somme}, ylabel={Fréquence}, ymin=0, ymax=7, xmin=1.5, xmax=12.5, xtick={2,3,4,5,6,7,8,9,10,11,12}, ytick={1,2,3,4,5,6,7}, area style, width=12cm, height=8cm, major grid style={line width=.2pt,draw=gray!50}, grid=major, bar width=0.8, title={Distribution des sommes}, nodes near coords, ] \addplot[fill=blue!40, draw=blue!80, opacity=0.8] coordinates { (2,1) (3,2) (4,3) (5,4) (6,5) (7,6) (8,5) (9,4) (10,3) (11,2) (12,1) }; \end{axis} \end{tikzpicture} \end{center} 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. \begin{center} \begin{tikzpicture} \begin{axis}[ xlabel={$x$}, ylabel={$f(x)$}, xmin=-1, xmax=8, ymin=-0.2, ymax=1.4, xtick={-1,0,1,2,3,4,5,6,7,8}, ytick={0,1}, width=10cm, height=6cm, grid=major, title={Fonction nulle puis constante égale à $1$ puis nulle}, samples=100, domain=-1:8, axis lines=middle, clip=false, ] % Plot the function \addplot[blue, thick, const plot] coordinates { (-1, 0) (0.999, 0) (1, 1) (6, 1) (6.001, 0) (8, 0) }; % Add markers to explicitly show the endpoints \addplot[only marks, mark=*, mark options={fill=red}] coordinates { (1, 1) (6, 1) }; \end{axis} \end{tikzpicture} \end{center} \begin{center} \includegraphics[width=14cm, trim={0 0 0 0}, clip]{2025-04-01T23:42:09,818697014+02:00.png} \captionof{figure}{numpy confirme ce résultat} \end{center} 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. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=300, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(2*x))}; \addlegendentry{$\cos(2x)$} \addplot[red, thick] {sin(deg(x))}; \addlegendentry{$\sin(x)$} \end{axis} \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=300, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(2*x))}; \addlegendentry{$\cos(2x)$} \addplot[red, thick] {sin(deg(-x))}; \addlegendentry{$\sin(-x)$} \end{axis} \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=300, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(2*x))*sin(deg(-x))}; \addlegendentry{$\cos(2x)\sin(-x)$} % Highlight positive areas in red \addplot [fill=red, opacity=0.3, domain=-2*pi:2*pi] {max(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle; % Highlight negative areas in blue \addplot [fill=blue, opacity=0.3, domain=-2*pi:2*pi] {min(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle; \end{axis} \end{tikzpicture} \end{center} On voit que c'est nul dans ce cas. \subsubsection{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. \begin{center} \begin{tikzpicture} \begin{axis}[ xlabel={$x$}, ylabel={$y$}, xmin=-0.5, xmax=23.5, ymin=-0.5, ymax=9, xtick={0,2,4,6,8,10,12,14,16,18,20,22}, width=12cm, height=8cm, grid=major, title={Moyennage par convolution}, legend pos=north west, ] % Plot the first dataset with connected lines \addplot[blue, thick] coordinates { (0,0) (1,0) (2,0) (3,0) (4,1) (5,2) (6,3) (7,4) (8,5) (9,6) (10,7) (11,8) (12,0) (13,4) (14,2) (15,3) (16,8) (17,1) (18,2) (19,4) (20,8) (21,0) (22,0) (23,0) }; % Plot the second dataset with connected lines \addplot[red, thick] coordinates { (0,0) (1,0.16666667) (2,0.5) (3,1) (4,1.66666667) (5,2.5) (6,3.5) (7,4.5) (8,5.5) (9,5) (10,5) (11,4.5) (12,4) (13,4.16666667) (14,3) (15,3.33333333) (16,3.33333333) (17,4.33333333) (18,3.83333333) (19,2.5) (20,2.33333333) (21,2) (22,1.33333333) (23,0) }; \legend{$g$, $g*\frac{f}{6}$} \end{axis} \end{tikzpicture} \end{center} 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. \subsubsection{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] = \sum\limits_{\substack{i\in\mathbb{Z} \\ 0\leq p-i\leq k-1\\ 0\leq i\leq n-1}} a_ib_{p-i}$ On reconnait la formule du produit des polynômes $P(X)=\sum_{i=0}^{n-1} a_i X^i$ et $Q(X)=\sum_{j=0}^{k-1} a_j X^j$ $(a*b)[p]$ est le coefficient du terme de degré $p$ dans le produit $PQ(X)=\sum_{j=0}^p(\sum\limits_{\substack{i\in\mathbb{Z} \\ 0\leq p-i\leq k-1\\ 0\leq i\leq n-1}} a_ib_{j-i})X^j$ Bon, de premier abord, ça n'a pas l'air de beaucoup nous aider, le produit de polynômes est $\mathcal{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 $\mathcal{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 $\mathcal{O}(n \cdot log(n))$ Les coefficients que l'on obtient quand on évalue en les $\omega_n^k = e^{-\frac{2ki\pi}{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. \subsubsection{Radix-2 decimation-in-time (DIT) - Factorisation Cooley-Tukey} Evaluer notre polynôme $P(x)=\sum_{i=0}^{n-1} a_i x^i$ en les $\omega_n^k$ revient à faire la multiplication matricielle suivante: \[ R=\begin{bmatrix} f_0 \\ f_1 \\ \vdots \\ f_{n-1} \end{bmatrix}=\begin{bmatrix} P(\omega_n^0) \\ P(\omega_n^1) \\ \vdots \\ P(\omega_n^{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} \] On reconnaît la matrice de Vandermonde transposée associée à l'évaluation en les \( \omega_n^k \) On se restreint au cas où $n=2^p$ ($p \in \mathbb{N}$) 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{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix}$ 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é, si on note, $D_{2^{p-1}} = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & \omega_{2^{p-1}} & 0 & \cdots & 0 \\ 0 & 0 & \omega_{2^{p-1}}^2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \omega_{2^{p-1}}^{2^{p-1}-1} \end{bmatrix}$ est: \begin{equation*} R=F_{2^p}\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}= \begin{pNiceArray}{cc|cc} \eqnmarkbox[blue]{mcoeffdiag1}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag2}{D_{2^{p-1}}} \\ \hline \eqnmarkbox[blue]{mcoeffdiag3}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag4}{-D_{2^{p-1}}} \end{pNiceArray} \begin{pNiceArray}{cc|cc} \eqnmarkbox[cyan]{mrecurrence1}{F_{2^{p-1}}} && \mathbf{0} \\ \hline \mathbf{0} && \eqnmarkbox[cyan]{mrecurrence2}{F_{2^{p-1}}} \end{pNiceArray} \begin{bmatrix} \eqnmarkbox[green]{mcoeffpairs1}{a_0} \\ \eqnmarkbox[green]{mcoeffpairs2}{a_2} \\ \vdots \\ \eqnmarkbox[green]{mcoeffpairs3}{a_{2^{p-1}-2}} \\ \eqnmarkbox[green]{mcoeffpairs4}{a_{2^{p-1}}} \\ \eqnmarkbox[magenta]{mcoeffimpairs1}{a_{1}} \\ \eqnmarkbox[magenta]{mcoeffimpairs2}{a_{3}} \\ \vdots \\ \eqnmarkbox[magenta]{mcoeffimpairs3}{a_{n-3}} \\ \eqnmarkbox[magenta]{mcoeffimpairs4}{a_{n-1}} \end{bmatrix} \end{equation*} \annotatetwo[yshift=1.5em]{above}{mcoeffdiag1}{mcoeffdiag2}{Blocs diagonaux, de l'ordre de $\mathcal{O}(n)$ opérations} \annotate[yshift=0.5em]{above}{mcoeffpairs1}{Coefficients d'indice pair} \annotate[yshift=0em]{below}{mcoeffimpairs4}{Coefficients d'indice impair} \annotate[yshift=0em]{below,left}{mrecurrence2}{Récurrence} On évalue la compléxité de l'algorithme.\\ On note $u_p$ sa complexité en fonction de $p$ et $C_n$ sa complexité en fonction de $n$. \begin{equation*} u_{p+1}=\eqnmarkbox[magenta]{cdiagmult}{A \cdot 2^{p+1}}+\eqnmarkbox[red]{crec}{2u_{p}} \end{equation*} \annotate[yshift=0.5em]{above,left}{cdiagmult}{Compléxité du produit sur la diagonale} \annotate[yshift=0em]{below,right}{crec}{Traitement des coefficients par récurrence} On factorise par la solution homogène, $\frac{u_{p+1}}{2^{p+1}}=A+\frac{u_{p}}{2^p}$ $\frac{u_{p}}{2^{p}}=u_{0}+A\cdot p$ $u_{p}=u_{0}\cdot2^{p}+A\cdot p \cdot 2^{p}$ Or $p=\log_2{n}$ Donc $C_n = u_{\log_2{n}} = u_{0} \cdot n+A \cdot \log_2{n} \cdot n = \mathcal{O}(n \log n)$ \subsubsection{IFFT} On peut montrer que: $F_{2^p}^{-1}=\frac{1}{2^p}\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \overline{\omega_n}^1 & \overline{\omega_n}^2 & \cdots & \overline{\omega_n}^{n-1} \\ 1 & \overline{\omega_n}^2 & \overline{\omega_n}^4 & \cdots & \overline{\omega_n}^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \overline{\omega_n}^{n-1} & \overline{\omega_n}^{2(n-1)} & \cdots & \overline{\omega_n}^{(n-1)(n-1)} \end{bmatrix}$ $A=\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}=F_{2^p}^{-1}R=F_{2^p}^{-1}\begin{bmatrix} f_0 \\ f_1 \\ \vdots \\ f_{n-1} \end{bmatrix}$ Le conjugué passe au produit et à la somme, donc aussi pour les matrices (prendre le conjugué d'une matrice c'est prendre le conjugué des termes de la matrice). $A=\overline{\overline{A}}=\overline{\overline{F_{2^p}^{-1}R}}=\overline{\overline{F_{2^p}^{-1}}\overline{R}}=\frac{1}{2^p}\overline{F_{2^p}\overline{R}}$ On peut donc utiliser la même technique, en prenant le conjugué avant d'appliquer un FFT et en le prenant après puis en renormalisant. \subsubsection{Passage du domaine temporel au domaine fréquentiel} Sans passer trop de temps dessus, le fait de multiplier par ces coefficients spécifiques, revient à décomposer en ondes sinusoidales de différentes fréquences et phases notre signal. Il y a beaucoup d'applications dans d'autres domaines. Mais ce qui est remarquable c'est que pour calculer plus rapidement notre résultat. On est passé par une autre réprésentation. Maintenant qu'on peut calculer des convolutions rapidement, revenons à notre problème. % TODO: https://fr.wikipedia.org/wiki/Transformation_de_Fourier_rapide \subsubsection{Corrélation croisée} % TODO: https://en.wikipedia.org/wiki/Cross-correlation Si on définit $\tilde{f}(t)=f(-t)$ La corrélation croisée de \(f\) et \(g\) est \( (g * \tilde{f})(t) = \int_{-\infty}^{\infty} \overline{f(x-t)} g(x) \, dx \) Ou sur une seule période pour une fonction périodique. On va donc faire glisser les fonctions dans le même sens et multiplier les valeurs de la fonction évaluée en ses points ensembles deux à deux puis faire la somme. Quand les fonctions sont similaires et bien alignées, alors, on multiplie les grandes valeurs ensemble et les petites valeurs ensemble et donc on a résultat grand, alors que quand elles ne le sont pas, on multiplie les grandes valeurs avec les petites valeurs et on aura un résultat plus petit. On prend le conjugué pour que si c'est des fonctions complexes, que quand elles s'alignent elles donnent une plus grande valeur au lieu d'une plus petite à cause de la valeur imaginaire. Et donc on ne regardera que la partie réele et à quel point elle est grande. Donc si on a une fonction $f$ et une fonction $g$ où $g$ est dans $f$ mais avec un certain retard $\tau$ alors $(g * \tilde{f})$ sera maximale en $\tau$. Il s'agit donc d'un produit scalaire entre les deux familles de nombres dans le cas discret pour chaque décalage. Le produit scalaire des deux fonctions étant donc maximal quand elles se "ressemblent" le plus. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(x))}; \addlegendentry{$g(x)=\cos(x)$} \addplot[red, thick] {sin(deg(x))}; \addlegendentry{$f(x)=\sin(x)$} \end{axis} \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {pi*sin(deg(x))}; \addlegendentry{$(g * \tilde{f})(t) = \pi \sin(t)$} \end{axis} \end{tikzpicture} \end{center} Ainsi, si on a un enregistrement de référence pour nos différentes ondes $P$ et $S$ alors on peut utiliser la corrélation croisée entre eux et le signal en temps réel pour détecter ces ondes et mesurer avec exactitude le temps d'arrivée de ces ondes. \subsection{Modélisation de la propagation des ondes sismiques} \subsubsection{Calcul du temps de propagation selon iasp91} On va utiliser un modèle de vitesse de propagation des différentes ondes sismiques selon la profondeur et tracer des rayons pour calculer le temps de propagation. On remarque que la courbure qu'on obtient est du à la réfraction. L'onde va emprunter le chemin le plus court, même si cela veut dire aller plus en profondeur d'abord. \noindent\begin{minipage}{.5\textwidth} \centering \includegraphics[width=3cm, trim={0 0 0 0}, clip]{2025-03-26T23:10:40,393866488+01:00-side.png} \captionof{figure}{TauPy} \end{minipage}% \begin{minipage}{.5\textwidth} \centering \includegraphics[width=4cm, trim={0 0 0 0}, clip]{IASP91.png} \captionof{figure}{iasp91} \end{minipage} \subsubsection{Tabulation et interpolation} \begin{center} \includegraphics[width=10cm, trim={12cm 4cm 12cm 8cm}, clip]{2025-03-24T00:28:53,070002973+01:00.png}\\ \captionof{figure}{Visualization des deux tables précalculées} \end{center} \subsection{Multilatération} Une fois qu'on a une table des différentes vitesses de propagation pour les différentes profondeurs et les ondes P et S ainsi que nos mesures, on peut résoudre l'équation avec nos mesures. Sauf qu'on a plusieurs mesures donc on veut une moyenne pour nos mesures pour avoir le meilleur résultat et nos équations ne sont pas linéaires comme vous pouvez le constater sur la diapositive précédente. Ainsi, on va utiliser une méthode pour trouver le minimum sur une fonction d'erreur. \subsubsection{Fonction d'erreur} Voici notre fonction d'erreur: \vspace*{7em} \begin{equation*} \text{E}(\eqnmarkbox[blue]{p1}{depth},\eqnmarkbox[blue]{p2}{lat},\eqnmarkbox[blue]{p3}{lon},\eqnmarkbox[blue]{p4}{epoch},\eqnmarkbox[green]{p5}{obs})=\sum_i \left(\frac{\eqnmarkbox[green]{p6}{(obs_iS-epoch)}-\eqnmarkbox[blue]{p7}{\text{S}(depth, \eqnmark[red]{p8}{\text{greatCircleAngle}(lat,lon,lat_i,lon_i))}}}{\eqnmarkbox[green]{p9}{obs_iS-epoch}}\right)^{\eqnmarkbox[magenta]{p11}{2}} \end{equation*}\vspace*{3em} \begin{equation*} +\eqnmarkbox[pink]{p10}{\left(\frac{(obs_iP-epoch)-\text{P}(depth, \text{greatCircleAngle}(lat,lon,lat_i,lon_i))}{obs_iP-epoch}\right)^2} \end{equation*} \annotate[yshift=-0.3em]{below}{p1}{Profondeur estimée} \annotatetwo[yshift=5.5em]{above}{p2}{p3}{Latitude et longitude estimée} \annotate[yshift=4.5em]{above}{p4}{Date de début du séisme estimée} \annotate[yshift=-0.4em]{below}{p5}{Tableau des observations} \annotate[yshift=2em]{above}{p6}{Temps de propagation avec la date de début du séisme estimée} \annotate[yshift=5em]{above,left}{p7}{Temps de propagation que l'on devrait observer compte tenu du modèle} \annotate[yshift=6.5em]{above,left}{p8}{Calcul de l'angle entre le seismographe et la position estimée du séisme} \annotate[yshift=-1em]{below,left}{p9}{Renormalisation pour éviter de favoriser les temps de propagation longs} \annotate[yshift=3em]{above,left}{p11}{On fait la moyenne quadratique pour avoir l'écart} \annotate[yshift=0em]{below,left}{p10}{Idem mais pour l'onde P} \vspace*{3em} \begin{center} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-31T10:14:11,496825900+02:00.png}\\ \captionof{figure}{Implémentation de la fonction d'erreur} \end{center} \subsubsection{Méthode de Nelder-Mead} La méthode d'optimisation que j'utilise est celle de Nelder-Mead. Elle permet de trouver un minimum local dans une fonction continue avec une recherche directe heuristique. Elle utilise ce qu'on appelle des simplexes, c'est la généralisation des triangles à plus de dimensions, en dimension 1 c'est une ligne, en dimension $3$ ça correspond à un tetraèdre, à chaque fois en dimension $n$, il $n+1$ côtés (il y a $4$ dimensions ici donc il faudra $n+1=5$ cotés au simplexe). L'approche la plus simple consiste à remplacer le plus mauvais point par un point reflétant le centroïde des $n$ points restants. Si ce point est meilleur que le meilleur point actuel, nous pouvons essayer de l'étirer de manière exponentielle le long de cette ligne. En revanche, si ce nouveau point n'est pas beaucoup mieux que la valeur précédente, nous traversons une vallée et nous rétrécissons le simplexe vers un meilleur point. \begin{center} \includegraphics[width=18cm, trim={0 0 0 0}, clip]{neldermead/An-iteration-of-the-Nelder-Mead-method-over-two-dimensional-space-showing-point-p-min.png} \captionof{figure}{Une itération de Nelder-Mead sur un espace de dimension 2} \end{center} \subsection{Magnitude sismique} On utilise l'echelle logarithmique de Richter pour quantifier la magnitude d'un séisme pour laquelle on a la formule suivante: \[ M_\mathrm{L} = \log_{10} \left[ \frac{A}{A_\mathrm{0}(\delta)} \right] \] \noindent où $A$ correspond à l'amplitude maximale mesurée (en m) par le sismographe et $A_\mathrm{0}(\delta)$ un coefficient de correction qui dépend de la distance ($\delta$) à l'epicentre et dont le calcul diffère selon les modèles employés (généralement on utilise une table de corrélation empirique). Pour les grands séismes on préfère utiliser la magnitude de moment qui rend mieux compte de l'énergie du séisme plutôt que l'échelle de Richter (magnitude locale), mais elle est plus compliquée à calculer et nécessite une étude approfondie du séisme (largeur et longueur de la faille, rigidité du milieu, déplacement sur la faille) et n'est donc pas adaptée à un système de prévention en temps réel. Il existe plusieurs formules et tables empiriques pour calculer la magnitude locale d'un séisme. On utilisera la formule empirique de Tsuboi (Université de Tokyo): \[ M_{\mathrm{L}} = \log_{10} A + 1.73 \log_{10} \Delta - 0.83 \] \noindent où \( A \) est l'amplitude maximale en micromètres et \( \Delta \) est la distance en kilomètres. Les différentes valeurs obtenues permettent de quantifier les dégats provoqués par les séismes: \begin{center} \begin{tabular}{|>{\columncolor{white}}l|l|l|p{5cm}|} \hline \rowcolor{gray!30} \textbf{Magnitude} & \textbf{Description} & \textbf{MMI Typique} & \textbf{Effets Moyens du Séisme} \\ \hline \cellcolor{mmi1}1.0 - 1.9 & Micro & I & Micro-séismes, non ressentis. Enregistrés par les sismographes. \\ \hline \cellcolor{mmi1}2.0 - 2.9 & Mineur & I & Légèrement ressenti par certaines personnes. Aucun dommage aux bâtiments. \\ \hline \cellcolor{mmi3}3.0 - 3.9 & Léger & II à III & Souvent ressenti, mais cause rarement des dégâts. Secousses perceptibles d’objets à l’intérieur. \\ \hline \cellcolor{mmi5}4.0 - 4.9 & Faible & IV à V & Secousses intérieures notables et bruits de cliquetis. Légèrement ressenti à l’extérieur. Dégâts minimes possibles. \\ \hline \cellcolor{mmi6}5.0 - 5.9 & Modéré & VI à VII & Peut endommager les bâtiments mal construits ; ressenti par tous. Peu ou pas de dégâts aux bâtiments solides. \\ \hline \cellcolor{mmi7}6.0 - 6.9 & Fort & VII à IX & Dégâts modérés aux structures solides ; dégâts sévères aux structures faibles. Ressenti sur de grandes régions. \\ \hline \cellcolor{mmi8}7.0 - 7.9 & Majeur & VIII ou plus & Dégâts majeurs et effondrements possibles. Dommages concentrés dans un rayon de 250 km. \\ \hline \cellcolor{mmi9}8.0 - 8.9 & Très fort & VIII+ & Destructions majeures à totales. Dommages sur des zones très vastes. Ressenti à très grande distance de l’épicentre. \\ \hline \cellcolor{mmi10}9.0 - 9.9 & Extrême & XII & Destruction quasi-totale, dégâts graves ou effondrement de tous les bâtiments. Modification du relief. \\ \hline \end{tabular} \end{center} \section{Toutes les pièces mises ensembles} \subsection{Fonctionnement} Ainsi, le programme fonctionne comme ceci: Il se connecte à tous les sismographes, commence à streamer leurs mesures. Il fait tourner en continu une corrélation croisée entre les signaux et plusieurs signaux de référence pour les ondes S et P. Si suffisamment de sismographes dans une zone géographique proche détectent des ondes P et S alors une détéction d'événement sismique est enregistrée. Ainsi, le programme met dans une liste les sismographes qui ont détecté cet évenement (donc leur position), et les temps d'arrivée des ondes P et S. À mesure que plus de détecteurs sismiques detectent cet évenement, ils sont rajoutés à la liste. À chaque fois que la liste est mise à jour, il calcule avec la méthode de Nelder-Mead la position (longitude, latitude, profondeur) et temps de séisme qui minimise l'erreur et renvoie celle-ci. Il calcule ensuite à partir de cette position et l'amplitude des mesures la magnitude du séisme est calculée. Toutes ces informations sont renvoyées en temps réel à l'utilisateur et peuvent être utilisées pour alimenter un système d'alerte des populations. Il renvoie une fois toutes les détections finies un compte rendu final avec ces paramètres ainsi qu'une carte qui montre la position du séisme, la fonction d'erreur superposée sur celle-ci selon la position à temps de début de séisme fixé, ainsi que les détecteurs qui ont permis leur détection. Par ailleurs, il y a une phase de calibration qui prend en compte les spéficités de chaque detecteur (sensitivité, bruit de fond, signal continu present à tout moment). \begin{center} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-21-53.png} \captionof{figure}{Les données de stations telles qu'elles sont stockées pour éviter de les recalculer} \end{center} \newpage Voici une carte des stations utilisées: \begin{center} \includegraphics[width=12cm, trim={0 0 0 6cm}, clip]{stations.png} \captionof{figure}{Une carte des stations utilisées} \end{center} Voici ce que renvoie le programme: \begin{center} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1423-40-06.png} \captionof{figure}{Phase de registration des stations} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1423-21-09.png} \captionof{figure}{Arrivée d'une onde} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-08-06.png} \captionof{figure}{Phase de calibration terminée} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1501-18-07.png} \captionof{figure}{Détection d'un séisme de magnitude 3} \end{center} Il y a même une interface web pour voir tout ça en temps réel (notamment les traces des sismographes). \begin{center} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-14-28.png} \captionof{figure}{Interface web} \end{center} \newpage \subsection{Résultats} Voici par exemple un séisme plutôt faible (M6.2) detecté par mon système. On constate que les positions prédites sont en accord avec celles données par l'USGS (United States Geological Survey) ainsi que les magnitudes. \begin{center} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{debugMaps/rs2025fwmrzv.png} \captionof{figure}{rs2025fwmrzv - Mon estimation - Il donne également une magnitude de M5.9} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-25T21:07:44,659295790+01:00-cropped.png} \captionof{figure}{rs2025fwmrzv - Estimation de l'USGS - M6.2 selon eux} \end{center} \newpage \section{Précision et erreur} \subsection{Cas réduit avec propagation rectiligne dans un millieu transparent, homogène, isotrope, linéaire à longues distances et 4 detecteurs} Dans certains cas, comme le GPS (cadre dans lequel on va se placer pour l'analogie), les détecteurs (pour le GPS ils émettent mais c'est symétrique) sont suffisamment loin par rapport à l'erreur pour pouvoir linéariser les équations régissant la détermination de la position (il faudrait encadrer la position de l'appareil avec des surfaces de sphères qui deviennent donc des plans à cette distance). Ainsi dans le cas de 4 détecteurs (pour avoir exactement une solution, quatre inconnues $(x,y,z,t)$ donc quatre équations, $t$ étant la date de l'envoi du ping par tous les satellites), on a une dépendance linéaire entre les temps d'arrivée et la position. Ainsi si on fixe la position du récepteur à l'origine et qu'on désigne la position des quatre satellites par les vecteurs $\vec{a}, \vec{b}, \vec{c}$ et $\vec{d}$. Soit $\vec{M} = (x, y, z)^\top$ notre mesure et $t$ le temps d'émission du signal par tous les satellites (de façon synchrone). \begin{figure}[h] \centering \tdplotsetmaincoords{80}{140} % angle de vue (élévation, azimuth) \begin{tikzpicture}[tdplot_main_coords, scale=2] % Origine (récepteur) \filldraw[black] (0,0,0) circle (1.5pt) node[below] {$\vec{M'}=0$}; % Axes optionnels \draw[->] (0,0,0) -- (1.5,0,0) node[left] {$x$}; \draw[->] (0,0,0) -- (0,1.5,0) node[right] {$y$}; \draw[->] (0,0,0) -- (0,0,1.5) node[above] {$z$}; % Satellite a \coordinate (A) at (2,1,2); \filldraw[blue] (A) circle (1.5pt); \draw[-{Stealth}, blue] (0,0,0) -- (A) node[right] {$\vec{a}$}; \draw[-{Stealth}, red, thick] (0,0,0) -- (0.666667,0.333333,0.666667) node[midway,right] {$\vec{u}_a$}; % Satellite b \coordinate (B) at (-1,2,2); \filldraw[blue] (B) circle (1.5pt); \draw[-{Stealth}, blue] (0,0,0) -- (B) node[right] {$\vec{b}$}; \draw[-{Stealth}, red, thick] (0,0,0) -- (-0.333333,0.666667,0.666667) node[midway,left] {$\vec{u}_b$}; % Satellite c \coordinate (C) at (1,-2,2); \filldraw[blue] (C) circle (1.5pt); \draw[-{Stealth}, blue] (0,0,0) -- (C) node[right] {$\vec{c}$}; \draw[-{Stealth}, red, thick] (0,0,0) -- (0.333333,-0.666667,0.666667) node[midway,right] {$\vec{u}_c$}; % Satellite d \coordinate (D) at (-1,-1,3); \filldraw[blue] (D) circle (1.5pt); \draw[-{Stealth}, blue] (0,0,0) -- (D) node[right] {$\vec{d}$}; \draw[-{Stealth}, red, thick] (0,0,0) -- (-0.3015113445,-0.3015113445,0.9045340337) node[midway,right] {$\vec{u}_d$}; \end{tikzpicture} \caption{Géométrie GPS : satellites et vecteurs unitaires $\vec{u}_i$ associés} \end{figure} \begin{align} \|\vec{M} - \vec{a}\|^2 &= v^2(\tau_a - t)^2 \\ \|\vec{M} - \vec{b}\|^2 &= v^2(\tau_b - t)^2 \\ \|\vec{M} - \vec{c}\|^2 &= v^2(\tau_c - t)^2 \\ \|\vec{M} - \vec{d}\|^2 &= v^2(\tau_d - t)^2 \end{align} On écrit $\vec{M} = \vec{M'} + \vec{\delta M}$, $t = t' + \delta t$, $\tau_a = \tau_a' + \delta\tau_a$, $\tau_b = \tau_b' + \delta\tau_b$, $\tau_c = \tau_c' + \delta\tau_c$, $\tau_d = \tau_d' + \delta\tau_d$. $\vec{M'}$ est la vraie position, donc $\vec{M'} = 0$ vu notre choix de coordonnées, et $t' = 0$ également ; $\tau_a'$, $\tau_b'$, $\tau_c'$ et $\tau_d'$ sont les temps de propagation entre les satellites et l'origine. \newpage On a donc en développant à l'ordre 1 et en identifiant : \begin{align} -2\langle \vec{\delta M}, \vec{a} \rangle &= 2v^2 \tau_a'(\delta\tau_a - \delta t) \\ -2\langle \vec{\delta M}, \vec{b} \rangle &= 2v^2 \tau_b'(\delta\tau_b - \delta t) \\ -2\langle \vec{\delta M}, \vec{c} \rangle &= 2v^2 \tau_c'(\delta\tau_c - \delta t) \\ -2\langle \vec{\delta M}, \vec{d} \rangle &= 2v^2 \tau_d'(\delta\tau_d - \delta t) \end{align} Ce qui donne en réarrangeant : \begin{align} \langle \vec{\delta M}, \vec{a} \rangle + v^2 \tau_a'(\delta\tau_a - \delta t) &= 0 \\ \langle \vec{\delta M}, \vec{b} \rangle + v^2 \tau_b'(\delta\tau_b - \delta t) &= 0 \\ \langle \vec{\delta M}, \vec{c} \rangle + v^2 \tau_c'(\delta\tau_c - \delta t) &= 0 \\ \langle \vec{\delta M}, \vec{d} \rangle + v^2 \tau_d'(\delta\tau_d - \delta t) &= 0 \end{align} On divise chaque ligne par la distance entre le récepteur au repos et le satellite correspondant ($v\tau_i'$). On note $\vec{u}_i$ le vecteur unitaire correspondant. \begin{align} \langle \vec{\delta M}, \vec{u}_a \rangle + v(\delta\tau_a - \delta t) &= 0 \\ \langle \vec{\delta M}, \vec{u}_b \rangle + v(\delta\tau_b - \delta t) &= 0 \\ \langle \vec{\delta M}, \vec{u}_c \rangle + v(\delta\tau_c - \delta t) &= 0 \\ \langle \vec{\delta M}, \vec{u}_d \rangle + v(\delta\tau_d - \delta t) &= 0 \end{align} Si on pose : \[ A = \frac{1}{v} \begin{pmatrix} u_{ax} & u_{ay} & u_{az} & -v \\ u_{bx} & u_{by} & u_{bz} & -v \\ u_{cx} & u_{cy} & u_{cz} & -v \\ u_{dx} & u_{dy} & u_{dz} & -v \end{pmatrix}, \quad X = \begin{pmatrix} \delta x \\ \delta y \\ \delta z \\ \delta t \end{pmatrix}, \quad Y = \begin{pmatrix} \delta\tau_a \\ \delta\tau_b \\ \delta\tau_c \\ \delta\tau_d \end{pmatrix} \] On a $AX = Y$. Ainsi si on choisit $Y$ comme variable aléatoire où chaque coordonnée est indépendante, centrée et de $\mathbb{E}[(\delta\tau_i)^2] = \sigma_i^2$, alors $X = A^{-1}Y$ est une variable aléatoire. Calculons l'espérance de sa norme au carré (on peut également calculer l'espérance du carré d'une coordonée en projettant et en méttant au carré). \[ \|X\|^2 =X^\top X = Y^\top {A^{-1}}^\top A^{-1} Y = Y^\top (AA^\top)^{-1} Y \] $(AA^\top)^{-1}$ étant symétrique positive et inversible (matrice de Gram associée à $A^\top$, i.e.\ aux vecteurs $\vec{u}_i$ complétés d'un $-1$), on peut donc la diagonaliser dans une BON par le théorème spectral. \[ Q = (\vec{v}_1 \mid \vec{v}_2 \mid \vec{v}_3 \mid \vec{v}_4) \] où $\vec{v}_1$, $\vec{v}_2$, $\vec{v}_3$ et $\vec{v}_4$ sont orthonormés. \begin{align} \|X\|^2 = X^\top X &= Y^\top {A^{-1}}^\top A^{-1} Y \\ &= Y^\top Q^\top \operatorname{diag}(\lambda_1, \lambda_2, \lambda_3, \lambda_4)\, Q\, Y \\ &= \sum_{i=1}^{4} \lambda_i \langle \vec{v}_i, Y \rangle^2 \end{align} On passe à l'espérance : \begin{align} \mathbb{E}(\|X\|^2) &= \sum_{i=1}^{4} \lambda_i \,\mathbb{E}\!\left(\langle \vec{v}_i, Y \rangle^2\right) \\ &= \sum_{i=1}^{4} \lambda_i \sum_{j=1}^{4} (v_i)_j^2 \,\sigma_j^2 \end{align} La dernière égalité étant donnée par indépendance. Si on suppose de plus que l'incertitude en temps ($\sigma_j$) est la même pour tous les détecteurs, égale à $\sigma_\tau$, alors cette égalité devient : \[ \mathbb{E}(\|X\|^2) = \operatorname{tr}\!\left((AA^\top)^{-1}\right)\sigma_\tau^2 \] On remarque donc que l'erreur explose quand les valeurs propres de $AA^\top$ sont très petites, i.e.\ quand les vecteurs $\vec{u}_i$ sont proches (dilution de précision géométrique). Si on regarde la forme quadratique induite par $Q(Y) = Y^\top (AA^\top)^{-1} Y$, l'ensemble $Q(Y) = C$ pour une constante $C > 0$ définit un ellipsoïde. Sur cette surface, plus $\|Y\|^2$ est petit, plus le résultat est sensible à l'erreur en temps de mesure selon la direction (dans l'espace des temps de mesure) de $Y$. \begin{figure}[h] \centering \tdplotsetmaincoords{70}{135} \begin{tikzpicture}[tdplot_main_coords, scale=1.5] % j'ai décidé de tourner par la matrice %[ 0.8636654, 0.4567464, -0.2132249; %-0.4446277, 0.8895858, 0.1046102; % 0.2374622, 0.0044575, 0.9713866 ] % Ellipsoide Q(Y)=C projeté en 3D (delta_tau_d fixé) % On visualise dans l'espace (delta_tau_a, delta_tau_b, delta_tau_c) % Axes \draw[-{Stealth}, gray] (0,0,0) -- (4,0,0) node[left] {$\delta\tau_a$}; \draw[-{Stealth}, gray] (0,0,0) -- (0,4,0) node[right] {$\delta\tau_b$}; \draw[-{Stealth}, gray] (0,0,0) -- (0,0,2) node[above] {$\delta\tau_c$}; % Ellipsoide (demi-axes 1/sqrt(lambda_i), ici lambda_1 petite vp -> grande erreur \def\b{1.9} \def\c{1.0} % petit axe -> grande vp -> petite erreur % Meridiens \foreach \ph in {0,30,...,150}{ \draw[blue!30, thin] plot[domain=0:360, samples=60] ({0.8636654*\a*cos(\x)*cos(\ph)+0.4567464*\b*cos(\x)*sin(\ph)-0.2132249*\c*sin(\x)}, {-0.4446277*\a*cos(\x)*cos(\ph)+0.8895858*\b*cos(\x)*sin(\ph)+0.1046102*\c*sin(\x)}, {0.2374622*\a*cos(\x)*cos(\ph)+0.0044575*\b*cos(\x)*sin(\ph)+0.9713866*\c*sin(\x)}); } % Paralleles \foreach \th in {-60,-30,...,60}{ \draw[blue!50, thin] plot[domain=0:360, samples=60] ({0.8636654*\a*cos(\th)*cos(\x)+0.4567464*\b*cos(\th)*sin(\x)-0.2132249*\c*sin(\th)}, {-0.4446277*\a*cos(\th)*cos(\x)+0.8895858*\b*cos(\th)*sin(\x)+0.1046102*\c*sin(\th)}, {0.2374622*\a*cos(\th)*cos(\x)+0.0044575*\b*cos(\th)*sin(\x)+0.9713866*\c*sin(\th)}); } % Cercles de l'ellipsoide sur les plans principaux \draw[blue, thick] plot[domain=0:360, samples=60] ({0.8636654*\a*cos(\x)+0.4567464*\b*sin(\x)}, {-0.4446277*\a*cos(\x)+0.8895858*\b*sin(\x)}, {0.2374622*\a*cos(\x)+0.0044575*\b*sin(\x)}); % ({\a*cos(\x)}, {\b*sin(\x)}, 0); \draw[blue, thick] plot[domain=0:360, samples=60] ({0.8636654*\a*cos(\x)-0.2132249*\c*sin(\x)}, {-0.4446277*\a*cos(\x)+0.1046102*\c*sin(\x)}, {0.2374622*\a*cos(\x)+0.9713866*\c*sin(\x)}); % ({\a*cos(\x)}, 0, {\c*sin(\x)}); \draw[blue, thick] plot[domain=0:360, samples=60] ({0.4567464*\b*cos(\x)-0.2132249*\c*sin(\x)},{0.8895858*\b*cos(\x)+0.1046102*\c*sin(\x)},{0.0044575*\b*cos(\x)+0.9713866*\c*sin(\x)}); % (0, {\b*cos(\x)}, {\c*sin(\x)}); % Vecteurs propres v_i (directions des axes de l'ellipsoide) \draw[-{Stealth}, red, very thick] (0,0,0) -- ({0.8636654*\a}, {-0.4446277*\a}, {0.2374622*\a}) node[left] {$\sqrt{\frac{C}{\lambda_1}}\vec{v}_1$}; \draw[-{Stealth}, orange, very thick] (0,0,0) -- ({0.4567464*\b},{0.8895858*\b} , {0.0044575*\b}) node[right] {$\sqrt{\frac{C}{\lambda_2}}\vec{v}_2$}; \draw[-{Stealth}, green!60!black, very thick] (0,0,0) -- ({-0.2132249*\c},{0.1046102*\c},{0.9713866*\c}) node[above] {$\sqrt{\frac{C}{\lambda_3}}\vec{v}_3$}; % Origine \filldraw[black] (0,0,0) circle (1.5pt) node[below left] {$0$}; % Légende \node at (0,-2.5,1) {($\delta\tau_d = 0$ fixé)}; \end{tikzpicture} \caption{Ellipsoïde $Q(Y) = Y^\top(AA^\top)^{-1}Y = C$ projeté dans l'espace $(\delta\tau_a, \delta\tau_b, \delta\tau_c)$. Les grands axes correspondent aux petites valeurs propres de $AA^\top$ : une faible norme de $Y$ dans ces directions suffit à atteindre la surface, signifiant une grande sensibilité à l'erreur.} \end{figure} \newpage \subsection{Méthode de Monte-Carlo pour les cas plus complexes} Pour les cas plus complexes où la position est déterminée par recherche de minimum ou encore que les distances sont suffisamment faibles pour ne pas pouvoir linéariser les équations, les calculs deviennent rapidement compliqués. Ainsi on peut opter pour une méthode de Monte-Carlo. On estime l'espérance de la norme carrée du vecteur position selon nos paramètres en faisant pleins de simulations en faisant varier aléatoirement les valeurs des mesures selon l'erreur qu'on veut. J'ai écrit un programme pour le faire dans le cas de mesures faites à la main avec la propagation d'un onde sonore (clappement de main) avec 4 microphones synchronisés et placés dans une petite salle (hypothèse de linéarisation non vérifiée). D'une part je calcule la position par la même méthode qu'avec les séismes et d'autre part je calcule l'erreur par Monte-Carlo selon la disposition des microphones avec la même erreur de mesure pour chacun (distribution uniforme par ailleurs). On voit bien l'effet de la dilution de précision géométrique. Voici une configuration écrasée, il y a une grande dilution de précision géométrique dans le plan engendré par les detecteurs. \begin{figure}[h] \centering \begin{minipage}{0.49\textwidth} \centering \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/1.png} \caption{Rouge: détecteurs, Vert: position mesurée, Bleu: référence (configuration écrasée)} \end{minipage} \hfill \begin{minipage}{0.49\textwidth} \centering \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/2.png} \caption{Erreur calculée par Monte-Carlo selon la position dans l'espace à erreur de mesure égale par détecteur (configuration écrasée)} \end{minipage} \end{figure} \newpage Lorsque la configuration donne des vecteurs moins colinéaires on a bien une meilleure précision. \begin{figure}[h] \centering \begin{minipage}{0.49\textwidth} \centering \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/3.png} \caption{Rouge: détecteurs, Vert: position mesurée, Bleu: référence (configuration favorable)} \end{minipage} \hfill \begin{minipage}{0.49\textwidth} \centering \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/4.png} \caption{Erreur calculée par Monte-Carlo selon la position dans l'espace à erreur de mesure égale par détecteur (configuration favorable)} \end{minipage} \end{figure} \newpage \section{Bonus} De plus ces expérimentations de traitements de signaux m'ont amené à construire un oscilloscope logiciel avec auto-correlation (pour détecter les périodes) et corrélation avec une onde sinusoidale de la période détectée pour toujours avoir la même phase au centre (\url{https://github.com/make-42/xyosc}). \begin{center} \includegraphics[width=5cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-42-57.png} \captionof{figure}{Mode canal unique avec auto-correlation} \includegraphics[width=5cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-41-26.png} \captionof{figure}{Mode XY} \includegraphics[width=5cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1512-44-34.png} \captionof{figure}{Mode spectrogramme avec détection de note} \end{center} \newpage \section{Bibliographie} \begin{thebibliography}{9} \bibitem{convolution_video} 3Blue1Brown, \textit{But what is a convolution?}, YouTube video, 2022.\\ \url{https://www.youtube.com/watch?v=KuXjwB4LzSA} \bibitem{fft_algorithm_video} Michael Pound, \textit{The Fast Fourier Transform Algorithm}, YouTube video, 2023.\\ \url{https://www.youtube.com/watch?v=toj_IoCQE-4} \bibitem{fft_matrix_factorizations} Charles Van Loan, \textit{The FFT Via Matrix Factorizations}, Lecture notes, 2010.\\ \url{https://www.cs.cornell.edu/~bindel/class/cs5220-s10/slides/FFT.pdf} \bibitem{eq_sources} R.J. Mitchell, \textit{Earthquake Sources}, Lecture notes, Western Washington University.\\ \url{https://www.geol.wwu.edu/rjmitch/L4_EQsources.pdf} \bibitem{seismic_waves} University of Hawaii, \textit{Compare, Contrast, and Connect: Seismic Waves and Determining Earth's Structure}.\\ \url{https://manoa.hawaii.edu/exploringourfluidearth/physical/ocean-floor/layers-earth/compare-contrast-connect-seismic-waves-and-determining-earth-s-structure} \bibitem{gfz-1-d-earth-models} Peter Bormann, \textit{Global 1-D Earth models (IASP91 tables)}, GFZ German Research Centre for Geosciences.\\ \url{https://gfzpublic.gfz-potsdam.de/rest/items/item_4031/component/file_4032/content} \bibitem{earthquake_data_centers} Yacine Boussoufa, \textit{Earthquake Data Centers}, GitHub repository.\\ \url{https://github.com/YacineBoussoufa/EarthquakeDataCenters} \bibitem{cfcs_lecture} School of Informatics, University of Edinburgh, \textit{Computational Foundations of Cognitive Science: Lecture 15 (Convolutions and Kernels)}.\\ \url{https://www.inf.ed.ac.uk/teaching/courses/cfcs1/lectures/cfcs_l15.pdf} \bibitem{penn_state_anova} Penn State Eberly College of Science, \textit{STAT 510: Lesson 8.2 - Cross Correlation Functions and Lagged Regressions}, Online course material.\\ \url{https://online.stat.psu.edu/stat510/lesson/8/8.2} \bibitem{nelder_mead} Jason Cantarella, \textit{Nelder-Mead Method}, Lecture notes.\\ \url{https://jasoncantarella.com/downloads/NelderMeadProof.pdf} \bibitem{broadband_magnitude} Tatsuhiko Hara, \textit{Determination of Broadband Moment Magnitude}, IISEE/BRI.\\ \url{https://iisee.kenken.go.jp/lna/download.php?f=2011082925678c01.pdf&n=T0-100-2007_Mwp-2-new.pdf&cid=T0-100-2007} \bibitem{obspy} Moritz Beyreuther, Robert Barsch, Lion Krischer, Tobias Megies, Yannik Behr and Joachim Wassermann,\\ \textit{ObsPy: A Python Toolbox for Seismology},\\ Seismological Research Letters, vol. 81, no. 3, pp. 530--533, 2010.\\ doi:10.1785/gssrl.81.3.530 \end{thebibliography} \end{document}