versión On-line ISSN 2411-9717versión impresa ISSN 0038-223X

J. S. Afr. Inst. Min. Metall. vol.114 no.3 Johannesburg mar. 2014

School of Mining and Petroleum Engineering, Department of Civil and Environmental Engineering, University of Alberta, Alberta, Canada

SYNOPSIS

Block model estimates are commonly calculated by the well-established technique of kriging. The mathematics are well established, but implementation details require site-specific considerations. Checking and measuring the performance of the estimates is important to ensure the block model is suitable for its intended purpose. There are different criteria for long-term planning and short-term planning. The implementation of kriging should be checked with cross-validation and assessed for conditional bias and departure from theoretical optimality with the calculation of kriging efficiency. A new expression of kriging efficiency, which compares the kriging variance with the theoretically optimal kriging variance, is developed to aid in this assessment. The assessment tools presented here can be applied in a number of situations; however, ordinary kriging with a reasonably large search performs well in most cases.

Keywords: block model, kriging, variance.

Introduction

The challenge of spatial prediction is common in the earth science and other disciplines including meteorology and hydrology. Estimating spatially related geological variables is difficult; estimates must be made given samples taken from only one billionth of the deposit volume or less. Early computerized kriging techniques were developed by Krige (1951), Sichel (1973) and others in South Africa. Matheron (1962, 1963) formalized the methods for mining resource and reserve estimation and named them kriging after the South African pioneer Danie Krige. There were parallel developments in forestry (Matérn, 1960), bathymetry (Switzer et al., 1963), meteorology (Gandin, 1965), and mathematics (Kolmogorov et al., 1962). Further discussion of the historical development of kriging is given by Cressie (1990).

The properties of kriging are well known; it is a best, linear, unbiased estimator that minimizes the variance of estimation errors; often referred to as the estimation variance or kriging variance (Journel and Huijbregts, 1978). In practice, many decisions are required to make a kriging estimate such as kriging type, search parameters, and data selection. The resulting estimates should then be checked and assessed. In this paper, we examine kriging and how necessary decisions should be made and assessed with metrics and measures of performance. There are three types of estimates, each requiring a different strategy and criteria for assessing results: (1) estimates for visualization and geological understanding, (2) estimates for interim planning, and (3) final estimates for reserve classification.

Review of kriging framework

Consider a regionalized random variable Z (a realization is referred to as the 'data'). The data is subset into reasonable geologic rock types or domains called stationary populations in geostatistics. The data is assembled:

The notation ua denotes a particular sample location within the domain. Estimates at unsampled locations are required. The data is usually at a nominal point scale, whereas the unsampled locations are volumes relevant for mining reserves and resources assessment. That is, they are blocks of a much larger size. The block estimate Z* at an unsampled location u is computed by a linear combination of n nearby data:

The notation V denotes the larger block volume. The weights applied to each of the nearby data are denoted by λi. The stationary mean is m, and the usual assumption is that this applies throughout the domain. The challenge is to choose the nearest n data to combine and to choose the weights. One evident choice is to retain 4 to 20 data and assign weights by inverse distance interpolation to the unsampled location. This can work well; but there are nuances of spatial structure and anisotropy that can be exploited for even better estimates.

Consider a separation vector h defined by an anisotropic distance and direction (two angles, usually azimuth and inclination from horizontal). The spatial variability for h is represented by the semivariogram y(h). The semivariogram (commonly referred to as the variogram) is defined as the half expected squared difference for the random variable at data points separated by h:

The variogram is experimentally estimated using the JV(h) data separated by a vector h within some defined tolerance:

The experimental variogram is fit with a parametric variogram model. This permits calculation of the closely related covariance function C(h) = σ2 - y(h) where σ2 is the variance of the data for the domain. The local weights for the estimate are optimized to minimize the variance of the errors of estimation (estimation variance):

Where C(ui,uj) is the covariance between two data locations ui and uj, C(ui;, ui) is the average covariance between the block centered on u and the data at ui and C(V,V) is the average covariance value between any two points within the block variance (Journel and Huijbregts, 1978). Minimizing this function by taking the partial derivatives and setting them equal to zero results in the simple kriging (named by Matheron) equations which are solved to determine the optimal weights:

It should be noted that for simple kriging, there is no constraint that the sum of the weights equals unity and that the simple kriging equation can be expressed differently to include the weight given to the mean; this expression can be found in any standard geostatistical text (Journel and Huijbregts, 1978). The corresponding minimized estimation variance (kriging variance) can then be calculated:

Choosing to constrain the weights to sum to unity leads to the ordinary kriging equations (so named by Matheron). The error variance is minimized subject to the unbiasedness constraint that the sum of the weights is equal to unity (Journel and Huijbregts, 1978). This relaxes the assumption of global stationarity since the mean is effectively estimated using the local search data. Other variants of kriging such as universal kriging or intrinsic random function kriging can also be considered for specific cases. Complete development of these kriging methods can be found in many geostatistical texts (Chiles and Delfiner, 2012; David, 1977; Journel and Huijbregts, 1978). For the estimation of mining variables in practice, ordinary kriging is used more often than simple kriging or alternative formulations, although integrated simulation approaches are also being adopted (Godoy and Dimitrakopoulos, 2011), and simple kriging is used in areas of wide data spacing.

Decisions for kriging and assessment of results

When kriging is used for estimation, there are a number of significant decisions to be made. The type of kriging, generally either simple or ordinary kriging, must be decided. A restricted search is often considered to reduce the computational burden of computing and inverting large covariance matrices and, in the case of ordinary kriging, to reduce reliance on the hypothesis of a stationary mean (Rivoirard, 1987). A restricted search for simple kriging may also be used to eliminate the presence of negative weights and therefore decrease the weight applied to the mean (Boyle, 2010). As such, the next step is data selection and the definition of search parameters. Search parameters commonly include a maximum range around the location being estimated to search for local data and a maximum number of local data to consider. Other search parameters could include a maximum number of data to use from the same drill-hole and a maximum number of data to use from each octant of the local neighbourhood (Deutsch and Journel, 1998). This data pooled together composes the kriging neighbourhood (Rivoirard, 1987; Vann et al., 2003).

Rivoirard (1987) suggests the use of the weight given to the mean in simple kriging and the slope of regression for ordinary kriging as diagnostics on the search parameters. If the weight given to the mean in simple kriging with a limited neighbourhood is large, then a larger neighbourhood (less restrictive search) can be used. Conversely, if the weight applied to the mean in simple kriging is small, then the local neighbourhood has a strong influence. From Rivoirard's analysis, it is almost always necessary to use a larger neighbourhood for ordinary kriging compared to simple kriging to reduce conditional bias. The slope of regression of the true values against the estimated values is used as a diagnostic for conditional bias. Ideally, the slope of this regression would be as close to unity as possible. This diagnostic depends heavily on the stationarity of the domain, as many of the diagnostics presented here do.

The large number of decisions required for kriging requires careful assessment of the resulting estimates. The most widely used method for assessing kriging results is cross-validation. Data locations are estimated with local data, excluding data from the same drill-hole. The estimated values are compared to the true values. Ideally, the estimates would correlate very strongly with the truth. In practice, the variance of the estimates is lower, so on average, high values are underestimated and low values are overestimated. A more robust cross-validation method is jackknifing, where a subset of the data is omitted from the beginning of modelling, and all calculations (estimation of the mean, variance, and variogram) are performed without this subset of the data.

Cross-validation and jackknifing are excellent assessment methods, but can be used only for locations where we know the true value. For the majority of locations where we do not know the true value, kriging metrics such as the slope of regression and kriging efficiency can be considered. The slope of regression provides a measure of the conditional bias in the estimates. The kriging efficiency as defined by Krige (1997) is akin to a goodness of fit. In this paper, a novel definition for kriging efficiency estimates the departure from the theoretical minimum squared estimation error. The results of cross-validation and these metrics can be used to guide decisions regarding data selection and kriging type.

Data selection

Search anisotropy and ties

The determination of search anisotropy for kriging is not trivial. The decision on search anisotropy parameters is guided by the variogram and data spacing. The variogram provides information on the spatial variability law of the domain, and the data spacing is used to choose a practical range. The search anisotropy should align with the variogram anisotropy, but be modified by the data spacing. As each of the variogram structures can have a different orientation, an effective anisotropy across all nested variogram structures should be considered. A reasonable distance for the search anisotropy should then be used, such as two data spacings. Other considerations can include restricting the search to a maximum number per octant or maximum number per drillhole. Imposing additional restrictions may be necessary for very complicated deposits or deposits with long-range geologic continuity.

Care should be taken in the case of evenly spaced drillhole data when a limited number of data are used for the kriging estimate. In this case the data used in the estimation (if they are equidistant) will depend on the kriging program and original order in the file. This can lead to inconsistencies between estimates if the data order in the file changes. This is a concern in mining, where establishing an audit trail and the reproducibility of estimates are critically important. One solution is to pre-sort the data along an arbitrary vector. The positions of data at tied distances are then chosen in order of the pre-sort vector. Alternatively, equidistant data could always be included if the situation arises.

A related issue when using only a subset of the data for kriging is the generation of discontinuous kriging surfaces. For some applications such as hill shading, hydrological modelling, and visualization, this is undesired (Gribov and Krivoruchko, 2004; Rivoirard and Romary, 2011). Gribov and Krivoruchko propose the use of a smoothing kernel which is applied to the weights of distance data. Rivoirard and Romary use a penalizing function to drive distance data weights to zero. These approaches are useful when the kriging application calls for it, but the resulting histogram of estimates should be checked.

Number of search data

The choice of the number of local search data to use for the kriging estimate is very important. As discussed earlier, the number of search data is commonly limited to reduce computational time and reliance on a global stationary mean. With increasing computer speed, this first concern is relatively minor. These kriging decisions were evaluated by the authors using a series of case studies similar to the evaluation of nonlinear kriging methods by Moyeed and Papritz (2002). A series of case studies were carried out to determine the effect of the number of search data on the mean squared error of the estimates. For each test case, kriging was carried out in cross-validation mode considering both simple and ordinary kriging with between 5 and 100 local data used to make each estimate. The mean squared error between the estimates and true value were then compared. In all cases, data from the same drill-hole as the estimate being made was excluded as is typical in cross-validation studies.

The first case study was a low-grade porphyry copper deposit. This data set includes copper data from 134 drillholes with an average spacing of 100 m. The vertical copper data was composited into 3 m sections. The copper assays are moderately skewed with a mean of 0.25% copper and standard deviation of 0.27%. The variograms for this deposit are approximately isotropic and well behaved with a nugget effect of 20% of the sill. Figure 1 plots the effect of the number of search data on simple and ordinary kriging results using this copper data. For a low number of search data, simple kriging performed better than ordinary kriging, but with a large number of search data, ordinary kriging was the better estimator.

The second case study is a set of bitumen data from the Athabasca oil sands. The percentage bitumen was sampled at 280 drill-holes with the vertical data composited into 3 m sections. Due to the stratified geology of the oil sands, this deposit displays very strong vertical to horizontal anisotropy of approximately 150:1. The bitumen histogram is approximately uniform with a mean of 7.7% bitumen. As with the copper case study, ordinary kriging performed better than simple kriging for a large number of search data (Figure 2).

The third and final case study considered was a zinc deposit with zinc assays from 367 drill-holes. For this deposit, the metal assays were skewed, and there was a moderate amount (approximately 3:1) of anisotropy between the horizontal and vertical directions observed in the variograms. For this deposit, the mean squared errors for ordinary or simple kriging with large numbers of search data were comparable (Figure 3).

Metrics

Given the large number of choices that must be made by the practitioner, we are motivated to consider measures and metrics of how kriging performs. Clearly, the number one measure is the minimized mean squared error variance (kriging variance); however, if search parameters are not used there will be less emphasis on local stationarity. In such cases other measures related to conditional bias and departure from theoretical optimality should be considered.

Slope of regression

Consider cross-plotting the unknown true value of a random variable Z with estimation volume V against the known estimate of the random variable Z* on the same volume V. The regression of the true values given the estimate is an indication of the conditional bias in the estimate (Rivoirard, 1987). Specifically, the slope of regression, b, approximates the conditional bias of the kriging estimate:

If the bivariate relationship between the truth and estimate was Gaussian, then the linear regression would be exactly the conditional expectation. Even in non-Gaussian settings, the slope of regression is a reasonable approximation of the conditional bias. The slope of regression is calculated using the estimation weights:

For simple kriging, the slope of regression is exactly 1. In the case of ordinary kriging where the Lagrange multiplier is used, the slope of regression is generally less than 1. Given this, we could say that ordinary kriging is conditionally biased unless the slope of regression is equal to 1.

Kriging efficiency

Kriging efficiency was first introduced by Krige (1997) as a measure of the efficiency of block estimates. The kriging efficiency is expressed as the kriging variance (σ2K) normalized by the variance of the true blocks (σ2) as a percentage (Equation [10]). For block kriging, we consider the estimates to be made on some volume V, and the corresponding variance of the blocks is equal to the average covariance between points within the blocks: C (V,V). We express Krige's kriging efficiency, KEDK, subscripted by Krige's initials:

A high efficiency means that the kriging variance is low, and the variance of the block estimates is approximately equal to the variance of the true block values. A low efficiency implies a high kriging variance relative to the block variance. The kriging variance varies from block to block, so the kriging efficiency will vary as well. There are a number of limiting cases discussed by Krige (1997). For perfect valuations, the efficiency is 100%:

Where all blocks are valued at the global mean (global estimate of all blocks is the only estimate made):

With no conditional bias for imperfect valuations:

With a conditional bias for imperfect valuations:

Krige also notes that the efficiency can be negative if the kriging variance is greater than the true block variance. When the estimation variance exceeds the block variance, Krige deems this a kriging anomaly and states that valuing the block with the mean would be more efficient assuming we know the mean accurately.

Statistical efficiency

Efficiency is a measure of the relative amount of effort to accomplish a task. If a different process can accomplish the same task with less effort, then it is more efficient. In statistics, the efficiency of a statistical quantity is defined differently depending on the property. The three most commonly considered are the efficiencies of (1) an estimator, (2) an experimental design matrix, and (3) a hypothesis test. For an unbiased estimator, the efficiency is defined as the minimum possible estimation variance divided by the variance of the estimator (Fisher, 1922). The minimum possible estimation variance is determined by the sample size; the minimum possible estimation variance will change depending on whether we have 10, 100, or 1 million samples. This minimum possible estimation variance is given by the Cramér-Rao bound, which states that the minimum possible variance of the estimator is the inverse of the Fisher information matrix (Rao, 1945). The most efficient estimator (if one exists) will have an efficiency of 1. Less efficient estimators will have efficiencies in [0,1).

Krige's definition of kriging efficiency differs from the classical definition in that the kriging efficiency is a measure of the variance in the true value that is not represented in the kriged estimate. This is more akin to a goodness-of-fit or proportion of variance explained. The absolute minimum variance would be the global simple kriging variance. If the definition of efficiency proposed by Krige were to be reworked as a statistical measure of efficiency, it should incorporate the global simple kriging variance which is the absolute minimum kriging variance. If the kriging variance is equal to the global simple kriging variance, then the estimator is efficient and will have an efficiency of 1. If the kriging variance is larger than the global simple kriging variance, then the efficiency will be less than 1. This could also be expressed as a percentage. The efficiency will vary depending on the block being estimated and data available. Although global simple kriging is the lowest variance linear estimator, there may exist lower variance nonlinear estimators. Restricting the definition of efficiency to linear, unbiased estimators, a new definition for the kriging efficiency is:

Global simple kriging is the most statistically efficient unbiased linear estimator. It does, however, force the model to rely heavily on a strong global stationarity assumption, which is unrealistic in many cases. For example, issues with an increase in the mean squared error and decrease in the actual slope of regression were encountered by Boyle (2010) when a very large search area was used. The reason for using ordinary kriging is to obtain a local estimate of the mean at the expense of introducing some conditional bias by limiting the search, and therefore reduce reliance on the assumption of a global stationarity mean throughout the whole domain.

Constraining the estimator will introduce a conditional bias. Constrained estimators are sometimes used in statistics because it is possible that an estimator with a small bias will have a smaller mean squared error (the mean squared error is equal to the variance plus the square of the bias). To compare efficiencies, the ratio of efficiencies is taken. This gives a measure of how much more efficient one estimator is compared to a second estimator.

For example, kriging with a more restrictive search radius will decrease the relative efficiency compared to kriging with a large search radius. The ratio of efficiencies would give a measure of how much variance is being incurred to decrease reliance on the global stationarity assumption. Ultimately, knowledge of how much the decision of stationarity could be relied upon would be necessary to determine an acceptable efficiency level (McLennan and Deutsch, 2004).

Application of kriging metrics to case studies

Recall the three case studies introduced earlier: (1) a low-grade copper porphyry deposit, (2) an oil sands deposit, and (3) a zinc deposit. We can consider calculating the kriging metrics discussed for the same range of search data. For simplicity, the global simple kriging variance is approximated using 100 search data. This approximation was found to be very reasonable with the data-sets.

Kriging metrics were calculated for the copper porphyry case study (Figure 4). The mean slope of regression, Krige's efficiency, and kriging efficiency are each given as the mean of the calculated theoretical values. The mean weight given to the simple kriging is also plotted for reference. The shape of Krige's efficiency, the theoretical slope of regression, and kriging efficiency are all similar, with the exception of the simple kriging slope. Recall that for simple kriging, the slope of regression will always be 1. A key aspect of this plot is that as the number of search data increases, the mean kriging efficiency (using the measure introduced in this paper) for both ordinary and simple kriging asymptotically approaches 1. This is because the estimator variances are approaching that of the theoretically best linear estimator: global simple kriging. This contrasts with what is observed for Krige's efficiency, which plateaus around a value of 0.3. This is because the kriging variance is not being compared to the best possible case but to the block variance instead.

For the oil sands case study, the kriging metrics follow a similar shape (Figure 5), although with a distinct elbow at the point (10-20 search data) where the cross-validated mean squared errors displayed an elbow (Figure 2). This elbow indicates that any additional search data used receive only very small weights.

Scatter plots of the estimated values versus the true values are shown for simple and ordinary kriging with 10 and 90 data in Figure 6. As expected, the observed conditional bias for simple kriging with any number of data is very low, as is the conditional bias with ordinary kriging and a large number of data. Note that the metrics presented here, including the slope of regression, are not the average of the theoretical values as they are in Figure 5, but are calculated using the known true values from cross-validation. As the number of search data used for simple kriging increases, the average weight to the simple kriging mean decreases, resulting in an increased standard deviation of the estimates.

Finally, the effect of the number of search data on kriging metrics was quantified for the zinc case study (Figure 7). As with the mean squared error, it can be observed that the simple kriging and ordinary kriging efficiencies quickly converge. In this test case, the slope of regression exceeds 1 for ordinary kriging; it is not common for the theoretical slope of regression to exceed 1, but this has been documented in the literature (Boyle, 2010).

Discussion

There are a number of useful observations that can be gathered from the test cases. Increasing the number of search data will generally decrease the mean squared error for simple and ordinary kriging. This is contrary to the notion that the search should be restricted to relax the assumption of global stationarity and make better estimates. Increasing the number of search data for ordinary kriging will implicitly increase the accuracy of the estimate of the local mean and improve the estimates made. In the case studies examined here, using a large number of data resulted in estimates with a lower mean squared error when examined with cross-validation. For the majority of the deposit where the truth is not known, kriging metrics, including kriging efficiency and the slope of regression, provide a reasonable proxy for assessing the departure from theoretical optimality. If increasing the number of search data results in a large increase in kriging efficiency, or the slope of regression, then increasing the number of search data will also be likely to reduce the mean squared error of the estimate. This is observed when the effects of the number of search data on the mean squared error and kriging metrics are compared for each of the case studies.

In both the oil sands and copper case studies, ordinary kriging with a large number of search data resulted in a lower observed mean squared error compared with simple kriging. In the zinc case study, the results of ordinary and simple kriging with a large number of data were comparable. Theoretically, simple kriging will always result in a lower mean squared error compared with ordinary kriging. The better performance of ordinary kriging relative to simple kriging is attributed to better estimation of a local mean compared to the global mean used by simple kriging. If the mean is not globally stationary in the domain, then using a local stationary mean with ordinary kriging will result in better estimates. This is supported by the plots of the mean weight given to the simple kriging mean in the case studies (Figures 4, 5, and 7). In the copper (Figure 4) and oil sands (Figure 5) case studies, the simple kriging mean receives weights of up to 0.1 when a very large search neighbourhood is used. In these cases simple kriging was outperformed by ordinary kriging. Conversely in the zinc case study (Figure 7), the simple kriging mean receives a negligible weight when a large search neighbourhood is used. In this case the results of simple and ordinary kriging were almost identical. Using a very large number of search data resulted in a very accurate estimate of the local mean. Increasing the number of search data did not decrease the accuracy of this estimate, as very distant data points received negligible weight.

Increasing the number of search data had the expected effect on the kriging metrics; it decreased the conditional bias, which was accessed through the slope of regression, and increased the kriging efficiency as the kriging variance approached the global simple kriging variance. When the number of search data is restricted to reduce smoothing and produce estimates that more closely match with the histogram predicted by volume-variance relationships, a conditional bias will be introduced into the estimates. For the purposes of estimating resource and reserves, this conditional bias is generally accepted since the cost of using a very smooth kriged map is often drastically under- or overestimating of reserves unless the estimates are post-processed with a local conditioning approach such as uniform conditioning.

Ultimately, the kriging search strategy depends on the estimate type. Estimate types can be broadly classified as (1) visualization and geological understanding, (2) interim estimates, or (3) final estimates. When the goal is visualization and gaining a level of geological understanding, a very smooth map is desired since this is easier to understand. For this purpose global simple kriging (to prevent any search artifacts) could be considered. An alternative approach would be to use local kriging with a smoothing method (Gribov and Krivoruchko, 2004; Rivoirard and Romary, 2011) applied.

Interim estimates are used for long-term planning. The goal is to anticipate the information effect for the future and consider volume variance relations. For interim estimates a restricted search to increase the variance of the estimates could be considered. This increases understanding of what variability can be expected in the future. This approach is straightforward and commonly implemented, but not the only available approach. Emery (2006) compared the use of Monte Carlo integration with the point-support multigaussian model and the commonly used discrete Gaussian model for use with ordinary multigaussian kriging. Either method could be used as an alternative to evaluate ore deposit reserves without introducing a bias.

Final estimates come down to the decision of ore or waste. These decisions must be made with minimum conditional bias and should have a minimum mean squared error. From the results of the case studies presented, the practitioner should consider the use of ordinary kriging with a large number of search data. A series of cross-validation case studies could be done with the data-set to confirm that this was the best approach. These results are similar to the results obtained by Boyle (2010), who demonstrated that a large number of samples increased the accuracy of estimates, although this improvement was marginal. Boyle emphasized that dividing the area into suitable domains is likely to be far more important than the selection of the precise number of data to use when estimating, a suggestion that we echo here. The incorporation of additional data to decrease estimation error and conditional bias must be done within a stationary domain.

Conclusions

There are a number of useful kriging measures, including Krige's efficiency, the slope of regression, and the statistically-based kriging efficiency proposed in this paper. The proposed kriging efficiency considers the estimation variance compared to the global simple kriging variance, which is the theoretically minimum possible estimation variance. In the case studies considered, the lowest mean squared error in cross-validation was normally for ordinary kriging with a very large number of search data. In this case, the kriging efficiency was close to 1. The kriging metrics used in this paper can be used to ascertain the departure from theoretical optimality of a kriging implementation. If the kriging efficiency or slope of regression is significantly lower than 1, then the mean squared error and conditional bias will likely be high. This may be acceptable if the goal of the estimate is to assess the proportion of the deposit above a certain cut-off grade, in which case the smoothing introduced by using a very large number of search data may be unacceptable. Ultimately, the decision of how many search data to use and search strategy depends on the type of estimate being made. The kriging metrics and guidelines introduced in this paper should assist the practitioner in making these decisions.

Fisher, R.A. 1922. On the mathematical foundations of theoretical statistics. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, vol. 222. pp. 309--368. [ Links ]

Gandin, L.S. 1965. Objective analysis of meteorological fields. (translated from Russian), Israel Program for Scientific Translations for the US Department of Commerce and the National Science Foundation. [ Links ]