Method of Soil Analysis
1.5 Geostatistics
1.5.1 Introduction
1.5.2 Using Geostatistical Methods
1 Dec. 2004
D1 Takeshi TOKIDA
1
1.5.1 Introduction

True understanding of the spatial variability in
the soil map is very limited.
 Distinct
boundary (too continuous or sudden change).
 Assumption of uniformity within a mapping unit is not
necessarily valid.


Spatial and temporal variability diversify our
environment. It’s Benefit!
However Soil variation can be problematic for
landscape management.
2
1.5.1 Introduction

There is a need to study surface variations
in a systematic manner.
Geostatistical methods are used in a
variety of disciplines.
e.g. mining, geology, and recently biological sciences also.
Numerous books have been published.
3
1.5.1.1 Geostatistical Investigations
Geostatistics is used to…



map and identify the spatial patterns of
given attributes across a landscape.
improve the efficiency of sampling
networks.
identify locations in need of remediation.


Disjunctive kriging→Probability map
predict future effects in the landscape.

Random field generation→Conditioned→Predict
4
1.5.2 Using Geostaitstical Methods
1.5.2.1 Sampling

Consider the appropriate sampling methodology (see
Section 1.4)
Analysis
Objective of
the study
Appropriate data
collection
The analysis of the data depend on the objective of the study and
appropriate data collection.
5
Table 1.5-1


If the Kolmogrov-Smirnov statistic is greater than the
critical value, the hypothesis of “not being normal” is
adopted.
If the distribution is completely normal, skew and
kurtosis values are 0.
6
Table 1.5-2
?
θ
θ
ρb
Na
ln(Na)
B
ln(B)
ρb
Na
1 -0.49
ln(Na)
0.06
1 -0.14
1
0.14
B
ln(B)
0.09
0.19
-0.25 -0.02 -0.16
0.67
0.58
0.53
1
0.57
0.83
1
0.76
1
Na values can be used to estimates B content at lower cost.
7
Randome function & realization

Observed data are a single realization of
the random field, Z(x).
Z(xα)
Realization
+
Assumptions, i.e.
stationarity
Random field
(Random function)
8
1.5.2.2 Spatial Autocorrelation
Only if a spatial
correlation exists,
geostatistical analysis
can be used.
 Fig. 1.5-1 A: No
spatial correlation
 Fig.1.5-1 B: spatially
correlated

Fig. 1.5-1
9
1.5.2.2.a Variogram
Variogram
Experimental variogram (Estimator)
How to create pairs?
10
Variogram model
Variance is undefined
Var(Z)
95%
Practical range
11
Important considerations when calculating the variogram 1
Between a lag interval, in this case
1.5 to 4.5, a wide range of actual
separation distance occurs.
Imprecision
compared with a situation where every
sampling pair has the same distance
A large number of pairs are used to
calculate a variogram value.
It is generally accepted that
30 or more pairs are sufficient to
produce a reasonable sample
variogram.
Fig. 1.5-3
12
Important considerations when calculating the variogram 2
Fig. 1.5-4



Width of the lag interval can affect the variance.
This is not the case.
The value for h (actual separation distance) is
affected by the lag width.
13
Fig. 1.5-1 & 1.5-5

Fig. 1.5-1
The variograms reproduce
spatial structure of
simulated random fields.
14
Example of variogram
Sill
Sill
Nugget effect
Range


Some information at the smaller scales (less than
48 m) has been lost.
For both attribute, the range is about 900 m.
15
1.5.2.2.d Directional Variograms
Geometric anisotropy
Fig. 1.5-7


Fig. 1.5-8
Often there is a preferred orientation with higher
spatial correlation in a certain direction.
For many situations, the anisotropic variogram
can be transformed into an isotropic variogram
by a linear transformation.
16
1.5.2.2.e Stationarity
A sample at a location
Impossible to determine the probability
distribution at the point!
Assumption:
The joint distribution do not depend on the
location.
A stationary Z(x) has the same joint probability distribution
for all locations xi and xi+h.
17
Second-Order Stationarity
Autocovariance
E[ Z (x  h)]  E[ Z (x)]  
cov[Z (x  h), Z (x)]  E[(Z (x  h)  E[ Z (x  h)])(Z (x)  E[ Z (x)])]
 E[(Z (x  h)   )(Z (x)   )]
 E[ Z (x  h) Z (x)]   2
 C (h)
1.5-3, 1.5-6
＊
C(0)
 (h)  C (0)  C (h)
(h)
1
2
1
1
 E[(Z (x  h))2 ]  E[ Z (x  h) Z (x)]  E[(Z (x))2 ]
2
2
1
1
 C (0)   2  C (h)   2  C (0)   2
2
2
 C (0)  C (h)
 (h)  E[(Z (x  h)  Z (x))2 ]

Sill
Nugget
effect

C(h)
h
Range
 
 

18
Intrinsic Stationarity (Hypothesis)
E[ Z (x  h)  Z (x)]  m(h)  0
var[Z (x  h)  Z (x)]  2 (h)
No Drift
Theoretical Variogram
1
 (h)  E[( Z (x  h)  Z (x)) 2 ]
2
Drift?
No Drift?
Fig. 1.5-9
If lim 
h 
 (h)
h
2
 0, the random field is stationary in terms of Intrinsic hypothesis.
19
1.5.2.2.c Integral Scale
1.5-5
A measure of the distance for which the attribute
is spatially correlated.
1.5-4
Autocorrelation function:
normalized form of the autocovariance function
20
1.5.2.3 Geostatistics and Estimation
Kriging produces a best linear unbiased
estimate of an atribute together with
estimation variance.
 Multivariate or cokriging: Superior
accuracy
 Powerful tool, useful in a wide variety of
investigations.

21
1.5.2.3.a Ordinary Kriging
We wish to estimate a value at xo using
the data values and combining them
linearly with the weiths: λi
xo
1.5-7
Z* should be unbiased:
n
E[ Z * (x 0 )  Z (x 0 )]   i EZ xi   
i 1
n
n

   i  1  0   i  1
i 1
 i 1

1.5-9
22
Derivation of equation 1.5-10
Z* should be best-linear, unbiased estimator.
Our goal is to reduce as much as possible the
variance of the estimation error.
First, rewrite the estimation variance



 E Z (x )  Z (x )  

var Z * (x 0 )  Z (x 0 )  E ( Z * (x 0 )  Z (x 0 ))  E Z * (x 0 )  Z (x 0 )
0
2
0
2
*
 
0
2
 n
 
 E   i Z (x i )  Z (x 0 )  
 
 i 1
 n
 n

n

2
 E   i Z (x i )   i Z (x i )   2 E  i Z (x i ) Z (x 0 )  E Z (x 0 ) 
 i 1

 i 1

 i 1

n

i 1
   EZ (x )Z (x ) 2  EZ (x )Z (x )  EZ (x ) 
n
j 1
n
i
j
i
j
i 1
2
i
i
0
0
23

Derivation of equation 1.5-10
Let’s rewrite the estimation variance in terms of the semivariogram.
We assume intrinsic hypothesis.
From the definition of the semivariogram we know:



1
 (x i  x j )  E Z (x i )  Z (x j )2
2
1
1
2
2
 E Z (x i )   E Z (x i ) Z (x j )  E Z (x j ) 
2
2

 








1
1
2
2
E Z (x i ) Z (x j )   (x i  x j )  E Z (xi )   E Z (x j ) 
2
2
1
1
2
2
EZ (xi ) Z (x 0 )   (x i  x o )  E Z (x i )   E Z (x o ) 
2
2




E Z (x0 )2   (x0  x0 )  E Z (x0 )
2




24
Derivation of equation 1.5-10
Just substitute:


n
var Z (x 0 )  Z (x 0 )  
*
i 1

2








E
Z
(
x
)
Z
(
x
)

2

E
Z
(
x
)
Z
(
x
)

E
Z
(
x
)
 i j
 i
i
j
i
0
0
n
n
j 1
i 1





2
2








(
x

x
)

E
Z
(
x
)

E
Z
(
x
)
 i j  i j 2
i
j

2
n

n

i 1
1
j 1


1


1
1

2
2 
 2 i   (x i  x o )  E Z (x i )   E Z (x o )  
2
2


i 1
n

  (x 0  x 0 )  E Z (x 0 ) 
n
2
n
 
    (x
i 1
i
j 1
1 n
 
2 i 1
j

n
i
 x j )  2 i (x i  x o )   (x 0  x 0 )
i 1

 


2








E
Z
(
x
)

E
Z
(
x
)


E
Z
(
x
)
 i j
 i
i
j
i
n
j 1
2
2
n
i 1

0
    EZ (x )  EZ (x ) 
n
n
i 1
j 1
2
i
j
i
2
j







n
n
n
n
n

2
2 
2 
2
  i   j E Z (xi )     j E Z (x j )      j E Z (x j )  1   i   2  j E Z (x j ) 
i 1
j 1
j 1
 i 1 
 j 1
 j 1
n
1
1
25

Derivation of equation 1.5-10
We define an objective function φ containing a term with the
Lagrange multiplier, 2β.


 n

 i ,    var Z (x 0 )  Z (x 0 )  2   i  1
 i 1

n
 
i 1
*
 n




(
x

x
)

2


(
x

x
)


(
x

x
)

2



1





i j
i
j
i
i
o
0
0
i
j 1
i 1
 i 1

n
To solve the optimization
problem we set the partial
derivatives to zero:
n
 i ,  
 0 for i  1,...n
i
 i ,  
0

26
Derivation of equation 1.5-10
Ordinary Kriging system
n
 2  j  (x i  x j )  2 (x i  x o )  2
j 1
n
   (x
j 1
j
n

j 1
j
i
 x j )     (x i  x o ) for i  1,2,...,n
equation 1.5-10
1
Example:
27
Derivation of equation 1.5-10
Kriging Variance


n
var Z (x 0 )  Z (x 0 )  
*
i 1
n
    (x
j 1
i
j
n
n
i
 x j )  2 i (x i  x o )   (x 0  x 0 )
i 1
n
  i  (x i  x o )    2 i (x i  x o )   (x 0  x 0 )
i 1
i 1
n
   (x
j 1
j
i
 x j )     (x i  x o )
n
  i (x i  x o )     (x 0  x 0 )
equation 1.5-12
i 1
Block Kriging
Estimation of an average value of a spatial attribute over a region.
Average variogram values
Variance
equation 1.5-13
equation 1.5-14
equation 1.5-15
28
1.5.2.3.b Validation
Cross validation
Little bias
Estimated kriging
variance is nearly equal
to the actual estimation
error.
29
1.5.2.3.c Examples
equation 1.5-18
Isotropic Case, Kriging Matrix.
equation 1.5-10
1.5-11
But we can’t find the values of a given attribute!
λ1=0.107, λ2=0.600, λ3=0.154, λ4=0.140
Note that the weight for point 1 is less than point 4, even though the
distance from the estimation site is almost the same.
30
Creating Maps Using Kriging
Directional variogram oriented
in 0°& 90°
Length of each ray is equal to the range
of the directional variogram.
Anisotropy ratio = major axes / minor axes
31
Creating Maps Using Kriging
Fig. 1.5-12 Based on
Anisotropic variogram
Fig. 1.5-13 Based on
isotropic variogram
32

Bootstrapping high dimensional vector: interplay between