GEOSTATISTICAL

INTERPOLATIONS
HISTOGRAM MODELING
GAUSSIAN ANAMORPHOSIS

2
Gaussian Anamorphosis
A continuous variable Z could be considered as the
transform of a standard Gaussian variable Y (mean 0,
variance 1), through an anamorphosis function Φ

3
The Gaussian Anamorphosis
Transforming a non Gaussian distribution into a Gaussian
distribution (“anamorphosis” means “transformation”)

F G
1 1

F( z) = G ( y)

0.5 0.5

0 0
0 10 20 z 30 40 Z -3 -2 -1 0 y 1 2 3 Y

Anamorphosis between an empirical distribution and a normal

distribution (cumulated distribution functions)

4
Point Anamorphosis

Experimental
Experimental
anamorphosis
anamorphosis
The anamorphosis links the
raw distribution to a
gaussian one. Raw values: Z

The experimental
anamorphosis is modeled Model
Model
by a continuous function. anamorphosis
anamorphosis

=Φ
Gaussian values: Y
Hermite polynomials
The experimental distribution is modeled by Hermite
polynomials:

=Φ

Note that : = =

=
Gaussian Anamorphosis Modeling Now

Modeling the histogram with:

- a possibility to disperse the data values and no use of
Hermite polynomials
- performing the Gaussian transform of variables
affected by zero-effect or of indicators which allows
the joint simulation of grade and domain indicators

7
Normal score transform

• Normal score transform of the data is required:

- For improving variograms from data with skewed
histogram;
- To change support and define distribution on different
supports;
- To perform geostatistical simulations

8
HISTOGRAM MODELING
VARIOGRAPHY FOR SKEWED DISTRIBUTIONS
Introduction
In case of skewed distributions or with outliers the
variograms may be difficult to interpret.
300

200
Y (m)

100 8000 217

1393
7000

Variogram : simu lognormal

2893
757 2672
0 6000 1461 21102438
1572 1934 2671
1297
1075 3102
0 100 200 743 1205 19402148 2385
2672 2844 2952
2495
5000 1706 3036 3192
3067
X (m) 2009
2299
2331 2795
2870
1110 1738 2825
2827
18072093 2358
2581 2975
4000
607
0.8 Nb Samples: 837 84473898
324 1550
Minimum: 0.0058 3000 957
0.7 Maximum: 1263.1195
Mean: 17.5658
Std. Dev.: 68.8139
0.6 2000
Frequencies

0.5
1000
0.4

0.3 0
0 10 20 30 40 50 60 70 80 90 100
0.2 Distance (m)
0.1

0.0
0 100 200
simu lognormal

Page 10
Introduction

The structure of correlation hidden by few outliers is

present as can be seen on the pairwise relative
variogram
Pairwise Relative Vario. : simu lognorm

1.00

0.75

0.50

0.25

0.00
0 10 20 30 40 50 60 70 80 90 100
Distance (m)

Page 11
Variographic analysis
In the gaussian discrete model the structure of
ℎ = , +ℎ

And the structure ℎ of

ℎ = , +ℎ

are related through the relationship:

ℎ = × ℎ + × ℎ ++ × ℎ

Page 12
Variographic analysis
Gaussian Anamorphosis with 30 polynomials
1500

1000
simu lognormal

500

1393 1738 1940

1205 1550
1297 1807
1706 2009
1.00 1934
1110 1461
1572
1075
0

Variogram : gaussian
957
898
0.75 757
-5 -4 -3 -2 -1 0 1 2 3 4 5 743
Gaussian values 607
0.50 473
324
217
100 0.25 84

90 32

80
0.00
0 10 20 30 40 50
70
Frequencies (%)

Distance (m)
60 Isatis
Grids/Blocks 1x1(sampling 10x10)
50 - Variable #1 : gaussian
Variogram : in 1 direction(s)
40 D1 :
Angular tolerance = 90.00
30 Lag = 2.0000m, Count = 25 lags, Tolerance = 50.00%
Model : 1 basic structure(s)
20 S1 - Spherical - Range = 30.0000m, Sill = 1

10

0
0 50 100 150 200 250 Variograms on the normal
simu lognormal

score transforms

Page 13
Workflow

Modeling the distribution using the gaussian anamorphosis

Calculating the experimental variogram on gaussian variables

Fitting the experimental gaussian variogram

Performing the transformation from gaussian to raw

variogram

Fitting the experimental raw variogram

Page 14
KRIGING
Limitations on variables
• The variables must be additive
• Data must have same (or similar) support
E.g. length, diameter, volume
The reason comes from the linearity of the kriging interpolator

Hint: Sometimes 2D estimation is preferable to 3D

estimation. Thickness is estimated together with
the accumulation (metal quantity)

16
Linear Interpolation
The estimated value is a weighted sum of values at sampled
locations:

Z = å λα Zα
Z0 * : estimated value
* Za: value at locations a
0 λ a : weight at locations a
α

Z3

Z2
λ3 λ2
Z0
λ1

Z1
17
Regression – Data Configuration
The variable at a given location can be estimated from data:

One data configuration <-> One regression

Good Interpolation Poor Interpolation Extrapolation

This spatial regression built from the

Least Square Method is the Kriging.
18
 Least square method
∗
−

• Unbiased:
∗
− =0
• Optimal: the error between the estimated value and the true but
unknown value is minimized.
∗
−
19
 Global unbiasedness

Kriging is a globally unbiased estimator, i.e. on average

the error between the estimated value Z* and true value
Z0 is zero:
∗
− =0
∗
⇒ =

20
 Global unbiasedness
Let’s introduce a constant value to the estimator:
∗
= +

From the global unbiasedness :

− − =0

⇒ = × 1−

1) If is considered as known everywhere : Simple Kriging, 1 − ∑

is the weight associated to the mean.
2) If needs to be re-estimated locally, the condition to apply on the
weights is ∑ =1 : Ordinary Kriging , then = 0.

21
 Simple Kriging

∗
= + × 1−

A weight is applied to the mean m,

Problem with SK: m is generally not known, it is replaced by an

estimate m*

22
What type of Kriging?

• Simple Kriging • Ordinary Kriging

- Order 2 stationarity ₋ Intrinsic stationarity
₋ Unknown mean
- Known mean
₋ sum of weights must
- No conditions on weights equal 1

∗ ∗
= + × =

Problem with SK: m is generally not =1

known, it is replaced by an
estimate m*.

23
Kriging optimality
• Kriging is said the Best Linear Unbiased Estimator
(BLUE)

• Kriging is the Best estimator because the weights are

calculated in order the variance of estimation is
minimized

∗
− → =0

• The development of the derivatives leads to the

kriging equations that involve the variogram model

24
Variance of error
The variance of error can be calculated using the
covariance:

∗
− = −

∗
− = + − .

∗
− = + − .

25
Minimisation of variance SK
To minimize the variance
∗
• Let the variance : − =
= + − 2.

• Calculate the value of the derivative = 0 i.e. minimum

=

⟹ =
Simple Kriging
System

26
Simple Kriging, covariances
In matrix notations the kriging system is:

⋯
⋮ ⋱ ⋮ × ⋮ = ⋮
⋯

Distances between data points, Distances between data and

pair by pair target point

The value of the minimized variance is called Kriging Variance :

∗
− = = −

27
 Minimization of variance OK

= + + .

To minimize , and estimate the average : ∑ =

The Lagrange technique converts a constrained minimization into an

unconstrained minimization, by adding one variable called the
Lagrange Multiplier m :

= + + . + −

28
Lagrange Multiplier method

: − ∗ = × −

⟺ =
is the Lagrange Multiplier, it takes negative values, but we
expect it to be close to 0.
29
Ordinary Kriging, 2nd order stationarity
Under second order stationarity, minimizing the estimation
variance gives:

Φ + =
=0
⇒ Ordinary
Φ Kriging System
=0 =1

With :
∗
− =
= − −

30
Ordinary Kriging, intrinsic stationarity
Intrinsic stationarity implies the use of the variogram, the
kriging system using the variogram is expressed:

Φ − + =−
=0
⇒ Ordinary
Φ Kriging System
=0 =1

With :
∗
− =
= −

31
Ordinary Kriging, covariances
In matrix notations the kriging system is:

⋯
⋮ ⋱ ⋮ ⋮ × ⋮ ⋮
=
⋯
…

Covariances between data points, Weights

Covariances between data
pair by pair and target point

32
 Ordinary Kriging (points)
In matrix notations the kriging system is:
− ⋯ − −
⋮ ⋱ ⋮ ⋮ × ⋮ ⋮
=
− ⋯ − −
…

Variogram values between data Weights

Variogram values between
points, pair by pair data and target point
The value of the minimized variance is called Kriging Variance :

− ∗
= = −

33
 KRIGING
PROPERTIES OF KRIGING WEIGHTS
 How Kriging Works
γ(h), C(h) Sill = 1 Range = 6
Y Z0
Z1 h0-1=2
1
4 h0-2=4
h0-3=4
3
ì 3 h 1 h3
ï
C(h) = í1- 2 r + 2 r3 if h < r
0.52
2 h1-2=3
1 h1-3=5 0.31 ïî 0 if h ³ r
Z2 Z3 h2-3=4 0.15
0.04
1 2 3 4 5
X 1 2 3 4 5 6 7 8 h

3 h3 1 h33 » 0.31
Z*0

Target-Data Between data

= l1Z1 + l 2Z2 + l3Z3 C1- 2 = C(3) = 1 - ´ + ´ 33
2 6r 2 6r
3 5 1 53 » 0.04
Kriging system: C1-3 = C(5) = 1 - ´ + ´ 3
2 6 2 6
ìå lb Cab + m = Cao "a C2 - 3
3 4 1 4 3 » 0.15
= C(4) = 1 - ´ + ´ 3
ïb 2 6 2 6
í
ï
î
åa
la = 1 3 2 1 23 » 0.52
C0 -1 = C(2) = 1 - ´ + ´ 3
2 6 2 6
3 4 1 4 3 » 0.15
C0 - 2 = C0 - 3 = C(4) = 1 - ´ + ´ 3
2 6 2 6

35
 How Kriging Works
3 h3 1 h33 » 0.31

Target-Data Between data

Z0 C1- 2 = C(3) = 1 - ´ + ´ 33
Y 2 6r 2 6r
Z1 h0-1=2
4 h0-2=4 3 5 1 53 » 0.04
C1-3 = C(5) = 1 - ´ + ´ 3
h0-3=4 2 6 2 6
3
3 4 1 4 3 » 0.15
C2 - 3 = C(4) = 1 - ´ + ´ 3
2 h1-2=3 2 6 2 6

1 h1-3=5 3 2 1 23 » 0.52
C0 -1 = C(2) = 1 - ´ + ´ 3
Z2 Z3 h2-3=4 2 6 2 6

1 2 3 4 5
X 3 4 1 4 3 » 0.15
C0 - 2 = C0 - 3 = C(4) = 1 - ´ + ´ 3
2 6 2 6

-1
C11-1 0.31 C1-3 1
C1-2 0.04 λ1 0.52
C0-1 λ1 C11-1 0.31 C1-3 1
C1-2 0.04 0.52
C0-1
C2-1 C1
0.31 C2-3 1
2-2 0.15 λ2 0.15
C0-2 λ2 C2-1 C1
0.31 C2-3 1
2-2 0.15 C0-2
0.15
x = Þ = x
C3-1 0.15
0.04 C3-2 C13-3 1 λ3 0.15
C0-3 λ3 C3-1 C
0.04 3-2 C1
0.15 3-3 1 C0-3
0.15
1 1 1 0 μ 1 μ 1 1 1 0 1

36
 How Kriging Works

λ1=0.65 Z
Y 0 l1 + l2 + l3 = 0.65 + 0.28 + 0.07 = 1
Z1 h0-1=2
4 h0-2=4 Z 0* = l1Z1 + l2 Z 2 + l3 Z 3
h0-3=4 = 0.65 ´ Z1 + 0.07 ´Z 2+0.28 ´ Z 3
3 λ2= 0.07 λ3= 0.28
2 h1-2=3
h1-3=5 s k2 = C( 0 ) - å la Cao - m
1 a
Z2 Z3 h2-3=4 = C( 0 ) - l1 ´ C( 2 ) - l2 ´ C( 4 ) - l3 ´ C( 4 ) - m
1 2 3 4 5
X = 1 - 0.68 ´ 0.52 - 0.07 ´ 0.15 - 0.28 ´ 0.15 + 0.16
= 0.75

-1
λ1 C11-1 0.31
C1-2 C 1-3 1
0.04 0.52 0.65

λ2 C2-1 C1
0.31 2-2 C 2-3 1
0.15 0.15 0.07
= x =
λ3 C3-1 0.15
0.04 C3-2 C13-3 1 0.15 0.28
μ 1 1 1 0 1 -0.16

37
Properties of Kriging: respects symmetry

Gaussian model with range a=200 m

100 m

200 m

Very continuous model à higher weight at closest points

38
Properties of Kriging: respects symmetry

Spherical model with range a= 200 m

100 m

200 m

Far points get higher weights

39
Properties of Kriging: respects symmetry

Pure Nugget effect model: no « local » estimation

100 m

200 m

Every data point receives the same weight à simple average

40
 Properties of Kriging: accounts for anisotropy

Spherical with isotropic

range a=200 m

Spherical with anisotropic

ranges ax=250 m and
ay=200 m

41
Properties of Kriging : Stable Position
Same model, different data configuration
σ 2k = 0.45 σ 2k = 0.68

Left configuration à more reliable estimate

Shorter distances between data à smaller spherical model

42
Properties of Kriging: declustering effect
Over-sampling

St-Dev : 0.919 St-Dev : 0.935

Two vs. three clustered samples à similar estimation

variance

43
 Negative weights
Theoretically weights can be negative, as the only
constraints applied on the unit sum:
• Negative weights are attributed to samples with very
weak relationship to the target
• They help recovering variability and prevent smoothing
• But they can lead to inconsistent results (negative
grade, above 100%, …)

44
 Properties of kriging: Screen Effect
Adding one extra sample

Spherical (1, 200 m)

The added data point screens

off the distant sample

Perfect screen effect à Rare in 3D

45
Influence of the Nugget Effect

Nugget effect à the weights are more homogeneous

à the negative weights disappear

46
Influence of nugget on estimation
variance
Kriging with a
Spherical Model

Map of Porosity
Kriging Error

Model

47
Influence of nugget on estimation
variance
Kriging with
Spherical Model plus
Nugget: the kriging
error increases

Map of Porosity
Kriging Error

Model

48
Comparing kriging to IDW
Inverse Distance Weighting
• Weights do not reflect the spatial structure of the data
- No reference to a spatial model
- Particularly problematic in high nugget situation

• Estimation error is not minimized

- Greater likelihood of conditional bias

• There is little possibility to control the weighting

- Continuous grade: n large -> close samples have bigger weights
- Less continuous grade: n small -> close samples have smaller weights

• Weights do not depend on the geometrical configuration

- No screening effect
- No declustering effect

49
 BLOCK KRIGING
INTRODUCTION : SUPPORT EFFECT
 Support effect

• Support: volume
- Point support
- SMU (block) support
- Panel support

• Mean similar at each support

• Variance decreases at larger
supports
• Area under curve above a
cut-off
• Is a physical property:
continuity of the
mineralization

51
Support effect – Krige’s relationship

Var Z1 > Var Z2 > Var Z3

Page 52
Support effect – Krige’s relationship

, is the
average point
variogram inside
the block.

It measures the
variability absorbed
inside the block.
= 10 = 10

̅, = 6.09 ̅, = 0.67

Page 53
Support effect – Krige’s relationship

, is the
average point
variogram inside
the block.

It measures the
variability absorbed
inside the block.
= 10 = 10

̅, = 6.09 ̅, = 0.67

= 3.91 = 9.33

For every estimation method the support effect is taken into account through
the block variance calculated by the Krige’s relationship:

= − ̅
( , )

Page 54
 Block Variance
ℎ : ( )= ( )− ̅
( , )

For one example block:

Long range: Short range:

̅
( , )→0 ̅, → ( )

( )→ ( ) ( )→
 How to calculate ?
̅ , is the average variogram value within one block.
To calculate it, the block is discretized:
To sample every distance, a random secondary
set of 3x3 points is used :

From ℎ we get the

variogram value ℎ .

Ex: discretization 3x3 Each variogram value between a regular point and a
random point is assessed. The average of these values
is ̅ , . For a 3x3 discretization :

, = ∑ ∑ ,

Page 56
Experimental ( , )
How many points of discretization ?

Guidelines:
• For a continuous variogram, a small number of discretisation
points is sufficient to converge.
• For an erratic variogram, a higher number is advised.
Support Change Consequences
l A sample taken inside a block may be quite different
from the block value.

l Point estimation methods (NN, IDW, Point Kriging) will

falsely increase the importance of a data close to the
block centroid, and create biased estimates.
BLOCK KRIGING
KRIGING SYSTEM
Block kriging

We want to estimate the average value of ∗ of the

block from the data :

∗ =

60
Block kriging

• Because of the support effect, the average block grade

is not simply the point grade of the block centroid

• The volume of the block is « v » (SMU, mining block,

Panel, etc…)

61
Block kriging
The variability between one data and one block is given
by a point-block variogram. It is derived automatically
for each configuration:

( )

( )

, = ,

62
 Ordinary Kriging (blocks)
In matrix notations the kriging system is:

− ⋯ − − ,
⋮ ⋱ ⋮ ⋮ × ⋮ ⋮
=
− ⋯ − − ,
…
Variogram values between data Weights Variogram values between
points, pair by pair data and target block

63
Block kriging variance
The Kriging block variance has the
following expression:

= , − , −

In terms of covariance, we have:

= − , − ,

ℎ = − ,

Note that: , =

64
 Example Block estimation
= − , −

Simple Kriging

= , = , =

Ordinary Kriging

=− = = , = , = −

Pure Nugget Effect Spherical (30m) 65

 Block kriging with no discretization

Estimating the block centroid with

no discretization, is point kriging

Legend
red squares are selected samples for the
estimation
values above are the kriging weights,
values below are data values

66
Block kriging with discretization

4x4 discretization
The weights are calculated differently

Legend
red squares are selected samples for the
estimation
values above are the kriging weights,
values below are data values

67
Block kriging with discretization

• Implicit assumption: the « sub-block » size is

comparable to the composite support size on which the
variogram has been established

• The variability is already averaged by regularizing along

the drill hole and the composite length

Support length (vertical) should equal block dimension (vertical)

68
 Block kriging with discretization
Example 1:
₋ long range variogram

₋ vertical composites of 2 m

₋ horizontal blocks of 6 m x 6 m with a vertical extent of 2 m

à discretization of 4x4x1 is adequate

Example 2:
₋ small range variogram

₋ vertical composites of 2 m

₋ horizontal blocks of 10 m x 10 m with a vertical extent of 4 m

à discretization of 6x6x2

69
Polygon kriging

Blocks do not need to be

regular-irregular polygons
(e.g. production blocks)
can also be kriged

The POINT-to-BLOCK
covariance must be
calculated for each
polygon, by discretization

70
 KRIGING
KRIGING PROPERTIES
Indications of kriging quality
Characteristics of estimated points or blocks to be checked:

• Conditional bias
slope of regression E[Z|Z*]

• Precision
estimation variance or standard deviation sk

• Adequacy of the data configuration

weight assigned to the mean in SK

• Smoothness
dispersion variance of the kriged estimates Var[Z*]

72
Kriging outputs
• Three linked random variables are used to define Kriging
outputs : Z, Z* and Z-Z*

• Their relationship is characterized by the variogram: ℎ =

− ∗

--> We have access to probabilistic parameters, although we

only have one value of
∗
Z*.
− = = +

= =

∗
= − −
∗
, = − −

∗ , ∗
, =
∗

73
Kriging system
∗

∗
−
= ∗

≥ ∗

∗
−
∗
,
∗
,
∗
∗
( , )
= ∗

∗ ∗ = 0
− ∗=
/

74
Kriging system -> refaire le nuage !
∗

∗
−
= ∗

≥ ∗

∗
−
∗
,
∗
,
∗
∗
( , )
= ∗

∗ ∗ = 0
− ∗=
/

75
 Synthetic case : a view of kriging outputs
For one block, a distribution is re-estimated using Kriging of 100 realizations of a
random function. Two drillholes sampling are tested:
Dense sampling
Loose sampling

Real histogram Estimated histogram Cross-plot

76
Conditional Unbiasedness
In addition to the global unbiasedness, a good interpolator should
be conditionally unbiased.

[
E Z0 Z = Z *
] *

77
Slope of regression
The linear regression between the true and estimated grades, with the
local condition of search ellipsoïd and variogram model

Slope = 0.95 and Slope = 0.67 and conditionally

conditionally unbiased unbiased ??

78
 Kriging variance
• It depends only on the variogram model and the data configuration

• It is not suited for calculation of confidence intervals, or for calculating

probabilities to exceed a cut-off, unless :
1

0.9
- The grade histogram is Gaussian and 0.8

0.7 Kriging
- The block(s) are well estimated and 0.6 standard

fr e q u e n c y
0.5 deviation
- The stationary assumption is fulfilled 0.4

0.3

0.2

0.1

0
0 10 20 30 40 50 60 70 80 90 100
Possible Z values

79
Weight of the mean
• In Simple Kriging
- The mean ‘m’ is considered to be known (and invariant)
- It is used determine weights
- A specific weight lm is assigned to the mean
- When local estimation is of poor quality the weight assigned to the mean
is high

• In Ordinary Kriging
- Local mean is considered to be unknown (local stationarity only)
- It is determined by the local search ellipsoid
- The weight of the mean can be calculated (by SK) as a test to quantify the
uncertainty of the OK

80
Variance of Z*
Kriging provides an estimate of the dispersion variance of the estimator.
It is a global parameter locally estimated.

Low variance of Z* means the estimate is smoothed (range < sampling distance)
High variance of Z* means the local estimate is good and allows high selectivity

∗ = − −
81
Lagrange Parameter

It is highly sensitive to the number of data and their configuration:

5 samples 20 samples 200 samples

is between 0 and -1 .

It can influence other parameters greatly especially if the

estimation is performed with few samples: Var Z* , Cov(Z,Z*)

82
Dispersion variance
l The dispersion variance (Var Z*) describes the smoothing
effect of Kriging:
- High variance : no smoothing effect, the estimate can potentially take
values between high and low grades.
- Low variance : smoothing effect, the estimate will be close to the average
grade.
∗ = − −

Estimate Z* Var Z*
Page 83
Smoothing effect of kriging
How is the variogram calculated on kriged values compared to the variogram of
the data?

The variogram computed on the kriging grid has a sill that is much lower than the
variogram computed on the real data. (2 times lower)
The range is highest on the variogram computed on the kriging estimation.

84
Kriging efficiency
If the KE is negative, it means that the quality of the average
estimation, characterized by is too poor, and that there is a clear
conditional bias.
∗ ∑
− + +
= = =

0.75

0.50

V KE
0.25
∗
Good estimates : ≈ 0.00

Bad estimates : ≈
-0.25

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

V Slope

85
KRIGING
NEIGHBORHOOD
Neighborhood
2 types of neighborhood can be used:
Unique:
- all data are used for each target estimate
- applicable only for small datasets

Moving:
Each time a new point/block needs to be kriged, a new set of
limited number of data receiving weights will be performed.

87
Unique Neighborhood

• All data used for each target estimate

• Kriging matrix is inverted once for all targets

• Provides maximum data for estimation

• Prevents artifacts from abrupt changes

• Impractical for many minerals projects due to memory

limitations and CPU time- kriging matrix is unworkable

88
Moving Neighborhood
• Compulsory for large datasets or when stationarity not fulfilled
at any scale

• Data outside a 2D ellipse or 3D ellipsoid is discarded

• Selection criteria (defined by user) define a subset of the data

inside the ellipse for estimation

• Each target will be kriged with a different set of data

• How to optimize?

89
Moving Neighborhood
Search ellipse parameters depend
upon:
• Stationarity extension
• Data density
• Variogram parameters
⁻ spatial scale
⁻ anisotropy ratio
⁻ azimuth

Kriging: a Linear Estimator

Z 0 = l1Z1 + l2 Z 2 + ... + ln Z n
Z1 , Z 2 ,..., Z n Data Points
l1 , l2 ,..., ln Weights

90
Moving neighborhood parameters
• Range: ellipsoid radius should be slightly greater than
variogram range
₋ High variance model (high nugget, short range) and/or scarce data: large
ellipsoid is advised
₋ High continuous model (low nugget, large range) and/or numerous data:
large ellipsoid is not required.
₋ Definition of ellipsoid size and shape should be consistent with the
stationary hypothesis

• Orientation of ellipsoid is often parallel to anisotropy

but it is not compulsory, it depends on model
orientation and data configuration

91
 Anisotropic distances

When neighborhood
ellipse is parallel to
variogram anisotropy:
Using anisotropic
distances during data
selection phase
increases selection of
data points in
direction of maximum
continuity

Using anisotropic distances Using euclidean distances

92
Moving neighborhood parameters
• Sectors force the search to be spread around the target, to avoid
selecting samples from just one direction or drill hole

• Minimum number of samples: common to need estimate for all

points/blocks, to be balanced against quality

• Maximum number of samples: more data is generally better but

can induce negative estimates

• The combination of range and maximum samples is important-

need to balance need to estimate (all) blocks against quality

93
 How to choose the model?
• The variogram model should not contradict the
naturalist interpretation

• The cross validation can help to choose the model:

- For a given data set the cross validation checks the

consistency of the model of spatial correlation with the data.
- The procedure is the following :
o Each data point is removed in turn
o It is then kriged using the other data.
o Statistics on the errors are then achieved.

94
COKRIGING
Ordinary Cokriging
• Principle: Estimate Z1 using a linear combination of Z1
and Z2 measured at scattered data points
Z *
1 = åa l 1a Z a
1 + åa l a2 Z a
2

• As for the Kriging the Least square method is apply so

as to find the optimal weights on each data:
λ Z2,z Z1,z
Target Z1,0
λ

Z1,x

Z2,y
96
Ordinary Cokriging
• Principle: Estimate Z1 using a linear combination of Z1
and Z2 measured at scattered data points
Z *
1 = åa l 1a Z a
1 + åa l a2 Z a
2

• The ordinary cokriging system is

ì å λ 1β C αβZ 1 + å λ 2β C αβZ 1 Z 2 + μ1 = C Z1
α 0 "a
ï β β
ï
ï åβ åβ 2 αβ
λ 1
β
C Z 1Z 2
αβ + λ β
C Z 2
+ μ2 = C Z 1Z
α 0
2
"a
í
ï åa 1 l a
= 1
ï
ï
î
åa
l a
2 = 0

97
 Ordinary Cokriging
Matrix notation
é l 11 ù é C 10Z 1 ù
é C Z1
L C Z1
C Z 1Z 2
L C Z 1Z 2
1 0ù ê ú ê ú
ê
11 1n 11 1n
ú ê M ú ê M ú
ê M M M M M Mú ê n ú ê C nZ 01 ú
l2
ê C nZ11 L C nnZ 1 C Z 1Z
L C nnZ 1 Z 2 1 0 ú ê n ú ê Z 1Z 2 ú
ú êl2
2

ê Z 1Z 2 ú ê C 10 ú
n1

ê C 11 L C 1Zn1 Z 2 C Z 2
L C 1Zn 2 0 1ú ê ú = ê ú
11
´ M M
ê M M M M M M ú ê ú ê Z 1Z 2 ú
ú ê l1
n
ê Z 1Z 2 Z 1Z Z 2 Z 2 ú êC n0 ú
ê C n1 L C nn
2
C n1 L C nn 0 1ú êm ú ê 1 ú
1
ê 1 L 1 0 L 0 0 0ú ê ú ê ú
ê ú êm 2 ú ê 0 ú
ë 0 L 0 1 L 1 0 0û ê ú
ê ú
ë û ë û
Covariances between samples and variables Unknown Covariances
weights between the
data and the
n ´ N equations + 2 target
10 samples, 10 variables Þ 102 equations regarding the
variable of
interest

98
No spatial correlation between
variables: self krigeable
• When variables Z1 and Z2 are not at all correlated,
- Correlation coefficient is 0
- Cross variogram is flat and equal to 0
• Therefore Z1 is self-krigeable
• There is no benefit in cokriging Z1 and Z2
Cross variogram

γZ1Z 2 (h) = 0 "h

Perfect correlation envelope
(negative or positive)

99
Perfect linear correlation: self krigeable
• When variables Z1 and Z2 are
very highly correlated
Z2 = a ´ Z1 + b

• The cross variogram is equal to

the perfect correlation
envelope (negative or positive)

Cross variogram
• Information is redundant à
Using Z2 in a cokriging of Z1 does
not improve the estimation of
Z1

b iju £ ± b uii b u
jj

" i, j
100
Simple and cross variograms
proportional: self krigeable
Variable Z1 is self-krigeable when Z1 and
Z2 are intrinsically correlated
Cross variogram
• The simple and cross variograms are
proportional
Proportional variograms differ only in their sills, the
ranges are equal

• Isotopic case: there is no benefit in cokriging

Z1 and Z2

• Heterotopic case: some benefit of cokriging

101
Example
Cokriging of two variables:
- U: high statistical variance – low continuity
- V: low statistical variance – high continuity

102
Weights
Weights on Main variable (U) Weights on Secondary variable (V)

U less continuous than V: positive weights in center

103
Weights
Weights on Main variable (V) Weights on Secondary variable (U)

V more continuous than U: negative weights in center

104
Weights
Weights on Main variable Weights on Secondary variable

Same continuity : intrinsic correlation, variables independant

105
Weights
Secondary variable Weight profile

+
Main
Sec. -
Less continuous than Main « Smooth »

+
Sec. Main
-
More continuous than Main « Wrinkle »

+
Sec. Main
-
Same continuity as Main No transformation

106
Weights
On the main variable:
Weights are independent of the statistical variance.
On the secondary variable:
Weights are proportional to the ratio of both statistical variances
(V1/V2) to correct for the variability around the secondary variable.
NB: Extreme difference in variance may lead to numerical
instabilities -> need to work with normalized variables

Example:
∗
= × 1.1 + × 0.8 + × 11401 + × 21971
, ,
= × 1.1 + × 0.8 + × × 11401 + × × 21971

107
Example
Two variables:
- U (variable to be estimated)
- V (secondary variable – always defined on 470 locations)
Two datasets:
Isotopic Heterotopic
(U defined for 470 data) (U defined for 282 data)

108
Example: Estimates vs. Reality
Kriging Cokriging
Heterotopic
Isotopic

109
Example: Quality of estimation
Kriging Cokriging
Heterotopic
Isotopic

110
Example: Quality of estimation
Kriging Cokriging
Heterotopic
Isotopic

111
Example: Correlations
Kriging Cokriging

Heterotopic
Reality
Isotopic

112
Useful features of Cokriging

• Cokriging is consistent
- Cokriging the sum is same as sum of cokriging
{Z1 ( x) + Z2 ( x)}CK = Z1 ( x)CK + Z2 ( x)CK

• The Cokriging system can be simplified

- (3) variables: Z 2 - Z1 = Z
one of these variables is redundant
- (3-1) variables are estimated

113
Kriging or cokriging ?

No

114
CROSS-VALIDATION
Cross validation
• A tool for evaluating kriging parameters and inputs
- Neighborhood parameters
- Data configurations

• Checks consistency of the model with data

• Procedure
- A data point is removed from the neighborhood
- The removed point location is kriged using the other data
- All points are removed and estimated in turn, returning a full
set of estimated values with corresponding true values
- Statistics are computed

116
Cross validation
å (Z )
N
Mean errors 1
• merror = i - Z i*
Not globally biased if 0 N i =1

å (Z )
• Variance of the error 1 N
* 2
s ² error = i -Z i
As close to 0 as possible N i =1

Mean standardized error 1 N

(Z - Z *
)
å
•
mstd -error = i i
The standardized
Should be as close to 0 as possible N i =1 s i* error should
have a normal
distribution
• Variance of the standardized
1 N
(Z - Z ) * 2
error s2 std -error
=
N
åi =1
i
s i*2
i

Should be as close to 1 as possible

117
Isatis Cross Validation Test Window
Panel output includes:

• True value

• Estimated value

• Standard deviation

118
 Cross Validation Graphical Outputs

Cross plot of true

Base map and estimated
values
Z i vs. Z i*

Check for
conditional bias

Histogram of
standardized errors Cross plot of

(Z i - Z *
i ) standardized
errors and
s i2* estimated values
(Z - Z )
i
*
i
vs. Z*i
σ i2*

119
Cross validation limitations
• Cross-validation allows to check and compare numerically
neighborhood parameters and model choices

• In 3D, the method is limited by the drill hole configuration:

- Often, only the drill hole direction component of the
variogram is adequately tested
- The cross-hole samples are too distant for testing the cross-
hole variogram component

• The standardized results are to be taken cautiously as they are

dependent on the sill value of the variogram

120