Direct sequential indicator simulation

Sequential Gaussian simulation (sgsim) and sequential indicator simulation (sisim) have emerged as powerful tools for stochastic imaging of Earth Science phenomena and are currently widely accepted fast simulation algorithms. In this paper, we will expand the sequential simulation toolbox by relying on a recent development in geostatistical theory (Journel, 1994) denoted as direct sequential simulation (dssim). This theorem states that in order to obtain variogram reproduction in the resulting simulations, any type of local conditional distribution can be used to simulate values, as long as its mean and variance identifies the kriging mean and kriging variance. This theorem is a clear extension of the Gaussian simulation paradigm where the kriging mean and variance are imposed on the standard normal distribution. dssim allows variogram reproduction under nonGaussian assumptions, however, it does not ensure reproduction of the global histogram of the data, neither does dssim allow, through its construction, the reproduction of indicator variograms. We propose a new direct sequential simulation method that solves these two shortcomings by imposing linear constraints on the successive local conditional distributions. We present our results on a case study. Introduction Stochastic simulation is a widely accepted tool in various areas of geostatistics. The goal of stochastic simulation is to reproduce geological texture in a set of equiprobable (=equally likely to be drawn) simulated realizations. Simulations are termed globally accuracte through the reproduction of one-, two-, or multiple-point statistics representative for the area under study. Its counterpart, kriging, is locally accurate in the minimum error variance sense, yet it provides inaccurate representations of spatial variability. Sequential simulation (Johnson, 1987), in its various geostatistical flavors (Isaaks, 1990, Gomez-Hernandez and Journel, 1993), has grown to be one of the most popular and computationally efficient tools for obtaining simulations that honor a certain histogram and a certain type of variogram(s) which have been obtained from the data. In mathematical terms, the most convenient method for simulation is sequential Gaussian simulation (sgsim, Deutsch and Journel, 1998) because all successive conditional distributions from which simulated values are drawn are Gaussian with parameters determined by the solution of a simple kriging system. However various limitations and shortcomings can be attributed to the sequential Gaussian simulation: sgsim relies on the assumption of multi-variate Gaussianity, an assumption that can never be fully checked in practice, yet always seems to be taken for granted. Multi-Gaussianity leads to simulated realizations that have maximally disconnected extremes (maximum entropy), a property that often conflicts with geological reality. sgsim requires a transformation into Gaussian space before simulation and a corresponding backtransformation after simulation is finished. However, often the primary variable to be simulated has to be conditioned to a secondary variable that is a linear or non-linear volume average of the primary variable (Journel, 1998). Normal-score transforms are non-linear transforms, hence they destroy the possible linear relation that exists between primary and secondary variable, or, they change the non-linearlity if that relation is non-linear. sgsim reproduces, by theory, only the normal score variogram, not the original variogram model. Usually reproduction of the normal score variogram entails reproduction of the original data variogram if the data histogram is not too skewed. However in case of high skewness, the reproduction of the variogram model after back-transformation is not guaranteed at all. It is therefore much more convenient to simulate directly in the space of the primary variable without transformation and without relying on the usual multi-Gaussian assumptions. Direct sequential simulation (dssim) (Journel, 1994, Caers, 2000) has been proposed for simulating directly in the original data space and does not rely on multi-Gaussian assumptions. In fact, dssim can produce realizations with low entropy characteristics (more connected extremes) than traditional sgsim. dssim relies on the important theoretical result (Journel, 1994) that, in order to reproduce a given covariance model, the successive conditional distributions used in the sequential path can be of any type as long as they identify the kriging mean and variance. There are however two important limitations to the current dssim algorithm: dssim does not ensure reproduction of the histogram, the only univariate statistics reproduced is the global sample mean and variance. Hence, a posterior identification of the histogram might be necessary, which might detroy the variogram reproduction. dssim does not allow through its construction the reproduction of indicator variograms as would be possible through the use of a sequential indicator simulation (Deutsch and Journel, 1998). Reproduction of indicator variograms is important when connectivity of extremes have been quantified from data and need to be reproduced in the simulated realizations. Note that by reproduction of a statistics (histogram/variogram), we mean reproduction on average, i.e. over a number of realizations. In this paper we propose a solution to these two shortcomings. Direct sequential simulation: a recall Sequential simulation The aim of sequential simulation, as it was originally constructed, is to reproduce desirable multivariate properties through the sequential use of conditional distributions. Consider a set of random variables defined at locations . The objective is to generate several joint realizations conditional to the available data and to some structural model such as the variogram. It can be shown that an -point multivariate distribution can be decomposed into a set of one-point conditional cdf’s as (1) where is the conditional cdf of given the set of original data values and the previous realizations . This decomposition allows us to generate an image by sequentially visiting each node. Sequential simulation, under a given multivariate distribution, amounts to read the decomposition (1) from left to right, i.e. the purpose of sequential simulation is to reproduce the properties of the given multivariate distribution. The simulation algorithm proceeds as follows, (see for example Goovaerts, 1998, p 380) Perform a transformation if neccesitated by the theory Define a random path visiting all nodes For each node , , do – model the conditional distribution of , given the original data values and all previously drawn values , – draw the simulated value from End loop Perform a backtransform to identify the target histogram (if needed). The last step ensures reproduction of the target univariate distribution deemed relevant for the field of observation. Under a multivariate Gaussian distribution, each conditional distribution is Gaussian. The purpose of sequential Gaussian simulation is to reproduce the properties of the multi-variate Gaussian distribution. In sequential Gaussian simulation the data has to be transformed into Gaussian space, hence sgsim reproduces only the normal score variogram, not the original space variogram However, relation (1) states that any product of a series of conditional distributions results in a specific multivariate distribution (not necessarily multi-variate Gaussian) that will be reflected in the resulting simulated realizations. Variogram reproduction is ensured if we constrain the local conditional distributions using the following theorem (Journel, 1994). Theorem 1 For the sequential simulation algorithm to reproduce a specific covariance model it suffices that all ccdf’s identify the simple kriging mean and variance derived from that covariance model The practical importance of this theorem for sequential simulation cannot be overstated. The conditional distribution can be of any type and need not be the same at each location. The distribution may have more than two parameters allowing identification of other moments or quantiles beyond the simple kriging mean and variance. In Caers (2000), it was shown that the choice of the type of local conditional distributions determines the multivariate properties of the resulting simulated realizations. Local distributions with long tails (exponential, power laws) tend to increase the connectivity of extreme values and result in highly non-Gaussian simulations while shorter tailed distributions such as the Gaussian or the uniform tend to enforce a high entropy property on the simulations. The need for back-transformation Direct sequential simulation only identifies the variogram. The histogram observed in the resulting realizations will typically depend on the type of local conditional distribution used, on the value of the global mean and global variance (as used in the simple kriging) and on the amount of conditioning data available. Theorem 1 does not ensure reproduction of the histogram of the data, hence a backtransformation should be applied to identify any target histogram. The GSLIB-program trans (Journel and Xu, 1996) can be used to perform this task. However it would be more convenient to have a method that reproduces histogram and variogram at the same time, without having to rely on a posterior identification of the histogram which might destroy the variogram reproduction. A discretized local conditional distribution model The model The dssim theory specifies that any type of local conditional distribution (ccdf) can be used. In order to obtain a flexible non-parametric model for the local conditional distribution, we propose to define the ccdf at each by a series of probabilities at predefined thresholds. These thresholds could be, for example, a series of quantile values. We define a set of thresholds for which there exist a series of conditional probability densities The following specifies how the are determined Theorem 1 enforces the following linear constraints on this series of local probabilities, for the kriging mean (2) and for the kriging variance we use the second moment (3) Also the following consistency conditions have to hold and (4) As a final constraint we force the local conditional distribution to be as close as possible (in least square sense) to the global histogram. This can be achieved by defining an objective function such that is minimal where are the global proportions obtained from the histogram, i.e If the derivative of the objective function is calculated with respect to the probabilities we get the evident result that (5) (A) 2-D Reference Data