Imagerie sismique quantitative de la marge convergente d’Equateur-Colombie

Mèthodes d’imagerie sismique 

 Introduction 

Les mèthodes de traitement classiquement utilisèes en sismique rèflexion multitrace (NMO+stack+migration post-stack ou migration pre-stack cinèmatique) fournissent une image gèomètrique des discontinuitès du sous-sol en exploitant exclusivement les propriètès cinèmatiques des arrivèes sismiques. Cependant, une partie importante de l’information vèhiculèe par les signaux sismiques n’est pas prise en considèration : il s’agit de la modification d’amplitude et de phase de ces signaux au cours de leur propagation. Des techniques d’interprètation qualitative des dèformations du signal (modèlisation) ont ètè dèveloppèes. Ces techniques sont difficilement exploitables dans le cas de milieux très hètèrogènes. En effet, beaucoup de sections sismiques restent difficiles à interprèter dans des contextes variès (industriels ou acadèmiques) et à des èchelles très diffèrentes (imageries de subsurface et crustale, èchelles de la longueur d’onde ou mèsoscopique). Citons par exemple des problèmatiques aussi varièes que la caractèrisation des rèservoirs de pètrole ou de gaz en exploration pètrolière, l’ètude du diapirisme salifère, des BSRs (Bottom Simulating Reflector), des zones de dècollement et de la zone sismogène dans les zones de subduction. Il s’agit de zones hètèrogènes à structure et à rhèologie complexe. Les signaux enregistrès prèsentent des anomalies locales d’amplitude èvidentes. Comme ces anomalies rèsultent des caractèristiques pètrophysiques des milieux traversès, l’inversion de ces signaux devrait permettre de caractèriser plusieurs propriètès physiques des milieux traversès. On voit donc l’importance de mettre au point de mèthodes d’inversion de la forme d’onde complète (tomographie en diffraction du champ d’onde fondèes sur des modèlisations asymptotyques et numèriques des ondes) pour ètudier ces milieux de manière à pouvoir extraire le plus d’informations possibles : gèomètrie correcte des rèflecteurs et leurs propriètès physiques (i.e., variations de vitesse, d’attènuation, de densitè). Pour l’imagerie de la zone superficielle et intermèdiaire j’ai choisi d’appliquer la mèthode de tomographie en diffraction rai+Born (alias migration/inversion avant-sommation en profondeur) qui est bien adaptèe aux dispositifs d’acquisition multitrace et à l’ètude des niveaux supèrieurs de la croûte terrestre. 28 Chapitre 3. Mèthodes d’imagerie sismique 

 Bref historique de la mèthode de tomographie en diffraction rais+Born

 La mèthode de tomographie en diffraction rais+Born, ou mèthode de migration avant-sommation en amplitude prèservèe, a ètè formulèe pour la première fois par [Beylkin, 1985]. A la suite de ces travaux, [Jin et al., 1992] ont proposè un nouvel algorithme d’inversion, pour des milieux èlastiques paramètrès par les impedances P et S et la densitè. Des rèsultats sur donnèes synthètiques èlastiques ont ètè prèsentès ainsi qu’une ètude thèorique de la rèsolution du problème inverse (sèparation des paramètres, influence du dispositif d’acquisition). Parrallèlement, [Lambarè et al., 1992] prèsentent une application de la mèthode à des donnèes rèelles de sismique rèflexion multitrace acquises en mer du Nord dans l’hypothèse d’un milieu acoustique mono-paramètre, dèmontrant ainsi la faisabilitè des applications de la mèthode pour l’imagerie de milieux hètèrogènes. [Forgues, 1996] a poursuivi ces travaux en vèrifiant la mèthode sur des milieux acoustiques reprèsentès par deux paramètres (e.g. vitesse de propagation des ondes P et densitè). Il a identifiè les problèmes apparaissant lors du traitement de donnèes rèelles dans le cas d’une inversion multi-paramètres : effets de la bande passante limitèe de la source, de l’ouverture limitèe des dispositifs d’acquisition et l’imprècision du modèle de vitesse de rèfèrence sur la sèparation des diffèrents paramètres. L’extension de la mèthode aux milieux visco-acoustiques et visco-èlastiques a ètè prèsentèe d’un point de vue thèorique par [Ribodetti et al., 1998] et sa validation sur des donnèes ultrasoniques enregistrèes en cuve acoustique a ètè prèsentèe par [Ribodetti et al., 2000a]. Par la suite, la mèthode a ètè ètendue aux dispositifs d’acquisition 3D pour l’imagerie des milieux acoustiques monoparamètre. L’enjeu ètait de dèmontrer que l’application de mèthodes d’imagerie 3D ètait accessible avec des moyens de calcul disponibles dans la plupart des laboratoires (station de travail). Des algorithmes 3D rapides, utilisant le formalisme rais+Born, ont ètè dèveloppès au Centre de Recherche en Gèophysique de l’Ecole ´ des Mines de Paris et ont donnè lieu à une application pionnière en milieu acadèmique sur des donnèes 3-D du champ pètrolier d’Oseberg [Thierry et al., 1999a]. Une application au modèle synthètique OVERTHRUST 3D a ègalement ètè prèsentèe et a permis de clarifier la relation entre la gèomètrie du dispositif d’acquisition et le pouvoir de rèsolution de la mèthode de migration [Lambarè et al., 2003, Operto et al., 2003]. L’application de cette approche à des milieux fortement hètèrogènes a ègalement donnè lieu à une nouvelle extension de la mèthode prenant en compte des phènomènes de propagation complexes impliquant l’apparition de caustiques dans le champ des rais et la prise en compte d’arrivèes multiples ([Thierry et al., 1999b, Operto et al., 2000, Xu, 2001, Xu and Lambarè, 2004, Xu et al., 2004]). Mon travail s’inscrit dans la continuitè de ces travaux dans la mesure o`u il constitue un des premières applications à vocation acadèmique de ces mèthodes pour l’ètude d’une zone de subduction (marge d’Equateur-Colombie) et pour des cibles gèologiques plus profondes que celles gènèralement explorèes par la sismique pètrolière. Des outils nècessaires en amont (analyse de vitesse) et en aval (post-traitement des images migrèes) de l’algorithme de migration ont ètè dèveloppès pour disposer d’une chaîne de traitement quantitative complète dans le domaine profondeur. Un objectif majeur de cette chaîne de traitement, qui la dèmarque des approches plus conventionnelles, est d’extraire une information quantitative sur des propriètès physiques des rèflecteurs.

 La chaîne de traitement dans le domaine profondeur

 La chaîne de traitements appliquèe aux donnèes MCS est schèmatisèe sur la figure 3.1 : les donnèes sont prè-traitèes, un macro-modèle de vitesse (obtenu par analyse de vitesse, conversion en vitesse quadratique moyenne avec la formule de Dix, passage en coupes vitesse profondeur, interpolation et lissage) est construit. Un tracè de rais ([Lambarè et al., 1996]) est utilisè pour calculer tous les paramètres (temps d’arrivèe, amplitude, vecteurs lenteurs) nècessaires à la mise en oeuvre de la migration ([Thierry et al., 1999b]). Estimation du macromodèle de vitesse L’utilisation des mèthodes de migration avant-sommation profondeur (tomographie en diffraction rai+Born) n’ont de sens que si l’on dispose d’un modèle de vitesse de rèfèrence (dècrivant les variations de vitesse basse frèquence) suffisament prècis. Beaucoup d’approches ont ètè dèveloppèes pour vèrifier et amèliorer les modèles de vitesse à partir des rèsultats de la migration (”migration velocity analysis”). Pour analyser la fiabilitè de l’image obtenue (gèomètrie et amplitudes correctes des rèflecteurs), des migrations partielles à angle de diffraction constant sont calculèes [Xu, 2001]. Chaque point du milieu est èclairè par un certain nombre de trace, chaque trace dèfinissant un angle de diffraction θ. La migration de chaque sous-ensemble de donnèes associè à un angle de diffraction θ fournit une image de la structure du sous-sol dans le domaine (x − z). Les images partielles (x − z) peuvent ˆetre retrièes à distance horizontale x constante de manière à produire des panneaux profondeur-angle (z −θ) dit panneaux ISO-X (Common Images Gathers=CIG). L’analyse des panneaux ISO-X permet d’èvaluer si le modèle de vitesse de rèfèrence est prècis ([Al-Yahya, 1987]). Dans ce cas, l’image redondante d’un mˆeme rèflecteur aux diffèrents angles doit ˆetre horizontale sur le panneau iso-X. Si les panneaux presentent des rèflecteurs courbès vers le bas (ou vers le haut) (Figure 3.1 (a) en bas, exemple d’un panneau), une correction du modèle de vitesse de rèfèrence doit ˆetre appliquèe pour horizontaliser les èvènements. Une telle approche a ètè implèmentèe durant cette thèse suivant l’approche proposèe par [Al-Yahya, 1987]. Cette mèthode permet de quantifier les imprècisions du macro-modèle de vitesse via une fonction de semblance Gamma qui prèsente des valeurs supèrieurs à 1, quand les panneaux ISO-X prèsentent des rèflecteurs courbès vers le bas, (Figure 3.1 (b) en bas) ; cette fonction indique que le macro-modèle a ètè localement sous-estimè. Il est possible de calculer une fonction de correction à appliquer au macro-modèle de vitesse (Figure 3.1 (c)). On peut voir l’effet de cette correction sur un log de vitesse, qui a fait agmenter la vitesse au toit de la croûte ocèanique (Figure 3.1 (c) à droite pour un log de vitesse situè en X = 85km cf. modèle de vitesse 2D en dessus). En rèpètant ces analyses 1D à intervalle règulier le long du modèle, on peut construire par interpolation une fonction de correction 2D qui sera appliquèe au modèle de vitesse. 30 Chapitre 3. Mèthodes d’imagerie sismique Fig. 3.1 – Schèma du traitement des donnèes MCS. Dans la procèdure de migration, une mèthode de correction du modèle de vitesse v0 a ètè intègrèe. Le segment du profil SIS44 correspondant au remplissage de la fosse (sèrie stratifièe) qui repose sur le toit de la croûte ocèanique est prèsentè. Cette mèthode intègre un outil de contrˆole sur la prècision de l’image migrèe en profondeur (gèomètrie et amplitude correctes des rèflecteurs) et de correction du macro-modèle de vitesse. Voir le texte pour les dètails. On procède itèrativement jusqu’à ce que les fonctions Gamma soient ≈ 1. Des tests de validation de l’outil d’analyse de vitesse seront prèsentès dans ce chapître. La Figure 3.1 (d) en bas montre une portion du profil  (comme nous le verrons plus en dètail dans la section §6.2) avant et après correction du modèle de vitesse. On constate que les rèflecteurs sont mieux focalisès et que le toit de la croûte a ètè dèplacè d’environ 2 Km. Si des donnèes OBS coincidentes au profil MCS sont disponibles, un modèle de vitesse peut ˆetre construit par tomographie des temps de trajet des ondes rèflechies et rèfractèes (voir le chapitre sur ”Imagerie de la structure profonde : inversion conjointe de donnèes de sismique rèflexion et rèflexion/rèfraction grand-angle”). Ce modèle peut ˆetre combinè à celui obtenu par analyse de vitesse de migration pour fournir un modèle amèliorè particulièrement aux profondeurs peu rèsolues par la sismisque MCS. Pouvoir de rèsolution de la migration et dèconvolution En raison de la bande passante limitèe de la source sismique et de l’extension limitèe du dispositif d’acquisition, l’image migrèe a une rèsolution limitèe qui ne facilite pas l’interprètation quantitative des rèflecteurs : par exemple, lorsque deux rèflecteurs dèlimitant le toit et la base d’une couche sont proches, la rèsolution limitèe des images sismiques ne permet pas de sèparer ces deux rèflecteurs et d’estimer les contrastes d’impèdance acoustique au toit et à la base de la couche. Un post-traitement de l’image migrèe implèmentè sous forme d’un problème inverse par optimisation globale a ètè dèveloppè pour dèterminer une famille de modèles de vitesse impulsionnels (modèle de vitesse qui serait obtenu avec un dispositif expèrimental idèal) qui soient compatibles avec l’image migrèe. Ces modèles impulsionnels localisent les interfaces entre les couches et quantifient les vitesses de propagation des ondes P dans les couches fournissant ainsi une double information de nature structurale et lithologique. Dans la mesure o`u l’objectif de ce post-traitement est d’amèliorer la rèsolution de l’information extraite des images migrèes, nous pouvons l’assimiler à une dèconvolution. Des applications prèliminaires de ce type de post-traitement ont ètè prèsentès dans [Ribodetti et al., 2004] dans le cas d’applications à des donnèes ultrasoniques. La mise en oeuvre de ce post-traitement nècessite l’ estimation du signal source que l’on a extrait de la rèflexion sur le fond de l’eau corrigèe du coefficient de rèflexion. Le post-traitement consiste schèmatiquement à dèfinir une famille de modèles (impulsionnels) gèologiques plausibles (gènèrès alèatoirement) et de leur appliquer le filtre rèsultant de la bande passante limitèe de la source (convolution). Le modèle filtrè reprèsente une image migrèe ”synthètique” que l’on compare à l’image migrèe ”observèe” obtenue par migration rai+Born des donnèes. Le ou les meilleurs modèles gèologiques sont ceux qui fournissent la meilleur correspondance entre les 2 images migrèes. Cette dèmarche a ètè automatisèe en la formulant sous forme d’un problème inverse rèsolu par une mèthode d’optimisation globale de recuit simulè (en anglais VFSA= Very fast Simulated Annealing). Les diffèrentes phases du post-traitement de l’image migrèe sont schèmatisèes sur la Figure 3.2. Source estimation Source wavelet s (t) Random Model m (z) Random Model m (t) Time conversion * Synthetic migrated trace gt (m) Depth conversion Synthetic migrated trace gz (m) Inversion VFSA Fig. 3.2 – Schèma du post-traitement de l’image migrèe. L’image migrèe est calibrèe. L’estimation de la source est effectèe pour la convoluer avec la famille de modèles impulsionnels du sous-sol qui sont gènèrès alèatoirement. La convolution fournit une trace synthètique qui est comparèe avec la trace migrèe. Le problème de la minimisation de la fonction coût entre la trace synthètique et la trace migrèe est effectuè via un algorithme de recuit simulè (VFSA) adaptè à notre cas. Notre objectif est d’obtenir, parmi la famille de modèles impulsionnels du sous-sol ceux qui, convoluès avec la source, s’approchent le plus de la trace migrèe. La famille de modèles impulsionnels obtenue par VFSA nous reinseigne sur la gèomètrie correcte des rèflecteurs sismiques et sur leur vitesse absolue. 

 La migration/inversion Rai+Born 

Principe cinèmatique de la migration avant-sommation profondeur La migration avant-sommation profondeur a pour objectif de fournir une image en profondeur des discontinuitès du sous-sol en positionnant l’ènergie des arrivèes rèfl´chies au niveau des rèflecteurs sur lesquels elles se sont rèflèchies. On anticipe aisèment que pour atteindre cet objectif nous devons disposer d’une information sur les variations de vitesse lentes du milieu qui sont utilisèes pour ”rètropropager” l’ènergie enregistrèe aux rècepteurs au niveau du rèflecteur. On distingue deux principales familles de migration avant-sommation profondeur : les migrations avant-sommation cinèmatique qui ne modèlisent que les temps de trajet des arrivèes rèflèchies pour les positionner en profondeur. Ces migrations ne fournissent pas d’information quantifièe sur une ou plusieurs propriètès physiques du rèflecteur telles que la rèflectivitè ou les contrastes d’impèdance acoustique. La deuxième catègorie de migration avant-sommation prend en compte les aspects cinèmatique (temps de trajet) et dynamique (amplitude) de la propagation des ondes sismiques via la rèsolution de l’èquation d’onde par diffèrentes approches numèriques possibles (thèorie des rais, diffèrences finies de l’èquation d’onde paraxiale (one-way) ou complète (two-way)). En formulant la migration dans le cadre de la thèorie des problèmes inverses, des informations sur les propriètès physiques des rèflecteurs, qui agissent sur l’amplitude et la forme des ondes rèflèchies, peuvent ˆetre obtenues. C’est cette deuxième catègorie de mèthodes appellèes migration/inversion ou migration en amplitude prèservèe qui est utilisèe dans cette thèse. Avant de dèvelopper les fondements thèoriques de la mèthode de migration/inversion rai+Born, je rappelle ci-après le principe d’imagerie cinèmatique sur lequel repose les algorithmes populaires de migration avant-sommation tels que la mèthode de Kirchhoff. Considèrons tout d’abord le problème direct qui consiste à calculer le temps d’arrivèe d’un trajet diffractè entre une source s et un rècepteur r via un point diffractant x. Pour cela, on calcule le temps entre la source et le diffractant tsx et le temps entre le point diffractant et le rècepteur txr. La somme de ces deux temps fournit le temps d’arrivèe de la diffraction (s − x − r) (Figure 3.3). Fig. 3.3 – Problème direct du problème de migration. Le temps d’arrivèe d’une diffraction (s − x − r) est calculè en additionnant le temps du trajet incident (s − x) plus celui du trajet diffractè (x − r). Considèrons maintenant le cas d’un rèflecteur continu. L’hodochrone (courbe temps- 34 Chapitre 3. Mèthodes d’imagerie sismique offset) de l’arrivèe rèflèchie sur ce rèflecteur peut ˆetre obtenue par l’enveloppe des diffractions gènèrèes par des points diffractants règulièrement rèpartis le long de l’interface. Ce rèsultat est une consèquence du principe de Huygens (Figure 3.4). Considèrons maintenant le problème inverse : des arrivèes rèflèchies sont observèes dans les donnèes et nous voulons reconstruire l’image a priori inconnue des rèflecteurs dont elles sont issues en faisant l’hypothèse que l’on dispose d’un macromodèle de vitesse. L’amplitude d’une arrivèe rèflèchie enregistrèe à un temps t va ˆetre ”ètalèe” en tous les points x de l’espace vèrifiant la condition d’imagerie : t = tsx + txr. Dans le cas d’un macromodèle de vitesse homogène, ces points de l’espace vèrifient l’èquation d’une ellipse appellèe isochrone. Dans l’approximation frèquence infinie, la migration d’une trace (qui Fig. 3.4 – L’hodochrone d’une rèflexion est formèe par l’enveloppe des diffractions èlèmentaires issues des diffractants discrètisant le rèflecteur. peut ˆetre vu comme un èchantillon de l’hodochrone de l’arrivèe rèflèchie de la mˆeme manière que le point diffractant ètait un èchantillon du rèflecteur continu dans le cadre du problème direct) va imager un point du rèflecteur. L’image continue du rèflecteur va ˆetre fournie par la superposition des images partielles fournies par la migration de chaque trace (Figure 3.5b). Cela relève encore du principe de Huygens appliquè au problème inverse cette fois. La couverture multiple fournie par les dispositifs de sismique rèflexion mulitrace va permettre d’additionner constructivement une information redondante provenant de plusieurs couples tir-capteur au niveau du rèflecteur tandis que l’ènergie ètalèe le long des isochrones à l’extèrieur des rèflecteurs va ˆetre additionnèe de manière destructive au cours de la sommation. Fig. 3.5 – Problème inverse. a) L’ènergie d’un èvènement est distribuèe en tout point de l’espace vèrifiant la condition d’imagerie t = tsx + txr. Ces points de l’espace dècrivent une ellipse appellèe isochrone ou courbes iso-temps. Des isochrones sont reprèsentèes pour diffèrents offsets source-rècepteur. b) L’image d’un rèflecteur continue est formèe par l’enveloppe de toutes les isochrones. La sommation est constructive sur le rèflecteur et destructive à l’extèrieur de celui-ci. Plus l’information redondate sommèe est volumineuse, meilleur sera le rapport signal sur bruit de l’image migrèe. De manière schèmatique, un algorithme de migration avant-sommation profondeur contient les ètapes èlèmentaires suivantes : 1) Pour chaque tir, calcul des temps de trajet entre le tir et chaque point de l’image migrèe, tsx. 2) Pour chaque capteur associè à ce tir, calcul des temps de trajet entre le capteur et chaque point de l’image migrèe, txr. 3) Pour chaque point x de l’image migrèe, 3.1. calcul de t = tsx + txr. 3.2. Extraction de l’amplitude a de la trace au temps t 3.3. addition de l’amplitude a pondèrèe de corrections dynamiques w à l’image migrèe m : m(x)=m(x)+w.a De nombreuses stratègies numèriques fondèes essentiellement sur des interpolations des grandeurs modèlisèes ont ètè dèveloppèes pour limiter le nombre d’opèrations effectuèes dans les algorithmes de migration [Thierry et al., 1999b, Thierry et al., 1999a].