5833 Simple Zonation and Principal Component Analysis for Speeding Up Porosity and Permeability Estimation from 4D Seismic an
M. Dadashpour* (NTNU), D. Echeverría Ciaurri (Stanford University), T. Mukerji (Stanford University), J. Kleppe (NTNU) & M. Landrø (NTNU)
SUMMARY
Simple zonation and a mathematical transformation based on principal component analysis are used in a distributed computing environment for the estimation of porosity and permeability from production and 4D-seismic data, in the form of zero offset amplitudes and amplitude versus offset gradients. The parameter estimation problem is formulated as a least-squares minimization. The GaussNewton technique is used for this purpose which requires gradient information and derivatives are estimated numerically. The results show that with only a bound constraint the solution space cannot be properly limited, and therefore that the solutions found might not be geologically consistent. Although simple gradzone analysis can speed up the optimization, the resulting spatial distributions of porosity and permeability are not acceptable from a geological point of view. Principal component analysis gives us the opportunity to not only accelerate the optimization but also, by incorporating some more complex spatial constraints, to force the parameter estimation to honour geology. By considering both distributed computing and parameter space reduction techniques it is possible to apply the methodology presented to larger problems with practical relevance.
71st EAGE Conference & Exhibition — Amsterdam, The Netherlands, 8 - 11 June 2009
Introduction The extensive computation time required for optimization is one of the main challenges in every parameter estimation problem. This becomes important when we encounter problems with a large number of unknown parameters and it will be even more critical when estimation of gradients is needed. Several methods have been introduced for speeding up the optimization. One efficient way for some optimization algorithms is to introduce a parallel computing framework (Ouenes et al. 1995, Ken-Ichi et al. 2007). Parallelization is common in global minimization algorithms but it can be also effective for gradient computation. The gradient of the data mismatch can be calculated by means of adjoints (Ruijian et al. 2003, Dong and Oliver 2008). This can reduce significantly the total time required for optimization, but unfortunately such algorithms require explicit knowledge of the simulator source code. This problem becomes more evident when 4D seismic computations are added to the forward model. Another method for speeding up the optimization is reparameterization of the reservoir to reduce the number of parameters. This reparameterization can be done by defining pilot-points (Marsily et al. 1984, Arenas et al. 2001) and then reduce the parameter space by kriging. It can also be done by simple zonation or by mathematical transforms. Gradzone analysis (Harb 2004) avoids the calculation of sensitivities at all grid nodes. Reynolds et al. (1996) use a subspace method for reducing the matrix problem size, and compare it with reparameterization based on spectral decomposition. Principal component analysis (Shah et al. 1978) and it generalization kernel principal component analysis (Sarma et al.2007) are two other methods for reducing the parameter space in an optimization problem. These methods are based on the singular value decomposition of a covariance matrix (or sensitivity matrix) that incorporates previous statistical knowledge of the parameters to estimate. This project considers two different reparameterization techniques based on gradzone and principal component analysis for generating reservoir descriptions conditioned to 4D seismic and production data. Gradzone Analysis Gradzone analysis (GZA) is a method that, based on sensitivities, groups parameters to define a region. This requires first and second derivative information. The GZA procedure is as follows (Gosselin and Berg 2001): 1. compute the Hessian matrix, 2. make a spectral decomposition for obtaining the eigenvectors, 3. define a first threshold, 4. group together all parameters which associated eigenvectors higher than that threshold (two different groups are generated), 5. define a next threshold , 6. repeat steps 4 and 5. The reduction of parameters depends on the number of defined thresholds. Each threshold defines two different groups of parameters. Principal Component Analysis Principal component analysis (PCA) is a multivariate statistical technique. It is based on a vector space transform that can be used to reduce data dimensionality by a covariance analysis between factors. The method generates a new set of parameters which are called principal components. Those new parameters are linear combinations of the original parameters. The PCA procedure for porosity and permeability estimation is as follows (Yadav 2006): 1. obtain a set of multiple porosity and permeability realizations for the reservoir by using geostatistical prior information (for example, by sequential Gaussian simulation), 2. estimate the covariance matrix of first for porosity,
71st EAGE Conference & Exhibition — Amsterdam, The Netherlands, 8 - 11 June 2009
3. calculate the eigenvectors and eigenvalues of these covariance matrix, 4. select the eigenvectors corresponding to the highest eigenvalues, 5. obtain a set of permeability realizations from the porosity realizations and empirical rock type knowledge (relation between porosity and permeability), 6. repeat steps 2 – 4 for the permeability realizations. Each eigenvector represents a parameter distribution pattern. These patterns are different to the available realizations. Only patterns associated to the highest eigenvalues are retained for analysis. Numerical results The approaches proposed are tested on a semi-synthetic reservoir from the Norne Field offshore Norway. The model is a two-phase (oil and water), two-dimensional black-oil reservoir which is discretized by 1014 (39×1×26) grid blocks with 739 active cells. It has one producer and one water injector which are located at column 6 and 33, respectively. The reservoir is subdivided into four different formations, Garn, Ile, Tofte, and Tilje, from top to base. Hydrocarbons are located in the Lower- to Middle-Jurassic sandstones. Different geological and dispositional environments gave rise to nine different rock types (with different properties). Since two parameters (porosity and horizontal permeability; vertical permeability is assumed to be equal to horizontal permeability) are attached to each grid cell, there are 1478 parameters to estimate. The objective function in this project represents the mismatch in seismic and production and can be written as follows J(θ) = w p M p (θ ) − M p θ * + w s M s (θ ) − M s θ* ,
( )
( )
(1)
where, J is the objective function. ws and wp are weighting factors for seismic and production mismatch, respectively, θ is the vector of unknown parameters (porosity and permeability) and Mp(θ), Mp(θ*), Ms(θ*) and Ms(θ*) represent simulated and observed production and seismic data respectively. The seismic data considered in this project are zero offset amplitudes and amplitude versus offset (AVO) gradients. The production data used are water and oil well production rates (WWPR and WOPR, respectively) and bottom hole pressures for injector and producer (WBHPI and WBHPP, respectively). Because this optimization problem can be formulated in terms of a least-squares minimization, the Gauss-Newton method can be applied to solve it. This optimization scheme requires gradient information. The fact that the commercial flow simulation software used in this project does not provide derivative information for the parameter estimation problem forces us to compute numerical and approximate derivatives. Though this approach in general might be very time consuming, when combined with parameter reduction technique and distributed computing, the whole procedure can be fairly efficient. The adjustment of the perturbation parameter in the gradient approximation may be problematic in some cases. In this work a selection of relative perturbation size of 1% for this parameter yields acceptable results. It should be noticed that (unlike when exact derivatives are used) the adjustment of the perturbation parameter in the gradient estimation can be beneficial in situations where noise is present in the data. Gradzones are defined based on parameter sensitivities with respect to the 4D seismic data. For this purpose the Hessian of the objective function (1) with wp = ws = 1 is computed. The eigenvectors of this matrix are ordered according to the eigenvalues, and as indicated above, gradzones can be defined by means of some thresholds. In this way, two sets of parameters have been created (49 and 92 parameters for porosity and permeability, respectively). In the first step of the PCA analysis porosity realizations are simulated. One thousand porosity realizations conditioned to local well data and variograms are obtained by sequential Gaussian simulation. Another one thousand permeability realizations are generated by empirical knowledge for the different rock types (relations between porosity and permeability). The number of principal components retained is such that 95% of the total variance is kept. This results in reducing the number of component to 100 for porosity and permeability, respectively.
71st EAGE Conference & Exhibition — Amsterdam, The Netherlands, 8 - 11 June 2009
An IBM P690, 32 CPU (Power 4X) with 32GB RAM operating on AIX (version 5.3) is used in this work. Each cost function evaluation (production and seismic simulations) takes 25 seconds. This means that the whole parameter estimation process, without considering any reparameterization technique and with only one processor, requires approximately 206 hours. If 20 processors are used in a distributed computing environment for simulation and gradient calculation (this approach will be denoted here by WREPA), the total time is reduced to 20 hours. If GZA and PCA reparameterization techniques are used (these strategies are denoted REGZA and REPCA) the time drops to 2 and 3 hours, respectively. The REPCA approach is applied in two different ways, REPCA-STG1 and REPCA-STG2. In REPCA-STG1 the principal component analysis are updated simultaneously while in REPCA-STG2 they are modified sequentially (i.e, first the component corresponding to the highest eigenvalue, then the second one, the third one and so on). In Figure 1 the performance of the four strategies above is compared (with respect to the cost function reduction). We conclude that all the approaches succeed in obtaining a cost function value acceptable for estimation purposes.
Table 1: Different reparameterization strategies
WREPA REGZA REPCA STG1
STG2
Figure 1: Objective function versus number of required forward simulation for the parameter estimation approaches considered.
Without REPArameterization REparameterization by GradZone Analysis REparameterization by Principal Component Analysis STrateGy 1 REparameterization by Principal Component Analysis STrateGy 2
Figures 2 and 3 represent real and estimated porosity and permeability distributions for all cases. The distributions corresponding to WREPA and REGZA are not acceptable from a geological point of view since the values for porosity and permeability are blocky and spatially disconnected. Therefore, these methods might not be able to detect (high permeable) thin layers. The low resolution of 4D seismic in vertical direction and the fact that WREPA and REGZA do not include complex geological constraints (unlike REPCA) can be two of the reasons of the malperformance mentioned above. The constraints considered in WREPA and REGZA are only bound constraints, clearly not enough to limit the solution space properly (from a geological point of view). By PCA the optimization is not only speeded up but also better conditioned thanks to additional continuity constraints. These constraints are implicitly included in the porosity variogram, in the well data (local conditioning), and in the relations between permeability and porosity. The solutions obtained using PCA are geologically acceptable.
71st EAGE Conference & Exhibition — Amsterdam, The Netherlands, 8 - 11 June 2009
Figure 2: Real and estimated porosity distributions for all the approaches considered.
Figure 3: Real and estimated permeability distributions for all the approaches considered
Conclusions
Gradzone and principal component analysis are incorporated in a distributed computing environment to speed up the parameter estimation. This process is formulated as a least-squares minimization problem and the Gauss-Newton method is used for that purpose. The parameters to estimate are porosity and permeability distributions, and the optimization cost function is a mismatch measurement for near and far offset time-lapse seismic and production data. If distributed computing and parameter reduction techniques are not used, the computational load in the complete process can be prohibitive for practical applications. Additionally, if some geological constraints are not included in the parameter estimation, the solutions obtained might not be geologically realistic. By principal component analysis we can take into account these geologic constraints of spatial continuity, and at the same time and in conjunction with a distributed computing framework, significantly speed up the optimization procedure. Since the seismic data considered do not have much vertical resolution and the approach introduced here tends to produce vertical spread, some very thin (high) permeable layers cannot be detected accurately. All in all, the good results obtained in a semi synthetic field case with the approach presented in this work encourage application of the methodology to parameter estimation of the real field
Acknowledgements
The authors thank the Norwegian University of Science and Technology (NTNU), the Integrated Operation Centre (CIO), the Research Council of Norway (NFR), and TOTAL for financial support. Authors are grateful to StatoilHydro for permitting the use of the reservoir model. We also thank Alexey Stovas for preparing the seismic forward model, and Jan Ivar Jensen, Knut Backe and Erlend Våtevik, all from NTNU, for their help during this study.
References
Arenas E, Kruijsdijk C, Oldenziel T 2001 Semi-Automatic History Matching Using the Pilot Point Method Including TimeLapse Seismic Data, SPE 71634, SPE Annual Technical Conference and Exhibition, New Orleans, USA, September 30 - October 3. Brun B, Gosselin O, and Barker J W 2001 Use of Prior Information in Gradient-Based History-Matching, SPE 66353, SPE Reservoir simulation symposium, Texas, USA, February 11–14.. Dong Y, and Oliver D S 2008 Reservoir Simulation Model Updates via Automatic History Matching With Integration of Seismic Impedance Change and Production Data, SPE 12550, International Petroleum Technology Conference, Kuala Lumpur, Malaysia, December 3-5. Gosselin O and Berg S V D 2001 Integrated History-Matching of Production and 4D Seismic Data, SPE 71599, SPE Annual Technical Conference and Exhibition, Louisiana, September 30 - October 3
71st EAGE Conference & Exhibition — Amsterdam, The Netherlands, 8 - 11 June 2009
Harb R 2004 History matching including 4-D seismic - An application to a field in the North Sea, Master thesis, NTNU Norway Marsily G, Levedan G, Boucher M, Fasanino G 1984 Interpretation of interference tests in a well field using geostatistical techniques to fit the permeability distribution in a reservoir model, Geostatistics for Natural Resources Characterization, 2 831. Nomura K, Kalia R K, Nakamo A, Vashishta P, and Landa J 2007 Parallel history matching and associated forcast at the centre for interactive smart oilfield technologies, J Supercomput, 41 109-117. Ouenes, Ahmed, Weiss, William 1995 Parallel Reservoir Automatic History Matching Using a Network of Workstations and PVM, SPE 29107, SPE Reservoir Simulation Symposium, San Antonio, Texas, USA, February 12-15. Reynolds A C, He N, Chu L, Oliver D S 1996 Reparameterization techniques for generating reservoir descriptions conditioned to variograms and well-test pressure data, Soc. Petrol. Eng. J., 1, 413 - 426. Ruijian L, Reynolds A C, Oliver D S 2003 History Matching of Three-Phase Flow Production Data, SPE 87336, SPE Journal, 4 328-340. Sarma P, Durlofsky L, Aziz K 2007 A new approach to automatic history matching using Kernel PCA, SPE 106176, SPE Reservoir Simulation Symposium, Houston, USA, February 26-28. Shah P C, Gavalas G R, and Seinfeld J H 1978 Error analysis in history matching: the optimum level of parameterization, Soc. Petrol. Eng. J., 18(6), 219–228. Yadav S 2006 History matching Using Face-Recognition Technique Based on Principal Component Analysis, SPE 102148, SPE Annual Technical Conference and Exhibition, Texas, September 24 - 27.
71st EAGE Conference & Exhibition — Amsterdam, The Netherlands, 8 - 11 June 2009