\documentclass[a4paper, 11pt]{article} \usepackage[a4paper, total={18cm, 24cm}]{geometry} \usepackage[T1]{fontenc} \usepackage[utf8]{inputenc} \usepackage[french]{babel} \usepackage{lmodern} \usepackage{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsthm} \usepackage{graphicx} \usepackage{color} \usepackage{xcolor} \usepackage{url} \usepackage{theorem} \usepackage{textcomp} \usepackage{listings} \usepackage{hyperref} %\usepackage{glossaries} \usepackage{parskip} \usepackage{amsmath, amsthm, amssymb} \usepackage{stmaryrd} \usepackage{graphicx} \usepackage{subfig} \usepackage{color} \usepackage{longtable} \usepackage{pgfplots} \usepackage{nicematrix} \usepackage[table]{xcolor} \usepackage{assets/texpackages/annotate-equations} \newenvironment{graytext}{\color{gray}}{\ignorespacesafterend} \graphicspath{ {./assets/} } % Define custom colors \definecolor{mmi1}{HTML}{FFFFFF} \definecolor{mmi2}{HTML}{BFCCFF} \definecolor{mmi3}{HTML}{A0E6FF} \definecolor{mmi4}{HTML}{80FFFF} \definecolor{mmi5}{HTML}{7AFF93} \definecolor{mmi6}{HTML}{FFFF00} \definecolor{mmi7}{HTML}{FFC800} \definecolor{mmi8}{HTML}{FF9100} \definecolor{mmi9}{HTML}{FF0000} \definecolor{mmi10}{HTML}{C80000} \definecolor{mmi11}{HTML}{A40000} \definecolor{mmi12}{HTML}{800000} \begin{document} \title{Traitement de signaux pour la détection de séismes et leur multilatération} \author{Dalibard Louis} \date{\today} \maketitle \tableofcontents \newpage \section{Séismes} \subsection{Introduction} Un séisme ou tremblement de terre est une secousse du sol résultant de la libération brusque d'énergie accumulée par les contraintes 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.\\ \includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{eq-ed-fault-labeled.png} \subsection{Ondes P et S} Les ondes sismiques, ou ondes élastiques, sont des mouvements vibratoires qui se propagent à travers un milieu matériel et peuvent le modifier irréversiblement si leur amplitude est suffisante. Elles sont engendrées par un événement initial, généralement un séisme. L'impulsion de départ déplace des atomes, qui en poussent d'autres à leur tour avant de reprendre leur place, ces déplacements oscillatoires se propageant ensuite de proche en proche. Un séisme émet des ondes dans toutes les directions de l'espace. La physique des milieux élastiques est régie par l'équation de Navier, qui implique l'existence de deux grands types d'ondes, mises en évidence expérimentalement : les ondes de volume qui traversent l'intérieur de la Terre et les ondes de surface qui se propagent dans une couche d'épaisseur limitée en suivant la surface terrestre. Les ondes de volume se subdivisent par ailleurs en ondes longitudinales (ondes P) et transversales (ondes S). \noindent\begin{minipage}{.5\textwidth} \centering \includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_compression_impulsion_1d_30_petit.png} \captionof{figure}{Ondes P} \end{minipage}% \begin{minipage}{.5\textwidth} \centering \includegraphics[width=3.5cm, trim={0 0 0 0}, clip]{Onde_cisaillement_impulsion_1d_30_petit.png} \captionof{figure}{Ondes S} \end{minipage} \subsection{Différence de vitesse de propagation} Les ondes P ou ondes primaires, appelées aussi ondes de compression (en anglais Pressure waves) ou ondes longitudinales. Le déplacement du sol qui accompagne leur passage se fait par des dilatations et des compressions successives. Ces déplacements du sol sont parallèles à la direction de propagation de l'onde. Elles se propagent dans tous les milieux et sont les plus rapides (autour de 6 km/s près de la surface), parcourant le chemin le plus court, même par noyau terrestre, et sont donc les premières enregistrées sur les sismogrammes. Elles sont responsables du grondement sourd que l'on peut entendre au début d'un tremblement de terre. Les ondes S ou ondes secondaires, appelées aussi ondes de cisaillement (Shear waves) ou ondes transversales. À leur passage, les mouvements du sol s'effectuent perpendiculairement au sens de propagation de l'onde. Ces ondes ne se propagent pas dans les milieux liquides, elles sont en particulier arrêtées par le noyau externe de la Terre. Leur vitesse est d'environ 4 km/s. Elles apparaissent en deuxième sur les sismogrammes. \begin{center} \includegraphics[width=7cm, trim={0 12.5mm 0cm 2cm}, clip]{2025-03-26T23:29:53,572838335+01:00.png}\\ \captionof{figure}{6 km/s (ondes P) vs 4 km/s (ondes S)} \end{center} On peut donc exploiter les temps d'arrivée de ces différentes ondes pour localiser un séisme. \section{Théorie} \subsection{Principe} Différentes étapes: \begin{enumerate} \item Acquisition de données en temps réel (SeedLink) \item Reconnaissance d'un séisme et mesure automatique des temps \item Calcul de la position et de la magnitude \begin{enumerate} \item Modélisation de la propagation des ondes sismiques \item Méthode numérique d'optimisation de fonction à plusieurs variables pour la multilatération \item Calcul de la magnitude \end{enumerate} \end{enumerate} \subsection{DSP} \subsubsection{Acquisition des données} On utilise le protocole Seedlink pour récupérer les données des sismographes en temps réel. \renewcommand{\arraystretch}{1.2} \begin{longtable}{|l|l|} \hline \textbf{Name} & \textbf{Host} \\ \hline AusPass & auspass.edu.au \\ BGR & eida.bgr.de \\ CISMID & www.cismid.uni.edu.pe \\ ENS & ephesite.ens.fr \\ GEOFON, GFZ & geofon.gfz-potsdam.de \\ GEONET & link.geonet.org.nz \\ Geoscience Australia & seis-pub.ga.gov.au \\ GSRAS (?) & 89.22.182.133 \\ Helsinki & finseis.seismo.helsinki.fi \\ Haiti & ayiti.unice.fr \\ ICGC & ws.icgc.cat \\ IDA Project & rtserve.ida.ucsd.edu \\ IFZ & data.ifz.ru \\ IPGP & rtserver.ipgp.fr \\ IRIS DMC & rtserve.iris.washington.edu \\ IRIS Jamaseis & jamaseis.iris.edu \\ ISNET - UNINA & 185.15.171.86 \\ LMU & erde.geophysik.uni-muenchen.de \\ NIGGG & 195.96.231.100 \\ NRCAN & earthquakescanada.nrcan.gc.ca \\ OBSEBRE & obsebre.es \\ OGS & nam.ogs.it \\ Oklahoma University & rtserve.ou.edu \\ ORFEUS & eida.orfeus-eu.org \\ PLSN (IGF Poland) & hudson.igf.edu.pl \\ Red Sìsmica de Puerto Rico & 161.35.236.45 \\ Red Sìsmica Baru & helis.redsismicabaru.com \\ RESIF & rtserve.resif.fr \\ SANET & 147.213.113.73 \\ RSIS & rsis1.on.br \\ SCSN-USC (South Carolina Seismic Network) & eeyore.seis.sc.edu:6382 \\ Seisme IRD & rtserve.ird.nc \\ Staneo & vibrato.staneo.fr \\ SNAC NOA & snac.gein.noa.gr \\ TexNet & rtserve.beg.utexas.edu \\ Thai Meteorological Department & 119.46.126.38 \\ UFRN (Universidade Federal do Rio Grande do Norte) & sislink.geofisica.ufrn.br \\ Unical Universita Della Calabria & www.sismocal.org \\ UNITS Università degli studi di Trieste & rtweb.units.it \\ UNIV-AG Université des Antilles & seedsrv0.ovmp.martinique.univ-ag.fr \\ Universidade de Évora & clv-cge.uevora.pt \\ Universidad de Colima & 148.213.24.15 \\ UPR & worm.uprm.edu \\ USGS & cwbpub.cr.usgs.gov \\ USP-IAG & seisrequest.iag.usp.br \\ \hline \end{longtable} \subsubsection{Extraction des temps d'arrivée} \begin{center} \includegraphics[width=15cm, trim={0 0 0 0}, clip]{R0ED0.EHZ-1743142835749.9944-cropped.png}\\ \captionof{figure}{Exemple d'un enregistrement de sismographe} \end{center} On a un signal comme ça, on peut facilement voir l'onde P et l'onde S. Mais comment faire pour extraire ces temps d'arrivée automatiquement ? Notamment quand on aura des centaines de mesures venant de partout dans le monde ? Pour cela on va faire un détour par des maths. \subsubsection{Convolution} Si on a deux fonctions, quelles opérations peut-on faire pour avoir une nouvelle fonction? La plus simple semble être les additionner. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] % Define the functions \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$}; \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$}; \addplot[green, thick] {sin(deg(x)) + cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) + g(x)$}; % Add legend \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) + g(x)$} \end{axis} \end{tikzpicture} \end{center} Une autre façon serait de les multiplier. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] % Define the functions \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$}; \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$}; \addplot[green, thick] {sin(deg(x)) * cos(deg(2*x))} node[right] {\footnotesize $h(x) = f(x) \cdot g(x)$}; % Add legend \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = f(x) \cdot g(x)$} \end{axis} \end{tikzpicture} \end{center} Mais il existe aussi une autre opération, la convolution. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] % Define the functions \addplot[blue, thick] {sin(deg(x))} node[right] {\footnotesize $f(x) = \sin x$}; \addplot[red, thick] {cos(deg(2*x))} node[right] {\footnotesize $g(x) = \cos 2x$}; \addplot[green, thick] {0} node[right] {\footnotesize $h(x) = (f*g)(x)$}; % Add legend \legend{$f(x) = \sin x$, $g(x) = \cos 2x$, $h(x) = (f*g)(x) = 0$} \end{axis} \end{tikzpicture} \end{center} Qui a de nombreuses applications dans le traitement d'images, la théorie de la probabilité, dans la solution d'équations différencielles (température), et comme nous le verrons (aussi la multiplication de polynômes). On définit formellement la convolution comme un produit sur l'espace des fonctions intégrables. Soit $f$ et $g$ deux fonctions intégrables sur $\mathbb{R}$. Pour tout $t \in \mathbb{R}$, \[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \] Dans le cas discret on utilisera une somme. Soit $f$ et $g$ deux fonctions de $\mathbb{Z}$ dans $\mathbb{C}$. Pour tout $n \in \mathbb{Z}$, \[ (f*g)[n]=\sum_{m=-\infty}^{+\infty} f[m] \cdot g[n-m] \] et dans le cas de fonction périodiques, il suffit de faire l'intégrale sur une période. Pour tout $t \in \mathbb{R}$, \[ (f*g)(t)=\int_{0}^{T} f(x) \cdot g(t-x) dx \] (on peut également généraliser cela à des fonctions discretès périodiques) \subsubsection{Propriétés algébriques de la convolution} \begin{itemize} \item Commutatif On remarquera que si \[ (f*g)(t)=\int_{-\infty}^{+\infty} f(x) \cdot g(t-x) dx \] Et on fait le changement de variable $u=t-x$ On a \[ (f*g)(t)=\int_{+\infty}^{-\infty} f(t-u) \cdot g(u) -du=\int_{-\infty}^{+\infty} f(t-u) \cdot g(u) du=(g*f)(t)\] \item Associatif \[ (f*g)*h= f*(g*h) \] C'est une conséquence du théorème de Fubini. \item Distributif \[ f*(g+h)= f*g+f*h \] Par linéarité de l'intégrale. \end{itemize} L'espace des fonctions intégrables muni de $*$ forme un demi-groupe commutatif (car pas d'élement neutre). Pourquoi une telle construction ? Déjà à quoi ça ressemble? Si on voulait compter le nombre de combinaisons de dés qui donnent un certain total, on les disposerait en liste, on metterait une d'elle à l'envers et on regarderait le nombre de dés qui s'alignent à chaque fois. \begin{center} \begin{tikzpicture}[scale=0.5] % Define Dice Faces \def\diceA{1,2,3,4,5,6} % Top Die (Static) \def\diceB{6,5,4,3,2,1} % Bottom Die (Sliding) % Define step spacing \def\stepSpace{3} % Loop through each step \foreach \step in {-5,-4,-3,-2,-1,0,1,2,3,4,5} { \pgfmathsetmacro{\sommedonnant}{int(7 + \step)} % Label step number \node[anchor=east] at (-0.5, -\step*3) {Somme donnant \sommedonnant:}; % Top Row - Fixed Dice \foreach \x [count=\i] in \diceA { \node[draw, minimum size=0.5cm] at (\i, -\step*3) {\x}; } % Bottom Row - Sliding Dice \foreach \x [count=\i] in \diceB { \pgfmathsetmacro{\sumval}{int(\i + \step)} \ifnum \sumval > 0 \ifnum \sumval < 7 \node[draw=red, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x}; \else \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x}; \fi \else \node[draw, minimum size=0.5cm] at (\i+\step, -\step*3-1) {\x}; \fi } } \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ xlabel={Somme}, ylabel={Fréquence}, ymin=0, ymax=7, xmin=1.5, xmax=12.5, xtick={2,3,4,5,6,7,8,9,10,11,12}, ytick={1,2,3,4,5,6,7}, area style, width=12cm, height=8cm, major grid style={line width=.2pt,draw=gray!50}, grid=major, bar width=0.8, title={Distribution des sommes}, nodes near coords, ] \addplot[fill=blue!40, draw=blue!80, opacity=0.8] coordinates { (2,1) (3,2) (4,3) (5,4) (6,5) (7,6) (8,5) (9,4) (10,3) (11,2) (12,1) }; \end{axis} \end{tikzpicture} \end{center} On obtient les fréquences d'apparition des différentes combinaisons. Nous pouvons également remarquer que si à la place de rajouter $1$ à chaque fois, on ponderait les probabilités et faisait le produit des probabilités qui se trouvent face à face, on aurait également le résultat pour le cas non équiprobable. On vient de faire la convolution entre cette fonction et elle même. \begin{center} \begin{tikzpicture} \begin{axis}[ xlabel={$x$}, ylabel={$f(x)$}, xmin=-1, xmax=8, ymin=-0.2, ymax=1.4, xtick={-1,0,1,2,3,4,5,6,7,8}, ytick={0,1}, width=10cm, height=6cm, grid=major, title={Fonction nulle puis constante égale à $1$ puis nulle}, samples=100, domain=-1:8, axis lines=middle, clip=false, ] % Plot the function \addplot[blue, thick, const plot] coordinates { (-1, 0) (0.999, 0) (1, 1) (6, 1) (6.001, 0) (8, 0) }; % Add markers to explicitly show the endpoints \addplot[only marks, mark=*, mark options={fill=red}] coordinates { (1, 1) (6, 1) }; \end{axis} \end{tikzpicture} \end{center} \begin{center} \includegraphics[width=14cm, trim={0 0 0 0}, clip]{2025-04-01T23:42:09,818697014+02:00.png} \captionof{figure}{numpy confirme ce résultat} \end{center} On remarque que si on le fait dans le cas discret, on a des valeurs en plus sur les bords que l'on peut ignorer. Faire la convolution revient donc à retourner une des listes des coefficients et le résultat est une fonction qui à un décalage $t$ entre les deux listes associe une somme des produits deux à deux des coefficients qui se retrouvent en face l'un de l'autre par ce décalage. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(2*x))}; \addlegendentry{$\cos(2x)$} \addplot[red, thick] {sin(deg(x))}; \addlegendentry{$\sin(x)$} \end{axis} \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(2*x))}; \addlegendentry{$\cos(2x)$} \addplot[red, thick] {sin(deg(-x))}; \addlegendentry{$\sin(-x)$} \end{axis} \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(2*x))*sin(deg(-x))}; \addlegendentry{$\cos(2x)\sin(-x)$} % Highlight positive areas in red \addplot [fill=red, opacity=0.3, domain=-2*pi:2*pi] {max(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle; % Highlight negative areas in blue \addplot [fill=blue, opacity=0.3, domain=-2*pi:2*pi] {min(0, cos(deg(2*x))*sin(deg(-x)))} \closedcycle; \end{axis} \end{tikzpicture} \end{center} On voit que c'est nul dans ce cas. \subsubsection{Moyennage} Si on applique la fonction $\frac{f}{6}$ à une fonction $g$ quelconque on observe une sorte de moyennage qui se fait entre 6 valeurs proches. \begin{center} \begin{tikzpicture} \begin{axis}[ xlabel={$x$}, ylabel={$y$}, xmin=-0.5, xmax=23.5, ymin=-0.5, ymax=9, xtick={0,2,4,6,8,10,12,14,16,18,20,22}, width=12cm, height=8cm, grid=major, title={Moyennage par convolution}, legend pos=north west, ] % Plot the first dataset with connected lines \addplot[blue, thick] coordinates { (0,0) (1,0) (2,0) (3,0) (4,1) (5,2) (6,3) (7,4) (8,5) (9,6) (10,7) (11,8) (12,0) (13,4) (14,2) (15,3) (16,8) (17,1) (18,2) (19,4) (20,8) (21,0) (22,0) (23,0) }; % Plot the second dataset with connected lines \addplot[red, thick] coordinates { (0,0) (1,0.16666667) (2,0.5) (3,1) (4,1.66666667) (5,2.5) (6,3.5) (7,4.5) (8,5.5) (9,5) (10,5) (11,4.5) (12,4) (13,4.16666667) (14,3) (15,3.33333333) (16,3.33333333) (17,4.33333333) (18,3.83333333) (19,2.5) (20,2.33333333) (21,2) (22,1.33333333) (23,0) }; \legend{$g$, $g*\frac{f}{6}$} \end{axis} \end{tikzpicture} \end{center} Ce processus est fait en deux dimensions sur des images avec differentes fonctions pour avoir des effets comme la détection de bord ou le floutage. \subsubsection{Transformée de Fourier discrète} Si on doit faire la convolution de deux listes $a$ et $b$ telles que $a=[a_0,a_1,...,a_{n-1}]$ et $b=[b_0,b_1,...,b_{k-1}]$, On a $(a*b)[p] = \sum\limits_{\substack{i\in\mathbb{Z} \\ 0\leq p-i\leq k-1\\ 0\leq i\leq n-1}} a_ib_{p-i}$ On reconnait la formule du produit des polynômes $P(X)=\sum_{i=0}^{n-1} a_i X^i$ et $Q(X)=\sum_{j=0}^{k-1} a_j X^j$ $(a*b)[p]$ est le coefficient du terme de degré $p$ dans le produit $PQ(X)=\sum_{j=0}^p(\sum\limits_{\substack{i\in\mathbb{Z} \\ 0\leq p-i\leq k-1\\ 0\leq i\leq n-1}} a_ib_{j-i})X^j$ Bon, de premier abord, ça n'a pas l'air de beaucoup nous aider, le produit de polynômes est $\mathcal{O}(n^2)$. Mais pour un polynôme, du fait de leur rigidité, on peut juste évaluer leur produit en un nombre fini de points pour retrouver les coefficients. Et donc l'idée si on avait un nombre discret de coefficients serait de les traiter comme les coefficients d'un polynôme, on les évalue un certain nombre de fois et on multiplie ensemble les évaluations, puis on retrouve les coefficients en résolvant le système de Vandermonde correspondant. Le problème c'est que cela semble prendre plus de temps, un pivot de Gauss est en $\mathcal{O}(n^3)$ Mais si on choisit d'évaluer nos valeurs en les racines de l'unité, on a un système plus simple à résoudre. On va pouvoir faire les opérations en $\mathcal{O}(n \cdot log(n))$ Les coefficients que l'on obtient quand on évalue en les $\omega_n^k = e^{-\frac{2ki\pi}{n}}$ valeurs sont appellés la transformation de Fourier discrète TFD (un équivalent discret de la transformation de Fourier continue). Il existe aussi d'autres algorithmes. \subsubsection{Radix-2 decimation-in-time (DIT) - Factorisation Cooley-Tukey} Evaluer notre polynôme $P(x)=\sum_{i=0}^{n-1} a_i x^i$ en les $\omega_n^k$ revient à faire la multiplication matricielle suivante: \[ R=\begin{bmatrix} f_0 \\ f_1 \\ \vdots \\ f_{n-1} \end{bmatrix}=\begin{bmatrix} P(\omega_n^0) \\ P(\omega_n^1) \\ \vdots \\ P(\omega_n^{n-1}) \end{bmatrix} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} \] On reconnaît la matrice de Vandermonde transposée associée à l'évaluation en les \( \omega_n^k \) On se restreint au cas où $n=2^p$ ($p \in \mathbb{N}$) pour simplifier l'explication même si en théorie ça ne change pas grand chose vu qu'il y a des algorithmes généraux possibles et on peut rajouter des 0 à la fin pour avoir une liste de taille $2^p$. On note la matrice de Vandermonde transposée, qui permet de calculer le DFT pour tous nos coefficients, $F_{2^p}=\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)(n-1)} \end{bmatrix}$ On peut permuter les lignes, à la fin, ce qu'on a c'est que, du fait des propriétés des racines de l'unité, si on note, $D_{2^{p-1}} = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & \omega_{2^{p-1}} & 0 & \cdots & 0 \\ 0 & 0 & \omega_{2^{p-1}}^2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \omega_{2^{p-1}}^{2^{p-1}-1} \end{bmatrix}$ est: \begin{equation*} R=F_{2^p}\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}= \begin{pNiceArray}{cc|cc} \eqnmarkbox[blue]{mcoeffdiag1}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag2}{D_{2^{p-1}}} \\ \hline \eqnmarkbox[blue]{mcoeffdiag3}{I_{2^{p-1}}} && \eqnmarkbox[blue]{mcoeffdiag4}{-D_{2^{p-1}}} \end{pNiceArray} \begin{pNiceArray}{cc|cc} \eqnmarkbox[cyan]{mrecurrence1}{F_{2^{p-1}}} && \mathbf{0} \\ \hline \mathbf{0} && \eqnmarkbox[cyan]{mrecurrence2}{F_{2^{p-1}}} \end{pNiceArray} \begin{bmatrix} \eqnmarkbox[green]{mcoeffpairs1}{a_0} \\ \eqnmarkbox[green]{mcoeffpairs2}{a_2} \\ \vdots \\ \eqnmarkbox[green]{mcoeffpairs3}{a_{2^{p-1}-2}} \\ \eqnmarkbox[green]{mcoeffpairs4}{a_{2^{p-1}}} \\ \eqnmarkbox[magenta]{mcoeffimpairs1}{a_{1}} \\ \eqnmarkbox[magenta]{mcoeffimpairs2}{a_{3}} \\ \vdots \\ \eqnmarkbox[magenta]{mcoeffimpairs3}{a_{n-3}} \\ \eqnmarkbox[magenta]{mcoeffimpairs4}{a_{n-1}} \end{bmatrix} \end{equation*} \annotatetwo[yshift=1.5em]{above}{mcoeffdiag1}{mcoeffdiag2}{Blocs diagonaux, de l'ordre de $\mathcal{O}(n)$ opérations} \annotate[yshift=0.5em]{above}{mcoeffpairs1}{Coefficients d'indice pair} \annotate[yshift=0em]{below}{mcoeffimpairs4}{Coefficients d'indice impair} \annotate[yshift=0em]{below,left}{mrecurrence2}{Récurrence} On évalue la compléxité de l'algorithme.\\ On note $u_p$ sa complexité en fonction de $p$ et $C_n$ sa complexité en fonction de $n$. \begin{equation*} u_{p+1}=\eqnmarkbox[magenta]{cdiagmult}{A \cdot 2^{p+1}}+\eqnmarkbox[red]{crec}{2u_{p}} \end{equation*} \annotate[yshift=0.5em]{above,left}{cdiagmult}{Compléxité du produit sur la diagonale} \annotate[yshift=0em]{below,right}{crec}{Traitement des coefficients par récurrence} On factorise par la solution homogène, $\frac{u_{p+1}}{2^{p+1}}=A+\frac{u_{p}}{2^p}$ $\frac{u_{p}}{2^{p}}=u_{0}+A\cdot p$ $u_{p}=u_{0}\cdot2^{p}+A\cdot p \cdot 2^{p}$ Or $p=\log_2{n}$ Donc $C_n = u_{\log_2{n}} = u_{0} \cdot n+A \cdot \log_2{n} \cdot n = \mathcal{O}(n \log n)$ \subsubsection{IFFT} On peut montrer que: $F_{2^p}^{-1}=\frac{1}{2^p}\begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & \overline{\omega_n}^1 & \overline{\omega_n}^2 & \cdots & \overline{\omega_n}^{n-1} \\ 1 & \overline{\omega_n}^2 & \overline{\omega_n}^4 & \cdots & \overline{\omega_n}^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \overline{\omega_n}^{n-1} & \overline{\omega_n}^{2(n-1)} & \cdots & \overline{\omega_n}^{(n-1)(n-1)} \end{bmatrix}$ $A=\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix}=F_{2^p}^{-1}R=F_{2^p}^{-1}\begin{bmatrix} f_0 \\ f_1 \\ \vdots \\ f_{n-1} \end{bmatrix}$ Le conjugué passe au produit et à la somme, donc aussi pour les matrices (prendre le conjugué d'une matrice c'est prendre le conjugué des termes de la matrice). $A=\overline{\overline{A}}=\overline{\overline{F_{2^p}^{-1}R}}=\overline{\overline{F_{2^p}^{-1}}\overline{R}}=\frac{1}{2^p}\overline{F_{2^p}\overline{R}}$ On peut donc utiliser la même technique, en prenant le conjugué avant d'appliquer un FFT et en le prenant après puis en renormalisant. \subsubsection{Passage du domaine temporel au domaine fréquentiel} Sans passer trop de temps dessus, le fait de multiplier par ces coefficients spécifiques, revient à décomposer en ondes sinusoidales de différentes fréquences et phases notre signal. Il y a beaucoup d'applications dans d'autres domaines. Mais ce qui est remarquable c'est que pour calculer plus rapidement notre résultat. On est passé par une autre réprésentation. Maintenant qu'on peut calculer des convolutions rapidement, revenons à notre problème. % TODO: https://fr.wikipedia.org/wiki/Transformation_de_Fourier_rapide \subsubsection{Corrélation croisée} % TODO: https://en.wikipedia.org/wiki/Cross-correlation Si on définit $\tilde{f}(t)=f(-t)$ La corrélation croisée de \(f\) et \(g\) est \( (g * \tilde{f})(t) = \int_{-\infty}^{\infty} \overline{f(x-t)} g(x) \, dx \) Ou sur une seule période pour une fonction périodique. On va donc faire glisser les fonctions dans le même sens et multiplier les valeurs de la fonction évaluée en ses points ensembles deux à deux puis faire la somme. Quand les fonctions sont similaires et bien alignées, alors, on multiplie les grandes valeurs ensemble et les petites valeurs ensemble et donc on a résultat grand, alors que quand elles ne le sont pas, on multiplie les grandes valeurs avec les petites valeurs et on aura un résultat plus petit. On prend le conjugué pour que si c'est des fonctions complexes, que quand elles s'alignent elles donnent une plus grande valeur au lieu d'une plus petite à cause de la valeur imaginaire. Et donc on ne regardera que la partie réele et à quel point elle est grande. Donc si on a une fonction $f$ et une fonction $g$ où $g$ est dans $f$ mais avec un certain retard $\tau$ alors $(g * \tilde{f})$ sera maximale en $\tau$. \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {cos(deg(x))}; \addlegendentry{$g(x)=\cos(x)$} \addplot[red, thick] {sin(deg(x))}; \addlegendentry{$f(x)=\sin(x)$} \end{axis} \end{tikzpicture} \end{center} \begin{center} \begin{tikzpicture} \begin{axis}[ axis lines = middle, xlabel = {$x$}, ylabel = {$y$}, samples=100, domain=-6:6, legend pos=north east, width=0.8\textwidth, height=.4\textheight ] \addplot[blue, thick] {pi*sin(deg(x))}; \addlegendentry{$(g * \tilde{f})(t) = \pi \sin(t)$} \end{axis} \end{tikzpicture} \end{center} Ainsi, si on a un enregistrement de référence pour nos différentes ondes $P$ et $S$ alors on peut utiliser la corrélation croisée entre eux et le signal en temps réel pour détecter ces ondes et mesurer avec exactitude le temps d'arrivée de ces ondes. \subsection{Modélisation de la propagation des ondes sismiques} \subsubsection{Calcul du temps de propagation selon iasp91} On va utiliser un modèle de vitesse de propagation des différentes ondes sismiques selon la profondeur et tracer des rayons pour calculer le temps de propagation. On remarque que la courbure qu'on obtient est du à la réfraction. L'onde va emprunter le chemin le plus court, même si cela veut dire aller plus en profondeur d'abord. \noindent\begin{minipage}{.5\textwidth} \centering \includegraphics[width=3cm, trim={0 0 0 0}, clip]{2025-03-26T23:10:40,393866488+01:00-side.png} \captionof{figure}{TauPy} \end{minipage}% \begin{minipage}{.5\textwidth} \centering \includegraphics[width=4cm, trim={0 0 0 0}, clip]{IASP91.png} \captionof{figure}{iasp91} \end{minipage} \subsubsection{Tabulation et interpolation} \begin{center} \includegraphics[width=10cm, trim={12cm 4cm 12cm 8cm}, clip]{2025-03-24T00:28:53,070002973+01:00.png}\\ \captionof{figure}{Visualization des deux tables précalculées} \end{center} \subsection{Multilatération} Une fois qu'on a une table des différentes vitesses de propagation pour les différentes profondeurs et les ondes P et S ainsi que nos mesures, on peut résoudre l'équation avec nos mesures. Sauf qu'on a plusieurs mesures donc on veut une moyenne pour nos mesures pour avoir le meilleur résultat et nos équations ne sont pas linéaires comme vous pouvez le constater sur la diapositive précédente. Ainsi, on va utiliser une méthode pour trouver le minimum sur une fonction d'erreur. \subsubsection{Fonction d'erreur} Voici notre fonction d'erreur: \vspace*{7em} \begin{equation*} \text{E}(\eqnmarkbox[blue]{p1}{depth},\eqnmarkbox[blue]{p2}{lat},\eqnmarkbox[blue]{p3}{lon},\eqnmarkbox[blue]{p4}{epoch},\eqnmarkbox[green]{p5}{obs})=\sum_i \left(\frac{\eqnmarkbox[green]{p6}{(obs_iS-epoch)}-\eqnmarkbox[blue]{p7}{\text{S}(depth, \eqnmark[red]{p8}{\text{greatCircleAngle}(lat,lon,lat_i,lon_i))}}}{\eqnmarkbox[green]{p9}{obs_iS-epoch}}\right)^{\eqnmarkbox[magenta]{p11}{2}} \end{equation*}\vspace*{3em} \begin{equation*} +\eqnmarkbox[pink]{p10}{\left(\frac{(obs_iP-epoch)-\text{P}(depth, \text{greatCircleAngle}(lat,lon,lat_i,lon_i))}{obs_iP-epoch}\right)^2} \end{equation*} \annotate[yshift=-0.3em]{below}{p1}{Profondeur estimée} \annotatetwo[yshift=5.5em]{above}{p2}{p3}{Latitude et longitude estimée} \annotate[yshift=4.5em]{above}{p4}{Date de début du séisme estimée} \annotate[yshift=-0.4em]{below}{p5}{Tableau des observations} \annotate[yshift=2em]{above}{p6}{Temps de propagation avec la date de début du séisme estimée} \annotate[yshift=5em]{above,left}{p7}{Temps de propagation que l'on devrait observer compte tenu du modèle} \annotate[yshift=6.5em]{above,left}{p8}{Calcul de l'angle entre le seismographe et la position estimée du séisme} \annotate[yshift=-1em]{below,left}{p9}{Renormalisation pour éviter de favoriser les temps de propagation longs} \annotate[yshift=3em]{above,left}{p11}{On fait la moyenne quadratique pour avoir l'écart} \annotate[yshift=0em]{below,left}{p10}{Idem mais pour l'onde P} \vspace*{3em} \begin{center} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-31T10:14:11,496825900+02:00.png}\\ \captionof{figure}{Implémentation de la fonction d'erreur} \end{center} \subsubsection{Méthode de Nelder-Mead} La méthode d'optimisation que j'utilise est celle de Nelder-Mead. Elle permet de trouver un minimum local dans une fonction continue avec une recherche directe heuristique. Elle utilise ce qu'on appelle des simplexes, c'est la généralisation des triangles à plus de dimensions, en dimension 1 c'est une ligne, en dimension $3$ ça correspond à un tetraèdre, à chaque fois en dimension $n$, il $n+1$ côtés (il y a $4$ dimensions ici donc il faudra $n+1=5$ cotés au simplexe). L'approche la plus simple consiste à remplacer le plus mauvais point par un point reflétant le centroïde des $n$ points restants. Si ce point est meilleur que le meilleur point actuel, nous pouvons essayer de l'étirer de manière exponentielle le long de cette ligne. En revanche, si ce nouveau point n'est pas beaucoup mieux que la valeur précédente, nous traversons une vallée et nous rétrécissons le simplexe vers un meilleur point. \begin{center} \includegraphics[width=18cm, trim={0 0 0 0}, clip]{neldermead/An-iteration-of-the-Nelder-Mead-method-over-two-dimensional-space-showing-point-p-min.png} \captionof{figure}{Une itération de Nelder-Mead sur un espace de dimension 2} \end{center} \subsection{Magnitude sismique} On utilise l'echelle logarithmique de Richter pour quantifier la magnitude d'un séisme pour laquelle on a la formule suivante: \[ M_\mathrm{L} = \log_{10} \left[ \frac{A}{A_\mathrm{0}(\delta)} \right] \] \noindent où $A$ correspond à l'amplitude maximale mesurée (en m) par le sismographe et $A_\mathrm{0}(\delta)$ un coefficient de correction qui dépend de la distance ($\delta$) à l'epicentre et dont le calcul diffère selon les modèles employés (généralement on utilise une table de correrlation empirique). On utilisera la formule empirique de Tsuboi (Université de Tokyo): \[ M_{\mathrm{L}} = \log_{10} A + 1.73 \log_{10} \Delta - 0.83 \] \noindent où \( A \) est l'amplitude en micromètres et \( \Delta \) est la distance en kilomètres. Les différentes valeurs obtenues permettent de quantifier les dégats provoqués par les séismes: \begin{center} \begin{tabular}{|>{\columncolor{white}}l|l|l|p{5cm}|} \hline \rowcolor{gray!30} \textbf{Magnitude} & \textbf{Description} & \textbf{MMI Typique} & \textbf{Effets Moyens du Séisme} \\ \hline \cellcolor{mmi1}1.0 - 1.9 & Micro & I & Micro-séismes, non ressentis. Enregistrés par les sismographes. \\ \hline \cellcolor{mmi1}2.0 - 2.9 & Mineur & I & Légèrement ressenti par certaines personnes. Aucun dommage aux bâtiments. \\ \hline \cellcolor{mmi3}3.0 - 3.9 & Léger & II à III & Souvent ressenti, mais cause rarement des dégâts. Secousses perceptibles d’objets à l’intérieur. \\ \hline \cellcolor{mmi5}4.0 - 4.9 & Faible & IV à V & Secousses intérieures notables et bruits de cliquetis. Légèrement ressenti à l’extérieur. Dégâts minimes possibles. \\ \hline \cellcolor{mmi6}5.0 - 5.9 & Modéré & VI à VII & Peut endommager les bâtiments mal construits ; ressenti par tous. Peu ou pas de dégâts aux bâtiments solides. \\ \hline \cellcolor{mmi7}6.0 - 6.9 & Fort & VII à IX & Dégâts modérés aux structures solides ; dégâts sévères aux structures faibles. Ressenti sur de grandes régions. \\ \hline \cellcolor{mmi8}7.0 - 7.9 & Majeur & VIII ou plus & Dégâts majeurs et effondrements possibles. Dommages concentrés dans un rayon de 250 km. \\ \hline \cellcolor{mmi9}8.0 - 8.9 & Très fort & VIII+ & Destructions majeures à totales. Dommages sur des zones très vastes. Ressenti à très grande distance de l’épicentre. \\ \hline \cellcolor{mmi10}9.0 - 9.9 & Extrême & XII & Destruction quasi-totale, dégâts graves ou effondrement de tous les bâtiments. Modification du relief. \\ \hline \end{tabular} \end{center} \section{Résultats} \subsection{Tests} Je 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. On constate que les positions prédites sont en accord avec celles données par l'USGS (United States Geological Survey) \begin{center} \includegraphics[width=18cm, trim={30cm 1cm 30cm 1cm}, clip]{debugMaps/rs2025flrsdg.png} \captionof{figure}{rs2025flrsdg - Mon estimation (avec les données RaspberryShake)} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{debugMaps/rs2025fwmrzv.png} \captionof{figure}{rs2025fwmrzv - Mon estimation (avec les données RaspberryShake)} \includegraphics[width=12cm, trim={0 0 0 0}, clip]{2025-03-25T21:07:44,659295790+01:00-cropped.png} \captionof{figure}{rs2025fwmrzv - Estimation de l'USGS} \end{center} \newpage \section{Bibliographie} \begin{thebibliography}{9} \bibitem{convolution_video} 3Blue1Brown, \textit{But what is a convolution?}, YouTube video, 2022.\\ \url{https://www.youtube.com/watch?v=KuXjwB4LzSA} \bibitem{fft_algorithm_video} Michael Pound, \textit{The Fast Fourier Transform Algorithm}, YouTube video, 2023.\\ \url{https://www.youtube.com/watch?v=toj_IoCQE-4} \bibitem{fft_matrix_factorizations} Charles Van Loan, \textit{The FFT Via Matrix Factorizations}, Lecture notes, 2010.\\ \url{https://www.cs.cornell.edu/~bindel/class/cs5220-s10/slides/FFT.pdf} \bibitem{eq_sources} R.J. Mitchell, \textit{Earthquake Sources}, Lecture notes, Western Washington University.\\ \url{https://www.geol.wwu.edu/rjmitch/L4_EQsources.pdf} \bibitem{seismic_waves} University of Hawaii, \textit{Compare, Contrast, and Connect: Seismic Waves and Determining Earth's Structure}.\\ \url{https://manoa.hawaii.edu/exploringourfluidearth/physical/ocean-floor/layers-earth/compare-contrast-connect-seismic-waves-and-determining-earth-s-structure} \bibitem{gfz-1-d-earth-models} Peter Bormann, \textit{Global 1-D Earth models (IASP91 tables)}, GFZ German Research Centre for Geosciences.\\ \url{https://gfzpublic.gfz-potsdam.de/rest/items/item_4031/component/file_4032/content} \bibitem{earthquake_data_centers} Yacine Boussoufa, \textit{Earthquake Data Centers}, GitHub repository.\\ \url{https://github.com/YacineBoussoufa/EarthquakeDataCenters} \bibitem{cfcs_lecture} School of Informatics, University of Edinburgh, \textit{Computational Foundations of Cognitive Science: Lecture 15 (Convolutions and Kernels)}.\\ \url{https://www.inf.ed.ac.uk/teaching/courses/cfcs1/lectures/cfcs_l15.pdf} \bibitem{penn_state_anova} Penn State Eberly College of Science, \textit{STAT 510: Lesson 8.2 - Cross Correlation Functions and Lagged Regressions}, Online course material.\\ \url{https://online.stat.psu.edu/stat510/lesson/8/8.2} \bibitem{nelder_mead} Jason Cantarella, \textit{Nelder-Mead Method}, Lecture notes.\\ \url{https://jasoncantarella.com/downloads/NelderMeadProof.pdf} \bibitem{broadband_magnitude} Tatsuhiko Hara, \textit{Determination of Broadband Moment Magnitude}, IISEE/BRI.\\ \url{https://iisee.kenken.go.jp/lna/download.php?f=2011082925678c01.pdf&n=T0-100-2007_Mwp-2-new.pdf&cid=T0-100-2007} \end{thebibliography} \end{document}