Kriging

Algorithm

Kriging can be seen as a point interpolation which requires a point map as input and returns a raster map with estimations and optionally an error map.

The estimations or predictions are calculated as weighted averages of known input point values, similar to the Moving Average operation.

The estimate to be calculated, i.e. an output pixel value , is a linear combination of weight factors (w_i) and known input point values (Z_i):
= S(w_i* Z_i)

In case the value of an output pixel would only depend on 3 input points, this would read:
= w₁* Z₁ + w₂* Z₂ + w₃* Z₃

Thus, to calculate one output pixel value , first, three weight factors w₁, w₂, w₃ have to be found (one for each input point value Z₁, Z₂, Z₃), then, these weight factors can be multiplied with the corresponding input point values, and summed.

In Moving average, the weight factors are simply determined by the distances of the input points towards an output pixel. In Kriging, however, the weight factors are calculated by finding the semi-variogram values for all distances between input points and by finding semi-variogram values for all distances between an output pixel and all input points; then a set of simultaneous equations has to be solved. When the spherical distance option is used, distances are calculated over the sphere using the projection of the coordinate system that is used by the georeference of the output raster map.

All semi-variogram values are calculated by using a user-specified semi-variogram model (based on the output of the Spatial correlation operation). The weight factors are calculated in such a way that the estimation error in each output pixel is minimized.

The optional error map contains the standard errors of the estimates.

Process Simple Kriging:

Find the valid input points:

input points which coordinates are undefined are ignored,
input points which value is undefined are ignored,
handle duplicates or coinciding points as specified by the user (no, average, first value).

Determine the distances between all valid input points (n) and find the semi-variogram value for these distances:

for each combination of 2 input points, i.e. a point pair, the distance between the points is determined,
for each combination of 2 input points, the distance value is substituted in the user-selected semi-variogram model, using the user-specified nugget, sill, and range parameters; this gives a semi-variogram value.
the semi-variogram values are filled out in matrix C (as in Equation 1 below),
matrix C is inverted as a preparation for calculations in step 4.

For the first output pixel, determine the distances of this pixel towards all input points, and find the semi-variogram value for these distances:

semi-variogram values are determined using the selected semi-variogram model and its parameters as above,
the semi-variogram values are filled out in vector D (as in Equation 1).

Calculate the weight factors (vector w):

by multiplying the inverted matrix C (result of step 2) with vector D (result of step 3).

The obtained weight factors apply to the current output pixel only.

Calculate the estimated or predicted values for this output pixel:

as the sum of the products of the weight factors and the input point values (Equation 4).

Optionally, calculate the error variance and standard error for this output pixel:

error variance: by multiplying vector w (result of step 4) with vector D (result of step 3),
standard error or standard deviation: as the square root of the error variance, according to Equation 5b.

Consider the next output pixel and repeat steps 3-7, until all output pixels are done.

Ordinary Kriging:

Find the valid input points:

input points which coordinates are undefined are ignored,
input points which value is undefined are ignored,
handle duplicates or coinciding points as specified by the user (no, average, first value).

For the first output pixel, determine the input points (n) which will make a contribution to the output value depending on the specified limiting distance and minimum and maximum number of points:

input points that are farther away from this output pixel than the specified limiting distance are ignored,
if the number of points found within the limiting distance is smaller than the specified minimum nr of points, assign the undefined value to this output pixel,
use only the specified maximum number of points within the limiting distance, and, in case more points are found within the limiting distance than the specified maximum number of points, use only the points that are nearest to this output pixel.

Determine the distances between all input points that will make a contribution to this output pixel (result of step 2), and find the semi-variogram value for these distances.

for each combination of 2 contributing input points, the distance between the points is determined,
for each combination of 2 contributing input points, the distance value is substituted in the user-selected selected semi-variogram function, using the user-specified nugget, sill, and range parameters; this gives a semi-variogram value.
the semi-variogram values are filled out in matrix C below (eq. 1).

Determine the distances of this output pixel towards all input points, and find the semi-variogram value for these distances:

semi-variogram values are determined using the selected semi-variogram function or model and its parameters as above,
the semi-variogram values are filled out in vector D (eq. 1).

Calculate the weight factors (vector w):

by first inverting matrix C (result of step 3),
by solving the set of simultaneous equations.

The obtained weight factors apply to the current output pixel only.

Calculate the estimated or predicted values for this output pixel:

as the sum of the products of the weight factors and the input point values (Equation 4).

Optionally, calculate the error variance and standard error for this output pixel:

error variance: by multiplying vector w (result of step 4) with vector D (result of step 3), according to Equation 5a.
standard error or standard deviation: as the square root of the error variance, according to Equation 5b.

Consider the next output pixel and repeat steps 2-8, until all output pixels are done.

Formulae to calculate weight factors:

The Kriging weight factors of n valid input points i (i = 1, ..., n) are found by solving the following matrix equation:

[ C ]

[ w ]

[ D ]

(1)

This matrix equation can be written as a set of n +1 simultaneous equations:

S_i ( w_i * g (h_ik) ) + l = g (h_pi)	for k = 1, ..., n	(2)
S_i w_i = 1		(3)

where:

h_ik	is the distance between input point i and input point k
h_pi	is the distance between the output pixel p and input point i
g(h_ik)	is the value of the semi-variogram model for the distance h_ik, i.e. the semi-variogram value for the distance between input points i and input point k
g(h_pi)	is the value of the semi-variogram model for the distance h_pi , i.e. the semi-variogram value for the distance between the output pixel p and input point i
w_i	is a weight factor for input point i
l	is a Lagrange multiplier, used to minimize possible estimation error

Matrix C thus contains the semi-variogram values for all combinations of valid input points that will make a contribution to the output pixel value.

Vector w thus contains the weight factors for all valid input points that will make a contribution to the output pixel value.

Vector D thus contains the semi-variogram for an output pixel and all combinations of valid input points.

Equation (3) guarantees unbiasedness of the estimates. The solutions w_i minimize the Kriging error variance s².

Formulae to calculate an estimate or predicted value for an output pixel:

= S_i(w_i* Z_i)

(4)

where:

	is the estimate or predicted value for one output pixel to be calculated
w_i	is the weight factor for input point i
Z_i	is the value of input point i

Formulae to calculate error variance and standard error:

The error variance is calculated as:

s² = S_i ( w_i * g (h_pi)) + l

(5a)

The standard error or standard deviation is the square root of the error variance, thus:

s = Ö ( S_i ( w_i * g (h_pi)) + l )

(5b)

where:

s²	is the error variance for the output pixel estimate
s	is he standard error or the standard deviation of the output pixel estimate
h_pi	is the distance between the output pixel p and input point i
g (h_pi)	is the value of the semi-variogram model for the distance h_pi , i.e. the semi-variogram value for the distance between the output pixel p and input point i
w_i	is a weight factor for input point i
l	is a Lagrange multiplier, used to minimize possible estimation error

Notes:

The contents of matrix C depends on the semi-variogram model selected by the user, its parameters nugget, sill and range, and the geometric distribution of the points within the limiting circle in the input map.
The contents of vector D is determined by the location of the estimated pixel value with respect to the surrounding input points (inside the limiting circle) and the semi-variogram.
The estimates are computed as linear combinations of the n point sample values with the weights w_i as coefficients (the w_i are found from equation (1)). Therefore the estimates are called 'linear predictors'.
Equation (3) guarantees unbiasedness of the estimates. The solutions w_i minimize the Kriging error variance s².
Equation (5) does not contain the sample attribute information. This means that the error variances solely depend on the spatial distribution of the samples and not on their measurement values (the attribute values).
In the case of Simple Kriging, it is assumed that all input points contribute in some way to the estimate in each pixel. The Kriging matrix has thus a constant value for all pixels estimated and needs to be inverted only once; however the right hand-side D keeps changing.
In Ordinary Kriging the number of points used (n <= N) and hence the size of the Kriging matrix (n+1) will change from pixel to pixel while calculating the output map(s). Hence it is theoretically possible that for each output pixel a new Kriging system of order n+1 has to be solved. The algorithm also takes care that for each new set of surrounding input points, this set is sorted according to distance from the estimated pixel in order to enable to select the closest points satisfying the max nr of points condition.

References:

Clark, I. 1979. Practical geostatistics. Applied Science Publishers, London. 129 pp.
Cressie, N.A.C. 1993. Statistics for spatial data. Wiley, New York. 900 pp.
Davis, J. C. 1973. Statistics and data analysis in geology. Wiley, New York. 646 pp.
Deutsch, C.V., and A.G. Journel. 1992. Geostatistical software library and user's guide. Oxford University Press, New York. 340 pp.
Isaaks, E. H., and R. M. Srivastava. 1989. An introduction to applied geostatistics. Oxford University Press, New York. 561 pp.
Journel, A. G. and Ch. J. Huijbregts. 1978. Mining geostatistics. Academic Press, London, 600 pp.
Krige, D.G. Two-dimensional weighted moving average trend surfaces for ore-valuation, in Proc. Symposium on Mathematical Statistics and Computer Applications in Ore Valuation: Journ. South African Inst. Mining and Metallurgy, Johannesburg, 1966, Mar. 7-8, pp. 13-38.
Matheron, G.F. Principles of geostatistics: Economic Geology, 1963, vol. 58, pp.1246-1266.
Meer, F. D. van der. Introduction to geostatistics. ITC lecture notes. 72 pp.
Olea, R.A. 1991. Geostatistical glossary and multilingual dictionary. Oxford University Press, New York.
Stein, A. 1998. Spatial statistics for soils and the environment. ITC lecture notes. 47 pp.