1\documentclass[a4paper, 11pt]{article}
2\usepackage[a4paper, total={18cm, 24cm}]{geometry}
3\usepackage[T1]{fontenc}
4\usepackage[utf8]{inputenc}
5\usepackage[french]{babel}
6\usepackage{lmodern}
7\usepackage{amsmath}
8\usepackage{amsfonts}
9\usepackage{amssymb}
10\usepackage{amsthm}
11\usepackage{graphicx}
12\usepackage{color}
13\usepackage{xcolor}
14\usepackage{url}
15\usepackage{theorem}
16\usepackage{textcomp}
17\usepackage{listings}
18\usepackage{hyperref}
19%\usepackage{glossaries}
20\usepackage{parskip}
21\usepackage{amsmath, amsthm, amssymb}
22\usepackage{stmaryrd}
23\usepackage{graphicx}
24\usepackage{subfig}
25\usepackage{color}
26\usepackage{longtable}
27\usepackage{pgfplots}
28\usepackage{nicematrix}
29\usepackage[table]{xcolor}
30\usepackage{tikz}
31\usepackage{tikz-3dplot}
32\usetikzlibrary{arrows.meta, 3d}
33\usepackage{assets/texpackages/annotate-equations}
34\newenvironment{graytext}{\color{gray}}{\ignorespacesafterend}
35\graphicspath{ {./assets/} }
36
37% Define custom colors
38\definecolor{mmi1}{HTML}{FFFFFF}
39\definecolor{mmi2}{HTML}{BFCCFF}
40\definecolor{mmi3}{HTML}{A0E6FF}
41\definecolor{mmi4}{HTML}{80FFFF}
42\definecolor{mmi5}{HTML}{7AFF93}
43\definecolor{mmi6}{HTML}{FFFF00}
44\definecolor{mmi7}{HTML}{FFC800}
45\definecolor{mmi8}{HTML}{FF9100}
46\definecolor{mmi9}{HTML}{FF0000}
47\definecolor{mmi10}{HTML}{C80000}
48\definecolor{mmi11}{HTML}{A40000}
49\definecolor{mmi12}{HTML}{800000}
50
51\begin{document}
52\title{Traitement de signaux pour la détection de séismes et leur multilatération}
53\author{Dalibard Louis}
54\date{\today}
55
56\maketitle
57\tableofcontents
58
59\newpage
60\section{Séismes}
61\subsection{Introduction}
62Un 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.\\
63\includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{eq-ed-fault-labeled.png}
64\subsection{Ondes P et S}
65Les 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.
66
67La 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).
68\noindent\begin{minipage}{.5\textwidth}
69\centering
70\includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_compression_impulsion_1d_30_petit.png} \captionof{figure}{Ondes P}
71\end{minipage}%
72\begin{minipage}{.5\textwidth}
73\centering
74\includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_cisaillement_impulsion_1d_30_petit.png}
75 \captionof{figure}{Ondes S}
76\end{minipage}
77\subsection{Différence de vitesse de propagation}
78Les 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.
79
80Les 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.
81\begin{center}
82 \includegraphics[width=7cm, trim={0 12.5mm 0cm 2cm}, clip]{2025-03-26T23:29:53,572838335+01:00.png}\\
83 \captionof{figure}{6 km/s (ondes P) vs 4 km/s (ondes S)}
84 \end{center}
85
86On peut donc exploiter les temps d'arrivée de ces différentes ondes pour localiser un séisme.
87
88\section{Théorie}
89\subsection{Principe}
90Différentes étapes:
91 \begin{enumerate}
92 \item Acquisition de données en temps réel (SeedLink)
93 \item Reconnaissance d'un séisme et mesure automatique des temps
94 \item Calcul de la position et de la magnitude
95 \begin{enumerate}
96 \item Modélisation de la propagation des ondes sismiques
97 \item Méthode numérique d'optimisation de fonction à plusieurs variables pour la multilatération
98 \item Calcul de la magnitude
99 \end{enumerate}
100 \end{enumerate}
101
102\subsection{DSP}
103\subsubsection{Acquisition des données}
104On utilise le protocole Seedlink pour récupérer les données des sismographes en temps réel.
105
106\renewcommand{\arraystretch}{1.2}
107
108
109\begin{longtable}{|l|l|}
110 \hline
111 \textbf{Name} & \textbf{Host} \\
112 \hline
113 AusPass & auspass.edu.au \\
114 BGR & eida.bgr.de \\
115 CISMID & www.cismid.uni.edu.pe \\
116 ENS & ephesite.ens.fr \\
117 GEOFON, GFZ & geofon.gfz-potsdam.de \\
118 GEONET & link.geonet.org.nz \\
119 Geoscience Australia & seis-pub.ga.gov.au \\
120 GSRAS (?) & 89.22.182.133 \\
121 Helsinki & finseis.seismo.helsinki.fi \\
122 Haiti & ayiti.unice.fr \\
123 ICGC & ws.icgc.cat \\
124 IDA Project & rtserve.ida.ucsd.edu \\
125 IFZ & data.ifz.ru \\
126 IPGP & rtserver.ipgp.fr \\
127 IRIS DMC & rtserve.iris.washington.edu \\
128 IRIS Jamaseis & jamaseis.iris.edu \\
129 ISNET - UNINA & 185.15.171.86 \\
130 LMU & erde.geophysik.uni-muenchen.de \\
131 NIGGG & 195.96.231.100 \\
132 NRCAN & earthquakescanada.nrcan.gc.ca \\
133 OBSEBRE & obsebre.es \\
134 OGS & nam.ogs.it \\
135 Oklahoma University & rtserve.ou.edu \\
136 ORFEUS & eida.orfeus-eu.org \\
137 PLSN (IGF Poland) & hudson.igf.edu.pl \\
138 Red Sìsmica de Puerto Rico & 161.35.236.45 \\
139 Red Sìsmica Baru & helis.redsismicabaru.com \\
140 RESIF & rtserve.resif.fr \\
141 SANET & 147.213.113.73 \\
142 RSIS & rsis1.on.br \\
143 SCSN-USC (South Carolina Seismic Network) & eeyore.seis.sc.edu:6382 \\
144 Seisme IRD & rtserve.ird.nc \\
145 Staneo & vibrato.staneo.fr \\
146 SNAC NOA & snac.gein.noa.gr \\
147 TexNet & rtserve.beg.utexas.edu \\
148 Thai Meteorological Department & 119.46.126.38 \\
149 UFRN (Universidade Federal do Rio Grande do Norte) & sislink.geofisica.ufrn.br \\
150 Unical Universita Della Calabria & www.sismocal.org \\
151 UNITS Università degli studi di Trieste & rtweb.units.it \\
152 UNIV-AG Université des Antilles & seedsrv0.ovmp.martinique.univ-ag.fr \\
153 Universidade de Évora & clv-cge.uevora.pt \\
154 Universidad de Colima & 148.213.24.15 \\
155 UPR & worm.uprm.edu \\
156 USGS & cwbpub.cr.usgs.gov \\
157 USP-IAG & seisrequest.iag.usp.br \\
158 \hline
159\end{longtable}
160\subsubsection{Extraction des temps d'arrivée}
161\begin{center}
162 \includegraphics[width=15cm, trim={0 0 0 0}, clip]{R0ED0.EHZ-1743142835749.9944-cropped.png}\\
163 \captionof{figure}{Exemple d'un enregistrement de sismographe}
164\end{center}
165On 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 ?
166Pour cela on va faire un détour par des maths.
167
168\subsubsection{Convolution}
169Si on a deux fonctions, quelles opérations peut-on faire pour avoir une nouvelle fonction?
170
171La plus simple semble être les additionner.
172
173\begin{center}
174\begin{tikzpicture}
175 \begin{axis}[
176 axis lines = middle,
177 xlabel = {$x$},
178 ylabel = {$y$},
179 samples=300,
180 domain=-6:6,
181 legend pos=north east,
182 width=0.8\textwidth,
183 height=.4\textheight
184 ]
185 % Define the functions
186 \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$};
187 \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$};
188 \addplot[green, thick] {sin(deg(x)) + cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) + g(x)$};
189
190 % Add legend
191 \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) + g(x)$}
192 \end{axis}
193\end{tikzpicture}
194\end{center}
195
196
197
198Une autre façon serait de les multiplier.
199
200\begin{center}
201\begin{tikzpicture}
202 \begin{axis}[
203 axis lines = middle,
204 xlabel = {$x$},
205 ylabel = {$y$},
206 samples=300,
207 domain=-6:6,
208 legend pos=north east,
209 width=0.8\textwidth,
210 height=.4\textheight
211 ]
212 % Define the functions
213 \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$};
214 \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$};
215 \addplot[green, thick] {sin(deg(x)) * cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) \cdot g(x)$};
216
217 % Add legend
218 \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) \cdot g(x)$}
219 \end{axis}
220\end{tikzpicture}
221\end{center}
222
223Mais il existe aussi une autre opération, la convolution.
224
225\begin{center}
226\begin{tikzpicture}
227 \begin{axis}[
228 axis lines = middle,
229 xlabel = {$x$},
230 ylabel = {$y$},
231 samples=300,
232 domain=-6:6,
233 legend pos=north east,
234 width=0.8\textwidth,
235 height=.4\textheight
236 ]
237 % Define the functions
238 \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$};
239 \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$};
240 \addplot[green, thick] {0} node[right] {\footnotesize $h(x) = (f*g)(x)$};
241
242 % Add legend
243 \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = (f*g)(x) = 0$}
244 \end{axis}
245\end{tikzpicture}
246\end{center}
247
248Qui 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).
249
250On définit formellement la convolution comme un produit sur l'espace des fonctions intégrables.
251
252Soit $f$ et $g$ deux fonctions intégrables sur $\mathbb{R}$.
253
254Pour tout $t \in \mathbb{R}$,
255\[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \]
256
257Dans le cas discret on utilisera une somme.
258
259Soit $f$ et $g$ deux fonctions de $\mathbb{Z}$ dans $\mathbb{C}$.
260
261Pour tout $n \in \mathbb{Z}$,
262\[ (f*g)[n]=\sum_{m=-\infty}^{+\infty} f[m] \cdot g[n-m] \]
263
264et dans le cas de fonction périodiques, il suffit de faire l'intégrale sur une période.
265Pour tout $t \in \mathbb{R}$,
266\[ (f*g)(t)=\int_{0}^{T} f(x) \cdot g(t-x) dx \]
267
268(on peut également généraliser cela à des fonctions discretès périodiques)
269
270\subsubsection{Propriétés algébriques de la convolution}
271\begin{itemize}
272\item Commutatif
273
274On remarquera que si \[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \]
275
276Et on fait le changement de variable $u=t-x$
277
278On 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)\]
279
280\item Associatif
281\[ (f*g)*h= f*(g*h) \]
282C'est une conséquence du théorème de Fubini.
283
284\item Distributif
285\[ f*(g+h)= f*g+f*h \]
286Par linéarité de l'intégrale.
287
288\end{itemize}
289
290L'espace des fonctions intégrables muni de $*$ forme un demi-groupe commutatif (car pas d'élement neutre).
291
292Pourquoi une telle construction ?
293
294Déjà à quoi ça ressemble?
295
296Si 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.
297
298\begin{center}
299 \begin{tikzpicture}[scale=0.5]
300 % Define Dice Faces
301 \def\diceA{1,2,3,4,5,6} % Top Die (Static)
302 \def\diceB{6,5,4,3,2,1} % Bottom Die (Sliding)
303
304 % Define step spacing
305 \def\stepSpace{3}
306
307 % Loop through each step
308 \foreach \step in {-5,-4,-3,-2,-1,0,1,2,3,4,5} {
309 \pgfmathsetmacro{\sommedonnant}{int(7 + \step)}
310 % Label step number
311 \node[anchor=east] at (-0.5, -\step*3) {Somme donnant \sommedonnant:};
312
313 % Top Row - Fixed Dice
314 \foreach \x [count=\i] in \diceA {
315 \node[draw, minimum size=0.5cm] at (\i, -\step*3) {\x};
316 }
317
318 % Bottom Row - Sliding Dice
319 \foreach \x [count=\i] in \diceB {
320 \pgfmathsetmacro{\sumval}{int(\i + \step)}
321 \ifnum \sumval > 0 \ifnum \sumval < 7 \node[draw=red, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x};
322 \else
323 \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x};
324 \fi
325 \else
326 \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x};
327 \fi
328 }
329 }
330 \end{tikzpicture}
331\end{center}
332\begin{center}
333 \begin{tikzpicture}
334 \begin{axis}[
335 xlabel={Somme},
336 ylabel={Fréquence},
337 ymin=0,
338 ymax=7,
339 xmin=1.5,
340 xmax=12.5,
341 xtick={2,3,4,5,6,7,8,9,10,11,12},
342 ytick={1,2,3,4,5,6,7},
343 area style,
344 width=12cm,
345 height=8cm,
346 major grid style={line width=.2pt,draw=gray!50},
347 grid=major,
348 bar width=0.8,
349 title={Distribution des sommes},
350 nodes near coords,
351 ]
352 \addplot[fill=blue!40, draw=blue!80, opacity=0.8] coordinates {
353 (2,1)
354 (3,2)
355 (4,3)
356 (5,4)
357 (6,5)
358 (7,6)
359 (8,5)
360 (9,4)
361 (10,3)
362 (11,2)
363 (12,1)
364 };
365 \end{axis}
366\end{tikzpicture}
367\end{center}
368
369On 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.
370
371On vient de faire la convolution entre cette fonction et elle même.
372\begin{center}
373 \begin{tikzpicture}
374 \begin{axis}[
375 xlabel={$x$},
376 ylabel={$f(x)$},
377 xmin=-1,
378 xmax=8,
379 ymin=-0.2,
380 ymax=1.4,
381 xtick={-1,0,1,2,3,4,5,6,7,8},
382 ytick={0,1},
383 width=10cm,
384 height=6cm,
385 grid=major,
386 title={Fonction nulle puis constante égale à $1$ puis nulle},
387 samples=100,
388 domain=-1:8,
389 axis lines=middle,
390 clip=false,
391 ]
392 % Plot the function
393 \addplot[blue, thick, const plot] coordinates {
394 (-1, 0)
395 (0.999, 0)
396 (1, 1)
397 (6, 1)
398 (6.001, 0)
399 (8, 0)
400 };
401
402 % Add markers to explicitly show the endpoints
403 \addplot[only marks, mark=*, mark options={fill=red}] coordinates {
404 (1, 1)
405 (6, 1)
406 };
407
408 \end{axis}
409\end{tikzpicture}
410\end{center}
411\begin{center}
412\includegraphics[width=14cm, trim={0 0 0 0}, clip]{2025-04-01T23:42:09,818697014+02:00.png}
413\captionof{figure}{numpy confirme ce résultat}
414\end{center}
415
416On remarque que si on le fait dans le cas discret, on a des valeurs en plus sur les bords que l'on peut ignorer.
417
418Faire 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.
419
420\begin{center}
421 \begin{tikzpicture}
422 \begin{axis}[
423 axis lines = middle,
424 xlabel = {$x$},
425 ylabel = {$y$},
426 samples=300,
427 domain=-6:6,
428 legend pos=north east,
429 width=0.8\textwidth,
430 height=.4\textheight
431 ]
432 \addplot[blue, thick] {cos(deg(2*x))};
433 \addlegendentry{$\cos(2x)$}
434 \addplot[red, thick] {sin(deg(x))};
435 \addlegendentry{$\sin(x)$}
436 \end{axis}
437 \end{tikzpicture}
438 \end{center}
439
440 \begin{center}
441 \begin{tikzpicture}
442 \begin{axis}[
443 axis lines = middle,
444 xlabel = {$x$},
445 ylabel = {$y$},
446 samples=300,
447 domain=-6:6,
448 legend pos=north east,
449 width=0.8\textwidth,
450 height=.4\textheight
451 ]
452 \addplot[blue, thick] {cos(deg(2*x))};
453 \addlegendentry{$\cos(2x)$}
454 \addplot[red, thick] {sin(deg(-x))};
455 \addlegendentry{$\sin(-x)$}
456 \end{axis}
457 \end{tikzpicture}
458 \end{center}
459
460 \begin{center}
461 \begin{tikzpicture}
462 \begin{axis}[
463 axis lines = middle,
464 xlabel = {$x$},
465 ylabel = {$y$},
466 samples=300,
467 domain=-6:6,
468 legend pos=north east,
469 width=0.8\textwidth,
470 height=.4\textheight
471 ]
472 \addplot[blue, thick] {cos(deg(2*x))*sin(deg(-x))};
473 \addlegendentry{$\cos(2x)\sin(-x)$}
474
475 % Highlight positive areas in red
476 \addplot [fill=red, opacity=0.3, domain=-2*pi:2*pi]
477 {max(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle;
478
479 % Highlight negative areas in blue
480 \addplot [fill=blue, opacity=0.3, domain=-2*pi:2*pi]
481 {min(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle;
482
483 \end{axis}
484 \end{tikzpicture}
485 \end{center}
486
487 On voit que c'est nul dans ce cas.
488
489\subsubsection{Moyennage}
490
491Si on applique la fonction $\frac{f}{6}$ à une fonction $g$ quelconque on observe une sorte de moyennage qui se fait entre 6 valeurs proches.
492\begin{center}
493 \begin{tikzpicture}
494 \begin{axis}[
495 xlabel={$x$},
496 ylabel={$y$},
497 xmin=-0.5,
498 xmax=23.5,
499 ymin=-0.5,
500 ymax=9,
501 xtick={0,2,4,6,8,10,12,14,16,18,20,22},
502 width=12cm,
503 height=8cm,
504 grid=major,
505 title={Moyennage par convolution},
506 legend pos=north west,
507 ]
508 % Plot the first dataset with connected lines
509 \addplot[blue, thick] coordinates {
510 (0,0) (1,0) (2,0) (3,0) (4,1) (5,2) (6,3) (7,4) (8,5) (9,6) (10,7) (11,8)
511 (12,0) (13,4) (14,2) (15,3) (16,8) (17,1) (18,2) (19,4) (20,8) (21,0) (22,0) (23,0)
512 };
513
514 % Plot the second dataset with connected lines
515 \addplot[red, thick] coordinates {
516 (0,0)
517 (1,0.16666667)
518 (2,0.5)
519 (3,1)
520 (4,1.66666667)
521 (5,2.5)
522 (6,3.5)
523 (7,4.5)
524 (8,5.5)
525 (9,5)
526 (10,5)
527 (11,4.5)
528 (12,4)
529 (13,4.16666667)
530 (14,3)
531 (15,3.33333333)
532 (16,3.33333333)
533 (17,4.33333333)
534 (18,3.83333333)
535 (19,2.5)
536 (20,2.33333333)
537 (21,2)
538 (22,1.33333333)
539 (23,0)
540 };
541
542 \legend{$g$, $g*\frac{f}{6}$}
543
544 \end{axis}
545 \end{tikzpicture}
546\end{center}
547
548Ce 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.
549
550\subsubsection{Transformée de Fourier discrète}
551Si 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}]$,
552
553On 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}$
554
555
556On 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$
557
558$(a*b)[p]$ est le coefficient du terme de degré $p$ dans le produit
559
560$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$
561
562Bon, de premier abord, ça n'a pas l'air de beaucoup nous aider, le produit de polynômes est $\mathcal{O}(n^2)$.
563
564Mais 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.
565
566Et 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.
567
568Le problème c'est que cela semble prendre plus de temps, un pivot de Gauss est en $\mathcal{O}(n^3)$
569
570Mais si on choisit d'évaluer nos valeurs en les racines de l'unité, on a un système plus simple à résoudre.
571
572On va pouvoir faire les opérations en $\mathcal{O}(n \cdot log(n))$
573
574Les 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.
575\subsubsection{Radix-2 decimation-in-time (DIT) - Factorisation Cooley-Tukey}
576Evaluer 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:
577
578\[
579R=\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} =
580\begin{bmatrix}
5811 & 1 & 1 & \cdots & 1 \\
5821 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\
5831 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\
584\vdots & \vdots & \vdots & \ddots & \vdots \\
5851 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}
586\end{bmatrix}
587\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}
588\]
589
590On reconnaît la matrice de Vandermonde transposée associée à l'évaluation en les \( \omega_n^k \)
591
592On 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$.
593
594On note la matrice de Vandermonde transposée, qui permet de calculer le DFT pour tous nos coefficients,
595
596$F_{2^p}=\begin{bmatrix}
5971 & 1 & 1 & \cdots & 1 \\
5981 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\
5991 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\
600\vdots & \vdots & \vdots & \ddots & \vdots \\
6011 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}
602\end{bmatrix}$
603
604On peut permuter les lignes, à la fin, ce qu'on a c'est que, du fait des propriétés des racines de l'unité,
605
606si on note, $D_{2^{p-1}} =
607\begin{bmatrix}
6081 & 0 & 0 & \cdots & 0 \\
6090 & \omega_{2^{p-1}} & 0 & \cdots & 0 \\
6100 & 0 & \omega_{2^{p-1}}^2 & \cdots & 0 \\
611\vdots & \vdots & \vdots & \ddots & \vdots \\
6120 & 0 & 0 & \cdots & \omega_{2^{p-1}}^{2^{p-1}-1}
613\end{bmatrix}$
614
615est:
616
617\begin{equation*}
618R=F_{2^p}\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}=
619\begin{pNiceArray}{cc|cc}
620 \eqnmarkbox[blue]{mcoeffdiag1}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag2}{D_{2^{p-1}}} \\
621 \hline
622 \eqnmarkbox[blue]{mcoeffdiag3}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag4}{-D_{2^{p-1}}}
623\end{pNiceArray}
624\begin{pNiceArray}{cc|cc}
625 \eqnmarkbox[cyan]{mrecurrence1}{F_{2^{p-1}}} && \mathbf{0} \\
626 \hline
627 \mathbf{0} && \eqnmarkbox[cyan]{mrecurrence2}{F_{2^{p-1}}}
628\end{pNiceArray}
629\begin{bmatrix}
630 \eqnmarkbox[green]{mcoeffpairs1}{a_0} \\
631 \eqnmarkbox[green]{mcoeffpairs2}{a_2} \\
632 \vdots \\
633 \eqnmarkbox[green]{mcoeffpairs3}{a_{2^{p-1}-2}} \\
634 \eqnmarkbox[green]{mcoeffpairs4}{a_{2^{p-1}}} \\
635 \eqnmarkbox[magenta]{mcoeffimpairs1}{a_{1}} \\
636 \eqnmarkbox[magenta]{mcoeffimpairs2}{a_{3}} \\
637 \vdots \\
638 \eqnmarkbox[magenta]{mcoeffimpairs3}{a_{n-3}} \\
639 \eqnmarkbox[magenta]{mcoeffimpairs4}{a_{n-1}}
640\end{bmatrix}
641\end{equation*}
642\annotatetwo[yshift=1.5em]{above}{mcoeffdiag1}{mcoeffdiag2}{Blocs diagonaux, de l'ordre de $\mathcal{O}(n)$ opérations}
643\annotate[yshift=0.5em]{above}{mcoeffpairs1}{Coefficients d'indice pair}
644\annotate[yshift=0em]{below}{mcoeffimpairs4}{Coefficients d'indice impair}
645\annotate[yshift=0em]{below,left}{mrecurrence2}{Récurrence}
646
647
648On évalue la compléxité de l'algorithme.\\
649On note $u_p$ sa complexité en fonction de $p$ et $C_n$ sa complexité en fonction de $n$.
650
651
652\begin{equation*}
653u_{p+1}=\eqnmarkbox[magenta]{cdiagmult}{A \cdot 2^{p+1}}+\eqnmarkbox[red]{crec}{2u_{p}}
654\end{equation*}
655\annotate[yshift=0.5em]{above,left}{cdiagmult}{Compléxité du produit sur la diagonale}
656\annotate[yshift=0em]{below,right}{crec}{Traitement des coefficients par récurrence}
657
658On factorise par la solution homogène,
659$\frac{u_{p+1}}{2^{p+1}}=A+\frac{u_{p}}{2^p}$
660
661$\frac{u_{p}}{2^{p}}=u_{0}+A\cdot p$
662
663$u_{p}=u_{0}\cdot2^{p}+A\cdot p \cdot 2^{p}$
664
665Or $p=\log_2{n}$
666
667Donc $C_n = u_{\log_2{n}} = u_{0} \cdot n+A \cdot \log_2{n} \cdot n = \mathcal{O}(n \log n)$
668
669\subsubsection{IFFT}
670On peut montrer que:
671$F_{2^p}^{-1}=\frac{1}{2^p}\begin{bmatrix}
6721 & 1 & 1 & \cdots & 1 \\
6731 & \overline{\omega_n}^1 & \overline{\omega_n}^2 & \cdots & \overline{\omega_n}^{n-1} \\
6741 & \overline{\omega_n}^2 & \overline{\omega_n}^4 & \cdots & \overline{\omega_n}^{2(n-1)} \\
675\vdots & \vdots & \vdots & \ddots & \vdots \\
6761 & \overline{\omega_n}^{n-1} & \overline{\omega_n}^{2(n-1)} & \cdots & \overline{\omega_n}^{(n-1)(n-1)}
677\end{bmatrix}$
678
679$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}$
680
681Le 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).
682
683$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}}$
684
685On 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.
686
687\subsubsection{Passage du domaine temporel au domaine fréquentiel}
688Sans 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.
689
690Maintenant qu'on peut calculer des convolutions rapidement, revenons à notre problème.
691
692% TODO: https://fr.wikipedia.org/wiki/Transformation_de_Fourier_rapide
693\subsubsection{Corrélation croisée}
694% TODO: https://en.wikipedia.org/wiki/Cross-correlation
695Si on définit $\tilde{f}(t)=f(-t)$
696
697La corrélation croisée de \(f\) et \(g\) est \( (g * \tilde{f})(t) = \int_{-\infty}^{\infty} \overline{f(x-t)} g(x) \, dx \)
698
699Ou sur une seule période pour une fonction périodique.
700
701On 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.
702
703Quand 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.
704
705On 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.
706
707Et donc on ne regardera que la partie réele et à quel point elle est grande.
708
709Donc 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$.
710
711Il 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.
712
713\begin{center}
714 \begin{tikzpicture}
715 \begin{axis}[
716 axis lines = middle,
717 xlabel = {$x$},
718 ylabel = {$y$},
719 samples=100,
720 domain=-6:6,
721 legend pos=north east,
722 width=0.8\textwidth,
723 height=.4\textheight
724 ]
725 \addplot[blue, thick] {cos(deg(x))};
726 \addlegendentry{$g(x)=\cos(x)$}
727 \addplot[red, thick] {sin(deg(x))};
728 \addlegendentry{$f(x)=\sin(x)$}
729 \end{axis}
730 \end{tikzpicture}
731 \end{center}
732
733 \begin{center}
734 \begin{tikzpicture}
735 \begin{axis}[
736 axis lines = middle,
737 xlabel = {$x$},
738 ylabel = {$y$},
739 samples=100,
740 domain=-6:6,
741 legend pos=north east,
742 width=0.8\textwidth,
743 height=.4\textheight
744 ]
745 \addplot[blue, thick] {pi*sin(deg(x))};
746 \addlegendentry{$(g * \tilde{f})(t) = \pi \sin(t)$}
747 \end{axis}
748 \end{tikzpicture}
749 \end{center}
750
751Ainsi, 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.
752
753\subsection{Modélisation de la propagation des ondes sismiques}
754\subsubsection{Calcul du temps de propagation selon iasp91}
755On 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.
756
757\noindent\begin{minipage}{.5\textwidth}
758\centering
759\includegraphics[width=3cm, trim={0 0 0 0}, clip]{2025-03-26T23:10:40,393866488+01:00-side.png}
760 \captionof{figure}{TauPy}
761\end{minipage}%
762\begin{minipage}{.5\textwidth}
763\centering
764 \includegraphics[width=4cm, trim={0 0 0 0}, clip]{IASP91.png}
765 \captionof{figure}{iasp91}
766\end{minipage}
767
768
769\subsubsection{Tabulation et interpolation}
770\begin{center}
771 \includegraphics[width=10cm, trim={12cm 4cm 12cm 8cm}, clip]{2025-03-24T00:28:53,070002973+01:00.png}\\
772 \captionof{figure}{Visualization des deux tables précalculées}
773\end{center}
774\subsection{Multilatération}
775Une 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.
776
777Sauf 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.
778
779Ainsi, on va utiliser une méthode pour trouver le minimum sur une fonction d'erreur.
780
781\subsubsection{Fonction d'erreur}
782Voici notre fonction d'erreur:
783\vspace*{7em}
784
785\begin{equation*}
786 \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}}
787 \end{equation*}\vspace*{3em}
788 \begin{equation*}
789 +\eqnmarkbox[pink]{p10}{\left(\frac{(obs_iP-epoch)-\text{P}(depth, \text{greatCircleAngle}(lat,lon,lat_i,lon_i))}{obs_iP-epoch}\right)^2}
790\end{equation*}
791\annotate[yshift=-0.3em]{below}{p1}{Profondeur estimée}
792\annotatetwo[yshift=5.5em]{above}{p2}{p3}{Latitude et longitude estimée}
793\annotate[yshift=4.5em]{above}{p4}{Date de début du séisme estimée}
794\annotate[yshift=-0.4em]{below}{p5}{Tableau des observations}
795\annotate[yshift=2em]{above}{p6}{Temps de propagation avec la date de début du séisme estimée}
796\annotate[yshift=5em]{above,left}{p7}{Temps de propagation que l'on devrait observer compte tenu du modèle}
797\annotate[yshift=6.5em]{above,left}{p8}{Calcul de l'angle entre le seismographe et la position estimée du séisme}
798\annotate[yshift=-1em]{below,left}{p9}{Renormalisation pour éviter de favoriser les temps de propagation longs}
799\annotate[yshift=3em]{above,left}{p11}{On fait la moyenne quadratique pour avoir l'écart}
800\annotate[yshift=0em]{below,left}{p10}{Idem mais pour l'onde P}
801\vspace*{3em}
802
803\begin{center}
804 \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-31T10:14:11,496825900+02:00.png}\\
805 \captionof{figure}{Implémentation de la fonction d'erreur}
806 \end{center}
807\subsubsection{Méthode de Nelder-Mead}
808La 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.
809
810Elle 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).
811
812L'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.
813
814\begin{center}
815\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}
816 \captionof{figure}{Une itération de Nelder-Mead sur un espace de dimension 2}
817\end{center}
818
819\subsection{Magnitude sismique}
820On utilise l'echelle logarithmique de Richter pour quantifier la magnitude d'un séisme pour laquelle on a la formule suivante:
821\[
822M_\mathrm{L} = \log_{10} \left[ \frac{A}{A_\mathrm{0}(\delta)} \right]
823\]
824
825\noindent
826où $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).
827
828Pour 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.
829
830Il existe plusieurs formules et tables empiriques pour calculer la magnitude locale d'un séisme.
831
832On utilisera la formule empirique de Tsuboi (Université de Tokyo):
833\[
834M_{\mathrm{L}} = \log_{10} A + 1.73 \log_{10} \Delta - 0.83
835\]
836
837\noindent
838où \( A \) est l'amplitude maximale en micromètres et \( \Delta \) est la distance en kilomètres.
839
840Les différentes valeurs obtenues permettent de quantifier les dégats provoqués par les séismes:
841
842\begin{center}
843\begin{tabular}{|>{\columncolor{white}}l|l|l|p{5cm}|}
844\hline
845\rowcolor{gray!30}
846\textbf{Magnitude} & \textbf{Description} & \textbf{MMI Typique} & \textbf{Effets Moyens du Séisme} \\
847\hline
848\cellcolor{mmi1}1.0 - 1.9 & Micro & I & Micro-séismes, non ressentis. Enregistrés par les sismographes. \\
849\hline
850\cellcolor{mmi1}2.0 - 2.9 & Mineur & I & Légèrement ressenti par certaines personnes. Aucun dommage aux bâtiments. \\
851\hline
852\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. \\
853\hline
854\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. \\
855\hline
856\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. \\
857\hline
858\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. \\
859\hline
860\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. \\
861\hline
862\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. \\
863\hline
864\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. \\
865\hline
866\end{tabular}
867\end{center}
868
869\section{Toutes les pièces mises ensembles}
870\subsection{Fonctionnement}
871Ainsi, 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.
872
873Il 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.
874
875Par 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).
876
877\begin{center}
878\includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-21-53.png}
879\captionof{figure}{Les données de stations telles qu'elles sont stockées pour éviter de les recalculer}
880\end{center}
881
882\newpage
883Voici une carte des stations utilisées:
884
885\begin{center}
886\includegraphics[width=12cm, trim={0 0 0 6cm}, clip]{stations.png}
887\captionof{figure}{Une carte des stations utilisées}
888\end{center}
889
890Voici ce que renvoie le programme:
891
892\begin{center}
893\includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1423-40-06.png}
894\captionof{figure}{Phase de registration des stations}
895
896\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1423-21-09-1.png}
897\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1423-21-09-2.png}
898\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1423-21-09-3.png}
899\captionof{figure}{Arrivée d'une onde}
900
901\includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-08-06.png}
902\captionof{figure}{Phase de calibration terminée}
903
904\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1501-18-07-1.png}
905\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1501-18-07-2.png}
906\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1501-18-07-3.png}
907\includegraphics[height=0.35cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1501-18-07-4.png}
908\captionof{figure}{Détection d'un séisme de magnitude 3}
909\end{center}
910
911Il y a même une interface web pour voir tout ça en temps réel (notamment les traces des sismographes).
912
913\begin{center}
914\includegraphics[width=12cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-14-28.png}
915\captionof{figure}{Interface web}
916\end{center}
917
918\newpage
919\subsection{Résultats}
920Voici par exemple un séisme plutôt faible (M6.2) detecté par mon système.
921
922On constate que les positions prédites sont en accord avec celles données par l'USGS (United States Geological Survey) ainsi que les magnitudes.
923
924\begin{center}
925
926 \includegraphics[width=12cm, trim={0 0 0 0}, clip]{debugMaps/rs2025fwmrzv.png}
927 \captionof{figure}{rs2025fwmrzv - Mon estimation - Il donne également une magnitude de M5.9}
928
929
930 \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-25T21:07:44,659295790+01:00-cropped.png}
931 \captionof{figure}{rs2025fwmrzv - Estimation de l'USGS - M6.2 selon eux}
932\end{center}
933\newpage
934
935\section{Précision et erreur}
936\subsection{Cas réduit avec propagation rectiligne dans un millieu transparent, homogène, isotrope, linéaire à longues distances et 4 detecteurs}
937Dans 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}$.
938
939Soit $\vec{M} = (x, y, z)^\top$ notre mesure et $t$ le temps d'émission du signal par tous les satellites (de façon synchrone).
940
941\begin{figure}[h]
942\centering
943\tdplotsetmaincoords{80}{140} % angle de vue (élévation, azimuth)
944\begin{tikzpicture}[tdplot_main_coords, scale=2]
945
946 % Origine (récepteur)
947 \filldraw[black] (0,0,0) circle (1.5pt) node[below] {$\vec{M'}=0$};
948
949 % Axes optionnels
950 \draw[->] (0,0,0) -- (1.5,0,0) node[left] {$x$};
951 \draw[->] (0,0,0) -- (0,1.5,0) node[right] {$y$};
952 \draw[->] (0,0,0) -- (0,0,1.5) node[above] {$z$};
953
954 % Satellite a
955 \coordinate (A) at (2,1,2);
956 \filldraw[blue] (A) circle (1.5pt);
957 \draw[-{Stealth}, blue] (0,0,0) -- (A) node[right] {$\vec{a}$};
958 \draw[-{Stealth}, red, thick] (0,0,0) -- (0.666667,0.333333,0.666667) node[midway,right] {$\vec{u}_a$};
959
960 % Satellite b
961 \coordinate (B) at (-1,2,2);
962 \filldraw[blue] (B) circle (1.5pt);
963 \draw[-{Stealth}, blue] (0,0,0) -- (B) node[right] {$\vec{b}$};
964 \draw[-{Stealth}, red, thick] (0,0,0) -- (-0.333333,0.666667,0.666667) node[midway,left] {$\vec{u}_b$};
965
966 % Satellite c
967 \coordinate (C) at (1,-2,2);
968 \filldraw[blue] (C) circle (1.5pt);
969 \draw[-{Stealth}, blue] (0,0,0) -- (C) node[right] {$\vec{c}$};
970 \draw[-{Stealth}, red, thick] (0,0,0) -- (0.333333,-0.666667,0.666667) node[midway,right] {$\vec{u}_c$};
971
972 % Satellite d
973 \coordinate (D) at (-1,-1,3);
974 \filldraw[blue] (D) circle (1.5pt);
975 \draw[-{Stealth}, blue] (0,0,0) -- (D) node[right] {$\vec{d}$};
976 \draw[-{Stealth}, red, thick] (0,0,0) -- (-0.3015113445,-0.3015113445,0.9045340337) node[midway,right] {$\vec{u}_d$};
977
978\end{tikzpicture}
979\caption{Géométrie GPS : satellites et vecteurs unitaires $\vec{u}_i$ associés}
980\end{figure}
981
982\begin{align}
983\|\vec{M} - \vec{a}\|^2 &= v^2(\tau_a - t)^2 \\
984\|\vec{M} - \vec{b}\|^2 &= v^2(\tau_b - t)^2 \\
985\|\vec{M} - \vec{c}\|^2 &= v^2(\tau_c - t)^2 \\
986\|\vec{M} - \vec{d}\|^2 &= v^2(\tau_d - t)^2
987\end{align}
988
989On é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$.
990
991$\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.
992
993\newpage
994On a donc en développant à l'ordre 1 et en identifiant :
995\begin{align}
996-2\langle \vec{\delta M}, \vec{a} \rangle &= 2v^2 \tau_a'(\delta\tau_a - \delta t) \\
997-2\langle \vec{\delta M}, \vec{b} \rangle &= 2v^2 \tau_b'(\delta\tau_b - \delta t) \\
998-2\langle \vec{\delta M}, \vec{c} \rangle &= 2v^2 \tau_c'(\delta\tau_c - \delta t) \\
999-2\langle \vec{\delta M}, \vec{d} \rangle &= 2v^2 \tau_d'(\delta\tau_d - \delta t)
1000\end{align}
1001
1002Ce qui donne en réarrangeant :
1003\begin{align}
1004\langle \vec{\delta M}, \vec{a} \rangle + v^2 \tau_a'(\delta\tau_a - \delta t) &= 0 \\
1005\langle \vec{\delta M}, \vec{b} \rangle + v^2 \tau_b'(\delta\tau_b - \delta t) &= 0 \\
1006\langle \vec{\delta M}, \vec{c} \rangle + v^2 \tau_c'(\delta\tau_c - \delta t) &= 0 \\
1007\langle \vec{\delta M}, \vec{d} \rangle + v^2 \tau_d'(\delta\tau_d - \delta t) &= 0
1008\end{align}
1009
1010On 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.
1011\begin{align}
1012\langle \vec{\delta M}, \vec{u}_a \rangle + v(\delta\tau_a - \delta t) &= 0 \\
1013\langle \vec{\delta M}, \vec{u}_b \rangle + v(\delta\tau_b - \delta t) &= 0 \\
1014\langle \vec{\delta M}, \vec{u}_c \rangle + v(\delta\tau_c - \delta t) &= 0 \\
1015\langle \vec{\delta M}, \vec{u}_d \rangle + v(\delta\tau_d - \delta t) &= 0
1016\end{align}
1017
1018Si on pose :
1019
1020\[
1021A = \frac{1}{v}
1022\begin{pmatrix}
1023u_{ax} & u_{ay} & u_{az} & -v \\
1024u_{bx} & u_{by} & u_{bz} & -v \\
1025u_{cx} & u_{cy} & u_{cz} & -v \\
1026u_{dx} & u_{dy} & u_{dz} & -v
1027\end{pmatrix},
1028\quad
1029X = \begin{pmatrix} \delta x \\ \delta y \\ \delta z \\ \delta t \end{pmatrix},
1030\quad
1031Y = \begin{pmatrix} \delta\tau_a \\ \delta\tau_b \\ \delta\tau_c \\ \delta\tau_d \end{pmatrix}
1032\]
1033
1034On a $AX = Y$.
1035
1036Ainsi 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é).
1037
1038\[
1039\|X\|^2 =X^\top X = Y^\top {A^{-1}}^\top A^{-1} Y = Y^\top (AA^\top)^{-1} Y
1040\]
1041
1042$(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.
1043
1044\[
1045Q = (\vec{v}_1 \mid \vec{v}_2 \mid \vec{v}_3 \mid \vec{v}_4)
1046\]
1047
1048où $\vec{v}_1$, $\vec{v}_2$, $\vec{v}_3$ et $\vec{v}_4$ sont orthonormés.
1049
1050\begin{align}
1051\|X\|^2 = X^\top X &= Y^\top {A^{-1}}^\top A^{-1} Y \\
1052&= Y^\top Q^\top \operatorname{diag}(\lambda_1, \lambda_2, \lambda_3, \lambda_4)\, Q\, Y \\
1053&= \sum_{i=1}^{4} \lambda_i \langle \vec{v}_i, Y \rangle^2
1054\end{align}
1055
1056On passe à l'espérance :
1057
1058\begin{align}
1059\mathbb{E}(\|X\|^2) &= \sum_{i=1}^{4} \lambda_i \,\mathbb{E}\!\left(\langle \vec{v}_i, Y \rangle^2\right) \\
1060&= \sum_{i=1}^{4} \lambda_i \sum_{j=1}^{4} (v_i)_j^2 \,\sigma_j^2
1061\end{align}
1062
1063La dernière égalité étant donnée par indépendance.
1064
1065Si 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 :
1066
1067\[
1068\mathbb{E}(\|X\|^2) = \operatorname{tr}\!\left((AA^\top)^{-1}\right)\sigma_\tau^2
1069\]
1070
1071On 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).
1072
1073Si 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$.
1074
1075\begin{figure}[h]
1076\centering
1077\tdplotsetmaincoords{70}{135}
1078\begin{tikzpicture}[tdplot_main_coords, scale=1.5]
1079
1080 % j'ai décidé de tourner par la matrice
1081 %[ 0.8636654, 0.4567464, -0.2132249;
1082 %-0.4446277, 0.8895858, 0.1046102;
1083 % 0.2374622, 0.0044575, 0.9713866 ]
1084
1085 % Ellipsoide Q(Y)=C projeté en 3D (delta_tau_d fixé)
1086 % On visualise dans l'espace (delta_tau_a, delta_tau_b, delta_tau_c)
1087
1088 % Axes
1089 \draw[-{Stealth}, gray] (0,0,0) -- (4,0,0) node[left] {$\delta\tau_a$};
1090 \draw[-{Stealth}, gray] (0,0,0) -- (0,4,0) node[right] {$\delta\tau_b$};
1091 \draw[-{Stealth}, gray] (0,0,0) -- (0,0,2) node[above] {$\delta\tau_c$};
1092
1093 % Ellipsoide (demi-axes 1/sqrt(lambda_i), ici lambda_1<lambda_2<lambda_3)
1094 % demi-axes : a1=2.0, a2=1.2, a3=0.7 pour illustrer lambda_1 petit
1095 \def\a{2.4} % grand axe -> petite vp -> grande erreur
1096 \def\b{1.9}
1097 \def\c{1.0} % petit axe -> grande vp -> petite erreur
1098
1099 % Meridiens
1100 \foreach \ph in {0,30,...,150}{
1101 \draw[blue!30, thin] plot[domain=0:360, samples=60]
1102 ({0.8636654*\a*cos(\x)*cos(\ph)+0.4567464*\b*cos(\x)*sin(\ph)-0.2132249*\c*sin(\x)},
1103 {-0.4446277*\a*cos(\x)*cos(\ph)+0.8895858*\b*cos(\x)*sin(\ph)+0.1046102*\c*sin(\x)},
1104 {0.2374622*\a*cos(\x)*cos(\ph)+0.0044575*\b*cos(\x)*sin(\ph)+0.9713866*\c*sin(\x)});
1105 }
1106
1107 % Paralleles
1108 \foreach \th in {-60,-30,...,60}{
1109 \draw[blue!50, thin] plot[domain=0:360, samples=60]
1110 ({0.8636654*\a*cos(\th)*cos(\x)+0.4567464*\b*cos(\th)*sin(\x)-0.2132249*\c*sin(\th)},
1111 {-0.4446277*\a*cos(\th)*cos(\x)+0.8895858*\b*cos(\th)*sin(\x)+0.1046102*\c*sin(\th)},
1112 {0.2374622*\a*cos(\th)*cos(\x)+0.0044575*\b*cos(\th)*sin(\x)+0.9713866*\c*sin(\th)});
1113 }
1114
1115 % Cercles de l'ellipsoide sur les plans principaux
1116 \draw[blue, thick] plot[domain=0:360, samples=60]
1117 ({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)});
1118 % ({\a*cos(\x)}, {\b*sin(\x)}, 0);
1119
1120 \draw[blue, thick] plot[domain=0:360, samples=60]
1121 ({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)});
1122 % ({\a*cos(\x)}, 0, {\c*sin(\x)});
1123
1124 \draw[blue, thick] plot[domain=0:360, samples=60]
1125 ({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)});
1126 % (0, {\b*cos(\x)}, {\c*sin(\x)});
1127 % Vecteurs propres v_i (directions des axes de l'ellipsoide)
1128 \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$};
1129 \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$};
1130 \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$};
1131
1132 % Origine
1133 \filldraw[black] (0,0,0) circle (1.5pt) node[below left] {$0$};
1134
1135 % Légende
1136 \node at (0,-2.5,1) {($\delta\tau_d = 0$ fixé)};
1137
1138\end{tikzpicture}
1139\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)$.
1140Les grands axes correspondent aux petites valeurs propres de $AA^\top$ :
1141une faible norme de $Y$ dans ces directions suffit à atteindre la surface,
1142signifiant une grande sensibilité à l'erreur.}
1143\end{figure}
1144\newpage
1145\subsection{Méthode de Monte-Carlo pour les cas plus complexes}
1146Pour 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.
1147
1148J'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).
1149
1150D'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).
1151
1152On voit bien l'effet de la dilution de précision géométrique.
1153
1154Voici une configuration écrasée, il y a une grande dilution de précision géométrique dans le plan engendré par les detecteurs.
1155
1156\begin{figure}[h]
1157 \centering
1158 \begin{minipage}{0.49\textwidth}
1159 \centering
1160 \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/1.png}
1161 \caption{Rouge: détecteurs, Vert: position mesurée, Bleu: référence (configuration écrasée)}
1162 \end{minipage}
1163 \hfill
1164 \begin{minipage}{0.49\textwidth}
1165 \centering
1166 \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/2.png}
1167 \caption{Erreur calculée par Monte-Carlo selon la position dans l'espace à erreur de mesure égale par détecteur (configuration écrasée)}
1168 \end{minipage}
1169\end{figure}
1170
1171\newpage
1172Lorsque la configuration donne des vecteurs moins colinéaires on a bien une meilleure précision.
1173
1174\begin{figure}[h]
1175 \centering
1176 \begin{minipage}{0.49\textwidth}
1177 \centering
1178 \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/3.png}
1179 \caption{Rouge: détecteurs, Vert: position mesurée, Bleu: référence (configuration favorable)}
1180 \end{minipage}
1181 \hfill
1182 \begin{minipage}{0.49\textwidth}
1183 \centering
1184 \includegraphics[width=\linewidth, trim={0 0 0 0}, clip]{montecarlo/4.png}
1185 \caption{Erreur calculée par Monte-Carlo selon la position dans l'espace à erreur de mesure égale par détecteur (configuration favorable)}
1186 \end{minipage}
1187\end{figure}
1188
1189
1190
1191\newpage
1192\section{Bonus}
1193De 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}).
1194
1195\begin{center}
1196\includegraphics[width=5cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-42-57.png}
1197\captionof{figure}{Mode canal unique avec auto-correlation}
1198\includegraphics[width=5cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1500-41-26.png}
1199\captionof{figure}{Mode XY}
1200\includegraphics[width=5cm, trim={0 0 0 0}, clip]{screenshots/Screenshotfrom2026-05-1512-44-34.png}
1201\captionof{figure}{Mode spectrogramme avec détection de note}
1202\end{center}
1203
1204\newpage
1205\section{Bibliographie}
1206\begin{thebibliography}{9}
1207
1208\bibitem{fft_matrix_factorizations}
1209Charles Van Loan,
1210\textit{The FFT Via Matrix Factorizations},
1211Lecture notes, 2010.\\
1212\url{https://www.cs.cornell.edu/~bindel/class/cs5220-s10/slides/FFT.pdf}
1213
1214\bibitem{eq_sources}
1215R.J. Mitchell,
1216\textit{Earthquake Sources},
1217Lecture notes, Western Washington University.\\
1218\url{https://www.geol.wwu.edu/rjmitch/L4_EQsources.pdf}
1219
1220\bibitem{seismic_waves}
1221University of Hawaii,
1222\textit{Compare, Contrast, and Connect: Seismic Waves and Determining Earth's Structure}.\\
1223\url{https://manoa.hawaii.edu/exploringourfluidearth/physical/ocean-floor/layers-earth/compare-contrast-connect-seismic-waves-and-determining-earth-s-structure}
1224
1225\bibitem{gfz-1-d-earth-models}
1226Peter Bormann,
1227\textit{Global 1-D Earth models (IASP91 tables)},
1228GFZ German Research Centre for Geosciences.\\
1229\url{https://gfzpublic.gfz-potsdam.de/rest/items/item_4031/component/file_4032/content}
1230
1231\bibitem{earthquake_data_centers}
1232Yacine Boussoufa,
1233\textit{Earthquake Data Centers},
1234GitHub repository.\\
1235\url{https://github.com/YacineBoussoufa/EarthquakeDataCenters}
1236
1237\bibitem{cfcs_lecture}
1238School of Informatics, University of Edinburgh,
1239\textit{Computational Foundations of Cognitive Science: Lecture 15 (Convolutions and Kernels)}.\\
1240\url{https://www.inf.ed.ac.uk/teaching/courses/cfcs1/lectures/cfcs_l15.pdf}
1241
1242\bibitem{penn_state_anova}
1243Penn State Eberly College of Science,
1244\textit{STAT 510: Lesson 8.2 - Cross Correlation Functions and Lagged Regressions},
1245Online course material.\\
1246\url{https://online.stat.psu.edu/stat510/lesson/8/8.2}
1247
1248\bibitem{nelder_mead}
1249Jason Cantarella,
1250\textit{Nelder-Mead Method},
1251Lecture notes.\\
1252\url{https://jasoncantarella.com/downloads/NelderMeadProof.pdf}
1253
1254\bibitem{broadband_magnitude}
1255Tatsuhiko Hara,
1256\textit{Determination of Broadband Moment Magnitude},
1257IISEE/BRI.\\
1258\url{https://iisee.kenken.go.jp/lna/download.php?f=2011082925678c01.pdf&n=T0-100-2007_Mwp-2-new.pdf&cid=T0-100-2007}
1259
1260\bibitem{obspy}
1261Moritz Beyreuther, Robert Barsch, Lion Krischer, Tobias Megies, Yannik Behr and Joachim Wassermann,\\
1262\textit{ObsPy: A Python Toolbox for Seismology},\\
1263Seismological Research Letters, vol. 81, no. 3, pp. 530--533, 2010.\\
1264doi:10.1785/gssrl.81.3.530
1265
1266\bibitem{fourieranalysis}
1267Elias M. Stein \& Rami Shakarchi. \textit{Fourier Analysis: An Introduction.}\\
12682003. (ISBN 978-0-691-11384-5)
1269
1270\bibitem{convolution_video}
12713Blue1Brown,
1272\textit{But what is a convolution?},
1273YouTube video, 2022.\\
1274\url{https://www.youtube.com/watch?v=KuXjwB4LzSA}
1275
1276\bibitem{fft_algorithm_video}
1277Michael Pound,
1278\textit{The Fast Fourier Transform Algorithm},
1279YouTube video, 2023.\\
1280\url{https://www.youtube.com/watch?v=toj_IoCQE-4}
1281
1282
1283\end{thebibliography}
1284
1285\end{document}