Problèmes direct et inverse – Déformation graduelle

Calage d’historique de production

 Le problème inverse consistant à intégrer les données de production tels les débits et les pressions est appelé, en ingénierie de réservoir, calage d’historique de production. Il consiste à identifier un modèle de réservoir, qui une fois donné en entrée au système d’équations d’écoulement, donne des réponses proches des données dynamiques connues. On distingue deux types de méthodes pour résoudre ce type de problème. Le premier regroupe les méthodes de type Monte Carlo ou probabilistes (Evensen, 1994, 2009). Ces méthodes sont fondées sur la génération d’un ensemble de réalisations de la variable recherchée, ce qui donne une distribution a priori. On résout le problème direct pour chaque réalisation jusqu’au temps où on a la première mesure. Chaque variable est alors mise à jour à partir de la théorie des filtres de Kalman. Ce processus est répété jusqu’à ce que toutes les données soient assimilées. Les filtres de Kalman ont beaucoup été utilisés. Cependant, il reste un certain nombre de questions à propos du surparamétrage, de la mauvaise quantification des incertitudes ou encore de la gestion des solutions non-physiques (Gu et Oliver, 2005). Le second type de méthodes regroupe les approches variationnelles (Jacquard, 1965). Dans ce cadre, la solution recherchée minimise une fonction que l’on appelle fonction objectif. Cette fonction quantifie l’écart entre les données dynamiques et les réponses simulées pour le modèle de réservoir. Ce processus revient à implémenter un algorithme d’optimisation itératif. A chaque itération, les paramètres du réservoir sont ajustés et on effectue une simulation d’écoulement. Ces paramètres sont modifiés jusqu’à ce que les réponses simulées correspondent aux données dynamiques.  Le but est de construire un modèle de réservoir qui respecte au mieux les données disponibles, aussi bien statiques (mesures aux puits,…) que dynamiques (débits, pression,…). La boucle de calage montre comment on contraint le modèle de réservoir. La description détaillée de cette boucle de calage est disponible dans Le Ravalec et al. (2012 OGST), ici nous ne faisons qu’un résumé.  Premièrement, on construit un modèle géologique, sur une grille fine, à partir des données statiques telles que les mesures sur carotte, les diagraphies ou les modèles de sédimentation ou de déposition.  La deuxième étape consiste à amener ce modèle vers une échelle grossière. On parle ici d’ « upscaling ». L’intérêt est de réduire le temps de calcul lié à la simulation d’écoulement. On travaille donc sur une grille de moindre résolution où les propriétés pétrophysiques sont moyennées (à l’aide d’une moyenne arithmétique, géométrique,…).  Dans un troisième temps, on donne ce modèle grossier en entrée au simulateur d’écoulement qui fournit alors des réponses en production (débits au cours de la production, pressions, « water-cut »,…). A partir du modèle géologique et des résultats en pression et saturation issus du simulateur d’écoulement, on construit un modèle petroélastique, qui nous permet de simuler les attributs sismiques (données déduites de la sismique).  La dernière étape de la boucle est le calcul de la fonction objectif, qui compare les données de production et les données sismiques simulées avec les données réelles.  Pour réduire cette fonction objectif on perturbe le modèle géologique construit à l’étape 1 à l’aide de techniques telles que la méthode de déformation graduelle. On répète cette chaine de simulation jusqu’à ce que la fonction objectif passe sous un seuil préalablement défini. La fonction objectif est définie par : Problème inverse – Calage d’historique de production .Elle mesure l’écart entre les réponses simulées g(v) et les données dynamiques dobs. g est l’opérateur qui permet de passer de l’espace des paramètres v aux réponses dynamiques : il s’agit ici du simulateur d’écoulement. CD est la matrice de covariance qui quantifie les incertitudes expérimentales et théoriques. Généralement, un modèle de réservoir est composé de millions de mailles. Donc, il y a au moins autant d’inconnues que de mailles : chaque maille doit contenir une valeur de porosité, de perméabilité, de saturation, etc. D’un autre côté, le nombre de données dynamiques est nettement moins grand, ce qui rend le problème mal posé. Pour que le problème soit « mieux » posé, on utilise différentes techniques de régularisation. On peut par exemple ajouter de l’information dans la fonction objectif . Le terme ajouté sur la droite évalue l’écart entre le modèle courant v et le modèle a priori vo. La matrice de covariance CV caractérise les incertitudes sur vo. Une autre technique de régularisation consiste à limiter l’espace de recherche. On peut alors utiliser des méthodes de paramétrisation géostatistiques telle que la méthode de déformation graduelle (Hu, 2000). Cette méthode, intégrée dans un processus d’optimisation, permet d’ajuster le paramètre de déformation de manière à minimiser la fonction objectif (cette méthode est détaillée section 2.3). Comme la méthode de déformation graduelle tient compte de l’information a priori, la fonction objectif est réduite au terme de différence entre les données. La méthode de déformation graduelle est rappelée plus loin. De manière plus générale le problème d’optimisation s’écrit : min 𝑣∈Ω 𝐽(𝑣) Equation 7 On rappelle que J est la fonction objectif et v les paramètres inconnus du modèles tels que la porosité, la perméabilité, etc. J est généralement fortement non linéaire dans le cadre du calage d’historique. Différentes méthodes d’optimisation sont utilisées pour résoudre ce problème : méthodes stochastiques (Gao et al., 2007) ou déterministes (de Montleau et al., 2006), méthodes globales (Bouzarkouna, 2012) ou locales (Bissell et al., 1994), etc. Les méthodes globales permettent un meilleur calage mais sont généralement trop couteuse en temps de calcul. Elles demandent un grand nombre d’évaluations de la fonction objectif et donc de simulations d’écoulement. Les méthodes locales sont souvent préférées. Nous détaillerons ici les 3 catégories en les expliquant rapidement Les méthodes de descente. Ces méthodes s’appuient sur le calcul l’estimation du gradient ou de la Hessienne de J. En effet, la première étape de ces algorithmes est de partir d’un point initial et de le mettre à jour en se déplaçant dans la direction de descente. C’est cette direction de descente qui est calculée à partir des dérivées au premier et second ordre. Cependant dans les cas pétroliers ces dérivées n’ont jamais d’expression analytique, il faut donc calculer des approximations. On peut estimer le gradient par état adjoint. Cela consiste à créer un problème adjoint à la simulation d’écoulement et à résoudre les deux problèmes. Cependant peu de simulateur Problèmes direct et inverse – Déformation graduelle 25 d’écoulement proposent un code adjoint. On peut estimer le gradient par différences finies, en perturbant chaque paramètres. Mais cela implique de faire une évaluation de la fonction objectif pour chaque perturbation. Enfin citons la méthode bien connue de Quasi-Newton, qui s’appuie sur l’estimation de la Hessienne pour la direction de descente. De même la Hessienne est approximée à partir des calculs du gradient à chaque itérations.  Les méthodes de recherche par région de confiance Ces méthodes sont basées sur la construction de modèles approchés de J dans une région de confiance, recalculée à chaque itération. Ces modèles, souvent quadratiques, sont très facile à optimiser. Ils peuvent être calculés à partir de l’estimation du gradient et de la Hessienne comme dans SQP (Sequential Quadratic Programing) ou par interpolation (Bissell, et al., 1994).  Les méthodes de recherche directe Ces méthodes ne nécessitent pas la connaissance des dérivées et sont très simples à implémentées. La première étape consiste à échantillonner la fonction objectif sur un certain nombre de points, puis de décider des prochaines évaluations à faire en se basant sur l’échantillonnage précédent. On retrouve deux types de méthodes : les méthodes par recherche directionnelle (Conn et al., 2009) et les méthodes par simplexe (Nelder et Mead, 1965). Une des principales préoccupations lors d’un calage d’historique de production est le temps de calcul, qui augmente avec le nombre d’itérations. Beaucoup de techniques ont été étudiées pour accélérer les techniques d’optimisation. L’approche développée ici se concentre sur la réduction des dimensions de l’espace de recherche, c’est-à-dire du nombre de paramètres inconnus. Comme expliqué précédemment, elle fait intervenir une paramétrisation multiéchelles des réalisations représentant les variations des propriétés pétrophysiques dans le réservoir. L’idée principale est de se donner la possibilité de changer les propriétés pétrophysiques à différentes échelles tout en préservant les liens entre les échelles. Pour résumer, on veut pouvoir changer les tendances d’une propriété à l’échelle grossière ou les variations plus fines de cette même propriété à une échelle plus fine. L’intérêt de cette approche est de gérer un nombre de paramètres qui dépend de l’échelle considérée. Ainsi, à l’échelle grossière, le nombre de paramètres inconnus est plus petit qu’à l’échelle fine. On cherche donc à déterminer d’abord les paramètres à l’échelle grossière. Puis, dans un deuxième temps, on envisage de raffiner le modèle de la propriété recherchée en variant les paramètres à l’échelle fine. Chaque fois qu’un modèle est raffiné, le nombre de paramètres augmente.

Déformation graduelle 

La méthode de déformation graduelle est une technique de paramétrisation géostatistique utilisée pour perturber une réalisation d’une variable aléatoire à l’aide d’un nombre réduit de paramètres tout en préservant la variabilité spatiale. Son schéma de base implique la combinaison de deux bruits blancs Gaussiens indépendants : () cos sin 1 2 z  z  z Equation 8 z1 et z2 sont deux bruits blancs Gaussiens indépendants. Quelle que soit la valeur du paramètre de déformation θ, on peut montrer que z est aussi un bruit blanc Gaussien (Hu, 2000). Comme la règle de déformation est périodique, θ est compris entre -1 et 1. Lorsque θ = 0, z est égal à z1 ; lorsque θ = ½, z est égal à z2. Le bruit blanc Gaussien correspond aux nombres aléatoires attribués aux mailles où on veut simuler des valeurs de la variable. Faire varier le paramètre de déformation permet de varier le bruit blanc Gaussien utilisé pour générer la réalisation. Ce processus permet de déformer la réalisation déduite du bruit blanc Gaussien. Cette technique de paramétrisation, une fois intégrée au processus d’optimisation ou au calage d’historique, est un outil qui permet de parcourir des chaînes de réalisations dans l’espace de recherche. Parmi toutes celles-ci, on peut identifier une réalisation qui permet la décroissance de la fonction objectif. La probabilité d’avoir un bon calage des données dynamiques est très faible si on se contente de parcourir une seule chaine. Pour être efficace, la recherche par déformation graduelle doit se faire sur Problèmes direct et inverse – Déformation graduelle 27 plusieurs chaines. Pour plus de détails sur la méthode de déformation graduelle, un ebook a été publié par Le Ravalec et al. (2014). Figure 8 Schéma d’une chaine d’optimisation à un paramètre de déformation graduelle θ Pour chaque chaine explorée, on a un processus d’optimisation. Au final, l’optimisation par déformation graduelle correspond à une séquence d’optimisations comme schématisée Figure 8. Une première chaine est construite à partir de deux bruits blancs Gaussiens (z1 et z2). On lance un premier processus d’optimisation qui donne la réalisation réduisant le plus la fonction objectif. On obtient donc un θoptim qui donne le bruit blanc Gaussien qui réduit le plus la fonction objectif dans cette première chaine. Comme une chaine seule ne représente qu’une toute petite partie de l’espace de recherche, il faut réitérer le processus. Le bruit blanc Gaussien correspondant à la réalisation optimale du processus d’optimisation précédent est utilisé pour mettre à jour z1 que l’on appelle ici dans la Figure 8, z3(θoptim), et un nouveau bruit blanc Gaussien est tiré au hasard pour z4. La combinaison de ces deux nouveaux bruits blancs Gaussiens donne une nouvelle chaine de réalisations qui une fois parcourue amène à une réalisation optimale qui réduit la fonction objectif autant que possible, et ainsi de suite. Dans ce contexte, on espère que le schéma multi-échelles apporte plus de flexibilité. Une réalisation dépend de deux bruits blancs Gaussiens : un à l’échelle fine et un à l’échelle grossière. Le bruit blanc Gaussien à l’échelle grossière permet de générer la tendance sur la grille grossière. Le bruit blanc Gaussien à l’échelle fine génère la réalisation sur la grille fine conditionnellement à la réalisation générée précédemment sur la grille grossière. En conséquence, la méthode de déformation graduelle peut être appliquée sur la grille grossière (Figure 9), sur la grille fine (Figure 10) ou sur les deux en même temps. La Figure 9 montre l’impact de la déformation graduelle sur la réalisation simulée sur la grille grossière et les modifications qui en résultent pour la réalisation sur la grille fine. La Figure 10, quant à elle, montre l’impact de la déformation graduelle sur la réalisation à l’échelle fine pour une réalisation grossière fixée. Pour ce dernier exemple, on remarque que la tendance reste la même à l’échelle fine, on joue seulement sur les variations autour de la tendance. L’avantage de cette méthode de simulation imbriquée est qu’il est possible de changer une réalisation à une échelle donnée, avec un nombre d’inconnues dépendant de cette échelle. Plus l’échelle est grossière, moins il y a de mailles dans la grille associée et donc, moins il y a d’inconnues. Le calage devrait s’en trouver facilité : les données peuvent être intégrées à la bonne échelle avec moins d’inconnues. 

Simulation séquentielle Gaussienne 

Dans ce chapitre nous développons un modèle de réservoir issu de la géostatistique à deux points. Nous utilisons donc un variogramme pour modéliser les propriétés pétrophysiques du réservoir. Après une brève bibliographie sur les modèles couramment utilisés, nous présentons la simulation séquentielle Gaussienne, puis nous détaillons l’approche que nous avons proposée pour l’étendre à deux échelles. Nous l’illustrons dans un premier pour des propriétés continues, puis nous expliquons comment procéder pour des variables discrètes. Enfin nous mettons en place un cas synthétique pour mettre en évidence les points fort de la méthode. 

Etat de l’art 

Les méthodes géostatistiques sont un ensemble de méthodes statistiques, dédiées à l’analyse de données définies par des coordonnées spatiales et temporelles. Jusqu’à la fin des années 1980, elles étaient essentiellement vues comme un moyen de décrire les variations spatiales à partir d’un variogramme et de prédire les valeurs des propriétés aux points non échantillonnés par krigeage (Vieira et al., 1983 ; Trangmar et al., 1985). Le krigeage, ou estimation, est un nom générique adopté par la communauté des géostatisticiens pour décrire une famille d’algorithmes de régression par moindre carré. Il existe de nombreuses méthodes de krigeage : simple, ordinaire, universel, ou avec tendance, cokrigeage, krigeage avec dérive externe, etc. Une présentation détaillée de ces méthodes mathématiques a été faite par Goovaerts (1997a, pp 125- 258). Le krigeage a tendance à lisser la variabilité spatiale de la propriété modélisée (Goovaerts, 1999). Au contraire, les méthodes de simulation stochastiques ne minimisent pas la variance de l’erreur, mais se concentrent sur la reproduction des propriétés statistiques, notamment la moyenne, la variance et la covariance (ou variogramme). La simulation stochastique est préférée au krigeage pour les applications où la variabilité spatiale de la propriété doit être préservée. La simulation peut être faite à partir d’une large palette de techniques qui dépendent du modèle de la fonction aléatoire, du nombre et du type d’informations disponibles, et des spécifications sur le temps de calcul. Une revue détaillée est disponible dans le livre de Goovaerts (1997a, pp 376-424), notamment sur les techniques de recuit simulé (Deutsch et Lewis, 1992), d’approche par champ de probabilité (Froidevaux, 1992) et des algorithmes booléens, que nous ne détaillerons pas ici. Nous allons plutôt détailler la simulation séquentielle, méthode dans laquelle s’inscrit la simulation séquentielle Gaussienne. La simulation séquentielle : au lieu de modéliser la fonction de répartition conditionnelle à N-points, on calcule une fonction de répartition conditionnelle à un point sur chaque nœud visité selon une séquence aléatoire. Pour s’assurer de reproduire le variogramme, chaque fonction de répartition conditionnelle est conditionnée non seulement aux données disponibles, mais aussi Algorithme à une échelle 30 aux points déjà simulés. Deux classes se distinguent : les techniques pour la simulation de fonctions aléatoires multi-Gaussiennes ou d’indicatrices. o La simulation séquentielle Gaussienne (Desbarats, 1996 ; Deutsch et Journel, 1998) : la grille est parcourue séquentiellement. Pour une maille donnée, la moyenne et la variance de la loi conditionnelle (supposée Gaussienne) sont obtenues par krigeage. Cette méthode est détaillée dans la section 3.2 o La simulation séquentielle par indicatrices (Goovaerts, 1997b ; Garcia et Froidevaux, 1997) : la fonction de distribution se détermine à partir du calcul de quelques points de la courbe de distribution (calcul discret de la probabilité pour un certain nombre de seuils convenablement choisis). La simulation non paramétrique utilise des indicatrices.

 Algorithme à une échelle

Nous cherchons à présent à simuler une réalisation d’une variable aléatoire stationnaire V(u), connue en n points ui, i ∈ à [1,n], de moyenne m, de variance σ 2 et de covariance C V (h). Les n valeurs connues de V sont aussi appelées données dures. Pour ce faire, nous appliquons la méthode de simulation séquentielle Gaussienne, identifiée par l’acronyme SGSim. Les principales étapes de cet algorithme sont présentées ci-dessous. i) Remplir la grille de départ avec un bruit blanc Gaussien. ii) Déterminer un chemin aléatoire pour parcourir toute la grille. iii) Pour chaque maille i visitée, vérifier qu’il n’y a ni donnée dure, ni valeur déjà simulée. Si c’est le cas :  chercher dans un voisinage proche toutes les données disponibles (les données dures et les valeurs déjà simulées)  calculer les moyenne m SK et variance σ 2 SK de la fonction de densité de probabilité (pdf) en utilisant un krigeage simple.Les coefficients λ sont les poids du krigeage. Ils sont obtenus par la résolution du système C V λ =B V . C V est la matrice de covariance calculée pour les mailles avec des données ou des valeurs déjà simulées. B V est le vecteur des covariances pour les distances entre la maille i visitée et les mailles voisines contenant des données ou des valeurs déjà simulées.  à partir de cette pdf et du nombre aléatoire assigné à la maille courante, on obtient une valeur que l’on attribue à la maille.