FR 2776393 A1 19990924 FR 9804084 A 19980320 La présente invention concerne une méthode optimrniséec de modélisation de la structure stratigraphique d'une zone souterraine permettant l'ajustement rapide de paramètres du modèle avec des données observées ou mesurées. La méthode selon l'invention trouve des applications dans la modélisation de bassins sédimentaires en vue notamment d'aider les ingénieurs de réservoir à mieux localiser les gisements souterrains d'hydrocarbures. Les bassins sédimentaires sont évolutifs au cours du temps, ils sont de géométrie variable sous l'effet de la compaction (le volume poreux diminue), de la subsidence (déformation du fond du bassin) et de l'érosion (enlèvement d'une partie des sédiments supérieurs du bassin). Dans un milieu marin, des sédiments viennent se poser sur le socle pour former un bassin sédimentaire. On cherche à estimer la hauteur du sédiment et les teneurs en lithologies sous l'influence de la compaction, de la subsidence et de l'érosion. Pour définir tous les paramètres intervenant dans la définition et la formation d'un tel bassin sédimentaire, on va utiliser les notation suivantes: Q: Domaine d'étude qui représente le bassin sédimentaire, Q c IR2. N: Nombre de lithologies, N > I. E',:Efficacité de l'eau, (unité [0]: dépendante de la lithologie 1. qt,l= 1...N: Flux sédimentaire de la lithologie I sur la frontière du domaine, il dépend du temps, (unité [L2T-I], Lest la longueur et Test le temps). v,,l = 1,.... N: Teneur en lithologie 1, elle dépend du temps, (unité: []). QL, = 1... N: Flux sédimentaire de la lithologie 1, (unité: [ L-2T-]). H(x,v,t): Hauteur du sédiment, (unité: [L]. Y D t: Porosité, (unité: [ó]). K,(h): Coefficient de diffusion, appelé aussi diffusivité, il dépend des variables d'espace et de temps, (unité: [ L2T-']). L'ensemble des équations décrivant le processus de sédimentation est le suivant: - Equation en h: donnée par une loi de conservation de la massc. - Equation en E,,,.; elle correspond à l'efficacité de l'eau (transport des sédiments par l'eau). - Equation liée au coefficient de diffusion K, (h), elle fait intervenir la bathymétrie et permet de définir le coefficient de diffusion relativement à trois zones différentes: continentale, marnage (zonc de marée) et marine. La valeur du coefficient de diffusion sera donc définie à partir des trois valeurs KI ( K,,,,.,) K2 ( K,,,,,.) et K3 ( Kn,er) (Fig. 1) dont le raccord entre la zone continentale ct marine est linéaire. - Equation liant le flux Q; à la hauteur h de la lithologie i: Qi =-K, vE,,,, Vh, elle est donnée par une loi d'écoulement similaire à la loi Darcy. Pour définir le modèle direct d'un bassin, il est connu de coupler, pour chaque lithologie, des équations de transport avec celles de conservation de la masse dans le but de définir en chaque point du bassin le taux d'érosion ou de sédimentation et la teneur en chaque lithologie. On représente d'abord le bassin sédimentaire par -Q domaine ouvert de IR2 de frontière F = rF ur2 et l'on doit trouver h,Qi et v;,i = 1...N tels que: -(v,'(Dh) + divQ, = gi sur 2x[O, T[ Q; =-K, v E 0 tel que K > oai, h sera solution de l'équation parabolique suivante LviVi2h+div Yfi(h)WhïJ= Y i (l) o (fi (h) -K;v, Eo,,) On connait dans le domaine des géosciences, des méthodes d'inversion permettant de contraindre un modèle de simulation initial d'une zone souterraine telle qu'un bassin sédimentaire, résultant par exemple d'une simulation géostatistique, par des données observées ou mesurées dans la zone, suivant un processus itératif automatique. Plusieurs techniques existent pour résoudre les problèmes d'inversion. Une technique connue de type essai-erreur consiste à appliquer un algorithme d'optimisation de type génétique calqué sur le processus de l'évolution par la sélection naturelle. Un exemple de cette technique appliquée à la modélisation stratigraphique de bassins sédimentaires, est décrite par exemple par: - Bornhotdt S., Optimization using Genetic Algorithms; Proceedings Numerical Experiments in Stratigraphy ((NES) 1996. La mise en oeuvre de cette technique est facile et ne nécessite pas de calcul de gradients. Cependant, dans certains tests numériques, du fait du choix aléatoire de la population initiale, on a besoin de plusieurs générations et d'une grande population pour diminuer la fonction critère, ce terme désignant l'erreur des moindres carrés entre les observations (valeurs de référence) et les prédictions (valeurs calculées). Une autre méthode connue consiste à utiliser un algorithme itératif de descente pour l'ajustement des paramètres, en supposant que la fonctioncritère respecte certaines formes. Une technique de ce type, appliquée à la modélisation stratigraphique de bassins sédimentaires, est décrite par exemple par - Lessenger M. et aI.: Forward and Inverse Simulation Models of Stratal Architecture and Facies Distributions in Marine Shelf to Coastal Plain Environments. Thesis, Colorado School, November 1993; ou par - Lessenger M. et ai.: Estimating Accuracy and Uncertaintv of Stratigraphic Predictions from Inverse Numerical Models; in Proceedings Numerical Experiments in Stratigraphy (NES). 1996; Par le brevet FR-A-2 744 224 du demandeur. on connaît également une méthode par essais et erreur de type manuelle qui peut être utilisée pour répondre à un problème d'ajustement des paramètres géologiques d'un processus stratigraphique. lorsqu'on a des connaissances géologiques fortes dans l'environnement des dépôts sédimentaires. La méthode selon l'invention permet une modélisation automatique de la structure stratigraphique d'une zone souterraine avec un ajustement rapide de paramètres du modèle avec des données observées ou mesurées. Elle comporte - la détermination par mesure ou observation d'une série de grandeurs représentatives de la structure stratigraphique de la zone en différents points (telles que des hauteurs de sédiments ou des données provenant d'interprétations d'essais de puits ou des mesures sismiques notamment) - l'implémentation d'un modèle initial direct dépendant d'un ensemble de paramètres représentatifs du processus de formation géologique de la zone - la détermination au moyen du modèle initial, des valeurs prises par les grandeurs aux différents points de mesure ou d'observation - le calcul au moyen du modèle direct. des gradients de chaquc grandeur par rapport aux paramètres à ajuster; et l'optimisation du modèle eni minimisant par itérations successives la fonction critère, jusqu'à réduction suffisante des écarts, par le calcul de la sensibilité du modèle direct aux paramètres à ajuster, de manière à produire une représentation stratigraphique de la zone souterraine. La méthode comporte par exemple l'utilisation d'un modèle de type diffusif. Dans le cas o le modèle direct sélectionné est dérivable par rapport aux paramètres, la méthode comporte par exemple l'utilisation d'un gradient de type analytique. La méthode peut aussi comporter l'utilisation de gradients par différences finies. La méthode peut comporter également l'utilisation par exemple d'un algorithme quasi-newtonien pour minimiser la fonction-critère. La méthode selon l'invention permet d'obtenir une convergence rapide du modèle stratigraphique choisi initialement avec les mesures et observations des hauteurs sédimentaires mesurées, ce qui conduit à minimiser les coûts des calculs d'ajustement du modèle. On utilise un dispositif permettant l'implémentation de la méthode qui comporte: - des moyens pour déterminer par mesure ou observation d'une séric de grandeurs représentatives de la structure stratigraphique de la zone en différents points; et - un ensemble de traitement informatique programmé pour réaliser une implémentation d'un modèle initial direct dépendant d'un ensemble de paramètres représentatifs d'un processus de formation géologique de la zone. cet ensemble de traitement incluant: des moyens pour déterminer au moyen du modèle initial, des valeurs prises par les grandeurs aux différents points de mesure ou d'observation - des moyens pour calculer à partir du modèle direct. des gradients de chaque grandeur par rapport aux paramètres à ajuster; - des moyens d'optimisation du modèle en minimisant par itérations successives une fonction critère dépendant des dits paramètres. jusqu'à réduction suffisante des écarts, par le calcul de la sensibilité du modèle direct aux paramètres à ajuster; et - des moyens pour produire une représentation stratigraphique de la zone souterraine. D'autres caractéristiques et avantages de la méthode selon l'invention. apparaîtront à la lecture de la description ci-après d'un exemple non limitatif de réalisation, en se référant aux dessins annexés o: - la Fig.l montre un exemple de variation du coefficient de diffusion dans une bordure littorale; - la Fig.2 montre un organigramme représentatif du processus d'optimisation d'un modèle stratigraphique - la Fig.3 montre, en fonction de h, la variation d'un terme source représentant par exemple les apports sédimentaires; - la Fig.4 représente un algorithme de calcul du gradient' - la Fig.5 représente un algorithme de calcul de paramètres optimaux; - la Fig.6 représente un algorithmc préféré utilisé pour minimiser unc ftonction- critère - la Fig.7 représente un organigramme du processus d'inversion utilisé pour optimiser le modèle stratigraphique; - les Fig.8A, 8B, représentent deux exemples de variation dans l'espace de la hauteur de sédiments respectivement à deux instants successifs - la Fig.8C représente la variation de la fonction-critère J par rapport au coefficient de diffusion K; - les Fig.9A, 9B représentent des exemples de variations dans l'espace de la hauteur de sédiments h, la première à un instant initial, la deuxième à trois instants successifs t, , t. et t,; - la Fig.9C représente un exemple de la variation en fonction des coefficients de diffusivité KI (K,,t.) et K2 (K,nr) de la fonction critère J - la Fig. 10 est un tableau comparatif des valeurs optimales atteintes par les paramètres qt I, Kfi1t,K--, par référence à un premier jeu de valeurs des paramètres de départ, avec K, =.1; - la Fig. I 1 est un tableau comparatif des valeurs optimales atteintes par les paramètres q -""a, KIt"'', K2/""t par référence à un deuxième jeu de valeurs des paramètres de départ, avec K3 ( Kh,,h,) =1; et - la Fig.12 est un tableau comparatif des valeurs optimales atteintes par les paramètres q tiK,, K/,,,, K",, par référence à un troisième jeu de valeurs des paramètres de départ, avec K3 ( Kl,,,,,) = 1. Description de la méthode Dans ses grandes lignes qui sont schématisées à la Fig.2, ct développées dans la suite de la description, la méthode comporte la succession d'étapes suivantes - Choisir un modèle direct (de type diffusif par exemple) pour simuler l'objet géologique réel - Choisir un jeu initial de paramètres qui vont permettre avec le modèle direct, d'aboutir aux grandeurs (telles que les hauteurs sédimentaires par exemple) que l'on cherche à faire coïncider avec les grandeurs réelles - Calculer en utilisant le modèle direct, les grandeurs (hauteurs sédimentaires par exemple) qui en résultent aux points de mesure et d'observation (tels que les puits forés au travers du milieu) et les comparer aux observations réelles aux mêmes points; - Calculer au moyen du modèle direct, des gradients de chaque grandeur par rapport aux paramètres à ajuster, que ce soit des gradients de type analytique si le modèle direct est dérivable, ou des gradients de type par différences finies en faisant varier les paramètres. Si les écarts entre le modèle et les observations sont grands, optimiser le modèle ainsi créé en minimisant de manière itérative la tfonctioncritère, de préférence par un algorithme quasi-newtonien, de préférence via une méthode connue d'optimisation dite BFGS. On va définir ci-après les termes dans lesquels se posent le problème d'optimisation en utilisant les notations suivantes h"h, pour désigner la hauteur observée. P = {p = (Pi)>o0 E IR"}, pour désigner l'ensemble des paramètres P à identifier; et m, pour désigner le nombre de paramètres à ajuster. L'ensemble Pd, des paramètres admissibles (P, c P) est défini par l'ensemble des paramètres qui sont sous contraintes. Si p est le paramètre L identifier, on note h(.:p) la solution du problème direct. J(p) qui représente l'écart entre les quantités calculées h(.;p) et les quantités observées, h"";- est appelée fonction objectif ou critère. Il s'agit dans ce contexte de trouver e P,/ tel que J(ô) L2(çi), 'fI2 =(fn f|)',2). Le paramètre p doit vérifier VJ(') = O. Le problème inverse est un problème de minimisation par rapport aux paramètres. Suivant la définition de l'ensemble des paramètres admissibles (égal ou inclus dans l'espace des paramètres inconnus), on a à traiter un problème de minimisation sans contraintes ou avec contraintes. Dans le cas o K(h) = cte, Vh et o K est le seul paramètre qu'on cherche à identifier, il faut, pour estimer l'optimum K, calculer 7h(; K), gradient de h par rapport à K, du fait que VJ = 2Vh(.;K).(h(.; K) - h""), et à trouver à l'aide des techniques d'optimisation, le paramètre K qui vérifie VJ(K) = 0. Dans le cas o la fonction K est donnée en fonction de K,,,.,. K, et K,,, ,,, de la Fig. 1, il faut calculer les gradients de h par rapport à K,,,, ,,, K,. et K,,,, et définir les valeurs optimales Ke,,.,, K,,, et Kh,,,h, par un procédé d'optimisation. Avec la méthode selon l'invention, on actualise les paramètres par le calcul du gradient. et pour cela. on peut utiliser plusieurs techniques de calcul différentes. Méthode de gradient pour le simulateur stratigraphique Cas d'une couche homogène Le paramètre K(h) est ici indépendant de la hauteur du bassin (on est soit en dessous de l'eau, soit au-dessus) et l'on considère que la porosité (D0 est gl, KI indépendante du temps. Si l'on pose: q = -, K= K' I'équation à résoudre est du type: dh i div(KEe,,,,,Vh) = q(h) (5) qui se réduit à d - d K(h) =q (6) dans le cas d'une couche monodimensionnelle et l'on utilise pour sa résolution un code de résolution numérique d'un type connu, que l'on désigne par 'codediffusion". Le terme source q peut correspondre à la production des carbonates ou des sédiments siliciclastiques; une représentation graphique possible est dans la Fig. 3. Méthode de dérivation du modèle direct D'après ce qui précède, pour résoudre le problème d'inversion lié à cette équation, il faut d'abord estimer Vh, qui représente le gradient de hI par rapport à p o p est le paramètre à identifier. ll Restriction au cas unidimensionnel L'équation satisfaite par le gradient Vh, est, on le montre, duI même type que celle du modèle direct. Dans le cas o q(Ih) = =,nction (h, s, l,),s1.s, E IR (cf. Fig.3), il faut estimer ah ah s - et - et chercher ensuite le minimum de J(K,sI s,). L'algorithme résumant ce qui précède est schématisé suLIr la Fig. 4: Formulation générale de la méthode de dérivation pour le calcul du gradient Soit U l'ensemble des inconnues et d l'ensemble des paramètres d'un modèle direct. On note U"' la solution approchée au temps t"'. Le problème (P) peut s'exprimer, après discrétisation. sous la forme suivante: f (U, U ", d) = ((8) que l'on résoud par une méthode classique de type Newton. La méthode du gradient pour l'inversion des paramètres consiste à dériver l'équation (8) par rapport au paramètre d; dans ce cas Y'+l = dz est solution du svstème linéaire suivant Yil yn+l=_( f_ Qf V(9) dUl+I & dU" ', Une application de cette méthode est décrite par exemple par Rahon, D. et al: "Gradient Method constrained Geological Bodies for History Matching. " SPE 36568: Society of Petroleum Engineers. 1996. Méthodes numériques de calcul du gradient Une première méthode consiste à dériver les équations discrétisées du modèle direct, et par exemple l'équation (5) dans le cas d'une lithologie unique. Une deuxième possibilité consiste à utiliser une méthode cide différences finies Fl -h pour approcher le gradient V,i(p) par ôi ou h = h(. p+A,). La perturbation Ap est obtenue en perturbant l'élément i'"" du vecteur p Ait, = (Po, P,....Pi +)... On Fait une première simulation avec le paramètre p pour obtenir h. puis on fait une nouvelle simulation avec le paramètre (p, p,p..., p, +...) pour évaluer h(.; p+ Ai,). Méthodes d'optimisation On connaît des méthodes pour l'optimisation d'une tfonction itérative qui se basent sur la recherche d'un point stationnaire c'est-à-dire le point qui annule le gradient de la fonction et parmi elle particulièrement la méthode dite du gradient conjugué pour les fonctions quadratiques et la méthode BFGS pour les fonctions quelconques. On choisit de préférence, la méthode dite BFGS qui est très faiblement sensible aux imprécisions dans la procédure de recherche unidimensionnelle. Cela permet d'utiliser des méthodes d'optimisation unidimensionnelle économiques qui ne nécessitent qu'un très petit nombre d'évaluations de la fonction à chaque itération. Pour la méthode BFGS. les approximations HA du hessien sont définies par Ho = I (l est la matrice identité) et par la suite récurrente suivante: H, =H +l [1,]e IJ 5C H, + H4& Hk+= Hk +] - (20) o 5. = p - p et y7 =VJ(pk])-VJ(pL). x' désigne le transposé du vecteur x (matrice-ligne) et x. v désigne le produit scalaire des vecteurs x et v. L'algorithme BFGS de minimisation du critère est schématisé à la Fig.6. Application de l'algorithme d'inversion On va utiliser l'algorithme BFGS décrit ci-dessus pour définir un modèle inverse pour l'identification des paramètres intervenant le modèle direct suivant dont la formulation mathématique. s'écrit, pour i = I....N, sous la forme suivante t V (vOih) +div(-K, v; E".Vh) = gi sur Q x [0, T[ ZIV i=l Vi 1 h=ho sur F x[0,T[ vi E,,, KiVh.i = fi sur r2, x[0, T[ h(x, y,0) = g(x, v) sur Q o Q représente le bassin sédimentaire, N le nombre de lithologies et Ki les diffusivités qui peuvent être constantes ou dépendantes de h comme dans la figure I de la section 1. Dans ce cas, l'ensemble P des paramètres qu'on voudrait ajuster contient les coefficients de diffusion Ki K.;i et K:,i (respectivement à terre, en bathymétrie et en mer) pour chaque lithologic i. les termes sources g; et les flux.fi au bord F,. Nous envisageons deux cas: I) Si Vi, gi est constantc indépendamment de h alors P = I Kl.,, K.,,,, Ka./:...:Kiv. K.; N. K2.N; gl.....g,v; 1l...ffN r' h. 2) Si les fonctions gi = fJction(h, s.i,s,,.i,) sont données par la loi de la Fig. 3, c'est-à- dire si les gi sont définies via les coefficients xj,, s2,. alors P= { Ki./. K,.. K2,:...:KtN. K.,N. K.,N: SI N I..... 2),N f.,,fN }. Organigramme du modèle inverse Ayant fixé les paramètres à identifier, la procédure à suivre est résumée dans l'organigramme de la Fig.7. Tests numériques On traite dans ce qui suit de l'inversion des paramètres du modèle diffusif dans le cas monolithologique (N = 1) sous des conditions aux limites de type Neumann et Dirichiet au moyen de l'équation aux dérivées partielles: dgh (E) --div(KE,,,,Vh) =q Inversion dans le cas linéaire On a réalisé à titre de test, I'inversion du coefficient de diffusion K(h)=cte, h en cherchant à caler le modèle numérique sur un autre modèle numérique (jouant le rôle du modèle réel) tel que h"h''(.,.) = h(.,.; K,),.......), K",,t...f.. = 12 A0 o h(... Kf,,pi,,) est la solution dc I'équation (E). On minimise enfin le critère J(K). on choisit K'=4 conmme point de départ de la méthode BFGS. On vérifie que le coefficient optimal trouvé coïncide exactement avec K,,,,i,,,,,,. Ceci est illustré sur les Fig.9A à 9C o la troisième représentation correspond au critère J en fonction de K. La convergence vers le minimum est ici obtenue au bout de trois itérations seulement. On observe que l'on obtient le nmême résultat dans le cas o le critère est évalué sur tous les points du domaine ou sur quelques uns seulement (c'est le cas si les observations l""' sont connues en quelques puits seulement). Inversion dans le cas non linéaire Dans cette partie, on s'intéresse à l'étude numérique de l'équation (E) dans le cas o K(h) est défini en fonction des trois paramètres KI. K3 et K, (cf. Fig.1).On fixe K3 et l'on inverse les paramètres q. Ki et K2.. Les Fig. 10B, I B correspondent à la représentation de quelques simulations directes et du critère J(K, K2) o K =. 1 et qc = 5. Le coefficient de diffusion initial K(h) et la condition initiale sont représentés respectivement sur les Fig. I 1 A. 1 OA. Les tables des Fig.l I et 12 correspondent à l'inversion du paramètre p = (q. K., K2) avec des observations calculées à l'aide du paramètre P",rine,,, = ( 10, 5.2.5). La table de la fig. 13 correspond à l'inversion du paramètre p = (q. KI, K2) avec des observations calculées à l'aide du paramètre Ppti,,, = (10. 5, 3). Le critère J(K,. K_,) qui correspond à ce cas, ainsi que la condition initiale h(.,.,O) sont représentés dans les Fig.14A-14C. On constate d'abord que le calcul du minimum de J dépend des valeurs initiales qo, Ky et K". Comme on véritfic sur Ics tables des Fig. I I- 13. la convergence vers l'optimum est assurée et ceci bien que les différents points de départ sont choisis très loin de l'optimum P On a décrit des exemples de mise en oeuvre dc la méthode o la grandeur considérée était une hauteur (épaisseur) de dépôts sédimentaires. On ne sortirait pas du cadre de l'invention toutefois en considérant d'autres grandeurs telles que des données provenant d'interprétations des essais de puits (valeurs de pression) ou des données géologiques (position et nature de lithologies) obtenues par interprétation de mesures sismiques par exemple. REVENDICATIONS I) Méthode pour modéliser automatiquement la structure stratigraphique d'une zone souterraine avec un ajustement rapide automatique de paramètres du modèle avec des données lithologiques observées ou mesurées, caractérisée en ce qu'elle comporte: - la détermination par mesure ou observation cd'une série de grandeurs représentatives de la structure stratigraphique de la zone en différents points; - I'implémentation d'un modèle initial direct dépendant d'un ensemble de paramètres représentatifs du processus de formation géologique de la zone, - la détermination au moyen du modèle initial. des valeurs prises par les grandeurs aux différents points de mesure ou d'observation - le calcul au moyen du modèle direct, des gradients de chaque grandeur par rapport aux paramètres à ajuster; et - lI'optimisation du modèle en minimisant par itérations successives une fonction critère, jusqu'à réduction suffisante des écarts, par le calcul de la sensibilité du modèle direct aux paramètres à ajuster. de manière à produire une représentation stratigraphique de la zone souterraine. 2) Méthode selon la revendication 1, caractérisée en ce qu'elle comporte l'utilisation d'un modèle de type diffusif. 3) Méthode selon la revendication I ou 2, caractérisée en ce qu'elle comporte l'utilisation d'un modèle direct sélectionné dérivable par rapport aux paramètres. et l'utilisation d'un gradient de type analytique. A0 4) Méthode selon l'une des revendications précdecntes, caractérisée en ce qu'elle comporte l'utilisation de gradients par différences finies. ) Méthode selon l'une des revendications précédentes. caractérisée en ce qu'elle comporte l'utilisation d'un algorithme quasi-newtonien pour minimiser la fonction-critère. 6) Dispositif pour modéliser automatiquement la structure stratigraphique d'une zone souterraine avec un ajustement rapide automatique de paramètres du modèle avec des données lithologiques observées ou mesurées. caractérisée en ce qu'il comporte: des moyens pour déterminer par mesure ou observation d'une série de grandeurs représentatives de la structure stratigraphique de la zone en différents points; et - un ensemble de traitement informatique programmé pour réaliser une implémentation d'un modèle initial direct dépendant d'un ensemble de paramètres représentatifs d'un processus de formation géologique de la zone, cet ensemble de traitement incluant - des moyens pour déterminer au moyen du modèle initial, des valeurs prises par les grandeurs aux différents points de mesure ou d'observation - des moyens pour calculer à partir du modèle direct. des gradients de chaque grandeur par rapport aux paramètres à ajuster - des moyens d'optimisation du modèle en minimisant par itérations successives une fonction critère dépendant des dits paramètres, jusqu'à réduction suffisante des écarts, par le calcul de la sensibilité du modèle direct aux paramètres à ajuster; et des moyens pour produire une représentation stratigraphique de la zone souterraine.