FR2852710A1 - Methode pour former rapidement un modele stochastique representatif de la distribution d'une grandeur physique dans un milieu heterogene par une selection appropriee de realisations geostatistiques - Google Patents

Abstract

- Méthode pour former rapidement un modèle stochastique 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, par une sélection appropriée de réalisations géostatistiques à combiner linéairement.- On utilise un processus itératif de déformation graduelle où l'on combine linéairement à chaque itération, une réalisation géostatistique initiale (y) du milieu, et un nombre (N-1) (N>1) autres réalisations (zi) i=1,....(N-1), indépendantes de la réalisation initiale (y), en imposant des contraintes aux coefficients de combinaison linéaire des réalisations (y) et (zi) i=1,...,(N-1), et l'on minimise une fonction objectif (J) mesurant l'écart entre un jeu de données simulées déduites de la dite combinaison au moyen d'un simulateur, et les dites données dynamiques. Pour minimiser rapidement la fonction objectif, on calcule le multiplicateur de Lagrange (λ), λ∈RN-1 associé à la contrainte portant sur les valeurs spécifiées des coefficients des réalisations (zi) i=1,...,(N-1). La valeur absolue λi de la i-ième composante du multiplicateur λ indique la sensibilité de la fonction objectif par rapport à la i-ième réalisation géostatistique (zi). Par ce moyen, avec un choix judicieux de coefficients initiaux, l'utilisateur peut donc déterminer à priori, si une réalisation donnée fera baisser la fonction objectif de manière significative (au premier ordre).- Applications notamment aux problèmes de calage d'historiques de production en ingénierie de réservoirs par exemple.

Description

La présente invention concerne une méthode pour former rapidement un

modèle stochastique représentatif de la distribution d'une grandeur physique telle que la perméabilité par exemple, dans un milieu hétérogène poreux, qui soit calé par rapport à des données dynamiques, par une sélection appropriée de réalisations géostatistiques combinées linéairement.

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 la 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 assurée dans le modèle par le biais d'une procédure inverse.

Les études intégrées en ingénierie de réservoirs ont principalement deux objectifs: - sur des champs matures, l'ingénieur réservoir cherche à quantifier les incertitudes liées aux prévisions de production, - sur des projets en cours de développement ou entrant dans une nouvelle phase de production, l'ingénieur réservoir veut pouvoir tester différents scénarii de production pour effectuer des études de risques.

Dans ce contexte, l'utilisation des géostatistiques aussi bien comme méthodes d'estimation que de simulation stochastiques est devenue courante. Les outils actuels de simulation géostatistique permettent de générer rapidement des modèles de réservoirs pétroliers contenant plusieurs millions de mailles. Les défis associés à l'utilisation de tels modèles sont principalement de deux ordres: [1] d'une part, il faut pouvoir intégrer les données disponibles pour mettre à jour le modèle de réservoir tout en assurant la conservation des propriétés géostatistiques du modèle géologique initial, [2] d'autre part, il faut résoudre le problème inverse associé à cette intégration de données dans des délais compatibles avec les contraintes économiques.

Dans les deux cas, la paramétrisation du modèle géologique joue un rôle primordial. Une approche classique consiste à réduire le nombre de paramètres et à ne prendre en compte que ceux ayant une sensibilité maximale. La méthodes des points pilotes, initialement présentée par Ramarao, B.S., LaVenue, A.M., de Marsily, G. & Marietta, M.G., "Pilot point methodology for automated calibration of an ensemble of conditionally simulated transmissivity fields: 1. Theory and computational experiments", Water Resources Research, 31(3): 475-493, March 1995 permet d'effectuer un calage d'historique sur un modèle de réservoir paramétré par les données de puits et un certain nombre de points pilotes spécifiés par l'utilisateur. La mise à jour du modèle est néanmoins localisée au voisinage des points pilotes.

La méthode de zonation du réservoir en principales unités ou zones ayant des propriétés pétrophysiques constantes a été initialement présentée par - Bissel, R., "Calculating optimal parameters for history matching", Proceedings of the 4th European Conference on the mathematical of oil recovery (ECMOR IV), 1994.

Elle permet d'effectuer un calage dhistorique pour peu que la zonation choisie soit correcte du point de vue géologique. Le point [1] n'est néanmoins pas respecté.

Une autre approche réside dans la paramétrisation multi-échelles pour laquelle le problème est résolu successivement sur des échelles de plus en plus fines. Elle présente l'inconvénient de conduire généralement à une sur-paramétrisation car tous les degrés de liberté de l'échelle inférieure sont utilisés alors que seuls certains seraient nécessaires pour expliquer les données.

Les approches multi-échelles adaptatives permettent de corriger cet inconvénient. Le concept d'indicateur de raffinement, présenté par: Chavent, G., &, Bissell, R., "Indicator for the refinement of parameterization", Proceedings of the International Symposium on Inverse Problems in Engineering Mechanics, Nagano, Japan, p. 185-190, 1998 permet d'identifier les degrés de liberté utiles pour expliquer les données tout en évitant le piège de la sur-paramétrisation.

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, Elle est aussi décrite dans la publication suivante: - Roggero, F., &, Hu, L.Y., (1998) "Gradual deformation of continuous geostatistical models for history matching", SPE 49004.

Elle permet d'effectuer un calage d'historique tout en conservant les propriétés géostatistiques initiales du modèle de réservoir. La paramétrisation se réduit alors aux paramètres de déformation graduelle à partir desquels l'utilisateur calcule les coefficients de la combinaison linéaire. Le modèle géologique étant paramétré par une combinaison linéaire de réalisations géostatistiques, le point [1] peut être satisfait par une contrainte spécifique sur les coefficients de cette combinaison linéaire.

Le point [2] implique deux autres conditions dès lors que l'on travaille avec la méthode des déformations graduelles. Il s'agit de savoir: [3] quel nombre N de réalisations géostatistiques est à considérer pour la combinaison linéaire; et [4] comment choisir de la manière la plus efficace possible les N réalisations géostatistiques optimales (décroissance la plus rapide de la fonction objectif du problème inverse considéré).

Une approche classique concernant le point [3] consiste à effectuer une décomposition en valeurs et vecteurs propres. Les différentes valeurs propres obtenues permettent alors de trouver un compromis entre l'incertitude obtenue sur les paramètres en fin de calage et le nombre de paramètres que l'on peut estimer à partir des données disponibles.

Suivant une autre approche, toujours dans le cadre d'une paramétrisation par la méthode des déformations graduelles, on considère que toute réalisation géostatistique apporte toujours une contribution, même minime, à la baisse de la fonction objectif (J). En conséquence, on combine linéairement par la méthode des déformations graduelles un certain nombre de réalisations géostatistiques optimales. Ces réalisations optimales sont elles-mêmes issues d'une combinaison linéaire de réalisations géostatistiques initiales dont les coefficients de combinaison sont choisis de telle manière à proposer une direction de recherche des déformations graduelles aussi proche que possible de la direction de descente donnée par les gradients. Cette approche a fait l'objet de la demande de brevet 02/13.632 du demandeur.

Concernant le point [4], aucune approche n'a été jusqu'à présent proposée permettant d'effectuer un choix à priori des réalisations (ou cartes) utilisées dans la combinaison linéaire dans un contexte de calage d'historiques. Seul une technique permettant un choix à priori des réalisations géostatistiques correspondant aux scenarii de production extrêmes dans un contexte de quantification des incertitudes a été proposée par: - Roggero, F., "Direct sélection of stochastic model realizations constrained to historical data", SPE 38731, 1997.

LA METHODE SELON L'INVENTION La méthode selon l'invention permet de former un modèle numérique stochastique de type Gaussien ou apparenté, donnant une image de la distribution d'une grandeur physique dans un milieu hétérogène poreux, qui soit calé par rapport à des données obtenues par des mesures effectuées dans le milieu ou des observations préalables, et caractéristiques du déplacement des fluides dans le milieu. Elle comporte un processus itératif de déformation graduelle o l'on combine linéairement à chaque itération, une réalisation géostatistique initiale (y) d'au moins une partie du milieu, et un nombre (N-1) (N>1) autres réalisations indépendantes de la réalisation initiale en imposant des contraintes aux coefficients de la combinaison linéaire, 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, et les dites données dynamiques, le processus itératif étant répété jusqu'à obtenir une réalisation optimale du modèle stochastique.

Dans le but de sélectionner les réalisations géostatistiques de façon à minimiser rapidement la dite fonction objectif et obtenir une convergence plus rapide du modèle, on génère un certain nombre N#, beaucoup plus grand que le nombre N, de réalisations et on sélectionne parmi elles, les réalisations ayant des indicateurs de plus forte valeur absolue jusqu'à avoir au total N cartes, l'indicateur associé à chaque réalisation correspondant à un produit scalaire entre cette réalisation et le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique initiale considérée.

Suivant un mode de mise en oeuvre, en plus de la réalisation géostatistique initiale (y) choisie arbitrairement, on sélectionne les (N1) autres réalisations à partir des (N-1) indicateurs de plus forte valeur absolue. Dans ce cas, on considère le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique initiale (y) conditionnée aux données.

Suivant un mode de mise en oeuvre, on sélectionne les N réalisations à partir des N indicateurs de plus forte valeur absolue. Dans, ce cas, on considère le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique de la variable aléatoire nulle conditionnée aux données statiques de puits.

On détermine les gradients géostatistiques par exemple par la méthode de l'état adjoint.

PRESENTATION SUCCINTE DES FIGURES

Les caractéristiques et avantages de la méthode selon l'invention, apparaîtront plus clairement à 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: - la figure 1 montre les principales étapes du problème direct lorsque l'on utilise les déformations graduelles comme paramétrisation du modèle géologique; - la figure 2 est une représentation géométrique de la méthode permettant de trouver le meilleure point d'initialisation pour l'algorithme d'optimisation dans le cas des indicateurs de raffinement pour les déformations graduelles; - la figure 3 montre les valeurs absolues des indicateurs de raffinement pour les déformations graduelles de 700 réalisations géostatistiques; - la figure 4 montre l'évolution de la fonction objectif J pour quatre boucles d'optimisation, les réalisations géostatistiques utilisées pour les déformations graduelles étant choisies ou non par les indicateurs de raffinement pour les déformations graduelles; - la figure 5 montre les valeurs absolues des indicateurs d'initialisation pour les déformations graduelles de 700 réalisations géostatistiques; et - la figure 6 montre l'évolution de la fonction objectif J pour quatre boucles d'optimisation, les réalisations géostatistiques utilisées pour les déformations graduelles étant choisies ou non par les indicateurs d'initialisation pour les déformations graduelles.

DESCRIPTION DETAILLEE DE LA METHODE

On considère que l'on effectue un calage d'historique de production sur un modèle de réservoir (à l'échelle géologique ou de simulation d'écoulement) paramétré à l'aide de la méthode des déformations graduelles. Cette méthode consiste à optimiser de manière itérative une combinaison linéaire de réalisations indépendantes du modèle stochastique jusqu'à ce que les contraintes dynamiques soient respectées (via une fonction objectif).

Les propriétés pétrophysiques (telles que la perméabilité ou la porosité) du modèle de réservoir sont modélisées par une fonction aléatoire Z. Considérons spécifiquement une fonction aléatoire gaussienne Z(x) centrée réduite. La méthode des déformations graduelles consiste à écrire une nouvelle réalisation z de Z comme une combinaison linéaire de N réalisations indépendantes zi de Z:

N N

z(p) =p,z, avec pi2 = 1 (1) i=1 i=1 La contrainte de normalité est automatiquement vérifiée lorsque l'on travaille en coordonnées sphériques {Io..... N- }. La nouvelle réalisation z est dès lors fonction de (N1) paramètres de déformation graduelle indépendants {0J, ...,ON-1} calculés par une relation: p=S(O) (2) Le concept des indicateurs de raffinement a été appliqué à la paramétrisation par déformations graduelles. Dans la pratique, un nombre N (généralement petit) de réalisations zi e Rm sont choisies au hasard, ou nm correspond au nombre de mailles géostatistiques (généralement grand). L'algorithme d'optimisation est initialiser avec des coordonnées sphériques 0 =O 0. Cela revient à initialiser les coefficients de la combinaison linéaire [Eq. (1)] à pl = 1, P2 = ... = PN = O0. L'utilisation des indicateurs de raffinement va permettre à l'utilisateur de réduire à une ou même zéro le nombre de réalisations géostatistiques choisies au hasard.

Cas 1: Choix au hasard d'une réalisation géostatistique zi Une fois que la réalisation initiale zl a été choisie, l'utilisateur choisit les (N-1) réalisations z2...zN utilisées dans les déformations graduelles parmi un jeu de N# (N# " N) réalisations Z2...ZN. L'idée est de générer les réalisations Z2...ZN# à partir d'un grand nombre de germes aléatoires et de ne prendre en compte pour la déformation graduelle que les (N-1) réalisations ayant les plus forts indicateurs de raffinement 2i,. que nous définissons maintenant.

Soit p = (p '...pN# ), et considérons le problème d'optimisation sous contraintes suivant: Trouverp e RN qui minimise J(z) avec z = pii i=1 sous les contraintes: (3) N# [A] pi = bi, i = 2...N, [B] p2 = 1 i=I o bi, i= 2...N# sont des nombres spécifiés de telle manière que 2b2 < 1. Toute solution locale p* de (3) satisfait la condition nécessaire de Lagrange associée il existe A *= (2.. .2.;) (multiplicateurs pour la contrainte [A]) et u* (multiplicateur pour la contrainte [B]) tels que OL *,2',/2') = 0 4 -(p,i "u)=O (4) ap o L est le Lagrangien définit par: L(p,c,#)= (Zi=lPizi)+ Ei (pi -bi)2i +(,# 152 1 i=2 i=i i (5) IVP =(pi..PNEU ,A #-i)' VI La formule (4) nous donne immédiatement les multiplicateurs de Lagrange: )i Z* = Z RÀ = (-,(z* ,zi i = 2...N#, avec z*= PiZi (6) aDz Rn i=1 La gradient géostatistique aJ/az(z* )e Rm correspond à la dérivée de la fonction objectif J par rapport à chaque cellule géostatistique de la réalisation z*. La solution p* de (3), - de même que la réalisation géostatistique associée z* et les multiplicateurs de Lagrange A*, ,u* dépendent du second membre b = (b2...bN. ) de la contrainte [1], de sorte que nous pouvons les noter par Pb, z*, < and /*. La valeur minimale (optimale) de la fonction objectif associée au second membre b est donc J = (Zb) (7) Un résultat bien connu d'optimisation sous contraintes nous dit que le multiplicateur de Lagrange i coïncide avec la dérivée de la fonction objectif optimale Jb par rapport au iième élément bi du second membre de la contrainte: b s- - (8) abi Pour appliquer ce résultat à notre problème, nous remarquons que pour le choix b2 = b3 =... = bN# = 0, l'espace des solutions du problème (4) contient uniquement deux points isolés p = ( 1,0,....0)e RN qui sont donc des solutions locales de (4)! Nous pouvons donc appliquer l'analyse précédente avec p* = (1,0,...,0) et z* = zl. Si nous notons J(z ) la valeur minimale de la fonction objectif J lorsque le i-ième élément du second membre b passe de bi = 0 à bi = bi, nous voyons qu'au premier ordre: J(z )-J(z*)=-i;Sbi (9) La i-ième composante 2R du multiplicateur de Lagrange 2 nous indique donc la sensibilité de la fonction objectif optimale Ja lorsque l'on prend en compte le i-ième degré de liberté c'est-à-dire que l'on utilise la i-ième réalisations zi pour la paramétrisation des déformations graduelles. Nous appellerons dorénavant ces multiplicateurs de Lagrange des indicateurs de raffinement pour les déformations graduelles. Afin de sélectionner, parmi les N# réalisations candidates, les N-1 qui doivent être associée à z1 pour utiliser les déformations graduelles, nous calculons les N# -1 indicateurs de raffinement 22...2N# par la formule 6. Ceci ce fait de manière extrêmement rapide, étant donné que chaque 2R correspond à un simple produit scalaire, une fois que le gradient géostatistique aJ/Dz(z* ) Rm a été calculé (cf. paragraphe suivant). Les indicateurs de raffinement sont alors triés par ordre décroissant selon leur valeur absolue et nous sélectionnons les (N-1) réalisations géostatistiques correspondant aux (N-l) indicateurs de raffinement de plus forte valeur absolue.

Cas 2: Eviter le choix au hasard de la réalisation géostatistique zl Supposons que les propriétés pétrophysiques (telles que la perméabilité ou la porosité) du modèle de réservoir soient modélisées par une distribution lognormale Y(x) de moyenne m et de variance 02. Cette réalisation Y(x) est liée à une distribution normale U(x) de moyenne m' et de variance a'2 à travers la relation: Y(x)=eU(x) (10) Les déformations graduelles utilisent une variable aléatoire centrée réduite Z(x). L'équation 10 peut alors être réécrite comme: Y(x) = e(m"+Z(x)) (11) Au lieu de calculer le gradient géostatistique aJ/az pour une réalisation z1 choisie au hasard (cas 1), l'utilisateur évalue le gradient aJ/az pour la variable aléatoire nulle Z - O conditionnée aux données statiques de puits. La motivation vient du fait que si une infinité de réalisations étaient prises en compte dans la combinaison linéaire [Eq. (1)], la réalisation résultante serait égale à la moyenne em c'est-à-dire l'équation (11) avec Z - 0. Ce gradient nous donne la sensibilité de la fonction objectif pour une distribution déterministe égale à la moyenne de la distribution conditionnée aux données statiques de puits.

Afin de sélectionner les N réalisations utilisées dans les déformations graduelles, l'utilisateur génère comme dans le cas 1 un grand nombre N# de réalisations géostatistiques et calcule les indicateurs d'initialisation pour les déformations graduelles: Ai = (-o),Zi) i 1...N# (12) az R-' Par définition, nous avons: J(SPilZi) - J(O) ayn aiepi, i = 1.

N# (13) Les réalisations zi ayant une valeur absolue lAil importante peuvent potentiellement faire décroître de manière forte la fonction objectif pour un Spi de signe correct. L'utilisateur trie par ordre décroissant de valeur absolue les indicateurs d'initialisation et choisit pour la déformation graduelle les N réalisations ayant la plus forte valeur absolue. Calcul du gradient géostatistique aJ/Dz L'étape préliminaire au calcul des indicateurs de raffinement [Eq. (6)] consiste à calculer le gradient géostatistique aJ/az. A cette fin, considérons tout d'abord les différentes étapes du problème direct lorsque les déformations graduelles sont utilisées comme paramétrisation du modèle géologique (Fig. 1). Les quatre étapes successives sont: [1] Déformation graduelle des N réalisations zi qui résulte dans la réalisation z, [2] Modélisation géologique: cas de distributions lognormales ou modèles en faciès. Le conditionnement aux données statiques de puits est effectué. [3] Changement d'échelle pour passer du modèle géologique au modèle de simulation d'écoulement si nécessaire, [4] Simulation d'écoulement et calcul de la fonction objectif J. Par dérivation composée, on peut écrire le gradient DJ/Dz comme: aD J a aKk (14) az aKak az i 2 3 Un des points essentiels de la méthodologie proposée est de calculer ces gradients par la méthode de l'état adjoint. Le calcul du gradient aJ/aK se fait par un état adjoint discret. Le second terme peut lui aussi être calculé par état adjoint si nécessaire. Le troisième terme correspond à la dérivation de l'étape de modélisation géologique et peut être facilement calculé de manière analytique...DTD: Avec cette approche, le coût supplémentaire induit par le calcul du gradient aJ/aK est similaire au coût informatique requis pour le calcul de la fonction objectif J - et, plus important, est indépendant du nombre de cellules géostatistiques, qui est très élevé.

Utilisation des indicateurs pour initialiser l'algorithme d'optimisation L'équation (9) nous indique que le signe des indicateurs de raffinement contient une information importante. Supposons qu'un indicateur donné ait un signe positif. Si l'on fixe à la réalisation associée un poids positif, cela aura tendance à faire décroître la fonction objectif (au premier ordre). La même analyse tient pour des indicateurs négatifs. L'utilisateur aura donc tout intérêt à initialiser l'algorithme d'optimisation avec des 0 donnant des coefficients de combinaison linéaire de même signe que l'indicateur sous considération.

Cas 1: Choix au hasard d'une réalisation géostatistique zl Une fois que la réalisation géostatistique zl a été choisie, l'utilisateur détermine les (N-1) réalisations z2... ZN utilisées pour les déformations graduelles en se basant sur les indicateurs de raffinement associés. L'utilisateur a calculé VZJ(zl) et, par simple produit scalaire, les composantes 2...AN de VPJ((1,O...O)). Il peut donc chercher sur la sphère E iPi2 = 1 un nouveau point initial dans la direction - VpJ: - Si Ri > O, J tend à décroître lorsque pi est augmenté. Dans ce cas, p = (1,0...0) est le meilleur point d'initialisation, - Si) < 0, l'utilisateur peut se déplacer dans la direction -7VpJ jusqu'à intersecter à nouveau la sphère (fig. 2): Pl= lEI, P2 = - Eú2.... PN = N ú 2 2 zzO&1,o**N =0 - 2N 2 2q_(15) P1 + + .3 = 1 N + + A Cas 2 Eviter le choix au hasard de la réalisation géostatistique zl L'utilisateur a calculé VzJ(z =O), il a donc à sa disposition VpJ(p= O)=(Al...AN). Il peut initialiser l'algorithme d'optimisation en suivant les règles suivantes {pij = O -.-A, j- 1 =...N e1 (16) 2: ±-. + p:12+ 2)2 Pl N (A +... + AN3 13 2852710 La transformation inverse S-l de l'équation (2) donne les paramètres de déformation graduelle {O....N-1} correspondants.

Extension de la méthodologie aux cas de modèles en faciès Les modèles en faciès correspondent à des modèles présentant des discontinuités au niveau de grandeurs physiques telles que la perméabilité par exemple, ce qui rend le terme [3] de l'équation (14) non dérivable. Une méthodologie proposée par: - Schaaf, T., Mezghani, M. and Chavent, G., 2002, "Direct conditioning of fine-scale facies models to dynamic data by combining graduai déformation and numerical upscaling techniques", Paper E-44: Proc. 8th European Conference on Mathematics of Oil Recovery (ECMOR VIII), 3-6 September 2002, Freiberg, Germany.

introduit le concept de faciès de transition qui permet d'avoir une fonction k = f(z) dérivable et donc de calculer le terme [3] de l'équation (14). La méthode selon l'invention s'applique donc à des modèles en faciès dès lors que l'on utilise le concept des faciès de transition décrit dans cette publication pour calculer le terme [3] de l'équation (14).

EXEMPLE

On décrit ici un exemple synthétique illustrant la méthode de façon nonlimitative.

On considère un réservoir avec un maillage cartésien contenant suivant les axes X, Y et Z respectivement 141, 141 et 4 mailles. La seule propriété pétrophysique considérée est la perméabilité. La perméabilité est modélisée par une distribution lognormale de moyenne égale à 100 mD et de variance égale à (100)2nD'. La porosité est constante dans le réservoir. On considère un cas synthétique de tests d'interférence, il s'agit donc d'une simulation d'écoulement monophasique. Les seules données issues du simulateur d'écoulement sont les pressions puits. Le réservoir est traversé par cinq puits: un puits producteur au centre du modèle et quatre puits observateurs disposés en croix autour du puits producteur.

Une réalisation générée à partir d'un germe aléatoire sert de modèle de référence. Une simulation d'écoulement de fluides est conduite sur ce modèle, ce qui nous donne des pressions puits de référence. La fonction objectif est formulée au sens des moindres 14 2852710 carrées. On ne considère pas d'étape de changement d'échelle. La simulation est directement conduite sur le modèle considéré. On cherche à résoudre notre problème inverse en paramétrant notre modèle par deux paramètres de déformation graduelle 0 c'est-à-dire que l'on considère une combinaison linéaire de trois réalisations géostatistiques zi.

Pour cette exemple synthétique, nous allons tester les deux approches décrites dans les cas 1 et 2.

Cas 1: Choix au hasard d'une réalisation géostatistique z1 On choisit une réalisation géostatistique z1 au hasard pour laquelle on calcule le gradient aJ/az. On génère par la suite un jeu de 700 réalisations géostatistiques pour lesquelles on calcule les indicateurs de raffinement correspondants (Fig. 3) Afin de démontrer l'utilité de la démarche proposée, nous allons considérer les quatre boucles d'optimisation suivantes (Fig. 4): - la boucle 1 correspond à l'évolution de la fonction objectif J pour une optimisation conduite avec la réalisation géostatistique z1 et deux autres réalisations choisies au hasard. L'algorithme d'optimisation est initialisé avec 0 = 0, - la boucle 2 correspond à l'optimisation conduite avec la réalisation z1 et les réalisations ayant les indicateurs de raffinement de plus forte valeur absolue dans le jeu considéré.

L'algorithme d'optimisation est initialisé avec O = 0, - la boucle 3 est similaire à la boucle 2 en ce qui concerne les réalisations choisies.

L'algorithme d'optimisation est initialiser en utilisant les valeurs des indicateurs (cf. paragraphe précédent), - la boucle 4 correspond à l'optimisation conduite avec la réalisation z1 et les réalisations ayant les indicateurs de raffinement de plus faible valeur absolue dans le jeu considéré.

L'algorithme d'optimisation est initialisé avec O = 0.

On vérifie que dans les cas 2 et 3 o l'on applique la méthode, la fonction coût décroît très rapidement.

2852710 L'utilisation des indicateurs de raffinement non seulement pour le choix à priori des réalisations géostatistiques mais aussi pour l'initialisation de l'algorithme d'optimisation apparaît comme optimal. La valeur finale de la fonction objectif dans le cas ou l'on utilise les indicateurs de raffinement pour les déformation graduelles sera inférieure ou égale au cas ou ces indicateurs ne sont pas utilisés.

Cas 2: Eviter le choix au hasard de la réalisation géostatistique z1 On choisit la réalisation z de la variable aléatoire Z nulle (sur laquelle est effectué le conditionnement aux données statiques de puits) pour laquelle on calcule le gradient aJ/az. On génère par la suite un jeu de 700 réalisations géostatistiques pour lesquelles on calcule les indicateurs d'initialisation correspondants (Fig. 5). Afin de démontrer l'utilité de la démarche proposée, nous allons considérer les quatreboucles d'optimisation suivantes (Fig. 6): - la boucle 1 correspond à l'évolution de la fonction objectif J pour une optimisation conduite avec la réalisation géostatistique z1 et deux autres réalisations choisies au hasard. L'algorithme d'optimisation est initialisé avec 0 = 0, - la boucle 2 correspond à l'optimisation conduite avec la réalisation z1 et les réalisations ayant les indicateurs d'initialisation de plus forte valeur absolue dans le jeu considéré.

L'algorithme d'optimisation est initialisé avec 0 = 0, - la boucle 3 est similaire à la boucle 2 en ce qui concerne les réalisations choisies.

L'algorithme d'optimisation est initialisé en utilisant les valeurs des indicateurs (cf. paragraphe précédent), - la boucle 4 correspond à l'optimisation conduite avec la réalisation, z1 et les réalisations ayant les indicateurs d'initialisation de plus faible valeur absolue dans le jeu considéré.

L'algorithme d'optimisation est initialisé avec 0 = 0.

On vérifie que dans les cas 2 et 3 o l'on applique la méthode, la fonction coût décroît aussi très rapidement.

16 2852710

Claims (6)

REVENDICATIONS

1) Méthode pour former un modèle numérique stochastique de type Gaussien ou apparenté, donnant une image de la distribution d'une grandeur physique dans un milieu hétérogène poreux, qui soit calé par rapport à des données obtenues par des mesures effectuées dans le milieu ou des observations préalables, et caractéristiques du déplacement des fluides dans le milieu, comportant un processus itératif de déformation graduelle o l'on combine linéairement à chaque itération, une réalisation initiale (y) et un nombre (N-1) (N>1) autres réalisations (yi) (avec i = 1,..., (N - i)), indépendantes de la réalisation initiale, en imposant des contraintes aux coefficients de la combinaison linéaire, 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, et les dites données dynamiques, le processus itératif étant répété jusqu'à obtenir une réalisation optimale du modèle stochastique, caractérisée en ce que, dans le but de sélectionner les réalisations géostatistiques permettant de minimiser rapidement la dite fonction objectif et d'obtenir une convergence plus rapide du modèle: - on génère un certain nombre N# (N# " N) de réalisations et on sélectionne parmi elles, (N-1) autres réalisations ayant des indicateurs de plus forte valeur absolue, l'indicateur associé à chaque réalisation correspondant à un produit scalaire entre elle et le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique initiale (y) laquelle est choisie arbitrairement et conditionnée aux dites données.

2) Méthode selon la revendication 1, caractérisée en ce que l'on sélectionne les (N-1) autres réalisations ayant des indicateurs de plus forte valeur absolue, l'indicateur associé à chaque réalisation correspondant à un produit scalaire entre elle et le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique de la variable aléatoire nulle conditionnée aux dites données, la réalisation géostatistique initiale (y) correspondant à l'indicateur de plus forte valeur absolue, et les autres réalisations (yi) correspondant aux indicateurs suivants de plus forte valeur absolue par ordre décroissant.

3) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on utilise les indicateurs pour initialiser un algorithme d'optimisation avec des paramètres initiaux 0 de déformation graduelle optimaux et non nuls.

17 2852710

4) Méthode selon l'une des revendications précédentes, caractérisée en ce que l'on détermine le gradient géostatistique par la méthode de l'état adjoint.

5) Méthode selon l'une des revendications précédentes, 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 le gradient géostatistique i.e. la dérivée de la fonction objectif par rapport aux composantes du bruit blanc gaussien.

6) Application de la méthode aux modèles en faciès en rendant continu le gradient géostatistique dès lors qu'on rajoute des faciès de transition aux facies réel du milieu.