FR 2846767 A1 20040507 FR 0213632 A 20021030 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 d'un milieu hétérogène poreux (tel qu'un gisement d'hydrocarbures par exemple) calé par rapport à des données dites dynamiques caractéristiques du déplacement des fluides dans le milieu telles que 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é 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. 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é. Ces réalisations forment des modèles numériques de réservoir. Les données disponibles sont, par exemple, des mesures ponctuelles de perméabilité ou porosité, un modèle de variabilité spatiale déterminé selon des mesures ponctuelles ou encore 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. La cohérence vis à vis des données 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. Pour ce faire, la technique la plus simple est celle de l'essai et de l'erreur. Cette approche consiste à tirer au hasard des réalisations jusqu'à l'obtention d'une réalisation satisfaisant aux données. Elle offre l'avantage de la conservation du modèle de variabilité spatiale, mais requiert un temps calcul prohibitif. Aussi n'est-elle que très rarement utilisée en pratique. Une alternative souvent préférée prend appui sur le calcul des gradients. Les méthodes de gradients permettent de modifier une réalisation initiale suivant une direction de recherche, cette dernière étant estimée à partir des gradients. Les modifications sont apportées itérativement jusqu'à ce que le calage des données soit jugé acceptable. Les méthodes de gradients sont attrayantes de par leur efficacité. Néanmoins, elles deviennent inadaptées dès lors que le nombre de paramètres, c'est-à-dire le nombre de valeurs de perméabilité et de porosité constituant le modèle numérique, est important. En outre, elles ne permettent pas de modifier les réalisations en respectant la structure spatiale du modèle stochastique. Plus récemment, une technique de paramétrage géostatistique, 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 itérative 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): Gradual déformation and itérative 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. Comme on va le rappeler plus en détail ci-après, la méthode de déformation graduelle permet de modifier graduellement une réalisation d'un modèle stochastique de type Gaussien ou de type apparenté au modèle Gaussien tout en respectant la structure spatiale du modèle. Optimisation stochastique Soit Jobs = (Jobs 1 , jObS 2 ,....... jobs) M l'ensemble des données sur le terrain et 1 =J1,J2,......,fM) les réponses correspondantes simulées pour une réalisation z d'un modèle stochastique Z donné. En général, les réponses f = (fl, f2,.....,fM ) sont obtenues en résolvant numériquement le problème direct. Ainsi, si z 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 Z 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 z. En ce sens, la minimisation de la fonction objectif est un problème à plusieurs variables. Soit N le nombre de mailles ou de composantes formant la réalisation z. N est souvent très grand (104 # 107). Il est donc très difficile de mener une optimisation directement par rapport aux composantes de z. De plus, la réalisation z, même modifiée, doit rester une réalisation de Z. La méthode de déformation graduelle permet de s'affranchir de ces difficultés. Recherche aléatoire à partir de la méthode de déformation graduelle On considère à présent un champ aléatoire Z qui peut être transformé en champ aléatoire Gaussien Y. 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 (figure la). 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 + ul sin t où t est le paramètre de déformation. Une autre technique de construction de chaînes de réalisations consiste à combiner la réalisation initiale avec non pas une, mais P réalisations complémentaires up (p=l,P) de Y (figure lb). Les coefficients de la combinaison sont tels que la somme de leurs carrés est égale à 1. 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 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 y(topt) 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 + ul 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. Ces dernières sont aussi indépendantes de yo. Cette formulation implique que la chaîne de réalisations correspond à une hyper ellipse de dimension N. 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. Le coté aléatoire de la méthode tient au fait qu'à chaque itération, la réalisation complémentaire est tirée au hasard. De fait, la direction de recherche dans laquelle on part depuis la réalisation optimisée à l'étape précédente est aléatoire. En effet, la direction de recherche, pour une chaîne donnée et à partir de la réalisation optimale définie précédemment, est : Cette direction de recherche dépend uniquement de ul. De plus, ul étant indépendante des réalisations complémentaires déjà générées ul,u2,...,ul~1 et aussi de yo, la direction de recherche au départ de chaque nouvelle chaîne est orthogonale à la tangente définie pour la chaîne précédente au même point (figure 2). Bien qu'il puisse paraître adéquat de sélectionner une direction de recherche orthogonale à cette tangente, il existe une infinité de directions orthogonales possibles. L'expérience montre que, après quelques itérations, les nouvelles directions de recherche ne contribuent plus significativement à la décroissance de la fonction objectif (figure 6). Il a été envisagé aussi de combiner la réalisation initiale non pas à une, mais à plusieurs réalisations complémentaires. Dans ce cas, le nombre de paramètres de déformation augmente : il égale le nombre de réalisations complémentaires mises en jeu lors d'une combinaison graduelle. Bien que le processus d'optimisation soit alors plus flexible, il faut pouvoir gérer plusieurs paramètres ce qui n'est pas forcément évident. Par ailleurs, un tel processus n'est pas nécessairement plus efficace, car il peut dépendre de l'exécution d'un plus grand nombré de simulations directes d'écoulement. La méthode selon l'invention La méthode selon l'invention permet de former plus rapidement un modèle numérique représentatif de la distribution d'une grandeur physique dans un milieu hétérogène poreux, tel qu'une zone souterraine (réservoirs pétroliers, aquifères, etc.), contraint par des données collectées dans le milieu (données dynamiques caractéristiques du déplacement des fluides dans le milieu, collectées par des mesures (dans des puits de production, d'injection ou d'observation par exemple) ou des observations préalables. Elle comporte un processus itératif de déformation graduelle où l'on combine linéairement à une réalisation initiale d'une partie au moins du modèle choisi du milieu, au moins une deuxième réalisation indépendante de la réalisation initiale, les coefficients de cette combinaison linéaire étant tels que la somme de leurs carrés est égale à 1, et l'on minimise une fonction objectif mesurant l'écart entre un jeu de données non linéaires déduites de la dite combinaison au moyen d'un simulateur du milieu, et les dites données géologiques et dynamiques par ajustement des coefficients de la combinaison, le processus itératif étant répété jusqu'à obtenir une réalisation optimale du modèle stochastique. La méthode est caractérisée en ce que l'on accélère la vitesse de déformation graduelle vers le modèle optimal représentatif du milieu, en sélectionnant comme deuxième réalisation à combiner avec la réalisation initiale, au moins une réalisation composite obtenue en sélectionnant au préalable une direction de descente définie en fonction des gradients de la fonction objectif par rapport à toutes les composantes de la dite réalisation initiale. On obtient la réalisation composite par exemple par combinaison linéaire d'un jeu de réalisations indépendantes du modèle, les coefficients de la combinaison étant calculés de façon que la direction de descente à partir de la réalisation initiale y soit le plus proche possible de celle définie par les gradients de la fonction objectif par rapport à toutes les composantes de la réalisation initiale. On conduit l'optimisation par exemple à partir d'un paramètre de déformation, qui contrôle la combinaison entre la réalisation initiale et la réalisation composite. Dans le cas où la dite combinaison n'affecte qu'une partie de la réalisation initiale, on applique le processus itératif de déformation graduelle à un bruit blanc Gaussien utilisé pour générer une réalisation gaussienne et on détermine les dérivées de la fonction objectif par rapport aux composantes du bruit blanc Gaussien. Suivant un mode de mise en #uvre, on combine la réalisation initiale avec un certain nombre M de réalisations composites, toutes obtenues par composition à partir de Pm réalisations indépendantes de Y, l'optimisation faisant intervenir M paramètres. En d'autres termes, la méthode comprend essentiellement un nouveau schéma de combinaison graduelle qui tient compte de l'information apportée par les gradients au point initial de toute chaîne de réalisations. Pour construire une chaîne, on part toujours d'une réalisation de départ et d'un ensemble de réalisations complémentaires, toutes indépendantes et issues du même modèle stochastique. Par contre, on ne combine pas directement la réalisation de départ avec les réalisations complémentaires. Les réalisations complémentaires offrent la possibilité d'explorer l'espace des réalisations dans des directions différentes. Ces directions ne sont pas équivalentes : certaines permettent de se diriger davantage vers l'optimum. A ce stade, on élabore une réalisation, dite composite, en combinant les réalisations complémentaires uniquement. Puis, on crée une chaîne de réalisations à partir de la réalisation de départ (initiale) et de cette réalisation composite. Cette chaîne, comme celle qui était proposée dans le cas de base de la déformation graduelle, peut être explorée à l'aide d'un seul paramètre de déformation. La réalisation composite est construite de façon à proposer une direction de recherche le long de la chaîne aussi proche que possible de la direction de descente donnée par les gradients. Comme mentionné plus haut, toutes les réalisations complémentaires ne sont pas équivalentes : la réalisation composite prend le meilleur de chaque réalisation complémentaire. La méthode permet d'arriver plus rapidement à la formation du 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ù : - les figures la, lb montrent des schémas de déformation graduelle (dits GDM) déjà connus ; - la figurelc montre le schéma de déformation graduelle (dit GBC) correspondant à la méthode selon l'invention ; - la figure 2 représente des chaînes de réalisations dans un espace euclidien à N dimensions, où la tangente au niveau de la réalisation optimisée pour une chaîne de réalisations l-1 (RCl-1) est orthogonale à la direction de recherche pour la réalisation de départ de la chaîne de réalisations suivante l (RC1 ); - la figure 3 montre la projection de la direction de descente v dans le sous-espace U défini par la base orthonormale formée par P réalisations indépendantes (Ul ,..., u p ) ; - les figures 4A à 4C montrent des exemples de distributions de perméabilité respectivement pour les réalisations de référence, initiale et contrainte aux données de pression; - les figure 5A à 5E montrent respectivement les variations en fonction du temps des pressions simulées respectivement aux cinq puits (BHP-OBS1, BHP-OBS2, BHP- PRO1, BHP-OBS3, BHP-OBS4) des figures 4, respectivement pour les distributions de perméabilité de référence (data), initiale et contrainte (match) ; et - la figure 6 montre l'évolution (OF) de la fonction objectif en fonction du nombre k de simulation d'écoulement exécutées, GDM1 correspondant à une optimisation menée en combinant une réalisation de départ et une seule réalisation complémentaire, GDMGBC3, à une optimisation menée en combinant la réalisation de départ et une réalisation composite construite à partir de trois réalisations complémentaires, GDMGBC10, à l'optimisation menée en combinant la réalisation de départ et une réalisation composite construite à partir de dix réalisations complémentaires, et GDMGBC30, à une optimisation menée en combinant la réalisation de départ et une réalisation composite construite à partir de trente réalisations complémentaires. Description détaillée de la méthode La méthode selon l'invention permet d'orienter à chaque itération la construction de la chaîne de réalisations afin d'atteindre une direction de recherche souhaitable. La technique retenue tire avantage de l'information apportée par la direction de descente définie par les gradients de la fonction objectif J. Elle peut être mise en #uvre au moyen d'un simulateur numérique d'un type connu des gens de l'art tel que les simulateurs ATHOS ou ECLIPSE. Recherche orientée par la direction de descente On considère à ce stade la direction de descente (évolution vers une valeur minimale qu'elle soit ou non locale) définie par les gradients de la fonction objectif J par rapport à toutes les composantes de la réalisation z. Ces gradients se déduisent des coefficients de sensibilité #fm / #zn : Le problème du calcul des coefficients de sensibilité a été largement abordé dans la littérature scientifique. On pourra, par exemple, consulter : - Sun, N. Z., Inverse problems in groundwater modelling, Kluwer Acad. Publ., Dordrecht, The Netherlands, 1994. En particulier, la technique de l'état adjoint permet de calculer l'ensemble de ces coefficients en résolvant le problème direct et son problème adjoint, dans un temps équivalent à deux fois celui requis pour le résolution du problème direct. Comme l'objectif est d'intégrer l'information fournie par les gradients dans la méthode de déformation graduelle, il faut se ramener à la réalisation Gaussienne qui résulte de la transformation de z. On distingue deux cas selon 1) que l'on dispose ou non de données ponctuelles de conditionnement, c'est-à-dire de mesures dez en certains points. Dans le cas où l'on ne dispose d'aucune donnée, si y est la réalisation Gaussienne obtenue en transformant z, alors : Dans le cas contraire, on suppose que z est connu en certains points a. Ainsi la réalisation Gaussienne s, déduite de la transformation de z, est une réalisation conditionnelle de Y. Pour obtenir s, on génère une réalisation non conditionnelle y de Y et on la conditionne aux valeurs connues en a par krigeage. La réalisation conditionnelle se déduit de : s = s*+ (y-y*) s* et y* sont calculées respectivement par krigeage à partir des données réelles et les données de y générées au niveau des points a. On montre alors que : Pour des propriétés physiques continues, dzn/dyn ou dzn I dsn se calculent à partir de la fonction anamorphose qui permet de transformer la réalisation.z en réalisation Gaussienne. Lorsque les propriétés physiques considérées sont catégorielles ou discrètes, ces dérivées n'existent pas. Les techniques de gradients ne peuvent alors s'appliquer. Ces différentes relations sont d'un intérêt direct si on se contente de pouvoir déformer globalement une réalisation. Si on souhaite par contre la déformer localement, il faut appliquer la méthode de déformation graduelle au bruit blanc Gaussien x utilisé pour générer y. Dans ce cas, une étape supplémentaire s'impose : le calcul des dérivées de la fonction objectif par rapport aux composantes du bruit blanc Gaussien. Pour illustrer le calcul de ces dérivées, on propose de se concentrer sur le cas particulier suivant lequel la réalisation Gaussienne non conditionnelle y est obtenue à partir du générateur FFT-MA décrit dans l'article publié par Le Ravalec et al. 2000 cité plus haut. Le principe de base de ce générateur FFT-MA (pour "FFT-Moving Average") est de transformer un bruit blanc Gaussien x en réalisation Gaussienne y corrélée à partir d'un produit de convolution : y=g*x La fonction g résulte de la décomposition de la fonction de covariance C telle que C = g * g', où gt est la transposée de g. Les dérivées de la fonction objectif par rapport aux composantes du bruit blanc Gaussien sont : L'expression discrète du produit de convolution conduit à y, = Y, 91-k Xk ce qui k implique #yl/#xn -gl-n. Si on introduit cette formulation dans les dérivées de la fonction objectif, on montre : Cette formule exprime le fait que la dérivée de la fonction objectif par rapport au nième composant du bruit blanc Gaussien est donnée par la nième composante du champ obtenu en convoluant toutes les dérivées de la fonction objectif par rapport aux composantes de la réalisation Gaussienne avec le kernel de la fonction de covariance. A partir du cadre mis en place pour le générateur FFT-MA, on procède comme suit pour déterminer ces dérivées. 1-calcul de la transformée de Fourier de #J/#y, ces dérivées étant fournies par mise en #uvre du simulateur numérique direct ; 2-multiplication de cette transformée de Fourier avec celle de g qui est fournit par FFT-MA; 3-calcul de la transformée de Fourier inverse du produit précédent. Le temps requis pour le calcul de ces dérivées est négligeable : il représente un temps additionnel de 2/3 par rapport à la simulation d'une réalisation Gaussienne par FFTMA. Quels que soient le générateur de réalisations et le simulateur numérique direct, on suppose, dans la suite, que l'on est capable de définir une direction de descente à partir des dérivées de la fonction objectif. Si l'optimisation de la fonction objectif est menée par rapport à cette direction de descente seulement, on détruit en général la cohérence de la réalisation vis à vis du modèle de variabilité spatiale. Dans la section qui suit, on intègre l'information apportée par les dérivées de la fonction objectif dans le schéma de combinaison graduelle. Prise en compte des dérivées de la fonction objectif dans le processus de déformation graduelle On considère la chaîne de réalisations yl(t) construite à partir de yo et d'une autre réalisation u de Y (figure la). A présent, au lieu d'utiliser une réalisation complémentaire u telle qu'elle, on tire au hasard P réalisations complémentaires up (p=1,2,...,P) de Y et on substitue à u une combinaison des p=1,P (figure le). Cette combinaison n'est pas quelconque : elle suit le mode de construction suivant : La réalisation u résultante est une réalisation de Y et est indépendante de yo. D'après l'équation 1, u est aussi la direction de recherche calculée pour la chaîne y1(t) au point de départ yo. On propose donc de construire u de sorte que la direction u soit aussi proche que possible de la direction de descente donnée par les gradients en yo. On définit au préalable l'espace V des vecteurs à N dimensions muni du produit scalaire : y@,net yj,n sont respectivement les nièmes composantes des vecteurs y, et yj# Soit U un sous-espace de V défini par la base orthonormale (u1,u2,...,uP). La direction de recherche dans U la plus proche de la direction de descente v est donnée par la projection ' de v dans U (figure 3) : En normalisant ce vecteur, on obtient la direction u souhaitée. Les coefficients de combinaison À de l'équation (2) sont donc : La réalisation u ainsi définie est appelée réalisation composite. On a considéré jusqu'ici la construction de chaînes de réalisations par combinaison la réalisation initiale avec une réalisation composite. Cette technique peut toutefois être généralisée à la construction de chaînes faisant intervenir un certain nombre M de réalisations composites, toutes obtenues par composition à partir de Pm réalisations de Y, ces #m=1M Pm réalisations étant indépendantes. L'optimisation fait alors intervenir M paramètres. Cette technique augmente le degré de liberté dans le processus d'optimisation, mais oblige à gérer M paramètres d'optimisation. Exemple numérique On construit un modèle synthétique de réservoir sur lequel on teste la méthode selon l'invention. Le réservoir synthétique de référence est représenté sur la figure 4A. Il s'agit d'un réservoir mono-couche comprenant 51x51 mailles de 10 m d'épaisseur et de 40 m de côté. La distribution de perméabilité est lognormale de moyenne 200 mD et d'écart-type 100 mD. Le logarithme du champ de perméabilité est caractérisé par un variogramme sphérique isotrope et une longueur de corrélation de 480 m. Les autres propriétés pétrophysiques sont constantes : la porosité est de 25 %, la compressibilité totale de 5.10-4 bar-1et la viscosité du fluide de 1 cP. Un puit de production BHP-PROl, de rayon 7. 85 cm dépourvu d'effet pariétal ("skin effect") se trouve au centre du réservoir : il est entouré par quatre puits d'observation (BHP-OBS1, BHP-OBS2, BHP-OBS3, BHP-OBS4) (figures 4). Une simulation numérique d'écoulement permet d'obtenir pour ce réservoir un jeu de données de référence comprenant les pressions aux cinq puits (figure 5). L'objet du problème inverse est de déterminer un modèle de réservoir cohérent avec les données de pression en supposant inconnue la distribution des perméabilités. Dans ce but, on lance quatre processus d'optimisation en partant de la même réalisation de départ (figure 4B). Pour chacun de ces processus, on considère un unique paramètre d'optimisation, à savoir le paramètre de déformation. Le premier processus (GDM1) reprend le schéma de déformation graduelle classique avec la construction d'une chaîne de réalisation en utilisant la réalisation de départ et une seule réalisation complémentaire. Les trois autres processus (GDMGBC3, GDMGBC10 et GDMGBC30) illustrent l'application de la méthode selon l'invention : les chaînes sont dans ce cas élaborées avec la réalisation de départ et une réalisation composite issue de la combinaison de 3,10 et 30 réalisations complémentaires. La réalisation composite est construite comme expliqué ci-dessus (Eq. 2) en intégrant l'information fournie par les gradients de la fonction objectif par rapport au bruit blanc Gaussien. Pour chacun des processus, on représente sur la figure 6 l'évolution de la fonction objectif (OF) en fonction du nombre k de simulation d'écoulement exécutées. On observe qu'en utilisant les gradients et en augmentant le nombre de réalisations complémentaires, la fonction objectif diminue plus vite pour un même nombre de simulations d'écoulements. REVENDICATIONS 1) Méthode pour former plus rapidement 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 calé par rapport à des données dynamiques obtenues par des mesures effectuées dans le milieu ou des observations préalables, et caractéristiques du déplacement des fluides dans celui-ci, comportant un processus itératif de déformation graduelle où l'on combine linéairement à chaque itération, une réalisation initiale (y) d'au moins une partie du milieu, et au moins une deuxième réalisation indépendante de la réalisation initiale, les coefficients ([alpha]i) de cette combinaison linéaire étant tels que la somme de leurs carrés est égale à 1, et l'on minimise une fonction objectif (J) mesurant l'écart entre un jeu de données non linéaires déduites de la dite combinaison au moyen d'un simulateur d'écoulement dans le milieu, et les dites données géologiques et dynamiques par ajustement des coefficients de la combinaison, le processus itératif étant répété jusqu'à obtenir une réalisation optimale du modèle stochastique, caractérisée en ce que l'on accélère la vitesse de déformation graduelle vers le modèle optimal représentatif du milieu, en sélectionnant comme deuxième réalisation à combiner avec la réalisation initiale, au moins une réalisation composite obtenue en sélectionnant au préalable une direction de descente définie en fonction des gradients de la fonction objectif par rapport à toutes les composantes de la dite réalisation initiale (y). 2) Méthode selon la revendication 1, caractérisée en ce que l'on obtient la réalisation composite par combinaison linéaire d'un jeu de P réalisations indépendantes du modèle, les coefficients de la combinaison étant calculés de façon à ce que la direction de descente à partir de la réalisation initiale y soit le plus proche possible de celle définie par les gradients de la fonction objectif par rapport à toutes les composantes de la réalisation initiale. 3) Méthode selon la revendication 1 ou 2, caractérisée en ce que l'on conduit l'optimisation à partir d'un paramètre de déformation, qui contrôle la combinaison entre la réalisation initiale et la réalisation composite. 4) Méthode selon la revendication 1 ou 2, caractérisée en ce que la dite combinaison n'affectant qu'une partie de la réalisation initiale, on applique le processus itératif de déformation graduelle à un bruit blanc Gaussien utilisé pour générer une réalisation gaussienne et on détermine les dérivées de la fonction objectif par rapport aux composantes du bruit blanc Gaussien. 5) Méthode selon la revendication 1 ou 2, caractérisée en ce que l'on combine la réalisation initiale avec un certain nombre M de réalisations composites, toutes obtenues par composition à partir de Pm réalisations indépendantes de Y, l'optimisation faisant intervenir M paramètres.