FR 2654515 A1 19910517 FR 9013631 A 19901026 La présente invention a trait à une méthode statistique dite de Monte Carlo pour l'estimation d'une variation dans le type d'une roche ou sa texture, c'est-à-dire le changement dans sa lithologie, le long d'une couche donnée ou d'un groupe apparenté de couches à l'intérieur d'une formation géologique sélectionnée, à partir de données sismiques recueillies le long d'un réseau de surveillance coïncidant avec des points de contrôle géologiques tels que des puits de forage. Comme cela est bien connu dans l'art de lletude sismique, on génère une impulsion acoustique sur ou près de la surface du sol. L'impulsion acoustique se propage suivant un front d'onde sphérique dans le sol. Le front d'onde est réfléchi de sur les couches souterraines jusqu'à la surface où les données réfléchies sont détectées par des géophones ou des hydrophones sous forme de signaux électriques numériques ou modulés en amplitude. Les signaux sont emmagasinés sous forme d'enregistrement en fonction du temps pour le traitement ultérieur pour procurer une coupe sismique d'une partie du sol. Dans l'exploration pétrolière, les données sismiques sont utilisées pour évaluer la profondeur, l'altitude, et la lithologie des formations rocheuses entre les puits et pour rechercher les zones probables pour un nouveau forage. L'amplitude d'un front d'onde qui est réfléchi depuis une interface entre deux couches, est déterminée par le coefficient de réflexivité qui est une fonction de la différence d'impédance acoustique entre les deux couches. L'impédance acoustique d'une couche est définie comme la vélocité caractéristique à l'intérieur de cette couche de roche (la vélocité d'intervalle) multipliée par la densité de la roche. L'impédance acoustique est mesurée en unités de longueur (pieds par seconde, par gramme, par centimètre cube. L'amplitude des réflexions sismiques reçues et enregistrées peut être mesurée et corrigée par rapport à leur étalement sphérique, aux artefacts instrumentaux et aux autres variables prévisibles, pour fournir des données "d'amplitude vraie. Les impédances acoustiques des couches respectives peuvent ensuite être calculées à partir des données d'amplitude. Un profil d'impédance acoustique est utilisé pour l'estimation du type et de l'étendue de certaines classes de roches et de leurs caractéristiques pétrophysiques, le long d'un réseau de surveillance sismique tel que celui forme; par des puits. Par exemple, des argiles schisteuses ont tendance en moyenne à être caractérisées par une impédance acoustique un peu plus grande que le grès. Toutefois, les couches de la roche ne sont pas nécessairement homogènes sur leur étendue horizontale et verticale. Des sables hétérogènes se classent dans des argiles schisteuses sablonneuses.Un changement dans la composition, la porosité, la perméabilité, ou le contenu en fluide de la masse rocheuse peut causer une variation concomitante des attributs sismiques tels que l'amplitude vraie, le temps de transit acoustique (vélocité réciproque) ou l'impédance acoustique. Donc, une valeur particulière d'un attribut sismique, associée å un horizon réfléchissant sélectionné, n'est pas le seul diagnostic des caractéristiques petrophysiques qui ont donné lieu à cette réflexion. L'impédance acoustique des sables, argiles schisteuses et autres types ou classes de roches peut être mesurée directement dans les puits de forage par l'usage combiné des enregistrements de gamma graphe, des enregistrements de neutrons et des enregistreurs de la vélocité continue conjointement avec des analyses de carottes. Des e;chantil- lons de roche provenant d'exploitations à ciel ouvert peuvent être étudies de la même façon. Etant donné une étude sismique proche, les caractéristiques pétrophysiques des types de roche peuvent être mis en corrélation directe avec les attributs sismiques. Par exemple, sur la fig. 1, un histogramme de frequence illustre les temps synthétiques de transit sismiques comme fonction des intervalles discrets ou des rangs de perméabilité pour un échantillon de grès de Berea. La dispersion des données dans le temps de transit a produit un recouvrement considérable des données observées. On pourrait s'attendre à ce qu'un histogramme d'autres variables pétrophysiques dans une seule couche telles que la porosité ou la lithologie, puisse illustrer une diffusion similaire sur une base surfacique. Un profil des attributs sismiques est donc ambigu ; le degré d'ambiguïté est proportionnel à l'aire de recouvrement des distributions de fréquence. On ne connaît bien les données de mesure pétrophysiques qu'aux points de contrôle, lesquels sont généralement assez dispersés sur une certaine zone. Des données sismiques peuvent être disponibles à des postes d'échantillonnage rapprochés à l'intérieur de cette zone mais les données sont ambigues quant aux propriétés précises de la roche. Etant donné deux sources de données imparfaites, l'objet de l'invention est de combiner les deux sources pour évaluer statistiquement les propriétés pétrophysiques vraies. L'idée-mère de l'invention est fondée en partie sur les techniques d'enrichissement d'image telles qu'appliquées pour la télédétection, technique décrite par J. BESAG ("On the Statistical Analysis of Dirty Pictures", Journal of the Royal Statistical Society, v. 48, N03,1986). Une méthode apparentée a été décrite par C.L. Farmer, "The Mathematical Generation of Reservoir Geology", présentée en juillet 1989 à IMA/SPE European Conférence on the Mathematics of Oil Recovery. L'invention a trait à une méthode statistique assistée par ordinateur dite de Monte Carlo pour combiner les mesures géologiques discrètes des propriétés de la roche avec les mesures continues des attributs sismiques, et pour la transformation des ces mesures physiques combinées, en une représentation de la meilleure estimation des classes de roche souterraine définies comme des rangs discrets des propriétés pétrophysiques ou des classes lithologiques. selon l'invention, il est propose une méthode pour modéliser les caractéristiques pétrophysiques d'une partie du sous-sol. Il est établi la présence d'une pluralité de classes de roche discrètes dans une couche du sol. Les données statistiques de transition spatiales entre les classes de roche d'une couche du sol sont mesurées entre des points de contrôle géologique discrets espaces. Des nombres entiers ou des teintes distinctes sont attribués à chaque classe de roche à des fins d'affichage et de représentation des propriétés pe;tro- physiques. Un réseau continu de stations sismiques est disposé entre les points de contrôle geolo giques. Des enregistrements d'attributs sismiques sélectionnés sont réalisés à chaque station sismique. Les attributs sont sélectionnes sur l'hypothèse qu'ils sont indicatifs des propriétés petrophysiques étudiées. Les moyennes et variances des attributs sélectionnés sont calculées pour chaque classe de roche ou de propriété petrophysique. Un modèle initial est créé qui se compose d'un réseau de pixels à N dimension où N peut être égal à 2 ou 3. Les dimensions du réseau de pixels sont en rapport avec le nombre des stations sismiques. Il est affecté au hasard à chaque pixel du modèle initial l'une d'au moins deux teintes possibles . Un pixel est choisi au hasard. L'énergie globale du modèle entier est calculée pour chaque teinte possible du pixel choisi. L'énergie globale est la somme de l'énergie probabiliste et de l'énergie antérieure. La teinte affectée au pixel choisi sera la teinte pour laquelle l'énergie totale est minimisée. Le processus est répété pour tous les pixels du modèle. Les teintes de pixel sont ensuite perturbées itérativement une pluralité de fois pour tous les pixels jusqu a ce que l'énergie globale soit minimisée et le modèle stabilisé. Selon un aspect de l'invention, les attributs sismiques pris en compte sont l'impédance acoustique, l'amplitude sismique vraie et le temps de transit sismique. Selon un autre aspect de l'invention, l'énergie probabiliste est une mesure de l'écart entre les attributs sismiques enregistrés et les attributs sismiques moyens calculés à partir des données sismiques pour chaque classe de roche. L'énergie antérieure est une mesure de l'encart entre les données statistiques de transition spatiales calculées à partir des données géologiques et les données statistiques de transition spatiales déterminées à partir du modèle. La présente invention sera mieux comprise par référence à la description détaillée et aux dessins dans lesquels - La fig. 1 est un graphe de distribution de fréquence des temps de transit sismiques synthétiques pour des rangs discrets de la perméabilité sur un échantillon du grès de Berea. - La fig. 2 est un organigramme d'un procédé assisté par ordinateur, conforme à l'invention - la fig. 3 illustre la distribution vraie bi-dimensionnelle des rangs de perméabilité mesurées sur l'échantillon ci-dessus - la fig. 4 est une représentation synthétique d'un certain attribut sismique, exprimée en termes de temps de transit sismique par pixel, dérivé des données mesurées de la perméabilité - la fig. 5 est le modèle initial aléatoire ; - la fig. 6 est une carte recons tituée à partir des données sismiques et géologiques après la première itération ; - la fig. 7 est la reconstitution après la troisième itération ; - la fig. 8 illustre l'apparence du modèle après la cinquième itération ; et - la fig. 9 est le modèle final après onze itérations. A partir de données expe;rimen- tales, telles que celles recueillies à partir d'une étude géologique de reconnaissance ou à partir de puits de forage, il est établi que pour une couche donnée de roche, il existe des variations dans les caractéristiques pétrophysiques ou les catégories de roche à l'intérieur de cette couche. En utilisant ces données, il est proposé un projet de classification pétrophysique statistique qui utilise à la fois les attributs sismiques tels que 1' impédance, les données d'amplitude vraie, ou les temps de transit sismiques et des informations connues relatives à la variation spatiale des propriétés pétrophysiques exprimées en termes de fréquences de transition entre les catégories pétrophysiques comme fonction de distance.Le schéma de classification statistique est basé sur la théorie de Bayes, théorie selon laquelle les probabilités sont associées à des évènements individuels indépendants entre eux et non pas simplement aux fréquences de ces évènements. L'estimation maximale à posteriori de la classe vraie de la roche est déduite ensuite de cette classification ; en cela, elle est la meilleure estimation de la classe de la roche comme provenant de faits observes. Etant donné une zone N-dimensionnelle R, divisée en un réseau de pixels représentant des stations d'observation des données, on représente une image pétrophysique du sous-sol par un vecteur Xtrue = (xltrue,..., xntrue) dans lequel xl représente la classe vraie de la roche au pixel i et n est le nombre total d'éléments d'image (pixels). Les variables Xi peuvent prendre seulement des valeurs discrètes l,2,...,c telles que par exemple 1 pour le sable et 2 pour l'argile schisteuse, bien que des valeurs additionnelles 3 et 4 puissent être introduites comme pour la pierre à chaux ou la dolomite ou, dans le cas d'espèce, des rangs de perméabilité à l'intérieur d'un seul type de la roche. Alternativement, pour des fins de représentation visuelle, on pourrait assigner des niveaux de gris aux pixels comme du noir, du blanc et des nuances de gris ou d'autres teintes du specte à la place des valeurs numériques. Un enregistrement sismique zi est associé à chaque pixel i. Le vecteur des données sismiques pour#l'image entière est défini par z = (z1,..., zn). Le nombre des pixels, n, est rapporté au nombre total des échantillons d'enregistrements sismiques. Il est à noter que les variables sismiques Zi sont continues. La taille physique de chaque pixel est par conséquent rapportée à l'espacement sismique de la station. Le problème est de déduire la distribution pétrophysique spatiale vraie xtrue , habituellement connue seulement aux points de contrôle sus-visés tels que puits de forage, à partir du vecteur sismique Z mesuré. L'estimation xest du vecteur pétrophysique vrai Xtrue est choisie comme celle donnant la probabilité maximale des données sismiques Z. Si on applique la règle de Bayes La probabilité de X étant donné les données sismiques Z doit être maximisée par rapport à X. L'estimation Xest est le maximum à posteriori (MAP) de l'estimation de la valeur vraie Xtrue Supposant la condition Gaussienne d'indpendance des enregistrements d'impédance sismique (ou des enregistrements d'autres attributs), la probabilité conditionnelle p(ZIX) peut être représentée par où E1 est l'énergie probabiliste définie par L'énergie probabiliste peut être considérée comme une mesure de l'écart entre les attributs sismiques enregistrés et les moyennes des attributs sismiques de chaque classe de roche.Dans (3), u(xi) et 02(xi) représentent respectivement la moyenne et la variance d'un attribut tel que l'impédance pour les pixels du type xi, qui peut être égal à 1, 2 ou tout autre nombre entier (noir, blanc ou autre teinte) selon la classe de la roche au pixel i. Les moyennes et variances peuvent être estimées à partir d'un graphe de distribution de fréquence du paramètre étudié , comme le montre la figure l qui est un diagramme de distribution de fréquence des temps de transit sismiques pour différents rangs de perméabilité. Dans (1), p(X) est le modèle de distribution antérieur qui est construit à partir de fréquences expérimentales de transition entre les classes pétrophysiques, en fonction des distances horizontale et verticale, u et v respectivement, qui peuvent par exemple être mesures en termes d'unités de pixels. Supposant un exemple bi-dimensionnel, r et s sont les valeurs des variables de la classe de roche Xi et x3 pour n'importe quelle paire de pixels i et j, séparés (en unités de pixels) par un vecteur de distance (u,v) dans le plan de la coupe sismique. Par la définition ci-dessus, r et s sont des nombres entiers égaux, par exemple 1 ou 2. Pour un vecteur distance donne (u,v), séparant une paire de pixels i et j, la quantité Nrs(u,v) est définie comme le nombre de paires de pixels i-j de limage, de telle sorte que le pixel i appartienne à la classe r et le pixel j appartienne à la classe s. En utilisant la définition cidessus pour Nrs(u,v), la distribution antérieure est donnée par p(x) exp(-Ep), (4) dans laquelle Dans l'équation (5), NtrUers représente les valeurs statistiques de transition prévues pour l'image pétrophysique vraie xtrue, Nrs indique les valeurs statistiques de transition associées à une classe arbitraire image X de la roche, tandis que c est le nombre total de teintes possibles ou de types de roche. La quantité Ep est appelée l'énergie antérieure, c'est-à-dire l'énergie qui est attribuée, comme étant connue, à une classe d'image de la roche, à priori. Ep est une mesure de l'écart entre les valeurs statistiques de transition spatiales dérivées de la géologie connue aux points de contrôle et les valeurs statistiques de transition spatiales du modèle X. Il est évident, à partir des équations (4) et (5), que la qualification basse énergie, ou haute probabilité, est affectée à des images pour lesquelles les valeurs statistiques de transition sont proches de leurs valeurs attendues Ntruers et vice-versa. Le paramètre de variance d2p dans (5) contrôle la grandeur de l'energie Ep. C'est un paramètre d'usage défini qui reflète le degré de precision de la valeur statistique Ntruere. Si d2px O, cela implique une grande confiance dans les données, tandis que si d2p#oo cela implique des données peu fiables. En pratique, les dimensions de 62p sont comparées avec la grandeur des valeurs statistiques de Ntruers. Dans ltéquation (5), la première sommation s'étend à toutes les distances u et v pour lesquelles les valeurs statistiques de transition Ntruere, sont disponibles à partir par exemple du noir au blanc, du blanc au noir, du blanc au blanc, etc.... Il doit être rappelé que, bien que des données sismiques soient habituellement disponibles à des incréments de distance régu- lièrement proches, u et v, les données réelles pour déterminer Ntruers peuvent être éparses et disponibles seulement aux points de contrôle discrets qui peuvent être largement dispersés. A Partir des définitions exponentielles pour p(X) et p(ZIX) dans (2) et (4), et à partir de (1), il est clair que trouver le maximum à priori de la valeur estimée xest des propriétés pétrophysiques vraies xtrue est équivalent à trouver un modèle pétrophysique X qui minimise l'énergie totale ou globale E E(X) = E1(X)#+ E#(X), (6) On arrive à la réduction minimale de la fonction objective E(X) en (6) par un processus itératif assisté par ordinateur montré de façon schématique dans l'organigramme de la fig. 2 qui sera utilisé pour expliquer le meilleur mode operatoire, en relation avec les fig. 3 à 9.Aux fins d'explication, on se référera aux "teintes" plutôt qu'aux nombres entiers, avec l'entier 1 associé à la teinte blanche (si des niveaux de gris sont utilisés) à titre d'exemple, mais pas de façon limitative. En référence maintenant à l'organigramme de la fig. 2, on cre;e d'abord un modèle pe;tro- physique initial N-dimensionnel qui peut être un réseau de pixels teintés au hasard. Dans la fig.2, N est égal à 2 à titre d'exemple non limitatif. Les pixels peuvent prendre n'importe laquelle de plusieurs teintes possibles. A partir des données sismiques et pétrophysiques connues, les moyennes Lz et les variances d21 pour chaque teinte (équation (3)) sont calculées et mémorisées. Similairement, les valeurs statistiques de transition spatiales Ntruers(u,v) sont générées à partir des données géologiques connues pour chaque teinte possible aux points de contrôle géologiques. Un pixel i est choisi au hasard à partir du modèle initial. L'energie probabiliste, E1 est calculée à partir de la moyenne et de la variance pour cette teinte du pixel. Les valeurs statistiques de transition spatiales pour le modèle Nrs(u,v) sont calculées pour chaque valeur u et v et chaque valeur possible de r et s, et comparée avec Ntruers pour calculer l'énergie antérieure Ep. Ep et E1 sont sommées pour déterminer l'énergie globale E(X). La teinte assignée au pixel choisi sera la teinte pour laquelle l'énergie globale est minimisée. On continue le processus pour tous les pixels du modèle. Des cycles d'itération globale additionnels sont exécutés en vue de perturber la teinte de pixel sur tous les pixels jusqu'à ce que l'énergie globale E(X) converge jusqu a un minimum et que le modèle devienne stable. En général moins de quinze cycles d'itération sont nécessaires pour obtenir cette convergence. La méthode de l'invention sera mieux comprise en référence aux figures 3 à 9 qui en illustrent non limitativement une forme préférée d'exécution. La fig. 3 illustre un réseau bidimensionnel de pixels qui représentent les noyaux effectifs prélevés dans un bloc de grès de Berea. La perméabilité de chaque noyau a été mesurée et on lui a assigné une teinte qui correspond à un intervalle discret, ou rang, de permeabilité , la mesure est exprimée en milli-darcy, comme le montre l'échelle dans le coin supérieur droit du réseau. A titre d'exemple, les noyaux respectifs servent de pluralité de points de contrôle géologiques à partir desquels on a calcule; les valeurs statistiques Ntruers. Sur la fig. 4, il est montré des rangs discrets d'un attribut sismique synthétique, à savoir le temps de transit sismique mesuré en ,us/m, pour chaque pixel. Dans cet exemple synthétique, les données du temps de transit forment un réseau continu des stations sismiques entre les points de contrôle géologiques. Les temps de transit sont indicatifs des classes de roche représentées par les rangs de perméabilité. On peut voir sur la fig. 1 qu'il y a un recouvrement considérable des rangs. La fig. 5 est le modèle initial tel que créé selon l'organigramme de la fig. 2. Les pixels sont teintés au hasard. La fig. 6 illustre le modèle après l'achèvement du premier cycle d'itération global. Au cours du processus, les teintes des 1163 pixels ont été remplacées. Après la première itération globale, l'energie globale calculée s'élève à 3305. Le pourcentage des pixels "disclassésQswélève à 57%. La fig. 7 est l'image après 3 itérations globales. La valeur de E est réduite à 1102 après que quelque 267 teintes de pixels aient été échangées. Le taux de "dis-classement"ha été réduit à 36%. Les résultats de la cinquième itération globale sont montrés sur la fig.8. Le modèle initial se stabilise en ce que la valeur de E, 982 n'a pas beaucoup changé par rapport à la valeur première de 1102. Au cours de la onzième itération globale de la fig.9, E a diminué à 950. Grâce à toutes ces opérations pratiques, le modèle initial est maintenant stabilisé ; il est improbable qu'une itération globale additionnelle donne lieu à des améliorations significatives. Après l'itération finale, le taux de déclassementssest de 33%. il doit être compris que l'exemple synthétique pre sente; est quelque peu idéalisé de façon à démontrer plus clairement, le processus global. Dans des conditions ef fectives de terrain, il est improbable que les points de contrôle géologiques soient espacés de façon aussi dense par rapport aux données sismiques, comme cela est illustré ici. Bien que l'on ait décrit et représenté une forme préférez d'application de la méthode de l'invention, il doit être compris que la portée de la présente invention n'est pas limitée à cette forme mais qu'elle s'étend à toute méthode comportant les caractéristiques énoncées plus haut. REVENDICATIONS l.-Methode pour modéliser un profil continu des propriétés pétrophysiques des couches du sou-sol, caractérisée en ce qu'elle comprend les étapes consistant à - établir la présence et mesurer les propriétés pétrophysiques d'une pluralité de classes de roche géologique dans une couche du sol, lesdites classes de roche incluant les rangs discrets de propri étés pétrophysiques, - assigner l'une d'au moins deux teintes possibles à chaque dite classe de roche, - mesurer les valeurs statistiques de transition spatiales entre lesdites classes de roche dans ladite couche du sous-sol entre une pluralité de points de contrôle géologiques espacés, - disposer un réseau continu de stations sismiques entre ladite pluralité des points de contrôle géologiques espacés, et enregistrer des attributs sismiques indicatifs des propriétés pétrophy siques desdites classes de roche, - calculer les moyennes et les variances desdits attributs sismiques pour chaque dite classe de roche, - créer un modèle initial des couches du sous-sol consistant en un réseau N-dimension nel de pixels et assigner au hasard une des deux au moins couleurs possibles à chacun desdits pixels dudit modèle initial, - choisir un pixel au hasard à partir dudit modèle, - calculer pour chaque teinte possible dudit pixel l'energie probable et l'énergie antérieure, - sommer l'énergie probable et l'énergie antérieure pour évaluer ltenergie globale, - choisir celle desdites au moins deux couleurs possibles de pixel qui minimise l'énergie globale, répéter les étapes de choix, calcul, sommation, et sélection pour tous les autres pixels du réseau qui forment le modèle initial ce par quoi est réalisé un premier cycle global d'itération, - exécuter des cycles additionnels d'itération globale en perturbant la teinte de pixel sur tous les pixels jusqu a ce que l'énergie globale converge en un minimum pour fournir un modèle final, et - afficher ledit modèle. 2.-Méthode selon la revendication 1, caractérisée en ce que ledit attribut sismique est l'impedance sismique. 3.-Méthode selon la revendication 1, caractérisée en ce que ledit attribut sismique est l'amplitude vraie. 4.-Méthode selon la revendication 1, caractérisée en ce que l'attribut sismique est le temps de transit sismique. 5.-Méthode selon la revendication 1, caractérisée en ce que ladite énergie proba ble est une mesure de l'écart entre les attributs sis miques enregistrés et les moyennes des attributs sis miques calculées pour chaque classe de roche à partir desdites données sismiques. 6.-Méthode selon la revendication 1, caractérisée en ce que l'énergie antérieure est une mesure de l'encart entre les valeurs statistiques de transition spatiales calculées à partir des données géologiques et les valeurs statistiques de transition spatiales déterminées à partir du modèle. 7.-Méthode selon la revendication 1, caractérisée en ce que la taille physique de chacun desdits pixels est rapportée à l'espacement des dites stations sismiques.