Geostatistical Ore Reserve Estimation – Quick Tutorial

The development of geostatistics as an ore reserve estimation methodology emerged in France in early 1960 from the work of Matheron (1962) and was based on original studies by D.G. Krige functioning the optimal assigning of weights to the neighbouring sample values used in estimating the grade of blocks in South African gold mines. Geostatistics can be defined as the application of the theory of “regionalized variable” to the evaluation of a mineral deposit, involving the study of the spatial relationship between sample values, thickness or any geological phenomena showing intrinsic dispersion.

The primary purpose of any natural resource estimation is to reliably estimate the overall ore resource and the distribution of insitu and recoverable tonnage and grade throughout a mineral deposit. Conventional methods e.g. area of influence, polygonal, other geometrical, distance weighting etc. may provide a global estimate of an ore body’s reserve. However, meaningful geostatistical reserve study with careful attention to geological controls on mineralization would provide an adequate indication of relative confidence in the block grade estimate.
Geostatistics aims at providing quantitative descriptions of natural variables distributed in space or in time and space. Examples of such variables include:

1. Ore body parameters in a mineral deposit
2. Depth and thickness of a geological layer
3. Porosity and permeability in a porous medium
4. Density of trees of a certain species in a forest
5. Soil properties in a region
6. Concentration of pollutants in a contaminated site.

Geostatistical estimation involves analysis of Regionalised variable (Rev). A Regionalised variable has magnitude (support or volume in space) as well as its spatial location (in reference of coordinate axes) and spatial correlation. It treats data to be distributed in space. The geostatistical estimation techniques are based on a set of theoretical concepts known as the theory of Regionalised Variables developed by Matheron (1971) based on empirical work carried out by Krige (1951, 1952 and 1962). Most regionalised variables, in ore reserve estimation, display two aspects; viz.

1. a random aspect, consisting of highly irregular and unpredictable variations and
2. a structured aspect, reflecting spatial characteristics of the regionalised phenomena.

The Geostatistical ore reserve estimation involves the following steps:

1. Geostatistical structural analysis through semi-variography that involves:
   i)  Construction of an experimental semi-variogaram; and
   ii) Fitting a model/mathematical Semi-variogram to it;
2. Deciding BLOCK dimension and delineating mineralized boundary on each slice;
3. Block kriged estimation employing (Ordinary) Kriging;
4. Mineral inventory estimation
5. Establishment of grade-tonnage relations: and
6. Test for Conditional non-bias & smoothing effect.

Semivariography

The function that measures the spatial variability among the sample values is called ‘semi-variance’. Semi-variograms are constructed by comparing a sample value with the remaining ones at constantly increasing distance called lag interval (Isaaks and Srivastava, 1989; Kim, 1991; and Armstrong, 1998). The mathematical formulation of a semi-variance function γ (h) is given as:

where Z (x_i) is the value of regionalised variable (e.g., grade) at a point xi in the space; Z(x_i+h) is the grade at another point at a distance ‘h’ known as lag distance/lag interval and ‘N’ is the number of sample pairs being considered.

Fig.1 One-dimensional data arrays illustrating how data are paired for the purpose of constructing an experimental semivariogram.

An experimental semi-variogram is a plot between lag-distance and semi-variance and it provides the following information regarding the characteristics of a mineral deposit:

1. A measure of continuity of mineralisation (C);
2. Range (a); A measure of the zone or area (range or zone) of influence.
3. Sill (C0+C);
4. Nugget (C0); In theory the semivariogram value at the origin (0 lag) should be zero. If it is significantly different from zero for lags very close to zero, then this semivariogram value is referred to as the nugget.
5. The nugget represents variability at distances smaller than the typical sample spacing, including measurement error.
6. A measure of trend, where grade (gx) = Trend (mx) + Residual (Rx); and
7. A measure of Anisotropy-(Geometric/elliptical and Zonal/stratified).

Fig.2 Characteristics of the Semivariogram.

Semi-variography with due consideration to deposit geology is able to quantify the characteristics of spatial continuity via nugget effect, range, sill and directional anisotropy, which in turn, provides an adequate model of geological influences that are used in reserve estimation. An experimental semi-variogram is first constructed using the raw data. Then a model semi-variogram is fitted to it by (i) Hand fit method (ii) Non-linear least squares fit method; and (iii) Point Kriging Cross-Validation (PKCV) method. This fitted model provides parametrisation of the characteristics of the deposit and it is utilized during kriging process.

To an experimental semi-variogram, various mathematical models may be fitted such as Spherical/Matheron model, Exponential model, De Wijsian/Logarithmic model, Linear model, parabolic model, Hole-effect model, Mixed/Nested Spherical model etc. Spherical/Matheron model is most frequently used as more than 95% of the mineral deposits conform to this model (Rendu, 1981).

The geostatistical procedure of estimating values of a regionalized variable using the information obtained from a model semi-variogram is Kriging. It is an optimal spatial interpolation technique. It is called BLUE as (i) it is BEST (because of minimum estimation variance); (ii) LINEAR (because of weighted arithmetic average); (iii) UNBIASED (since the weight sum to unity) Estimator.

Let G* be the kriged estimate of the average value of grade G of the samples having values g₁, g_2,, g₃, …..,g_n. Let w₁, w₂, w₃,….. w_n be the weightage given to each of the values respectively such that Σw_i = 1; and G* = Σw_ig_i. Thus the estimation becomes unbiased; the mean error is zero for a large number of estimated values and the estimated variance is minimum. The Kriging variance is given as σk² = Σ (g_i – G*)². To make kriging variance minimum, a function called Lagrange multiplier (λ), is used for optimal solution of the kriging system.

POINT KRIGING

Point kriging is a method of estimation or interpolation of a point by a set of neighbouring sample points applying the theory of regionalised variables where the sum of weight coefficients sum to unity and produce a minimum variance of error.
Expressed mathematically, kriged estimate is given as:
P*= Σw_i s_i
where P* = the estimate of true value at a point ‘p’.
w_i = weight coefficients of the individual samples
s_i = individual sample values at sample points, s_i.
and kriging variance, σk² = Σ w_iγ(s_i,p) + λ
where λ = Lagrangian multiplier and γ(s_i,p) = average  semi-variance among
samples and the point to be estimated.

BLOCK KRIGING

It is a method of estimation of a block of ground with the help of surrounding sample values using the theory of Regionalised variables. The kriged estimate G* of a block is mathematically expressed as:
G* = Σwi gi

where G*= estimated value of the block using a set of sample S_i.; w_i = weight  coefficients and n = number of samples used for estimation of block.
Kriging variance, σk² = Σ (g_i – G*)², is mathematically given as: σk²  = Σ w_i γ (S_i, V) – γ (V, V) + λ,
where
γ (V, V) = Average semi-variance within block V.
γ (S_i, V) = Average semi-variance between sample, S_i and whole of the block, V and
λ = Lagrange multiplier, introduced minimization processes to balance the number of equations with the number of unknown coefficients.
The weight coefficients, wi, and the Lagrangian multiplier, are computed from the matrix form of kriging equation, which is given as:

Fig.3 The block and data illustrating how various average gamma values are obtained in order to krige a single block. Lines extending from Sample 3 illustrate some of the pairs of data used to estimate the average gamma value between s3 and all other data points; lines extending between discretetized block locations (+signs) illustrate how pairs are selected to determine the average gamma value within the block; lines from sample s7 to discretized block locations (+signs) illustrate how paired values are obtained to determine the average gamma value between a sample and the block.

The methods of kriging described here viz., Point Kriging and Block Kriging belong to linear geostatistics. The non-linear geostatistics deals with Lognormal Kriging (Rendu, 1979 and 1981; Sinclair and Blackwell, 2002), Disjunctive Kriging (Matheron, 1976) and Multi Gaussian Kriging (Verly, 1983), while the non-parametric geostatistics includes Indicator Kriging (Journel, 1983; Sinclair and Blackwell, 2002) and Probability Kriging (Sulivan, 1984). Additionally, there are other models such as Universal Kriging (Kriging in the presence of trend, Journel and Huijbregts, 1978; Deutsch and Journel, 1997), Co-Kriging (Kriging of one variable based on the correlation of it with the other variable, Journel and Huijbregts, 1978); Polygonal kriging and Blast hole Kriging (David, 1988; Kim, 1993).

Outline of steps for performing Block Kriging

The entire mineralized body is divided into regularly spaced horizontal sections, by projecting the sample data from the (vertical) cross sections earlier constructed. The vertical height or gap between the sections is kept at length equalling the vertical lift or bench height as per the method of mining. In each of the horizontal sections, the mineralized boundary delineated, is divided into smaller grids based on selective mining unit (SMU). Usually at least one fourth of the drill spacing (for a square grid) is taken as the side of a grid. Each slice forms a set of X and Y arrays of blocks with constant Z values (X-Easting, Y-Northing, Z- Elevation). The arrays of blocks are then kriged slice by slice, producing kriged estimates and kriging variance for each of them.

At first step, the following input parameters are required for block Kriging:

1. A minimum of 4 and a maximum of 16 samples to krige a block, ;
2. The radius of search for sample points around a block centre should be within the range of influence,
3. The semi-variogram parameters: nugget variance (C0), transition variance or continuity and range (a),
4. The ratio of anisotropy in case of anisotropic semi-variogram model, and
5. The dimension of blocks to be kriged and block coordinates.

The next steps that follow include:

1. Computation of average variability of sample values contained within the dimensions of small blocks;
2. Selection of nearest samples lying within the radius of search;
3. Counting the number of the samples. If found insufficient with reference to a minimum specified to krige a block, the next block is taken up and procedure is repeated from step ii);
4. Establishing kriging matrices and computation of weight coefficient;
5. Multiplication of weight coefficient by their respective sample values to provide kriged estimates. Kriging variance is calculated from the sum of the products of the weight coefficient and their respective sample-block variances. An extra constant, called lagrange multiplier is added to minimise the kriging variance; and
6. Move to next block and repeat the procedure from step (ii).
   The individual slices are then averaged to produce a global estimate of kriged mean together with associated variance.

Mineral inventory

The kriged blocks are then stacked slice-wise one below the other from top to bottom provide a 3D array of regularly spaced gridded blocks with estimated values. Development of such an array of blocks provides an estimate of the total stock of mineral in place and is called mineral inventory. It provides the number of blocks kriged, total tonnage, kriged estimate and kriging variance in respect of each slice and finally provides total tonnage of the deposit as a whole together with mean kriged estimate and mean kriging variance.

Fig.4 A three-dimensional array of blocks designed to approximate the geometry of an ore deposit. The block size is commonly taken as the selective mining unit (SMU), the smallest volume for which selection as ore or waste is possible and thus the smallest volume for which an average grade must be estimated.

Grade-tonnage relations

The next step of the geostatistical evaluation is to produce a series of grade-tonnage estimates at various hypothetical cutoff grades. Generally, a greater tonnage is associated with a lower grade. Progressively higher grades are worked out by increasing the degree of selectivity. A numerical approach involving a step-wise integration of the block grade frequency distribution over a range of grades has been employed to calculate: (i) quantities of ore, metal and waste; (ii) average grade of ore and waste; and (iii) waste to ore ratio at various hypothetical cut-off grades. A plot of these relationships provided grade-tonnage curves.

Advantages of Geostatistics:

Grade control problems to keep mill feed grade fluctuation within predetermined limits can be addressed by this technique.

1. Design of optimum drilling patterns
2. Directional variability
3. Volume- variance relationship
4. 2D and 3D block model with prediction accuracy
5. Ore evaluation for mine planning and production
6. Calculation of the error of estimation of the volume or tonnage of an ore body