Kriging in a global neighborhood

The kriging estimator is usually computed in a moving neighborhood; only the data near the point to be estimated are used. This moving neighborhood approach creates discontinuities in mapping applications. An alternative approach is presented here, whereby all points are estimated using all the available data. To solve the resulting large linear system the kriging estimator is expressed in terms of the inverse of the covariance matrix. The covariance matrix has the advantage of being positive definite and the size of system which can be solved without encountering numerical instability is substantially increased. Because the kriging matrix does not change, the estimator can be written in terms of scalar products, thus avoiding the more time-consuming matrix multiplications of the standard approach. In the particular case of a covariance which is zero for distances greater than a fixed value (the range), the resulting banded structure of the covariance matrix is shown to lead to substantial computational savings in both run time and storage space. In this case the calculation time for the kriging variance is also substantially reduced. The present method is extended to the nonstationary case.     

Correcting the Smoothing Effect of Estimators: A Spectral Postprocessor

The postprocessing algorithm introduced by Yao for imposing the spectral amplitudes of a target covariance model is shown to be efficient in correcting the smoothing effect of estimation maps, whether obtained by kriging or any other interpolation technique. As opposed to stochastic simulation, Yao's algorithm yields a unique map starting from an original, typically smooth, estimation map. Most importantly it is shown that reproduction of a covariance/semivariogram model (global accuracy) is necessarily obtained at the cost of local accuracy reduction and increase in conditional bias. When working on one location at a time, kriging remains the most accurate (in the least squared error sense) estimator. However, kriging estimates should only be listed, not mapped, since they do not reflect the correct (target) spatial autocorrelation. This mismatch in spatial autocorrelation can be corrected via stochastic simulation, or can be imposed a posteriori via Yao's algorithm.     

Application of FFT-based Algorithms for Large-Scale Universal Kriging Problems

Looking at kriging problems with huge numbers of estimation points and measurements, computational power and storage capacities often pose heavy limitations to the maximum manageable problem size. In the past, a list of FFT-based algorithms for matrix operations have been developed. They allow extremely fast convolution, superposition and inversion of covariance matrices under certain conditions. If adequately used in kriging problems, these algorithms lead to drastic speedup and reductions in storage requirements without changing the kriging estimator. However, they require second-order stationary covariance functions, estimation on regular grids, and the measurements must also form a regular grid. In this study, we show how to alleviate these rather heavy and many times unrealistic restrictions. Stationarity can be generalized to intrinsicity and beyond, if decomposing kriging problems into the sum of a stationary problem and a formally decoupled regression task. We use universal kriging, because it covers arbitrary forms of unknown drift and all cases of generalized covariance functions. Even more general, we use an extension to uncertain rather than unknown drift coefficients. The sampling locations may now be irregular, but must form a subset of the estimation grid. Finally, we present asymptotically exact but fast approximations to the estimation variance and point out application to conditional simulation, cokriging and sequential kriging. The drastic gain in computational and storage efficiency is demonstrated in test cases. Especially high-resolution and data-rich fields such as rainfall interpolation from radar measurements or seismic or other geophysical inversion can benefit from these improvements.     

Kriging in the Presence of Locally Varying Anisotropy Using Non-Euclidean Distances

A stationary specification of anisotropy does not always capture the complexities of a geologic site. In this situation, the anisotropy can be varied locally. Directions of continuity and the range of the variogram can change depending on location within the domain being modeled. Kriging equations have been developed to use a local anisotropy specification within kriging neighborhoods; however, this approach does not account for variation in anisotropy within the kriging neighborhood. This paper presents an algorithm to determine the optimum path between points that results in the highest covariance in the presence of locally varying anisotropy. Using optimum paths increases covariance, results in lower estimation variance and leads to results that reflect important curvilinear structures. Although CPU intensive, the complex curvilinear structures of the kriged maps are important for process evaluation. Examples highlight the ability of this methodology to reproduce complex features that could not be generated with traditional kriging.     

Mononodal indicator variography—Part I: Theory

For any distribution of grades, a particular cutoff grade is shown here to exist at which the indicator covariance is proportional to the grade covariance to a very high degree of accuracy. The name “mononodal cutoff” is chosen to denote this grade. Its importance for robust grade variography in the presence of a large coefficient of variation—typical of precious metals—derives from the fact that the mononodal indicator variogram is then linearly related to the grade variogram yet is immune to outlier data and is found to be particularly robust under data information reduction. Thus, it is an excellent substitute to model in lieu of a difficult grade variogram. A theoretical expression for the indicator covariance is given as a double series of orthogonal polynomials that have the grade density function as weight function. Leading terms of this series suggest that indicator and grade covariances are first-order proportional, with cutoff grade dependence being carried by the proportionality factor. Kriging equations associated with this indicator covariance lead to cutoff-free kriging weights that are identical to grade kriging weights. This circumstance simplifies indicator kriging used to estimate local point-grade histograms, while at the same time obviating order relations problems.     

Inner product matrices, kriging, and nonparametric estimation of variogram

Two important problems in the practical implementation of kriging are: (1) estimation of the variogram, and (2) estimation of the prediction error. In this paper, a nonparametric estimator of the variogram to circumvent the problem of the precise choice of a variogram model is proposed. Using orthogonal decomposition of the kriging predictor and the prediction error, a method for selecting, what may be considered, a statistical neighborhood is suggested. The prediction error estimates based on this scheme, in fact, reflects the true prediction error, thus leading to proper coverage for the corresponding prediction interval. By simulations and a reanalysis of published data, it is shown that the proposals made in this paper are useful in practice.     

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.     

Conditional Spectral Simulation with Phase Identification

Spectral simulation is used widely in electrical engineering to generate random fields with a given covariance spectrum. The algorithms used are fast particularly when based on Fast Fourier Transform (FFT). However, because of lack of phase identification, spectral simulation only generates unconditional realizations. Local data conditioning is obtained typically by adding a simulated kriging residual. This conditioning process requires an additional kriging at each simulated node thus forfeiting the speed advantage of FFT. A new algorithm for conditioning is proposed whereby the phase values are determined iteratively to ensure approximative data reproduction while reproducing the frequency spectrum, that is, the covariance model. A case study is presented to demonstrate the algorithm.     

Kriging without negative weights

Under a constant drift, the linear kriging estimator is considered as a weighted average ofn available sample values. Kriging weights are determined such that the estimator is unbiased and optimal. To meet these requirements, negative kriging weights are sometimes found. Use of negative weights can produce negative block grades, which makes no practical sense. In some applications, all kriging weights may be required to be nonnegative. In this paper, a derivation of a set of nonlinear equations with the nonnegative constraint is presented. A numerical algorithm also is developed for the solution of the new set of kriging equations.     

A sensitivity analysis for universal kriging

Universal kriging is an interpolation method for producing contour maps from irregularly spaced sample points, taking into account the trend (or polynomial drift), which is of known form. It assumes a known covariance model to express correlation of points a short distance apart. A sensitivity analysis examines how the fitted surface will change for given perturbations in the covariance model. We develop a simple approximate analysis in preference to exact analysis and show that it is adequate for small perturbations. For large data sets, a dramatic reduction in computer time is possible using approximate analysis. Possible extensions of this work are noted.     

Improvement of Fourier-Based Unconditional and Conditional Simulations for Band Limited Fractal (von Kármán) Statistical Models

We evaluate the performance and statistical accuracy of the fast Fourier transform method for unconditional and conditional simulation. The method is applied under difficult but realistic circumstances of a large field (1001 by 1001 points) with abundant conditioning criteria and a band limited, anisotropic, fractal-based statistical characterization (the von Kármán model). The simple Fourier unconditional simulation is conducted by Fourier transform of the amplitude spectrum model, sampled on a discrete grid, multiplied by a random phase spectrum. Although computationally efficient, this method failed to adequately match the intended statistical model at small scales because of sinc-function convolution. Attempts to alleviate this problem through the covariance method (computing the amplitude spectrum by taking the square root of the discrete Fourier transform of the covariance function) created artifacts and spurious high wavenumber content. A modified Fourier method, consisting of pre-aliasing the wavenumber spectrum, satisfactorily remedies sinc smoothing. Conditional simulations using Fourier-based methods require several processing stages, including a smooth interpolation of the differential between conditioning data and an unconditional simulation. Although kriging is the ideal method for this step, it can take prohibitively long where the number of conditions is large. Here we develop a fast, approximate kriging methodology, consisting of coarse kriging followed by faster methods of interpolation. Though less accurate than full kriging, this fast kriging does not produce visually evident artifacts or adversely affect the a posteriori statistics of the Fourier conditional simulation.     

Is Lognormal Kriging Suitable for Local Estimation?

Lognormal kriging was developed early in geostatistics to take account of the often seen skewed distribution of the experimental mining data. Intuitively, taking the distribution of the data into account should lead to a better local estimate than that which would have been obtained when it is ignored. In practice however, the results obtained are sometimes disappointing. This paper tries to explain why this is so from the behavior of the lognormal kriging estimator. The estimator is shown to respect certain unbiasedness properties when considering the whole working field using the regression curve and its confidence interval for both simple or ordinary kriging. When examined locally, however, the estimator presents a behavior that is neither expected nor intuitive. These results lead to the question: is the theoretically correct lognormal kriging estimator suited to the practical problem of local estimation?     

Robust kriging—A proposal

Geological data frequently have a heavy-tailed normal-in-the-middle distribution, which gives rise to grade distributions that appear to be normal except for the occurrence of a few outliers. This same situation also applies to log-transformed data to which lognormal kriging is to be applied. For such data, linear kriging is nonrobust in that (1)kriged estimates tend to infinity as the outliers do, and (2)it is also not minimum mean squared error. The more general nonlinear method of disjunctive kriging is even more nonrobust, computationally more laborious, and in the end need not produce better practical answers. We propose a robust kriging method for such nearly normal data based on linear kriging of an editing of the data. It is little more laborious than conventional linear kriging and, used in conjunction with a robust estimator of the variogram, provides good protection against the effects of data outliers. The method is also applicable to time series analysis.     

Bayesian Maximum Entropy Analysis and Mapping: A Farewell to Kriging Estimators?

The Bayesian Maximum Entropy (BME) method of spatial analysis and mapping provides definite rules for incorporating prior information, hard and soft data into the mapping process. It has certain unique features that make it a loyal guardian of plausible reasoning under conditions of uncertainty. BME is a general approach that does not make any assumptions regarding the linearity of the estimator, the normality of the underlying probability laws, or the homogeneity of the spatial distribution. By capitalizing on various sources of information and data, BME introduces an epistemological framework that produces predictive maps that are more accurate and in many cases computationally more efficient than those derived by traditional techniques. In fact, kriging techniques can be derived as special cases of the BME approach, under restrictive assumptions regarding the prior information and the data available. BME is a more rigorous approach than indicator kriging for incorporating soft data. The BME formulation, in fact, applies in a spatial or a spatiotemporal domain and its extension to the case of block and vector random fields is straightforward. New theoretical results are presented and numerical examples are discussed, which use the BME approach to account for important sources of knowledge in a systematic manner. BME can be useful in practical situations in which prior information can be used to compensate for the limited amount of measurements available (e.g., preliminary or feasibility study levels) or soft data are available that can be combined with hard data to improve mapping significantly. BME may be then viewed as an effort towards the development of a more general framework of spatial/temporal analysis and mapping, which includes traditional geostatistics as its limiting case, and it also provides the means to derive novel results that could not be obtained by traditional geostatistics.     

On the influence of kriging parameters on the cartographic output—A study in mapping subglacial topography

Geostatistics provides a suite of methods, summarized as kriging, to analyze a finite data set to describe a continuous property of the Earth. Kriging methods consist of moving window optimum estimation techniques, which are based on a least-squares principle and use a spatial structure function, usually the variogram. Applications of kriging techniques have become increasingly wide-spread, with ordinary kriging and universal kriging being the most popular ones. The dependence of the final map or model on the input, however, is not generally understood. Herein we demonstrate how changes in the kriging parameters and the neighborhood search affect the cartographic result. Principles are illustrated through a glaciological study. The objective is to map ice thickness and subglacial topography of Storglaciären, Kebnekaise Massif, northern Sweden, from several sets of radio-echo soundings and hot water drillings. New maps are presented.     

Kriging Regionalized Positive Variables Revisited: Sample Space and Scale Considerations

Frequently, regionalized positive variables are treated by preliminarily applying a logarithm, and kriging estimates are back-transformed using classical formulae for the expectation of a lognormal random variable. This practice has several problems (lack of robustness, non-optimal confidence intervals, etc.), particularly when estimating block averages. Therefore, many practitioners take exponentials of the kriging estimates, although the final estimations are deemed as non-optimal. Another approach arises when the nature of the sample space and the scale of the data are considered. Since these concepts can be suitably captured by an Euclidean space structure, we may define an optimal kriging estimator for positive variables, with all properties analogous to those of linear geostatistical techniques, even for the estimation of block averages. In this particular case, no assumption on preservation of lognormality is needed. From a practical point of view, the proposed method coincides with the median estimator and offers theoretical ground to this extended practice. Thus, existing software and routines remain fully applicable.     

An application of the deterministic variogram to design-based variance estimation

When estimating the mean value of a variable, or the total amount of a resource, within a specified region it is desirable to report an estimated standard error for the resulting estimate. If the sample sites are selected according to a probability sampling design, it usually is possible to construct an appropriate design-based standard error estimate. One exception is systematic sampling for which no such standard error estimator exists. However, a slight modification of systematic sampling, termed 2-step tessellation stratified (2TS) sampling, does permit the estimation of design-based standard errors. This paper develops a design-based standard error estimator for 2TS sampling. It is shown that the Taylor series approximation to the variance of the sample mean under 2TS sampling may be expressed in terms of either a deterministic variogram or a deterministic covariance function. Variance estimation then can be approached through the estimation of a variogram or a covariance function. The resulting standard error estimators are compared to some more traditional variance estimators through a simulation study. The simulation results show that estimators based on the new approach may perform better than traditional variance estimators.