Comparaison des données simulées avec des données réelles (IFSTTAR)

Télécharger le fichier pdf d’un mémoire de fin d’études

Description des méthodes envisageables

Les fondations dont on cherche la géométrie sont situées dans la très proche surface. Ceci restreint le champ des méthodes d’imagerie utilisables à l’imagerie électrique, au radar de proche surface et à certaines méthodes sismiques que nous allons décrire.

Méthodes électriques

Sondage électrique Le sondage électrique a pour but d’estimer la résistivité d’un sol. La démarche de mesure consiste à faire circuler un courant électrique dans le milieu étudié à l’aide de deux électrodes puis à mesurer la différence de potentiel induite avec deux électrodes de mesure, les électrodes étant classiquement placées en surface. Il existe plusieurs stratégies de positionnement des électrodes selon la zone que l’on cherche à son-der. Une démarche classique est la recherche d’un profil vertical de résistivité. On centre chacune des paires d’électrodes sur la verticale du milieu que l’on cherche à caractériser. Dans cette configuration, on réalise différentes mesures de potentiel en faisant varier l’espacement des électrodes de mesure, ce qui correspond à augmenter la profondeur de la région caractérisée. En faisant l’hypothèse que le milieu est de résistivité homogène et connaissant les positions des électrodes, on peut déduire la valeur de cette résistivité pour chaque espacement entre les électrodes de mesure. Cependant, l’hypothèse d’ho-mogénéité étant généralement fausse, on procède alors à une inversion pour retrouver la véritable résistivité du milieu à partir des résistivités apparentes mesurées. On obtient alors un profil vertical de résistivité. En extrapolant, on conçoit qu’en translatant l’en-semble du dispositif d’acquisition à la surface du milieu, on obtient des mesures de la résistivité apparente du milieu à la verticale d’autres positions. La démarche de mesure précédente permet donc théoriquement de réaliser des sections 2D voire des descriptions 3D de la résistivité du milieu selon les déplacements du dispositif d’acquisition à la sur-face du milieu et en faisant des mesures pour plusieurs espacements entre les électrodes pour caractériser plus profondément le milieu. Cependant, il serait trop long et fastidieux de le faire manuellement c’est pourquoi on a recours au panneau électrique.
Panneau électrique Le panneau électrique consiste à répartir uniformément un grand nombre d’électrodes sur un segment à la surface d’une section 2D ou d’une portion 3D du milieu dont on souhaite caractériser la résistivité. Pour les profils 2D, on utilise en général plusieurs dizaines d’électrodes – typiquement une cinquantaine – et au lieu de déplacer les électrodes, un appareil désigne par un jeu de commutations les deux électrodes d’injection et les deux électrodes de mesure parmi l’ensemble des électrodes installées. La stratégie de choix des électrodes permet de régler la sensibilité de la mesure  à l’orientation des structures enfouies, à la profondeur et influe sur le rapport signal sur bruit (Samouëlian, 2005). De plus, le système de mesure automatique exploite la rapidité de mesure pour reproduire plusieurs fois chaque mesure et ainsi améliorer la précision du résultat. Notons que les mesures que l’on réalise sont identiques à celles que l’on aurait pu faire en déplaçant manuellement les électrodes mais avec un gain de temps et de qualité de mesure considérable.
Les méthodes électriques sont classiquement utilisées pour caractériser des terrains de l’échelle à métrique à décamétrique, le réglage de la résolution se faisant par le choix de l’espacement des électrodes de la dizaine de mètres à la dizaine de centimètres. Les applications de ce type de méthode d’imagerie électrique sont nombreuses : elle est uti-lisée en archéologie pour identifier les zones intéressantes d’un site (Osella et al., 2005), en hydrologie pour quantifier la salinité de l’eau (Amidu et Dunbar, 2008) ainsi qu’en géotechnique pour réaliser l’investigation de sites (Denis et al., 2002; Sudha et al., 2009) ou encore pour imager des fissures dans des structures en béton (Lataste et al., 2003).
Dans le cas des fondations de pylônes électriques en béton, deux difficultés se posent à l’utilisation de méthode électrique. D’une part, la méthode doit être utilisable quelque soit la nature du sol, en particulier si le sol est conducteur comme c’est le cas des sols argileux ( conductivité de l’ordre de la 10−1 à 10− 2S/m ). Or la conductivité électrique du béton étant plus faible ( de l’ordre de 10−4 S/m ), les lignes de courant contournent la fondation de telle sorte que d’un point de vue électrique, la fondation est court-circuitée par le sol plus conducteur qui l’entoure ; dans un environnement plus conducteur que le béton, l’information du panneau électrique 2D est donc erronée par ces importantes fuites de courant en dehors de la section imagée. D’autre part les fondations de pylônes électriques contiennent des câbles de cuivre assurant la mise à la terre de la fondation et qui, par leur forte conductivité (∼ 6.106S/m), vont canaliser les lignes de courant et donc diminuer la sensibilité de l’imagerie électrique à la présence du béton.
En résumé, la diversité des milieux dans lesquels se trouvent les fondations de pylônes électriques est importante et il peut se trouver des milieux très conducteurs ou contenant des éléments métalliques susceptibles de perturber la localisation des contrastes de ré-sistivité et donc la caractérisation de la géométrie de la fondation. Il est donc nécessaire d’envisager d’autres méthodes.

Méthode électromagnétique

GPR Une alternative intéressante est l’utilisation des ondes électromagnétiques. Le radar à pénétration de sol, ground penetrating radar (GPR) en anglais, est un outil clas-siquement utilisé pour la caractérisation de la subsurface (Mari et al., 1998). Son principe consiste à générer une onde électromagnétique avec une antenne émettrice et à enregis-trer le champ d’onde résultant de l’interaction de l’onde incidente avec le sol en utilisant une antenne réceptrice. L’onde incidente va être perturbée par l’impédance électroma- gnétique du milieu Z = 0 r/ 0 r , où 0 est la perméabilité magnétique du vide et vaut 4π10−7H/m, 0 est la permittivité diélectrique du vide et vaut 1/(36π)10−9F/m, r est la perméabilité magnétique relative du milieu et r est la permittivité diélectrique relative du milieu. Cette impédance va notamment créer des réflexions aux endroits où il y a de forts contrastes. Les sols dans lesquels nous allons travailler sont dépourvus de métaux et de minéraux si bien que l’on peut simplifier l’expression de leur impédance :
Z = 0/ 0 r. La démarche d’acquisition la plus courante consiste à maintenir un dé-port constant entre les antennes. En effet, avec cette configuration d’acquisition, la forte directivité des antennes aussi bien à l’émission qu’à la réception favorise l’enregistre-ment des ondes réfléchies ; d’autre part, la simplicité de mise en oeuvre de cette méthode et la rapidité de l’interprétation sont des atouts importants. La fréquence centrale des antennes utilisées est choisie en haute fréquence, au delà de la fréquence de coupure séparant la zone basse fréquence – fortement atténuante et dispersive – de la zone haute fréquence – modérément atténuante et propagative. La gamme de fréquence centrale des antennes d’émission utilisées en GPR s’étend de la dizaine de M Hz au GHz. La réso-lution – au sens du quart de la longueur d’onde – que l’on peut attendre de l’imagerie radar dépend de la vitesse des ondes dans le milieu considéré. En régime propagatif et dans un sol homogène, la vitesse dépend principalement de la permittivité diélectrique du sol r selon la loi v = c/√ r où c est la vitesse de la lumière dans le vide. Le tableau 1.1 présente la résolution verticale que l’on peut attendre de l’interprétation d’une ac-quisition radar effectuée avec des antennes de fréquence centrale d’émission de 500M Hz en négligeant l’influence de la puissance d’émission, c’est-à-dire en acceptant une erreur de l’ordre de 5% entre le fréquence émise et la fréquence dominante enregistrée.

Imagerie des temps de première arrivée

Historiquement, l’imagerie des temps de première arrivée a été introduite pour amé-liorer la localisation des séismes (Aki et Lee, 1976). Cette méthode consiste à perturber un modèle acoustique du milieu étudié de telle sorte que les temps mis par le front d’onde le plus rapide pour arriver sur les récepteurs dans le milieu synthétique soit le même que dans le milieu réel. On suppose alors que le milieu synthétique a les mêmes propriétés que le milieu réel. Cette démarche permet donc d’obtenir une image des propriétés du milieu grâce aux ondes émises par les séismes, mais on peut aussi appliquer la même démarche en utilisant une source active comme un pot vibrant ou un marteau. L’utili-sation de cette méthode est très répandue en sismologie que ce soit à l’échelle globale (Montelli et al., 2004b,a) ou locale (Thurber et al., 1997). En revanche, cette méthode n’est pas utilisée en géotechnique car, en l’absence de réflecteur, le gradient de vitesse au voisinage de la subsurface est trop petit pour redresser suffisamment les rais et faire en sorte que l’on puisse les enregistrer sur le profil d’acquisition.
L’inversion des temps de première arrivée n’utilise qu’une information partielle des sismogrammes, une première idée consiste alors à utiliser la forme d’onde de la première arrivée (Sheng et al., 2006). Dans le prolongement de cette démarche, on peut envisa-ger de prendre en compte l’ensemble des arrivées en temps et en amplitude résultant de l’interaction du front d’onde incident avec le milieu, autrement dit l’intégralité des sismogrammes. Cette démarche s’appelle l’imagerie de la forme d’onde complète – full waveform inversion, (FWI) en anglais (Tarantola, 1984a).

Imagerie sismique de la forme d’onde complète

Il semble assez avantageux d’utiliser l’ensemble des phénomènes enregistrés sur les sismogrammes plutôt que la première arrivée uniquement. Plus on considère un grand nombre de phénomènes dans lesquels s’expriment les propriétés du milieu et meilleure devrait être la caractérisation du sol. Cependant, différentes arrivées sont imbriquées et peuvent rendre les données très complexes ; il est donc important de définir une stratégie pour utiliser efficacement tout le signal. Précisons donc le principe de l’imagerie de la forme d’onde complète.
Principe de la FWI La FWI requiert d’utiliser une source active et de placer des géophones dans le milieu étudié. Lorsque la source émet des ondes dans ce milieu, le signal reçu par les géophones est caractéristique de la source (signal temporel, bande de fréquence et position), du milieu (propriétés mécaniques) et des géophones (positions et bande passante). Pour une acquisition sismique réelle, les positions de la source et des géophones sont connues ; si on dispose d’une estimation initiale des propriétés du milieu et de l’amplitude du signal source en fonction du temps, on peut alors faire une simulation numérique des sismogrammes et comparer les données synthétiques ainsi obtenues avec les données réelles.
Dans l’hypothèse où la modélisation numérique est suffisamment réaliste, c’est-à-dire capable de simuler les phénomènes de la physique des ondes se produisant dans le milieu étudié, les différences observées entre les données réelles et synthétiques sont unique-ment dues à une mauvaise estimation des propriétés du milieu ou du signal source. On peut alors ajuster ces propriétés et le signal source pour améliorer la ressemblance entre les sismogrammes réels et les sismogrammes synthétiques sur toute la durée du sismo-gramme. Cet ajustement peut se faire par tâtonnement, mais une démarche heuristique peut prendre du temps et présente une grande incertitude de succès. Pour automatiser cet ajustement on introduit un critère de distance entre les sismogrammes réels et syn-thétiques et l’on va chercher à le minimiser par une méthode de minimisation de fonction. Pour ce faire, une première idée est d’estimer la sensibilité du critère à une perturbation de chacun des paramètres physiques pour trouver la combinaison de paramètres qui fait décroître le plus le critère. Appliquée telle quelle, cette démarche est très coûteuse et des techniques plus astucieuses conduisent au même résultat pour un coût de calcul bien plus faible (Lailly, 1983; Tarantola, 1984a; Pratt et al., 1998). Or cette direction n’est pas nécessairement celle de la solution puisque la sensibilité du critère à une perturbation du modèle a été estimée au voisinage d’un modèle courant. On réitère donc la démarche jusqu’à ce qu’une petite perturbation du modèle ne puisse plus entraîner une diminu-tion du critère ; on dit alors que le critère est dans un minimum local que l’on considère comme le résultat de l’imagerie. Une fois l’inversion réalisée, de deux choses l’une :
• soit le critère est suffisamment faible et l’on considère que l’inversion a convergé vers le minimum global du problème, ou à proximité ; les sismogrammes synthétiques ressemblant alors au sismogrammes réels, on considère que le milieu reconstruit est proche du milieu réel.
• soit le critère est relativement élevé ce qui est symptomatique d’une mauvaise estimation des propriétés du milieu reconstruit.
Notons que ce que l’on vient de décrire pour une seule position de la source – et donc un seul jeu de sismogrammes – est envisageable pour plusieurs jeux de sismogrammes obtenus en déplaçant la source. L’intérêt de changer la position de la source peut s’illus-trer par l’analogie d’une pièce fermée dans laquelle se trouve un objet éclairé par une source ponctuelle de lumière : pour une position de la source, l’ombre projetée sur les murs apporte une information partielle sur la forme de l’objet, de même qu’un jeu de sismogrammes obtenus pour une position de la source sismique est porteur de la signa-ture du milieu ; en changeant la position de la source, on change l’éclairement de l’objet et, grâce à l’ombre supplémentaire obtenue dans cette deuxième position de la source, on a une meilleure idée de sa forme. Ainsi, en remarquant qu’un même jeu de données peut correspondre à des milieux différents, on conçoit qu’il est intéressant d’augmenter le nombre de récepteurs et le nombre de positions de tir pour restreindre le nombre de milieux reconstruits ambigus avec le milieu solution ; mathématiquement, cela signifie que l’on augmente le nombre de données permettant d’identifier les paramètres inconnus du milieu.
D’après le principe que nous venons de présenter, deux points essentiels doivent être dé-crits : la démarche de modélisation et la stratégie d’inversion. Ces thèmes feront l’objet des chapitres 2 et 3 respectivement.
Il est intéressant ici de connaître l’historique de l’inversion de la forme d’onde sismique : ceci permet de mettre en évidence que cette méthode dont la théorie a bientôt 30 ans est relativement jeune en termes d’application sur des données réelles.
Historique L’imagerie par inversion de données est un domaine des mathématiques appliquées qui requiert des ressources informatiques importantes. C’est pourquoi les ap-plications ont émergé conjointement au développement de l’informatique et de la mo-délisation de phénomènes physiques de plus en plus complexes. L’inversion de données a été appliquée avec succès en imagerie médicale dès le début des années 70 (Herman et al., 1973). Contrairement à la tomographie à rayons X, la propagation d’onde sis-mique ne peut pas être décrite par un modèle linéaire et il a fallu attendre le début des années 80 pour que des sismologues envisagent la mise en oeuvre d’une méthode d’ima-gerie utilisant l’intégralité de la forme d’onde des sismogrammes (Tarantola et Valette, 1982; Tarantola, 1984a). Les premières applications d’imagerie en milieu 2D – supposé invariant selon une direction horizontale – ont été initiées au milieu des années 80 sur des données synthétiques (Gauthier et al., 1986) et réelles (Mora, 1987) au prix d’un fort investissement en temps de calcul. A cette époque, l’imagerie par inversion de la forme d’onde était une méthode demandant des moyens informatiques sensiblement plus importants que les autres et pour un gain modéré ; cependant Pratt (1990) fit des tests synthétiques dans le cadre d’acquisition en transmission entre puits pour comparer les performances de cette méthode avec la tomographie de diffraction et conclut sur la su-périorité de la tomographie de diffraction comparée à une seule itération de la méthode d’inversion de la forme d’onde. Une particularité de sa démarche est le fait de travailler dans le domaine fréquentiel pour, principalement, limiter la quantité de données trai-tées puisqu’alors il inverse seulement une fréquence au lieu d’inverser l’amplitude des sismogrammes pour l’ensemble des instants. Dans la continuité de cette démarche, Pratt et al. (1996) appliquent une inversion séquentielle de certaines fréquences des données de façon croissante à des données sismologiques synthétiques acquises en surface ; selon son modèle de départ de l’inversion, il montre le potentiel intérêt de l’inversion de la forme d’onde à l’échelle crustale appliquée à des données où l’éclairement du milieu est grand pour une acquisition en surface. Depuis, de nombreuses applications sur données réelles ont été menées avec succès en particulier dans la caractérisation de structures géologiques (Dessa et al., 2004; Operto et al., 2006) et plus récemment en exploration pétrolière (Sirgue et al., 2009; Prieux et al., 2011).
En imagerie de subsurface, les applications géotechniques de l’inversion de la totalité de la forme d’onde sont très peu répandues. On pourra cependant citer les travaux de Bre-taudeau et al. (2010) concernant la détection de cavités et de Smithyman et al. (2009) sur la détection d’objets enfouis. Si les succès de la FWI en géophysique sont avérés, il est important de noter que son application à notre cadre d’étude va rencontrer un certain nombre d’obstacles.

Difficultés prévues

Ce travail de thèse est réalisé dans la continuité de la thèse de Gelis (2005) sur la FWI et de Magnin (2008) concernant l’identification de la source sismique appropriée à la problématique de l’imagerie des fondations de pylônes électriques. Concernant l’outil d’imagerie, plusieurs obstacles vont devoir être surmontés :
• La structure que l’on cherche à imager présente un contraste fort avec l’encais-sant, le rapport entre les propriétés de la fondation et celles de l’encaissant étant de l’ordre de dix dans certains milieux. Or la méthode d’imagerie est basée sur l’hy-pothèse de petites perturbations du milieu à chaque itération : plus le contraste est grand et plus le chemin de convergence est long, ce qui fragilise la convergence vers le bon milieu, l’inversion étant alors plus susceptible de se trouver piégée dans un minimum local correspondant à un modèle éloigné de la solution.
• une difficulté de la sismique de proche subsurface est la complexité du champ d’onde. En effet, celui-ci consiste en la superposition des ondes de surface et des ondes de volume. Or ces ondes de volume sont de bien plus faible énergie que les ondes de surface alors que seules les ondes de volume sont porteuses de la signature de la fondation. Il faudra donc veiller à appliquer l’inversion en privilégiant les ondes de volume.
• l’atténuation devrait jouer un rôle important dans l’imagerie. En effet ce para-mètre influe fortement sur les données : si on choisit de ne pas l’inverser, il faudra avoir une estimation digne de confiance de sa distribution spatiale dans le milieu. Si on décide de l’inverser, il faudra être conscient que sa reconstruction dépend beaucoup de la qualité de la reconstruction du milieu en Vp, la vitesse des ondes de compression.
• Certes nous cherchons un objet fortement contrasté par rapport au sol et dont la localisation approximative est connue. Cependant le sol lui même peut être inhomogène et le sera vraisemblablement, en particulier à cause de la zone altérée de plusieurs dizaines de centimètres que l’on trouve classiquement à la surface de tout milieu et donc autour de la fondation. La méthode d’imagerie va donc devoir estimer les propriétés physiques de l’ensemble du milieu.
• L’objet à imager est une fondation de pylône, un objet contenu dans un parallélépi-pède de 4 mètres de hauteur et de section horizontale 2m*2m. Or nous envisageons d’imager cette structure 3D avec un outil d’imagerie 2D pour des raisons de temps de calcul. Il va donc falloir s’assurer que cette hypothèse simplificatrice, consistant à négliger les variations des propriétés du milieu selon la direction per-pendiculaire à la section, ne perturbe pas trop l’imagerie ; en particulier, nous supposons que les diffractions des sommets de la fondation seront faiblement éner-gétiques. D’autre part, nous devrons déployer le dispositif d’acquisition en consé-quence, c’est à dire travailler avec une seule ligne de sources et récepteurs situés dans le plan de symétrie de la fondation.
• Le faible éclairement de la partie inférieure de la fondation rend difficile la connaissance de sa profondeur, nous espérons que l’information donnée par l’impact écho compensera ce déficit d’illumination.
Le travail proposé dans ce document consiste en une étude de l’applicabilité de l’in-version de la forme d’onde à notre problématique. Pour ce faire il a été nécessaire d’avoir des interactions scientifiques avec des partenaires spécialistes de l’acquisition de données et de l’inversion. Nous présentons ces partenaires dans la partie qui suit.

Compétences et moyens mis en oeuvre

Partenaires

Les collaborateurs qui sont intervenus sur ce projet sont :
• RTE qui a formulé la problématique et s’est intéressé au déroulement global du projet. Nos collaborateurs de RTE étaient Gwennou Le Mignon, Paul Penserini et Xavier Waymel.
• L’antenne de Nantes de l’Institut Français des Sciences et Technologies des Trans-ports, de l’Aménagement et des Réseaux (IFSTTAR) – anciennement Laboratoire Central des Ponts et Chaussées (LCPC) – qui s’est chargé de la réalisation de ma-quettes reproduisant des fondations en modèle réduit dans le but de permettre le test des algorithmes de reconstruction dans des milieux contrôlés. Nos collabora-teurs d’IFSTTAR étaient Philippe Cotte (acquisitions sur site), Olivier Durand (acquisitions sur maquette et sur site), Donatienne Leparoux (acquisitions sur ma-quette) et Anaëlle Luczak (acquisition sur site).
• TerraSeis s’est occupé de mettre en place le protocole d’acquisition de données : son rôle fut notamment de caractériser la diversité des terrains dans lesquels peuvent se trouver les fondations et de mettre en place une acquisition sur site pour obtenir des données réelles destinées à l’inversion de la forme d’onde. Notre collaborateur à TerraSeis était Olivier Magnin.
• La branche Recherche et Développement d’Electricité De France (EDF R&D) s’est intéressée à la centralisation des idées proposées et développées par les partenaires en charge de la modélisation et de l’inversion. Nos collaborateurs d’EDF ont été Alexandre Girard (inversion), Jean-Marie Hénault (impact écho), Nicolas Paul (in-version), Laurent Ulpat (impact écho) et Guy d’Urso (impact écho).
• Apside a participé à la modélisation en forte implication avec EDF R&D. Nos collaborateurs d’Apside étaient Laurie Cham-Lan et Steven Kerzale.
• L’Institut de Recherche en Communication et Cybernétique de Nantes (IRCCyN) qui a travaillé sur l’inversion de données avec introduction d’a priori. Nos collabo-rateurs à L’IRCCyN étaient Jérôme Idier et Denis Vautrin.
• L’Institut des Sciences de la Terre (ISTerre) de Grenoble et le laboratoire Grenoble Image Parole Signal et Automatique (GIPSA-Lab) ont travaillé sur l’acquisition, la modélisation du problème direct, le traitement des données et sur l’inversion des données sismiques. Les chercheurs qui ont encadré mon travail ont été mes directeurs de thèse, Jean Virieux et Jérôme Mars, ainsi que Romain Brossier. Enfin, des doctorants et post-doctorants de l’ISTerre ont participé à la mise en oeuvre de plusieurs acquisitions sismiques sur le site test de Grenoble. Il s’agissait de Bastien Dupuy, Ganghui Hu, François Lavoué, et Ludovic Métivier.
Un tel développement de compétence requiert des interactions fréquentes pour synchro-niser les avancées des différentes équipes : ainsi, une réunion technique entre EDF-R&D, l’IRCCyN et l’ISTerre/GIPSA-Lab a été organisée tous les six mois au cours de ma thèse pour diffuser l’information de l’avancement des volets Expérimental et Numérique du projet. Dans le courant du mois de juin 2009, je suis allé sur le lieu de travail des partenaires EDF-R&D à Chatou pour discuter de l’implémentation du problème direct en différences finies et par la méthode Galerkin discontinu. Début juillet 2009, le IFST-TAR et TerraSeis ont réalisé une campagne de mesure sur des fondations implantées sur le campus de St Martin d’Hères avec l’aide du doctorant de l’IRCCyN, Denis Vautrin, et de deux doctorants d’ISTerre Bastien Dupuy et moi-même. Fin Juillet 2009, je me suis rendu à Nantes pour observer la procédure d’acquisition sur maquette réalisée à IFSTTAR. Les 23 et 24 Mars 2010, Romain Brossier, Philippe Cotte, Steven Kerzalé, Guy D’Urso, Jean Virieux et moi-même avons procédé à une caractérisation du terrain au voisinage de fondations du site du Chesnoy. Ensuite les 5 et 6 Octobre 2010, Romain Brossier, Philippe Cotte, Olivier Magnin et moi-même avons mis en place une acquisition sismique de caractérisation du sol sur deux sites d’Hostun pendant que Guy d’Urso et Jean-Marie Hénault procédaient à des acquisitions d’impact écho. Parallèlement à ces réunions techniques, des réunions semestrielles étaient organisées avec tous les interve-nants : elles étaient l’occasion pour RTE de constater l’avancement global du projet.

Moyens mis en oeuvre

Ma thèse a été financée par une bourse ministérielle associée à un monitorat délivrée par le Ministère de l’Enseignement Supérieur et de la Recherche. Les frais logistiques (déplacements sur les lieux de conférence, compte informatique, logiciels) ont été pris en charge par RTE. Les ressources informatiques nécessaires à ma contribution au projet ont été :
• les serveurs de calcul R2d2 et Gofree de l’OSUG gérés par Françoise Roch et Bruno Breznik
• le serveur de calcul Vargas de l’IDRIS
• le serveur de post traitement 6po de l’OSUG et le serveur de l’équipe Risques du laboratoire ISTerre gérés par Patrick Fulconis

Cadre Physique

La modélisation des ondes sismiques consiste à décrire la propagation des vitesses particulaires et éventuellement des contraintes en chaque point du milieu. D’autre part, la modélisation requiert de faire des hypothèses simplificatrices sur le terrain et il faut justifier la pertinence de ces approximations dans notre cadre de travail.
Dans notre cas, la modélisation que nous utilisons est une modélisation en éléments finis Galerkin Discontinus 2D P-SV dans le domaine fréquentiel. Détaillons les hypothèses associées à ce choix et son intérêt dans le cadre de l’imagerie des fondations de pylônes électriques.

Hypothèse d’un milieu 2D

On considère que le milieu de propagation est à deux dimensions : on suppose ainsi qu’il est invariant suivant l’une des directions horizontales. Cette hypothèse est bien sûr fausse dans le cas d’un milieu qui contient une fondation de pylône électrique, cependant le temps de calcul des simulations 3D est prohibitif : par exemple une modélisation de 160 tirs dans l’un des milieux tests de Grenoble avec une méthode de type éléments finis discontinus en 3D dans le domaine temporel demande 40000h de calcul contre 45h en 2D dans le domaine fréquentiel. L’étude de l’applicabilité de l’imagerie sismique à la caractérisation géométrique des fondations de pylônes requiert donc de travailler en 2D, comme ont pu le faire Smithyman et al. (2009) en géotechnique ou Prieux et al. (2010) en prospection pétrolière sur des données réelles. La figure 2.1 illustre un exemple de géométrie 3D de fondation (cf a et b) et le résultat d’imagerie que nous souhaitons idéalement dans ce cadre, en l’occurrence la section de la fondation dans son plan de symétrie (cf c). Ce point conditionne la géométrie d’acquisition : les sources et les récepteurs doivent alors être situés dans ce plan compte tenu de la modélisation 2D. Cette comparaison devrait mettre en évidence l’influence de l’effet 3D, en particulier celui des diffractions sur les coins de la fondation situés en dehors de la section imagée. Cependant, l’existence d’une différence n’est pas rédhibitoire : ce n’est pas parce que le modèle ne peut pas reproduire certains phénomènes, que l’inversion ne peut pas converger. En effet, il est par exemple possible d’utiliser un modèle acoustique pour inverser des données élastiques, l’idée étant que les phénomènes élastiques ne peuvent pas être interprétés comme des phénomènes acoustiques et que l’inversion se contente alors d’expliquer ce qu’elle est en mesure d’expliquer, à savoir la propagation des ondes de compression. Cependant, la présence de ces phénomènes non pris en compte par la modélisation réduit la résolution de l’imagerie. C’est pourquoi, lorsque c’est possible, on élimine les arrivées correspondant à des phénomènes physiques que l’on ne modélise pas ; ainsi Operto et al. (2006) annule les arrivées correspondant à des phénomènes de cisaillement avant de mettre en oeuvre une inversion acoustique. Dans le cas de données 3D, il est difficile d’identifier clairement des composantes liées à la géométrie 3D du milieu qu’une modélisation 2D ne pourrait pas reproduire. Il existe cependant des traitements de données corrigeant l’atténuation géométrique et l’extension du dispositif d’acquisition – supposé infini dans la direction transversale en 2D – qui permettent de convertir des données 3D en données 2D ; nous abordons ces corrections au chapitre 5. Précisons que ces traitements de données ne peuvent bien sûr pas corriger le caractère 3D d’une structure enfouie ce qui peut s’avérer problématique si cette structure a une faible extension transversale.
Abordons les hypothèses physiques faites sur le milieu de propagation.

Hypothèse d’un milieu visco-élastique

Les lois physiques Dans le cadre des ondes mécaniques élastiques, les grandeurs qui se propagent sont des ondes d’allongements relatifs et de contraintes dans le milieu. Elles se propagent en raison du couplage entre les champs d’allongement relatif et de contrainte qui intervient dans la loi de Hooke ainsi que dans le principe fondamental de la dynamique.
L’équation (2.1) est l’expression de la loi de Hooke en un point M d’un milieu à 2 dimensions. σxx λ + 2 0 xx σxz 0 0 2 xz σzz = λ λ + 2 0 zz (2.1)
où σxx et σzz sont respectivement les contraintes en compression horizontale et verticale au point M et σxz est la contrainte en cisaillement au point M. xx et zz sont respective-ment les allongements relatifs horizontal et vertical au point M et xz est l’allongement relatif de cisaillement au point M. Les propriétés du milieu sont décrites par les coeffi-cients de Lamé λ et . La définition des allongements relatifs est précisée en (2.2).xx = ∂ux (2.2)

Hypothèse d’un milieu atténuant

Le modèle scalaire d’atténuation classiquement utilisé en milieu visco-élastique dans le domaine temporel est celui de Boltzmann tnσ(t) = ψ(t − τ ) ˙(τ )dτ (2.9)
où σ est la contrainte, ˙ est la dérivée temporelle du déplacement relatif et ψ est la fonction de relaxation de la contrainte. L’interprétation physique de cette relation est que la proportionnalité entre déformation relative et contrainte n’est plus satisfaite en raison de la mémoire du milieu. Par propriété de l’opération de convolution, on peut écrire l’expression (2.10) équivalente à l’expression (2.9). −∞ ψ(t − τ ) (τ )dτσ(t) = ˙ (2.10)
On peut alors considérer que ˙ est la réponse en contrainte à une impulsion de défor- ψ mation. Notons que l’on retrouve la loi de Hooke dans le cas limite où ˙ , ψ = Cδ(t) δ(t) étant la fonction Dirac et C un coefficient de déformation élastique. Dans notre cas, nous faisons l’hypothèse que l’atténuation du milieu est indépendante de la fréquence et on note qp et qs les coefficients de qualité des ondes P et S respectivement. En conclusion, le milieu modélisé est décrit par 5 paramètres dépendant de l’espace : Vp, Vs, ρ, qp et qs.
Justification d’une modélisation visco-élastique Les ondes les plus énergétiques qui se propagent dans un volume sont en général les ondes de compression. Il est donc essentiel de les modéliser dans la perspective de l’imagerie d’un milieu. Dans certaines cir-constances, l’utilisation de ces seules ondes permet d’imager efficacement des structures géologiques. C’est le cas de certaines acquisitions marines dans lesquelles des capteurs sont placés sur le fond océanique (Operto et al., 2006; Kamei et al., 2011). Cependant, lorsque le milieu est complexe, il peut être intéressant voire nécessaire de considérer d’autres phénomènes dans la perspective de l’imagerie. Dans le cas des fondations de pylônes électriques, nous choisissons une description visco-élastique du milieu qui prend en compte les ondes de compression, les ondes de cisaillement et leur atténuation. La mo-délisation des ondes de cisaillement est cruciale dans le cas de l’imagerie de très proche surface (<5m) car elle permet la modélisation des ondes de surface qui comme nous le verrons dominent les sismogrammes. Elle permet également d’expliquer les conversions d’onde qui vont avoir lieu à l’interface du milieu et de la fondation. D’autre part, les inhomogénéités du milieu que l’on rencontre au voisinage de la surface entraînent une atténuation importante des ondes qu’il est donc pertinent de prendre en compte dans la modélisation.
Notons qu’il existe d’autres types de modèles prenant en compte des phénomènes supplémentaires comme la modélisation de propagation d’ondes en milieux anisotropes ou en milieux poreux. Pourquoi ne pas les avoir choisies ? Une première raison est que les fondations ont normalement été placées dans des milieux suffisamment consolidés pour assurer un transfert efficace des efforts, il n’est donc pas a priori nécessaire d’envisager la porosité du milieu. Ensuite, ces modélisations décrivent une physique plus complexe : elles prennent donc en compte plus de paramètres et elles sont plus coûteuses en temps de calcul. Enfin, les paramètres physiques les plus pertinents à inverser dans ces types de milieu ne sont pas encore bien définis.
Nous avons présenté le système d’équations d’onde 2D PS-V en temporel ainsi que les ondes classiques qu’il permet de modéliser et nous allons maintenant nous intéresser à sa résolution. Ce système d’équations peut être résolu semi-analytiquement avec certaines conditions de distribution de propriétés du milieu ; ainsi Garvin (1956) et Bouchon et al. (1989) ont proposé respectivement la solution de l’équation d’onde en milieu homogène avec surface libre puis en milieu stratifié. Cependant il n’y a pas de méthode analytique donnant une solution dans le cas général. Pour trouver les champs d’ondes générés dans un milieu arbitrairement hétérogène, il est nécessaire de recourir à des méthodes numé-riques.
Présentons maintenant la résolution numérique de l’équation d’onde élastique.

Cadre Numérique de la résolution de l’équation d’onde

Les lois physiques et le système d’équations qui en découle sont formulées localement : pour trouver le champ d’onde créé par une source dans un milieu, les méthodes numé-riques procèdent à une discrétisation de l’espace physique en petites régions. On conçoit que la façon dont on effectue cette discrétisation spatiale va influer sur la précision de la modélisation. Dans une première partie, nous aborderons la méthode des différences finies – démarche la plus courante en modélisation visco-élastique à l’heure actuelle -puis celle des éléments finis Galerkin-Discontinus, méthode que nous utilisons. Nous pré-senterons leur principe en précisant leurs avantages et inconvénients respectifs. Ensuite, nous savons qu’il est possible de représenter de façon équivalente un signal temporel par sa représentation dans le domaine fréquentiel grâce à la transformée de Fourier. Cette représentation est cruciale pour notre problématique d’imagerie ; nous présenterons son incidence sur la formulation du problème direct.

Les méthodes de modélisation

Dans cette partie, nous présentons le principe des différences finies et des éléments finis Galerkin-Discontinus ainsi que leur pertinence dans la problématique de l’imagerie des fondations de pylônes électriques.
Pour illustrer leur principe, nous travaillerons sur l’équation de propagation acoustique en 1D comme le suggèrent Hesthaven et Warburton (2008) et Etienne (2011). L’équation d’onde s’exprime alors de la façon suivante : ∂tu(x, t) + A ∂xu(x, t) = 0 (2.11)
où le vecteur u désigne les champs physiques de vitesse particulaire et de contrainte u = (v, σ)t (2.12) et A désigne la matrice d’impédance du milieu : 1 A = E 0 0 ρ où E est le module d’Young du matériau et ρ sa masse volumique.
Dans ce cadre présentons les différentes méthodes de discrétisation de cette équation.
Différences finies La méthode des différences finies classiques est basée sur l’approxi-mation de la dérivée d’une quantité par des différences de cette quantité. Ainsi la dérivée temporelle définie par l’équation (2.13) est approchée par la formule (2.14). ∂tu(x, t) = lim u(x, t) − u(x, t − τ ) τ →0 τ u(x, t) − u(x, t − Δt) Δt (2.13) (2.14)
Δt étant le pas d’échantillonnage temporel choisi. Ce terme est constant et il influence de façon critique la précision de la modélisation. L’expression (2.14) est appelée dérivée temporelle décentrée d’ordre 2 : cette approximation est la plus utilisée en modélisation sismique en temps (Virieux, 1986; Saenger et al., 2000). De la même manière il est possible d’approximer la dérivée spatiale d’un champ par la relation suivante : ∂xu(x, t) = lim u(x + ξ, t) − u(x, t) ξ→0 ξ u(x + Δx, t) − u(x − Δx, t) 2Δx (2.15)(2.16) où Δx est le pas d’échantillonnage spatial choisi. Ce terme est constant et il influence de façon critique la précision de la modélisation. L’expression 2.16 est appelée dérivée spatiale centrée d’ordre 2. En utilisant le développement de Taylor, on montre que cette expression a une précision de l’ordre de (Δx)3 et qu’il est possible d’améliorer cette précision en utilisant des développements à des ordres supérieurs. L’expression suivante par exemple permet d’avoir une erreur de l’ordre de (Δx)5 : ∂tu(x, t) = 8 u(x + Δx, t) − u(x − Δx, t) − u(x + 2Δx, t) − u(x − 2Δx, t) + (Δx)5 12Δx
Ainsi, pour améliorer la précision d’une modélisation, on peut diminuer la taille du maillage ou bien augmenter l’ordre d’interpolation des opérations de dérivée ce qui, dans les deux cas, entraîne une augmentation du coût de calcul. Ce point est commun à toutes les méthodes de modélisation.
En différences finies, les deux grandeurs Δx et Δt sont en général constantes : on ne les ajuste pas localement. D’autre part, il existe un troisième point qui influe sur la qualité de la modélisation des différences finies : il s’agit de la façon dont on répartit les propriétés et les champs localement, au niveau de chaque cellule, ce que l’on appelle le stencil. Le stencil le plus simple consiste à attribuer à chaque point de la grille les propriétés élastiques du milieu et les champs de vitesse-contrainte. Cependant, des tests numériques montrent que ce n’est pas le choix le plus efficace en termes de calcul pour une même précision et d’autres stencils ont été proposés dans lesquels chaque champ n’est pas exprimé en chaque point dans le but de nuancer le caractère ponctuel de la formulation des différences finies alors que la physique est continue.
Les différence finies ont eu du succès en modélisation sismique dès le début des années 80 (Virieux et Madariaga, 1982; Virieux, 1986) en raison de leur efficacité : pour une précision donnée, c’est la méthode la plus rapide. Concernant la façon de considérer les propriétés locales du milieu dans les équations, plusieurs solutions existent : le stencil en quinconce (Virieux, 1984) présente de bonnes performances ainsi que le stencil tourné en quinconce – rotated staggered grid (Saenger et al., 2000) qui est moins efficace que le stencil non tourné en termes de dispersion numérique mais qui permet de modéliser plus précisément les milieux avec de fortes hétérogénéités. Les différences finies sont basées sur un développement en séries de Taylor de chaque champ considéré et ce type de développement suppose que les champs ont des propriétés de dérivabilité d’autant plus grandes que le développement est fait à un ordre élevé. Or, si les propriétés du milieu sont localement très contrastées – ce qui peut se produire au voisinage de la surface libre, dans le cas d’une faille ou encore d’une cavité – cette hypothèse n’est pas vérifiée. La modélisation par différences finies est donc une méthode adaptée aux milieux faiblement hétérogènes. Le fait de raffiner le maillage nuance la mauvaise modélisation des hétérogénéités, mais le fait de ne pas pouvoir ajuster localement la dimension de la grille est une contrainte qui alors devient critique. En effet, il peut se trouver de petites régions qui nécessitent une maille bien plus petite que d’autres, ce qui entraine un pas de maillage excessivement fin dans le reste du milieu.
Dans la problématique des fondations, ce point est particulièrement critique pour la modélisation des ondes de Rayleigh. Ces ondes de surface requièrent une densité de maillage allant de 15 à 60 cellules par longueurs d’onde, selon que l’on considère une topographie plate ou complexe (Bohlen et Saenger, 2006) c’est-à-dire 1.5 à 6 fois supé-rieure à celle réclamée par les ondes de volume en milieu homogène. Le problème de la topographie complexe ou des hétérogénéités pentées par rapport au maillage est qu’il est difficile de décrire un segment penté avec une grille cartésienne, le cas le plus défavorable étant une pente de 45°. Lorsque la topographie n’est pas plate, le maillage cartésien des différences finies ne peuvent l’approcher que par une succession de marches d’escalier qui sont source de diffractions parasites si l’on ne choisit pas des marches suffisamment petites devant la longueur d’onde. Notons qu’il existe des propositions alternatives à la description de la surface sous la forme de marche d’escalier (Lombard et al., 2008), mais elles ne sont pas appliquées en raison de leur complexité de mise au point. La contrainte d’un pas de maillage uniforme est donc un inconvénient majeur des différences finies comparées aux autres méthodes dans le cas de milieux réalistes. Une solution consiste-rait à implémenter les différences finies avec un pas de maillage variable spatialement, cependant ce n’est pas trivial mathématiquement et cela fait actuellement l’objet de re-cherches (Bon et al., 2009; Tarrass et al., 2011).
Dans le cadre de notre problématique d’imagerie des fondations de pylônes, on renonce donc à utiliser les différences finies car elles sont trop coûteuses lorsqu’il s’agit de modé-liser des surfaces complexes.

Table des matières

Introduction générale
Présentationde la problématique 
1.1 Objectif
1.2 Le choix de l’imagerie sismique
1.2.1 Description des méthodes envisageables
1.2.1.1 Méthodes électriques
1.2.1.2 Méthode électromagnétique
1.2.1.3 Méthodes sismiques
1.2.1.4 Imagerie des temps de première arrivée
1.2.1.5 Imagerie sismique de la forme d’onde complète
1.2.2 Difficultés prévues
1.3 Compétences et moyens mis en oeuvre
1.3.1 Partenaires
1.3.2 Moyens mis en oeuvre
2 Modèle direct 
2.1 Description physique
2.1.1 Cadre Physique
2.1.1.1 Hypothèse d’un milieu 2D
2.1.1.2 Hypothèse d’un milieu visco-élastique
2.1.1.3 Equation de propagation
2.1.1.4 Hypothèse d’un milieu atténuant
2.1.2 Cadre Numérique de la résolution de l’équation d’onde
2.1.2.1 Les méthodes de modélisation
2.1.2.2 Discrétisation temps-fréquence
2.1.2.3 Les limites physiques du milieu
2.2 Comparaison des données simulées avec des données réelles (IFSTTAR)
2.2.1 L’acquisition sismique sur maquette
2.2.2 La maquette IFSTTAR
2.2.3 Simulations
2.2.3.1 Cadre de simulation
2.2.3.2 Résultats de simulation
2.2.4 Conclusion
3 Inversion 
3.1 Principe de l’inversion
3.2 Fonction coût
3.3 Inversion de la source
3.3.1 Une source méconnue
3.3.2 Démarche d’estimation
3.3.3 Intérêts d’inverser la source
3.4 Inversion du milieu par FWI
3.4.1 Méthode du gradient
3.4.2 Le gradient conjugué
3.4.3 La méthode de Newton
3.4.3.1 Méthode de Newton dans un milieu à une dimension
3.4.3.2 Méthode de Newton en dimension N
3.4.4 Calcul du gradient
3.4.4.1 Expression du gradient
3.4.4.2 Interprétation physique du gradient
3.4.5 Calcul du hessien
3.4.5.1 Expression et Interprétation du hessien
3.4.5.2 Les méthodes de Quasi-Newton
3.4.5.3 L’algorithme L-BFGS
3.4.6 Line search
3.5 Choix de la Norme
3.6 Régularisation
3.6.1 Parcours des données en fréquence
3.6.2 Fenêtrages temporels
3.6.3 Le milieu initial
3.6.4 Fonction coût régularisée
3.6.4.1 Calcul de l’inverse de la matrice de covariance
3.6.4.2 Mise en oeuvre numérique
3.7 Conclusion
4 Imagerie à partir de données synthétiques 2D 
4.1 Choix des paramètres imagés
4.1.1 Le choix classique de V p et V s
4.1.2 Test du paramètre {ln(V p), ln(V s)}
4.1.2.1 Démarche de changement de paramètres
4.1.2.2 Cas d’un petit modèle en transmission
4.1.2.3 Fondation 4n avec un rapport fondation-encaissant R = 5 118
4.2 Résolution maximale de l’imagerie
4.3 Imagerie dans le milieu synthétique de Grenoble
4.3.1 Cadre de la FWI
4.3.2 Inversion dans le milieu sans surface libre
4.3.3 Inversion dans le milieu avec surface libre
4.3.3.1 Cas d’une acquisition avec des sources et des récepteurs en surface et sur les côtés
4.3.3.2 Cas d’une acquisition avec des récepteurs en surface et des sources en surface et sur les côtés
4.3.3.3 Cas d’une acquisition avec des sources en surface et des récepteurs en surface et sur les côtés
4.3.3.4 Cas d’une acquisition uniquement en surface
4.3.4 Choix influençant l’inversion
4.3.4.1 Cadre
4.3.4.2 Résultat
4.3.4.3 Influence de l’offset d’acquisition
4.3.4.4 Influence du maillage
4.3.4.5 Influence de l’intertrace
4.3.4.6 Influence de la gamme de fréquences
4.3.5 Migration dans le milieu de Grenoble avec surface libre
4.3.5.1 Cadre et mise en oeuvre
4.3.5.2 Conclusion
4.3.6 Imagerie avec une illumination favorable
4.3.6.1 Cadre de la FWI
4.3.6.2 Résultats de la FWI
4.3.6.3 Migration avec une acquisition réaliste et favorable
4.3.6.4 Conclusion
4.4 Imagerie dans le milieu synthétique de la maquette IFSTTAR
4.4.1 Cadre d’inversion
4.4.2 Résultats d’inversion dans le milieu avec surface libre
4.4.2.1 Inversion séquentielle
4.4.2.2 Inversion séquentielle en figeant la proche surface
4.4.2.3 Inversion multifréquencielle
4.4.2.4 Inversion avec une stratégie Bunks
4.4.3 Migration dans le milieu avec surface libre
4.4.3.1 Données
4.4.3.2 Résultats
4.5 Conclusions
5 Imagerie à partir de données réalistes 
5.1 Un milieu 3D et une modélisation 2D
5.1.1 Données synthétiques 3D
5.1.1.1 Outil de Modélisation
5.1.1.2 Cadre de modélisation
5.1.1.3 Sismogrammes
5.1.2 Données synthétiques 2D
5.1.2.1 Cadre
5.1.2.2 Sismogrammes
5.1.3 Comparaison et correction 3D-2D
5.1.3.1 Comparaison des données
5.1.3.2 Correction des données
5.1.4 Inversion des données synthétiques 3D
5.1.4.1 Démarche d’inversion
5.1.4.2 Inversion monofréquentielle de données
5.1.4.3 Inversion multifréquentielle de données
5.1.4.4 Conclusion
5.2 Imagerie à partir des données de la maquette IFSTTAR
5.2.1 Données
5.2.1.1 Données « maquette »
5.2.1.2 Données modélisées
5.2.2 Imagerie
5.2.2.1 Imagerie par inversion de la forme d’onde (FWI)
5.2.2.2 Imagerie par migration (RTM)
5.2.3 Conclusion
5.3 Imagerie à partir de données de terrain
5.3.1 Données
5.3.1.1 Données de terrain
5.3.1.2 Comparaison entre les données synthétiques et les don- nées de terrain
5.3.2 Imagerie basée sur l’inversion de la forme d’onde (FWI)
5.3.3 Imagerie basée sur la migration par retournement temporel (RTM)
5.4 Conclusion
Conclusion et Perspectives 
Bibliographie 

Télécharger le rapport complet