FR 2902196 A1 20071214 FR 0605123 A 20060609 PROCEDE DE TRAITEMENT SISMIQUE DU TYPE METTANT EN' îUVRE UNE MIGRATION PROFONDEUR AVANT ADDITION La présente invention est relative à un procédé de traitement sismique du type 5 mettant en oeuvre une migration profondeur avant addition. Plus précisément, elle concerne un procédé de ce type mettant en oeuvre une étape préalable de modulation des points de tirs. 10 La migration profondeur avant addition est une étape centrale du traitement sismique. Elle consiste à focaliser les événements sismiques enregistrés en temps, en des réflexions indexées en profondeur. Une des manières les plus précises d'effectuer cette étape est la migration par point de tir, qui consiste à propager numériquement une onde incidente représentant la source sismique et une onde réfléchie. On 15 initialise en surface à la profondeur z=0, et pour chaque fréquence f, l'onde descendante dn(f,x,y,z) par une représentation synthétique de la source et l'onde montante un(fx,y,z) par l'onde enregistrée par les capteurs sismiques. La propagation numérique propage ces ondes de proche en proche à travers des couches d'épaisseur Az permettent d'obtenir les ondes descendantes et montantes dn(f,x,y,z) et un(fx,y,z) 20 pour toute profondeur z sur la grille. Le résultat de la migration est une réflectivité calculée en sommant sur toutes les fréquences et tous les tirs, la cross-corrélation des ondes descendantes et montantes: (1) r(x,y,z)= Edn(.f,x,y,z)un(.f,x,y,z) f,n L'inconvénient de cette méthode est que l'ensemble des points de tirs de 25 l'acquisition sismique un(f,x,y,z) est représenté par une seule image r(x,y,z). On souhaite en général avoir un résultat de migration qui soit une collection rt(x,y,z) telle que : (2) r(x, y, z) = E ri (x, y, z) Ceci permet de faire des analyses vitesses en vérifiant que les arrivées sismiques ont le même temps d'arrivée sur toutes les collections. On peut aussi faire des traitements pour atténuer les arrivées multiples en rejetant les événements présentant une courbure dans le plan profondeur-i de la collection. ETAT DE LA TECHNIQUE La migration par point de tir permet de produire une collection par point de tir où rn(x,y,z) est la réflectivité fournie par le point de tir n. Mais ceci n'est pas commode au vu du grand nombre de points de tirs (105 à 106) d'une acquisition moderne. En général, on souhaite des collections comportant quelques dizaines de composantes. La migration par onde plane (Schultz et Claerbout (1978), Rietveld et al (1992), Duquet et al (2001)) est une méthode où des tirs composés sont construits à partir des tirs individuels, chaque tir composé étant obtenus en sommant les tirs individuels après application de retards linéaires proportionnels à une valeur de lenteur p donnée (la lenteur est l'inverse de la vitesse), ou de façon équivalente en construisant le tir composé d'indice m en pondérant dans le domaine fréquentiel le tir d'indice n se trouvant à l'abscisse xä par C,,,,,(f)=exp(-2jitfp,nxn). Cette superposition linéaire des ondes descendantes et montantes est suivie d'une migration analogue à la migration par point de tir, les points de tir composés remplaçant les points de tirs individuels. Cette technique a deux inconvénients : elle ne permet pas de retrouver exactement la même image que la migration par point de tir et ne donne pas de critère sur le nombre de valeurs de p à migrer. Romero et al (2000) utilisent des tirs composés obtenus à partir d'autres pondérations qui ne dépendent pas de la position de la source xn (phases linéaires, aléatoires, chirp, chirp modifié). Ces méthodes souffrent des défauts des méthodes précédentes et en outre ne permettent pas d'obtenir des collection interprétables. BUTS ET PRESENTATION GENERALE DE L'INVENTION Un but de l'invention est de proposer un procédé de traitement avec migration basé 30 sur une composition de tirs permettant d'obtenir exactement la même image qu'une migration par points de tir. Un autre but de l'invention est de proposer un procédé de traitement de ce type avec un nombre de tirs composés à migrer minimal. Un autre but encore de l'invention est de proposer un procédé de traitement avec migration permettant d'obtenir des collections analogues à celles obtenues dans le cas d'une migration par onde plane. L'invention propose notamment un procédé de traitement sismique du type par migration, selon lequel : on combine linéairement des ondes descendantes émises par différentes sources sismiques en différents points de tirs, ainsi que des ondes montantes enregistrées par différents capteurs sismiques, - on propage fictivement les ondes composées ainsi obtenues afin d'obtenir des ondes descendantes et montantes migrées pour différentes profondeurs, - on détermine en différentes profondeurs au moins une caractéristique du sous-sol en fonction des ondes montantes et descendantes ainsi propagées ; caractérisé en ce que les ondes composées descendantes (respectivement montantes) sont calculées par une combinaison linéaire des ondes descendantes (resp. montantes) selon laquelle lesdites ondes descendantes (resp. montantes) sont pondérées par les coefficients d'une matrice de modulation spatiale fonction de la position des points de tirs. DESCRIPTION DES FIGURES - la figure 1 représente dans un plan f-k un exemple de tirs composés à migrer dans le 25 cas d'une migration par modulation des points de tirs conforme à un mode de mise en oeuvre de l'invention ; - la figure 2: représente dans un plan f-k un exemple de tirs composés à migrer pour une migration par ondes planes ; - les figures 3 et 4 représentent les images fournies par la migration par modulations 30 de tirs (figure 3) et pour la collection (figure 4), dans le cas d'un jeu de données synthétiques dit modèle Sigsbee et d'un mode de mise en oeuvre de l'invention ; - la figure 5 illustre schématiquement les étapes d'un exemple de mise en oeuvre de l'invention. DESCRIPTION D'UN OU PLUSIEURS MODES DE MISE EN îUVRE DE L'INVENTION Principes généraux 10 Condition d'exactitude d'une migration par point de tir composé: La migration par point de tir étant la plus exacte des migrations, on peut appeler condition d'exactitude une condition garantissant que l'image obtenue par une autre migration est la même que celle obtenue par une migration par point de tir. Nous allons obtenir une telle condition pour une migration par tir composé. 15 dn(f,x,y,z=0) and un(f,x,y,z=0), n=O,N-1 étant les ondes descendantes et montantes à la profondeur z=0 correspondant aux N tirs individuels d'une acquisition, la migration par point de tir composé définit une matrice de composition C(f) de dimension (M,N) permettant d'obtenir les ondes composées descendantes et montantes Dm(f,x,y,z=0) et Um(f,,x,y,z=0): Nù1 20 (3) Um( f , x, y, z = 0) = ECmn(f)un(f,x,y,z=0),m=O,M-1 n=0 la même relation linéaire reliant Dm à dn. On effectue ensuite une migration par point de tir de ces deux ondes composées. Pour cela, les deux ondes sont extrapolées en profondeur et la linéarité de l'extrapolation garantit que la relation linéaire (3) est valable pour toute profondeur z, entre Um(f,x,y,z) et un(f,x,y,z), et de manière similaire 25 entre Dm(f x,y,z) et dn(f,x,y,z). Ceci peut s'écrire en notation matricielle, en utilisant les vecteurs d(f)=[dI(,,x,y,z),••,dN I(f,x,y,z)]T et D(f)=[Dl(f,x,y,z),••,DM-I(f,x,y,z)]T et la matrice C(f) d'éléments Cmn(f): (4) D(f) = C(f)d(f),U(f) = C(f)u(.f) L'image obtenue par la migration par tir composé est:5 (5) M -1 R(x,y,z)= 1Dm(f,x,Y,z)Um(f,x,y,z)=1D*(f)U(f)=>d*(f)C*(f)C(f f n=0 f f Si la matrice C(f) est unitaire pour tout f, alors C*(f)C(f) (6) R(x, y,z) = Id*(f)C*(f)C(f)u(f) = Id*(f)u(f) = r(x,Y,z)=Id f f ce qui signifie qu'une migration par tir composé avec une matrice de composition unitaire produit la même image que la migration par point de tir. 10 Migration par modulations de points de tirs : On considère une migration par point de tir composé, où le nombre de tirs composés M est égal au nombre de tirs d'origine N et où la matrice de composition est , pour toute les fréquences, la matrice de la transformée de Fourier discrète spatiale: 1 mn \ 15 (7) Cmn (f) = Cmn = ,/N exp 2 jrc N ) Cette matrice (qui ne dépend pas dei) est unitaire de sorte que la migration fournit exactement la même image que la migration par point de tir. L'égalité des images découle directement du théorème de Parseval. En règle générale, les positions des tirs sont espacés régulièrement sur une grille de 20 pasdx. xn ndx est alors la position du tir d'indice n. Comme il y a N tirs, la longueur totale de l'ensemble des tirs est xm NAx. On définit le pas en nombre d'onde Ak =1/xmax, qui permet de définir une grille de nombre d'ondes km mAk. Avec ces notations on obtient : 1 (8) Cmn = ,/N exp(2Jrkmxn ) 25 ce qui montre que le tir composé d'indice m est obtenu par sommation des tirs individuels après modulation avec le nombre d'onde km. 5 (5) Si à partir de N tirs individuels, on choisit de migrer M dn+kM(.f,x,y,z) Eun+kM(.f,x,Y,z) f n=O\k / k L'expression (9) montre que l'image obtenue est la même que celle de la migration 10 par point de tir si les ondes descendantes et montantes dn et un, sont décorrélées pour I n-n'j M. Ceci est le cas si M est égal au nombre de tir individuels N. Nombre de tirs composés à migrer Il est cependant d'usage dans une migration par point de tir de définir une distance 15 d'illumination xu1, et de considérer qu'il est inutile de propager les ondes descendantes et montantes en dehors d'un intervalle [-x111,x111] centré sur la position de la source, les ondes étant nulles ou négligeables en dehors de cet intervalle. Dans ce cas, une distance de périodisation x,nax 2x11/ est suffisante pour garantir l'exactitude de la migration par modulations de points de tirs puisque, du fait qu'une 20 onde descendante dn et une onde montante un+M ne se recouvrent pas, leur cross-corrélation ne contribue pas à l'image. Le nombre de tirs composés à migrer pour la migration par modulations de points de tir peut encore être diminué. Définissons en effet un paramètre de lenteur maximale pmax, qui peut être calculé à partir d'une vitesse de référence v et d'un angle 25 maximum Ornai par p,nax sinOmax/v. Par défaut, v est la vitesse de surface et Ornai -90 degrés. Nous définissons un nombre de composantes à migrer par fréquence Mf tel que knax (Mf/2)Ak soit la plus petite des deux valeurs : le Nyquist spatial kNyg 1/2dx et fp,nax. Ceci veut dire que les composantes en k en dehors de l'intervalle [- fp,nax,fp,nax] ne sont pas migrées. Ceci ne dégrade pas l'image puisque cela équivaut à réduire le nombre de tirs individuels en les rééchantillonnant sur la plus grande grille satisfaisant le critère de Nyquist spatial pour la fréquence donnée. Ceci fait que la migration par modulations de points de tirs peut être aussi exacte tout en étant plus économique que la migration par point de tir. Obtention des collections indexées par p au cours d'une migration par modulation des points de tirs: L'algorithme décrit ci-dessus, qui peut être appelé migration par modulations de tirs, permet d'obtenir la même migration qu'une migration par point de tir pour un coût inférieur. Il va maintenant être décrit un procédé pour obtenir, au cours d'une telle migration, des collections indexées par une lenteur p. L'intervalle [ pmax,pmax] est divisé en Np intervalles de largeur Op. Au cours de la migration du tir composé correspondant à la fréquence et au nombre d'onde (f km), la contribution à l'image de ce tir composé est accumulée dans l'intervalle en p correspondant à p=km/f. Ceci est équivalent à reéchantillonner d'une grille régulière de pas Ak vers une grille régulière de pas fAp, ce qui peut être fait de différentes manières comme une interpolation linéaire, en sinus cardinal, ou une interpolation en bande limitée, etc... Ce reéchantillonnage peut s'effectuer à partir de pondérations précalculées 22.i(f,km), i=O,Np 1 telles que Ei %i(f,km)=l qui sont utilisées durant la migration pour actualiser les différentes collections par: Ri(x,y,z)= E Â,i(f,km)Rf,m(x,Y,z) f,m (10) où Rfm(x,y,z) est la contribution à l'image du tir composé correspondent à (f,km) : R f,m (x, y, z) = Dm ( f , x, y, z)Um (.Î, x, y, z) (11) et Ri(x,y,z) l'élément de la collection en p correspondant à l'indice i. La migration par modulation des points de tirs permet donc d'obtenir des collections 30 en p analogues à la migration par onde plane. Le procédé pour les obtenir est cependant différent. La migration par modulations de tirs utilise une composition de tirs indépendantes de la fréquence et correspondant à une modulation alors que la migration par ondes planes utilise une composition de tirs dépendant de la fréquence et correspondant à des retards. Une étape propre à la migration modulée, la distribution dépendante de la fréquence de la contribution courante à l'image parmi les collections permet à celle-ci d'obtenir les collections indexées en p. L'avantage de la migration par modulation des points de tirs est que le nombre de tir à migrer est plus petit que la migration par onde plane. Ceci est illustré par les figures 1 et 2, chaque croix correspondant à un tir composé. Dans ce cas, le coût de la migration par modulation des points de tirs est de 42% celui de la migration par onde plane. La migration par modulation des points de tirs peut utiliser une matrice de composition dépendant de la fréquence, en utilisant par exemple une distance de périodisation xmax dépendant de la fréquence xmax. Exemple de mise en oeuvre On a représenté sur la figure 5 les étapes principales (modulation, combinaison linéaire, puis calcul des collections) d'un exemple de mise en oeuvre d'une migration par points de tirs du type de celles déjà décrites. La migration par modulations de tirs a deux paramètres principaux: la distance de périodisation xmax et la lenteur maximum Pmax. N étant le nombre de tirs et Ax la distance entre tirs, la matrice de composition est: Cmn = 1 exp 2 jIcm ndx ,m E ù Mf , Mf ,n /N xmax / 2 2 _ (12) avec: Mf=min xmax 2fPmaxxmax Ax / (13) où Mf est le nombre de tirs composés à migrer pour la fréquence f, le tir composé d'indice m correspondant au nombre d'onde km=m/xmax. Pour le calcul des collection en lenteurs, après avoir choisi l'intervalle Ap entre les collections : Ri(x,Y,Z)= E Ài(.f,km)Rf,m(x,Y,Z) f,m (14) où Rfm(x,y,z) est la contribution à l'image du tir composé et 24km) les coefficients d'un filtre qui reéchantillonne une grille de pas 1/xma, vers une grille de pas fAp. Autre exemple de mise en oeuvre En variante notamment, au lieu de synthétiser les tirs modulés à la surface, ces tirs peuvent être synthétisées en profondeur. Cette procédure, décrite dans le cas général par Rietveld et al (1992) s'applique au cas de la migration par modulations de tirs. On place à une profondeur zo du milieu L des sources fictives aux abscisses xi et on leur affecte une amplitude égale à l'élément de la matrice de composition Cmi (f), m correspondant à un nombre d'onde km donné (les sources fictives peuvent aussi être placées à des profondeurs différentes). Pour chaque m, on calcule le champ montant généré par ces sources fictives en profondeur aux positions des N sources réelles (les positions des tirs individuels, habituellement à la profondeur z=0 et aux abscisses xn). On obtient des ondelettes Smn(f). On applique alors comme matrice de composition des sources en surface le complexe conjugué Cmn(f=Smn(f)*•Les propriétés de quasiunitarité des matrices de propagation font que l'onde descendante du tir composé aura la forme voulue d'une suite de tir modulés lorsqu'il aura été extrapolé à la profondeur zo. Ceci permet d'avoir des collections où la complexité de la propagation entre la surface et la profondeur zo a été enlevée. Les sources peuvent aussi être placées à des profondeurs différentes. Cas 3D Le procédé décrit peut être généralisé au traitement 3D. Les sources sont indexées de n=O,N-1 et se trouvent sur une position (xn,yn) sur une 30 grille de pas Ax en x et Ay en y. Définition des paramètres xmax et ymax, pxmax, Apx,py,,,ax,Opy. Pour chaque fréquence, calcul de MM(f) et My(f) : Mx (.Î) = min 2fpxmax , xx , M y (Î) =min 2fpymax ' YDyx / (15) Calcul des nombres d'onde kxä et ky,ny, indexés par mx,my avec mx [-Mx(f)/2,Mx(f)/2] et my=[-My(f)/2,My(f )l2] : xmax Ymax (16) Calcul des ondes descendantes et montantes des tirs composés à partir de celles des tirs individuels par : Nù1 Dmx,my(f,x,Y,z=O) E Cmx,my,ndn(f,x,Y,z=O) n=0 Nù1 Umx,my (f,x,y,z = 0) = Cmx,my,nun(f,xY,Z =0) n=0 M(f) Mx (.Î) M(f) M y (.f ) mx 2 2 ' my 2 2 avec : 1 f Cmx,my,n = exp[2j7r(kxmxxn +kymyYn ,I N (18) Pour tout mx,my, extrapolation de D,,,x,my et U,,,my pour toutes les valeurs de z sur la grille de pas 4z . Calcul de la contribution à l'image du tir composé mx,my : kx = mx , kv = my mx (17) R f,mx,my (x, y, z) = Dmx,my ( f,x,Y,z)Umx,my (f,x,Y,z) (19) Actualisation des collections en vecteur lenteur p=(px,py) indexées en ix,iy correspondant apz ixApx, py iyLpy : R1x,zy(x, y, z) = ,iy (f,kxmx,kYmy )Rf,mx,my (x,y,z) f,mx,my (20) ou XL1,iy(f,kxmx,kymy) est un filtre qui interpole une grille à 2 dimensions au pas Ak,,,Aky vers une grille au pas fApxffApy. Références Schultz, P. S., Claerbout, J.F. [1978]: Velocity estimation and downward continuation by wavefront synthesis, Geophysics, 43, 691-714. 15 Rietveld, W. E. A., Berkhout A. J., Wapenaar C.P.A. [1992]: Optimum seismic illumination of hydrocarbon reservoirs, Geophysics, 57, 1334-1345. Ober, C. C., Romero, L. A., Ghiglia, D. C [2000]: Method of migrating seismic records, US Patent 6021094, delivré le 1 er Février 2000. Lailly, P., Duquet B., Ehinger, A [2000] : Méthode pour réaliser en 3D avant 20 sommation, une migration de données sismiques, Brevet 2784195, publié le 7 Avril 2000. Duquet, B., Lailly, P. and Ehinger, A. [2001] 3D plane wave migration of streamer data. SEG 71" Annual International Meeting, Expanded Abstracts 1033-1036. Romero, L. A., Ghiglia, D. C., Ober, C. C. and Morton, S. A. [2000] Phase encoding 25 of shot records in prestack migration, Geophysics, 65, 426-436. 10 Revendications 1. Procédé de traitement sismique du type par migration, selon lequel : on combine linéairement des ondes descendantes émises par différentes sources sismiques en différents points de tirs, ainsi que des ondes montantes enregistrées par différents capteurs sismiques, on propage fictivement les ondes composées ainsi obtenues afin d'obtenir des ondes descendantes et montantes migrées pour différentes profondeurs, on détermine en différentes profondeurs au moins une caractéristique du sous-sol en fonction des ondes montantes et descendantes ainsi propagées; caractérisé en ce que les ondes composées descendantes (respectivement montantes) sont calculées par une combinaison linéaire des ondes descendantes (resp. montantes) selon laquelle lesdites ondes descendantes (resp. montantes) sont pondérées par les coefficients d'une matrice de modulation spatiale fonction de la position des points de tirs. 2. Procédé de traitement sismique selon la revendication 1, selon lequel, dans le cas d'une ligne de points de tirs, une pondération appliquée à une onde descendante ou montante est exp(2 jrckn,x,,où km est un nombre d'onde défini par mAkx, avec m nombre entier correspondant à l'indice de l'onde composée calculée et Akx = =1 /Xmax, où Xmax est une distance de périodisation. 3. Procédé selon la revendication 1, selon lequel, dans le cas d'une grille de points de tirs selon deux dimensions x et y, une pondération appliquée à une onde 30 descendante ou montante pour une position de tirs x,,,yä est :exp 2pr kx,,,xxä +kyntyytt où kmx et kmy sont deux nombres d'ondes définis par 5 k,,,=mxAkx et kmy = myAky , avec mx et my deux nombres entiers correspondant aux indices de l'onde composée calculée, et 10 Akx = 1 /Xmax, où Xmax est une distance de périodisation. Akx = 1 /ymax, où ymax est une distance de périodisation. 4. Procédé selon l'une des revendications 2 ou 3, selon lequel la ou les distances de périodisation dépendent de la fréquence. 15 5. Procédé selon la revendication 2, selon lequel pour chaque fréquence f, le nombre Mf d'ondes descendantes ou montantes composées que l'on calcule et sur lesquelles on applique le traitement de migration est Mf tel que kmax=(Mj/2)Ak est la plus petite des deux valeurs que constituent d'une part le Nyquist spatial 20 kNyq=1/2Ax et d'autre partfpmax, pmax étant un paramètre de lenteur maximale. 6. Procédé selon la revendication 3, selon lequel pour chaque fréquence f, les indices mx et my sont respectivement compris dans les intervalles [-Mx(f)/2,M;(f)/2] et [-My(f)/2,My(f)/2], où Mxf et Myf sont tels que kxmax (Mxj/2)Ak 25 est la plus petite des deux valeurs que constituent d'une part le Nyquist spatial kNyq=1/2Ax et d'autre part fpmax, pmax étant un paramètre de lenteur maximale, et kyxmas (Mxj/2)Ak est la plus petite des deux valeurs que constituent d'une part le Nyquist spatial kNyq=1/24x et d'autre part fpmax, pmax étant un paramètre de lenteur 30 maximale. 13 10 15 20 7. Procédé selon l'une des revendications précédentes, selon lequel des collections migrées sont calculées en distribuant la contribution à l'image des points de tirs composés dans des intervalles fonction de paramètres de lenteurs. 8. Procédé de traitement sismique selon lequel, pour déterminer des ondes composées en surface, on affecte à une pluralité de sources fictives, considérées à une profondeur donnée par rapport à la surface ou à différentes profondeurs, une amplitude égale à l'élément d'une matrice de composition (Cml (f)), correspondant à une matrice de modulation correspondant elle-même à une pondération telle que définie dans les revendications 2 à 6, pour chaque m, (ou chaque couple (mx, my)) on calcule le champ montant généré par ces sources fictives en profondeur aux positions des sources réelles et on en déduit différentes ondelettes de surface (Smn({I) on met en oeuvre le procédé selon l'une des revendications précédentes en appliquant comme matrice de composition, à la place de la matrice de modulation, le complexe conjugué des ondelettes de surface ainsi déterminées (Cmn(f)=smn(f) )•