Multivariable variogram and its application to the linear model of coregionalization

In this article, we present the multivariable variogram, which is defined in a way similar to that of the traditional variogram, by the expected value of a distance, squared, in a space withp dimensions. Combined with the linear model of coregionalization, this tool provides a way for finding the elementary variograms that characterize the different spatial scales contained in a set of data withp variables. In the case in which the number of elementary components is less than or equal to the number of variables, it is possible, by means of nonlinear regression of variograms and cross-variograms, to estimate the coregionalization parameters directly in order to obtain the elementary variables themselves, either by cokriging or by direct matrix inversion. This new tool greatly simplifies the procedure proposed by Matheron (1982) and Wackernagel (1985). The search for the elementary variograms is carried out using only one variogram (multivariable), as opposed to thep(p + 1)/2 required by the Matheron approach. Direct estimation of the linear coregionalization model parameters involves the creation of semipositive definite coregionalization matrices of rank 1.     

Conditional simulation of intrinsic random functions of order k

Conditional simulation of intrinsic random functions of orderk is a stochastic method that generates realizations which mimic the spatial fluctuation of nonstationary phenomena, reproduce their generalized covariance and honor the available data at sampled locations. The technique proposed here requires the following steps: (i) on-line simulation of Wiener-Levy processes and of their integrations; (ii) use of the turning-bands method to generate realizations in Rⁿ; (iii) conditioning to available data; and (iv) verification of the reproduced generalized covariance using generalized variograms. The applicational aspects of the technique are demonstrated in two and three dimensions. Examples include the conditional simulation of geological variates of the Crystal Viking petroleum reservoir, Alberta, Canada.     

Automatic Modeling of (Cross) Covariance Tables Using Fast Fourier Transform

Covariance models provide the basic measure of spatial continuity in geostatistics. Traditionally, a closed-form analytical model is fitted to allow for interpolation of sample Covariance values while ensuring the positive definiteness condition. For cokriging, the modeling task is made even more difficult because of the restriction imposed by the linear coregionalization model. Bochner's theorem maps the positive definite constraints into much simpler constraints on the Fourier transform of the covariance, that is the density spectrum. Accordingly, we propose to transform the experimental (cross) covariance tables into quasidensity spectrum tables using Fast Fourier Transform (FFT). These quasidensity spectrum tables are then smoothed under constraints of positivity and unit sum. A backtransform (FFT) yields permissible (jointly) positive definite (cross) covariance tables. At no point is any analytical modeling called for and the algorithm is not restricted by the linear coregionalization model. A case study shows the proposed covariance modeling to be easier and much faster than the traditional analytical covariance modeling, yet yields comparable kriging or simulation results.     

On a controversial method for modeling a coregionalization

This paper reviews two alternative approaches for modeling the (cross) variograms in a coregionalization: (1) fitting the traditional linear model of coregionalization. or (2) deducing the crossvariogram model as a linear combination of prior direct (auto) variogram models while checking the Cauchy-Schwarz inequalities. We show that the second approach has no practical advantage over the traditional one and may not be valid if more than two variables are involved. In such case. Cauchy-Schwarz inequalities are necessary but not sufficient conditions for validity of a coregionalization model.     

Tests of Significance for Structural Correlations in the Linear Model of Coregionalization

In the linear model of coregionalization (LMC), when applicable to the experimental direct variograms and the experimental cross variogram computed for two random functions, the variability of and relationships between the random functions are modeled with the same basis functions. In particular, structural correlations can be defined from entries of sill matrices (coregionalization matrices) under second-order stationarity. In this article, modified t-tests are proposed for assessing the statistical significance of estimated structural correlations. Their specific aspects and fundamental differences, compared with an existing modified t-test for global correlation analysis with spatial data, are discussed via estimated effective sample sizes, in relation to the superimposition of random structural components, the range of autocorrelation, the presence of correlation at another structure, and the sampling scheme. Accordingly, simulation results are presented for one structure versus two structures (one without and the other with autocorrelation). The performance of tests is shown to be related to the uncertainty associated with the estimation of variogram model parameters (range, sill matrix entries), because these are involved in the test statistic and the degrees of freedom of the associated t-distribution through the estimated effective sample size. Under the second-order stationarity and LMC assumptions, the proposed tests are generally valid.     

Universal cokriging under intrinsic coregionalization

Under the intrinsic coregionalization model if both primary and secondary measurements are available at all sample locations, the conventional geostatistical wisdom is that cokriging provides exactly the same solution as univariate kriging on the primary process alone. However, recent eamples have been given where nonzero secondary cokriging weights have accurred under this spatial dependence structure. This note identifies the conditions under which secondary information is useful under the assumption of intrinsic coregionalization. An illustration is given using a dataset of plutonium and americium concentrations collected from a region of the Nevada Test Site.     

A Fixed-Path Markov Chain Algorithm for Conditional Simulation of Discrete Spatial Variables

The Markov chain random field (MCRF) theory provided the theoretical foundation for a nonlinear Markov chain geostatistics. In a MCRF, the single Markov chain is also called a “spatial Markov chain” (SMC). This paper introduces an efficient fixed-path SMC algorithm for conditional simulation of discrete spatial variables (i.e., multinomial classes) on point samples with incorporation of interclass dependencies. The algorithm considers four nearest known neighbors in orthogonal directions. Transiograms are estimated from samples and are model-fitted to provide parameter input to the simulation algorithm. Results from a simulation example show that this efficient method can effectively capture the spatial patterns of the target variable and fairly generate all classes. Because of the incorporation of interclass dependencies in the simulation algorithm, simulated realizations are relatively imitative of each other in patterns. Large-scale patterns are well produced in realizations. Spatial uncertainty is visualized as occurrence probability maps, and transition zones between classes are demonstrated by maximum occurrence probability maps. Transiogram analysis shows that the algorithm can reproduce the spatial structure of multinomial classes described by transiograms with some ergodic fluctuations. A special characteristic of the method is that when simulation is conditioned on a number of sample points, simulated transiograms have the tendency to follow the experimental ones, which implies that conditioning sample data play a crucial role in determining spatial patterns of multinomial classes. The efficient algorithm may provide a powerful tool for large-scale structure simulation and spatial uncertainty analysis of discrete spatial variables.     

Characterization and Quantification of Uncertainty in?Coregionalization Analysis

Coregionalization analysis has been presented as a method of multi-scale analysis for multivariate spatial data. Despite an increasing use of this method in environmental and earth sciences, the uncertainty associated with the estimation of parameters in coregionalization analysis (e.g., sills and functions of sills) is potentially high and has not yet been characterized. This article aims to discuss the theory underlying coregionalization analysis and assess the robustness and limits of the method. A theoretical framework is developed to calculate the ergodic and fluctuation variance-covariance matrices of least-squares estimators of sills in the linear model of coregionalization. To adjust for the positive semidefiniteness constraint on estimated coregionalization matrices, a confidence interval estimation procedure for sills and functions of sills is presented. Thereafter, the relative importance of uncertainty measures (bias and variance) for sills and structural coefficients of correlation and determination is assessed under different scenarios to identify factors controlling their uncertainty. Our results show that the sampling grid density, the choice of the least-squares estimator of sills, the positive semidefiniteness constraint, the presence of scale dependence in the correlations, and the number and range of variogram models, all affect the level of uncertainty, sometimes through multiple interactions. The asymptotic properties of variogram model parameter estimators in a bounded sampling domain impose a theoretical limit to their accuracy and precision. Because of this limit, the uncertainty was found to be high for several scenarios, especially with three variogram models, and was often more dependent on the ratio of variogram range to domain extent than on the sampling grid density. In practice, in the coregionalization analysis of a real dataset, the circular requirement for sill estimates in the calculation of uncertainty measures makes the quantification of uncertainty very problematic, if not impossible. The use of coregionalization analysis must be made with due knowledge of the uncertainty levels and limits of the method.     

Combination of Dependent Realizations Within the Gradual Deformation Method

Gradual deformation is a parameterization method that reduces considerably the unknown parameter space of stochastic models. This method can be used in an iterative optimization procedure for constraining stochastic simulations to data that are complex, nonanalytical functions of the simulated variables. This method is based on the fact that linear combinations of multi-Gaussian random functions remain multi-Gaussian random functions. During the past few years, we developed the gradual deformation method by combining independent realizations. This paper investigates another alternative: the combination of dependent realizations. One of our motivations for combining dependent realizations was to improve the numerical stability of the gradual deformation method. Because of limitations both in the size of simulation grids and in the precision of simulation algorithms, numerical realizations of a stochastic model are never perfectly independent. It was shown that the accumulation of very small dependence between realizations might result in significant structural drift from the initial stochastic model. From the combination of random functions whose covariance and cross-covariance are proportional to each other, we derived a new formulation of the gradual deformation method that can explicitly take into account the numerical dependence between realizations. This new formulation allows us to reduce the structural deterioration during the iterative optimization. The problem of combining dependent realizations also arises when deforming conditional realizations of a stochastic model. As opposed to the combination of independent realizations, combining conditional realizations avoids the additional conditioning step during the optimization process. However, this procedure is limited to global deformations with fixed structural parameters.     

Deriving Constraints on Small-Scale Variograms due to Variograms of Large-Scale Data

The application of kriging-based geostatistical algorithms to integrate large-scale seismic data calls for direct and cross variograms of the seismic variable and primary variable (e.g., porosity) at the modeling scale, which is typically much smaller than the seismic data resolution. In order to ensure positive definiteness of the cokriging matrix, a licit small-scale coregionalization model has to be built. Since there are no small-scale secondary data, an analytical method is presented to infer small-scale seismic variograms. The method is applied to estimate the 3-D porosity distribution of a West Texas oil field given seismic data and porosity data at 62 wells.     

Characterization and Quantification of Uncertainty in Coregionalization Analysis

Coregionalization analysis has been presented as a method of multi-scale analysis for multivariate spatial data. Despite an increasing use of this method in environmental and earth sciences, the uncertainty associated with the estimation of parameters in coregionalization analysis (e.g., sills and functions of sills) is potentially high and has not yet been characterized. This article aims to discuss the theory underlying coregionalization analysis and assess the robustness and limits of the method. A theoretical framework is developed to calculate the ergodic and fluctuation variance-covariance matrices of least-squares estimators of sills in the linear model of coregionalization. To adjust for the positive semidefiniteness constraint on estimated coregionalization matrices, a confidence interval estimation procedure for sills and functions of sills is presented. Thereafter, the relative importance of uncertainty measures (bias and variance) for sills and structural coefficients of correlation and determination is assessed under different scenarios to identify factors controlling their uncertainty. Our results show that the sampling grid density, the choice of the least-squares estimator of sills, the positive semidefiniteness constraint, the presence of scale dependence in the correlations, and the number and range of variogram models, all affect the level of uncertainty, sometimes through multiple interactions. The asymptotic properties of variogram model parameter estimators in a bounded sampling domain impose a theoretical limit to their accuracy and precision. Because of this limit, the uncertainty was found to be high for several scenarios, especially with three variogram models, and was often more dependent on the ratio of variogram range to domain extent than on the sampling grid density. In practice, in the coregionalization analysis of a real dataset, the circular requirement for sill estimates in the calculation of uncertainty measures makes the quantification of uncertainty very problematic, if not impossible. The use of coregionalization analysis must be made with due knowledge of the uncertainty levels and limits of the method.     

On the Equivalence of the Cokriging and Kriging Systems

Simple cokriging of components of a p-dimensional second-order stationary random process is considered. Necessary and sufficient conditions under which simple cokriging is equivalent to simple kriging are given. Essentially this condition requires that it should be possible to express the cross-covariance at any lag series h using the cross-covariance at |h|=0 and the auto-covariance at lag series h. The mosaic model, multicolocated kriging and the linear model of coregionalization are examined in this context. A data analytic method to examine whether simple kriging of components of a multivariate random process is equivalent to its cokriging is given     

Conditional Indicator Simulation: Application to a Saskatchewan uranium deposit

Mineral deposits frequently exhibit a mixture of rock types in which each type can be identified by a characteristic metal concentration. Such a mixture can be correctly simulated by first reproducing the spatial and geometric configuration of the various rock types in the deposit. Then the grades for each rock type can be jointly simulated and filled in according to their specific coregionalization characteristics. The method of Conditional Indicator Simulation and an uranium—arsenic joint simulation are presented with a detailed, step-by-step application to the Midwest deposit, a high grade uranium deposit in northern Saskatchewan.     

Conditional simulation of intrinsic random functions of orderk

Conditional simulation of intrinsic random functions of orderk is a stochastic method that generates realizations which mimic the spatial fluctuation of nonstationary phenomena, reproduce their generalized covariance and honor the available data at sampled locations. The technique proposed here requires the following steps: (i) on-line simulation of Wiener-Levy processes and of their integrations; (ii) use of the turning-bands method to generate realizations in Rⁿ; (iii) conditioning to available data; and (iv) verification of the reproduced generalized covariance using generalized variograms. The applicational aspects of the technique are demonstrated in two and three dimensions. Examples include the conditional simulation of geological variates of the Crystal Viking petroleum reservoir, Alberta, Canada.     

T -distributed Random Fields: A Parametric Model for Heavy-tailedWell-log Data1

Histograms of observations from spatial phenomena are often found to be more heavy-tailed than Gaussian distributions, which makes the Gaussian random field model unsuited. A T-distributed random field model with heavy-tailed marginal probability density functions is defined. The model is a generalization of the familiar Student-T distribution, and it may be given a Bayesian interpretation. The increased variability appears cross-realizations, contrary to in-realizations, since all realizations are Gaussian-like with varying variance between realizations. The T-distributed random field model is analytically tractable and the conditional model is developed, which provides algorithms for conditional simulation and prediction, so-called T-kriging. The model compares favourably with most previously defined random field models. The Gaussian random field model appears as a special, limiting case of the T-distributed random field model. The model is particularly useful whenever multiple, sparsely sampled realizations of the random field are available, and is clearly favourable to the Gaussian model in this case. The properties of the T-distributed random field model is demonstrated on well log observations from the Gullfaks field in the North Sea. The predictions correspond to traditional kriging predictions, while the associated prediction variances are more representative, as they are layer specific and include uncertainty caused by using variance estimates.     

The Correlation Bias for Two-Dimensional Simulations by Turning Bands

The turning bands method (TBM) generates realizations of isotropic Gaussian random fields by summing contributions from line processes. We consider two-dimensional simulations and study the correlation bias attributable to the use of only a finite number L of lines. Our analytical and numerical results confirm that the maximal bias is of order 1/L, and that L = 64 lines suffice for excellent covariance reproduction. The notorious banding observed in simulations with an insufficient number of lines is a related but different phenomenon and depends strongly on the choice of the line simulation technique. Clear-cut recommendations for the number of lines necessary to avoid the effect can only be based on practical experience with the specific code at hand.     

Blending-Based Stochastic Simulator

This paper presents a new method of constructing random functions whose realizations can be evaluated efficiently. The basic idea is to blend, both stochastically and linearly, a limited set of independent initial realizations previously generated by any chosen simulation method. The blending stochastic coefficients are determined in such a way that the new random function so generated has the same mean and covariance functions as the random function used for generating the initial realizations.     

Modeling soil variability as a random field

The observed variability in the spatial distribution of soil properties suggests that it is natural to describe such distribution as a random field. One of the ways to study engineering problems in such a stochastic setting is by the use of the Monte-Carlo simulation procedure. Application of this technique requires the capability to generate a large number of realizations of a given random field. A numerical procedure for the generation of such realizations in two-dimensional space is introduced as a finite difference approximation of a stochastic differential equation. The equation used was that suggested by Heine (1955). The resulting procedure is essentially similar to other autoregressive procedures used for the same purpose (Whittle, 1954; Smith and Freeze, 1979). However, contrary to these procedures, the present one is defined in terms of physically significant parameters:r ₀, the autocorrelation distance;, the discretization size; and the standard deviation, . Formulating the simulation procedure in terms of the physically significant parameters (r ₀,, ) greatly simplifies the task of generating realizations that are compatable with a given soil deposit.     

Two Artifacts of Probability Field Simulation

Probability field simulation is being used increasingly to simulate geostatistical realizations. The method can be faster than conventional simulation algorithms and it is well suited to integrate prior soft information in the form of local probability distributions. The theoretical basis of probability field simulation has been established when there are no conditioning data; however, no such basis has been established in presence of conditioning data. Realizations generated by probability field simulation show two severe artifacts near conditioning data. We document these artifacts and show theoretically why they exist. The two artifacts that have been investigated are (1) local conditioning data appear as local minima or maxima of the simulated values, and (2) the variogram model in range of conditioning data is not honored; the simulated values have significantly greater continuity than they are supposed to. These two artifacts are predicted by theory. An example flow simulation study is presented to illustrate that they affect more than the visual appearance of the simulated realizations. Notwithstanding the flexibility of the probability field simulation method, these two artifacts suggest that it be used with caution in presence of conditioning data. Future research may overcome these limitations.     

Two-dimensional simulation by turning bands

Journel (1974) developed the turning-bands method which allows a three-dimensional data set with specified covariance to be obtained by the simulation of several one-dimensional realizations which have an intermediate covariance. The relationship between the threedimensional and one-dimensional covariance is straightforward and allows the one-dimensional covariance to be obtained immediately. In theory a dense uniform distribution of lines in three-dimensional space is required along which the one-dimensional realizations are generated; in practice most workers have been content to use the fifteen axes of the regular icosahedron. Many mining problems may be treated in two dimensions, and in this paper a turning-bands approach is developed to generate two-dimensional data sets with a specified covariance. By working in two dimensions, the area on which the data is simulated may be divided as finely as desired by the lines on which the one-dimensional realizations are first generated. The relationship between the two-dimensional and one-dimensional covariance is derived as a nontrivial integral equation. This is solved analytically for the onedimensional covariance. The method is applied to the generation of a two-dimensional data set with spherical covariance.