Note: Descriptions are shown in the official language in which they were submitted.
CA 02284386 1999-09-29
1
METHODE POUR RÉALISER EN 3D AVANT SOMMATION,
UNE MIGRATION DE DONNÉES SISMIQUES
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
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
nûgration il faut, pour
réaliser la migration des données associées à un déport, effectuer autant de
modélisation
CA 02284386 1999-09-29
la
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,
CA 02284386 1999-09-29
2
- 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,
1o surtout quand le champ de vitesses varie latéralement (ce qui complique les
calculs de
temps d'arrivée requis pour la mise en ceuvre 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).
CA 02284386 1999-09-29
3
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:5 i<_
Ns, la réception par des récepteurs sismiques placés en des positions des
signaux
sismiques renvoyés par la zone en réponse à chacun de ces champs d'onde, et
l'enregistrement des différents signaux reçus par chaque récepteur sismique
sous la forme
l0 de traces sismiques d(t) dépendant du temps.
Elle est caractérisée en ce que, pour un modèle de vitesses donné, elle
comporte
l'ensemble des étapes suivantes
~
a) on définit un vecteur de lenteur p dont les deux composantes pi et p2
peuvent
chacune prendre une suite de valeurs préalablement définie ;
-~ -~
b) on définit pour un vecteur de lenteur p donné et lieu d'émission Si donné,
une
-->
fonction de décalage temporel to ( p , i)
~
c) on applique la fonction de décalage temporel to( p, i) à chaque champ
d'onde
~
élémentaire (associé au lieu d'émission Si ) et on forme un premier champ
d'ondes
composite en surface par superposition spatio-temporelle des différents champs
d'onde
élémentaires ainsi décalés ;
d) on applique un décalage temporel to ( p, i) à chaque trace sismique d~(t)
repérée
i
par le couple (i,j) et on forme un champ de traces composite en surface par
superposition
spatio-temporelle des différentes traces sismiques ainsi décalées
CA 02284386 1999-09-29
4
e) on effectue une migration du champ de traces composite en utilisant comme
champ d'ondes le champ d'ondes composite, ceci en modélisant la propagation de
champ
d'ondes composite ainsi que la retropropagation de champ de traces composite,
et en
combinant de façon adaptée les deux champs composites ainsi 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 pl
et p2 du dit vecteur p ; et
~
g) pour toute valeur fixée de la deuxième composante pZ 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 pZ, 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 pZ. On effectue ainsi un
sommation après
migration.
Suivant un mode de mise en ceuvre, 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.
CA 02284386 1999-09-29
La méthode proposée présente de grands avantages :
1) Elle permet de réaliser une 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
5 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 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 ceuvre à 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.1 montre de façon très schématique, une disposition de lieux
d'émission S; et de
lieux de réception sismique R;' 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 ;
CA 02284386 1999-09-29
6
- 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 METHODE
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 1:9 i<_ Ns, la réception par
des récepteurs
sismiques placés en des positions des signaux sismiques renvoyés par les
discontinuités de la zone et l'enregistrement des différents signaux reçus par
chaque
récepteur sismique sous la forme d'une trace sismique d(t). On admet que
chaque source
sismique émet le même signal W(t), car on peut toujours se ramener à une telle
situation
par un pré-traitement adapté des 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 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 vecteur
CA 02284386 1999-09-29
7
~
lenteur p(homogène à l'inverse d'une vitesse, comme il est bien connu). Dont
les
composantes pl et P2 mesurent respectivement les composantes suivant ces deux
axes de ce
vecteur lenteur p . Pour un vecteur lenteur p donné et une source placée en un
lieu
~
d'émission Si donné, on définit une fonction de décalage temporel :
-~ -~ -3 -~
(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
pl
lo et un ensemble fini P2 de valeurs prises par le paramètre P2. L'ensemble
E=P1xP2 sera
~
l'ensemble des valeurs (vectorielles) prises par le vecteur p. L'ensemble P1
des valeurs
que prendra le paramètre pl pourra par exemple être construit en
échantillonnant
l'intervalle [-plmin, plmax] avec un pas d'échantillonnage régulier Spl. La
valeur à donner
au pas Spldépend notamment de la précision souhaitée et de l'espacement entre
les profils.
Une valeur typique est : Spl = 2.5 . 10-5 s/m. Les valeurs à donner à plmin et
pimax
dépendent de la complexité de la structure dans la direction orthogonale au
profil. Des
valeurs typiques sont :-plmin = plmax = 2.5 . 10-4 s/m. L'ensemble P2 des
valeurs que
prendra le paramètre pz pourra par exemple être construit en échantillonnant
l'intervalle [-
pZmin, p2max] avec un pas d'échantillonnage régulier Spz. La valeur à donner à
Sp2 dépend
2o notamment de la précision souhaitée et de la finesse avec laquelle on veut
suivre
l'évolution, lorsque pZ varie, des événements dans les collections point
image. Une valeur
typique est : 8p2 =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
= pzmax = 2.5 . 10-4 s/m
A) Etapes de traitement :
CA 02284386 1999-09-29
8
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).S( x"p - Si ) où le vecteur x 2D repère la position d'un point quelconque
à la surface
du sol et S est une masse de Dirac définie sur R' 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) WP (x2D,t)=~S(x2D-Si)W(t-to(P, i))
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~ i (t) associée à
chaque
couple (lieu d'émission Si, lieu de réception R~) repéré par le couple (i,j),
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) D' ( x 2D, t) _~ S( x 2D - R 1) d (t - ta( P+ i))
P. J i
3) Modélisation de la propaaation 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 prenver
--~
champ d'ondes propagé (dépendant de l'espace et du temps) W-~ ( x 3D, t), ceci
pour tout
P.
temps t (la solution est périodique en temps) et pour 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
CA 02284386 1999-09-29
9
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 champ 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, t), ceci
pour tout
P.
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=
5) 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-4 ( X 3D) = J W~ ( X 3D, t) D-' ( X 3D, t) dt.
P. 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 Pt. 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
:
(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.
CA 02284386 1999-09-29
B) Sommation après migration
Il 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 :
5 (6) M(X 3D) = E M p2 ( X 3D)
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
1o 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.(X2D-So))
p
Cette variante met en évidence une certaine similitude des traitements mis en
oyuvre
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 ceuvre, 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
CA 02284386 1999-09-29
11
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, 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.
