Automatisation du calcul sur le serveur de L’ESGT 

Surcharge atmosphérique

L’atmosphère constitue le deuxième réservoir d’eau après l’océan. Les variations de pression atmosphérique causées par les mouvements de masses d’air provoquent alors les surcharges dites atmosphériques. Les déplacements qui en résultent sont les plus importants dans les zones où la latitude est importante et durant l’hiver. DACH et al. (2011) ont démontré que la zone la plus affectée par cette surcharge se trouve en Russie. Ainsi les régions des pôles sont celles qui subissent le plus de variabilité de pression durant une année, à l’inverse de l’équateur où l’énergie solaire reste à peu près constante durant l’année. Les déplacements de la croûte terrestre sont inversement proportionnels à la surcharge. Ainsi si la pression augmente, la croûte terrestre s’abaissera (PETROV et BOY, 2004). Cette surcharge entraine surtout un déplacement de la composante verticale. Théoriquement, les effets de la charge atmosphérique sont d’ordre sub-centimétrique (COLLILIEUX, 2008).

Surcharge hydrologique

Les effets de charge hydrologiques sont essentiellement dus à la masse d’eau présente dans les nappes phréatiques. Au cours du temps cette masse varie et par conséquent appuie plus ou moins sur la surface terrestre (COLLILIEUX, 2008). La variation du bilan entre les précipitations (pluie ou neige), le ruissellement et l’évapotranspiration, cause une redistribution de masse d’eau continentale, et se traduit ainsi par une flexure sensible de la croûte terrestre (BOEREZ, 2013).
Un autre facteur entraîne ces variations, il s’agit de l’énergie solaire. La période du cycle hydrologique est généralement d’un an. Ainsi les surcharges hydrologiques sont principalement annuelles. Les différents régimes climatiques présents sur le globe font que cette surcharge diffère aussi spatialement (montagnes, bassins fluviaux, …). L’humidité du sol et l’enneigement constituent les deux principaux réservoirs de l’eau continentale. Leur variation est avant tout saisonnière et génère des déformations de la croûte principalement verticales. Ces déformations sont de type élastique et peuvent atteindre 30 mm verticalement (BEVIS et al., 2005).

Surcharge océanique

Les effets de charge océanique sont en partie causés par les marées océaniques qui correspondent aux oscillations du niveau de la mer. Cependant ces effets sont corrigés dans les séries GPS, mais il existe des surcharges océaniques non maréales qui ne sont pas corrigées dans les séries et qui doivent être modélisées au mieux dans les modèles de surcharge. Les variations de masse d’eau induites par ces marées agissent sur l’écorce terrestre qui subit alors un mouvement de flexion. A l’inverse des stations côtières, les stations situées à quelques centaines de kilomètres des côtes ne subissent que très peu de déformations. Les plus forts déplacements verticaux ont été constatés par Van Dam et al. (2012) le long des côtes asiatiques, australiennes, nord et sud-américaines mais aussi le long de la mer du Nord et près des îles Aléoutiennes. Ce type de surcharge a une forte interaction avec la surcharge atmosphérique. L’océan étant toujours en contact avec l’atmosphère, la pression qu’il engendre sur les fonds marins dépend fortement des variations de pression atmosphérique. Il convient ainsi de différencier les modèles océaniques en fonction de la réponse des océans aux variations de l’atmosphère (VAN DAM et al., 2012).

Préparation des données

Les données utilisées ici (figure 3) sont celles employées par A. ASRI (2014) lors de son TFE afin de vérifier nos résultats. Il ne s’agit pas des séries temporelles que j’ai téléchargées précédemment. L’objectif de cette préparation des données est d’enlever la tendance dans le signal et de combler le manque de données (espaces blancs) à certains endroits de la série temporelle (figure 3).
Ensuite, j’ai créé une fonction Python qui permet de détecter la présence de trous dans la série temporelle. Pour cela, j’ai d’abord calculé le temps écoulé entre les deux premières mesures ce qui va correspondre à la fréquence d’échantillonnage de notre série.
Je suis parti du principe qu’il n’y avait pas de trous entre ces deux valeurs. Puis la fonction calcule le temps qui s’écoule entre chaque mesure et l’ajoute dans un tableau si celui-ci est inférieur à deux fois le premier temps calculé. Une fois tous les calculs réalisés, elle effectue la moyenne des valeurs du tableau créé, ce qui détermine le pas final de la série.
Une fois le pas de la série temporelle déterminé, la fonction calcule de nouveau le temps écoulé entre chaque mesure et enregistre dans deux tableaux la date de début du trou et la date de fin du trou si celui-ci est supérieur à trois fois le pas déterminé précédemment, ainsi que les indices de ces valeurs. La fonction renvoie donc le pas, un tableau comprenant les dates de début des trous, un tableau comprenant les dates de fin des trous et un tableau comprenant les indices des valeurs.
Par la suite, j’ai créé une seconde fonction qui va permettre de combler les trous dans nos séries temporelles. Cette fonction demande en paramètres les éléments obtenus précédemment, le signal et les coefficients calculés. Elle va ainsi pouvoir calculer les dates entre le début et la fin d’un trou et calculer une valeur correspondante à chaque nouvelle date en utilisant les coefficients du modèle AR. Les nouvelles valeurs sont obtenues à partir desvaleurs précédentes. Les valeurs qui sont calculées sont ensuite utilisées dans les calculs.

Retranscription du calcul SSA sous python

Objectifs et intérêt de la SSA

L’objectif de la SSA est de décomposer une série temporelle en composantesoscillatoires et en bruits (CLAESSEN et GROTH). L’intérêt de la SSA est de pouvoir ensuite obtenir facilement et rapidement, à l’aide d’un outil d’analyse, la corrélation entre la série temporelle et, par exemple, les signaux ou modèles de surcharge. Cette analyse peut se faire entre chaque modèle de surcharge (atmosphérique, océanique ou hydrologique) ou à partir de la somme des trois.

Le calcul

La retranscription du calcul de la SSA de Matlab à Python n’a pas été toujours évidente. La plupart des fonctions utilisées par A. ASRI existent sous Python. Le plus gros problème que j’ai pu rencontrer a été de gérer les « types » de variable sous Python, entre les listes, les tableaux ou encore les matrices, sachant que Matlab travaille uniquement sur un type matriciel. L’intérêt de cette retranscription est d’avoir le programme de SSA dans le même langage que celui de la MSSA. En ce qui concerne le calcul j’ai suivi les mêmes étapes qu’A. ASRI

Etape 1

L’utilisateur doit dans un premier temps définir le paramètre M. Il s’agit du retard de covariance, c’est l’étape la plus délicate puisque cela revient en fait à tester une certaine périodicité et donc à bien connaître les effets liés à la série étudiée (ASRI, 2014). Afin de choisir au mieux ce paramètre, j’ai effectué des tests qui seront expliqués par la suite (Voir § III.4). Ainsi, le paramètre M définit la « fenêtre » d’étude. Celui-ci peut être rentré directement par l’utilisateur ou calculé par le programme. En effet, il peut être calculé en réalisant la transformée de Fourier sur le signal. Ainsi, à partir des pics fréquentiels, le programme va déterminer la fréquence prédominante et la période correspondante ce qui donne la valeur de M. Une fois M choisi, on découpe le signal en morceaux de longueur M de sorte que :

Application et résultats

Comme A.ASRI (2014) a pu le faire, j’ai appliqué la SSA sur les données GPS de la station BRAZ (figure 13) afin de valider la retranscription sous python du programme de calcul pour pouvoir confronter mes résultats à ceux qu’ils a obtenus.
Après un premier calcul, on peut observer sur la figure 13, le signal d’origine en rouge et les 5 premières composantes que la méthode SSA fait ressortir. Je n’ai représenté que les 5 premières composantes par commodité, sachant que plus on affiche les dernières composantes plus celles-ci ont une amplitude faible et représentent du bruit. Le nombre de composantes calculées est égal à M. Ce que l’on constate en premier lieu et qu’A. ASRI (2014) a observé aussi, est que les composantes 3 et 4 sont identiques ainsi que leurs spectres fréquentiels (figure 14). On aurait pu le deviner à l’aide du diagramme des logarithmes des valeurs propres des composantes. Les valeurs propres des composantes 3 et 4 ayant la mêmeamplitude, on peut en déduire que ces dernières sont identiques et doivent donc êtresommées.

Automatisation du calcul sur le serveur de l’ESGT

La retranscription du calcul de SSA de Matlab à Python avait aussi un autre objectif qui était de rendre son utilisation plus commode ainsi que son automatisation sur le serveur de calcul de l’ESGT pour de nombreuses stations.
Après un entretien avec F. DURAND, Ingénieur d’études et spécialiste en calcul scientifique dans le laboratoire du GeF, nous avons opté pour la création d’un programme en langage linux (.ksh) pour lancer les calculs. C’est ce programme que l’utilisateur va devoir ouvrir et modifier (figure 17).
Les différents paramètres à rentrer par l’utilisateur sont :
– Le ou les noms des stations à calculer (ex : P155, P521)
– L’intervalle (en heure) entre deux mesures GPS (ex : 24)
– La valeur de M (Il suffit de mettre 0 si l’on souhaite que le programme le calcule) (ex : 60)
– Le taux de reconstruction du signal : Ce paramètre va permettre de limiter le nombre de composantes calculées afin que celles ne représentant que du bruit nesoient pas stockées et ne ralentissent pas le calcul (ex : 0.8)
– Les composantes géographiques dans lesquelles on souhaite réaliser le calcul (ex : North, Up)

La MSSA

Le principe

Le principe de la MSSA reste à peu de chose près le même que celui de la SSA.
L’objectif est toujours de décomposer une série temporelle en composantes oscillatoires et en bruits, sauf que cette fois-ci, les données en entrée ne se composent pas d’un seul signal mais de plusieurs. On peut traduire la MSSA comme une combinaison de plusieurs SSA (CLAESSEN et GROTH). Dans notre cas, nous allons utiliser les 3 signaux (East, North et Up) des séries temporelles de la station considérée.
Avant que tous les calculs de la MSSA soient réalisés, les mêmes traitements préalables sur les données GPS ont été faits. C’est-à-dire :
– Suppression de la tendance
– Détection des trous dans les séries GPS et remplissage de ces derniers à l’aide d’un modèle autorégressif.
– Sélection des données GPS et des modèles sur un même laps de temps.

Le calcul

Le calcul de la MSSA a été effectué sur des séries temporelles GNSS disponibles sur le site internet du NGL. En ce qui concerne les différents modèles de surcharge, j’ai utilisé ceux fournis par J-P BOY se trouvant aussi sur le cluster du GeF.
En ce qui concerne le calcul, j’ai suivi les étapes décrites par D. CLAESSEN et A. GROTH.
Ce qui différencie la MSSA de la SSA est que la matrice D n’est plus composée d’un signal mais de tous les signaux que l’on souhaite. Ainsi le nombre de lignes de la matrice D sera équivalent au nombre de valeurs composant la série temporelle et le nombre de colonnes correspond au retard M multiplié par le nombre de signaux étudiés, ici 3.

Reconstruction du signal à partir des composantes

Une fois les composantes obtenues, nous pouvons reconstruire le signal d’origine.
Cette étape nous est nécessaire car nous allons par la suite comparer les composantes reconstruites entres elles, c’est-à-dire que l’on peut recalculer les 3 signaux, East, North, Up. Pour cela, nous avons besoin des composantes non sommées, obtenues par la MSSA et des vecteurs propres calculés précédemment.
Dans un premier temps, nous devons construire la matrice Z. Celle-ci contient les composantes principales mais on rajoute un retard.

Choix de la valeur de M

Comme en SSA, le choix de la valeur de M est primordial. Afin d’observer l’influence de ce dernier sur les résultats, j’ai créé trois signaux synthétiques (figure 21), East, North et Up correspondant à des sinusoïdes simples légèrement bruitées ayant une période de un an mais sans unité en amplitude. D’après les premières observations que l’on peut faire sur les séries temporelles GPS d’origine, on se rend bien compte qu’une périodicité d’un an est présente, c’est pourquoi j’ai choisi cette période. Le principe va être d’effectuer le calcul de la MSSA sur ces signaux avec différentes valeurs de M. Ainsi nous allonspourvoir constater l’impact du M sur nos résultats.
J’ai choisi d’effectuer les calculs avec les valeurs de M suivantes : 31, 183, 340, 380, 740, 1 000 et 1 200. Ces valeurs n’ont pas été choisies de manière aléatoires puissent qu’elles correspondent respectivement à environ 1 mois, 6 mois, un peu moins d’un an, un peu plus d’un an, deux ans et demi et un peu plus de deux ans et demi.
Pour un M = 31 (figure 22), soit un mois, la composante principale numéro une (période de 1 an, amplitude de 10 mm) détectée ressemble fortement aux signaux, et la composante 2 a aussi une période de 1 an mais a une amplitude très faible (0,04 mm). La composante 2 constitue donc le bruit rajouté au signal. Suite au calcul des composantes, nous allons chercher à reconstruire les signaux d’origine. Pour ce faire, nous calculons les composantes reprojetées. Nous allons nous intéresser ici uniquement à la direction Up. Dans ce premier cas, les composantes reprojetées RC1_up et RC2_up sont différentes, et RC1_up représente très bien le signal Up d’origine mais possède encore le « bruit » du signal d’origine.
A l’issue du calcul de la MSSA, nous allons comparer les composantes reprojetées entre les données GPS et les modèles. Or, que la valeur de M soit égale à 183 ou à 1200 (figure 24), les composantes reprojetées, RC1 seront presque identiques (légère différence d’amplitude de 0.03 mm qui n’est pas significative).
Après analyses de ces résultats, on peut en déduire qu’il y a une limite dans la quantité d’information que l’on peut tirer d’un signal par le processus MSSA. Suite à ces tests, j’en ai conclu que le fait d’augmenter la valeur M n’a pas un grand intérêt, il doit dans notre cas avoir une valeur au minimum de 183 pour l’étude des effets de surcharge d’une période de un an. En outre, plus la valeur de M est grande, plus le temps de calcul est long (moins de 5 minutes pour M = 183 à plus de 40 minutes pour M = 1200).
Nous allons maintenant appliquer la MSSA sur des séries temporelles GPS réelles.

Applications et résultats

Cette fois ci, les calculs ont été réalisés sur les séries temporelles GPS de la station POVE obtenues sur le site du NGL. De la même manière que pour la SSA, nous allons pouvoir observer les composantes que la MSSA fait ressortir et calculer les spectres de ces dernières. Sachant que cette fois-ci, la condition qui réalise la sommes des composantes identiques a directement été appliquée.

L’outil d’analyse

Avant d’utiliser l’outil d’analyse développé, nous allons le présenter et expliquer comment l’interpréter et ainsi voir son intérêt dans notre cas. Afin de disposer d’un résultat graphique permettant une interprétation visuelle assez aisée des comparaisons entre les différents signaux (observations et modèles), nous avons opté pour un diagramme de Taylor (TAYLOR, 2001). En effet ce dernier offre l’avantage de regrouper un grand nombre d’informations sur un seul graphique. Ce type de diagramme est notamment utilisé pour comparer les modèles de changement climatique. Il permet de représenter à la fois les corrélations et les écarts-types par rapport à une donnée de référence. De plus, ce type de diagramme n’a pas encore été utilisé en géodésie, ce qui apporte une analyse nouvelle dans ce domaine.

Présentation du diagramme de Taylor

La représentation de l’outil d’analyse se fera donc sous forme d’un diagramme de Taylor (figure 27).

Interprétation du diagramme

Dans cette partie, nous allons voir quelques exemples simples de réalisation de diagramme de Taylor avec des données et des modèles choisis afin de comprendre son interprétation.
Comme premier exemple, nous avons créé des données synthétiques représentant un sinus compris entre 0 et 4π (signal data) ainsi que deux modèles correspondant au même sinus, à savoir des signaux en phase mais avec une amplitude différente. Le premier modèle représente la fonction 1.3 × sin(?)et le second modèle la fonction0.8 × sin(?). Avec ces deux modèles, nous allons pouvoir voir l’influence de la différence d’amplitude entre les données et les modèles sur le diagramme.

Analyses

Choix des stations

Suite au développement du programme de MSSA et de l’outil d’analyse, nous avons sélectionné un sous réseau de stations en Amérique du Sud pour appliquer la méthode etanalyser les signaux de surcharge inclus dans les séries temporelles.
J’ai retenu 13 stations pour mes analyses sachant que le calcul est possible pour un plus grand nombre de stations, pourvu que l’on possède les séries temporelles GPS ainsi que les modèles sur une période de temps suffisamment longue (minimum 5 ans) pour que les calculs et les analyses puissent se faire. Deux d’entre elles, PRGU et UNPA, ont une caractéristique qui rend leur analyse différente. En effet, PRGU, présente un saut dans les données. Nous allons donc nous intéresser à l’influence de ce dernier sur le calcul et les résultats. La station UNPA quant à elle est située extrêmement loin des trois rivières principales d’Amérique du sud (Amazone, Orénoque, Paraná). Ainsi, on observera si l’ajout du modèle de surcharge hydrologique contenant l’impact de ces trois rivières aura uneinfluence sur les résultats.
Afin d’expliquer au mieux les résultats obtenus, nous avons superposé les stations sur différentes cartes. La première carte possède un fond topographique (figure 31) (WESSEL et al., 2013) et la seconde un fond des différents climats (carte de Köppen (HUFTY, 2001)). L’objectif est de pouvoir interpréter les résultats. La seconde carte se trouve en annexe 1.

Conclusion

Suite au travail réalisé ainsi qu’aux différentes analyses que j’ai pu faire, nous pouvons observer la difficulté à comprendre les phénomènes qui engendrent les mouvements de la croûte terrestre. Cependant, grâce aux travaux de certains scientifiques, comme J-P BOY, qui cherchent à modéliser au mieux les effets de surcharge selon leur type et le lieu étudié, on commence à comprendre l’importance de chaque effet et la nécessité de leur prise en compte dans nos calculs de positions.
Notre étude sur l’Amérique du sud a permis de mettre en avant l’influence considérable de l’hydrographie continentale et plus particulièrement des grand fleuves, l’Amazone, l’Orénoque et le Paraná. Dans cette région où le climat est extrêmement humide avec une saison des pluies, les phénomènes de surcharges hydrographiques sont parmi les plus forts au monde, c’est pourquoi leur prise en compte est extrêmement importante.
Le calcul de la MSSA et l’outil d’analyse nous ont montré leur efficacité et leur intérêt afin de comprendre les phénomènes géodynamiques de la croûte terrestre.
Néanmoins, nous avons pu observer les limites de nos outils. Si le manque de données est trop important, si la longueur de la série temporelle est trop courte ou encore si des sauts sont présents, les résultats peuvent êtres faussés comme on a pu le voir sur la station UNPAou les calculs peuvent tout simplement être impossibles.
Afin de poursuivre cette étude, il serait pertinent de chercher maintenant à prendre en compte les différentes difficultés que l’on a pu rencontrer, notamment les sauts. De plus, il faudrait approfondir les recherches et les analyses sur beaucoup plus de stations ou sur des signaux autres que de l’annuel (semi-annuel, interannuel), afin de confirmer nos conclusions et comparer les données avec de nouveaux modèles. Il serait pertinent de réaliser des analyses spatiales suites aux résultats obtenus. La prise en compte des composantes 3D de la MSSA et non des composantes reprojetées comme on a pu le faire serait une autre analyse à réaliser afin de comprendre encore mieux les phénomènes géodynamiques de la Terre.

Table des matières
Remerciements 
Liste des abréviations 
Table des matières 
Introduction
I LES DONNEES
I.1 PRESENTATION
I.1.1 Zone d’étude
I.1.2 Les séries GPS
I.1.3 Les modèles de surcharges
I.2 PREPARATION DES DONNEES
II RETRANSCRIPTION DU CALCUL SSASOUS PYTHON
II.1 OBJECTIFS ET INTERET DE LA SSA
II.2 LE CALCUL
II.3 APPLICATION ET RESULTATS
II.4 AUTOMATISATION DU CALCUL SUR LE SERVEUR DE L’ESGT
III LA MSSA 
III.1 LE PRINCIPE
III.2 LE CALCUL
III.3 RECONSTRUCTION DU SIGNAL A PARTIR DES COMPOSANTES
III.4 CHOIX DE LA VALEUR DE M
III.5 APPLICATIONS ET RESULTATS
IV L’OUTIL D’ANALYSE
IV.1 PRESENTATION DU DIAGRAMME DE TAYLOR
IV.2 INTERPRETATION DU DIAGRAMME
V ANALYSES
V.1 CHOIX DES STATIONS
V.2 LES COMPARAISONS REALISEES
V.3 ANALYSE DES RESULTATS
Conclusion 
Bibliographie 
Table des annexes
Annexe 1 Carte des climats (classification de Köppen) d’Amérique du sud avec les stations étudiées
Annexe 2 Résultats sur la station SAGA
Annexe 3 Résultats sur la station NAUS
Annexe 4 Résultats sur la station BELE
Annexe 5 Résultats sur la station KOUG
Annexe 6 Résultats sur la station BOAV
Annexe 7 Résultats sur la station MAPA
Annexe 8 Résultats sur la station SALU
Annexe 9 Résultats sur la station BRFT
Annexe 10 Résultats sur la station IMPZ
Annexe 11 Résultats sur la station BRAZ
Liste des figures ..

Lire le rapport complet