FR 2862765 A1 20050527 FR 0313694 A 20031120 Hu, L.Y., Blanc, G. and Noetinger, B., 2001, Graduai deformation and iterative calibration of sequential stochastic simulations, Math. Geol., 33(4). Hu, L.Y., 2003, History matching of object-based stochastic reservoir models. Paper SPE 81503. - Le Ravalec, M., Noetinger, B. and Hu, L.Y., 2000, The FFT moving average (FFF1'-MA) generator: an efficient numerical method for generating and conditioning Gaussian simulations, Math. Geol., 32(6). Le Ravalec, M., Hu, L.Y. and Noetinger, B., 2001, Stochastic Reservoir modeling constrained to dynamic data: local calibration and inference of the structural 10 parameters, SPE Journal, March. - Roggero F. and Hu, L.Y., 1998, Graduai deformation of continuous geostatistical models for history matching, Paper SPE 49004, - Caers, J., 2002, Geostatistical history matching under training-image based geological model constraints, Paper SPE 77429 - Journel, A., 2002, Combining knowledge from diverse information sources: an alternative to Bayesian analysis, Math. Geol., 34(5) - Hoffman, B.T., Caers, J., 2003, Geostatistical history matching using a regional probability perturbation method, Paper SPE 84409. 1.2 Méthode de déformation graduelle Différentes méthodes basées sur un schéma de déformation graduelle, ont fait l'objet de brevets ou demandes de brevets suivants du demandeur: FR 2.780.798 (US 18 95), FR 2.795.841, FR 2.821.946, FR 02/13.632 ou FR 03/02.199, ainsi que de plusieurs des publications rappelées plus haut. Les réalisations (images géostatistiques) de modèles stochastiques gaussiens, ou apparentés, sont souvent utilisées pour représenter la distribution spatiale de certaines propriétés physiques, comme la perméabilité ou la porosité, dans les réservoirs souterrains. On recourt alors couramment aux méthodes inverses pour contraindre ces réalisations à des données dont elles dépendent de façon non linéaire. C'est le cas notamment en hydrologie 2862765 3 ou en ingénierie pétrolière. Ces méthodes reposent sur la minimisation d'une fonction objectif, encore appelée fonction coût, qui mesure l'écart entre les données mesurées sur le terrain et les réponses correspondantes simulées numériquement pour les réalisations représentant le milieu que l'on cherche à caractériser. Le but affiché est d'identifier les réalisations associées aux valeurs de fonctions objectifs les plus faibles, c'est-à-dire les réalisations les plus cohérentes au regard des données. La méthode de déformation graduelle a été introduite par les brevets précités et les publication suivantes précédemment citées: Hu, 2000a; Le Ravalec et al., 2000. Cette technique de paramétrisation géostatistique permet de modifier graduellement les réalisations à partir d'un nombre restreint de paramètres. Elle est particulièrement adaptée aux problèmes de minimisation, car appliquée aux réalisations, elle induit une variation continue et régulière de la fonction objectif. De fait, la minimisation peut être conduite à partir des techniques les plus avancées, à savoir les techniques de gradients. La méthode de déformation graduelle s'est avérée efficace pour contraindre des modèles de réservoirs pétroliers à des données de production (cf. Roggero et Hu, 1998; Le Ravalec et al., 2001). Par la suite, la méthode de déformation graduelle, initialement mise en place pour des modèles gaussiens, a été étendue à tout type de modèles stochastiques (cf. Hu, 2000b; Hu Blanc and Noetinger, 2001; Hu, 2003). 1.2-1 Rappels La méthode de déformation graduelle s'applique par exemple aux fonctions aléatoires multigaussiennes. Soient, par exemple, deux fonctions aléatoires indépendantes, Yl(x) et Y2(x), multigaussiennes et stationnaires d'ordre 2 où .x représente le vecteur position. On suppose que ces deux fonctions ont mêmes moyennes et variances, à savoir 0 et 1, et même fonction de covariance. On construit alors une nouvelle fonction aléatoire Y(t) en combinant Yl et Y2 suivant l'expression: Y(t) = Y, cos(t)+ Y2 sin(t). (1) On peut montrer que quel que soit t, Y a même moyenne, variance et modèle de covariance que Yl et Y2. En outre, Y(t) est aussi une fonction aléatoire multigaussienne, car elle est la somme de deux fonctions aléatoires multigaussiennes. Suivant ce principe de combinaison, on peut, à partir de deux réalisations indépendantes yl et y2 de YI et Y2, construire une chaîne de réalisations y(t) qui dépendent uniquement du paramètre de déformation t. L'idée de base des processus de minimisation par déformation graduelle est d'explorer cette chaîne de réalisations et de déterminer le paramètre de déformation qui fournit la réalisation la plus compatible avec les données mesurées sur le terrain, c'est-à-dire les pressions, les débits de production, les temps de percée, etc. Comme l'exploration d'une unique chaîne de réalisation ne permet pas en général d'identifier une réalisation fournissant une fonction objectif suffisamment faible, on itère le processus de recherche. On combine alors la réalisation optimale déterminée pour la 1ere chaîne avec une nouvelle réalisation indépendante de Y2 et on en déduit une nouvelle chaîne de réalisations dont l'exploration peut fournir une réalisation qui diminue davantage la fonction objectif, etc. Au lieu de combiner deux réalisations indépendantes, on peut généraliser la relation (1) en combinant plusieurs réalisations indépendantes afin d'augmenter le degré de liberté pour la déformation d'une réalisation comme décrit dans les publications (cf. Roggero et Hu, 1998), et (Hu, 2000) et les brevets du demandeur déjà cités. L'algorithme de déformation graduelle précédemment décrit dans le cadre gaussien s'applique de façon plus générale à tout type de modèles stochastiques selon le principe suivant défini dans la publication déjà citée référencée (Hu, 2000b). Toute simulation stochastique peut être considérée comme une opération déterministe qui transforme une série de nombres U1, U2, U3,... UN,... indépendants les uns des autres et répartis selon la loi uniforme, en une série de nombres structurés selon le modèle stochastique. Ces nombres uniformes peuvent être transformés en variables gaussiennes Y1, Y2, Y3,...YN,... mutuellement indépendantes de moyenne nulle et de variance égale à l'unité, par la relation: 1.. =G-'(U;) i=1,2,3... où G-1 désigne la fonction de répartition gaussienne standardisée. Ceci posé, on va appliquer l'algorithme de déformation graduelle développé dans un cadre gaussien, au vecteur gaussien Y = (Y1, Y2, YN...), pour construire une chaîne continue de réalisations du vecteur uniforme U = (U1, U2, UN...). Etant donné deux réalisations indépendantes yl et y2 de Y, on définit la chaîne de réalisations u(t) du vecteur U que l'on obtient selon la relation: u(t) = G(y, cost + y2 sin t) (Ibis) Pour chaque t, u(t) est une réalisation du vecteur U. En conséquence, on peut ajuster le paramètre t comme dans le cas gaussien de façon à caler une réalisation à des données non linéaires (dynamiques). Cette procédure est itérée jusqu'à l'obtention un calage satisfaisant. L'application de l'algorithme de déformation graduelle (1) aux vecteurs U = (U1, U2, UN...) et Y = (Y1, Y2, YN...) (appelés aussi bruits blancs) permet non seulement d'étendre la méthode de déformation à tout type de modèle stochastiques (modèles continus, en faciès ou en objets, etc.) mais aussi à la déformation locale ou par rapport aux paramètres structuraux d'un modèle comme développé dans la publication (Hu, 2000) déjà citée. 1.3 Méthode de perturbation de probabilités La méthode de perturbation de probabilités (dite PPM pour Probability Perturbation Method ) est décrite dans le document Caers (2002) déjà cité. Elle s'inscrit aussi dans le cadre de la construction d'un modèle d'un milieu hétérogène contraint par des données dynamiques par optimisation d'une fonction objectif. Il s'agit d'une technique de paramétrisation, différente de celle de déformation graduelle, permettant aussi de modifier graduellement les réalisations à partir d'un nombre restreint de paramètres. Elle est basée sur une formule approximative de calcul de probabilité conditionnelle (cf. Journel, 2002 déjà cité) et un algorithme de paramétrisation de probabilités discrètes (Caers, 2002). Jusqu'à présent, cette technique est limitée aux modèles en faciès (modèles catégoriels) générés par simulation séquentielle. 1.3-1 Rappels La méthode de perturbation de probabilités s'applique aux modèles en faciès (modèles catégoriels) générés par simulation séquentielle. Soit I(x) la fonction d'indicatrice du faciès K d'un modèle stochastique: 1 xE K I(x) 0 xE K Désignons par A l'événement où faciès K a lieu au noeud x, c'est-à-dire A=[I(x)=1 j; par B l'événement des faciès dans un voisinage du noeud x défini par l'algorithme de simulation séquentielle; et par D l'événement des données de production. Selon une formule approximative proposée par Journel (2002), on peut calculer la probabilité de A sachant B et D en fonction de la probabilité marginale de A, la probabilité conditionnelle de A sachant B et la probabilité conditionnelle de A sachant D: P(A/ B, D) = [1- P(A)]P(Al B)P(AI D) [1 P(A)]P(A/ B)P(A/ D) + P(A)[1 P(A/ B)] [1 P(Al D)] P(A) est connue a priori, P(A/B) est calculée par l'algorithme de simulation 10 séquentielle, alors que P(A/D) est définie par: P(A/ D) = (1 r)i(0) (x) + rP(A) (4) où i(0) (x) désigne l'indicatrice du faciès K de la réalisation initiale du modèle stochastique et r le paramètre de perturbation qui varie entre 0 et 1. Quand r=0, on a P(A/ D) = i(0)(x) , ce qui donne P(A/B,D) = i(0) (x). Par conséquent, on retrouve la réalisation initiale lors d'une génération d'une réalisation du modèle stochastique selon la probabilité P(A/B,D) et son complémentaire. Par contre, quand r=1, on a P(A/ D) = P(A) , on en déduit P(A/ B, D) = P(A/ B) . Par conséquent, la simulation selon P(A/B,D) conduit à une réalisation indépendante de celle initiale j (0) (x) . On s'aperçoit qu'en fixant le germe aléatoire et en faisant varier le paramètre r de 0 à 1, on modifie de façon graduelle la réalisation initiale i(0) (x) vers une autre réalisation i(l) (x) indépendante de la première. Comme dans le cas de la méthode de déformation graduelle, nous pouvons ainsi minimiser une fonction objectif par rapport au paramètre r afin d'améliorer le calage d'un modèle aux données de production. Nous pouvons aussi itérer cette procédure en partant d'une réalisation obtenue précédemment et ainsi de suite jusqu'à l'obtention d'un calage satisfaisant. (2) (3) La technique précédente de paramétrisation' selon la méthode de perturbation de probabilités constitue donc une alternative à la méthode de déformation graduelle. Néanmoins, cette alternative est limitée dans le cas (au cas ?) des modèles en faciès générés par simulation séquentielle et elle nécessite de recourir à la formule approximative (3). LA METHODE SELON L'INVENTION La méthode selon l'invention a pour objet la formation d'un modèle stochastique optimal représentatif de la répartition spatiale dans une zone souterraine hétérogène, des grandeurs physiques, telles que la perméabilité et la porosité, qui soit contraint par des données dynamiques mesurées, représentatives de déplacements de fluides dans le milieu, et en imposant une continuité dans l'évolution de réalisations successives du modèle vers une réalisation optimale. Elle comporte une optimisation par un processus itératif de minimisation d'une fonction objectif mesurant à chaque itération, l'écart entre les données dynamiques mesurées et des données dynamiques simulées au moyen d'un simulateur d'écoulement, obtenues à partir d'une réalisation interpolée entre une réalisation de référence (la réalisation initiale ou celle obtenue à l'issue de l'itération précédente) et une autre réalisation indépendante, par ajustement d'un paramètre, le processus d'ajustement itératif étant poursuivi jusqu'à l'obtention d'une réalisation optimale du modèle stochastique bien ajustée aux données dynamiques et/ou par déformation graduelle du modèle en effectuant une combinaison linéaire de réalisations indépendantes, les coefficients de cette combinaison étant tels que la somme de leurs carrés est égale à 1. La méthode se distingue en ce que, à chaque itération, on construit une réalisation intermédiaire entre la réalisation de référence et une autre réalisation indépendante selon un algorithme de perturbation de distributions continues. Suivant un mode de mise en oeuvre, on construit la réalisation intermédiaire, en changeant un paramètre (r) qui modifie la distribution de probabilité. Suivant un mode de mise en oeuvre, on construit la réalisation intermédiaire, en changeant un paramètre (t) commandant une déformation graduelle du modèle. La méthode peut s'appliquer par exemple à la formation de modèles générés par 30 simulations séquentielles gaussiennes ou à partir de bruits blancs gaussiens. En d'autres termes, la méthode selon l'invention constitue une extension de la méthode de perturbation de probabilités discrètes (limitée aux modèles catégoriels générés par simulation séquentielle). Elle peut s'appliquer aux bruits blancs gaussiens sans qu'aucune approximation ne soit nécessaire. Par conséquent, la méthode s'applique à tout type de modèle stochastique (y compris les modèles en faciès ou en objet). Quand elle est appliquée à un modèle en faciès via un bruit blanc, elle se distingue de la méthode de perturbation de probabilités décrite précédemment, en ce qu'elle ne nécessite aucune approximation. Quand elle est appliquée à tout type de modèle stochastique via des bruits blancs, la méthode constitue une alternative à la méthode de déformation graduelle (applicable aussi à tout type de modèles stochastiques). PRESENTATION SUCCINCTE DES FIGURES Les caractéristiques et avantages de la méthode selon l'invention, apparaîtront plus clairement à la lecture de la description détaillée ciaprès en se référant aux dessins 15 annexés où : - la figure 1 montre un exemple de perturbation de distribution pour plusieurs valeurs du paramètre de perturbation r en cas d'une simulation séquentielle gaussienne; et - la figure 2 montre un exemple d'évolution de deux nombres gaussiens vers deux autres nombres gaussiens en fonction du paramètre de perturbation r en cas de perturbation 20 d'un bruit blanc gaussien. DESCRIPTION DETAILLEE DE LA METHODE La méthode selon l'invention permet on l'a vu d'optimiser un modèle stochastique représentatif de la répartition spatiale de grandeurs telles que la perméabilité, la porosité, etc. dans une zone souterraine hétérogène, qui soit contraint par des données dynamiques mesurées telles que des données représentatives de déplacements de fluides dans le milieu, obtenues par des tests de puits forés au travers de la zone par exemple, et en imposant une continuité dans l'évolution de réalisations (états ou images géostatistiques) du modèle. Elle comporte l'application d'une technique d'optimisation itérative suivant laquelle, à partir d'une réalisation initiale du modèle, on minimise à chaque itération une fonction objectif mesurant l'écart entre des données dynamiques effectivement mesurées et des données dynamiques simulées, au moyen d'un simulateur d'écoulement, sur une réalisation interpolée, entre une réalisation de référence (réalisation initiale en début du processus ou réalisation courante obtenue à la fin de l'itération précédente) et une autre réalisation indépendante. Cette minimisation est conduite en ajustant un paramètre de perturbation, selon un algorithme que l'on va décrire, de perturbation de distributions continues. Méthode de perturbation de distributions continues La méthode existante de perturbation de probabilités discrètes selon la relation (4) rappelée plus haut, est étendu, dans le cadre de la présente invention, aux modèles continus. Soit Y(x) une fonction aléatoire (non nécessairement gaussienne) admettant la fonction de répartition marginale F(y). Désignons par A l'événement où Y(x) est inférieure à un certain seuil y au noeud x, que l'on traduit par l'expression A= {Y(x) On écrit: P(A) = F(y) P(A/ B) = FB (y) et on définit: P(A/ D) = FD (y) = (1 r)ly,yo + rF(y) (5) où yo est la réalisation initiale de Y(x) au noeud x, et où : 1 y > yo 0 y yo est une fonction d'indicatrice de la réalisation yo. Moyennant la formule approximative (3), on peut obtenir la fonction de répartition FB,D (y) de Y(x) sachant B et D en fonction de F(y), FB (y) et FD (y) . On vérifie que quand r=0, FB D (y) = 1y,yo. Par conséquent, on retrouve la réalisation initiale lors d'une (6) génération d'une réalisation du modèle stochastique selon la loi FB D (y) . Quand r=1, on obtient FB D (y) = FB (y) . Par conséquent, la simulation selon FB D (y) conduit à une réalisation de Y(x) indépendante de la réalisation initiale yo (x) . On s'aperçoit qu'en fixant le germe aléatoire et en faisant varier le paramètre r de 0 à 1, on modifie de façon graduelle 5 la réalisation initiale yo (x) vers une autre réalisation gaussienne yl (x) indépendante de la première. Application aux simulations séquentielles gaussiennes La méthode de perturbation des distributions continues s'applique en particulier aux modèles continus générés par simulation séquentielle gaussienne. Soit maintenant Y(x) une 10 fonction aléatoire multigaussienne. Nous avons: P(A) = G(y) P(A/B)=FB(y)=G(y m' P(A/ D) = FD (y) = (1 r)1y,yp + rG(y) (7) où G est la fonction de répartition gaussienne centrée et réduite, m et 6-2 sont 15 respectivement le krigeage et la variance du krigeage de Y(x) à partir des données B. Moyennant la formule (3), nous avons: r[l G(y)]G y rG(y) (1 r)G y m y yo (8) 2862765 11 simulation de loi mono-variable. La figure 1 montre un exemple de la distribution (8) pour plusieurs valeurs du paramètre de perturbation r. Application aux bruits blancs gaussiens Dans le cas particulier où Y(x) est un bruit blanc gaussien (fonction (ou vecteur) 5 aléatoire gaussienne sans corrélation spatiale), l'événement B n'apporte aucune information à l'événement A et la fonction de répartition (8) se simplifie: rG(y) Y _ Yo Notons que cette distribution s'obtient aussi sans passer par la formule approximative (3). En effet, comme l'événement B n'apporte aucune information à 10 l'événement A, nous avons directement: rG(y) Y Yo On peut montrer que pour toute valeur de r, l'interpolation, selon l'algorithme de perturbation (9), entre deux bruits blancs gaussiens indépendantes reste un bruit blanc gaussien. La méthode de perturbation de distributions continues appliquée aux bruits blancs n'est plus une méthode approximative, contrairement à son application à d'autres situations et aussi à la méthode de perturbations de probabilités (Caers, 2002) . Pour une valeur de r fixe, l'algorithme de simulation de la loi (9) est le suivant: 1. générer un bruit blanc gaussien standard yo(x), 2. générer un bruit blanc uniforme u(x) si rG[ yo (x)] 1 r + rG[ yo (x)] , alors y(x) = G 1 (u(x)+r l\ On montre que y(x) ainsi généré est un bruit blanc gaussien. La figure 2 montre l'évolution de deux nombres gaussiens vers deux autres nombres gaussiens en fonction du paramètre de perturbation r selon l'algorithme précédent. La paramétrisation selon la relation (9) permet donc de modifier un bruit blanc initial vers un autre bruit blanc indépendant du premier. Cette modification respecte toutes les caractéristiques du bruit blanc gaussien. Comme explicité plus haut dans le rappel des mécanismes de la méthode de déformation graduelle, on peut appliquer la paramétrisation selon la relation (9) en tant qu'alternative à la paramétrisation selon la méthode de déformation graduelle, pour modifier, de façon cohérente, globalement ou localement, tout type de simulation stochastique. Procédure de calage itératif Comme dans le cas de la méthode de déformation graduelle, la méthode de perturbation de distributions continues s'intègre naturellement dans les processus de minimisation. La procédure d'optimisation se traduit par l'exploration successive de chaînes de réalisations continues fabriquées à partir d'une réalisation initiale ou de la réalisation optimale déterminée pour la chaîne de réalisations formée à la fin de l'itération précédente, et d'une seconde réalisation générée aléatoirement pour chaque chaîne. Ce processus de recherche est itéré jusqu'à ce que la fonction objectif soit jugée suffisamment faible et que l'image géostatistique de la zone souterraine qui en résulte, corresponde au mieux aux données dynamiques mesurées. Utilisation conjointe à la méthode de déformation graduelle La méthode selon l'invention peut être utilisée conjointement à la méthode de défolination graduelle. En effet, à partir d'une réalisation initiale ou de la réalisation optimale obtenue à la fin de l'itération précédente, nous pouvons construire une nouvelle chaîne de réalisations en fonction des deux paramètres suivants: le paramètre r de la relation (9) qui modifie la distribution de probabilité et le paramètre t de la relation (Ibis) qui modifie le nombre uniforme pour échantillonner cette distribution. Nous pouvons ensuite ajuster simultanément ces deux paramètres pour améliorer le calage aux données de production, ce qui donne un degré de liberté de plus par rapport à l'utilisation seule de la méthode de perttirbation de distribution ou de la méthode de déformation graduelle. Comme dans le précédent, le processus de recherche est itéré jusqu'à ce que le calage aux données de production soit jugée suffisamment satisfaisant. REVENDICATIONS 1) Méthode pour former un modèle stochastique optimal représentatif de la répartition spatiale dans une zone souterraine hétérogène, des grandeurs physiques, telles que la perméabilité et la porosité, qui soit contraint par des données dynamiques mesurées, représentatives de déplacements de fluides dans le milieu, et en imposant une continuité dans l'évolution de réalisations successives du modèle vers une réalisation optimale, dans laquelle on réalise l'optimisation par un processus itératif de minimisation d'une fonction objectif mesurant à chaque itération, l'écart entre les données dynamiques mesurées et des données dynamiques simulées au moyen d'un simulateur d'écoulement, obtenues à partir d'une réalisation interpolée entre une réalisation de référence et une autre réalisation indépendante, par ajustement d'un paramètre, le processus d'ajustement itératif étant poursuivi jusqu'à l'obtention d'une réalisation optimale du modèle stochastique bien ajustée aux données dynamiques, caractérisée en ce que, à chaque itération, on construit une réalisation intermédiaire entre la réalisation de référence et une autre réalisation indépendante suivant un processus de perturbation de distributions continues et/ou par déformation graduelle du modèle en effectuant une combinaison linéaire de réalisations indépendantes, les coefficients de cette combinaison étant tels que la somme de leurs carrés est égale à 1. 2) Méthode selon la revendication 1 dans laquelle on construit la réalisation 20 intermédiaire, en changeant un paramètre (r) qui modifie la distribution de probabilité. 3) Méthode selon la revendication 1 dans laquelle on construit la réalisation intermédiaire, en changeant un paramètre (t) commandant une déformation graduelle du modèle. 4) Application de la méthode selon la revendication 1 à la formation de modèles 25 générés par simulations séquentielles gaussiennes. 5) Application de la méthode selon la revendication 1 à la formation de modèles générés à partir de bruits blancs gaussiens.