FR 2765692 A1 19990108 FR 9708512 A 19970704 La présente invention concerne une méthode pour modéliser en 3D l'impédance d'un milieu hétérogène à partir d'ondes réfléchies par les discontinuités du milieu, en réponse à des ondes émises dans le milieu. Les ondes utilisées peuvent être des ondes élastiques émises par une source sismique, ou des ondes électromagnétiques émises par un radar par exemple, etc. La méthode selon l'invention trouve des applications notamment pour l'étude du sous-sol dans le but par exemple de rechercher des gisements d'hydrocarbures. On connaît différentes approches pour convertir des données sismiques en impédances acoustiques, qui diffèrent par le mode de paramétrisation du modèle io d'impédance, par la définition des fonctions objectif et par la faculté d'introduire de l'information a priori. La paramétrisation peut être concise, l'impédance étant constante ou linéaire suivant une direction verticale. dans un certain nombre de couches de géométrie inconnue. Les paramètres inconnus sont ici non seulement les impédances acoustiques mais aussi les emplacements des interfaces ainsi que leur nombre. Une méthode pour modéliser un paramètre physique tel que l'impédance acoustique d'un milieu hétérogène est décrit par exemple dans le brevet EP 0.354.1 12 du demandeur. La paramétrisation peut aussi être plus dense, le champ d'impédance étant discrétisé par utilisation d'une grille 3D. La fonction objectif sismique consiste en une norme d'un type connu L' par exemple ou L2, de l'écart entre les données synthétiques et les données observées. Les données synthétiques résultent de la convolution d'une ondelette optimale donnée, calibrée aux emplacements des puits avec la série des coefficients de réflexion calculés à partir de la diagraphie d'impédance. La quantité concernée dans la fonction objectif géologique peut être l'impédance elle-même ou certaines de ses dérivées spécialement dans la direction des couches. Les types de contraintes appliquées à ces quantités peuvent être une norme de type L' ou le carré d'une norme L2, de l'écart entre le modèle optimal et le modèle a priori, aussi bien que des égalités ou inégalités entre la quantité recherchée et des valeurs numériques connues. Les méthodes d'optimisation les plus pratiquées sont les méthodes de gradient conjugé ou des techniques dites de "recuit simulé". Par Brac J. et al; "Inversion with A Priori Information: an Approach to Integrated Stratigraphic Interpretation"; in Sheriff R.E.Ed.. Reservoir Geophysics, Soc. Expl. Geophys., Investigations in Geophysics 7, il est connu de choisir comme fonction objectif géologique, le carré de la norme L2 de l'écart d'impédance et de l'écart du gradient d'impédance, après projection sur la direction du pendage vu dans le plan vertical des lignes sismiques. De cette façon, l'inversion d'un bloc 3D consiste I0 en plusieurs inversions 2D multi-traces appliquées ligne sismique après ligne sismique, même si le modèle a priori est élaboré en 3D. La fiabilité et la précision des modèles de réservoir dépendent du degré d'intégration entre les données sismiques et géologiques. Un modèle de réservoir réaliste doit combiner à un stade aussi précoce que possible du processus de réalisation, des données sismiques après sommation (post-stack) relativement exactes, des données diagraphiques de puits peu nombreuses latéralement mais précises verticalement, ainsi que des données interprétatives fournies par une étude géologique régionale. L'inversion stratigraphique après sommation (poststack) est un pas important vers l'intégation des géosciences. Dans la suite de la description, le terme général volume "feuilleté" ou volume pourvu d'un "feuilletage", utilisé pour définir le milieu à modéliser, est pris dans son sens géométrique. Il s'agit d'un volume constitué d'un ensemble de feuilles ou surfaces connexes et disjointes dont la réunion constitue le volume. La méthode selon l'invention a pour objet de déterminer un modèle optimal en trois dimensions (3D) de l'impédance présentée par un milieu hétérogènes à partir de données enregistrées correspondant à des ondes émanant du milieu hétérogène. en réponse à des ondes qui y sont transmises et d'un modèle a priori en 3D de ce milieu hétérogène, connaissant la position de différentes discontinuités d'impédance du milieu, et après interprétation stratigraphique. La méthode est caractérisée dans sa définition la plus générale, en ce qu'elle comporte en combinaison: - la construction d'un modèle géométrique en 3D du milieu hétérogène, comprenant plusieurs volumes feuilletés (pourvus d'un feuilletage) avec définition (positionnement) de leurs différentes feuilles; - la construction d'un modèle a priori d'impédance à partir de ce modèle géométrique et d'une pluralité de mesures d'impédances faites à différentes profondeurs du milieu hétérogène; - la sélection d'un modèle de covariance le long des feuilles des dits feuilletages; et - la formation du modèle optimal par inversion des données enregistrées en utilisant le modèle d'impédance a priori. La méthode peut être utilisée notamment pour obtenir un modèle optimal en 3D de l'impédance d'une zone du sous-sol, les différents volumes feuilletés étant dans ce cas des unités sédimentaires, les mesures d'impédance étant obtenues à différentes profondeurs d'au moins un puits au travers de la zone, et les feuilles des différents volumes feuilletés étant constitués par des isochrones de dépôt dans chacune de ces unités sédimentaires. Le modèle d'impédance a priori peut être formé par exemple par interpolation simple ou par krigeage le long d'isochrones de dépôt, des valeurs d'impédance connues au(x) puits. Les données enregistrées peuvent être obtenues par exemple à partir d'ondes élastiques ou d'ondes électromagnétiques. Suivant un mode préféré de mise en oeuvre, on sélectionne de préférence un modèle de covariance exponentiel à 2D isotrope ou anisotrope. La formation du modèle optimal peut être obtenu par exemple par détermination sur l'ensemble du modèle 3D du minimum d'une fonction objectif globale comprenant un terme relatif aux données enregistrées constitué par le carré de la norme L2 de l'écart entre des données synthétiques et les données enregistrées et un terme d'impédance constitué d'une part par le carré de la norme L2 de l'écart d'impédance et d'autre part, du carré de la norme L2 du gradient de l'écart d'impédance, après projection de ce gradient sur le plan tangent à une feuille locale des volumes feuilletés (un plan de pendage dans le cas d'une zone du sous-sol par exemple). La méthode selon l'invention constitue une approche flexible pour estimer les impédances dans la bande des fréquences utilisées, contraintes par les données 3D enregistrées, les logs d'impédance (mesurés dans les puits par exemple), et les informations connues du milieu (les données géologiques tant stratigraphiques que structurales par exemple). D'autres caractéristiques et avantages de la méthode selon l'invention, apparaîtront à la lecture de la description ci-après d'un exemple non limitatif de mise en oeuvre, en se référant aux dessins annexés o: - les Fig. lA, lB montrent une section d'impédance i transversale intersectant un puits après respectivement une inversion 2D multiple et une inversion 3D; - la Fig.2 montre les incohérences entre lignes voisines d'une section d'impédance i, qui peuvent exister après une inversion multi-2D; - la Fig.3 illustre bien en quoi l'inversion 3D contribue à atténuer les incohérences visibles sur la Fig.2; - la Fig.4 montre différentes courbes (a), (b), (c) de variation du coefficient de corrélation entre les reflectivités associées à deux traces d'impédance adjacentes, correspondant respectivement aux cas des Fig.2, Fig.3 et au cas de la Fig.5 ci- après, les corrélations faibles vers la gauche de la figure étant dues à la présence de failles; - la Fig.5 montre une coupe horizontale à la base du bloc considéré après une inversion 3D o les paramètres d'incertitude ont été choisis en faveur de l'information sismique; - les Fig.6A, 6B montrent à la fois les traces d'impédance inversées et les logs (ou diagraphies) d'impédance obtenues dans des puits; et - la Fig.7 montre un organigramme illustratif de la méthode selon l'invention. La méthode comporte on l'a vu dans un premier temps la construction d'un modèle géométrique en 3D comprenant plusieurs volumes feuilletés qui. dans l'exemple de mise en oeuvre décrit, sont des unités sédimentaires dont les feuilles sont des isochrones de dépôt. On construit ensuite un modèle a priori d'impédance à partir de ce modèle géométrique et d'un ensemble de mesures d'impédances faites par exemple le long d'un ou plusieurs puits au travers de la zone explorée. La mise en ceuvre de ces deux étapes est décrite par exemple dans: - Joint EAEG/EAPG/EAGO Multidisciplinary Workshop on "Developing New Reservoirs in Europe", 1994, St Petersburg; Extended Abstracts Book. Une étape essentielle de la méthode selon l'invention consiste à minimiser une fonction objectif globale calculée sur l'ensemble du modèle 3D au moyen d'une méthode de gradient conjugué. Cette fonction objectif globale est constituée de deux termes. Le premier, qui est le terme sismique, est constitué par le carré de la norme L2 de l'écart entre les données synthétiques et les données observées. Le second, qui est le terme d'impédance, est constitué d'une part par le carré de la norme L2 de l'écart d'impédance et d'autre part par le carré de la norme L2 du gradient de l'écart d'impédance, après projection sur le plan du pendage. Cette méthode d'inversion stratigraphique 3D donne, en une seule passe, un bloc d'impédance 3D optimal compatible à la fois avec le bloc sismique 3D et le modèle 3D a priori. Inversion stratigraphique 3D complète La méthode d'inversion 3D complète est basée sur les deux hypothèses pratiques suivantes couramment utilisées dans le cadre d'études relatives aux réservoirs, i.e.: - la fonction qui, aux traces d'impédance dans le domaine temps, associe les traces sismiques, est calculée par convolution d'une ondelette donnée avec les coefficients de réflexion de chaque trace; - la structure géologique comporte des pendages modérés, des failles de faible rejet et de faibles variations latérales de vitesse. - On admet en outre que lrondelette sismique est constante sur l'ensemble du bloc 3D considéré et que le bruit altérant les traces sismiques est un bruit blanc non corrélé d'une trace à l'autre. - On considère également que les incertitudes relatives au modèle d'impédance a priori sont modélisées au moyen d'une covariance exponentielle 2D. De nombreux modèles de covariance peuvent être utilisés. Certains sont définis simplement par une variance et une longueur de corrélation, mais le choix d'un opérateur de covariance exponentielle est particulièrement intéressant dans la mesure o la norme associée à son inverse est connue et facile à calculer. La fonction objectif sismique La fonction objectif sismique mesure l'écart entre chaque trace sismique et la trace synthétique liée à cette dernière Qs - 2s 2 E Y R*W)ijk ijk) Q Xj,(Rw) s) (1) 2 S i. k Dans cette équation, as = SRMS/p est l'incertitude relative à l'amplitude sismique, SRMS étant la valeur efficace des données sismiques et p le rapport signal- bruit. R désigne la trace du coefficient de réflexion et W lrondelette obtenue après calibrage des données de puits et des données sismiques. Sbks est le i-ième c ijk échantillon de la j-ième trace de la k-ième section ou ligne du bloc sismique 3D. La fonction objectif d'impédance Dans chaque unité sédimentaire, on considère que la fonction de covariance 2D le long de chaque surface de corrélation est la suivante: cov(oZ(x),ÈZ(y)) = C(x, y) = -2 exp-- Y1 (2) i o $Z(x) est l'écart d'impédance par rapport au modèle a priori à l'emplacement x et âZ(y), l'écart correspondant en y, G2 étant la variance et X. la longueur de corrélation. L'opérateur de covariance étant: C: ---> (p = C() = f C(x, y)(y) dS(y) (3) Comme le savent les spécialistes, l'inverse de cet opérateur est H: (p -> p = H(q) = 27r j-- AC - (4) et la norme associée est: y c-' = (x) () dS(x) (5) -1 27['24iy dS+ VJ *dS (6) les symboles A, V désignant respectivement le laplacien et le gradient. Alors que l'inversion 2D implique comme l'on sait un modèle de covariance ID le long des lignes de corrélation, l'inversion 3D implique un modèle de covariance 2D le long des surfaces de corrélation (les isochrones de dépôt) déduites de l'interprétation stratigraphique et structurale. La discrétisation de I(PI1 donne la fonction objectif géologique. Le contenu fréquentiel des impédances inversées Les fréquences contenues dans le modèle a priori qui sont supérieures à celles fournies par la bande passante sismique ne sont pas contrôlées dans le processus d'inversion. C'est la raison pour laquelle elles sont généralement rejetées. Les informations basse fréquence après inversion ne proviennent pas des données sismiques mais du modèle d'impédance a priori. Dans de nombreux cas, ces informations basse fréquence provenant des puits présentent une importance cruciale en matière de prédiction de propriétés quantitatives des réservoirs à partir d'impédances inversées. Ainsi, des informations basse fréquence supplémentaires telles que les vitesses sismiques pourraient être utilisées lors de la construction du modèle d'impédance a priori. Choix des paramètres d'incertitude La variance c2 mesure la variabilité de l'écart d'impédance par rapport au modèle a priori, le long des isochrones de dépôt. La longueur de corrélation X contrôle la continuité latérale de cet écart d'impédance le long des surfaces stratigraphiques. En pratique, on peut évaluer une valeur pour chaque unité géologique soit à partir d'informations géologiques antérieures, soit à partir d'analyses variographiques de valeurs d'impédance aux puits, le long d'isochrones de dépôt. Leurs valeurs peuvent aussi être choisies selon la confiance accordée par l'interprétateur au modèle a priori (nombre de puits, environnement géologique connu aux puits). Résultats numériques L'efficacité de la méthode selon l'invention a été comparée à celle qu'apporte la méthode connue d'inversion 2D multitrace définie dans le brevet EP 0.354.112 précité, sur la base de données d'un bloc sismique 3D acquis en mer du Nord. Ce bloc s'étend sur 8 km2 entre 2,0 et 2,8 s. Le maillage est constitué de 35 lignes avec 130 traces et 201 échantillons temps, ce qui correspond à environ 900 000 noeuds. Les pas d'échantillonnage sont de 75 m entre les lignes, 25 m entre les traces et 4 ms en temps. Neuf diagraphies d'impédance ont été choisis pour la construction du modèle a priori. L'inversion 3D comparée à l'inversion 2D multitrace En utilisant le même modèle 3D a priori, on a réalisé à la fois des inversions multitrace 2D et des inversions 3D sur ce jeu de données. Les Fig. 1 A, lB, relatives à la coupe transversale traversant le puits (trait continu blanc sur la Fig.2), et les Fig.2 et 3, relatives à la coupe effectuée le long d'un horizon sismique, montrent l'amélioration de la continuité le long des surfaces stratigraphiques apportée par l'inversion 3D. F représente ici une faille. Cette observation est confirmée par la Fig.4, o les courbes (a), (b), ( c) représentent le coefficient de corrélation entre les réflectivités associées aux traces d'impédance optimale des lignes 8 (trait pointillé blanc sur la Fig.2) et 9 et de même abscisse x, respectivement après une inversion multi-2D, après une inversion 3D et après une inversion 3D, les données sismiques étant favorisées. Pour les deux inversions (a) et (b), les paramètres sont les suivants: as = 30 %, cy = 250 g.cm'3. m.s'1 et X = 500 m, ce qui favorise les informations géologiques. L'effet des paramètres d'incertitude Une autre inversion 3D a été réalisée avec des paramètres d'incertitude choisis en faveur des données sismiques au lieu du modèle d'impédance a priori (as = 30 %, (y = 500 g.cm-3.m.s-1 et X. = 100 m). Les données sismiques présentent une variation latérale beaucoup plus importante que celle du modèle a priori le long des horizons. De ce fait, en raison de l'influence accrue des données sismiques, le résultat présenté sur la Fig. 5 est plus "rugueux" (bruité latéralement) que celui des Fig. 2 et 3. Avec des paramètres favorisant les données sismiques, la continuité latérale de l'impédance optimale est faible. Avec des paramètres favorisant les informations géologiques, la continuité latérale est augmentée seulement le long de l'axe transversal avec une inversion multitrace 2D et elle est augmentée à la fois le long de l'axe transversal et de l'axe longitudinal avec une inversion 3D. Pour les résultats obtenus à la fois avec l'inversion 2D multitrace et l'inversion 3D, comparaison a été faite entre la trace d'impédance à l'emplacement du puits et la trace d'impédance déduite de diagraphies de puits. On peut noter que ce puits a été utilisé pour construire le modèle a priori. La Fig.6 montre que la trace d'impédance obtenue après inversion 3D est légèrement meilleure que celle obtenue après inversion 2D multitrace. Sur cette figure, les symboles représentent respectivement Tc, le toit du crétacé, - Bc, la base du crétacé, - Bb. la base du Brent, - a), la diagraphie d'impédance, - b), l'impédance après inversion multi-2D, et - c), I'impédance après inversion 3D. Pour réaliser l'étape d'inversion des données enregistrées en utilisant le modèle d'impédance a priori, on peut utiliser un logiciel d'inversion connu et par exemple le logiciel INTERWELL d'inversion stratigraphique 2D ou 3D, comme le montre l'organigramme de la Fig.7. La méthode selon l'invention a été décrite dans une application à la modélisation de l'impédance du sous-sol. Il est bien évident qu'elle peut s'appliquer à tout autre paramètre physique sur les données enregistrées. D'une façon plus générale, on ne sortirait pas du cadre de l'invention en appliquant la méthode selon l'invention à la réalisation de tomographies dans des domaines très différents notamment dans le domaine de l'imagerie médicale ou du contrôle non destructif. Il REVENDICATIONS 1) Méthode pour obtenir un modèle optimal en trois dimensions (3D) de l'impédance présentée par un milieu hétérogène, à partir de données enregistrées correspondant à des ondes émanant du milieu hétérogène en réponse à des ondes qui y sont transmises et d'un modèle a priori en 3D du milieu hétérogène. connaissant la position de différentes discontinuités d'impédance du volume, et après interprétation stratigraphique, caractérisée en ce qu'elle comporte - la construction d'un modèle géométrique en 3D comprenant plusieurs volumes feuilletés avec définition des différentes feuilles de chacun d'eux, - la construction d'un modèle a priori d'impédance à partir de ce modèle géométrique et d'une pluralité de mesures d'impédance faites à différentes profondeurs du milieu hétérogène; - la sélection d'un modèle de covariance le long des feuilles des volumes feuilletés; et - la formation du modèle optimal par inversion des données enregistrées en utilisant le modèle d'impédance a priori. 2) Méthode selon la revendication 1, caractérisé en ce que le milieu hétérogène est une zone du sous-sol, les différents volumes feuilletés sont des unités sédimentaires, les mesures d'impédance sont obtenues à différentes profondeurs d'un puits au travers de la zone, et les feuilles des différents volumes feuilletés sont des isochrones de dépôt dans chacune de ces unités. 3) Méthode selon l'une des revendications 1 ou 2, caractérisée en ce que les ondes émises dans le milieu sont des ondes élastiques. 4) Méthode selon l'une des revendications 1 ou 2, caractérisée en ce que les ondes émises dans le milieu sont des ondes électromagnétiques. ) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on sélectionne un modèle de covariance exponentiel isotrope ou anisotrope. 6) Méthode selon l'une des revendications précédentes, caractérisée en ce que la formation du modèle optimal comporte la détermination sur l'ensemble du modèle 3D du minimum d'une fonction objectif globale comprenant un terme relatif aux données enregistrées constitué par le carré de la norme L2 de l'écart entre des données synthétiques et les données enregistrées, et un terme d'impédance constitué d'une part par le carré de la norme L2 de l'écart d'impédance. et d'autre part. du carré de la norme L2 du gradient de l'écart d'impédance, après projection de ce gradient sur le plan tangent à la feuille locale des volumes feuilletés. 7) Méthode selon les revendications 2 et 6, caractérisé en ce que le plan tangent à une feuille locale des volumes feuilletés, est un plan de pendage. 8) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on sélectionne un modèle de covariance exponentielle à 2D.