FR 3011939 A1 20150417 FR 1360016 A 20131015 La présente invention concerne le domaine de l'industrie pétrolière, plus particulièrement le domaine de la caractérisation et l'exploitation des réservoirs pétroliers ou encore de la caractérisation et de la surveillance de sites de stockage géologique de gaz. Le domaine de la surveillance de sites de stockage de gaz à effet de serre tels que le dioxyde de carbone (002) est en particulier concerné par cette invention. La présente invention propose une méthode permettant d'évaluer la taille moyenne des hétérogénéités du sous-sol, qu'elles soient liées à la matrice rocheuse ou bien à la répartition de fluides dans le sous-sol, à partir d'enregistrements d'ondes sismiques, quelle que soit la longueur d'onde considérée et y compris dans le cas de milieux viscoélastiques. Les milieux géologiques se caractérisent par de nombreuses échelles d'hétérogénéités imbriquées. Chacune de ces hétérogénéités est susceptible de modifier à son échelle l'écoulement des fluides dans le sous-sol. Aussi une meilleure connaissance de la répartition spatiale des hétérogénéités dans le sous-sol et de leurs dimensions caractéristiques contribuent à une meilleure connaissance des performances du réservoir, tant dans le domaine pétrolier que dans le domaine de l'environnement (stockage géologique de gaz à effet de serre par exemple), et à l'optimisation de son exploitation. État de la technique Dans le cadre de la présente invention, au moins deux types possibles d'hétérogénéités d'une formation souterraine sont considérés. Le premier type est dû aux variations spatiales des propriétés physiques d'une formation souterraine. Elles sont liées au fait que, même en l'absence de fluides dans les pores de la roche constituant la formation, la matrice de cette roche n'est en général pas parfaitement homogène et peut présenter des variations spatiales. Le deuxième type d'hétérogénéités est dû à la répartition spatiale des fluides dans une formation souterraine. En effet, si l'on suppose la matrice rocheuse fixée et si la formation est en cours d'exploitation, par injection de fluides venant se substituer aux fluides initialement en place dans cette formation, le réseau poreux de la formation est rempli de différents fluides dont la répartition spatiale n'est le plus souvent pas parfaitement homogène. Les documents suivants seront cités au cours de la description : Bourbié, T., Coussy, 0. and B. Zinszner, 1987, Acoustics of Porous Media, Editions Technip, Paris. Mukerji, T. and Mavko, G., 1994, Scale-dependent dynamic wave propagation in heterogeneous media: Il. theory. SEG Technical Program Expanded Abstracts 1994: pp. 1151-1154. Molyneux, J.B., et D.R. Schmitt, 1999, First-break timing: Arrivai onset times by direct correlation, Geophysics, 65(5), 1492-1501. Rasolofosaon, P. et B. Zinszner, 2004, "Laboratory petroacoustics for seismic monitoring feasibility study", The Leading Edge, 23(3), 252-258. Sato, H. et M. Fehler, 1998, Seismic Wave Propagation and Scattering in the Heterogeneous Earth: Springer-Verlag and American lnstitute of Physics Press). Virieux, J., et V. Farra, 1991, Ray tracing in 3-D complex isotropic media: An analysis of the problem, Geophysics, 56, 2057-2069. Wu, R. S., et K. Aki, 1988, Introduction: Seismic waves scattering in three-dimensionally heterogeneous earth, in: Wu, R-S., et Aki, K. (eds.). Scattering and attenuation of seismic waves: Pure and Applied Geophysics,128: 1-6. Zinszner, B. and F.M. Pellerin, 2007, "A geoscientist's guide to petrophysics", Editions Technip, Paris. La vitesse de propagation des ondes sismiques fait partie des attributs couramment utilisés pour caractériser un milieu donné. Dans le cas d'un milieu homogène, la vitesse de propagation des ondes sismiques est un invariant quel que soit le point de mesure considéré. Dans un milieu hétérogène, la vitesse de propagation des ondes va non seulement varier en fonction du point de mesure, mais il est aussi possible de mesurer des vitesses de propagation des ondes différentes pour un même point suivant la méthode d'estimation de la vitesse utilisée. En effet, pour peu que les longueurs d'onde sismique utilisées soit petites par rapport à la taille des hétérogénéités, un phénomène de dispersion de cheminement (aussi connue sous le nom de "path dispersion" en anglais) va se produire à la rencontre de ces hétérogénéités, une partie de l'énergie de l'onde pouvant prendre un chemin plus rapide et ainsi générer une arrivée précoce dans les enregistrements sismiques. Dans un tel cas, on observe alors une différence entre la vitesse dite d'impetus (aussi connue sous le nom de "first break velocity" en anglais, ou encore de "vitesse de première arrivée") et par exemple la vitesse de phase (Rasolofosaon et Zinszner, 2004 ; Figure 1). Toutefois, l'estimation quantitative de la taille des hétérogénéités n'a à ce jour jamais été réalisée à partir de l'observation de la différence entre vitesse de phase et d'impetus. Un certain nombre de références théoriques et expérimentales existent sur la propagation des ondes élastiques en milieux hétérogènes géologiques, et ce en vue de leur caractérisation (par exemple Wu et Aki, 1985 ; Sato et Fehler, 1998). Ces travaux, principalement dans le domaine de la sismologie, sont basés sur l'étude de la "queue" (ou coda) du signal sismique. En effet, après propagation dans un milieu hétérogène, supposé d'extension infinie, le signal sismique se divise en deux composantes, à savoir une composante majeure cohérente correspondant à l'onde directe, suivie d'arrivées plus tardives, étalées dans le temps et incohérentes, résultant des rebonds de l'onde sismique sur les différentes hétérogénéités du milieu de propagation. D'après Sato et Fehler (1998), K. Aki a été le premier à étudier le lien entre l'étude des codas sismiques et la diffraction induite par les hétérogénéités de la lithosphère terrestre. Toutefois, les parties tardives des enregistrements sismiques sont souvent contaminées par des réflexions multiples dans les différentes couches du sous-sol. Ceci peut fausser considérablement le traitement des signaux sismiques et aboutir à une caractérisation notablement erronée des hétérogénéités du sous-sol. Parmi les études basées sur les premières arrivées sismiques et non sur la queue du signal, un grand nombre de théories existent sur la propagation dans les milieux aléatoires divers (gaussiens, exponentiels, etc) donnant, entre autres, l'écart entre la vitesse à très basse fréquence VLF et la vitesse à très haute fréquence VHF . Dans le cas d'un milieu élastique aléatoire gaussien, présentant des fluctuations de vitesse et non de densité, par exemple Mukerji et Mavko (1994) ont montré le résultat suivant : - .'I (1) (VHF -VLF )1 Vo KSI 2 ( 1 (DIMA 1)D 2 où Vo désigne la vitesse moyenne dans le milieu considéré autour de laquelle existent des fluctuations aléatoires dont l'écart quadratique moyen normalisé par Vo vaut KSI , D désigne la distance entre les deux points de mesure entre lesquels les vitesses VLF et VHF sont calculées, DIM désigne la dimension de l'espace du problème considéré (il vaut donc 1, 2 ou 3 suivant que la propagation d'onde sismique est respectivement uni-, bi- ou tri-dimensionnelle) et A est la longueur de corrélation sous cette hypothèse gaussienne, assimilable à la taille moyenne des hétérogénéités. Cette technique permettant de travailler sur la partie initiale des enregistrements sismiques, la moins contaminée par les arrivées multiples internes tardives, assure une meilleure estimation des vitesses. Toutefois, les descriptions théoriques de la propagation des ondes élastiques dans les milieux aléatoires hétérogènes sont limitées à des milieux parfaitement élastiques non intrinsèquement dissipatifs, en d'autre termes non viscoélastiques. Ce caractère intrinsèquement dissipatif, au sens viscoélastique, est pourtant une caractéristique universelle des matériaux naturels, en particulier des milieux géologiques tels que les roches. L'objet de la présente invention est une méthode pour estimer la taille moyenne des hétérogénéités propres à une formation souterraine, ou bien la taille moyenne des hétérogénéités de répartition en fluides induites par une substitution de fluide dans cette formation, et ce quelle que soit l'échelle d'investigation. De plus, la présente invention permet de tenir compte des effets viscoélastiques. La méthode selon l'invention Ainsi, la présente invention concerne une méthode pour estimer la taille moyenne des hétérogénéités d'une formation souterraine, liées à la matrice rocheuse de ladite formation ou bien à la répartition de fluides dans ladite formation, au moyen de mesures d'ondes sismiques. La méthode comporte les étapes suivantes : - on détermine au moins un groupe de points d'intérêt au sein de ladite formation pour lesquels on dispose d'une information sismique associée ; - on calcule la distance D de propagation des ondes sismiques entre lesdits points de chaque groupe ; - on détermine la vitesse d'impetus VFB entre lesdits points de chaque groupe ; - on détermine la vitesse de phase VPH entre lesdits points de chaque groupe ; - on évalue la taille moyenne HS des hétérogénéités réparties entre lesdits points de chaque groupe, à partir de ladite distance D, de ladite vitesse d'impetus VFB et de ladite vitesse de phase VPH. Selon l'invention, on peut effectuer une étape supplémentaire pour estimer l'ordre de grandeur des fluctuations moyennes des vitesses sismiques KSI dans ladite formation. Dans ce cas, ladite taille moyenne HS des hétérogénéités réparties entre lesdits points de chaque groupe pourra être calculée par la formule suivante : D(DIM -1),F HS ,-- avec D correspondant à ladite distance entre lesdits points d'un même groupe, DIM à la dimension de l'espace du problème considéré, VFB à la ladite vitesse d'impetus entre lesdits points de chaque groupe, VPH à ladite vitesse de phase entre lesdits points de chaque groupe, et KSI auxdites fluctuations moyennes des vitesses sismiques dans ladite formation. Selon l'invention, lorsque les mesures d'ondes sismiques peuvent être réalisées : - au moyen d'ondes ultrasonores : lesdits points d'intérêt d'un même groupe peuvent correspondre à la position d'une source et à la position d'un récepteur ; - par diagraphies soniques : lesdits points d'intérêt d'un même groupe peuvent correspondre à deux récepteurs distincts ; - par de la sismique de puits en mode transmission et à déport nul : lesdits points d'intérêts d'un même groupe peuvent correspondre à deux récepteurs distincts ; {2 VFB -VPH 1 VPH (KSI)` , - par de la sismique de puits en mode réflexion : lesdits points d'intérêt d'un même groupe peuvent correspondre à deux points situés sur deux réflecteurs sismiques différents ; - par de la sismique de puits à puits : lesdits points d'intérêt d'un même groupe peuvent correspondre à la position d'une source et à la position d'un récepteur ; - par de la sismique réflexion de surface : lesdits points d'intérêt peuvent être situés sur deux réflecteurs sismiques différents. Selon l'invention, dans le cas d'une estimation de la taille des hétérogénéités liées à la répartition de fluides induite par une substitution de fluide dans la formation étudiée, les fluctuations moyennes des vitesses sismiques KSI peuvent correspondre à la variation de vitesse induite par cette substitution de fluide. Selon une application particulière de l'invention, le CO2 constitue le fluide de substitution. D'autres caractéristiques et avantages de la méthode selon l'invention apparaîtront à la 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. Présentation succincte des Figures La Figure 1 illustre la différence observée entre les temps d'arrivée TVFB et TVPH respectivement associés à la vitesse d'impetus VFB et à la vitesse de phase VPH dans le cas d'un milieu hétérogène. La Figure 2 illustre l'étape de sélection des points d'intérêt de la présente invention ; 2A : dans le cas d'une acquisition sismique de surface à déport nul (source S et récepteur R confondus), les points d'intérêt d'un même groupe seront au nombre de deux (M1 et M2) et localisés sur deux réflecteurs REF1 et REF 2 différents (M1 sur REF1 et M2 sur REF2) ; 2B : dans le cas d'une acquisition sismique de surface à déport non nul (source S et récepteur R distincts), les points d'intérêt seront au nombre de trois (M1, M2 et M3), deux localisés sur un premier réflecteur (M1 et M3 sur REF1) et un troisième localisé sur un deuxième réflecteur (M2 sur REF2) ; 2C : dans le cas d'une sismique de puits en mode transmission à déport nul (source S en surface très proche du puits W et série de récepteurs R dans le puits W), les points d'intérêt d'un même groupe sont au nombre de deux et correspondent à deux positions récepteur du dispositif (R1 et R2) ; 2D : dans le cas d'une sismique de puits en mode réflexion à déport non nul (source S dans le puits W et série de récepteurs R en surface), les points d'intérêt seront au nombre de 3 (M1, M2 et M3), deux localisés sur un premier réflecteur (M1 et M3 sur REF1) et un troisième localisé sur un deuxième réflecteur (M2 sur REF2) ; 2E : dans le cas d'une sismique de puits à puits en mode transmission (série de sources S dans un puits WS et séries de récepteurs R dans un puits WR), les points d'intérêt d'un même groupe sont au nombre de deux et correspondent à une position source et à une position récepteur du dispositif. La Figure 3 montre un abaque, établi grâce à la présente invention dans le cas tridimensionnel, donnant les variations de l'écart relatif (VFB - VPH)/VPH en pourcentage entre la vitesse d'impetus VFB et la vitesse de phase VPH , en fonction de la distance D de propagation de l'onde sismique rapportée à la taille moyenne HS des hétérogénéités, pour différentes valeurs de l'amplitude KSI des fluctuations de vitesses sismiques, allant de 1% à 10% par pas de 1%. Le Tableau 1 résume les résultats obtenus par la méthode sur quatre roches différentes, deux de type calcaire (dits d'Euville et de Brauvilliers) et deux de type grès (dits de Fontainebleau et de Meule). Les différentes colonnes du Tableau 1 listent, de gauche à droite, le nom et le type de roche, l'écart relatif (VFB - VPH)/ VPH , exprimé en pourcentage, entre la vitesse d'impetus VFB et la vitesse de phase VPH à forte pression de confinement (typiquement 30 à 40MPa), l'amplitude KSI des fluctuations de vitesses ultrasonores dans la roche, la distance D de propagation de l'onde sismique rapportée à la taille moyenne HS des hétérogénéités, la taille moyenne HS des hétérogénéités déterminée à partir de la Figure 3 ou de l'équation 3, et la taille moyenne HS des hétérogénéités déterminée par méthodes indépendantes, soit à partir d'images de tomographie de densité par rayons X ou par des lames minces, comme illustré par la Figure 4. La Figure 4 montre les images de tomographie de densité par rayons X et une photographie de lame mince des quatre roches choisies pour illustrer la méthode et dont les caractéristiques sont résumées dans le Tableau 1. Description détaillée de la méthode Les définitions suivantes sont utilisées au cours de la description de l'invention : vitesse d'impetus : il s'agit de la vitesse mesurée à partir du pointé d'impetus, qui correspond au premier signal cohérent ayant une amplitude supérieure au niveau de bruit moyen (par exemple Molyneux et al. 1999). vitesse de phase : il s'agit de la vitesse à laquelle un point particulier (tel un pic ou un creux) de l'onde se propage dans l'espace (par exemple Bourbié et al., 1987). ondes sismiques = on désigne par ce terme des ondes élastiques caractérisées par une longueur d'onde quelconque, qui pourra aller du domaine ultrasonore (longueurs d'ondes de l'ordre du millimètre) à celui de la sismologie (longueurs d'ondes de l'ordre de plusieurs kilomètres). taille moyenne des hétérogénéités = on désigne par ce terme un paramètre dimensionnel représentatif du diamètre moyen des hétérogénéités d'une formation, que ce soient des hétérogénéités liées à la matrice rocheuse de la formation, ou bien des hétérogénéités liées à la répartition des fluides après substitution de fluide dans cette formation. L'objet de la présente invention est une méthode pour estimer la taille moyenne d'hétérogénéités propres à une formation souterraine ou bien induites par une substitution de fluide dans cette formation, et ce quelle que soit l'échelle d'investigation sismique considérée et, qui plus est, permettant de tenir compte des effets viscoélastiques. Cette invention est basée sur l'exploitation de la différence entre vitesse de phase et vitesse d'impetus mesurées entre des points situés dans la formation. La présente invention comporte au moins les étapes suivantes : a) on détermine au moins un groupe de points d'intérêt au sein de la formation pour lesquels on dispose d'une information sismique associée ; b) on calcule la distance D de propagation des ondes sismiques entre les points de chaque groupe ; c) on détermine la vitesse d'impetus VFB entre les points de chaque groupe ; d) on détermine la vitesse de phase VPH entre les points de chaque groupe ; e) on évalue la taille moyenne HS des hétérogénéités réparties entre les points de chaque groupe, à partir de la distance D, de la vitesse d'impetus VFB et de la vitesse de phase VPH. Définition des groupes de points d'intérêts Il s'agit lors de cette première étape de définir les groupes de points entre lesquels on souhaite estimer la taille moyenne des hétérogénéités. Ces points, localisés au sein ou à la frontière de la formation d'intérêt, doivent être caractérisés par une information sismique directe ou indirecte. Cela signifie que ces points doivent avoir été traversés par une onde sismique, soit en tant qu'émetteurs ou récepteurs de cette onde, ou bien en tant que points se trouvant sur le passage de cette onde (en réflexion ou bien en transmission). Selon un mode réalisation de la présente invention, la méthode d'investigation utilisée correspondra à des mesures de laboratoire ultrasonores. Un dispositif de mesure ultrasonore en laboratoire comporte toujours au moins deux transducteurs pouvant jouer le rôle d'émetteur comme celui de récepteur. Pour ce mode de réalisation, le groupe de points d'intérêt sera typiquement constitué de deux points, l'un étant un point émetteur du dispositif et l'autre un point récepteur de ce même dispositif. Selon un autre mode réalisation de la présente invention, la méthode d'investigation utilisée correspondra à des diagraphies soniques. Les outils de diagraphies soniques comportent au minimum un émetteur et deux récepteurs, mais il existe aussi des dispositifs avec des émetteurs monopole (axisymétrique) ou encore dipôle (orientés) et une centaine de récepteurs. Pour ce mode de réalisation, l'étape a consistera à choisir deux points ayant correspondu à la position de deux récepteurs distincts pour une position émetteur donnée. Selon un autre mode de réalisation de la présente invention, la méthode d'investigation utilisée correspondra à la sismique réflexion de surface. Pour ce mode de réalisation, un générateur émet des ondes sismiques depuis la surface (ou légèrement sous la surface) et les signaux réfléchis par les discontinuités géologiques, aussi appelées réflecteurs sismiques, sont enregistrés par un réseau de capteurs placés également en surface (ou légèrement sous la surface). Pour ce mode de réalisation, les points d'intérêt d'un même groupe seront localisés sur deux réflecteurs sismiques différents, correspondant par exemple au toit et à la base de la formation à étudier. Dans le cas de données sismiques à déport nul (aussi connues sous le nom de données "stack" ou "poststack" en anglais ou encore de données sismiques après sommation), et dans le cas où les réflecteurs sismiques ne présentent pas de pendage, les ondes sismiques se propagent de façon verticale vers le bas puis vers le haut après réflexion. Dans ce cas, les points d'intérêts d'un même groupe seront au nombre de deux et seront localisés à la verticale l'un de l'autre, l'un correspondant à un point de transmission au toit de la formation et l'autre au point de réflexion sur la base de la formation (Figure 2A). Dans le cas de données sismiques multi-déports (aussi connues sous le nom de données "prestack" en anglais, ou encore de données avant sommation), les groupes de points seront constitués de trois points, deux étant placés sur le réflecteur sismique correspondant au toit de la formation et qui voient l'onde se transmettre, et le troisième point placé sur le réflecteur sismique correspondant à la base de la formation et qui voit l'onde se réfléchir (Figure 2B). Dans le cas où la formation étudiée et son encaissant génèrent des propagations d'ondes sismiques complexes (par exemple lorsque les réflecteurs sismiques sont non plans et horizontaux et/ou lorsque les couches géologiques sont caractérisées par des variations internes de vitesse), il pourra être nécessaire de faire appel à une méthode de simulation de la propagation des ondes sismiques telle que le tracé de rayons (par exemple Virieux et Farra, 1991) pour localiser précisément dans la formation les points d'intérêt de chaque groupe. A noter qu'une telle méthode nécessite une représentation en vitesse de la formation étudiée et son encaissant. Selon un autre mode de réalisation de la présente invention, la méthode d'investigation utilisée correspondra à la sismique de puits. On peut distinguer trois types de géométrie d'acquisition pour la sismique de puits : - sources en surface avec récepteurs dans le puits (Figure 2C), - sources dans le puits avec récepteurs en surface (Figure 2D), sources dans un puits avec récepteurs dans un autre puits (Figure 2E). Par ailleurs, on peut s'intéresser en sismique de puits aussi bien aux arrivées sismiques directes (mode dit "en transmission") qu'aux arrivées réfléchies (mode dit "en réflexion"). Dans le cas de la sismique de puits en mode transmission et dans le cas d'un déport nul ou quasi-nul (cas des PSV ou "Profils Sismiques Verticaux"), les groupes de points seront constitués de deux points, correspondant à deux récepteurs distincts du dispositif (Figure 2C). Dans le cas de la sismique de puits en mode réflexion, les groupes de points seront comme pour la sismique réflexion de surface localisés sur deux réflecteurs sismiques d'intérêt et seront au nombre de deux ou trois suivant le déport considéré (Figure 2D). De même que dans le cas de la sismique réflexion de surface, on pourra faire appel à une méthode de simulation de la propagation des ondes sismiques pour localiser précisément dans la formation les points d'intérêt sur les réflecteurs sismiques. Dans le cas de la sismique de puits à puits en mode transmission, les deux points d'un groupe correspondront à une position source et à une position récepteur du dispositif (Figure 2E). A noter que dans le cas où la formation étudiée et son encaissant génèrent des propagations d'ondes sismiques complexes (par exemple lorsque les réflecteurs sismiques sont non plans et horizontaux et/ou lorsque les couches géologiques sont caractérisées par des variations internes de vitesse), il pourra être nécessaire de faire appel à une méthode de simulation de la propagation des ondes sismiques pour localiser précisément dans la formation les points d'intérêt de chaque groupe. Calcul de la distance D pour chaque groupe de points Lors de cette seconde étape, il s'agit de calculer la distance D parcourue par l'onde sismique entre les points de chacun des groupes précédemment définis. Lorsque l'approximation sera justifiée, on considèrera que la trajectoire de l'onde sismique entre les points d'un groupe donné correspond à des segments de droite. On estimera alors la distance D comme étant la somme des longueurs de chacun de ces segments. A noter que dans le cas de la sismique réflexion (de puits ou de surface), l'onde sismique effectuant un trajet descendant (entre émission et réflexion) puis montant (entre réflexion et réception), dans un cas où les points d'intérêt d'un groupe sont au nombre de deux (cas de la sismique réflexion à déport nul), la distance D correspondra à deux fois la distance entre les points d'intérêt d'un même groupe. Dans les cas où l'approximation de la trajectoire de l'onde sismique entre les points d'intérêt d'un même groupe par des segments de droite est trop erronée, pour peu que l'on dispose d'une représentation en vitesses du sous-sol et d'un simulateur de la propagation des ondes comme par exemple un tracé de rayons (par exemple Virieux et al., 1988), on pourra plus précisément calculer la trajectoire de l'onde sismique entre les points d'un groupe donné puis la distance D. On aura par exemple recours à une méthode précise d'évaluation de la trajectoire de l'onde sismique dans le cas de la sismique réflexion 3D appliquée à des milieux dont les vitesses de propagation des ondes sismiques varient latéralement, alors qu'une approximation par segments de droite de la trajectoire de l'onde sismique sera majoritairement suffisante dans le cas de diagraphies acoustiques. Détermination de la vitesse d'impetus VFB pour chaque groupe de points Il s'agit dans cette étape de déterminer la vitesse d'impetus ou encore la vitesse de première arrivée entre les points d'intérêt d'un même groupe de points. Il s'agit d'une mesure bien connue. On pourra par exemple procéder selon la méthode définie dans Molyneux et al. (1999). Détermination de la vitesse de phase VPH Il s'agit de déterminer la vitesse de phase VPH entre les points d'intérêt d'un même groupe de points. Il s'agit là aussi d'une mesure bien connue. On pourra par exemple utiliser la méthode de la détermination de la vitesse de phase par rapport spectral telle que décrite dans Bourbié et al. (1987) ou encore Rasolofosaon et Zinszner (2004). Évaluation de la taille moyenne HS des hétérogénéités Tous les matériaux naturels, en particulier les milieux géologiques tels les roches, sont plus ou moins hétérogènes et intrinsèquement dissipatifs, au sens viscoélastique. De manière à intégrer au mieux ces deux caractères et développer la présente méthode, les hypothèses suivantes sont faites : on suppose que le milieu analysé est hétérogène et dissipatif. Si fc désigne approximativement la fréquence centrale des signaux sismiques mesurés, on suppose que la vitesse VPH considérée dans l'étape (4) est une estimation de la vitesse moyenne de phase autour de la fréquence fc incluant tous les effets viscoélastiques ainsi que les effets d'hétérogénéités (mais homogénéisés jusqu'à cette fréquence) dans le milieu viscoélastique homogène équivalent. Il est à noter que l'on suppose, d'une part, que toutes les hétérogénéités induisant des variations de vitesses jusqu'à la fréquence fc peuvent être intégrées dans un modèle viscoélastique homogène équivalent, et d'autre part, que l'on fait totalement abstraction de la vitesse de phase à très basse fréquence, qui est d'ailleurs inaccessible à la mesure dans la plupart des cas ; on suppose que des fluctuations de vitesse, inhérentes à tous les milieux naturels, existent autour de la vitesse de référence, choisie comme correspondant à la vitesse VPH. Ces fluctuations de vitesses vont induire une augmentation de vitesse entre la vitesse à la fréquence fc et la vitesse à très haute fréquence ; en découplant les effets dus à la présence d'hétérogénéités des effets viscoélastiques (négligés entre la fréquence fc et la fréquence infinie), et en assimilant la vitesse d'impetus VFB à la vitesse à très haute fréquence, l'écart de vitesse entre VPH et VFB peut être interprété comme le saut de vitesse dû aux hétérogénéités depuis les très basses fréquences jusqu'aux très hautes fréquences, et non uniquement de la fréquence fc à la fréquence infinie. En effet le saut de vitesse dû aux hétérogénéités des très basses fréquences à la fréquence fc est supposé nul, ou plus exactement a déjà été intégré dans le modèle viscoélastique homogène équivalent. Par conséquent, si l'on fait le parallèle avec l'équation (1), le premier membre de l'équation (1) peut être identifié à : (VFB - VPH )/VPH (2) le paramètre D de l'équation (1) s'identifie immédiatement au paramètre D tel que défini à l'étape b de la présente invention, et la longueur de corrélation A de l'équation (1) s'identifie à la taille moyenne HS des hétérogénéités. En négligeant les fluctuations de densité et si l'on suppose l'amplitude KSI des fluctuations de vitesses sismiques connue, la taille moyenne HS des hétérogénéités réparties entre les points d'intérêts d'un groupe est finalement donnée par la relation : D(DIM -1).jr HS Exploitation des résultats En vue d'une exploitation systématique de l'information contenue dans la différence entre vitesse de phase et vitesse d'impetus mesurées entre les points d'un même groupe, on peut grâce à l'équation (3) établir un abaque très général auquel l'homme du métier pourra se référer pour estimer la taille moyenne des hétérogénéités de la formation souterraine d'intérêt. Un tel abaque est présenté en Figure 3. Cet abaque présente les variations de l'écart relatif entre les vitesses de phase et d'impetus (VFB - VPH)/VPH , exprimé en pourcentage, en fonction de la distance D de propagation de l'onde sismique rapportée à la taille moyenne HS des hétérogénéités, et ce pour différentes valeurs de l'amplitude KSI des fluctuations de vitesses sismiques. 12VFB -VPH VPH (KSI)2 (3) En pratique cet abaque est utilisé de la manière suivante. On part des mesures de la vitesse d'impetus VFB et de la vitesse de phase VPH , et par conséquent de leur écart relatif (VFB - VPH)/VPH , exprimé en pourcentage, sur l'axe des ordonnées de la Figure 3. Ensuite, un intervalle de valeurs de l'amplitude KSI des fluctuations des vitesses sismiques, exprimée en pourcentage, est supposé connu (voir ci-dessous). L'intervalle de valeurs de l'amplitude KSI des fluctuations de vitesses sismiques étant fixé, l'intervalle de valeurs du rapport D/HS est lu sur l'axe des abscisses de l'abaque de la Figure 3 à partir de la valeur mesurée du rapport (VFB - VPH)/ VPH et lue sur l'axe des ordonnées de l'abaque de la Figure 3. Enfin, D ayant été établi lors de l'étape b, on déduit de l'intervalle de valeurs du rapport D/HS , un intervalle de valeurs pour la taille moyenne HS des hétérogénéités de la formation considérée. Calibration des fluctuations de vitesse KSI La mise en oeuvre de la présente invention nécessite d'avoir une estimation de l'amplitude KSI des fluctuations de vitesses sismiques. Cette estimation de l'ordre de grandeur des fluctuations des vitesses sismiques dans la formation d'intérêt pourra résulter d'études antérieures sur la formation en question. Qu'il s'agisse de réservoirs pétroliers ou bien de sites de stockage géologique de gaz, les ordres de grandeur des vitesses de propagation des ondes dans les réservoirs sont souvent bien connus de l'homme du métier. Mais cet ordre de grandeur pourra aussi être déterminé de façon plus précise, notamment à une échelle appropriée à celle de la méthode d'investigation utilisée pour la mise en oeuvre de l'invention. Dans ce cas, cette estimation pourra faire l'objet d'une étape supplémentaire. Ainsi une valeur de KSI pourrait être obtenue par calibration sur des cas de référence pour lesquels on dispose à la fois de mesures indépendantes du rapport (VFB - VPH)/VPH et d'estimation de la taille moyenne HS des hétérogénéités. Par exemple dans le cas de mesures de laboratoire dans le domaine ultrasonore sur un échantillon de roche donné, on pourra par exemple estimer de manière indépendante la taille des hétérogénéités par tomographie de densité par rayons X ou par des lames minces ((Zinszner et Pellerin, 2007) ou (Bourbié et al.,1987)) pour les échantillons de roches en question. Puis à partir de cette estimation de la taille moyenne des hétérogénéités et des mesures de D, VPH et VFB telles que décrites dans les étapes b à d de la présente invention, et à l'aide de l'abaque décrit précédent, on pourra en déduire une estimation des fluctuations de vitesse KSI dans la formation d'intérêt. Dans le cas de diagraphies soniques, on pourra par exemple estimer la vitesse sismique entre chaque couple de récepteurs du dispositif et ce pour différentes positions de source localisées au sein de la formation. Une estimation raisonnable de KSI sera alors l'écart quadratique moyen sur cet intervalle des fluctuations de la vitesse sismique autour de la moyenne des vitesses ainsi estimées. Dans le cas de la sismique réflexion 3D et si par exemple on dispose en plus de données de sismique de puits sans déport (de type Profil Sismique Vertical, noté aussi PSV ci-après), on pourra procéder par exemple de la manière décrite ci-après. A partir des enregistrements PSV des ondes sismiques directes en des points d'enregistrement répartis en profondeur à intervalle fin par rapport à la longueur d'onde dominante de la sismique 3D, on pourra faire une mesure de la vitesse d'impetus entre chaque couple de récepteurs successifs. Cette vitesse pourra être considérée comme l'équivalent d'une donnée de pseudo-diagraphie sonique, mais à l'échelle de la sismique de puits (légèrement plus fine que l'échelle de la sismique 3D). Une estimation raisonnable de KSI sera alors, comme dans le cas des diagraphies soniques, l'écart quadratique moyen des fluctuations de la vitesse d'impetus autour de la moyenne de la vitesse d'impetus sur l'intervalle choisi. Dans le cas d'un suivi d'exploitation de réservoir, de stockage de gaz d'hydrocarbures ou bien de gaz acides, l'hétérogénéité est supposée être essentiellement due à la répartition non-uniforme des fluides saturant (en négligeant l'hétérogénéité de la roche elle-même). Ainsi, une valeur raisonnable de KSI pourra être l'écart relatif entre la vitesse sismique moyenne mesurée dans une zone saturée du fluide en place avant l'exploitation, et la vitesse sismique moyenne mesurée dans une zone saturée du deuxième fluide venant se substituer au fluide d'origine au cours de l'exploitation du réservoir. Exemples Les caractéristiques et avantages de la méthode selon l'invention, apparaîtront plus clairement à la lecture des exemples d'application ci-après. Le Tableau 1 résume les mesures effectuées au laboratoire aux fréquences ultrasonores sur quatre échantillons de roches, deux de type calcaire (dits d'Euville et de Brauvilliers) et deux de type grès (dits de Fontainebleau et de Meule). Les mesures de vitesse d'impetus VFB et la vitesse de phase VPH ont été réalisées à forte pression de confinement (typiquement 30 à 40MPa) de manière à assurer que les effets observés ne sont pas dus à la présence de défauts mécaniques (microfissures, joints de grains, micro-fractures, etc), mais à des hétérogénéités intrinsèques aux échantillons. La distance de propagation D , correspondant à la distance entre source et récepteur du dispositif, ou encore à la taille de l'échantillon de roche, est de 5 cm. La taille moyenne HS des hétérogénéités a été par ailleurs estimée de manière indépendante par d'autres mesures, soit à partir d'images de tomographie de densité par rayons X, soit par des lames minces (voir Figure 4 et dernière colonne du Tableau 1). Dans le cas du calcaire d'Euville, l'écart relatif (VFB - VPH)/ VPH estimé vaut environ 6%. Pour des valeurs de l'amplitude KSI des fluctuations de vitesses ultrasonores dans la roche comprise entre 4% et 5%, le rapport D/HS déduit de l'abaque de la Figure 3 est compris entre 13 et 21. Connaissant D, on en déduit alors la taille moyenne des hétérogénéités, comprise entre 2,4 mm et 3,8 mm. Ceci est tout à fait en accord avec la taille moyenne des hétérogénéités de densité mesurée par tomographie à rayons X, à savoir 2 mm à 5 mm (voir Figure 4 et dernière colonne du Tableau 1). Le second exemple est le calcaire de Brauvilliers. L'écart relatif (VFB - VPH)/ VPH mesuré est compris entre 2% et 7%, suivant les différents échantillons de cette roche. Pour une estimation de KSI entre 4% et 5%, le rapport D/HS déduit de l'abaque de la Figure 3 est compris entre 4 et 25. D valant 5 cm, la taille moyenne estimée des hétérogénéités est alors comprise entre 2 mm et 12,5 mm, ce qui est cohérent avec la taille moyenne des hétérogénéités de densité mesurée par tomographie à rayons X, à savoir 4 mm à 7 mm (voir Figure 4 et dernière colonne du Tableau 1). Dans le cas du grès de Fontainebleau, l'écart entre la vitesse d'impetus VFB et la vitesse de phase VPH à forte pression de confinement est pratiquement négligeable. Par conséquent, pour toute valeur de l'amplitude KSI des fluctuations de vitesses ultrasonores dans la roche, le rapport D/HS entre la distance D de propagation de l'onde sismique et la taille moyenne HS des hétérogénéités est infiniment grand. Par conséquent, la taille moyenne HS des hétérogénéités attendue est très petite. Ceci est en accord avec la taille moyenne des hétérogénéités de cette roche mise en évidence sur les lames minces, à savoir des grains de taille moyenne de l'ordre de 0,25 mm (voir Figure 4 et dernière colonne du Tableau 1). Le dernier cas considéré est le grès à Meule. L'écart relatif (VFB - VPH)/ VPH mesuré vaut environ 1%. Pour des valeurs de l'amplitude KSI des fluctuations de vitesses ultrasonores dans la roche comprise entre 1% et 2%, intervalle des valeurs de KSI considéré pour tous les grès analysés, le rapport D/HS déduit de l'abaque de la Figure 3 est compris entre 14 et 57. La distance de propagation valant 5 cm, la taille moyenne estimée des hétérogénéités est comprise entre 0,9 mm et 4 mm. Ceci est tout à fait en accord avec la taille moyenne des hétérogénéités de densité mesurée par tomographie à rayons X, à savoir 1 mm à 3 mm (voir Figure 4 et dernière colonne du Tableau 1). Roche VFB-VPH KSI D/HS HS HS (mm) (%) (mm) par mesures indépendantes VPH (%) Calcaire 6% 4 à 5% 13 à 21 2,4mm - 3,8mm 2mm - 5mm Euville Calcaire 2 à 7% 4 à 5% 4 à 25 2mm - 12,5mm 4mm - 7mm Brauvilliers Grès 0% 1 à 2% Grand Petit 0,25mm Fontainebleau Grès 1% 1 à 2% 0,9 à 4 1mm - 3mm Meule 14 à 57 Tableau 1 16 REVENDICATIONS1. Méthode pour estimer la taille moyenne des hétérogénéités d'une formation souterraine, S liées à la matrice rocheuse de ladite formation ou bien à la répartition de fluides dans ladite formation, au moyen de mesures d'ondes sismiques, caractérisée en ce que l'on réalise les étapes suivantes : a. on détermine au moins un groupe de points d'intérêt au sein de ladite formation pour lesquels on dispose d'une information sismique associée ; 10 b. on calcule la distance D de propagation des ondes sismiques entre lesdits points de chaque groupe ; c. on détermine la vitesse d'impetus VFB entre lesdits points de chaque groupe ; d. on détermine la vitesse de phase VPH entre lesdits points de chaque groupe ; e. on évalue la taille moyenne HS des hétérogénéités réparties entre lesdits points de 15 chaque groupe, à partir de ladite distance D, de ladite vitesse d'impetus VFB et de ladite vitesse de phase VPH. 2. Méthode selon la revendication 1, dans laquelle on effectue une étape supplémentaire pour estimer l'ordre de grandeur des fluctuations moyennes des vitesses sismiques KSI dans ladite formation. 20 3. Méthode selon la revendication 2, dans laquelle ladite taille moyenne HS des hétérogénéités réparties entre lesdits points de chaque groupe est calculée par la formule suivante : {2 VFB-VPH 1 VPH(KSI)2 HS D(DLY/ avec D correspondant à ladite distance entre lesdits points d'un même groupe, DM à la 25 dimension de l'espace du problème considéré, VFB à la ladite vitesse d'impetus entre lesdits points de chaque groupe, VPH à ladite vitesse de phase entre lesdits points de chaque groupe, et KSI auxdites fluctuations moyennes des vitesses sismiques dans ladite formation. 17 3011939 4. Méthode selon l'une des revendications 2 à 3, dans laquelle lesdites mesures sont réalisées au moyen d'ondes ultrasonores, et en ce que lesdits points d'intérêt d'un même groupe correspondent à la position d'une source et à la position d'un récepteur. 5. Méthode selon l'une des revendications 2 à 3, dans laquelle lesdites mesures sont réalisées par diagraphies soniques, et en ce que lesdits points d'intérêt d'un même groupe correspondent à deux récepteurs distincts. 6. Méthode selon l'une des revendications 2 à 3, dans laquelle lesdites mesures sont réalisées par de la sismique de puits en mode transmission à déport nul, et en ce que lesdits points d'intérêts d'un même groupe correspondent à la position de deux récepteurs distincts. 7. Méthode selon l'une des revendications 2 à 3, dans laquelle lesdites mesures sont réalisées par de la sismique de puits en mode réflexion, et en ce que lesdits points d'intérêt d'un même groupe sont situés sur deux réflecteurs sismiques différents. 8. Méthode selon l'une des revendications 2 à 3, dans laquelle lesdites mesures sont réalisées par de la sismique de puits à puits, et en ce que lesdits points d'intérêt d'un même groupe correspondent à la position d'une source et à la position d'un récepteur. 9. Méthode selon l'une des revendications 2 à 3, dans laquelle lesdites mesures sont réalisées par de la sismique réflexion de surface, et en ce que lesdits points d'intérêt sont situés sur deux réflecteurs sismiques différents. 10. Méthode selon l'une des revendications 2 à 9, dans laquelle lesdites mesures sont réalisées après une substitution de fluide dans ladite formation, et en ce que les fluctuations moyennes des vitesses sismiques KSI correspondent à la variation de vitesse induite par cette substitution de fluide. 11. Méthode selon la revendication 10, dans laquelle le fluide de substitution est du CO2.