Abstract:
A method of gradual deformation of representations or realizations, generated by sequential simulation, not limited to a Gaussian stochastic model of a physical quantity z in a meshed heterogeneous medium, in order to adjust the model to a set of data relative to the structure or the state of the medium which are collected by previous measurements and observations. The method comprises applying a stochastic model gradual deformation algorithm to a Gaussian vector with N mutually independent variables which is connected to a uniform vector with N mutually independent uniform variables by a Gaussian distribution function so as to define realizations of the uniform vector, and using these realizations to generate representations of the physical quantity z that are adjusted to the data.

Description:
BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present invention relates to a method of gradual deformation of representations or realizations, generated by sequential simulation, of a model of a heterogeneous medium which is not limited to a Gaussian stochastic model, based on a gradual deformation algorithm of Gaussian stochastic models. 
     2. Description of the Prior Art 
     In French patent application 98/09,018 a method is described which gradually deforms a stochastic (Gaussian type or similar) model of a heterogeneous medium such as an underground zone, constrained by a set of parameters relative to the structure of the medium. This method comprises drawing a number p (p=2. for example) of realizations (or representations) independent of the model or of at least part of the selected model of the medium from all the possible realizations and one or more iterative stages of gradual deformation of the model by performing one or more successive linear combinations of p independent initial realizations and then composite realizations are successively obtained possibly with new draws, etc., the coefficients of this combination being such that the sum of their squares is 1. 
     Gaussian or similar models are well-suited for modelling continuous quantity fields but are therefore ill-suited for modelling zones crossed by fracture networks or channel systems for example. 
     The most commonly used geostatistical simulation algorithms are those referred to as sequential simulation algorithms. Although sequential simulation algorithms are particularly well-suited for simulation of Gaussian models, these algorithms are not in principle limited to this type of model. 
     A geostatistical representation of an underground zone is formed for example by subdividing thereof by a network with N meshes and by determining a random vector with N dimensions Z=(Z 1 , Z 2 , . . . Z N ) best corresponding to measurements or observations obtained on the zone. As shown for example by Johnson, M. E.; in “Multivariate Statistical Simulation”; Wiley &amp; Sons, New York, 1987, this approach reduces the problem of the creation of an N-dimensional vector to a series of N one-dimensional problems. Such a random vector is neither necessarily multi-Gaussian nor stationary. Sequential simulation of Z first involves the definition of an order according to which the N elements (Z 1 , Z 2 , . . . Z N ) of vector Z are generated one after the other. Apart from any particular case, it is assumed that the N elements of Z are generated in sequence from Z 1  to Z N . To determine a value of each element Z 1 , (i=1, . . . , N), the following operations have to be carried out: 
     a) building the distribution of Z i  conditioned by (Z 1 , Z 2 . . . Z i-1 ) F c (Z i )=P(Z i ≦/Z 1 , Z 2 , . . . Z i-1 ); and 
     b) determining a value of Zi from distribution F c  (Z i ). 
     In geostatistical practice, sequential simulation is frequently used to generate multi-Gaussian vectors and non-Gaussian indicator vectors. The main function of sequential simulation is to determine conditional distributions F c  (Z i ) (i=1, . . . , N). Algorithms and softwares for estimating these distributions are for example described in: 
     Deutsch, C. V. et al. “GSLIB (Geostatistical Software Library) and Useres Guide”; Oxford University Press, New York, Oxford 1992. 
     Conceming determining the values from distribution F c  (Z i ), there also is a wide set of known algorithms. 
     The inverse distribution method is considered by means of which a realization of Z i ;  i =F c   −1 (u l ) is obtained, where u 1  is taken from a uniform distribution between 0 and 1. A realization of vector Z therefore corresponds to a realization of vector U whose elements U 1 , U 2 , . . . , U N , are mutually independent and evenly distributed between 0 and 1. 
     It can be seen that a sequential simulation is an operation S which converts a uniform vector U=(U 1 , U 2 , . . . U N ) to a structured vector Z=(Z 1 , Z 2 , . . . , Z N ): 
     
       
           Z=S ( U )  (1). 
       
     
     The problem of the constraint of a vector Z to various types of data can be solved by constraining conditional distributions F c (Z i ) (i=1, . . . , N) and/or uniform vector U=(U 1 , U 2 , . . . , U N ). 
     Recent work on the sequential algorithm was focused on improving the estimation of conditional distributions F c  (z i ) by geologic data and seismic data integration. An article by Zhu, H. et al: “Formatting and Integrating Soft Data: Stochastic imaging via the Markov-Bayes Algorithm” in Soares, A., Ed. GeostaUstics Troia 92, vol.l: KiuwerAcad. Pubi., Dordrecht, The Netherlands, pp.1-12, 1993 is an example. 
     However, this approach cannot be extended to integration of non-linear data such as pressures from well tests and production records, unless a severe linearization is imposed. Furthermore, since any combination of uniform vectors U does not give a uniform vector, the method for gradual deformation of a stochastic model developed in the aforementioned patent application cannot be directly applied within the scope of the sequential technique described above. 
     The method according to the invention thus allows making the two approaches compatible, that is to extend the formalism developed in the aforementioned patent application to gradual deformation of realizations, generated by sequential simulation, of a model which is not limited to a Gaussian stochastic model. 
     SUMMARY OF THE INVENTION 
     The method allows gradual deformation of a representation or realization, generated by sequential simulation, of a model which is not limited to a Gaussian stochastic model of a physical quantity z in a heterogeneous medium such as an underground zone, in order to constrain the model to a set of data collected in the medium by previous measurements and observations relative to the state or the structure thereof. 
     The method of the invention comprises applying an algorithm of gradual deformation of a stochastic model to a Gaussian vector (Y) having a number N of mutually independent variables that is connected to a uniform vector (U) with N mutually independent uniform variables by a Gaussian distribution function (G), so as to define a chain of realizations u(t) of vector (U), and using these realizations u(t) to generate realizations z(t) of this physical quantity that are adjusted in relation to the (non-linear) data. 
     According to a first embodiment, the chain of realizations u(t) of uniform vector (U) is defined from a linear combination of realizations of Gaussian vector (Y) whose combination coefficients are such that the sum of their squares is one. 
     According to another embodiment, gradual deformation of a number n of parts of the model representative of the heterogeneous model is performed while preserving the continuity between these n parts of the model by subdividing uniform vector (U) into a number n of mutually independent subvectors. 
     The method of the invention finds applications in modeling underground zones which generates representations showing how a certain physical quantity is distributed in an underground zone (permeability z for example) and which is compatible in the best manner with observed or measured data: geologic data, seismic records, measurements obtained in wells, notably measurements of the variation with time of the pressure and of the flow rate of fluids from a reservoir, etc. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     Other features and advantages of the method according to the invention will be clear from reading the description hereafter of a non limitative example, with reference to the accompanying drawings wherein: 
     FIG. 1 shows the medial layer of a realization of a facies model used as a reference, generated by sequential simulation of indicatrices, 
     FIG. 2 shows the variation with time of the pressure obtained in a well test for the reference model, 
     FIGS. 3A to  3 E respectively show five initial realizations of the medial layer of a reservoir zone, constrained only by the facies along the well, 
     FIGS. 4A to  4 E respectively show, for these five realizations, the bottomhole pressure curves in the reference model compared with those corresponding to the initial models, 
     FIGS. 5A to  5 E respectively show five realizations of the medial layer of the facies model conditioned to the facies along the well and adjusted in relation to the pressure curve obtained by well tests, 
     FIGS. 6A to  6 E respectively show, for the five realizations, the bottomhole pressure curves in the reference model compared with those corresponding to the adjusted models, 
     FIGS. 7A to  7 E respectively show how the objective functions respectively corresponding to these five examples vary with the number of iterations, 
     FIGS. 8A to  8 E show the gradual deformations generated by an anisotropy coefficient change on a three-facies model generated by sequential simulation of indicatrices, and 
     FIGS. 9A to  9 E show the local gradual deformations of a three-facies model, generated by sequential simulation of indicatrices. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     A study zone is considered that is subdivided by an N-mesh grid. Realizations or representations of a stochastic model of a certain physical quantity z representing for example the permeability of the formations in the zone are attempted to be made. The model that is sought must adjust to data obtained by measurements or observations at a certain 
     Adjustment of a stochastic model to non-linear data can be expressed as an optimization problem. The quantity f obs=  (f 1   obs , f 2   obs , f 3   obs . . . f p   obs ) designates the vector of the nonlinear data observed or measured in the studied medium (the reservoir zone), and the quantity f=(f 1 , f 2 , f 3  . . . f p ) is the corresponding vector of the responses of the stochastic model of the permeability Z=(Z 1 , Z 2 , . . . , Z N ,). The problem of constraining the stochastic model of Z by observations consists in generating a realization z of Z which reduces to a rather low value an objective function that is defined as the sum of the weighted rms errors of the responses of the model in relation to the observations or measurements in the reservoir zone, i.e.:        O   =       1   2                       ∑     i   =   1     p                       ω   i                       (       f   1     -     f   1   obs       )     2                                  
     where ω i , represents th weight assigned to respons f i  Functions f i  (i=1, 2, . . . , p) and objective function O are functions of vector Z. This presents an optimization problem of dimension N. 
     In order to extend the formalism developed in the aforementioned patent application to the gradual deformation of realizations generated by, but not necessarily limited to Gaussian sequential simulation, a starting point is from a Gaussian vector with N variables Y i , with i=1, 2, . . . , N, mutually independent, of zero mean and of variance equal to 1, and N mutually independent uniform variables U i , U 2 , U 3 , . . . U N  are defined by: 
     
       
         U i =G(Y i )∀i=1, 2, . . . , N 
       
     
     where G represents the standardized Gaussian distribution function. 
     Assuming this to be the case, the gradual deformation algorithm developed within a Gaussian frame is applied to the Gaussian vector Y=(Y 1 , Y 2 , . . . , Y N ) in order to build a continuous chain of realizations of uniform vector U (U 1 , U 2 , . . . , U N ). Given two independent realizations Y A  and Y b  of Y, the chain of realizations u(t) of vector U obtained with the following relation is defined: 
     
       
           u ( t )= G ( y   a , cos t+y b , sin  t )  (2). 
       
     
     For each t, u(t) is a realization of vector U. A vector z(t) which is, for each t, a realization of random vector Z is then obtained by sampling of the conditional distribution F c  (z i ) (i=1, 2, . . . , N) using the elements of vector u(t). Parameter t can consequently be adjusted as in the Gaussian case so as to adjust z(t) to non-linear data. This procedure is iterated until satisfactory adjustment is obtained. 
     Adjustment of a Facles Model to Pressure Data Obtained by Means of Well Tests 
     In order to illustrate application of the stochastic optimization method defined above, adjusting a stochastic reservoir model to pressure data obtained by means of well tests is attempted. Building of the reservoir model is derived from a real oil formation comprising three types of facies: two reservoir facies of good quality (facies  1  and  2 ) and a reservoir facies of very bad quality (facies  3 ). Table 1 defines the petrophysical properties of the three facies: 
     
       
         
               
               
               
               
               
               
             
               
               
               
               
               
               
             
           
               
                   
                   
               
               
                   
                 K x  (md) 
                 K y  (md) 
                 K z  (md) 
                 φ (%) 
                 c t  (10 −5  bar −1 ) 
               
               
                   
                   
               
             
             
               
                   
               
             
          
           
               
                 Facies 1 
                 10 
                 10 
                 10 
                 17 
                 2.1857 
               
               
                 Facies 2 
                 1 
                 1 
                 1 
                 14 
                 2.0003 
               
               
                 Facies 3 
                 0.1 
                 0.1 
                 0.001 
                 9 
                 1.8148 
               
               
                   
               
             
          
         
       
     
     In order to represent the specific facies distribution of the oil formation, a binary realization is first generated to represent fades  3  and its complement. Then, in the complementary part of facies  3 , another binary realization independent of the first one is generated to represent facies  1  and  2 . The formation is discretized by means of a regular grid pattern of 60×59×15 blocks 15 m×15 m×1.5m in size. An exponential variogram model is used to estimate the conditional distributions. The main anisotropy direction is diagonal in relation to the grid pattern. The ranges of the variogram of facies  3  in the three anisotropy directions are 300 m, 80 m and 3 m respectively. The ranges of the variogram of facies  1  and  2  in the three anisotropy directions are 150 m, 40 m and 1.5 m respectively. The proportions of facies  1 ,  2 ,  3  are 6%, 16% and 78% respectively. 
     A well test has been carried out by means of a finite-difference well test, simulator as described by: Blanc, G. et al: “Building Geostatistical Models Constrained by Dynamic Data—A Posteriori Constraints” in SPE 35478, Proc. NPF/SPE European 3D Reservoir Modelling Conference, Stavenger, Norway, 1996. 
     The medial layer of a realization used as the reference model for this validation can be seen in FIG.  1 . The section of the well that has been drilled runs horizontally through the medial layer of the reservoir model along axis x. The diameter of the well is 7.85 cm, the capacity of the well is zero and the skin factors of facies  1 ,  2  and  3  are 0, 3 and 50 respectively. The synthetic well test lasts for 240 days with a constant flow rate of 5 m 3 /day so as to investigate nearly the entire oil field. FIG. 2 shows the pressure variation with time. 
     The objective was to build realizations of the facies model constrained by the facies encountered along the well and by the pressure curve obtained during well testing. The objective function is defined as the sum of the rms differences between the pressure responses of the reference model and the pressure responses of the realization. Since the dynamic behavior of the reservoir model is mainly controlled by the contrast between the reservoir facies of good and bad quality, the binary realization used to generate facies  1  and  2  has been fixed first and only the binary realization used to generate facies  3  has been deformed for pressure data adjustment. 
     The pressure responses resulting from the well tests for the five realizations of FIGS. 3A to  3 E are different from that of the reference model, as shown in FIGS. 4A to  4 E. Starting respectively from these 5 independent realizations, by using the iterative adjustment method above, after several iterations, five adjusted are obtained realizations (FIGS. 5A to SE) for which the corresponding pressure curves are totally in accordance with those of the reference model, as shown in FIGS. 6A to  6 E. 
     Gradual Deformation in Relation to the Structural Parameters 
     In many cases, sufficient data for deducing the structural parameters of the stochastic model: mean, variance, covariance function, etc., are not available. These structural parameters are often given in terms of a priori intervals or distributions. If their values are wrong, it is useless to seek a realization adjusted to non-linear data, It is therefore essential for applications to be able to perform a gradual deformation of a realization with simultaneous modification of random numbers and structural parameters. The sequential simulation algorithm defined by equation (1) makes possible changing, simultaneously or separately, a structural operator S and a uniform vector U. FIGS. 8A to  8 E show the gradual deformations obtained for a fixed realization of uniform vector U when the anisotropy coefficient is changed. 
     Local or Regionalized Gradual Deformation 
     When the observations are spread out over different zones of a studied formation, an adjustment using global deformation would be ineffective because the improvement obtained in a zone could deteriorate the improvement in another zone. It is therefore preferable to apply gradual deformations zone by zone. Consider a subdivision of vector U into a certain number n of mutually independent subvectors U 1 , U 2 , . . . , U n , which allows performing their gradual deformation individually. Separate application of the gradual deformation algorithm to each subvector U 1 , U 2 . . . , U n  allows obtaining a function of dimension n of uniform vector U:          U                   (       t   1     ,     t   2     ,   …              ,     t   n       )       =       [             U   1                     (     t   1     )                   U   2                     (     t   2     )               ⋮               U   n                     (     t   n     )             ]     =     [           G                   (         Y   a   1                   cos                   t   1       +       Y   b   1                   sin                   t   1         )                 G                   (         Y   a   2                   cos                   t   2       +       Y   b   2                   sin                   t   2         )               ⋮             G                   (         Y   a   n                   cos                   t   n       +       Y   b   n                   sin                   t   n         )             ]                              
     where Y i   a  and Y i   b  for any i=1, 2, . . . , n, are independent Gaussian subvectors. For a set of realizations of V i   a , and v i   b , a problem of optimization of n parameters t 1 , t 2 , . . . , t n  is solved to obtain a realization that improves or maintains the data adjustment This procedure can be iterated until satisfactory adjustment is obtained. 
     Gradual local deformations thus allow significant improvement of the adjustment speed in all the cases where measurements or observations are spread out over different zones of the medium. 
     The effect of this gradual local deformation on the three-facies model of FIGS. 9A to  9 E can be clearly seen therein where only the enclosed left lower part is affected. 
     The method according to the invention can be readily generalized to gradual deformation of a representation or realization of any stochastic model since generation of a realization of such a stochastic model always comes down to generation of uniform numbers.