A Non-Homogeneous Model for Kriging Dosimetric Data

Abstract

This paper deals with kriging-based interpolation of dosimetric data. Such data typically show some inhomogeneities that are difficult to take into account by means of the usual non-stationary models available in geostatistics. By including the knowledge of suspected potential sources (such as pipes or reservoirs) better estimates can be obtained, even when no hard data are available on these sources. The proposed method decomposes the measured dose rates into a diffuse homogeneous term and the contribution from the sources. Accordingly, two random function models are introduced, the first associated with the diffuse term, and the second with the unknown and unmeasured source contribution. Estimation of the model parameters is based on cross-validation quadratic error. As a bonus, the resulting model makes it possible to estimate the source activity. The efficiency of this approach is demonstrated on data simulated according to the physical equations of the system. The method is available to researchers through an R-package provided by the authors upon request.

Appendix A: Impact of Covariance Bias on Kriging

Here, we demonstrated that using the biased covariance

$$\begin{aligned} \tilde{K}\,=\, K+\varPhi ww^{t}\varPhi ^{t}, \end{aligned}$$

affects neither the drift nor the dose estimates:

• Drift coefficient estimation. Based on the $\tilde{K}$ matrix, the weights $\tilde{\varLambda }$ are obtained by

  $$\begin{aligned} (K+\varPhi ww^{t}\varPhi ^{t})\,\tilde{\varLambda }\,=\,-\varPhi \tilde{\mu }. \end{aligned}$$

  Since $\varPhi ^{t}\tilde{\varLambda }=I$ (universality condition), this gives

  $$\begin{aligned} \tilde{\varLambda }=-K^{-1}\varPhi \,(\tilde{\mu }+ww^{t}), \end{aligned}$$

  (5)

  and finally

  $$\begin{aligned} \varPhi ^{t}K^{-1}\varPhi (\tilde{\mu }+ww^{t})\,=\,-I. \end{aligned}$$

  But the Lagrange multiplier $\mu$ of the true kriging system satisfies $\varPhi ^{t}K^{-1}\varPhi \,\mu =-I$, so provided that the solution is unique, we identify $\tilde{\mu }+ww^{t}=\mu$. Reporting in (5), we check that indeed $\tilde{\varLambda }=\varLambda$. We conclude that despite the bias on the covariance, we get the correct kriging estimate of the drift. However, note that the apparent kriging variance

  $$\begin{aligned} \tilde{\varLambda }^{t}\tilde{K}\tilde{\varLambda }=\varLambda ^{t}\tilde{K}\varLambda =\varLambda ^{t}(K+\varPhi ww^{t}\varPhi ^{t})\varLambda =\varLambda ^{t}K\varLambda +ww^{t}=\sigma _{K}^{2}+ww^{t}, \end{aligned}$$

  is not the correct one.

• Interpolation. The perturbed kriging system is now

  $$\begin{aligned} \left( \begin{array}{c@{\quad }c} K+\varPhi ww^{t}\varPhi ^{t} &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) \,\left( \begin{array}{c} \tilde{\varLambda }\\ \tilde{\mu } \end{array}\right) \ =\ \left( \begin{array}{c} K_{0}+\varPhi ww^{t}\varPhi _{0}^{t}\\ \varPhi _{0} \end{array}\right) . \end{aligned}$$

  We decompose the perturbed kriging matrix as follows:

  $$\begin{aligned} \left( \begin{array}{c@{\quad }c} K+\varPhi ww^{t}\varPhi ^{t} &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) \,=\,\left( \begin{array}{c@{\quad }c} K &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) +\left( \begin{array}{c@{\quad }c} \varPhi ww^{t}\varPhi ^{t} &{} 0\\ 0 &{} 0 \end{array}\right) . \end{aligned}$$

  Then, provided that $K^{-1}$ and $(\varPhi ^{t}K^{-1}\varPhi )^{-1}$ exist, we can write the perturbed kriging weights as the solution of

  $$\begin{aligned} \left[ I+\left( \begin{array}{c@{\quad }c} K &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) ^{-1}\left( \begin{array}{c@{\quad }c} \varPhi ww^{t}\varPhi ^{t} &{} 0\\ 0 &{} 0 \end{array}\right) \right] \,\left( \begin{array}{c} \tilde{\varLambda }\\ \tilde{\mu } \end{array}\right)= & {} \left( \begin{array}{c} \varLambda \\ \mu \end{array}\right) +\, R, \end{aligned}$$

  (6)

  where $\varLambda$ is the unperturbed kriging weights, and the vector R is given by

  $$\begin{aligned} R\ =\ \left( \begin{array}{c@{\quad }c} K &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) ^{-1}\left( \begin{array}{c} \varPhi ww^{t}\varPhi _{0}^{t}\\ 0 \end{array}\right) . \end{aligned}$$

  It is not difficult to establish the following expression for the inverse of the kriging matrix (this formula is based on Schur complements)

  $$\begin{aligned} \left( \begin{array}{c@{\quad }c} K &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) ^{-1}=\left( \begin{array}{c@{\quad }c} K^{-1}-K^{-1}\varPhi (\varPhi ^{t}K^{-1}\varPhi )^{-1}\varPhi ^{t}K^{-1}\ &{} \ K^{-1}\varPhi (\varPhi ^{t}K^{-1}\varPhi )^{-1}\\ \\ (\varPhi ^{t}K^{-1}\varPhi )^{-1}\varPhi ^{t}K^{-1} &{} -(\varPhi ^{t}K^{-1}\varPhi ) \end{array}\right) . \end{aligned}$$

  But since

  $$\begin{aligned} (K^{-1}-K^{-1}\varPhi (\varPhi ^{t}K^{-1}\varPhi )^{-1}\varPhi ^{t}K^{-1})\,\varPhi ww^{t}\varPhi _{0}^{t}=0, \end{aligned}$$

  we get

  $$\begin{aligned} R=\left( \begin{array}{c} 0\\ ww^{t}\varPhi _{0}^{t} \end{array}\right) . \end{aligned}$$

  We also have

  $$\begin{aligned} \left( \begin{array}{c@{\quad }c} K &{} \varPhi \\ \varPhi ^{t} &{} 0 \end{array}\right) ^{-1}\left( \begin{array}{c@{\quad }c} \varPhi ww^{t}\varPhi ^{t}\ &{}\ 0\\ 0 \ &{} \ 0 \end{array}\right) \ =\ \left( \begin{array}{c@{\quad }c} 0 \ &{} \ 0\\ ww^{t}\varPhi ^{t}\ &{}\ 0 \end{array}\right) . \end{aligned}$$

  It follows that the first line of (6) implies $\tilde{\varLambda }=\varLambda$.

The conclusion of this discussion is that the observed covariance can be used in place of the true covariance of $Y_{0}$ in the interpolation of the dose rate and in the estimation of the drift coefficients. However the kriging variance will not be reliable, and we are not aware of any practical way to correct it.