CA2519184C - 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

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 Download PDF

Info

Publication number
CA2519184C
CA2519184C CA2519184A CA2519184A CA2519184C CA 2519184 C CA2519184 C CA 2519184C CA 2519184 A CA2519184 A CA 2519184A CA 2519184 A CA2519184 A CA 2519184A CA 2519184 C CA2519184 C CA 2519184C
Authority
CA
Canada
Prior art keywords
realization
realizations
geostatistical
initial
objective function
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CA2519184A
Other languages
English (en)
Other versions
CA2519184A1 (fr
Inventor
Thomas Schaaf
Guy Chavent
Mokhles Mezghani
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
IFP Energies Nouvelles IFPEN
Original Assignee
IFP Energies Nouvelles IFPEN
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by IFP Energies Nouvelles IFPEN filed Critical IFP Energies Nouvelles IFPEN
Publication of CA2519184A1 publication Critical patent/CA2519184A1/fr
Application granted granted Critical
Publication of CA2519184C publication Critical patent/CA2519184C/fr
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/17Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

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 (.lambda.), .lambda. .isin. 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 .lambda.i de la i-ième composante du multiplicateur .lambda. indique la sensibilité de la fonction objectif par rapport à la i-ième réalisation géostatistique (zi).

Description

METHODE POUR FORMER RAPIDEMENT UN MODELE STOCHASTIQUE
REPRESENTATIF DE LA DISTRIBUTION D'UNE GRANDEUR PHYSIQUE
DANS UN MILIEU HETEROGENE PAR UNE SELECTION APPROPRIEE DE

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 pet mé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 pelmé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 4e1 European Conference on the mathematical of ou l recovery (ECMOR IV), 1994.
Elle permet d'effectuer un calage d'historique 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 petmettent de corriger cet inconvénient. Le concept d'indicateur de raffinement, présenté par:
- Chavent, G., 8c, Bissel!, R., "Indicator for the refinement of para.meterization", 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éfoimation 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 paranié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 [I] 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 par.amè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 FR 2 846 767 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 à
5 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 selection of stochastic model realizations constrained to historical data", SPE 38731, 1997.
LA METHODE SELON L'INVENTION
L'invention vise une méthode pour tester au moins un scénario de production d'un milieu poreux hétérogène, par formation d'un modèle numérique stochastique de type Gaussien, donnant une image de la distribution d'une grandeur physique dans ledit milieu hétérogène poreux, incluant la sélection d'une réalisation géostatique initiale dudit modèle stochastique de type Gaussen, et modification de la réalisation géostatique initiale, 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, par un processus itératif de déformation graduelle de la réalisation géostatique initiale, caractérisé en ce qu'on réalise les étapes suivantes:
- on génère un nombre N# de réalisations indépendantes de la réalisation géostatique initiale;
- on associe à chacune des N# réalisations un indicateur de raffinement indiquant la sensibilité de ladite fonction objectif à l'utilisation de la réalisation associée dans ledit processus itératif de déformation graduelle et l'indicateur de raffinement étant un produit scalaire entre ladite réalisation associée et le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique initiale (y) , dans lequel l'on détermine le gradient géostatistique par la méthode de l'état adjoint;
- on sélectionne parmi lesdites N# réalisations, n réalisations correspondant aux n indicateurs de raffinement de plus forte valeur absolue;

5a - on combine linéairement à chaque itération, une réalisation géostatistique initiale (y) et au moins une autre réalisation (yi), indépendante de la réalisation initiale, en imposant des contraintes aux coefficients de la combinaison linéaire;
- on calcule, au moyen d'un simulateur, les données dynamiques déduites de la combinaison linéaire pour chaque itération et on calcule et on minimise une fonction objectif (J) mesurant l'écart entre le jeu de données calculées;
- on optimise par itération la combinaison linéaire pour chaque itération jusqu'à
obtenir une réalisation optimale du modèle stochastique de type Gaussien, correspondant à une valeur minimale de la fonction objectif; et - on teste le scénario de prodùction du milieu poreux hétérogène par utilisation de la réalisation optimale.
De préférence, la méthode 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 (yj) (avec i= 1, ..., (N-1)), 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 ladite combinaison au moyen d'un simulateur, et lesdites données dynamiques, le processus itératif étant répété jusqu'à obtenir une réalisation optimale du modèle stochastique.

De préférence, dans le but de sélectionner les réalisations géostatistiques de façon à minimiser rapidement ladite fonction objectif et obtenir une convergence plus 5b 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, (N-1) autres réalisation ayant des indicateurs de plus forte valeur absolue, 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 (y) laquelle est choisie arbitrairement et conditionnée auxdites données.
Suivant un mode de mise en oeuvre, en plus de la réalisation géostatistique initiale (y) choisie arbitrairement, on sélectionne les (N-1) 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 peiniettant 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:

z(p)= avecIp =1 1=1 (1) La contrainte de normalité est automatiquement vérifiée lorsque l'on travaille en coordonnées sphériques {0/,...,ON_1 }. La nouvelle réalisation z est dès lors fonction de (N-1) paramètres de déformation graduelle indépendants {0,,...,ON_1 calculés par une relation :

p = S(0) (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 lem 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. Cela revient à initialiser les coefficients de la combinaison linéaire [Eq. (1)] à pl =- 1, p, = ...= p N = O. 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 z1 Une fois que la réalisation initiale z1 a été choisie, l'utilisateur choisit les (N-1) réalisations z,...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 21 , que nous définissons maintenant.

Soit p =(pi,...,p N#), et considérons le problème d'optimisation sous contraintes suivant :

Trouver p* E RN# qui minimise J(z) avec z =1 p izi N#

, sous les contraintes :

(3) [A] pi = bi, i =2...N# , [B] E P
12 =1N#
-où b 1 , i = 2...N# sont des nombres spécifiés de telle manière que EN# b2 < 1 . Toute i=2 i solution locale p* de (3) satisfait la condition nécessaire de Lagrange associée : il existe .1.* = (il.....1; 1,, ) (multiplicateurs pour la contrainte [A]) et ,u*
(multiplicateur pour la contrainte [B]) tels que:

ap (4) où L est le Lagrangien définit par:

I L(p,i1)= J Epizi +(pi ¨ bi pli \ ( Ao 2 V p = (pi...p Aii,), Vii, --= (22...2N, ), V,u \
(5) La formule (4) nous donne immédiatement les multiplicateurs de Lagrange :

2: = -( f \ al ¨v* ),z, i=2...N# , avec z*
=I p:: zi N#
(6) La gradient géostatistique ai 1 az(z.')E Rnm 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 f, If - dépendent du second membre b = (b2...b N,) de la contrainte [1], de sorte que nous pouvons les noter par pb* , zb* , .1.: and 1.2b* . La valeur minimale (optimale) de la fonction objectif associée au second membre b est donc :

`1"1: J(Zb*) (7) Un résultat bien connu d'optimisation sous contraintes nous dit que le multiplicateur de Lagrange .1; coïncide avec la dérivée de la fonction objectif optimale Jb* par rapport au i-ième élément bi du second membre de la contrainte :

ai*
abi=-/Ç (8) Pour appliquer ce résultat à notre problème, nous remarquons que pour le choix b2 = b3 = bN# = O, l'espace des solutions du problème (4) contient uniquement deux points isolés p en qui sont donc des solutions locales de (4)! Nous pouvons donc appliquer l'analyse précédente avec p* = (1,0,...,0) et z* z1. 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 = cSbz, nous voyons qu'au premier ordre :

J(z ¨ (z* = ¨2; Mi (9) La i-ième composante il, du multiplicateur de Lagrange il nous indique donc la sensibilité

de la fonction objectif optimale Jb* 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 à z/ 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 2 correspond à un simple pioduit scalaire, une fois que le gradient géostatistique a,1 laz(z*)E R"ni 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-1) indicateurs de raffinement de plus forte valeur absolue.

Cas 2 : Eviter le choix au hasard de la réalisation géostatistique z1 Supposons que les propriétés pétrophysiques (telles que la perméabilité bu la porosité) du modèle de réservoir soient modélisées par une distribution lognormale 17(x) de moyenne ni et de variance c72. Cette réalisation 114 est liée à une distribution normale U(x) de moyenne in' et de variance (y/2 à travers la relation :

11.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 :

17(4= e('x)) (11) Au lieu de calculer le gradient géostatistique a//az pour une réalisation z1 choisie au hasard (cas 1), l'utilisateur évalue le gradient aflaz pour la variable aléatoire nulle Z 0 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éfonnations 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 :

A, =(¨(0),z1) , =1...N# az R"

(12) Par définition, nous avons :

J(c5p1z1)¨ AlOpi, i =1...N#

(13) Les réalisations zi ayant une valeur absolue A importante peuvent potentiellement faire décroître de manière forte la fonction objectif pour un 8p1 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 ailaz L'étape préliminaire au calcul des indicateurs de raffinement [Eq. (6)]
consiste à calculer le gradient géostatistique aJiaz . 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 af/az comme :

aJ DK a (14) az aK ak aZ

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 ai yaK 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.

Avec cette approche, le coût supplémentaire induit par le calcul du gradient ailaK 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 6' 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 zi Une fois que la réalisation géostatistique z1 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é V zJ(z1) et, par simple produit scalaire, les composantes 21...2N de V pf((1,0...0)) . Il peut donc chercher sur la sphère µ71N 2 p, =1 un nouveau point initial dans la direction - Vp J:

- Si 21 > 0, J tend à décroître lorsque p1 est augmenté. Dans ce cas, p=
(1,0...0) est le meilleur point d'initialisation, - Si 21 < 0 , l'utilisateur peut se déplacer dans la direction - Vp./.
jusqu'à intersecter à

nouveau la sphère (fig. 2):

{p1 = 1- al, P2 = - g/12, - =, Pi,/ = -"e2N
8 = 221 (15) Pi2 + ===+ PN2 =1 22 1 + ...- F r N

Cas 2 : Eviter le choix au hasard de la réalisation géostatistique z1 L'utilisateur a calculé VzJ(z = 0), il a donc à sa disposition Vp J(p = 0) =
(A/ ...AN ). Il peut initialiser l'algorithme d'optimisation en suivant les règles suivantes :

{ p i = 0 - eA i, j =1...N

pi2 + ...+ pN2 =1 8 = (42 + . .
. + A N2 )1/2 (16) La transformation inverse S' de l'équation (2) donne les paramètres de déformation graduelle {0 1,...,0 N_I} 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 gradual defoimation 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 peimet 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 non-limitative.
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)2 m12. 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 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 8 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 zi On choisit une réalisation géostatistique z1 au hasard pour laquelle on calcule le gradient ailaz . 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 zi 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 0 = 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 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 très rapidement.

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 zi 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 aJiaz . 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 quatre boucles 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.

Claims (6)

1. Méthode pour tester au moins un scénario de production d'un milieu poreux hétérogène, par formation d'un modèle numérique stochastique de type Gaussien, donnant une image de la distribution d'une grandeur physique dans ledit milieu hétérogène poreux, incluant la sélection d'une réalisation géostatique initiale dudit modèle stochastique de type Gaussien, et modification de la réalisation géostatique initiale, 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, par un processus itératif de déformation graduelle de la réalisation géostatique initiale, caractérisé en ce qu'on réalise les étapes suivantes:
- on génère un nombre N# de réalisations indépendantes de la réalisation géostatique initiale;
- on associe à chacune des N# réalisations un indicateur de raffinement indiquant la sensibilité de ladite fonction objectif à l'utilisation de la réalisation associée dans ledit processus itératif de déformation graduelle et l'indicateur de raffinement étant un produit scalaire entre ladite réalisation associée et le gradient géostatistique de la fonction objectif par rapport à la réalisation géostatistique initiale (y), dans lequel l'on détermine le gradient géostatistique par la méthode de l'état adjoint;
- on sélectionne parmi lesdites N# réalisations, n réalisations correspondant aux n indicateurs de raffinement de plus forte valeur absolue;
- on combine linéairement à chaque itération, une réalisation géostatistique initiale (y) et au moins une autre réalisation (y i), indépendante de la réalisation initiale, en imposant des contraintes aux coefficients de la combinaison linéaire;
- on calcule, au moyen d'un simulateur, les données dynamiques déduites de la combinaison linéaire pour chaque itération et on calcule et on minimise une fonction objectif (J) mesurant l'écart entre le jeu de données calculées;

- on optimise par itération la combinaison linéaire pour chaque itération jusqu'à
obtenir une réalisation optimale du modèle stochastique de type Gaussien, correspondant à une valeur minimale de la fonction objectif; et - on teste le scénario de production du milieu poreux hétérogène par utilisation de la réalisation optimale.
2. Méthode selon la revendication 1, dans laquelle la réalisation géostatique initiale est choisie arbitrairement et conditionnée auxdites données.
3. Méthode selon la revendication 1, dans laquelle ladite réalisation géostatistique initiale (y) est la réalisation, parmi lesdites n réalisations, dont l'indicateur a la plus forte valeur absolue, les (n-1) autres réalisations (yj) correspondant aux indicateurs de raffinement suivants de plus forte valeur absolue.
4. Méthode selon l'une quelconque des revendications 1 à 3, caractérisée en ce que l'on utilise les indicateurs de raffinement pour initialiser un algorithme d'optimisation avec des paramètres initiaux 0 de déformation graduelle optimaux et non nuls.
5. Méthode selon l'une quelconque des revendications 1 à 4, caractérisée en ce que ladite 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. Méthode selon l'une quelconque des revendications 1 à 5, dans laquelle ladite grandeur physique est discontinue, et on calcule ledit gradient géostatistique de la fonction objectif en utilisant une technique de faciès de transition.
CA2519184A 2003-03-18 2004-03-16 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 Expired - Fee Related CA2519184C (fr)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
FR03/03276 2003-03-18
FR0303276A FR2852710B1 (fr) 2003-03-18 2003-03-18 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
PCT/FR2004/000644 WO2004086280A2 (fr) 2003-03-18 2004-03-16 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

Publications (2)

Publication Number Publication Date
CA2519184A1 CA2519184A1 (fr) 2004-10-07
CA2519184C true CA2519184C (fr) 2013-05-14

Family

ID=32922242

Family Applications (1)

Application Number Title Priority Date Filing Date
CA2519184A Expired - Fee Related CA2519184C (fr) 2003-03-18 2004-03-16 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

Country Status (6)

Country Link
US (1) US7558715B2 (fr)
EP (1) EP1606727A2 (fr)
CA (1) CA2519184C (fr)
FR (1) FR2852710B1 (fr)
NO (1) NO20054259L (fr)
WO (1) WO2004086280A2 (fr)

Families Citing this family (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2869116B1 (fr) * 2004-04-14 2006-06-09 Inst Francais Du Petrole Methode pour construire un modele geomecanique d'une zone souterraine destine a etre couple a un modele de reservoir
FR2869421B1 (fr) * 2004-04-27 2006-06-02 Inst Francais Du Petrole Methode de reconstruction d'un modele stochastique, representatif d'un milieu heterogene poreux, pour ameliorer son calage par les donnees de production
US7805248B2 (en) * 2007-04-19 2010-09-28 Baker Hughes Incorporated System and method for water breakthrough detection and intervention in a production well
FR2919932B1 (fr) * 2007-08-06 2009-12-04 Inst Francais Du Petrole Methode pour evaluer un schema de production d'un gissement souterrain en tenant compte des incertitudes
WO2009145960A1 (fr) 2008-04-17 2009-12-03 Exxonmobil Upstream Research Company Outil de support de decision basee sur une optimisation robuste, utilise dans la planification de developpement de reservoir
EP2291761A4 (fr) * 2008-04-18 2013-01-16 Exxonmobil Upstream Res Co Outil d'aide à la décision à base de processus de décision de markov pour une planification de développement de réservoir
BRPI0910333A2 (pt) 2008-04-21 2015-10-06 Exxonmobil Upstream Res Co métodos para planejamento de desenvolvimento de um reservatório, para suporte à decisão que considera o desenvolvimento de recursos de petróleo, para otimização do planejamento de desenvolvimento com base em computador para um reservatório de hidrocarboneto, e para produzir hidrocarbonetos, e, produto de programa de computador
AU2009341851B2 (en) 2009-03-11 2015-07-16 Exxonmobil Upstream Research Company Gradient-based workflows for conditioning of process-based geologic models
AU2009341852B2 (en) 2009-03-11 2015-08-20 Exxonmobil Upstream Research Company Adjoint-based conditioning of process-based geologic models
EP2499567A4 (fr) 2009-11-12 2017-09-06 Exxonmobil Upstream Research Company Procédé et appareil de modélisation et de simulation de réservoirs
FR2953039B1 (fr) * 2009-11-26 2012-01-13 Inst Francais Du Petrole Methode d'exploitation d'un gisement petrolier par reconstruction de modele de reservoir
WO2012021292A1 (fr) 2010-08-09 2012-02-16 Conocophillips Company Procédé de mise à l'échelle supérieure de réservoir avec une transmissibilité préservée
BR112013006777B1 (pt) * 2010-09-27 2021-10-26 Total Sa Método de simulação de fenômenos de carstificação em uma região cárstica e dispositivo de simulação de fenômenos de carstificação em uma região cárstica
US8942966B2 (en) * 2010-10-20 2015-01-27 Conocophillips Company Method for parameterizing and morphing stochastic reservoir models
US9488047B2 (en) 2011-04-04 2016-11-08 Conocophillips Company Reservoir calibration parameterization method
CA2887319C (fr) * 2012-10-05 2019-09-17 Total Sa Procede pour determiner une region karstique
FR3004270B1 (fr) 2013-04-05 2015-05-01 Storengy Methode de determination d'un modele cale pour un reservoir souterrain de fluide

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5838634A (en) * 1996-04-04 1998-11-17 Exxon Production Research Company Method of generating 3-D geologic models incorporating geologic and geophysical constraints
US6549879B1 (en) * 1999-09-21 2003-04-15 Mobil Oil Corporation Determining optimal well locations from a 3D reservoir model
US6480790B1 (en) * 1999-10-29 2002-11-12 Exxonmobil Upstream Research Company Process for constructing three-dimensional geologic models having adjustable geologic interfaces
FR2821946B1 (fr) * 2001-03-07 2003-05-16 Inst Francais Du Petrole Methode pour deformer graduellement une repartition intiale d'objets dans un milieu heterogene, generee par simulation d'un modele stochastique de type objet, pour l'adapter au mieux a des contraintes physiques imposees
FR2823877B1 (fr) * 2001-04-19 2004-12-24 Inst Francais Du Petrole Methode pour contraindre par des donnees dynamiques de production un modele fin representatif de la repartition dans le gisement d'une grandeur physique caracteristique de la structure du sous-sol

Also Published As

Publication number Publication date
WO2004086280A2 (fr) 2004-10-07
WO2004086280A3 (fr) 2004-10-28
US7558715B2 (en) 2009-07-07
CA2519184A1 (fr) 2004-10-07
FR2852710B1 (fr) 2005-04-29
NO20054259D0 (no) 2005-09-15
FR2852710A1 (fr) 2004-09-24
EP1606727A2 (fr) 2005-12-21
NO20054259L (no) 2005-12-16
US20060241925A1 (en) 2006-10-26

Similar Documents

Publication Publication Date Title
CA2519184C (fr) 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
EP1760492B1 (fr) Méthode pour mettre à jour un modèle géologique de réservoir a l&#39;aide de données dynamiques
JP5797262B2 (ja) 成熟炭化水素産地をシミュレートするための生産シミュレータ
Schwenk et al. The life of a meander bend: Connecting shape and dynamics via analysis of a numerical model
van der Westhuysen Modeling of depth‐induced wave breaking under finite depth wave growth conditions
EP2400320B1 (fr) Procédé d&#39;exploitation d&#39;un gisement pétrolier à partir d&#39;un modèle de réservoir déformé graduellement au moyen de cosimulations
CA2445858C (fr) Methode pour former plus rapidement un modele stochastique representatif d&#39;un reservoir heterogene souterrain, contraint par des donnees dynamiques
EP2343576B1 (fr) Méthode d&#39;exploitation d&#39;un gisement pétrolier à partir d&#39;une construction d&#39;une carte de faciès représentative du gisement
Cumming et al. Multiwell deconvolution
CA2821099C (fr) Procede d&#39;exploitation d&#39;un reservoir geologique a partir d&#39;un modele de reservoir cale par le calcul d&#39;une loi analytique de distribution conditionnelle de parametres incertains du modele
FR3023641A1 (fr)
WO2009027598A1 (fr) Procede, programme et systeme informatique de conciliation de donnees de modele de reservoir d&#39;hydrocarbure
EP2963235A1 (fr) Procede d&#39;exploitation d&#39;un gisement petrolier a partir d&#39;une technique de positionnement des puits a forer
de Graaf et al. Hyper‐resolution continental‐scale 3‐D aquifer parameterization for groundwater modeling
FR3034894A1 (fr)
Rwechungura et al. Application of particle swarm optimization for parameter estimation integrating production and time lapse seismic data
EP2770162B1 (fr) Procédé d&#39;exploitation d&#39;un réservoir géologique au moyen d&#39;un modèle de réservoir calé et cohérent vis à vis des propriétés d&#39;écoulement
CN113219553A (zh) 一种烃源岩总有机碳含量预测方法及装置
CA2711926A1 (fr) Procede, programme et systeme informatique de construction d&#39;un modele geologique 3d
CN112666610B (zh) 地层有机碳含量计算方法及系统
Zabalza-Mezghani et al. Constraining reservoir facies models to dynamic data-impact of spatial distribution uncertainty on production forecasts
Landa et al. Sensitivity Analysis of Petrophysical Properties Spatial Distributions, and Flow Performance Forecasts to Geostatistical Parameters Using Derivative Coefficients
Baville Stratigraphic correlation uncertainty: On the impact of the sediment transport direction in computer-assisted multi-well correlation (Phd, Université de Lorraine)
Lopez–Sabater et al. NEURAL–NETWORK–BASED ALGORITHMS OF HYDRAULIC ROUGHNESS FOR OVERLAND FLOW
Yadavalli et al. Use of pressure transient data to estimate permeability variograms

Legal Events

Date Code Title Description
EEER Examination request
MKLA Lapsed

Effective date: 20160316