5.4 Kriging
Kriging is a method devised by geostatisticians to provide the best local estimate of the value of the mean value of a regionalized variable. Typically this was an ore grade and was motivated by the desire to extract the most value from the ore deposit with the minimum amount of capital investment. The technique and theory of geostatistics has grown since those early days into a field dedicated to finding the best linear unbiased estimator (BLUE) of the unknown characteristic being studied.
Simply put, a regionalized variable is a variable that can be said to be distributed in space. This space is not limited to the three-dimensional kind of space that we move around in every day, but can be extended to include time, parameter space, property space, etc. This definition "distributed in space" is purely descriptive and makes no probabilistic assumptions. It merely recognizes the fact that properties measured in space follow an irregular pattern that cannot be described by a mathematical function. Nevertheless, at every point in a space x it has a value. A regionalized variable, then, seems to have two contradictory characteristics:
Hence we are dealing with a naturally occurring property (variable) that has characteristics intermediate between a truly random variable and completely deterministic variable. In addition, this variable (property) can have what is known as a "drift" associated with it. These drifts are generally handled with trend surface analysis and can be analyzed for and subtracted out of the data much the same way an offset can be subtracted out of a data set.
While one can use ANOVA to determine at what level a trend surface is "significant", there is still something of an art to determining the correct variogram model to use. Using ANOVA as a guide one can fit trend surfaces to the data of every increasing order, eventually your ANOVA will tell you that at some specified level of significance, a trend surface of that order does not provide a statistically significant increase in the fit to the data. That's where you stop, at order minus one.
The weights and neighborhood of the trend surface analysis is dependent upon the semivariance of the data, i.e. dependent upon the structure function the data displays. This interdependency between trend surface and semivariance means that there is no unique solution/combination of trend and semivariance, hence the art. The degree of spatial continuity of the data (regionalized variable) is given by the semivariogram (see section 5.3) and some of the types of models used is covered in section 5.3.2.
To explain what kriging is I am going to concentrate on the simplest form of kriging, punctual kriging. Consider that you want to find:
That is to say, find the best linear weighted (unbiased) estimate for property (P). In addition suppose you also want to know what the estimation error is as well:
There is a way to do this by requiring that the weights wi sum to one, this will result in an unbiased estimate if there is no trend. You can then calculate the error variance as:
It seems only logical that the closer a data point is to the grid point you wish to estimate the more
weight it should carry. These weights used (wi) and the error of estimate (
) are related to YP
through the semivariogram. So, if we had three data points from which to estimate one grid
point, we would have:
for the estimate and:
for the weights. The question that remains is: how do we find the best set of wi's? Consider Fig. 5.4.1, here we have three data (control) points and from them we wish to make a best linear unbiased estimate of the field at grid point P.
Using the semivariogram we can create the following sets of equations:
where 
is the semivariance over the distance h between control points i and j and 
is
the semivariance over the distance h between the control point i and the grid point P. With Eqn.
5.4.5 we have three unknowns and four equations and to force Eqn 5.4.5 to always be true we
add a slack variable 
resulting in a matrix set of equations like:
This yields the wi and one more equation:
yields the error of estimate.
Now I want to point out that something really cool is happening here. If you stop and think about it you may wonder why the weights should apply to both the data points and the semivariances. We shouldn't have any problem considering equation 5.4.4, after all it's just the best linearly weighted combination of the surrounding data points. But what about equation 5.4.6? Why should these also be true? Well, strictly they aren't, not until you add the slack variable (lambda) that allows equation 5.4.5 to always be true. What insight does this give you into the nature of regionalized variables? I'll let you ponder that for a while.
Now it sometimes happens that you don't want to or can't remove the trend surface prior to kriging. It is still possible to come up with a best linear unbiased estimate of your grid points using Universal Kriging, the matrix you form is even more complicated than the one in Eqn. 5.4.7 and is covered in Davis, Chapter 5.
This example is taken from Davis, Chapter 5 and addresses only the concept of punctual kriging. Suppose you wanted to dig a well and wanted a good estimate of the depth of the water table before you began digging. Suppose further that you had three wells already dug distributed about your proposed site (p) much in the same fashion as are the control points in Fig. 5.4.1. Given the data in the table below, you can use punctual kriging to make a best linear unbiased estimate of the water table elevation at your proposed site.
| X1 Coordinate (km) | X2 Coordinate (km) | Water Elevation (m) | |
| Well 1 | 3.0 | 4.0 | 120 |
| Well 2 | 6.3 | 3.4 | 103 |
| Well 3 | 2.0 | 1.3 | 142 |
| Your site p | 3.0 | 3.0 | ??? |
From this information and a structure analysis (semivariogram) you can fill out the equations in equation 5.4.7 and then solve for the water table elevation at your site. The semivariogram analysis revealed a linear semivariogram out to 20 km with an intercept of zero and a slope of 4 m2/km. So the matrices in equation 5.4.7 look like:
The numbers on the lefthand side of the equation come from the semivariances between control points calculated by knowing the distance between them and the linear semivariogram model. The numbers on the righthand side of the equation come from knowing the distance between the proposed site and each control point and the linear semivariogram model.
If the condition number of the matrix on the lefthand side isn't too bad, you can invert directly to solve for the W's and lambda, otherwise you can use SVD to solve for the answer. Either way in this case you get a column vector of:
which when multiplied through equation 5.4.4 gives an estimate of 125.3 m and using equation 5.4.8 yields an error estimate of 5.28 m2. The square root of this number is represents one standard deviation (2.30 m), which represents the bounds of 68% confidence. So plus or minus two times this standard deviation yields the depth of the water table at your proposed site p with a 95% confidence.
MatLab's answers are a little different from the ones in Davis (1986), but I attribute that to the fact that Davis does not use singular value decomposition to invert his matrices (see Davis, 1986, Chap 3).
GoTo Next Section