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{assets/texpackages/annotate-equations}
31\newenvironment{graytext}{\color{gray}}{\ignorespacesafterend}
32\graphicspath{ {./assets/} }
33
34% Define custom colors
35\definecolor{mmi1}{HTML}{FFFFFF}
36\definecolor{mmi2}{HTML}{BFCCFF}
37\definecolor{mmi3}{HTML}{A0E6FF}
38\definecolor{mmi4}{HTML}{80FFFF}
39\definecolor{mmi5}{HTML}{7AFF93}
40\definecolor{mmi6}{HTML}{FFFF00}
41\definecolor{mmi7}{HTML}{FFC800}
42\definecolor{mmi8}{HTML}{FF9100}
43\definecolor{mmi9}{HTML}{FF0000}
44\definecolor{mmi10}{HTML}{C80000}
45\definecolor{mmi11}{HTML}{A40000}
46\definecolor{mmi12}{HTML}{800000}
47
48\begin{document}
49\title{Traitement de signaux pour la détection de séismes et leur multilatération}
50\author{Dalibard Louis}
51\date{\today}
52
53\maketitle
54\tableofcontents
55
56\newpage
57\section{Séismes}
58\subsection{Introduction}
59Un séisme ou tremblement de terre est une secousse du sol résultant de la libération brusque d'énergie accumulée par les contraintes exercées sur les roches. Cette libération d'énergie se fait par rupture le long d'une faille, généralement préexistante. Plus rares sont les séismes dus à l'activité volcanique ou d'origine artificielle (explosions par exemple). Le lieu de la rupture des roches en profondeur se nomme le foyer ; la projection du foyer à la surface est l'épicentre du séisme. Le mouvement des roches près du foyer engendre des vibrations élastiques qui se propagent, sous la forme de paquets d'ondes sismiques, autour et au travers du globe terrestre.\\
60\includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{eq-ed-fault-labeled.png}
61\subsection{Ondes P et S}
62Les 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.
63
64La 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).
65\noindent\begin{minipage}{.5\textwidth}
66\centering
67\includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_compression_impulsion_1d_30_petit.png} \captionof{figure}{Ondes P}
68\end{minipage}%
69\begin{minipage}{.5\textwidth}
70\centering
71\includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_cisaillement_impulsion_1d_30_petit.png}
72 \captionof{figure}{Ondes S}
73\end{minipage}
74\subsection{Différence de vitesse de propagation}
75Les 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.
76
77Les 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.
78\begin{center}
79 \includegraphics[width=7cm, trim={0 12.5mm 0cm 2cm}, clip]{2025-03-26T23:29:53,572838335+01:00.png}\\
80 \captionof{figure}{6 km/s (ondes P) vs 4 km/s (ondes S)}
81 \end{center}
82
83On peut donc exploiter les temps d'arrivée de ces différentes ondes pour localiser un séisme.
84
85\section{Théorie}
86\subsection{Principe}
87Différentes étapes:
88 \begin{enumerate}
89 \item Acquisition de données en temps réel (SeedLink)
90 \item Reconnaissance d'un séisme et mesure automatique des temps
91 \item Calcul de la position et de la magnitude
92 \begin{enumerate}
93 \item Modélisation de la propagation des ondes sismiques
94 \item Méthode numérique d'optimisation de fonction à plusieurs variables pour la multilatération
95 \item Calcul de la magnitude
96 \end{enumerate}
97 \end{enumerate}
98
99\subsection{DSP}
100\subsubsection{Acquisition des données}
101On utilise le protocole Seedlink pour récupérer les données des sismographes en temps réel.
102
103\renewcommand{\arraystretch}{1.2}
104
105
106\begin{longtable}{|l|l|}
107 \hline
108 \textbf{Name} & \textbf{Host} \\
109 \hline
110 AusPass & auspass.edu.au \\
111 BGR & eida.bgr.de \\
112 CISMID & www.cismid.uni.edu.pe \\
113 ENS & ephesite.ens.fr \\
114 GEOFON, GFZ & geofon.gfz-potsdam.de \\
115 GEONET & link.geonet.org.nz \\
116 Geoscience Australia & seis-pub.ga.gov.au \\
117 GSRAS (?) & 89.22.182.133 \\
118 Helsinki & finseis.seismo.helsinki.fi \\
119 Haiti & ayiti.unice.fr \\
120 ICGC & ws.icgc.cat \\
121 IDA Project & rtserve.ida.ucsd.edu \\
122 IFZ & data.ifz.ru \\
123 IPGP & rtserver.ipgp.fr \\
124 IRIS DMC & rtserve.iris.washington.edu \\
125 IRIS Jamaseis & jamaseis.iris.edu \\
126 ISNET - UNINA & 185.15.171.86 \\
127 LMU & erde.geophysik.uni-muenchen.de \\
128 NIGGG & 195.96.231.100 \\
129 NRCAN & earthquakescanada.nrcan.gc.ca \\
130 OBSEBRE & obsebre.es \\
131 OGS & nam.ogs.it \\
132 Oklahoma University & rtserve.ou.edu \\
133 ORFEUS & eida.orfeus-eu.org \\
134 PLSN (IGF Poland) & hudson.igf.edu.pl \\
135 Red Sìsmica de Puerto Rico & 161.35.236.45 \\
136 Red Sìsmica Baru & helis.redsismicabaru.com \\
137 RESIF & rtserve.resif.fr \\
138 SANET & 147.213.113.73 \\
139 RSIS & rsis1.on.br \\
140 SCSN-USC (South Carolina Seismic Network) & eeyore.seis.sc.edu:6382 \\
141 Seisme IRD & rtserve.ird.nc \\
142 Staneo & vibrato.staneo.fr \\
143 SNAC NOA & snac.gein.noa.gr \\
144 TexNet & rtserve.beg.utexas.edu \\
145 Thai Meteorological Department & 119.46.126.38 \\
146 UFRN (Universidade Federal do Rio Grande do Norte) & sislink.geofisica.ufrn.br \\
147 Unical Universita Della Calabria & www.sismocal.org \\
148 UNITS Università degli studi di Trieste & rtweb.units.it \\
149 UNIV-AG Université des Antilles & seedsrv0.ovmp.martinique.univ-ag.fr \\
150 Universidade de Évora & clv-cge.uevora.pt \\
151 Universidad de Colima & 148.213.24.15 \\
152 UPR & worm.uprm.edu \\
153 USGS & cwbpub.cr.usgs.gov \\
154 USP-IAG & seisrequest.iag.usp.br \\
155 \hline
156\end{longtable}
157\subsubsection{Extraction des temps d'arrivée}
158\begin{center}
159 \includegraphics[width=15cm, trim={0 0 0 0}, clip]{R0ED0.EHZ-1743142835749.9944-cropped.png}\\
160 \captionof{figure}{Exemple d'un enregistrement de sismographe}
161\end{center}
162On 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 ?
163Pour cela on va faire un détour par des maths.
164
165\subsubsection{Convolution}
166Si on a deux fonctions, quelles opérations peut-on faire pour avoir une nouvelle fonction?
167
168La plus simple semble être les additionner.
169
170\begin{center}
171\begin{tikzpicture}
172 \begin{axis}[
173 axis lines = middle,
174 xlabel = {$x$},
175 ylabel = {$y$},
176 samples=100,
177 domain=-6:6,
178 legend pos=north east,
179 width=0.8\textwidth,
180 height=.4\textheight
181 ]
182 % Define the functions
183 \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$};
184 \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$};
185 \addplot[green, thick] {sin(deg(x)) + cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) + g(x)$};
186
187 % Add legend
188 \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) + g(x)$}
189 \end{axis}
190\end{tikzpicture}
191\end{center}
192
193
194
195Une autre façon serait de les multiplier.
196
197\begin{center}
198\begin{tikzpicture}
199 \begin{axis}[
200 axis lines = middle,
201 xlabel = {$x$},
202 ylabel = {$y$},
203 samples=100,
204 domain=-6:6,
205 legend pos=north east,
206 width=0.8\textwidth,
207 height=.4\textheight
208 ]
209 % Define the functions
210 \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$};
211 \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$};
212 \addplot[green, thick] {sin(deg(x)) * cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) \cdot g(x)$};
213
214 % Add legend
215 \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) \cdot g(x)$}
216 \end{axis}
217\end{tikzpicture}
218\end{center}
219
220Mais il existe aussi une autre opération, la convolution.
221
222\begin{center}
223\begin{tikzpicture}
224 \begin{axis}[
225 axis lines = middle,
226 xlabel = {$x$},
227 ylabel = {$y$},
228 samples=100,
229 domain=-6:6,
230 legend pos=north east,
231 width=0.8\textwidth,
232 height=.4\textheight
233 ]
234 % Define the functions
235 \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$};
236 \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$};
237 \addplot[green, thick] {0} node[right] {\footnotesize $h(x) = (f*g)(x)$};
238
239 % Add legend
240 \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = (f*g)(x) = 0$}
241 \end{axis}
242\end{tikzpicture}
243\end{center}
244
245Qui 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).
246
247On définit formellement la convolution comme un produit sur l'espace des fonctions intégrables.
248
249Soit $f$ et $g$ deux fonctions intégrables sur $\mathbb{R}$.
250
251Pour tout $t \in \mathbb{R}$,
252\[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \]
253
254Dans le cas discret on utilisera une somme.
255
256Soit $f$ et $g$ deux fonctions de $\mathbb{Z}$ dans $\mathbb{C}$.
257
258Pour tout $n \in \mathbb{Z}$,
259\[ (f*g)[n]=\sum_{m=-\infty}^{+\infty} f[m] \cdot g[n-m] \]
260
261et dans le cas de fonction périodiques, il suffit de faire l'intégrale sur une période.
262Pour tout $t \in \mathbb{R}$,
263\[ (f*g)(t)=\int_{0}^{T} f(x) \cdot g(t-x) dx \]
264
265(on peut également généraliser cela à des fonctions discretès périodiques)
266
267\subsubsection{Propriétés algébriques de la convolution}
268\begin{itemize}
269\item Commutatif
270
271On remarquera que si \[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \]
272
273Et on fait le changement de variable $u=t-x$
274
275On 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)\]
276
277\item Associatif
278\[ (f*g)*h= f*(g*h) \]
279C'est une conséquence du théorème de Fubini.
280
281\item Distributif
282\[ f*(g+h)= f*g+f*h \]
283Par linéarité de l'intégrale.
284
285\end{itemize}
286
287L'espace des fonctions intégrables muni de $*$ forme un demi-groupe commutatif (car pas d'élement neutre).
288
289Pourquoi une telle construction ?
290
291Déjà à quoi ça ressemble?
292
293Si 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.
294
295\begin{center}
296 \begin{tikzpicture}[scale=0.5]
297 % Define Dice Faces
298 \def\diceA{1,2,3,4,5,6} % Top Die (Static)
299 \def\diceB{6,5,4,3,2,1} % Bottom Die (Sliding)
300
301 % Define step spacing
302 \def\stepSpace{3}
303
304 % Loop through each step
305 \foreach \step in {-5,-4,-3,-2,-1,0,1,2,3,4,5} {
306 \pgfmathsetmacro{\sommedonnant}{int(7 + \step)}
307 % Label step number
308 \node[anchor=east] at (-0.5, -\step*3) {Somme donnant \sommedonnant:};
309
310 % Top Row - Fixed Dice
311 \foreach \x [count=\i] in \diceA {
312 \node[draw, minimum size=0.5cm] at (\i, -\step*3) {\x};
313 }
314
315 % Bottom Row - Sliding Dice
316 \foreach \x [count=\i] in \diceB {
317 \pgfmathsetmacro{\sumval}{int(\i + \step)}
318 \ifnum \sumval > 0 \ifnum \sumval < 7 \node[draw=red, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x};
319 \else
320 \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x};
321 \fi
322 \else
323 \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x};
324 \fi
325 }
326 }
327 \end{tikzpicture}
328\end{center}
329\begin{center}
330 \begin{tikzpicture}
331 \begin{axis}[
332 xlabel={Somme},
333 ylabel={Fréquence},
334 ymin=0,
335 ymax=7,
336 xmin=1.5,
337 xmax=12.5,
338 xtick={2,3,4,5,6,7,8,9,10,11,12},
339 ytick={1,2,3,4,5,6,7},
340 area style,
341 width=12cm,
342 height=8cm,
343 major grid style={line width=.2pt,draw=gray!50},
344 grid=major,
345 bar width=0.8,
346 title={Distribution des sommes},
347 nodes near coords,
348 ]
349 \addplot[fill=blue!40, draw=blue!80, opacity=0.8] coordinates {
350 (2,1)
351 (3,2)
352 (4,3)
353 (5,4)
354 (6,5)
355 (7,6)
356 (8,5)
357 (9,4)
358 (10,3)
359 (11,2)
360 (12,1)
361 };
362 \end{axis}
363\end{tikzpicture}
364\end{center}
365
366On 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.
367
368On vient de faire la convolution entre cette fonction et elle même.
369\begin{center}
370 \begin{tikzpicture}
371 \begin{axis}[
372 xlabel={$x$},
373 ylabel={$f(x)$},
374 xmin=-1,
375 xmax=8,
376 ymin=-0.2,
377 ymax=1.4,
378 xtick={-1,0,1,2,3,4,5,6,7,8},
379 ytick={0,1},
380 width=10cm,
381 height=6cm,
382 grid=major,
383 title={Fonction nulle puis constante égale à $1$ puis nulle},
384 samples=100,
385 domain=-1:8,
386 axis lines=middle,
387 clip=false,
388 ]
389 % Plot the function
390 \addplot[blue, thick, const plot] coordinates {
391 (-1, 0)
392 (0.999, 0)
393 (1, 1)
394 (6, 1)
395 (6.001, 0)
396 (8, 0)
397 };
398
399 % Add markers to explicitly show the endpoints
400 \addplot[only marks, mark=*, mark options={fill=red}] coordinates {
401 (1, 1)
402 (6, 1)
403 };
404
405 \end{axis}
406\end{tikzpicture}
407\end{center}
408\begin{center}
409\includegraphics[width=14cm, trim={0 0 0 0}, clip]{2025-04-01T23:42:09,818697014+02:00.png}
410\captionof{figure}{numpy confirme ce résultat}
411\end{center}
412
413On remarque que si on le fait dans le cas discret, on a des valeurs en plus sur les bords que l'on peut ignorer.
414
415Faire 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.
416
417\begin{center}
418 \begin{tikzpicture}
419 \begin{axis}[
420 axis lines = middle,
421 xlabel = {$x$},
422 ylabel = {$y$},
423 samples=100,
424 domain=-6:6,
425 legend pos=north east,
426 width=0.8\textwidth,
427 height=.4\textheight
428 ]
429 \addplot[blue, thick] {cos(deg(2*x))};
430 \addlegendentry{$\cos(2x)$}
431 \addplot[red, thick] {sin(deg(x))};
432 \addlegendentry{$\sin(x)$}
433 \end{axis}
434 \end{tikzpicture}
435 \end{center}
436
437 \begin{center}
438 \begin{tikzpicture}
439 \begin{axis}[
440 axis lines = middle,
441 xlabel = {$x$},
442 ylabel = {$y$},
443 samples=100,
444 domain=-6:6,
445 legend pos=north east,
446 width=0.8\textwidth,
447 height=.4\textheight
448 ]
449 \addplot[blue, thick] {cos(deg(2*x))};
450 \addlegendentry{$\cos(2x)$}
451 \addplot[red, thick] {sin(deg(-x))};
452 \addlegendentry{$\sin(-x)$}
453 \end{axis}
454 \end{tikzpicture}
455 \end{center}
456
457 \begin{center}
458 \begin{tikzpicture}
459 \begin{axis}[
460 axis lines = middle,
461 xlabel = {$x$},
462 ylabel = {$y$},
463 samples=100,
464 domain=-6:6,
465 legend pos=north east,
466 width=0.8\textwidth,
467 height=.4\textheight
468 ]
469 \addplot[blue, thick] {cos(deg(2*x))*sin(deg(-x))};
470 \addlegendentry{$\cos(2x)\sin(-x)$}
471
472 % Highlight positive areas in red
473 \addplot [fill=red, opacity=0.3, domain=-2*pi:2*pi]
474 {max(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle;
475
476 % Highlight negative areas in blue
477 \addplot [fill=blue, opacity=0.3, domain=-2*pi:2*pi]
478 {min(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle;
479
480 \end{axis}
481 \end{tikzpicture}
482 \end{center}
483
484 On voit que c'est nul dans ce cas.
485
486\subsubsection{Moyennage}
487
488Si on applique la fonction $\frac{f}{6}$ à une fonction $g$ quelconque on observe une sorte de moyennage qui se fait entre 6 valeurs proches.
489\begin{center}
490 \begin{tikzpicture}
491 \begin{axis}[
492 xlabel={$x$},
493 ylabel={$y$},
494 xmin=-0.5,
495 xmax=23.5,
496 ymin=-0.5,
497 ymax=9,
498 xtick={0,2,4,6,8,10,12,14,16,18,20,22},
499 width=12cm,
500 height=8cm,
501 grid=major,
502 title={Moyennage par convolution},
503 legend pos=north west,
504 ]
505 % Plot the first dataset with connected lines
506 \addplot[blue, thick] coordinates {
507 (0,0) (1,0) (2,0) (3,0) (4,1) (5,2) (6,3) (7,4) (8,5) (9,6) (10,7) (11,8)
508 (12,0) (13,4) (14,2) (15,3) (16,8) (17,1) (18,2) (19,4) (20,8) (21,0) (22,0) (23,0)
509 };
510
511 % Plot the second dataset with connected lines
512 \addplot[red, thick] coordinates {
513 (0,0)
514 (1,0.16666667)
515 (2,0.5)
516 (3,1)
517 (4,1.66666667)
518 (5,2.5)
519 (6,3.5)
520 (7,4.5)
521 (8,5.5)
522 (9,5)
523 (10,5)
524 (11,4.5)
525 (12,4)
526 (13,4.16666667)
527 (14,3)
528 (15,3.33333333)
529 (16,3.33333333)
530 (17,4.33333333)
531 (18,3.83333333)
532 (19,2.5)
533 (20,2.33333333)
534 (21,2)
535 (22,1.33333333)
536 (23,0)
537 };
538
539 \legend{$g$, $g*\frac{f}{6}$}
540
541 \end{axis}
542 \end{tikzpicture}
543\end{center}
544
545Ce 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.
546
547\subsubsection{Transformée de Fourier discrète}
548Si 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}]$,
549
550On 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}$
551
552
553On 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$
554
555$(a*b)[p]$ est le coefficient du terme de degré $p$ dans le produit
556
557$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$
558
559Bon, de premier abord, ça n'a pas l'air de beaucoup nous aider, le produit de polynômes est $\mathcal{O}(n^2)$.
560
561Mais 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.
562
563Et 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.
564
565Le problème c'est que cela semble prendre plus de temps, un pivot de Gauss est en $\mathcal{O}(n^3)$
566
567Mais si on choisit d'évaluer nos valeurs en les racines de l'unité, on a un système plus simple à résoudre.
568
569On va pouvoir faire les opérations en $\mathcal{O}(n \cdot log(n))$
570
571Les 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.
572\subsubsection{Radix-2 decimation-in-time (DIT) - Factorisation Cooley-Tukey}
573Evaluer 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:
574
575\[
576R=\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} =
577\begin{bmatrix}
5781 & 1 & 1 & \cdots & 1 \\
5791 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\
5801 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\
581\vdots & \vdots & \vdots & \ddots & \vdots \\
5821 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}
583\end{bmatrix}
584\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}
585\]
586
587On reconnaît la matrice de Vandermonde transposée associée à l'évaluation en les \( \omega_n^k \)
588
589On 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$.
590
591On note la matrice de Vandermonde transposée, qui permet de calculer le DFT pour tous nos coefficients,
592
593$F_{2^p}=\begin{bmatrix}
5941 & 1 & 1 & \cdots & 1 \\
5951 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\
5961 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\
597\vdots & \vdots & \vdots & \ddots & \vdots \\
5981 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)}
599\end{bmatrix}$
600
601On peut permuter les lignes, à la fin, ce qu'on a c'est que, du fait des propriétés des racines de l'unité,
602
603si on note, $D_{2^{p-1}} =
604\begin{bmatrix}
6051 & 0 & 0 & \cdots & 0 \\
6060 & \omega_{2^{p-1}} & 0 & \cdots & 0 \\
6070 & 0 & \omega_{2^{p-1}}^2 & \cdots & 0 \\
608\vdots & \vdots & \vdots & \ddots & \vdots \\
6090 & 0 & 0 & \cdots & \omega_{2^{p-1}}^{2^{p-1}-1}
610\end{bmatrix}$
611
612est:
613
614\begin{equation*}
615R=F_{2^p}\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}=
616\begin{pNiceArray}{cc|cc}
617 \eqnmarkbox[blue]{mcoeffdiag1}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag2}{D_{2^{p-1}}} \\
618 \hline
619 \eqnmarkbox[blue]{mcoeffdiag3}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag4}{-D_{2^{p-1}}}
620\end{pNiceArray}
621\begin{pNiceArray}{cc|cc}
622 \eqnmarkbox[cyan]{mrecurrence1}{F_{2^{p-1}}} && \mathbf{0} \\
623 \hline
624 \mathbf{0} && \eqnmarkbox[cyan]{mrecurrence2}{F_{2^{p-1}}}
625\end{pNiceArray}
626\begin{bmatrix}
627 \eqnmarkbox[green]{mcoeffpairs1}{a_0} \\
628 \eqnmarkbox[green]{mcoeffpairs2}{a_2} \\
629 \vdots \\
630 \eqnmarkbox[green]{mcoeffpairs3}{a_{2^{p-1}-2}} \\
631 \eqnmarkbox[green]{mcoeffpairs4}{a_{2^{p-1}}} \\
632 \eqnmarkbox[magenta]{mcoeffimpairs1}{a_{1}} \\
633 \eqnmarkbox[magenta]{mcoeffimpairs2}{a_{3}} \\
634 \vdots \\
635 \eqnmarkbox[magenta]{mcoeffimpairs3}{a_{n-3}} \\
636 \eqnmarkbox[magenta]{mcoeffimpairs4}{a_{n-1}}
637\end{bmatrix}
638\end{equation*}
639\annotatetwo[yshift=1.5em]{above}{mcoeffdiag1}{mcoeffdiag2}{Blocs diagonaux, de l'ordre de $\mathcal{O}(n)$ opérations}
640\annotate[yshift=0.5em]{above}{mcoeffpairs1}{Coefficients d'indice pair}
641\annotate[yshift=0em]{below}{mcoeffimpairs4}{Coefficients d'indice impair}
642\annotate[yshift=0em]{below,left}{mrecurrence2}{Récurrence}
643
644
645On évalue la compléxité de l'algorithme.\\
646On note $u_p$ sa complexité en fonction de $p$ et $C_n$ sa complexité en fonction de $n$.
647
648
649\begin{equation*}
650u_{p+1}=\eqnmarkbox[magenta]{cdiagmult}{A \cdot 2^{p+1}}+\eqnmarkbox[red]{crec}{2u_{p}}
651\end{equation*}
652\annotate[yshift=0.5em]{above,left}{cdiagmult}{Compléxité du produit sur la diagonale}
653\annotate[yshift=0em]{below,right}{crec}{Traitement des coefficients par récurrence}
654
655On factorise par la solution homogène,
656$\frac{u_{p+1}}{2^{p+1}}=A+\frac{u_{p}}{2^p}$
657
658$\frac{u_{p}}{2^{p}}=u_{0}+A\cdot p$
659
660$u_{p}=u_{0}\cdot2^{p}+A\cdot p \cdot 2^{p}$
661
662Or $p=\log_2{n}$
663
664Donc $C_n = u_{\log_2{n}} = u_{0} \cdot n+A \cdot \log_2{n} \cdot n = \mathcal{O}(n \log n)$
665
666\subsubsection{IFFT}
667On peut montrer que:
668$F_{2^p}^{-1}=\frac{1}{2^p}\begin{bmatrix}
6691 & 1 & 1 & \cdots & 1 \\
6701 & \overline{\omega_n}^1 & \overline{\omega_n}^2 & \cdots & \overline{\omega_n}^{n-1} \\
6711 & \overline{\omega_n}^2 & \overline{\omega_n}^4 & \cdots & \overline{\omega_n}^{2(n-1)} \\
672\vdots & \vdots & \vdots & \ddots & \vdots \\
6731 & \overline{\omega_n}^{n-1} & \overline{\omega_n}^{2(n-1)} & \cdots & \overline{\omega_n}^{(n-1)(n-1)}
674\end{bmatrix}$
675
676$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}$
677
678Le 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).
679
680$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}}$
681
682On 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.
683
684\subsubsection{Passage du domaine temporel au domaine fréquentiel}
685Sans 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.
686
687Maintenant qu'on peut calculer des convolutions rapidement, revenons à notre problème.
688
689% TODO: https://fr.wikipedia.org/wiki/Transformation_de_Fourier_rapide
690\subsubsection{Corrélation croisée}
691% TODO: https://en.wikipedia.org/wiki/Cross-correlation
692Si on définit $\tilde{f}(t)=f(-t)$
693
694La corrélation croisée de \(f\) et \(g\) est \( (g * \tilde{f})(t) = \int_{-\infty}^{\infty} \overline{f(x-t)} g(x) \, dx \)
695
696Ou sur une seule période pour une fonction périodique.
697
698On 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.
699
700Quand 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.
701
702On 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.
703
704Et donc on ne regardera que la partie réele et à quel point elle est grande.
705
706Donc 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$.
707
708\begin{center}
709 \begin{tikzpicture}
710 \begin{axis}[
711 axis lines = middle,
712 xlabel = {$x$},
713 ylabel = {$y$},
714 samples=100,
715 domain=-6:6,
716 legend pos=north east,
717 width=0.8\textwidth,
718 height=.4\textheight
719 ]
720 \addplot[blue, thick] {cos(deg(x))};
721 \addlegendentry{$g(x)=\cos(x)$}
722 \addplot[red, thick] {sin(deg(x))};
723 \addlegendentry{$f(x)=\sin(x)$}
724 \end{axis}
725 \end{tikzpicture}
726 \end{center}
727
728 \begin{center}
729 \begin{tikzpicture}
730 \begin{axis}[
731 axis lines = middle,
732 xlabel = {$x$},
733 ylabel = {$y$},
734 samples=100,
735 domain=-6:6,
736 legend pos=north east,
737 width=0.8\textwidth,
738 height=.4\textheight
739 ]
740 \addplot[blue, thick] {pi*sin(deg(x))};
741 \addlegendentry{$(g * \tilde{f})(t) = \pi \sin(t)$}
742 \end{axis}
743 \end{tikzpicture}
744 \end{center}
745
746Ainsi, 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.
747
748\subsection{Modélisation de la propagation des ondes sismiques}
749\subsubsection{Calcul du temps de propagation selon iasp91}
750On 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.
751
752\noindent\begin{minipage}{.5\textwidth}
753\centering
754\includegraphics[width=3cm, trim={0 0 0 0}, clip]{2025-03-26T23:10:40,393866488+01:00-side.png}
755 \captionof{figure}{TauPy}
756\end{minipage}%
757\begin{minipage}{.5\textwidth}
758\centering
759 \includegraphics[width=4cm, trim={0 0 0 0}, clip]{IASP91.png}
760 \captionof{figure}{iasp91}
761\end{minipage}
762
763
764\subsubsection{Tabulation et interpolation}
765\begin{center}
766 \includegraphics[width=10cm, trim={12cm 4cm 12cm 8cm}, clip]{2025-03-24T00:28:53,070002973+01:00.png}\\
767 \captionof{figure}{Visualization des deux tables précalculées}
768\end{center}
769\subsection{Multilatération}
770Une 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.
771
772Sauf 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.
773
774Ainsi, on va utiliser une méthode pour trouver le minimum sur une fonction d'erreur.
775
776\subsubsection{Fonction d'erreur}
777Voici notre fonction d'erreur:
778\vspace*{7em}
779
780\begin{equation*}
781 \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}}
782 \end{equation*}\vspace*{3em}
783 \begin{equation*}
784 +\eqnmarkbox[pink]{p10}{\left(\frac{(obs_iP-epoch)-\text{P}(depth, \text{greatCircleAngle}(lat,lon,lat_i,lon_i))}{obs_iP-epoch}\right)^2}
785\end{equation*}
786\annotate[yshift=-0.3em]{below}{p1}{Profondeur estimée}
787\annotatetwo[yshift=5.5em]{above}{p2}{p3}{Latitude et longitude estimée}
788\annotate[yshift=4.5em]{above}{p4}{Date de début du séisme estimée}
789\annotate[yshift=-0.4em]{below}{p5}{Tableau des observations}
790\annotate[yshift=2em]{above}{p6}{Temps de propagation avec la date de début du séisme estimée}
791\annotate[yshift=5em]{above,left}{p7}{Temps de propagation que l'on devrait observer compte tenu du modèle}
792\annotate[yshift=6.5em]{above,left}{p8}{Calcul de l'angle entre le seismographe et la position estimée du séisme}
793\annotate[yshift=-1em]{below,left}{p9}{Renormalisation pour éviter de favoriser les temps de propagation longs}
794\annotate[yshift=3em]{above,left}{p11}{On fait la moyenne quadratique pour avoir l'écart}
795\annotate[yshift=0em]{below,left}{p10}{Idem mais pour l'onde P}
796\vspace*{3em}
797
798\begin{center}
799 \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-31T10:14:11,496825900+02:00.png}\\
800 \captionof{figure}{Implémentation de la fonction d'erreur}
801 \end{center}
802\subsubsection{Méthode de Nelder-Mead}
803La 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.
804
805Elle 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).
806
807L'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.
808
809\begin{center}
810\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}
811 \captionof{figure}{Une itération de Nelder-Mead sur un espace de dimension 2}
812\end{center}
813
814\subsection{Magnitude sismique}
815On utilise l'echelle logarithmique de Richter pour quantifier la magnitude d'un séisme pour laquelle on a la formule suivante:
816\[
817M_\mathrm{L} = \log_{10} \left[ \frac{A}{A_\mathrm{0}(\delta)} \right]
818\]
819
820\noindent
821où $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 correrlation empirique).
822
823
824On utilisera la formule empirique de Tsuboi (Université de Tokyo):
825\[
826M_{\mathrm{L}} = \log_{10} A + 1.73 \log_{10} \Delta - 0.83
827\]
828
829\noindent
830où \( A \) est l'amplitude en micromètres et \( \Delta \) est la distance en kilomètres.
831
832Les différentes valeurs obtenues permettent de quantifier les dégats provoqués par les séismes:
833
834\begin{center}
835\begin{tabular}{|>{\columncolor{white}}l|l|l|p{5cm}|}
836\hline
837\rowcolor{gray!30}
838\textbf{Magnitude} & \textbf{Description} & \textbf{MMI Typique} & \textbf{Effets Moyens du Séisme} \\
839\hline
840\cellcolor{mmi1}1.0 - 1.9 & Micro & I & Micro-séismes, non ressentis. Enregistrés par les sismographes. \\
841\hline
842\cellcolor{mmi1}2.0 - 2.9 & Mineur & I & Légèrement ressenti par certaines personnes. Aucun dommage aux bâtiments. \\
843\hline
844\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. \\
845\hline
846\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. \\
847\hline
848\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. \\
849\hline
850\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. \\
851\hline
852\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. \\
853\hline
854\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. \\
855\hline
856\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. \\
857\hline
858\end{tabular}
859\end{center}
860
861\section{Résultats}
862\subsection{Tests}
863Je n'ai pas eu le temps d'implémenter la partie qui mesure automatiquement les temps ni la magnitude. Cependant, j'ai fait des mesures manuelles que j'ai donné à l'algorithme de multilatération.
864
865On constate que les positions prédites sont en accord avec celles données par l'USGS (United States Geological Survey)
866
867\begin{center}
868\includegraphics[width=18cm, trim={30cm 1cm 30cm 1cm}, clip]{debugMaps/rs2025flrsdg.png}
869 \captionof{figure}{rs2025flrsdg - Mon estimation (avec les données RaspberryShake)}
870
871 \includegraphics[width=12cm, trim={0 0 0 0}, clip]{debugMaps/rs2025fwmrzv.png}
872 \captionof{figure}{rs2025fwmrzv - Mon estimation (avec les données RaspberryShake)}
873
874
875 \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-25T21:07:44,659295790+01:00-cropped.png}
876 \captionof{figure}{rs2025fwmrzv - Estimation de l'USGS}
877\end{center}
878\newpage
879\section{Bibliographie}
880\begin{thebibliography}{9}
881
882\bibitem{convolution_video}
8833Blue1Brown,
884\textit{But what is a convolution?},
885YouTube video, 2022.\\
886\url{https://www.youtube.com/watch?v=KuXjwB4LzSA}
887
888\bibitem{fft_algorithm_video}
889Michael Pound,
890\textit{The Fast Fourier Transform Algorithm},
891YouTube video, 2023.\\
892\url{https://www.youtube.com/watch?v=toj_IoCQE-4}
893
894\bibitem{fft_matrix_factorizations}
895Charles Van Loan,
896\textit{The FFT Via Matrix Factorizations},
897Lecture notes, 2010.\\
898\url{https://www.cs.cornell.edu/~bindel/class/cs5220-s10/slides/FFT.pdf}
899
900\bibitem{eq_sources}
901R.J. Mitchell,
902\textit{Earthquake Sources},
903Lecture notes, Western Washington University.\\
904\url{https://www.geol.wwu.edu/rjmitch/L4_EQsources.pdf}
905
906\bibitem{seismic_waves}
907University of Hawaii,
908\textit{Compare, Contrast, and Connect: Seismic Waves and Determining Earth's Structure}.\\
909\url{https://manoa.hawaii.edu/exploringourfluidearth/physical/ocean-floor/layers-earth/compare-contrast-connect-seismic-waves-and-determining-earth-s-structure}
910
911\bibitem{gfz-1-d-earth-models}
912Peter Bormann,
913\textit{Global 1-D Earth models (IASP91 tables)},
914GFZ German Research Centre for Geosciences.\\
915\url{https://gfzpublic.gfz-potsdam.de/rest/items/item_4031/component/file_4032/content}
916
917\bibitem{earthquake_data_centers}
918Yacine Boussoufa,
919\textit{Earthquake Data Centers},
920GitHub repository.\\
921\url{https://github.com/YacineBoussoufa/EarthquakeDataCenters}
922
923\bibitem{cfcs_lecture}
924School of Informatics, University of Edinburgh,
925\textit{Computational Foundations of Cognitive Science: Lecture 15 (Convolutions and Kernels)}.\\
926\url{https://www.inf.ed.ac.uk/teaching/courses/cfcs1/lectures/cfcs_l15.pdf}
927
928\bibitem{penn_state_anova}
929Penn State Eberly College of Science,
930\textit{STAT 510: Lesson 8.2 - Cross Correlation Functions and Lagged Regressions},
931Online course material.\\
932\url{https://online.stat.psu.edu/stat510/lesson/8/8.2}
933
934\bibitem{nelder_mead}
935Jason Cantarella,
936\textit{Nelder-Mead Method},
937Lecture notes.\\
938\url{https://jasoncantarella.com/downloads/NelderMeadProof.pdf}
939
940\bibitem{broadband_magnitude}
941Tatsuhiko Hara,
942\textit{Determination of Broadband Moment Magnitude},
943IISEE/BRI.\\
944\url{https://iisee.kenken.go.jp/lna/download.php?f=2011082925678c01.pdf&n=T0-100-2007_Mwp-2-new.pdf&cid=T0-100-2007}
945
946\bibitem{obspy}
947Moritz Beyreuther, Robert Barsch, Lion Krischer, Tobias Megies, Yannik Behr and Joachim Wassermann,\\
948\textit{ObsPy: A Python Toolbox for Seismology},\\
949Seismological Research Letters, vol. 81, no. 3, pp. 530--533, 2010.\\
950doi:10.1785/gssrl.81.3.530
951
952\end{thebibliography}
953
954\end{document}