19 Kriging
19.1 Theory and equations of kriging
Kriging is an interpolator in which to the unknown geospatial point \(\left({\mathbf u}\right)\) is assigned a value \(\left(z({\mathbf u})\right)\) given by the linear combination \(\left(\lambda_{\alpha}\right)\) of measured values \(\left(z({\mathbf u}_{\alpha})\right)\) in the neighbourhood of the unknown point \(\left({\mathbf u}\right)\).
\(\hat{z}({\mathbf u}) = \displaystyle\sum _{\alpha=1}^{n({\mathbf u})} \lambda_{\alpha}({\mathbf u}) z({\mathbf u}_{\alpha})\) \(\;\;\) kriging interpolation
Kriging is just a weighted moving average.
Weights \(\left(\lambda_{\alpha}\left(\mathbf{u}\right)\right)\) assigned to measured points are unknown and are calculated by solving the kriging system of linear equations, given in matrix form in the following equations:
\(\lambda_{\alpha}(\mathbf{u}) = K^{-1} k\) \(\;\;\) where \(\;\;\) \(K^{-1} = (K^{T}K)^{-1}K^T\) \(\;\;\) kriging system
where:
\(\lambda_{\alpha}(u) = \begin{bmatrix} \lambda_1({\mathbf u}) \\ \vdots \\ \lambda_n(u) \end{bmatrix}\) \(\;\;\) vector of weights
\(K = \begin{bmatrix} C(u_1-u_1) & \cdots & C(u_1,u_n) \\ \vdots & \ddots & \vdots \\ C(u_n-u_1) & \cdots & C(u_n-u_n) \end{bmatrix}\) \(\;\;\) variance-covariance matrix, calculated using measurements at point locations
\(k = \begin{bmatrix} C(u_1-u) \\ \vdots \\ C(u_n-u) \end{bmatrix}\) \(\;\;\) measured vs unknown point covariance, calculated from the variogram model
We have to constrain the weights to sum up to one (non-bias condition1):
\(\sum_{\alpha=1}^{n({\mathbf u})} \lambda_{\alpha}({\mathbf u}) = 1\)
Indeed, in the particular case when all data values are equal to a constant, the estimated value should also be equal to that constant (se Eq. 11.1, Wackernagel).
The \(n({\mathbf u})\) weights \(\lambda_{\alpha}({\mathbf u})\) are calculated to ensure that the estimator is unbiased (the non-bias condition) and that the estimation variance is minimal. To satisfy the non-bias condition, the expected error must be null:
\(E\left[\hat Z\left(\mathbf u\right) - Z\left(\mathbf u\right) \right] = 0\).
The estimation variance (also called the variance of the error) is:
\(\sigma^2 = var\left[ \hat{Z}\left(\mathbf{u}\right) \right] = E\left\{ \left[] \hat{Z}\left(\mathbf{u}\right) - Z\left(\mathbf{u}\right) \right]^2 \right\}\)
which is minimum setting the partial derivative to zero:
\(\partial \left( E\left\{ \left[ \hat{Z}\left(\mathbf{u}\right) - Z\left(\mathbf{u}\right) \right]^2 \right\} + Lagrangian \right) = 0\)
This procedure provides the kriging system shown above.
19.1.1 Random Function and Regionalized Variable
| RAF | Z(x) | → | Z(x₀) | The Random Variable’s |
| ⇓ | ↓ | ↓ | realization is a | |
| REV | z(x) | → | z(x₀) | Regionalized Value |
19.1.2 Stationarity
If the regionalised variable \(z({\mathbf u})\) is stationary it is also intrinsic (but the reverse is false), and the variogram is bounded:
\(\gamma(h) = C(0) - C(h)\)
- \(C(0)\) is the covariance for a zero distance or variance
- \(\gamma(h)\)
19.1.2.1 Types of Stationarity for Random Functions
stationarity of random function \({\mathbf Z}({\mathbf u})\): it’s distribution is (all moments are) invariant under translation
(weak or) second order stationarity of \({\mathbf Z}({\mathbf u})\): the first two moments (mean and covariance) are invariant (i.e. constant) under translation
the intrinsic hypothesis assumes that the increments of the \({\mathbf Z}({\mathbf u})\) function (i.e. the increments \(Z(x + h) – Z(x)\)) are weakly stationary.
(See M. Armstrong, Basic Linear Geostatistics)
19.1.3 Ordinary kriging
For the ordinary kriging (OK) system we can also write (see H. Wackernagel Multivariate Geostatistics; P. Goovaerts Geostatistics for Natural Resources Evaluation):
\(\begin{bmatrix} \gamma(u_1-u_1) & \cdots & \gamma(u_1,u_n) & 1 \\ \vdots & \ddots & \vdots & \vdots \\ \gamma(u_n-u_1) & \cdots & \gamma(u_n-u_n) & 1 \\ 1 & \cdots & 1 & 0 \end{bmatrix} \times \begin{bmatrix} \lambda_1(u) \\ \vdots \\ \lambda_n(u) \\ \mu_{OK}(u) \end{bmatrix} = \begin{bmatrix} \gamma(u_1-u) \\ \vdots \\ \gamma(u_n-u) \\ 1 \end{bmatrix}\) \(\;\;\) OK system using dissimilarity
where \(\mu_{OK}(u)\) is the Lagrange parameter2 and under the constraint:
\(\sum_{\alpha=1}^{n({\mathbf u})} \lambda_{\alpha}({\mathbf u}) = 1\)
That is, in matrix form:
\({\mathbf\lambda}_{\alpha}(\mathbf{u}) = {\mathbf\Gamma}^{-1} \; {\mathbf\gamma}\)