FR 3024243 A1 20160129 FR 1457140 A 20140724 1 La présente invention concerne le domaine de l'industrie pétrolière, plus particulièrement le domaine de l'exploration et l'exploitation des réservoirs pétroliers ou de sites de stockage géologique de gaz. Une technique largement utilisée dans l'industrie pétrolière pour la recherche et l'évaluation de réserves souterraines est la prospection sismique. La prospection sismique consiste en général en trois étapes : l'acquisition de données sismiques, le traitement de ces données, et enfin l'interprétation des données traitées, alors appelées image sismique. L'étape d'acquisition des données sismiques met généralement en oeuvre le principe de la sismique réflexion. La sismique réflexion consiste à émettre une onde ou plusieurs ondes (par explosif ou vibration en sismique terrestre, par canon à air ou à eau en sismique marine) et à enregistrer les signaux représentant les variations d'amplitude des ondes s'étant propagées dans l'eau ou le sous-sol, et s'étant au moins partiellement réfléchies au niveau d'au moins une limite (appelée aussi interface) de couche géologique caractérisée par un contraste d'impédance sismique. L'enregistrement des ondes s'étant ainsi réfléchies est effectué par des capteurs d'accélération (sismomètres), de vibration (géophones) ou de pression (hydrophones), ou par une combinaison de capteurs élémentaires des types précédents (par exemple des capteurs multi-composantes ou des OBC ("Ocean Bottom Cable")) localisés à des positions géographiques données. Le signal enregistré, typiquement pendant quelques secondes, par un capteur donné est appelé trace sismique. Les données sismiques correspondent à une collection d'une ou plusieurs traces sismiques, issues de capteurs localisés à différentes positions spatiales, formant un volume en deux ou trois dimensions (une de temps, une ou deux d'espace), voire en quatre dimensions si l'on inclut des acquisitions sismiques répétitives (acquises dans une même région à différentes périodes temporelles). La distance entre une source et un capteur sismiques est appelée déport (ou "offset" en anglais). Les données sismiques enregistrées au cours d'une expérience de sismique réflexion sont dites multi-déports ou encore avant sommation ("prestack" en anglais), c'est-à-dire que le signal émis par une source donnée est enregistré par plusieurs capteurs situés à différents déports. Les données sismiques enregistrées, dites brutes, sont bien souvent inexploitables. Suivant la qualité et les caractéristiques des données enregistrées, sont appliquées différentes étapes de traitement sismique, telles que les corrections d'amplitudes, la 3024243 2 déconvolution, les corrections statiques, le filtrage du bruit (aléatoire ou cohérent ), la correction NMO (« Normal Move Out » en anglais, « Courbure Normale » en français), le stack (« sommation » en français, conduisant à une section à déport nul ou après sommation des données sismiques avant sommation), ou encore la migration. Ces étapes de traitement, 5 nécessitant des calculs souvent très complexes et très longs, sont réalisées sur ordinateur. Les données sismiques résultantes sont alors appelées image sismique. Ces images sismiques sont le plus souvent représentées sur un ordinateur, par un maillage ou grille, chaque maille correspondant à une position latérale et verticale (la direction verticale correspondant au temps ou à la profondeur selon que le traitement a abouti à une image 10 temps ou à une image profondeur) au sein de la formation étudiée, et étant caractérisée par une amplitude sismique. Si le traitement sismique appliqué aux données sismiques enregistrées est optimal, l'amplitude sismique en une maille donnée d'une image sismique doit refléter l'amplitude de l'onde sismique ayant subi une unique réflexion (dite réflexion primaire) à la position de la maille dans la formation étudiée. 15 Les images sismiques sur lesquelles sont conduites les interprétations par un géologue ou bien un géophysicien doivent être suffisamment précises pour rendre compte des caractéristiques de la formation étudiée, notamment en terme de géométrie des couches géologiques et des failles, mais aussi en terme d'amplitudes sismiques, qui informent sur les 20 propriétés pétro-physiques de la formation étudiée. En effet, de ces informations sont construites des représentations de la formation étudiée, appelées modèles géologiques, qui permettent de déterminer de nombreux paramètres techniques relatifs à la recherche, l'étude ou l'exploitation d'un réservoir, d'hydrocarbures par exemple. 25 Une étape particulièrement délicate du traitement sismique consiste à filtrer le bruit engendré par des ondes parasites cohérentes, appelées réflexions multiples, qui subissent un ou plusieurs rebonds dans la couche d'eau ou entre au moins deux interfaces géologiques. La Figure 1 illustre le trajet d'une réflexion primaire P1 arrivant en un point A et le trajet de deux réflexions multiples M1 et M2 arrivant en un point B. Les réflexions multiples 30 ont connu une part de trajet en commun avec la réflexion primaire P1, mais ont, en plus, subi un rebond entre deux interfaces limitant en haut et en bas une couche géologique. On parle de multiples de première ordre lorsque l'onde sismique a subi un seul rebond, de second ordre si elle a subi deux rebonds, etc. Une réflexion multiple arrivant en B est enregistrée avec un certain retard par rapport la réflexion primaire arrivant en A, correspondant au 3024243 3 temps de propagation pour effectuer le rebond entre les deux interfaces. Suivant les vitesses de propagation dans les couches, il est possible qu'une réflexion multiple enregistrée en B interfère avec des réflexions primaires associées avec des réflecteurs sismiques plus profonds, comme par exemple avec la réflexion primaire P2 présentée en Figure 1. Par ces 5 interférences parasites, les réflexions multiples peuvent masquer ou encore fausser les informations utiles contenues dans les réflexions primaires, en modifiant la géométrie des réflexions primaires ou encore leur amplitude sismique. Avant toute interprétation d'une image sismique, il est donc nécessaire que les 10 données sismiques enregistrées subissent un traitement approprié permettant de filtrer, c'est-à-dire d'éliminer ou au moins d'atténuer, les réflexions multiples parasites. État de la technique 15 Les documents suivants seront cités au cours de la description : Nowak, E. J. and lmhof, M. G., Amplitude preservation of Radon-based reflexion multiple-removal filters, Geophysics, 2006, 71, V123-V126. 20 Chaux, C., Duval, L., Pesquet, J-C, Image Analysis Using a Dual-Tree M-Band Wavelet Transform, IEEE Transactions on Image Processing, August 2006, Volume 15, Issue 8, p. 2397-2412. Pica, A.; Poulain, G.; David, B.; Magesan, M.; Baldock, S.; Weisser, T.; Hugonnet, P. & 25 Herrmann, P. 3D surface-related multiple modeling, The leading Edge, 2005, 24, 292-296. Ventosa, S., Le Roy, S., Huard, I., Pica, A., Rabeson, H., Ricarte, P., Duval, L., Adaptive multiple subtraction with wavelet-based complex unary Wiener filters, Geophysics, 2012. 30 Il existe principalement deux familles de méthodes classiquement utilisées pour filtrer les réflexions multiples contenues dans des données sismiques : 3024243 4 1) Une famille basée sur l'usage d'une méthode de décomposition (telle que la transformée de Fourier, la transformée de Radon ou encore la transformée en ondelettes), basée sur l'hypothèse que les données sismiques correspondent à une sommation de composantes, chaque composante ayant des caractéristiques propres (par exemple une 5 gamme de fréquences particulières couplée à une gamme d'orientations dans l'espace particulière). Une composante est définie par une fonction mathématique, dépendante du type de décomposition choisi, et par un coefficient de décomposition. Dans le cas de l'application d'une transformée en ondelettes directionnelles visant à atténuer des réflexions multiples, les valeurs des coefficients correspondant aux réflexions primaires sont 10 conservées, les valeurs des autres coefficients, correspondants aux réflexions multiples que l'on cherche à éliminer, étant annulées ou mis à une valeur très petite par rapport à celles des coefficients correspondant aux réflexions primaires. Une transformation inverse est ensuite mise en oeuvre afin de recomposer les données sismiques, qui sont alors au moins partiellement atténuées des réflexions multiples. Une telle technique est décrite par exemple 15 dans Nowak and lmhof (2006). 2) Une famille basée sur l'usage d'un filtre adaptatif des réflexions multiples. Cette méthode consiste à prédire un ou plusieurs modèles de réflexions multiples, puis à le ou les soustraire aux données sismiques. Plus précisément, à partir d'un ou plusieurs 20 modèles de réflexions multiples, cette méthode consiste à estimer un ou plusieurs filtres adaptatifs, ayant un nombre de coefficients limité. Différentes méthodes d'obtention des coefficients de filtre sont connues. Une telle technique est par exemple décrite dans le document FR2994746. Un mode particulier d'application du filtrage adaptatif consiste à minimiser l'écart quadratique entre les données sismiques et le modèle de réflexions 25 multiples, en faisant une hypothèse d'orthogonalité entre les réflexions primaires et les réflexions multiples. Il n'existe en général pas de solution unique et satisfaisante quelles que soient les données sismiques considérées. Des méthodes alternatives de filtrage des réflexions 30 multiples peuvent être obtenues par combinaison des deux familles de méthodes citées ci- dessus. Une telle technique est par exemple décrite dans le document Ventosa et al. (2012). Contrairement aux autres méthodes connues, cette technique permet de compenser les défauts des méthodes de filtrage adaptatif classiques, qui ne préservent pas au mieux le 3024243 5 contenu fréquentiel, notamment les basses et les hautes fréquences des données sismiques. Ainsi cette méthode consiste à filtrer les données et les modèles dans une bande de fréquence, puis à calculer leur coefficient d'inter-corrélation dans cette bande de fréquence, et enfin à soustraire aux données originales le produit de chaque modèle, filtré 5 dans la bande de fréquence considérée, par le coefficient de corrélation obtenu. Cependant, la méthode décrite dans Ventosa et al. (2012) ne permet pas d'assurer une cohérence latérale optimale du filtrage de multiples en présence de bruits importants, notamment en cas de pendages importants des réflexions primaires ou des réflexions multiples. 10 De façon générale, les méthodes de l'art antérieur ne permettent pas un filtrage optimal des données sismiques, c'est-à-dire garantissant une élimination complète des réflexions multiples, tout en préservant les caractéristiques des réflexions primaires, telles que les amplitudes et les fréquences. En particulier, les méthodes d'acquisition sismique actuelles, telle que la technologie BroadSeisTM développée par CGG (France), permettant 15 d'acquérir des données sismiques ayant un très large contenu fréquentiel, il apparait essentiel que le filtrage des réflexions multiples garantisse au mieux la préservation du contenu fréquentiel des données sismiques enregistrées. La présente invention est une méthode alternative de filtrage des réflexions multiples 20 présentes dans des données sismiques, combinant une méthode de décomposition et une méthode de filtrage adaptatif, suivies d'une recombinaison pondérée des composantes filtrées. En particulier, la présente invention a pour objectif un meilleur filtrage des réflexions multiples, de mieux conserver l'amplitude des réflexions primaires, notamment quand ces dernières sont d'amplitude faible par rapport aux réflexions multiples et aux bruits aléatoires. 25 Le procédé selon l'invention Ainsi, la présente invention concerne un procédé pour construire une image sismique filtrée de réflexions multiples, à partir d'un enregistrement de données sismiques comprenant 30 des réflexions primaires et des réflexions multiples, à partir d'au moins un modèle desdites réflexions multiples. Le procédé comporte au moins les étapes suivantes : 3024243 6 a) on applique une méthode de décomposition pour décomposer, selon N directions de décomposition, lesdites données sismiques en un ensemble de N composantes desdites données sismiques ; b) pour au moins un desdits modèles desdites réflexions multiples, on applique 5 ladite méthode de décomposition pour décomposer, selon lesdites N directions de décomposition, ledit modèle de réflexions multiples en un ensemble de N composantes dudit modèle ; c) pour chaque direction de décomposition et pour au moins un desdits modèles desdites réflexions multiples, on calcule une concentration relative entre ladite 10 composante desdites données sismiques dans ladite direction et ladite composante dudit modèle dans ladite direction ; d) à partir d'au moins un desdits modèles desdites réflexions multiples, on applique une méthode de filtrage adaptatif desdites réflexions multiples sur chacune desdites N composantes desdites données sismiques, et on obtient un ensemble 15 de N composantes desdites données sismiques filtrées ; e) on calcule une recombinaison pondérée desdites composantes desdites données sismiques filtrées, à partir d'une pondération calculée pour chaque direction de décomposition en fonction desdites concentrations relatives calculées dans ladite direction de décomposition. 20 Selon un mode de réalisation de la présente invention, ladite méthode de décomposition peut être une transformée en ondelettes. Selon un mode particulier de réalisation de la présente invention, ladite méthode de 25 décomposition peut être une transformée en ondelettes M-bandes en arbre dual. Selon un mode de réalisation de la présente invention, ledit filtrage adaptatif peut être un filtre de Wiener unaire en trame d'ondelettes complexes. 3024243 7 Selon un mode préféré de mise en oeuvre de la présente invention, ladite concentration relative entre un signal S et un signal S' bi-dimensionnels peut être calculée selon la formule suivante : s21 51,52 SNIS'l 2 2 , S2 )1 CR(S(Si , S2), S'(Sp S2 )) S.j7S2 SI,S2 s2 )12 SI,S2 5 où si est le nombre d'échantillons dans une direction de l'espace et s2 est le nombre d'échantillons dans l'autre direction de l'espace. Selon un mode de réalisation de la présente invention, ladite recombinaison pondérée peut être calculée de la façon suivante : 10 SR = aIen.1Y (0 , , SDF , ,0) n=1,N où SDF, est ladite composante dans ladite direction n (n=1,N) desdites données sismiques filtrées, a est une constante, en est ladite pondération pour ladite direction de décomposition n, et D' est une inverse de ladite méthode de décomposition. 15 Selon un autre mode de réalisation de la présente invention, ladite recombinaison pondérée peut être calculée de la façon suivante : SR = crI D' (0,.., n .S D n , ,0) n , N où SDF, est ladite composante dans ladite direction n (n=1,N) desdites données sismiques 20 filtrées, a est une constante, en est ladite pondération pour ladite direction de décomposition n, et D' est une inverse de ladite méthode de décomposition. 3024243 8 Selon un mode de réalisation de la présente invention dans lequel on connait de manière approchée ladite inverse de ladite méthode de décomposition, on peut ajouter auxdites composantes desdites données sismiques, un résidu correspondant à la différence entre lesdites données sismiques et le résultat de ladite inverse approchée. 5 Selon un mode particulier de mise en oeuvre de la présente invention, ladite pondération pour ladite direction de décomposition n (n=1, N) peut être calculée selon la formule suivante : len en = i=1,I I 10 où en est une pondération pour ladite direction n et ledit modèle de réflexions multiples i (i=1,1). Selon un mode particulier de mise en oeuvre de la présente invention, ladite pondération pour ladite direction de décomposition n (n=1, N) peut être calculée selon la 15 formule suivante : e ==rilax(e) n n où en est une pondération pour ladite direction n et ledit modèle de réflexions multiples i (i=1,1). 20 Préférentiellement, ladite pondération pour ladite direction de décomposition n et ledit modèle de réflexions multiples i (avec i./,/ et n=1,N) peut être calculée selon la formule suivante : 1 W(x) , avec x = min(Cen, ) , CR où W est une fonction croissante, CRn est ladite concentration relative entre ladite 25 composante desdites données sismiques dans ladite direction n et ladite composante dudit modèle i dans ladite direction n. 3024243 9 Avantageusement, ladite fonction croissante W peut être définie telle que pour un x donné, W(x) = x' avec p . Préférentiellement, on peut choisir p égal à 4. 5 On peut définir un procédé d'exploitation d'une formation souterraine en réalisant les étapes suivantes : on construit une image sismique filtrée des réflexions multiples au moyen du procédé tel que décrit selon l'une des revendications 1 à 13 ; on construit un modèle géologique représentatif de la formation étudiée à 10 partir d'au moins l'image sismique ainsi déterminée ; on détermine un schéma d'exploitation optimal du réservoir à partir du modèle géologique ainsi déterminé ; on exploite ledit réservoir en mettant en oeuvre ledit schéma d'exploitation optimal. 15 En outre, l'invention concerne un produit programme d'ordinateur téléchargeable depuis un réseau de communication et/ou enregistré sur un support lisible par ordinateur et/ou exécutable par un processeur, comprenant des instructions de code de programme pour la mise en oeuvre du procédé selon la description ci-dessus, lorsque ledit programme 20 est exécuté sur un ordinateur. Présentation succincte des Figures D'autres caractéristiques et avantages du procédé selon l'invention, apparaîtront à la 25 lecture de la description ci-après d'exemples non limitatifs de réalisations, en se référant aux figures annexées et décrites ci-après. La Figure 1 présente un dispositif d'acquisition de données sismiques ainsi que des exemples de trajectoire de réflexions primaires et de réflexions multiples générées par ce 30 dispositif. 3024243 10 Les Figures 2A et 2B montrent un exemple de données sismiques avant sommation et le modèle de réflexions multiples correspondant. 5 La Figure 3 présente un ensemble d'ondelettes différentes, chaque ondelette étant caractérisée par une gamme de fréquences et une gamme d'orientations dans l'espace particulières. Les Figures 4A, 4B, 4C, 4D illustrent les coefficients de décomposition des données 10 présentées en Figure 2A, utilisant 4 ondelettes orientées différentes, et permettant d'obtenir des directions préférentielles sur quatre composantes SDi,SD,, SD3, et SD4. Les Figures 4E, 4F, 4G, 4H illustrent les coefficients de décomposition du modèle de multiples présenté en Figure 2B, utilisant 4 ondelettes orientées différentes, et permettant 15 d'obtenir des directions préférentielles sur quatre composantes MD,,MD2, MD3, et MD4. La Figure 5A montre les données sismiques déjà présentées en Figure 2A, avant filtrage du bruit aléatoire. La Figure 5B montre le résultat obtenu après application de la méthode décrite dans Ventosa et al. (2012). La Figure 5C montre le résultat obtenu après 20 application du procédé selon l'invention. Description détaillée du procédé Les définitions suivantes sont utilisées au cours de la description de l'invention : 25 - modèle de multiples : il s'agit d'un modèle approché des réflexions multiples contenues dans des données sismiques. Il existe plusieurs procédés d'obtention de ces modèles de réflexions multiples. Un procédé consiste à obtenir une version approchée de la réflexion primaire (par exemple par un filtrage), puis à la convoluer 30 avec elle-même ou encore avec la trace sismique initiale (voir par exemple Pica et al. 3024243 11 (2005)). Un modèle de multiples peut être également obtenu par résolution de l'équation des ondes dans le milieu considéré. En général, les modèles de multiples sont des approximations satisfaisantes des réflexions multiples. Ils peuvent toutefois être décalés sur l'axe vertical (axe des temps ou des profondeurs), avoir des 5 amplitudes et/ou un spectre en fréquence différents par rapport aux vraies réflexions multiples. A cause de ces imprécisions, le traitement d'atténuation des réflexions multiples a recours à une adaptation des modèles de multiples aux données sismiques, par filtrage adaptatif, également nommé recalage. A noter qu'un modèle de multiples peut avoir une pertinence limitée, en étant par exemple représentatif d'une 10 réflexion multiple pour une gamme de déports limitée, ou encore pour un ordre de multiples donné. On peut alors avoir recours à plusieurs modèles de multiples pour simuler, complètement, une réflexion multiple enregistrée dans des données sismiques. 15 - méthode de décomposition en ondelettes: il s'agit d'une méthode permettant de décomposer des données sismiques en différentes composantes représentées par des ondelettes. Une formule pour obtenir de telles ondelettes est donnée ci-après : ---- oxi) I un 1 rH x où les coordonnées x et y correspondent à la position de l'échantillon considéré dans 20 les données sismiques, le couple (u0,v0) définit les coordonnées d'un vecteur dans une direction donnée, la norme de ce vecteur définissant la fréquence, et les termes as et cf), correspondent aux largeurs d'enveloppe dans les directions x et y. Chaque choix de paramètres fournit une ondelette directionnelle particulière. La Figure 3 présente différentes ondelettes directionnelles, chacune étant caractérisée par une enveloppe, 25 une fréquence principale et une orientation particulière dans l'espace. - énergie d'un signal : dans le cas d'un signal bi-dimensionnel S caractérisé par si échantillons dans une direction et s2 échantillons dans l'autre direction, l'énergie d'un signal est définie par la formule suivante : 30 IIS(si,s2)2. si is2 3024243 12 L'objet de la présente invention est un procédé pour construire une image sismique filtrée des réflexions multiples présentes dans des enregistrements sismiques, à partir d'un ou plusieurs modèles de multiples. Cette invention utilise une décomposition à la fois des données sismiques et du modèle de multiples selon des directions de décomposition 5 préférentielles, et une recombinaison pondérée de chacune des composantes sismiques filtrées selon les composantes du modèle de multiples correspondantes. Les données sismiques peuvent avoir été enregistrées par des dispositifs d'acquisition bi-dimensionnels ou tri-dimensionnels et être organisées selon tout type de collection (i.e. en collection tir commun, récepteur commun, etc). Les données sismiques peuvent correspondre 10 indifféremment à des collections avant sommation ou après sommation, avant migration ou après migration. La présente invention peut être appliquée à tout stade du traitement sismique, et préférentiellement avant l'étape d'interprétation de l'image sismique résultant de l'application de l'ensemble du traitement sismique. L'invention nécessite de disposer d'au moins un modèle de réflexions multiples. 15 La présente invention comporte au moins les étapes suivantes : a) on applique une méthode de décomposition pour décomposer, selon N directions de décomposition, lesdites données sismiques en un ensemble de N 20 composantes desdites données sismiques ; b) pour au moins un desdits modèles desdites réflexions multiples, on applique ladite méthode de décomposition pour décomposer, selon lesdites N directions de décomposition, ledit modèle de réflexions multiples en un ensemble de N composantes dudit modèle ; 25 c) pour chaque direction de décomposition et pour au moins un desdits modèles desdites réflexions multiples, on calcule une concentration relative entre ladite composante desdites données sismiques dans ladite direction et ladite composante dudit modèle dans ladite direction ; d) à partir d'au moins un desdits modèles desdites réflexions multiples, on applique 30 une méthode de filtrage adaptatif desdites réflexions multiples sur chacune desdites N composantes desdites données sismiques, et on obtient un ensemble de N composantes desdites données sismiques filtrées ; 3024243 13 e) on calcule une recombinaison pondérée desdites composantes desdites données sismiques filtrées, à partir d'une pondération calculée pour chaque direction de décomposition en fonction desdites concentrations relatives calculées dans ladite direction de décomposition. 5 Les principales étapes de la présente invention sont décrites ci-après dans le cas de données sismiques S après sommation à deux dimensions, mais le procédé peut tout aussi bien être appliqué à des données sismiques tridimensionnelles. Les principales étapes de la présente invention sont déclinées ci-après dans le cas d'un seul modèle de multiples M. Le 10 cas de plusieurs modèles de multiples est décliné dans une variante. a) Décomposition des données sismiques en N composantes Il s'agit dans cette étape de décomposer les données sismiques S en N composantes selon une méthode de décomposition D . Le but de cette étape est de mieux distinguer 15 certains traits caractéristiques des données sismiques. On obtient ainsi une série de composantes SI), avec n=1,N associées à N directions de décomposition différentes. Selon un mode de réalisation de la présente invention, on choisit une méthode de décomposition D permettant de décomposer des données sismiques selon différentes 20 bandes de fréquence et différentes gammes d'orientation dans l'espace. Dans ce cas, une direction de décomposition est un couple formé par une bande de fréquence et une gamme d'orientation dans l'espace. Selon un mode de réalisation de la présente invention, on utilise une méthode de 25 décomposition en ondelettes M-bandes en arbre dual, telle que décrit dans Chaux et al. (2006). Il s'agit d'ondelettes directionnelles particulières, obtenues en choisissant un ensemble de filtres, dit banc de filtres primai, chaque filtre étant défini pour une bande de fréquence propre, et l'ensemble des filtres permettant de couvrir une gamme de fréquence prédéfinie. On calcule ensuite un autre ensemble de filtres déduit des précédents, dit banc 30 de filtres dual. Cet ensemble est obtenu en décalant chaque filtre d'un demi-coefficient par des techniques d'interpolation classiques ou par calcul dans le domaine de Fourier. Ces 3024243 14 filtres sont ensuite appliqués séparément sur les directions horizontales et verticales, et combinés pour fournir des directions diagonales. En pratique, l'expert définit P bandes de fréquences, et l'on construit le banc de filtres dual comme décrit précédemment de façon à obtenir les 4P décompositions selon les différentes gammes d'orientations dans l'espace et 5 bandes de fréquence. Selon un mode préféré de réalisation de la présente, P est compris entre quatre à huit. A noter que pour ce type d'ondelettes, les enveloppes sont définies par le choix des gammes de fréquence et des orientations dans l'espace. Selon un mode de réalisation de l'invention, on utilise une méthode de décomposition 10 pour laquelle on connait exactement au moins une inverse D' , c'est-à-dire telle que D' = D . Une inverse de la méthode de décomposition, dite aussi méthode de recombinaison, permet de recombiner les différentes composantes SD', avec n=1,1V, de sorte à produire les données sismiques initiales S. 15 Selon un mode de réalisation de l'invention, on utilise une méthode de décomposition pour laquelle on connait une inverse D' de manière approchée, dite aussi pseudo-inverse D' ou encore méthode de recombinaison approchée. Dans ce cas, de sorte à produire les données sismiques initiales S, on ajoute au préalable aux différentes composantes SD', avec n=1,N, le résidu correspondant à la différence entre les données sismiques enregistrées S et 20 le résultat de la pseudo-inverse choisie D'(SDn) . Le résultat de la méthode de décomposition décrite dans Chaux et al (2006), appliquée aux données sismiques présentées en Figure 2A, est donné en Figures 4A à 4D. Ainsi, les données sismiques S ont été décomposées en utilisant 4 ondelettes directionnelles 25 différentes, permettant d'obtenir quatre composantes SDL SD2 SD3, et SD4 orientées selon 4 directions dans l'espace différentes, plus précisément, une composante isotrope à basse fréquence (Figure 4A), et trois composantes directionnelles, couvrant approximativement les plages angulaires suivantes : -23° à +23° (Figure 4B), 22° à 68° (Figure 4C), 67° à 113° (Figure 4D). 30 3024243 15 b) Décomposition du modèle de multiples en N composantes Il s'agit dans cette étape de décomposer le modèle de multiples M en N composantes, avec la même méthode de décomposition D et selon les mêmes N directions que celles 5 définies à l'étape a) de décomposition des données sismiques. On obtient ainsi une série de composantes MD' avec n=1,N. Le résultat de la méthode de décomposition décrite dans Chaux et al (2006), 10 appliquée, selon les même directions que celles utilisées pour produire les Figure 4A à 4D, au modèle de réflexions multiples présenté en Figure 2B, est donné en Figures 4E à 4H. Ainsi, le modèle de réflexions multiples M a été décomposé en quatre composantes MD/, MD2 MD3, et MD4, orientées selon les mêmes directions que celles définies pour procéder à la décomposition des données sismiques. 15 c) Calcul de la concentration relative des N composantes Au cours de cette étape, il s'agit de calculer, pour chaque direction de décomposition n (avec n=1,N) définie lors des précédentes étapes, la concentration relative entre la composante SD' issue des données sismiques dans la direction considérée et la 20 composante MD' issue du modèle de multiples dans la même direction. La concentration relative entre deux signaux correspond au rapport de la concentration de chacun des signaux. 25 Selon un mode de réalisation de l'invention, la concentration C d'un signal S bi- dimensionnel défini par si échantillons dans une direction et s2 échantillons dans l'autre direction est définie de la façon suivante : s21 SI,S2 \IIIS(Si,S2)2 5.'172 -1 La concentration ainsi définie est donc comprise entre 0 et 1. Plus précisément, un signal 30 très concentré, par exemple tel qu'une seule valeur de ce signal soit non nulle, a une c(s(si, s2)) 3024243 16 concentration C égale à 1. Par contre, un signal très peu concentré, par exemple pour lequel toutes les valeurs du signal sont égales entre elles, aura une concentration C égale à 0. La concentration permet notamment de qualifier la qualité d'un traitement visant à concentrer l'énergie d'un signal donné de manière optimale. Par exemple, un signal sinusoïdal est 5 converti par transformation de Fourier en un pic fréquentiel bien localisé. Son indice de concentration est alors maximal. Selon un mode de réalisation de l'invention, la concentration relative CR entre un signal S et un signal S' tous deux bi-dimensionnels et caractérisés par si échantillons dans une 10 direction et s2 échantillons dans l'autre direction s'écrit alors : IIS(sps2)1 5..j72 J/1S(S S2 )2 p 51,52 CR(S(Si S2), S'(S1, S2 )) 1S" (Si 5115'1 2 11151,52 La concentration relative CR, pour une direction n donnée (avec n=1,N) entre une composante SD' issue des données sismiques dans cette direction n et la composante MD' 15 issue du modèle de multiples dans la même direction indique la qualité d'une méthode de décomposition à concentrer de façon équivalente les données sismiques et le modèle de multiples pour cette direction n. Ainsi, si la concentration relative CR, pour une direction n donnée est supérieure à 1, le signal sismique est plus concentré que le modèle de multiples pour la direction considérée. Et inversement, si la concentration relative est inférieure à 1, 20 alors le modèle de multiples est plus concentré que le signal sismique pour cette même direction. d) Filtrage adaptatif des N composantes sismiques Cette étape fait appel à une méthode F de filtrage adaptatif des réflexions primaires et 25 des réflexions multiples, à partir d'un modèle de multiples. Plus précisément, au cours de cette étape, une méthode de filtrage adaptatif F est appliquée à chacune des composantes SD, issues des données sismiques et calculées à l'étape a). On obtient ainsi N composantes sismiques SDF, filtrées des réflexions multiples de la façon suivante : 51,52 s2)12 51,52 3024243 17 SDFT, = F(SD,,,MD') avec n=1,N. Ainsi, le filtrage adaptatif, est appliqué non pas aux données sismiques elles-mêmes mais sur les composantes sismiques issues de la décomposition des données sismiques 5 selon des fréquences et des orientations dans l'espace privilégiées. De cette façon, les composantes sismiques étant plus simples et/ou plus homogènes que les données sismiques originelles, les paramètres du filtrage adaptatif sont plus simples à estimer, et par voie de conséquence, le filtrage adaptatif est plus efficace. Plus précisément, le fait d'appliquer un filtre adaptatif sur chaque composante issue des données sismiques a pour 10 but de mieux éliminer les réflexions multiples et de mieux préserver le contenu fréquentiel, composante par composante. Selon un mode de réalisation de l'invention, la méthode de filtrage adaptatif F est choisie pour sa qualité de préservation du contenu fréquentiel, notamment des basses et hautes fréquences des réflexions primaires. 15 Selon un mode préféré de réalisation de l'invention, on utilisera une méthode de filtrage adaptatif de réflexions multiples par filtres de Wiener unaires en trame d'ondelettes complexes mono-dimensionnelles telle que décrit dans Ventosa et al. (2012). Cette méthode permet en effet de préserver au mieux le contenu fréquentiel, notamment les basses et les 20 hautes fréquences des données sismiques. e) Recombinaison pondérée des N composantes sismiques filtrées Cette étape a pour objectif de recombiner les N composantes sismiques SDF' filtrées des réflexions multiples obtenues lors de l'étape d), la recombinaison étant pondérée en 25 tenant compte des concentrations relatives CR' (avec n=1,N) calculées à l'étape c). Ainsi, l'objectif de cette étape est de recombiner les différentes composantes sismiques filtrées, en tenant compte de l'efficacité de la méthode de décomposition appliquée aux données sismiques et au modèle de multiples. 30 Selon un mode de réalisation de la présente invention, pour réaliser la recombinaison des différentes composantes sismiques filtrées, on a recours à une inverse D' de la méthode de décomposition D telle que définie à l'étape a) 3024243 18 Selon un mode de réalisation de la présente invention, pour réaliser la recombinaison des différentes composantes sismiques filtrées, on a recours à une pseudo-inverse D', ou inverse approchée, de la méthode D telle que définie à l'étape a). 5 Selon un mode de réalisation de l'invention, on privilégie les composantes concentrées de manières équivalentes, en utilisant une pondération maximale quand la concentration relative dans la direction considérée est proche de 1, et minimale dans le cas contraire. 10 Selon un mode privilégié de réalisation de l'invention, on considère une fonction W croissante et l'on définit une pondération pour une direction donnée n par la valeur de cette fonction en un point x défini par le minimum entre la concentration relative calculée pour 1 cette direction et son inverse, soit en -W (x) avec x = min(CR , ).Dans ce cas, le point n CR n d'évaluation x de la fonction W est nécessairement compris entre 0 et 1. De plus, W étant 15 choisie croissante, on donne ainsi un poids plus important aux valeurs de x proches de 1, c'est-à-dire aux composantes équitablement concentrées. Selon un mode de réalisation de l'invention, on choisit une fonction W telle que pour une valeur x donnée, W(x) = x" avec p O. 20 Selon un mode préféré de réalisation de l'invention, on choisit une fonction W telle que pour une valeur x donnée, W(x) - xP avec p = 4. Selon un mode de réalisation de l'invention, on calcule une estimation SR des données 25 sismiques filtrées des réflexions multiples, résultant des meilleures combinaisons des N composantes sismiques SDFn filtrées par la méthode F et respectant la concentration relative des différentes composantes issues de la méthode de décomposition D , selon la formule suivante : SR = celen.D' (0,..,SDFn,...,0) n=1,N 3024243 19 où a est une constante. Selon un mode de réalisation de l'invention, a est choisie de telle manière que si l'on affecte des valeurs nulles au modèle de multiples M, SR est caractérisé par la même énergie que les données sismiques S en entrée. Selon un mode particulier de réalisation de l'invention, a est constant et égal à 1. De cette façon, on recombine, selon une 5 inverse de la méthode D , les composantes SDF' filtrées des réflexions multiples, en les annulant toutes sauf une, et en affectant la pondération calculée à la recombinaison de la composante en question. Selon un autre mode de réalisation, on calcule une estimation SR des données 10 sismiques filtrées des réflexions multiples de la façon suivante : SR = al D' (0,..,en.SDFn,...,0) n=1,N où a est une constante. Selon un mode de réalisation de l'invention, a est choisie de telle manière que si l'on affecte des valeurs nulles au modèle de multiples M, SR est caractérisé par la même énergie que les données sismiques S en entrée. Selon un mode particulier de 15 réalisation de l'invention, a est constant et égal à 1. De cette façon, on obtient une image sismique filtrée des réflexions multiples, grâce à un filtrage adaptatif appliqué sur des données sismiques décomposées selon des directions préférentielles, suivi d'une recombinaison des composantes filtrées tenant compte de 20 l'efficacité de la méthode de décomposition. Variantes Cas de plusieurs modèles de multiples 25 Dans le cas où l'on dispose de plusieurs modèles de multiples M' avec i./,/ on répète pour chaque modèle : l'étape b) de décomposition du modèle de multiples décrite ci-dessous : on obtient ainsi une série de décompositions MD;, avec i./,/ et n=1,N ; l'étape c) de calcul de la concentration relative : on calcule ainsi les 30 concentrations relatives CR;, pour chaque direction n et pour chaque modèle de multiples i, avec i./,/ et n=1,N. 3024243 20 Puis on applique, pour composante n (avec n=1,N), la méthode de filtrage F en tenant compte des différents modèles de multiples. Plus précisément on calcule les N composantes sismiques filtrées de la façon suivante : SDFn = (SIDn,MIDni MIDni ) avec n=1,N. 5 Puis, on répète, pour chaque modèle de multiples, le calcul de pondération décrit dans l'étape e) de recombinaison pondérée des composantes sismiques filtrées. On obtient ainsi une série de pondérations en' pour toutes les valeurs de i et de n, avec i./,/ et n=1,N. Puis, pour chaque direction n (avec n=1,N), on calcule une pondération composite en à 10 partir de l'ensemble des pondérations en (avec i=/,/) calculées pour chaque modèle de multiples. Selon un mode de réalisation de la présente invention, on détermine une pondération composite en en calculant la moyenne des pondérations individuelles en' soit en . 15 Selon un mode de réalisation préféré de la présente invention, afin de privilégier, pour une direction n donnée (avec n=1,N), la concentration relative entre les données sismiques S et au moins l'un des modèles de multiples A/11, on attribue à la pondération composite en la valeur maximale parmi la série de pondérations en (avec i=/,/) calculées pour chaque modèle de multiples. 20 Puis, après le calcul de la pondération composite en, on calcule une estimation SR des données sismiques filtrées des réflexions multiples par recombinaison pondérée des composantes sismiques filtrées SDF' comme décrit précédemment dans l'étape e). 25 Exploitation d'une formation souterraine En outre, l'invention concerne un procédé d'exploitation d'une formation souterraine, dans lequel on réalise les étapes suivantes : - on construit une image sismique filtrée des réflexions multiples au moyen du procédé tel que décrit précédemment ; 3024243 21 on construit un modèle géologique représentatif de la formation étudiée à partir d'au moins l'image sismique ainsi déterminée ; on détermine un schéma d'exploitation optimal du réservoir à partir du modèle géologique ainsi déterminé ; 5 on exploite ledit réservoir en mettant en oeuvre ledit schéma d'exploitation optimal. A partir d'un modèle géologique construit sur la base d'au moins les informations issues d'une image sismique obtenue lors des étapes précédentes, les spécialistes peuvent déterminer plusieurs schémas d'exploitation correspondant à différentes configurations 10 possibles d'exploitation du réservoir souterrain : emplacement des puits producteurs et/ou injecteurs, valeurs cibles pour les débits par puits et/ou pour le réservoir, le type d'outils utilisés, les fluides utilisés, injectés et/ou récupérés, etc. Pour chacun de ces schémas, il convient de déterminer leurs prévisions de production. Ces prévisions de production probabilistes peuvent être obtenues au moyen d'un logiciel de simulation d'écoulement ainsi 15 qu'au moyen du modèle numérique de réservoir calé. Une simulation de réservoir est une technique permettant de simuler les écoulements de fluides au sein d'un réservoir au moyen d'un logiciel appelé simulateur d'écoulement. Par exemple, le logiciel PumaFlow (IFP Énergies nouvelles, France) est un simulateur d'écoulement. 20 On définit un ou plusieurs schémas d'exploitation possibles adaptés au modèle géologique étudié. Pour chacun de ces schémas, on détermine les réponses par simulation. A partir des prévisions de productions probabilistes définies pour chaque schéma d'exploitation, les spécialistes peuvent par comparaison choisir le schéma d'exploitation qui 25 leur semble le plus pertinent. Par exemple : - en comparant le maximum du volume d'huile récupéré, on peut déterminer le schéma de production susceptible de fournir le maximum de récupération ou d'être le plus rentable ; - en comparant l'écart type du volume d'huile récupéré, on peut déterminer le schéma 30 de production le moins risqué. 3024243 22 On exploite alors le réservoir selon le schéma d'exploitation défini par exemple en forant de nouveaux puits (producteur ou injecteur), en modifiant les outils utilisés, en modifiant les débits et/ou la nature de fluides injectés, etc. 5 Produit programme d'ordinateur L'invention concerne, par ailleurs, un produit programme d'ordinateur téléchargeable depuis un réseau de communication et/ou enregistré sur un support lisible par ordinateur et/ou exécutable par un processeur. Ce programme comprend des instructions de code de programme pour la mise en oeuvre du procédé tel que décrit ci-dessus, lorsque le 10 programme est exécuté sur un ordinateur. Exemples d'application Le procédé selon l'invention est appliqué sur un cas de données sismiques réelles bi- 15 dimensionnelles après sommation dans le cas de la présence d'un bruit aléatoire. On observe sur la Figure 5A que les données sismiques en question correspondent à une superposition de réflexions sismiques primaires et de réflexions multiples, les réflexions multiples ayant une amplitude très forte par rapport aux réflexions primaires. On peut aussi observer que réflexions primaires et réflexions multiples sont fortement perturbées par un 20 bruit aléatoire important. La Figure 5B présente le résultat de l'application d'un filtrage adaptatif selon l'art antérieur (décrit dans Ventosa et al. (2012)) aux données de la Figure 5A. On observe que le filtrage adaptatif selon l'art antérieur fait correctement apparaître la réflexion primaire d'intérêt dans une grande partie de la zone considérée, mais que la continuité de cette réflexion primaire est perdue dans la partie la plus à gauche de la Figure 25 5B. Or il s'agit de la partie où le réflecteur ayant généré la réflexion d'intérêt est le plus penté. La Figure 5C présente le résultat du procédé selon l'invention appliqué aux données sismiques présentées en Figure 5A. Le procédé selon l'invention a été mis en oeuvre en utilisant la technique décrite dans Chaux et al. (2006) pour la méthode de décomposition et la technique décrite dans Ventosa et al. (2012) pour le filtrage adaptatif. On observe très 30 clairement que le procédé selon l'invention permet d'extraire de manière plus nette la réflexion primaire d'intérêt, en particulier dans la partie gauche de la Figure 5C. Cet avantage est en particulier obtenu par le fait que le procédé selon l'invention applique un filtrage adaptatif sur des données sismiques décomposées selon des orientations dans 3024243 23 l'espace privilégiées, et que la pondération par la concentration relative renforce les ressemblances entre les données sismiques et le modèle de multiples, au détriment des bruits ou perturbations non désirés, produisant une image sismique plus précise des réflexions primaires. Ainsi, même dans le cas de données sismiques présentant un fort bruit aléatoire, la présente invention permet d'améliorer le filtrage adaptatif de réflexions multiples contenues dans des données sismiques, en opérant de manière sélective selon des bandes de fréquence et des gammes d'orientations dans l'espace privilégiées. 5 10 REVENDICATIONS1) Procédé pour construire une image sismique filtrée de réflexions multiples, à partir d'un enregistrement de données sismiques comprenant des réflexions primaires et des réflexions multiples, à partir d'au moins un modèle desdites réflexions multiples, caractérisé en ce que l'on réalise au moins les étapes suivantes : a) on applique une méthode de décomposition pour décomposer, selon N directions de décomposition, lesdites données sismiques en un ensemble de N composantes desdites données sismiques ; b) pour au moins un desdits modèles desdites réflexions multiples, on applique ladite méthode de décomposition pour décomposer, selon lesdites N directions de décomposition, ledit modèle de réflexions multiples en un ensemble de N composantes dudit modèle ; c) pour chaque direction de décomposition et pour au moins un desdits modèles desdites réflexions multiples, on calcule une concentration relative entre ladite composante desdites données sismiques dans ladite direction et ladite composante dudit modèle dans ladite direction ; d) à partir d'au moins un desdits modèles desdites réflexions multiples, on applique une méthode de filtrage adaptatif desdites réflexions multiples sur chacune desdites N composantes desdites données sismiques, et on obtient un ensemble de N composantes desdites données sismiques filtrées ; e) on calcule une recombinaison pondérée desdites composantes desdites données sismiques filtrées, à partir d'une pondération calculée pour chaque direction de décomposition en fonction desdites concentrations relatives calculées dans ladite direction de décomposition. 2) Procédé selon la revendication 1, dans lequel ladite méthode de décomposition est une transformée en ondelettes. 3024243 25 3) Procédé selon l'une des revendications précédentes, dans lequel ladite méthode de décomposition est une transformée en ondelettes M-bandes en arbre dual. 4) Procédé selon l'une des revendications précédentes, dans lequel ledit filtrage 5 adaptatif est un filtre de Wiener unaire en trame d'ondelettes complexes. 5) Procédé selon l'une des revendications précédentes, dans lequel ladite concentration relative entre un signal S et un signal S' bi-dimensionnels est calculée selon la formule suivante : s21 10 51,52 SNIS'l 2 2 , S2 )1 CR(S(Si , S2), S'(Sp S2 )) SI,S2 s2 )12 SI,S2 où si est le nombre d'échantillons dans une direction de l'espace et s2 est le nombre d'échantillons dans l'autre direction de l'espace. 6) Procédé selon l'une des revendications précédentes, dans lequel ladite 15 recombinaison pondérée est calculée de la façon suivante : SR =aIen. D' (0 , , SDF , ,0) n=1,N où SDF, est ladite composante dans ladite direction n (n=1,1V) desdites données sismiques filtrées, a est une constante, en est ladite pondération pour ladite direction de décomposition n, et D' est une inverse de ladite méthode de décomposition. 20 72 7) Procédé selon l'une des revendications 1 à 5, dans lequel ladite recombinaison pondérée est calculée de la façon suivante : 3024243 26 SR = crI D' (0,..,en.SDF',...,0) n=1 , N où SDF', est ladite composante dans ladite direction n (n=1,N) desdites données sismiques filtrées, a est une constante, en est ladite pondération pour ladite direction 5 de décomposition n, et D' est une inverse de ladite méthode de décomposition. 8) Procédé selon l'une des revendications 6 à 7, dans lequel on connait de manière approchée ladite inverse de ladite méthode de décomposition, et dans lequel on ajoute auxdites composantes desdites données sismiques, un résidu correspondant à 10 la différence entre lesdites données sismiques et le résultat de ladite inverse approchée. 9) Procédé selon l'une des revendications 6 à 8, dans lequel ladite pondération pour ladite direction de décomposition n (n=1,N) est calculée selon la formule suivante : 15 en = i =1 , I I où en est une pondération pour ladite direction n et un modèle de réflexions multiples i 10) Procédé selon l'une des revendications 6 à 8, dans lequel ladite pondération pour 20 ladite direction de décomposition n (n=1,N) est calculée selon la formule suivante : max(e ) " =1,/ " où en est une pondération pour ladite direction n et un modèle de réflexions multiples i 3024243 27 11) Procédé selon l'une des revendications 9 à 10, dans lequel ladite pondération pour ladite direction de décomposition n et ledit modèle de réflexions multiples i (avec i=i,/ et n=1,N) est calculée selon la formule suivante : 1 n = W(x) , avec x = min(CR , ) , 5 où W est une fonction croissante, CR'' est ladite concentration relative entre ladite composante desdites données sismiques dans ladite direction n et ladite composante dudit modèle i dans ladite direction n. 12) Procédé selon la revendication 11, dans lequel ladite fonction croissante W est 10 définie telle que pour un x donné, W(x) = x' avec p O. 13) Procédé selon la revendication 12, dans lequel p vaut 4. 14) Procédé d'exploitation d'une formation souterraine, caractérisé en ce qu'on réalise les étapes suivantes : 15 on construit une image sismique filtrée des réflexions multiples au moyen du procédé tel que décrit selon l'une des revendications 1 à 13 ; on construit un modèle géologique représentatif de la formation étudiée à partir d'au moins l'image sismique ainsi déterminée ; on détermine un schéma d'exploitation optimal du réservoir à partir du modèle 20 géologique ainsi déterminé ; on exploite ledit réservoir en mettant en oeuvre ledit schéma d'exploitation optimal. 15) Produit programme d'ordinateur téléchargeable depuis un réseau de communication 25 et/ou enregistré sur un support lisible par ordinateur et/ou exécutable par un processeur, comprenant des instructions de code de programme pour la mise en oeuvre du procédé selon l'une des revendications précédentes, lorsque ledit programme est exécuté sur un ordinateur.