FR 2784195 A1 20000407 FR 9812337 A 19981001 La présente invention concerne une méthode pour réaliser une migration avant sommation d'événements sismiques enregistrés pour l'imagerie d'une partie d'une zone souterraine. La méthode selon l'invention permet par exemple de réaliser une migration profondeur de type 3D avant sommation (" 3D prestack depth migration "), pour un modèle de vitesses donné, de manière à imager les interfaces géologiques ou hétérogénéités diverses d'une partie de la subsurface. Art antérieur La migration avant sommation est une méthode classique de traitement de o0 données sismiques. De façon générale la technique consiste, connaissant la valeur d'un champ d'ondes à une profondeur connue, par exemple en surface, ainsi qu'un modèle de la répartition des vitesses de propagation des ondes dans la zone, à modéliser la propagation du champ source et la retropropagation des données de réflexion enregistrées et à chercher des cohérences de phase entre ces deux champs modélisés. On distingue trois grands types de migration avant sommation - la migration par point de tir: le champ source est l'état vibratoire engendré par le point de tir et les données de réflexion sont la réponse de la subsurface à ce champ source; - la migration par ondes planes aussi appelée migration par angle d'éclairage commun: le champ source est l'onde plane considérée et les données de réflexion sont la réponse de la subsurface à ce champ source - la migration par déport (" offset "): le champ source est celui émis par un point de tir et les données de réflexion sont les enregistrements réalisés par le (ou les) capteur(s) associé(s) à ce point de tir ayant l'offset considéré; dans une telle migration il faut, pour réaliser la migration des données associées à un déport, effectuer autant de modélisation de propagation et de retropropagation d'ondes qu'il y a de points de tir et sommer les résultats obtenus pour chaque point de tir. Des exemples d'utilisation de ce type de techniques sont décrits par exemple dans: - Claerbout, J.F., 1985; Imaging the Earth's interior; Blackwell Publications, -Duquet, B., 1996; Amélioration de l'imagerie Sismique de Structures Géologiques Complexes; thèse, Université Paris 13, ou - Whitmore, N.D., Felinsky, W.F., Murphy, G.E. and Gray, S.H., 1993; The Application of Common Offset and Common Angle Pre-stack Depth Migration in the North Sea, 55th Mtg., EAGE, Expanded abstract. Les implémentations classiques basées sur l'intégrale de Kirchhoff (ou des versions plus élaborées de cette technique, elle-même basée sur des techniques asymptotiques hautes fréquences), ont pour inconvénient majeur d'être généralement très coûteuses en temps de calcul, en raison du volume aussi bien des données à traiter que des résultats, surtout quand le champ de vitesses varie latéralement (ce qui complique les calculs de temps d'arrivée requis pour la mise en oeuvre de cette méthode). Pour des raisons économiques, on est souvent amené à limiter les volumes de données, (par décimation) et/ou la quantité de résultats produits (volume imagé de taille réduite, échantillonnage grossier des résultats). Définition de la méthode La méthode selon l'invention permet de réaliser la migration de données sismiques pour l'imagerie d'une partie d'une zone souterraine, les données sismiques étant obtenues à l'issue d'une série de Ns cycles de sismique réflexion, qui comprend chacun l'émission successive de champs d'onde élémentaires défini chacun par l'association d'un signal sismique W(t) et d'un lieu d'émission dans une série de lieux d'émission Si avec 1 modélisés en tout point de la zone à imager. f) On répète les étapes c) à e) pour toutes les valeurs prises par les composantes pi et p2 du dit vecteur p; et -.-) g) pour toute valeur fixée de la deuxième composante P2 du dit vecteur p, on somme les résultats de ces différentes combinaisons de manière à obtenir une image migrée associée à cette valeur fixée de p2, réalisant ainsi une migration avant sommation. Suivant un mode de mise en oeuvre, on peut sommer les résultats obtenus en g) pour toutes les valeurs prises par le paramètre p2. On effectue ainsi un sommation après migration. Suivant un mode de mise en oeuvre, il est possible de réaliser directement le sommation après migration, sans particulariser l'étape g). La méthode peut comporter en outre une mise a jour des vitesses par analyse des déformations obtenues lorsque l'on fait varier la deuxième coordonnée P2 du -. vecteur p Suivant un mode de mise en oeuvre, on peut former une image migrée d'une partie de la zone à imager en exploitant le phénomène de conversion d'ondes, par définition d'au moins une partie du champ de vitesse en ondes P et en ondes S (en appliquant par exemple au préalable un pré-traitement adapté aux données, visant à séparer les différents types d'événements sismique). On peut utiliser les étapes a) à g) pour la détermination du gradient d'une fonction coût intervenant dans un problème inverse de sismique. Il est possible également de remplacer une migration profondeur par une migration temps. La méthode proposée présente de grands avantages 1) Elle permet de réaliser tune migration à un prix (coût de calcul) avantageux car indépendant du volume de résultats calculé et du nombre de traces sismiques enregistrées, contrairement aux méthodes classiques de type Kirchhoff. Seul le volume de la zone dans lequel les ondes se propagent, ayant une incidence sur le coût de calcul. On obtient ainsi des images volumiques en prenant en compte l'ensemble des traces sismiques à un coût avantageux. On estime que cette méthode diminue 0o d'un facteur de l'ordre de quelques dizaines le temps de calcul nécessaire à la réalisation de la dite migration 3D avant sommation. 2) La méthode s'applique pour des modèles de vitesses d'une complexité arbitraire pour peu que la notion de migration avant sommation conserve un sens. Elle s'applique sans rencontrer aucune des limitations propres aux techniques asymptotiques hautes fréquences (optique géométrique) habituellement utilisées pour réaliser les migrations 3D avant sommation. La méthode peut être mise en oeuvre à l'aide d'outils classiques de modélisation de la propagation et de la retropropagation des ondes, que l'on trouve décrits par exemple dans la référence Duquet B. déjà citée. Par application de la méthode, on peut obtenir des images migrées élémentaires associées à une valeur donnée d'un paramètre et la somme de ces images (" post migration stack "), et ceci aussi bien dans le domaine profondeur que dans le domaine temps. D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après d'un exemple non limitatif de réalisation, en se référant aux dessins annexés o - la Fig.I montre de façon très schématique, une disposition de lieux d'émission S, et de lieux de réception sismique R,J d'une zone souterraine explorée o l'on cherche à restituer des horizons In de la subsurface; - la Fig.2 montre une section sismique extraite d'un bloc de données 3D migrées profondeur après sommation, obtenue dans le cadre de l'exploration sismique d'une structure salifère de mer du Nord; - la Fig.3 montre une section sismique, coïncidant avec la section sismique de la Fig. 2, extraite d'un bloc de données 3D migrées profondeur avant sommation par la méthode selon l'invention, en utilisant le modèle de vitesse utilisé pour l'obtention de la Fig. 2; et - la Fig.4 montre en quatre points différents de la surface de la zone souterraine, les images 3D migrées en profondeur avant sommation, par la méthode selon l'invention lorsque la deuxième composante P2 du vecteur p varie sur son domaine de définition, l'analyse de la déformation des événements observés sur ces images permettant la mise à jour du modèle de vitesse. DESCRIPTION DETAILLEE DE LA MIETHODE Définitions, notations et hypothèses requises La méthode selon l'invention permet de réaliser la migration de données sismiques pour l'imagerie d'une zone souterraine M, les données sismiques étant obtenues à l'issue d'une série de Ns cycles sismiques, chacun d'eux (Fig. 1) comportant l'émission d'un signal sismique W(t) depuis un lieu d'émission Si avec l données, tel qu'une déconvolution. On désignera par profil un ensemble de lieux d'émission alignés et l'on supposera que toutes les acquisitions de signaux sont effectuées en utilisant des profils d'émission parallèles. On supposera que, pour chaque source, les récepteurs -'4 sont localisés sur le profil o elle se trouve (repérée par le lieu d'émission Si). En fait, la méthode s'avère robuste vis à vis de la précision avec laquelle ces hypothèses sont vérifiées. Les hypothèses précédentes sont faites principalement pour pouvoir introduire les notations ci-dessous. On définit un repère orthonormé dont le premier axe est perpendiculaire à la direction du profil et le deuxième axe parallèle à cette direction et un vecteur un -4 vecteur lenteur p (homogène à l'inverse d'une vitesse, comme il est bien connu). Dont les composantes pi et P2 mesurent respectivement les composantes suivant ces -. -4 deux axes de ce vecteur lenteur p. Pour un vecteur lenteur p donné et une source -4 placée en un lieu d'émission Si donné, on définit une fonction de décalage temporel -..4 -"9 -'9 (1) to(p, i)= p.(Si -So) o So représente un point quelconque du domaine d'acquisition. Traitement Sélection du vecteur p On commence par choisir un ensemble fini Pl de valeurs prises par le paramètre pi et un ensemble fini P2 de valeurs prises par le paramètre p2. L'ensemble E=PxP2 sera l'ensemble des valeurs (vectorielles) prises par le vecteur p.-4 E=PlxP2 sera l'ensemble des valeurs (vectorielles) prises par le vecteur p. L'ensemble P1 des valeurs que prendra le paramètre pi pourra par exemple être construit en échantillonnant l'intervalle [-p1min, pimax] avec un pas d'échantillonnage régulier Spi. La valeur à donner au pas Spidépend notamment de la précision souhaitée et de l'espacement entre les profils. Une valeur typique est: 8p = 2.5. 10-5 s/m. Les valeurs à donner à p1min et pimax dépendent de la complexité de la structure dans la direction orthogonale au profil. Des valeurs typiques sont: - plmin = pimax = 2.5. 10-4 s/m. L'ensemble P2 des valeurs que prendra le paramètre P2 pourra par exemple être construit en échantillonnant l'intervalle [-p2min, p2max] avec un pas d'échantillonnage régulier 5p2. La valeur à donner à ap2 dépend i0 notamment de la précision souhaitée et de la finesse avec laquelle on veut suivre l'évolution, lorsque p2 varie, des événements dans les collections point image. Une valeur typique est: 5p2 =2.5. 10-5 s/m. Les valeurs à donner à p2min et p2max dépendent de la complexité de la structure dans la direction des profils. Des valeurs typiques sont: -p2min = p2max = 2.5. 10-4 s/m A) Etapes de traitement: 1) Création d'un premier champ d'ondes composite en surface (champ de sources) L'ondelette W(t) attachée à chacune source Si définit un champ d'onde élémentaire W (t).6( x 2D - Si) o le vecteur x 2D) repère la position d'un point quelconque à la surface du sol et 8 est une masse de Dirac définie sur R2 et centrée à l'origine. Après avoir retardé du temps to( p, i) le champ d'onde élémentaire associé, on forme un premier champ d'ondes composite en surface par superposition spatio- temporelle des différentes champs d'onde élémentaires ainsi décalées. On définit ainsi la fonction: (2) W X (x 2D, t) = Y ( X 2D - Si) W (t - to( p, i)) pi 2) Création d'un deuxième champ d'ondes composite en surface (champ de traces) Apres avoir retardé du temps to( p, i) la trace sismique d (t) associée à chaque couple (lieu d'émission Si, lieu de réception R J) repéré par le couple (i,j), i on forme un deuxième champ d'ondes composite en surface par superposition spatio- temporelle des différentes traces sismiques ainsi décalées. Avec les mêmes notations, on définit ainsi la fonction (3) Do (x2D,t)= x 2D - R) d (tto( p, i)) P. i,j j i Io 3) Modélisation de la propagation du champ source composite) On modélise la propagation du champ source composite en recherchant des solutions périodiques de période T de l'équation des ondes, en utilisant comme distribution de vitesses celle définie par le modèle de vitesses considéré. On obtient ainsi un premier champ d'ondes propagé (dépendant de l'espace et du temps) W- ( x 3D, t), ceci pour tout temps t (la solution est périodique en temps) et pour P. tout point de la partie de subsurface que l'on cherche à imager, point dont la position est repérée par le vecteur x 3D. La période T est choisie comme il est habituel dans les algorithmes classiques de migration (c'est à dire de l'ordre de la durée d'enregistrement). 4) Modélisation de la rétropropagation du cham de traces composite On modélise la rétropropagation du champ de traces composite n recherchant des solutions périodiques de période T de l'équation des ondes en utilisant comme distribution de vitesses celle définie par le modèle de vitesses considéré. On obtient ainsi un champ d'ondes rétropropagé (dépendant de l'espace et du temps) D- (x 3D, p. t), ceci pour tout temps t et pour tout point de la partie de la subsurface que l'on cherche à imager, point dont la position est repérée par le vecteur x 3D- ) Recherche de cohérence de phase On recherche alors une éventuelle cohérence de phase (par exemple par des calculs d'intercorrélation) entre le premier champ d'ondes composite propagé et le deuxième champ de traces composite rétropropagé, ceci en tout point de la zone souterraine (subsurface) que l'on cherche à imager, point dont la position est repérée par le vecteur x 3D. Pour ce faire, on évalue, toujours dans le cas d'un calcul par intercorrélation, la quantité: - T --, - (4) M- (X3D)= | W--- ( X 3D, t) D- (x3D, t) dt. P- o p P. La composante p2 ayant une valeur donnée dans P2, on somme les résultats calculés en (4) lorsque le paramètre pi parcourt Pl. Cette sommation permet de définir, pour tout point image repéré par le vecteur x 3D, la fonction Mp2 ( x 3D) par la formule: --4 --- (5) M p2( x 3D) = M- ( X 3D) pi e Pi P, La quantité (5) s'interprète comme la valeur au point x 3D de la superposition sur les différents profils d'acquisition des images migrées associées à une onde cylindrique ayant pour axe le profil d'acquisition et une pente définie par la lenteur p2. B) Sommation après migration I est possible en outre de réaliser une sommation après migration en sommant les contributions (5) obtenues pour les différentes valeurs du paramètre P2 appartenant à P2, selon la formule: (6) M(X3D)= E Mp2(X3D) P. E P2 VARIANTES 1) Une première variante des étapes A-3 et A-4 décrites ci-dessus consiste à résoudre à l'étape A-3 un problème de valeurs initiales et en A-4 un problème de valeurs finales au lieu de rechercher des solutions périodiques (procédures d'imagerie classiques connues sous le nom de migration par renversement du temps). 2) Une autre variante consiste à remplacer l'étape A-1 de création du premier champ d'ondes composite par la création d'une onde plane dont les lenteurs suivant la direction des profils et suivant la direction orthogonale sont respectivement p2 et pl. On aboutit ainsi à la formule: (7) W (X3D) = W (t- p.(x 2D-So)) P Cette variante met en évidence une certaine similitude des traitements mis en oeuvre dans la méthode et de ceux mis en oeuvre dans la méthode classique de migration par ondes planes. Une différence essentielle qui fait l'originalité de la méthode réside dans la possibilité de réaliser une procédure de migration par ondes planes même dans le cas o l'acquisition ne permet pas la synthèse de la réponse de la subsurface à une excitation ondes planes, réponse qu'il est indispensable de connaître dans les algorithmes connus de migration par ondes planes. Autres applications Suivant un mode de mise en oeuvre, il est possible de former une image migrée d'une partie de la subsurface en exploitant le phénomène de conversion d'ondes, par définition d'au moins une partie du champ de vitesses en ondes P et en ondes S, ceci après éventuellement appliqué un prétraitement adapté aux données, io visant à séparer les différents types d'événements sismiques. Les étapes précédemment définies peuvent aussi servir pour le calcul du gradient d'une fonction coût intervenant dans un problème inverse de sismique. REVENDICATIONS 1) Méthode pour réaliser une migration avant sommation d'événements sismiques pour l'imagerie d'une partie d'une zone souterraine à partir d'une série d'un nombre Ns de cycles de sismique réflexion, qui comprend l'émission successive de champs d'onde élémentaires défini chacun par l'association d'un signal sismique W(t) et d'un lieu d'émission défini dans une série de lieux d'émission Si avec i modélisés en tout point de la zone à imager. f) On répète les étapes c) à e) pour toutes les valeurs prises par les composantes pi et p2 du dit vecteur p; et g) pour toute valeur fixée de la deuxième composante p2 du dit vecteur p, on somme les résultats de ces différentes combinaisons de manière à obtenir une image migrée associée à cette valeur fixée de p2, réalisant ainsi une migration avant sommation. 2) Méthode selon la revendication précédente, caractérisée en ce qu'elle comporte en outre la formation d'une image finale représentative d'une sommation après migration des traces obtenues pour l'ensemble des valeurs de la deuxième -..4 composante P2 du vecteur p 3) Méthode selon la revendication 1, caractérisée en ce que l'on utilise les étapes a) à g) pour réaliser directement un sommation après migration. 4) Méthode selon la revendication 1, caractérisée en ce qu'elle comporte en outre une mise ajour des vitesses par analyse des déformations obtenues lorsque l'on fait varier la deuxième coordonnée P2 du vecteur p. fait varier la deuxième coordonnée p2 du vecteur p. ) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on forme une image migrée d'une partie de la zone à imager en exploitant le phénomène de conversion d'ondes, par définition d'au moins une partie du champ de vitesse en ondes P et en ondes S. 6) Méthode selon la revendication 5, caractérisé en ce que l'on applique un pré-traitement adapté aux données, visant à séparer les différents types d'événements sismiques. 7) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on utilise les dites étapes a) à g) pour la détermination du gradient d'une fonction io coût intervenant dans un problème inverse de sismique. 8) Méthode selon l'une des revendications précédentes, caractérisée par le fait que l'on remplace une migration profondeur par une migration temps.