FR 2851670 A1 20040827 FR 0302199 A 20030221 La présente invention concerne une méthode pour former plus rapidement un modèle numérique stochastique de type gaussien ou apparenté, représentatif de la distribution spatiale d'une grandeur physique (comme la perméabilité par exemple) dans un milieu hétérogène poreux (tel qu'un gisement d'hydrocarbures par exemple) calé par rapport à des données dites statiques et dynamiques incertaines. Les données dites statiques, correspondent à des observations sur la grandeur physique étudiée elle-même. Lorsqu'elle est certaine, la donnée statique est une valeur exacte. Dans le cas contraire, elle est définie par une loi de probabilité. Les données dites dynamiques, sont caractéristiques du déplacement des fluides dans le milieu : ce sont par exemple des données de production (pressions obtenues à partir d'essais de puits, débits, etc.). La méthode selon l'invention trouve des applications dans la modélisation de zones souterraines où il s'agit de générer des représentations montrant comment est distribuée une certaine grandeur physique dans une zone du sous-sol (la perméabilité, la porosité, les faciès notamment), compatibles au mieux avec des données observées ou mesurées, dans le but par exemple d'en favoriser l'exploitation. ETAT DE LA TECHNIQUE L'optimisation dans un contexte stochastique consiste à déterminer des réalisations d'un modèle stochastique qui satisfont un ensemble de données observées sur le terrain, dites statiques ou dynamiques suivant leurs natures. En ingénierie de réservoir, les réalisations à identifier correspondent à des représentations, dans le champ réservoir, de la distribution de propriétés de transport telles que la perméabilité ou porosité ou encore de la distribution des faciès, chaque faciès correspondant à une famille de propriétés de transport. Ces réalisations forment des modèles numériques de réservoir. Les données statiques disponibles sont, par exemple, des observations ponctuelles de perméabilité, porosité ou faciès et un modèle de variabilité spatiale déterminé selon des mesures ponctuelles. Les données ponctuelles peuvent être définies par des lois de probabilité plutôt que des valeurs exactes. Par exemple, en un point donné, une valeur de porosité peut être caractérisée par une loi de probabilité normale de moyenne 0. 20 et de variance 0. 03. Les données dynamiques sont des données directement liées aux écoulements des fluides dans un réservoir souterrain, c'est à dire des pressions, des temps de percée, des débits, etc. Ces dernières sont souvent non linéairement reliées aux propriétés physiques à modéliser. Une réalisation tirée au hasard n'est pas en général en adéquation avec l'ensemble des données collectées. Intégration des données statiques La cohérence vis à vis des données statiques est intégrée dans le modèle à partir de techniques de krigeage - Journel, A. G., and Huijbregts, C.J., "Mining geostatistics", Académie Press, San Diego, CA, 1978. L'approche générale consiste à générer une réalisation non conditionnelle et à la corriger de sorte qu'elle honore les observations ponctuelles et la structure spatiale. Dans ce cadre, les observations ponctuelles sont supposées exactes (par exemple, en un point donné, la perméabilité vaut 150mD). Soit une réalisation y gaussienne et non conditionnelle du modèle stochastique Y. La réalisation corrigée s'obtient comme suit : y c (x) = y dK (x) + [y(x) - y K (x)] yc est la réalisation gaussienne corrigée et donc conditionnelle. ydK et yK sont obtenues par krigeage à partir respectivement des observations ponctuelles et des valeurs de y simulées aux points d'observation. Lorsque que les observations ponctuelles ne correspondent pas à des valeurs exactes, mais sont définies par des lois de probabilité, on peut utiliser soit une approche bayésienne qui s'avère extrêmement coûteuse en temps calcul, soit une technique de krigeage. Dans ce dernier cas, cas sur lequel nous nous concentrons ici, une étape préliminaire s'impose afin de transformer la loi de probabilité en valeur ponctuelle. Cette transformation ne peut être quelconque : les valeurs ponctuelles en résultant doivent vérifier la structure spatiale et les lois de probabilité dont elles sont issues. Clairement, la relation entre lois de probabilité et valeurs ponctuelles n'est pas unique. Des méthodes de transformation itératives ont été proposées pour les réalisations en faciès, où la loi de probabilité est une loi uniforme, par: - Freulon, X., and Fouquet, C. de, "Conditioning a Gaussian model with inequalities", in Geostatistics Troia'92, A. Soares, ed., Kluwer Academic Publishers, Dordrecht, The Netherlands, 1993. - Le Ravalec-Dupin, M., "Conditioning truncated Gaussian realizations to static data", 8th Ann. Conf. Int. Ass. Math. Geol., Berlin, Germany, 15-20 Sept, 2002. Ces techniques convergent vers la fonction de densité de probabilité souhaitée, mais très lentement. En pratique, on ne cherche pas à échantillonner complètement cette fonction. On se contente de déterminer un jeu unique de valeurs ponctuelles respectant la structure spatiale et les lois de probabilités voulues et on l'utilise pour obtenir la réalisation corrigée yc. C'est à ce jeu unique, et donc invariable, que se réfère constamment le processus d'optimisation pour élaborer des réalisations contraintes simultanément par les données statiques et dynamiques. Ce point est fondamental dans le processus d'optimisation. Intégration des données dynamiques La cohérence vis à vis des données dynamiques est intégrée dans le modèle par le biais d'une procédure inverse: - Tarantola, A., "Inverse problem theory - Methods for data fitting and model parameter estimation", Elsevier Science Publishers, 1987. Récemment, une technique de paramétrage géostatistique qui simplifie le problème inverse, a été introduite pour contraindre, par déformation graduelle, les réalisations stochastiques à des données dont elles dépendent de manière non linéaire. Elle a fait l'objet des brevets FR 2. 780.798 et FR2. 795.841 du demandeur, et des publications suivantes, notamment : - Hu, L. Y., 2000, Graduai déformation and iterative calibration of Gaussian-related stochastic models : Geology, Vol.32, No.l. - Le Ravalec, M. et al. 2000, The FFT moving average (FFT-MA) generator: An efficient numerical method for generating and conditioning Gaussian simulations : Math. Geology, Vol.32, No. 6. - Hu, L. Y., Blanc, G. and Noetinger, B. (2001): Graduai déformation and iterative calibration of sequential stochastic simulations. Math. Geology, Vol. 33, No.4. Cette méthode a été appliquée avec succès à divers cas notamment à partir de données provenant de champs d'exploitation pétrolière, comme décrit dans les documents suivants : - Roggero, F. et al. 1998, Graduai déformation of continuous geostatistical models for history matching, Paper SPE 49004 : SPE Annual Technical Conference and Exhibition, New Orleans. - Hu, L. Y. et al. 1998, Constraining a réservoir facies model to dynamic data using a graduai déformation method, Paper B-01: Proc. 6th European Conference on Mathematics of Oil Recovery (ECMOR VI), 8-11 September 1998, Peebles, Scotland. Supposons que fobs = (f1obs, f2obs, .., fMobs) représente l'ensemble des données dynamiques recueillies sur le terrain et f = (f1, f2,.....,fM) les réponses correspondantes simulées pour une réalisation yc déjà contrainte par les données statiques comme explicité plus haut. En général, les réponses f (f1, f2,...,fM) sont obtenues en résolvant numériquement le problème direct. Ainsi, si yo représente un champ de perméabilité, les données f peuvent être des mesures de pression. Dans ce cas, elles sont simulées à partir d'un simulateur d'écoulement. L'objectif d'une optimisation stochastique est de produire des réalisations de Y qui réduisent les différences entre les données observées et les réponses correspondantes simulées numériquement. Ces différences sont mesurées par la fonction objectif suivante : Les coefficients #m sont des poids attribués aux données fmobs. Les fm sont des fonctions de la réalisation yo discrétisée sur un grand nombre de mailles. En ce sens, la minimisation de la fonction objectif est un problème à plusieurs variables. Soit N le nombre de mailles formant la réalisation yo. N est souvent très grand (104 # 107). Il est donc très difficile de mener une optimisation directement par rapport aux composantes de yo. De plus, la réalisation yo, même modifiée, doit rester une réalisation de Y. Le paramétrage par déformation graduelle permet de s'affranchir de ces difficultés. La technique de déformation graduelle permet de construire une chaîne continue de réalisations en combinant une réalisation initiale yo de Y avec une autre réalisation u1, dite complémentaire, de Y, u1 étant indépendante de yo. Les coefficients de combinaison sont par exemple cos(t) et sin(t) et la réalisation combinée vérifie la relation : y(t)= y0 cos t + u1 sin t oùt est le paramètre de déformation. Dès lors que la chaîne est élaborée, on peut l'explorer en variant le paramètre de déformation t et tenter d'identifier parmi toutes les réalisations de cette chaîne celle qui, après intégration des données statiques par krigeage, minimise la fonction objectif. Cette minimisation se fait par rapport à t. Le paramétrage suivant la méthode de déformation graduelle permet de réduire le nombre de dimensions du problème de N à 1, où N est le nombre de valeurs constituant le champ à contraindre. De plus, la somme des coefficients de combinaison au carré étant 1, la réalisation optimisée est encore une réalisation de Y: elle suit le même modèle de variabilité spatiale que toutes les réalisations de Y. Cependant, si l'on restreint l'exploration de l'espace des réalisations à une unique chaîne, on limite sévèrement nos possibilités de diminuer suffisamment la fonction objectif. Il faut donc répéter la procédure décrite ci-dessus, mais avec de nouvelles chaînes de réalisations. Ces chaînes de réalisations sont construites successivement en combinant une réalisation initiale qui est ici la réalisation optimale déterminée à l'itération précédente, avec une réalisation complémentaire de Y, à chaque fois tirée au hasard. Ainsi, à l'itération l, la chaîne continue de réalisation s'écrit : yl(t) = yl-1 cos t + ui sin t . yl-1 est la réalisation optimale définie à l'itération l-1 et les ul sont des réalisations indépendantes de Y. En minimisant la fonction objectif par rapport à t, on améliore, ou au moins préserve, le calage des données à chaque fois qu'une nouvelle chaîne de réalisations est explorée. Cette procédure de recherche de minimum itérative est poursuivie tant que le calage des données n'est pas satisfaisant. A ce jour, les valeurs ponctuelles intégrées par krigeage dans le modèle pour rendre compte des données statiques restent constantes tout au long du processus d'optimisation, même lorsqu'elles correspondent à des valeurs incertaines. Dans ce dernier cas, une telle hypothèse peut freiner significativement le processus d'optimisation et empêcher la minimisation de la fonction objectif. La méthode selon l'invention La méthode selon l'invention a pour objet de former un modèle numérique stochastique de type gaussien ou apparenté, représentatif de la distribution d'une grandeur physique dans un milieu hétérogène poreux (réservoirs pétroliers, aquifères, etc. ), qui soit ajusté par rapport à des données dynamiques, caractéristiques du déplacement des fluides dans le milieu, et des données statiques locales (par exemple, des valeurs de porosité définies par des lois de probabilité normales ou encore des valeurs précisant la nature des faciès observés définies par des lois de probabilité uniformes) mesurées (par diagraphies par exemple) en un certain nombre de points de mesure le long de puits au travers du milieu (puits de production, d'injection ou d'observation par exemple), et présentant une certaine marge d'incertitude. Elle comporte une optimisation du modèle par le biais d'un processus itératif de déformation où l'on forme à chaque itération une réalisation combinée obtenue par combinaison linéaire d'une part d'une réalisation initiale, représentant au moins une partie du milieu, et d'au moins une deuxième réalisation indépendante du même modèle stochastique, et une minimisation d'une fonction objectif mesurant l'écart entre des données dynamiques réelles et les données dynamiques simulées au moyen d'un simulateur d'écoulement, pour la réalisation combinée, par ajustement des coefficients de combinaison, le processus d' ajustement itératif étant poursuivi jusqu'à l'obtention d'une réalisation optimale du modèle stochastique. La méthode se distingue pour l'essentiel en ce que : a) on transforme les données statiques locales en pseudo-données ponctuelles en accord avec des lois de probabilité et un modèle de variabilité spatiale ; et b) on ajuste les pseudo-données en respectant les lois de probabilités dont elles sont issues, ainsi que le modèle de variabilité spatiale, par le biais d'un processus itératif où l'on combine un premier bruit blanc gaussien associé au jeu de pseudo-données à un deuxième bruit blanc gaussien. Suivant un mode préféré, on ajuste le modèle par un processus de déformation graduelle en imposant que la somme des carrés des coefficients de combinaison entre les réalisations soit égale à 1 et on ajuste alors les pseudo-données par un processus itératif où la somme des carrés des coefficients de la dite combinaison vaut également 1. Les lois de probabilités sont par exemple des lois normales ou des lois uniformes. On peut conduire l'ajustement itératif à partir de deux paramètres de déformation par exemple, avec un premier paramètre qui contrôle la combinaison entre la réalisation initiale et la deuxième réalisation, et un deuxième paramètre qui contrôle la combinaison entre le bruit blanc gaussien initial et le deuxième bruit blanc gaussien. Il est possible également de conduire l'optimisation à partir d'un unique paramètre quand les coefficients de combinaison sont identiques pour la combinaison de réalisations et la combinaison de bruits blancs gaussiens. Suivant un autre mode de mise en #uvre, l'optimisation peut également être conduite par une méthode de points pilotes. En d'autres termes, la méthode comprend une nouvelle technique de transformation qui associe à la loi de probabilité caractérisant une mesure locale incertaine une pseudodonnée ponctuelle en accord avec la loi de probabilité dont elle est issue et le modèle de variabilité spatiale. L'avantage de cette approche est sa rapidité et sa complète cohérence vis à vis de la méthode de déformation graduelle. De fait, il devient possible de mettre en #uvre des processus d'optimisation au cours desquels on déforme graduellement la réalisation et le jeu de pseudo-données ponctuelles. En variant ce ou ces deux paramètre(s), on explore une chaîne de réalisations respectant toutes la structure spatiale requise et une chaîne de jeux de données respectant tous la structure spatiale et les lois de probabilités requises. Il s'agit alors d'identifier la réalisation et le jeu de pseudo-données qui minimisent la fonction objectif. La méthode permet d'arriver plus rapidement à la formation d'un modèle numérique représentatif du milieu. Présentation des figures D'autres caractéristiques et avantages de la méthode selon l'invention apparaîtront plus précisément à la lecture de la description ci-après d'un exemple non limitatif d'application, en se référant aux dessins annexés où : - la figure 1 montre à gauche (h(x,y)# N(FK, #K).N(m,#) et à droite les fonctions densité et de répartition en résultant. La fonction densité est obtenue à partir de la diagonale de h(x,y) ; - la figure 2 montre deux jeux (t = 0. 0 et t = 0. 5) de pseudo-données ponctuelles issues de la transformation de lois de probabilité uniformes (cas des gaussiennes tronquées) et les pseudo-données qui en sont déduites par déformation graduelle (t = 0. 1, 0. 2, 0.3, 0.4) ; - la figure 3 montre le processus d'optimisation développé pour construire des réalisations contraintes aux données statiques (de type lois de probabilité) et dynamiques par déformation graduelle ; la figure 4 montre les distributions de faciès pour le réservoir de référence, le réservoir pris comme point de départ de l'optimisation et les réservoirs obtenus en fin d'optimisation en ayant modifié ou non les pseudo-données déduites de la conversion des observations de faciès aux puits ; la figure 5 montre les flux fractionnaires simulés aux puits producteurs et les pressions au puits injecteur ; et - la figure 6 montre l'évolution de la fonction objectif pour des optimisations menées sans modifier les pseudo-données aux puits (Conventional Matching) et en les modifiant (New Matching). Description détaillée de la méthode La méthode selon l'invention permet à chaque itération du processus de recherche du minimum de modifier graduellement la réalisation elle-même, mais aussi les pseudodonnées issues du processus de transformation des lois de probabilités traduisant l'incertitude aux points de mesure. Cette dernière propriété donne plus de flexibilité au processus d'optimisation et permet de tendre vers le minimum plus rapidement. Transformation des lois de probabilités associées aux points de mesure en données ponctuelles Une étape préliminaire au conditionnement repose sur la transformation de l'information fournie sous forme de lois de probabilité en pseudo-données ponctuelles. Cette transformation doit permettre de respecter la structure spatiale et les lois de probabilité dont les données sont issues. Elle doit en outre être compatible avec la méthode de déformation graduelle. Pour cette raison, on utilisera pour la fonction de répartition normale standard et son inverse, notées G et les approximations analytiques rappelées par - Deutsch et Journel, GSLIB - Geostatistical software library and user's guide, Oxford univ. Press, 1992. Ces fonctions assurent une relation bijective entre un bruit blanc gaussien initial et un jeu de données transformées. Les résultats obtenus sont précis jusque 5 décimales. Considérons N points xl,l#[1,N] pour lesquels on dispose de mesures. Ces N mesures sont incertaines et définies, de fait, à partir de lois de probabilité Pl#[1,N]. Le cadre général mis en place pour le processus de transformation proposé repose sur la technique de simulation séquentielle. L'algorithme de transformation procède comme suit: 1) On tire au hasard un bruit blanc gaussien zl,l#[1,N] dont chaque composante est associée à un point #l,l#[1,N]. 2) On convertit ce bruit blanc gaussien en nombres uniformes indépendants 3) On définit un chemin aléatoire visitant chacun des points n+1, 4) A l'itération standards on suppose que les n lois de probabilités des n points visités précédemment ont été transformées en n valeurs normales Fl,l#[1,n]. Au point xn+1, on détermine alors à partir des n valeurs déjà définies l'estimateur de krigeage FK et l'écart type associé #K. 5) On détermine la fonction de répartition H correspondant à la loi de probabilité suivante : où N(FK,#K) est la loi de probabilité normale de moyenne FK et d'écart type #K. On estime ensuite Fn+1 = H-1[un+1]. On ajoute ensuite cette valeur au jeu de données transformées. 6) On se rend sur le point suivant défini à partir du chemin aléatoire, et on répète les étapes 4 et 5. On procède ainsi jusqu'à ce que les N lois de probabilités aient été transformées. Quelques cas particuliers a) Loi de probabilité uniforme On suppose que les lois de probabilité Pi#[1,N] correspondent à des distributions uniformes sur des intervalles [An+i, Bn+1]. A l'étape 5, on procède alors comme suit: Ce cas correspond à la méthode des gaussiennes tronquées. La méthode selon l'invention permet alors de contraindre la réalisation gaussienne sous-jacente à la réalisation gaussiennes tronquée à des pseudo-données représentatives des observations de faciès en certains points. En d'autres termes, on contraint ainsi la réalisation en faciès aux faciès observés aux puits. b) Loi de probabilité normale On suppose que la loi de probabilité décrivant l'incertitude sur la mesure est une loi normale de moyenne m et de déviation standard #. Il faut alors déterminer la loi de probabilité suivante : où a est une constante de normalisation, ainsi que la fonction de répartition H correspondante, ce qui se peut se faire à partir de techniques numériques (cf. Figure 1). Ce cas convient à la description de tout champ dans la mesure où il est ramené par anamorphose à un champ gaussien. Applications La Figure 2 illustre le processus de transformation proposé dans le cas de lois de probabilité uniformes. On suppose qu'en 100 points séparés de 1 m, on dispose de l'information "en ce point, l'attribut modélisé appartient à l'intervalle 1, à l'intervalle 2 ou à l'intervalle 3". En outre, la structure spatiale est caractérisée par un variogramme exponentiel avec une longueur de corrélation de 20 m. En appliquant l'algorithme exposé plus haut, on transforme ces observations en données ponctuelles (courbe t = 0. 0). Toutes respectent les intervalles dont elles sont extraites ainsi que le modèle de variabilité spatiale. En prenant comme point de départ un nouveau bruit blanc gaussien, on obtient un nouveau jeu de données ponctuelles (courbet = 0. 5). L'avantage de l'algorithme de transformation décrit est sa compatibilité avec la méthode de déformation graduelle. On peut en effet appliquer le formalisme de la déformation graduelle pour combiner deux bruits blancs gaussiens zl,l#[1,N] qui fournissent chacun un jeu de pseudo-données en accord avec le modèle de variabilité spatiale et les lois de probabilité liées aux points de mesure. Un exemple est présenté sur la Figure 2 dans le cadre d'une loi de probabilité uniforme. On combine les bruits blancs gaussiens ayant conduit aux courbes t = 0. 0 et t = 0. 5. En variant le paramètre de déformation t, on obtient d'autres jeux de données ponctuelles respectant la structure spatiale et les lois uniformes requises. Conditionnement par les données statiques et dynamiques L'avantage de l'algorithme de transformation présenté dans la section précédente est son entière compatibilité avec la méthode de déformation graduelle. Or, cette méthode est une technique de paramétrage très pratique dans le cadre de l'optimisation stochastique : elle permet de déformer la distribution d'un attribut à l'aide d'un petit nombre de paramètres tout en préservant la structure spatiale. Nous intégrons la technique de transformation exposée plus haut dans le processus d'optimisation, selon le schéma décrit sur la Figure 3. On considère au départ un bruit blanc gaussien initial et un bruit blanc gaussien complémentaire pour la réalisation ; de même pour les observations statiques. Dans les deux cas, on combine suivant les principes de déformation graduelle le bruit initial au bruit complémentaire, cette combinaison étant contrôlée par un paramètre de déformation t. En réalité, on pourrait prendre deux paramètres de déformation différents, un pour la réalisation et un pour les observations statiques. On suppose ici ces deux paramètres identiques pour se retrouver dans le cadre d'un problème à une dimension. Ensuite, pour la réalisation, le bruit blanc gaussien issu de la combinaison graduelle est transformé en réalisation gaussienne structurée (composant FFTMA). En parallèle, pour les observations statiques, le bruit blanc gaussien issu de la combinaison graduelle est transformé en pseudo-données ponctuelles en accord avec les lois de probabilité observées. Ces pseudodonnées sont alors utilisées pour contraindre la réalisation fournie par le composant FFTMA. Une simulation d'écoulement est exécutée pour la réalisation contrainte qui fournit des données de production (par exemple, des pressions, des temps de percée, des débits, etc. ). Ces données simulées sont ensuite comparées aux données dynamiques observées sur le terrain par le biais de la fonction objectif. Durant le processus d'optimisation, on modifie le paramètre t afin de réduire la fonction objectif. Puis, cet technique de recherche est poursuivie plus avant avec de nouveaux bruits blancs gaussiens complémentaires. Exemple numérique On construit un modèle synthétique de réservoir sur lequel on teste la méthode selon l'invention. On considère un réservoir en faciès pour lequel on a observé la nature des faciès aux puits. Ce réservoir est simulé par la méthode des gaussiennes tronquées : les faciès sont alors intimement liés à lois de probabilités uniformes, c'est à dire à des intervalles dont la largeur dépend de la proportion des dits faciès. La technique des gaussiennes tronquées consiste à appliquer des seuillages suivant les intervalles définis à une réalisation gaussienne standard continue pour la transformer en réalisation discrète. Aux points d'observation, les valeurs de la réalisation gaussienne continue sont contraintes par les intervalles caractérisant les faciès observés. Le réservoir synthétique de référence est représenté sur la Figure 4. Il s'agit d'un réservoir mono-couche comprenant 100x100 mailles de 10 m d'épaisseur et de 10 m de côté. Le réservoir comprend trois faciès : 25 % de faciès 1,35 % de faciès 2 et 40 % de faciès 3. Leurs perméabilités sont 300 mD, 200 mD et 50 mD, respectivement. La réalisation gaussienne utilisée pour obtenir la réalisation tronquée est caractérisée par un variogramme anisotrope, stable, d'exposant a = 1.5: h est la distance, 1, la longueur de corrélation et # l'écart type. La longueur de corrélation vaut 50 m suivant l'axe principal (1;1;0) et 20 m suivant l'axe perpendiculaire (- 1;1;0). La porosité est constante et vaut 0,4. On a un puits où on injecte de l'eau au centre et quatre puits producteurs dans les coins. Le réservoir est supposé saturé initialement en huile. Les courbes de perméabilités relatives pour l'huile et l'eau suivent des lois de Corey avec un exposant de 2. Le rapport de mobilité est 1. L'historique de production pour ce réservoir de référence est présenté sur la Figure 5. Les faciès observés aux puits sont donnés dans le tableau 1. Table 1. Faciès observés aux puits. Well Facies INJ1 3 PROD1 2 PROD2 2 PROD3 2 PROD4 1 L'objet du problème inverse est de déterminer un modèle de réservoir cohérent avec les données dynamiques et les observations de faciès aux puits en supposant inconnue la distribution des faciès. Dans ce but, on lance deux processus d'optimisation en partant de la même réalisation de départ (figure 4). Pour chacun de ces processus, on considère un unique paramètre d'optimisation, à savoir le paramètre de déformation. Le premier processus est basé sur une approche classique : les valeurs gaussiennes représentatives des observations de faciès aux puits sont constantes pendant l'optimisation. Pour le deuxième processus, on teste l'approche proposée en variant ces valeurs gaussiennes. On observe qu'en rendant possible les variations au niveau des valeurs gaussiennes représentatives des observations de faciès aux puits (lois de probabilité uniformes), on accélère significativement la minimisation de la fonction objectif (Figure 6). On a décrit jusqu'ici une mise en #uvre de la méthode dans le cadre préféré d'un processus de déformation graduelle des réalisations du modèle qui impose que les coefficients des combinaisons de réalisations qui sont effectuées, soient tels que la somme des carrés des coefficients de combinaison soit égale à 1. On ne sortirait pas du cadre général de l'invention toutefois en appliquant la technique de transformation qui vient d'être décrite, associant à une loi de probabilité caractérisant une mesure locale incertaine, une pseudo-donnée ponctuelle en accord avec la loi de probabilité dont elle est issue et un modèle de variabilité spatiale, à une autre approche d'optimisation du modèle stochastique connue en soi, telle que la méthode dite des points pilotes. REVENDICATIONS 1) Méthode pour former un modèle numérique stochastique de type gaussien ou apparenté, représentatif de la distribution d'une grandeur physique dans un milieu hétérogène poreux, qui soit ajusté par rapport à des données dynamiques, caractéristiques du déplacement des fluides dans le milieu, et des données statiques locales mesurées en un certain nombre de points de mesure le long de puits au travers du milieu, et présentant une certaine marge d'incertitude, dans laquelle on optimise le modèle par le biais d'un processus itératif de déformation où l'on forme à chaque itération une réalisation combinée obtenue par combinaison linéaire d'une part d'une réalisation initiale, représentant au moins une partie du milieu, et d'au moins une deuxième réalisation indépendante du même modèle stochastique, et on minimise une fonction objectif mesurant l'écart entre des données dynamiques réelles et les données dynamiques simulées au moyen d'un simulateur d'écoulement, pour la réalisation combinée, par ajustement des coefficients de combinaison, le processus d' ajustement itératif étant poursuivi jusqu'à l'obtention d'une réalisation optimale du modèle stochastique, caractérisée en ce que : a) on transforme les données statiques locales en pseudo-données ponctuelles en accord avec des lois de probabilité et un modèle de variabilité spatiale ; et b) on ajuste les pseudo-données en respectant les lois de probabilités dont elles sont issues, ainsi que le modèle de variabilité spatiale, par le biais d'un processus itératif où l'on combine un premier bruit blanc gaussien associé au jeu de pseudo-données à un deuxième bruit blanc gaussien. 2) Méthode selon la revendication 1, caractérisée en ce que l'on ajuste le modèle par un processus de déformation graduelle en imposant que la somme des carrés des coefficients de combinaison entre les réalisations soit égale à 1 et on ajuste les pseudodonnées par un processus itératif où la somme des carrés des coefficients de la dite combinaison vaut également 1. 3) Méthode selon la revendication 2, caractérisée en ce que les lois de probabilités sont des lois normales ou des lois uniformes. 4) Méthode selon la revendication 2 ou 3, caractérisée en ce que l'on conduit l'ajustement itératif à partir de deux paramètres de déformation, un premier paramètre qui contrôle la combinaison entre la réalisation initiale et la deuxième réalisation, et un deuxième paramètre qui contrôle la combinaison entre le bruit blanc gaussien initial et le deuxième bruit blanc gaussien. 5) Méthode selon la revendication 2 ou 3, caractérisée en que l'optimisation est conduite à partir d'un unique paramètre quand les coefficients de combinaison sont identiques pour la combinaison de réalisations et la combinaison de bruits blancs gaussiens. 6) Méthode selon la revendication 1, dans laquelle l'optimisation est conduite par une méthode de points pilotes.