 ### Geostatistics - cracks in the foundations ?

Matheronian geostatistics, developed originally over 30 years ago, has been a spectacularly successful set of tools for estimation of mineral resources and reserves. However, practical success does not necessarily imply a secure theoretical basis. Both quantum mechanics and Christianity can be said to work well in practice, yet many consider the theoretical foundations of both to be problematic.

Geostatistics, as defined by the Regionalised Variable Theory of Georges Matheron, is rooted in the linear statistical model. Effectively it defined one basic estimation method, known as linear kriging (named after Professor Danie Krige, one of the key pioneers of the field). This is essentially a simple weighted moving average method. It is just like the much simpler inverse- power-of-distance moving average, except that weighting factors applied to the data are not simple functions of distance, but are computed from the interpreted covariance relationships among the data points and the point or block to be estimated. Nevertheless, it remains a linear estimation method, and very early it was proved rigorously that - subject to certain assumptions - linear kriging is the best linear unbiased estimator - ‘BLUE'. A nice by-product of linear kriging is that besides the estimate itself it also gives an estimation variance - an indication of the quality of the estimate for each point or block.

If geostatisticians had stopped there, all might have been fine. There were problems, however. The first and biggest is that the spatial covariance relationships are actually unknown for any real data set. They must be estimated by use of a device known as the semivariogram (more commonly known as the variogram), which can be interpreted according to one or other conventional model in order to calculate a set of model covariance relationships. The interpretation of variograms has been found in practice to be almost as difficult as Freud's interpretation of dreams. Geostatisticians have developed a number of ‘rules of thumb' to help, but it remains a highly subjective black art, though the quality of variogram models can, to some extent, be checked by cross-validation (‘jack-knifing') methods in which values at known points are estimated from those around them - and the set of estimates compared with the set of original data.

The estimation variance which in the early days was proclaimed by geostatisticians as one of the unique advantages of their method, has more recently fallen into disrepute. First, it has come to be realised by users that it is derived purely from the variogram and the spatial distribution of data points, and is not directly related to actual data values. This means that it has no dependence on local data properties which might deviate from those represented in the global variogram (for example, in a high-grade pocket in a gold deposit, the estimation variance would be just the same as in a low-grade ‘background' area - this is clearly unrealistic, since the predictability of grades will be expected to be much poorer in higher-grade areas. Second, when using any of the nonlinear methods it has no physical meaning at all, since it cannot in general be back-transformed to data units in the way that the estimates themselves can.

Another problem is that of stationarity - the assumption that the mean grade in a deposit, for example, is the same everywhere. Of course the nature of many economic mineral deposits is such that this often will not be true - the miner is mostly interested in ore bodies which have rich areas that are economically mineable and can be distinguished from poor areas that are not. Linear kriging does not work very well where there is an element of systematic variation over the region to be modelled. Although globally it remains unbiased, nevertheless the local estimates tend not to be so good. An early attempt to resolve this problem was known as universal kriging. It incorporated polynomial functions of the spatial coordinates into the kriging equations. These polynomial functions were exactly the same as those used in classical trend surface analysis, and the idea was to define a deterministic surface on which the geostatistical regionalised variable was superimposed. The problem with this method was that the semivariogram itself was sensitive to the form of the deterministic surface. Therefore, it required a number of iterations of kriging, variogram computation, and model fitting, to converge towards a consistent solution. More recent work has been based on the theory of ‘generalised covariances' in which progressively higher order derivatives of the data have been used in attempts to define the surface form. These methods have found little practical application because of their complexity, and the inherent instability of their solutions. Furthermore, the resulting kriging system is no longer linear and thus loses its ideal ‘BLUE' properties.

Another problem in geostatistics is that of the data distribution. Although regionalised variable theory does not require normal (gaussian) distribution of the data, it does assume normal distribution of error terms. Unfortunately, in most precious and base metal deposits, the data distribution is far from normal - and the error distribution itself cannot then be assumed to be normal. Even if the error distribution is normal, the use of linear kriging on any data set will tend to under-estimate the peaks and over-estimate the troughs - it is after all a moving average method with an essential smoothing effect. For positively skewed distributions such as those of gold deposits, this smoothing of the peaks can become serious and result in serious under- estimation of higher grade areas. Various ways have been developed to try to counter this problem. Lognormal kriging was developed to handle data that were lognormally distributed, as many gold deposits appeared to be. This was, of course a nonlinear method and hence lost its ‘BLUE' properties. However, worse than that, it was found that the data distribution in most deposits of gold, or anything else, were sufficiently different from lognormal to cause substantial errors in the estimates obtained from lognormal kriging. Another family of methods that has been developed uses empirical transformations to convert data to a normally distributed form before kriging, and then back-convert afterwards. Although the kriging operation then works fine, the use of the nonlinear transformation and back-transformation again entirely destroys the ‘BLUE' properties. Furthermore, these methods are difficult for the practising mining geologist or engineer to understand and to use.

A method which has been evolved for practical use on data populations with any distribution is known as indicator kriging. This is easier to understand. A cutoff grade is defined. All data below the cutoff are assigned a value of 1, all data above the cutoff are assigned 0. Variogram modelling and kriging is done on this set of ones and zeroes. If this is repeated for a selection of different cutoff grades, grade estimates can be obtained by combining the set of estimated indicators at each point or block. Although this method has been conspicuously successful in practice, it is again nonlinear and hence non-BLUE. The practice of using what was developed as a linear statistical model for continuous data, to operate instead on discrete data (ones and zeroes), is in my view analogous to using a carpenter's chisel to break rocks. You can do it, but it's not very good for the chisel. There have also been a number of practical problems found in applying indicator kriging, mostly fixed by ad hoc work-arounds. It is very much an engineer's method - it works well as long as you don't think too hard about its theoretical underpinning.

Possibly the root cause of all of these problems - and the necessity to develop the wide range of esoteric variants of the geostatistical method - is the limitation imposed by Matheron's regionalised variable theory. Brilliant though it was, it was rooted in the ‘standard' linear statistical model that evolved in the pre-computer age. The proliferation of geostatistical methods is reminiscent of the mediaeval epi-cycles progressively built into the Ptolemaic model for planetary motion to account for discrepancies from the assumed circular orbits. Maybe the time is ripe to search for a more comprehensive - and more objective - theory for spatial data estimation.

Stephen Henley
Matlock, England
steve@silicondale.co.uk
www.SiliconDale.com