Introduction to Geostatistics
!
"
#
"
$
%
&')(
*+
,.-
/10
23
4
*+
+5
6
7
3
4
8
3
9
:<;
=
>@?
;
A
;
B
CED
F
>HG
>
I
F
>
IKJ
L
;
M
M
V
WYX
Z
Centre de Géostatistique - EMP
Founded 1968 by:
Georges M
ATHERON
Director:
Jean-Paul C
HILÈS
Permanent staff:
14 research scientists
Funding (salaries):
60% by contract research,
,→
this conditions an application driven research
V
WYX
CG: Main application fields
Petroleum exploration, mining
Environmental sciences, climatology
Health: epidemiology
Fisheries, demography . . .
Software products:
Isatis
, Heresim. . .
sold by Geovariances International (
Also:
Bioinformatics group (SVM, kernel methods)
V
WYX
Geostatistics worldwide
Other groups:
Stanford (petroleum), Trondheim (petroleum), Calgary
(mining, petroleum), Brisbane (mining), Johannesburg
(mining), Valencia (hydrogeology),. . .
Main meetings:
International Geostatistics Conference:
1st in Rome (1975), . . . , 7th in Banff (2004)
−→
2008: Santiago de Chile
geoENV (european geostatistics conference for
environmental applications):
1st in Lisbon (1996),. . . , 5th in Neuchatel (2004)
−→
2006: Greece
Software:
R (
−→
V
WYX
Geostatistics
definition
V
WYX
Geostatistics
is an application of
the Theory of
Regionalized Variables
(usually considered as realizations of
Random Functions
)
to geology and mining (fifties)
to natural phenomena in general (seventies)
(re-)integrated mainstream statistics (nineties)
V
WYX
Concepts
Variogram: description of the spatial/temporal correlation of
a phenomenon
Kriging: optimal linear prediction method for estimating
values of a phenomenon at any location of a region
(
→
D. G. K
RIGE
)
Conditional Simulation: stochastic simulation of
realizations, conditional upon the data.
V
WYX
Basic Statistics
concepts
W
V
WYX
Center of mass
Seven weights
w
are hanging on a bar whose own weight is
negligible:
5
6
7
8
weight w
elementary weight v
center of mass
W
V
WYX
Center of mass
The weights
w
are suspended at points:
z = 5, 5.5, 6, 6.5, 7, 7.5, 8,
The mass
w(z)
of the weights is
w(z) = 3, 4, 6, 3, 4, 4, 2.
The location
z
where the bar, when suspended, stays in
equilibrium is:
z =
1
P
k
w(z
k
)
7
X
k
=1
z
k
w(z
k
)
W
V
WYX
Z
Center of mass
Defining normed weights:
p(z
k
) =
w(z
k
)
P
k
w(z
k
)
with
P
k
p(z
k
) = 1
, we can write:
z =
7
X
k
=1
z
k
p(z
k
)
W
V
WYX
Z
Z
Center of mass
The weights
w(z
k
)
are subdivided into
n
elementary weights
v
α
:
5
6
7
8
weight w
elementary weight v
center of mass
with corresponding normed weights
p
α
= 1/n
:
z =
n
X
α
=1
z
α
p
α
=
1
26
26
X
α
=1
z
α
= 6.4
W
V
WYX
Z
Center of mass
The average squared distance to the center of mass
dist
2
(z) =
1
n
n
X
α
=1
(z
α
− z)
2
= .83
gives an indication about the dispersion of the around the center
of mass
z
.
W
V
WYX
Z
Histogram
5
6
7
8
mean m*
The mean value
m
?
of data
z
α
is equivalently,
m
?
=
1
n
n
X
α
=1
z
α
W
V
WYX
Z
Histogram
The average squared deviation from the mean is the variance
s
2
=
1
n
n
X
α
=1
(z
α
− m
?
)
2
Its square-root is called the standard deviation.
The normalized weights
p(z
k
)
are the frequencies of the
occurence of the values
z = 5, 5.5, 6, 6.5, 7, 7.5, 8
.
n
is the number of samples.
W
V
WYX
Z
Cumulative histogram
An alternate way to represent the frequencies of the values
z
is
to cumulate them from left to right:
Z
FREQUENCY
CUMULATIVE
1
2
3
4
5
6
7
8
W
V
WYX
Z
Probability distribution
Suppose we draw randomly values
z
from a set of values
Z
.
We call
Z
a random variable and
z
its realizations,
z ∈ R
.
The mathematical idealization of the cumulative histogram
is the probability distribution function
F (z)
defined as:
F (z) = P (Z < z)
The probability
P (Z < z)
indicates the theoretical frequency
of drawing a realization lower than a given value
z
.
W
V
WYX
Z
Probability density
We shall only consider differentiable distribution functions.
The derivative of the probability distribution function is the
probability density
p(z)
:
F (dz) = p(z) dz
Properties:
0 ≤ p(z) ≤ 1
Z
p(z) dz = 1
W
V
WYX
Z
Expected value
The idealization of the concept of mean value is the
mathematical expectation:
E
Z
=
Z
z ∈ R
z p(z) dz = m.
The expectation is a linear operator.
Let
a
,
b
be constants:
E
a
= a,
E
b Z
= b E
Z
,
so that
E
a + b Z
= a + b E
Z
W
V
WYX
Z
Variance
The second moment of the random variable
Z
is:
E
Z
2
=
Z
z ∈ R
z
2
p(z) dz
The variance
σ
2
is defined as:
var(Z) = E
h
Z − E
Z
2
i
= E
(Z − m)
2
= σ
2
Alternate expression:
multiplying out we get
var(Z) = E
Z
2
+ m
2
− 2mZ
and, as the expectation is a linear operator,
var(Z) = E
Z
2
− E
Z
2
W
V
WYX
Covariance
Covariance
σ
ij
between
Z
i
and
Z
j
:
cov(Z
i
, Z
j
) = E
Z
i
− E
Z
i
· Z
j
− E
Z
j
= E
(Z
i
− m
i
) · (Z
j
− m
j
)
= σ
ij
where
m
i
and
m
j
are the means of the random variables.
Covariance of
Z
i
with itself:
σ
ii
= E
(Z
i
− m
i
)
2
= σ
2
i
Correlation coefficient:
ρ
ij
=
σ
ij
q
σ
2
i
σ
2
j
W
V
WYX
Z
Linear regression
V
WYX
Regression line
●
● ●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
1
z
2
z
2
z
1
z * = a
+ b
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
1
2
m*
m*
V
WYX
Optimal regression line
Two variables with experimental covariance:
s
12
=
1
n
n
X
α
=1
(z
α
1
− m
?
1
) · (z
α
2
− m
?
2
)
The regression line is:
z
?
1
= a z
2
+ b
with slope
a
and intercept
b
.
Minimizing the quadratic distance:
dist
2
(a, b)
=
1
n
n
X
α
=1
(z
α
1
− a z
α
2
− b)
2
we get
a =
s
12
s
2
2
b = m
?
1
− a m
?
2
V
WYX
Optimal regression line
z
?
1
=
s
12
s
2
2
(z
2
− m
?
2
) + m
?
1
= m
?
1
+
s
1
s
2
r
12
(z
2
− m
?
2
)
At the minimum the squared distance is:
dist
2
min
(a, b) = s
2
1
1 − (r
12
)
2
V
WYX
Multiple linear regression
W
V
WYX
Multivariate data set
The data matrix
Z
with
n
samples of
N
variables:
Samples
V ariables
z
11
. . . z
1i
. . . z
1N
..
.
..
.
..
.
z
α
1
. . . z
αi
. . . z
αN
..
.
..
.
..
.
z
n
1
. . . z
ni
. . . z
nN
W
V
WYX
Matrix of means
Define a matrix
M
with the same dimension
n × N
as
Z
,
replicating
n
times in its columns the mean value of each
variable:
M
=
m
?
1
. . . m
?
i
. . . m
?
N
..
.
..
.
..
.
m
?
1
. . . m
?
i
. . . m
?
N
..
.
..
.
..
.
m
?
1
. . . m
?
i
. . . m
?
N
W
V
WYX
Centered variables
A matrix
Z
c
of centered variables is obtained by subtracting
M
from the raw data matrix:
Z
c
=
Z
−
M
W
V
WYX
Variance-covariance matrix
The matrix
V
of experimental variances and covariances is:
V
=
1
n
Z
>
c
Z
c
=
var(z
1
)
. . . cov(z
1
, z
j
) . . . cov(z
1
, z
N
)
..
.
. ..
..
.
cov(z
i
, z
1
)
. . .
var(z
i
)
. . . cov(z
i
, z
N
)
..
.
. ..
..
.
cov(z
N
, z
1
) . . . cov(z
N
, z
j
) . . .
var(z
N
)
=
s
11
. . . s
1j
. . . s
1N
..
.
. ..
..
.
s
i
1
. . .
s
ii
. . .
s
iN
..
.
. ..
..
.
s
N
1
. . . s
N j
. . . s
N N
W
V
WYX
Multiple linear regression
For a regression of
z
0
on the
N
variables from
n
samples
we have the matrix equation
z
?
0
= m
0
+ (Z − M) a
The squared distance between
z
0
and the hyperplane is:
dist
2
(a) =
1
n
(z
0
− z
?
0
)
>
(z
0
− z
?
0
)
= var(z
0
) + a
>
Va
− 2 a
>
v
0
,
where
v
0
is the vector of covariances
between
z
0
and
z
i
, i = 1, . . . , N
.
W
V
WYX
Z
Minimizing the squared distance
The minimum is found for:
∂ dist
2
(a)
∂a
= 0
⇐⇒
2 Va − 2 v
0
= 0
⇐⇒
Va
= v
0
This system of linear equations:
var(z
1
)
. . . cov(z
1
, z
N
)
..
.
. ..
..
.
cov(z
N
, z
1
) . . .
var(z
N
)
a
1
..
.
a
N
=
cov(z
0
, z
1
)
..
.
cov(z
0
, z
N
)
has exactly one solution,
if the determinant of
V
is different from zero.
The squared distance at the minimum is:
dist
2
min
(a) = var(z
0
) − a
>
v
0
W
V
WYX
Simple kriging
W
V
WYX
Spatial data
Data points
x
α
and the estimation point
x
0
in a spatial domain
D
0
x
x
α
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
❍
D
W
V
WYX
Translation invariance
The expectation and the covariance are both assumed
translation invariant over the domain,
i.e. for any vector
h
between points
x
and
x
+h
:
E
Z(x+h)
= E
Z(x)
= m
cov
Z(x+h), Z(x)
= C(h)
The expectation
E
Z(x)
has the same value
m
at any point
x
of the domain
D
.
The covariance between any pair of locations
depends only on the vector
h
.
W
V
WYX
Known mean
We assume the mean
m
is known
and build the estimator:
Z
?
(x
0
) = m +
n
X
α
=1
w
α
Z(x
α
) − m
i.e.
Z
?
(x
0
) − m =
n
X
α
=1
w
α
Z(x
α
) − m
which is implicitly without bias:
E
h
Z
?
(x
0
) − m
i
=
n
X
α
=1
w
α
E
h
Z(x
α
) − m
i
= 0
W
V
WYX
Simple kriging equations
The kriging equations with known mean are simple:
n
X
β
=1
w
SK
β
C(x
α
−x
β
) = C(x
α
−x
0
)
for
α= 1, . . . , n
i.e.
the linear combination of weights with
the covariances between a data point
and the other data points
=
the covariance between that data point
and the point to estimate.
The variance of the Simple Kriging estimate is:
σ
2
SK
= σ
2
−
n
X
α
=1
w
SK
α
C(x
α
−x
0
)
W
V
WYX
Simple kriging: a multiple linear regression
Simple kriging is a multiple linear regression between spatial
random variables.
Like:
Va
= v
0
,
we have:
Cw
= c
0
Writing out the equation system:
var(Z(x
1
))
. . . cov(Z(x
1
), Z(x
N
))
..
.
. ..
..
.
cov(Z(x
N
), Z(x
1
)) . . .
var(Z(x
N
))
w
1
..
.
w
N
=
cov(Z(x
0
), cov(Z(x
1
))
..
.
cov(Z(x
0
), Z(x
N
))
W
V
WYX
Regionalized variables
and random function
V
WYX
The concept of a Random Function
Consider a domain
D
with points
x
:
●
x
D
Let
Z(x)
be a random variable at a location
x
∈ D
.
The family of random variables
n
Z(x); x ∈ D
o
is called a Random Function.
V
WYX
Regionalized Variable
The regionalized variable
z(x)
is the spatial variable of interest
(“reality”).
Data does not generally allow a deterministic reconstruction of
the regionalized variable.
The regionalized variable
z(x)
is considered as a realization
(draw) of a random function
Z(x)
.
For a given data set, different realizations containing the data
are equally plausible to represent the regionalized variable.
V
WYX
Z
Epistemological Problem
We possess data about only one realization:
how can we specify the random function?
Objective quantities
that describe the regionalized variable and
conventional parameters
that are constitutive of the model have
to be distinguished.
The quantities are estimated from data,
but the parameters are chosen.
−→
G
“Estimating and Choosing”
V
WYX
Variogram
definition
V
WYX
The Variogram
The vector
x
=
x
1
x
2
: coordinates of a point in 2D.
Let
h
be the vector separating two points:
●
●
h
D
x
x
β
α
We compare sample values
z
at a pair of points with:
z(x + h) − z(x)
2
2
V
WYX
The Variogram Cloud
Variogram values are plotted against distance in space:
●●
●● ●
●
●
●●
●
●
●
●
●
●
● ● ●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
2
2
●
●
|h|
| z(t+h) - z(t) |
V
WYX
The Experimental Variogram
Averages within distance (and angle) classes
k
are computed:
●●
●● ●
●
●
●●
●
●
●
●
●
●
● ● ●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
1
2
3
4
5
6
7
8
9
h
h
h h h h h h
h
|h|
γ∗
(h )
k
V
WYX
The Theoretical Variogram
A theoretical model is fitted:
γ
(h)
|h|
V
WYX
Intrinsic Hypothesis
The first two moments of the
increments
are assumed stationary
(translation-invariant):
the expectation does not depend on
x
E
h
Z(x+h) − Z(x)
i
= 0
the variance depends only on
h
var
h
Z(x+h)−Z(x)
i
= 2 γ(h)
This type of stationarity is called
intrinsic
.
,→
The stationarity of the increments does not imply the
stationarity of
Z
.
V
WYX
Definition of the Variogram
By the intrinsic hypothesis:
γ(h) =
1
2
E
h
Z(x+h) − Z(x)
2
i
Properties
- zero at the origin
γ(0) = 0
- positive values
γ(h) ≥ 0
- even function
γ(h) = γ(−h)
Regionalized variable
Behavior at the origin
←→
continuous and differentiable
←→
not differentiable
←→
discontinuous
V
WYX
Variogram and Covariance Function
The covariance function is defined as:
C(h) = E
h
Z(x) − m
·
Z(x+h) − m
i
where
stationarity of the first two moments
of
Z
is assumed.
A variogram can be constructed from any covariance function:
γ(h) = C(0) − C(h)
Conversely, however, only if the variogram is
bounded
does a
corresponding covariance function
C(h)
exist.
The variogram characterizes a larger class of random functions.
This is why it is preferred in geostatistics.
V
WYX
Variogram
examples
V
WYX
Z
Power variogram
γ(h) = |h|
p
,
0 < p ≤ 2
Power model
DISTANCE
VARIOGRAM
-4
-2
0
2
4
0
1
2
3
4
5
p=1.5
p=1
p=0.5
V
WYX
Spherical covariance function
C(h) =
3
2
|h|
a
−
1
2
|h|
3
a
3
|h|≤a
V
WYX
Exponential covariance function
C(h) = exp
−
|h|
a
V
WYX
Gaussian covariance function
C(h) = exp
−
|h|
2
a
2
V
WYX
Cardinal sine covariance function
C(h) =
sin
|h|
a
|h|
a
V
WYX
Geometric anisotropy
of the variogram
W
V
WYX
Geometric anisotropy
In practice the range of the variogram may change depending
on the direction:
h
1
h
2
1
h’
2
h’
Correction:
rotation
h
0
= Qh
of angle
θ
where
Q
=
cos θ sin θ
− sin θ cos θ
linear transformation of the coordinates
h
0
= (h
0
1
, h
0
2
)
W
V
WYX
Rotation in 3D
In 3D the rotation is obtained by a composition of elementary
rotations:
Q
=
cos θ
3
sin θ
3
0
− sin θ
3
cos θ
3
0
0
0
1
1
0
0
cos θ
2
sin θ
2
0
− sin θ
2
cos θ
2
0
cos θ
1
sin θ
1
0
− sin θ
1
cos θ
1
0
0
0
1
where
θ
1
,
θ
2
,
θ
3
are Euler’s angles.
W
V
WYX
2D example: Ebro river vertical section
-10.
-5.
0.
Ebro river (Kilometer)
-6.
-5.
-4.
-3.
-2.
-1.
0.
Depth (Meter)
185 Hydrolab Surveyor III conductivity measurements
W
V
WYX
2D conductivity variogram model
D1
D2
M1
M2
0.
1.
2.
3.
4.
5.
6.
Distance (D1: km; D2: m)
0.
250.
500.
750.
1000.
1250.
Variogram : CONDUCTIVITY
Experimental variogram for D1=horizontal, D2=vertical.
Anisotropic cubic variogram model in both directions (M1, M2).
Abscissa scale: kilometers for D1 and meters for D2.
W
V
WYX
Z
Behavior at the origin
of the variogram
V
WYX
Ebro river:
water samples
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
Distance from mouth (km)
-6.
-5.
-4.
-3.
-2.
-1.
0.
Depth (m)
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
-6.
-5.
-4.
-3.
-2.
-1.
0.
Depth (m)
47 water samples (top)
185 conductivity values (bottom)
V
WYX
Nitrate variogram: which behavior at origin?
D1
D2
M1
M2
CUBIC
0.
1.
2.
3.
Lag (RED: m; BLACK: km)
0.
1000.
2000.
3000.
Variogram: NITRATE
D1
D2
M1
M2
EXPONENTIAL
0.
1.
2.
3.
Lag (RED: m; BLACK: km)
0.
1000.
2000.
3000.
Variogram: NITRATE
Nitrate experimental variogram with two alternate models.
V
WYX
Cubic variogram: conditional simulations
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
-4.
-3.
-2.
-1.
0.
Depth (m)
>=124.8
117
109.2
101.4
93.6
85.8
78
70.2
62.4
54.6
46.8
39
31.2
23.4
15.6
7.8
<0
M
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
Distance from mouth (km)
-4.
-3.
-2.
-1.
0.
Depth (m)
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
-4.
-3.
-2.
-1.
0.
Depth (m)
V
WYX
Exponential model: conditional simulations
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
-4.
-3.
-2.
-1.
0.
Depth (m)
>=124.8
117
109.2
101.4
93.6
85.8
78
70.2
62.4
54.6
46.8
39
31.2
23.4
15.6
7.8
<0
M
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
-4.
-3.
-2.
-1.
0.
Depth (m)
-15.0
-12.5
-10.0
-7.5
-5.0
-2.5
Distance from mouth (km)
-4.
-3.
-2.
-1.
0.
Depth (m)
V
WYX
Kriging of the mean
of a random function
V
WYX
Spatially Correlated Data
Sample locations
x
α
in a geographical domain:
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
With spatial correlation we need to consider that:
sample points have a different number of immediate
neighbors,
distances to neighboring points play a role.
How should samples be weighted in an optimal way?
V
WYX
Estimation of the Mean Value
Using the formula of the arithmetic mean:
M
?
=
1
n
n
X
α
=1
Z(x
α
)
all samples get the same weight:
1
n
We rather need an estimator:
M
?
=
n
X
α
=1
w
α
Z(x
α
)
with weights
w
α
reflecting the spatial correlation.
V
WYX
Stationary random function
We assume translation-invariance of mean and covariance:
∀ x ∈ D :
E
h
Z(x)
i
= m;
∀ x
α
, x
β
∈ D :
C(x
α
, x
β
) = C(x
α
−x
β
).
The estimation error in our statistical model:
M
?
|
{z
}
estimated value
−
m
| {z }
true value
should be zero, on average:
E
h
M
?
− m
i
= 0
V
WYX
No bias
No bias is obtained using weights of unit sum:
n
X
α
=1
w
α
= 1
Consider:
E
h
M
?
− m
i
= E
h
n
X
α
=1
w
α
Z(x
α
) − m
i
=
n
X
α
=1
w
α
E
h
Z(x
α
)
i
|
{z
}
m
−m
= m
n
X
α
=1
w
α
| {z }
1
−m
=
0
V
WYX
Z
Variance of the estimation error
The variance
σ
2
E
of the estimation error is:
var(M
?
− m) = E
h
(M
?
− m)
2
i
−
E
h
M
?
− m
i
2
|
{z
}
0
= E
h
M
?2
− 2 mM
?
+ m
2
i
=
n
X
α
=1
n
X
β
=1
w
α
w
β
E
h
Z(x
α
) Z(x
β
)
i
−2 m
n
X
α
=1
w
α
E
h
Z(x
α
)
i
|
{z
}
m
+m
2
⇒
σ
2
E
=
n
X
α
=1
n
X
β
=1
w
α
w
β
C(x
α
− x
β
)
V
WYX
Minimal estimation variance
We want weights
w
α
that produce a minimal estimation variance:
minimum of
var(M
?
− m)
subject to
n
X
α
=1
w
α
= 1
The objective function
ϕ
has
n+1
parameters:
ϕ(w
1
, . . . , w
n
, µ) = var(M
?
− m) − 2 µ
n
X
α
=1
w
α
− 1
with
µ
a Lagrange multiplier. Setting partial derivatives to zero:
∀α :
∂ϕ(w
1
, . . . , w
n
, µ)
∂w
α
= 0,
∂ϕ(w
1
, . . . , w
n
, µ)
∂µ
= 0
V
WYX
Kriging equations
The method of Lagrange yields the equations for
the optimal weights
w
KM
α
of the kriging of the mean:
n
X
β
=1
w
KM
β
C(x
α
− x
β
) − µ
KM
= 0
for
α = 1, . . . , n
n
X
β
=1
w
KM
β
= 1
The variance at the minimum:
σ
2
KM
= µ
KM
is equal to the Lagrange multiplier.
V
WYX
Case of no autocorrelation
When the covariance model is a pure nugget-effect:
C(x
α
− x
β
) =
σ
2
if
x
α
= x
β
0
if
x
α
6= x
β
the kriging of the mean system simplifies to:
w
KM
α
σ
2
= µ
KM
for
α = 1, . . . , n
n
X
β
=1
w
KM
β
=
1
The solution weights are all equal:
w
KM
α
=
1
n
⇒
M
?
=
1
n
n
X
α
=1
Z(x
α
)
the arithmetic mean!
µ
KM
= σ
2
KM
=
1
n
σ
2
V
WYX
Ordinary Kriging
at a point in the domain
V
WYX
Estimation at a Point
Sample locations
x
α
(dots)
in a domain
D
:
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
❍
x
0
We wish to estimate a value
Z
?
at a point
x
0
.
V
WYX
Ordinary kriging
The estimate
Z
?
is a weighted average of data values
Z(x
α
)
:
Z
?
(x
0
) =
n
X
α
=1
w
α
Z(x
α
)
with
n
X
α
=1
w
α
= 1
The weights
w
OK
α
of the Best Linear Unbiased Estimator (BLUE)
are solution of the system:
n
X
β
=1
w
OK
β
γ(x
α
−x
β
) + µ
OK
=
γ(x
α
−x
0
)
∀α
n
X
β
=1
w
OK
β
=
1
Minimal variance:
σ
2
OK
= µ
OK
+
n
X
α
=1
w
OK
α
γ(x
α
−x
0
)
V
WYX
Cross-validation
leaving one out and reestimating it
V
WYX
Cross-validation
Comment:
the sound way to cross-validate is to leave out half
of the data locations and to re-estimate them from the other
half : this requires many data! For that reason it is often done in
the following way (implemented in sotware packages). . .
A data value
Z(x
α
)
is left out and a value
Z
?
(x
[α]
)
is estimated at
location
x
α
by ordinary kriging.
The notation
[α]
means that the sample at
x
α
has not been used
for estimating
Z
?
(x
[α]
)
.
The difference between the data value and the estimated value:
Z(x
α
) − Z
?
(x
[α]
)
gives an indication of how well the data value fits into the
neigborhood of the surrounding data values.
V
WYX
Average cross-Validation error
If the average of the cross-validation errors is not far from zero:
1
n
n
X
α
=1
Z(x
α
) − Z
?
(x
[α]
)
!
∼
= 0
then there is no systematic bias.
A negative (positive) average error represents
systematic overestimation (underestimation).
V
WYX
Z
Standardized cross-validation error
The kriging standard deviation
σ
K
represents the error predicted
by the model.
Dividing the cross-validation error by
σ
K
allows to compare the
magnitudes of both errors:
Z(x
α
) − Z
?
(x
[α]
)
σ
Kα
V
WYX
Average squared Standardized Errors
If the average of the squared standardized cross-validation
errors is not far from one:
1
n
n
X
α
=1
Z(x
α
) − Z
?
(x
[α]
)
!
2
σ
2
Kα
∼
= 1
then the actual estimation error is equal on average to the error
predicted by the model.
This quantity gives an idea of the adequacy of the model and of
its parameters.
V
WYX
Mapping with kriging
on a regular grid
with irregularly spaced data
W
W
V
WYX
Kriging for interpolation
Kriging is an estimation method.
It is not the quickest method to make an interpolation on a
regular grid for generating a map.
Its advantages are:
Kriging integrates the knowledge
gained from analysing the spatial structure:
the variogram.
Kriging interpolates exactly: when a sample value is
available at the location-of-interest, the kriging solution is
equal to that value.
Kriging provides an indication of the estimation error: the
kriging variance.
W
W
V
WYX
Generating a map
A regular grid is defined by the computer and
at each node of this grid a value is kriged.
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
●
x
0
Afterwards a graphical representation of this grid is performed,
as a raster of colour squares, as an isoline map, as a bloc
diagram. . .
W
W
V
WYX
Moving Neighborhood
If all data are used:
this is called a unique neighborhood.
Using a subset of close data points:
a moving neighborhood.
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
❍
●
x
0
To choose the size of the moving neighborhood, the range of the
variogram can give an indication.
W
W
V
WYX
Kriging weights
The shape of the kriging weights
V
WYX
Kriging weights
Nugget-effect model
σ
2
OK
= 1.25
●
●
●
❍
25%
●
25%
25%
25%
L
V
WYX
Isotropic variogram
Spherical model with range a/L = 2
σ
2
OK
= .84
●
●
●
❍
●
9.4%
9.4%
40.6%
40.6%
L
Gaussian model with range a/L = 1.5
σ
2
OK
= .30
●
●
●
❍
●
49.8
49.8%
0.2%
0.2%
L
V
WYX
Geometric anisotropy
Spherical with isotropic range
●
●
❍
25%
25%
●
●
25%
25%
L
Spherical with horizontal a/L = 1.5 and vertical a/L = .75
●
●
❍
●
●
L
17.6%
17.6%
32.4%
32.4%
V
WYX
Z
Relative position of samples
σ
2
OK
= .45
σ
2
OK
= .48
❍
●
●
●
33.3%
33.3%
33.3%
❍
●
●
25.9%
37.1%
●
37.1%
The left configuration gives a more reliable estimate.
V
WYX
The screen effect
Spherical model with range a/L = 2
σ
2
OK
= 1.14
34.4%
●
A
65.6%
●
B
❍
σ
2
OK
= 0.87
●
A
●
B
●
C
49.1%
48.2%
2.7%
❍
Adding the sample C screens off the sample B.
V
WYX
Nested variogram
and corresponding linear model
of the random function
V
WYX
Nested Variogram Model
A nested variogram
γ(h)
is composed of
a sum of elementary variograms
γ
u
(h)
with
u = 0, . . . , S
:
γ(h) = γ
0
(h) + . . . + γ
S
(h) =
S
X
u
=0
γ
u
(h)
Each variogram
γ
u
(h)
is build up with a normed variogram
g
u
(h)
multiplied with a coefficient
b
u
(sill, slope):
γ(h) =
S
X
u
=0
b
u
g
u
(h)
V
WYX
Example: Arsenic in soil (Loire, France)
310.
315.
320.
325.
330.
335.
340.
270.
275.
280.
285.
river
Loire
35
×
25 km
2
region. Dots are proportional to sample value.
V
WYX
Example: Nested Variogram Model
A nugget-effect (nug) and two spherical (sph) structures:
γ(h) =
b
0
nug
(h) +
b
1
sph
(h,
a
1
) +
b
2
sph
(h,
a
2
)
0.
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
0.0
0.5
1.0
nugget
long range
short range
(h)
γ
h
V
WYX
Nested Covariance Function
C(h) =
S
X
u
=0
C
u
(h) =
S
X
u
=0
b
u
ρ
u
(h)
where
ρ
u
(h)
are correlation functions.
The
ρ
u
(h)
characterize the spatial correlation
at
different scales
of index
u
.
The coefficents
b
u
represent a decomposition of
the
total variance
σ
2
into variances at different spatial scales:
C(0) = σ
2
=
S
X
u
=0
b
u
V
WYX
Regionalization Model
Z(x)
built up with
uncorrelated components
Y
u
(x)
of zero mean, with covariance functions
C
u
(h)
.
Example:
Z(x) = Y
1
(x) + Y
2
(x) + m
with
Y
1
⊥Y
2
The covariance function of
Z(x)
is nested:
C(h) = C
1
(h) + C
2
(h)
V
WYX
Linear Model with
S + 1
components
Z(x) =
S
X
u
=0
Y
u
(x)
+ m
with
Y
u
⊥Y
v
for
u 6= v
Corresponding
nested
covariance model:
C(h) =
S
X
u
=0
C
u
(h)
=
S
X
u
=0
b
u
ρ
u
(h)
Can components
Y
u
be
extracted
from samples
Z(x
α
)
?
V
WYX
Z
Kriging Spatial Components
A component
Y
1
(x)
at
x
0
is estimated from
n
data:
Y
?
1
(x
0
) =
n
X
α
=1
w
α
Z(x
α
)
“No bias” with
n
X
α
=1
w
α
= 0
:
this filters the mean
m
Minimizing the “estimation variance”:
n
X
β
=1
w
1
β
C(x
α
−x
β
) − µ
1
=
C
1
(x
α
−x
0
)
for
α= 1, . . . , n
n
X
β
=1
w
β
=
0
V
WYX
Z
Z
Example: Short-range Component of As
310.
315.
320.
325.
330.
335.
340.
270.
275.
280.
285.
V
WYX
Z
Example: Long-range Component of As
310.
315.
320.
325.
330.
335.
340.
270.
275.
280.
285.
V
WYX
Z
Demographic application
fertility data
W
W
W
V
WYX
Z
Demographic application: fertility 1990
5°
4°
3°
2°
1°
0°
1°
2°
3°
4°
5°
6°
7°
8°
43°
44°
45°
46°
47°
48°
49°
50°
51°
250 km
Communes de France
Nb of women per "commune"
Mean annual fertility ’90
0
50
100
150
100
500
5000
10000
25000
50000
5e+05
FERT500 class
Data provided by INSEE (
W
W
W
V
WYX
Z
Variograms: class 100-500 women / commune
D1
D2
D3
D4
0.
100.
200.
300.
Distance (km)
0.
25.
50.
75.
100.
Variogram : FERT500
D1
M1
0.
100.
200.
300.
400.
Distance (km)
0.
10.
20.
30.
40.
50.
60.
70.
80.
90.
100.
110.
Variogram : FERT500
long range
short range
W
W
W
V
WYX
Z
Kriging: short range effect only
D1
M1
0.
100.
200.
300.
400.
Distance (km)
0.
10.
20.
30.
40.
50.
60.
70.
80.
90.
100.
110.
Variogram : FERT500
short range
-100.
0.
100.
200.
300.
400.
500.
600.
700.
800.
UTM (Km)
4700.
4800.
4900.
5000.
5100.
5200.
5300.
5400.
5500.
5600.
UTM (Km)
Fert500 estimation
>=9.96
9.18
8.4
7.62
6.84
6.06
5.28
4.5
3.72
2.94
2.16
1.38
0.6
-0.18
-0.96
-1.74
-2.52
-3.3
-4.08
-4.86
-5.64
-6.42
-7.2
-7.98
-8.76
-9.54
-10.32
-11.1
-11.88
-12.66
-13.44
-14.22
<-15
W
W
W
V
WYX
Z
Kriging
-100.
0.
100.
200.
300.
400.
500.
600.
700.
800.
UTM (Km)
4700.
4800.
4900.
5000.
5100.
5200.
5300.
5400.
5500.
5600.
UTM (Km)
FERT500long_range
>=60
59.2188
58.4375
57.6562
56.875
56.0938
55.3125
54.5312
53.75
52.9688
52.1875
51.4062
50.625
49.8438
49.0625
48.2812
47.5
46.7188
45.9375
45.1562
44.375
43.5938
42.8125
42.0312
41.25
40.4688
39.6875
38.9062
38.125
37.3438
36.5625
35.7812
<35
-100.
0.
100.
200.
300.
400.
500.
600.
700.
800.
UTM (Km)
4700.
4800.
4900.
5000.
5100.
5200.
5300.
5400.
5500.
5600.
UTM (Km)
Fert500 estimation
>=65.2
64.1
63
61.9
60.8
59.7
58.6
57.5
56.4
55.3
54.2
53.1
52
50.9
49.8
48.7
47.6
46.5
45.4
44.3
43.2
42.1
41
39.9
38.8
37.7
36.6
35.5
34.4
33.3
32.2
31.1
<30
Long range effect
Short + long range
W
W
W
V
WYX
Z
Kriging with external drift
V
WYX
Z
Drift
Translation-invariant drift:
polynomials, trigonometric functions
External Drift:
an auxiliary variable known everywhere
V
WYX
Z
Z
External drift method
1.
Auxiliary variable
s(x)
known everywhere in the domain
D
.
2.
The relation to the variable of interest is linear:
E
h
Z(x)
i
= b
0
+ b
1
s(x)
V
WYX
Z
Z
Z
External Drift method
Z
?
(x
0
) =
n
X
α
=1
w
α
Z(x
α
)
⇒
E
h
Z
?
(x
0
)
i
=
n
X
α
=1
w
α
b
0
+ b
1
s(x
α
)
V
WYX
Z
Z
Constraint: no bias
The constraint
n
X
α
=1
w
α
= 1
has the effect that the coefficients
b
0
and
b
1
are filtered out:
E
h
Z
?
(x
0
)
i
= E
h
Z(x
0
)
i
=⇒
b
0
+ b
1
n
X
α
=1
w
α
s(x
α
) = b
0
+ b
1
s(x
0
)
=⇒
n
X
α
=1
w
α
s(x
α
) = s(x
0
)
V
WYX
Z
Z
Interpolation of external drift
This second constraint:
n
X
α
=1
w
α
s(x
α
) = s(x
0
)
generates weights
w
α
which interpolate exactly
s(x)
.
V
WYX
Z
Z
Kriging System with linear and external drift
n
X
β
=1
w
β
C(x
α
−x
β
) + µ
0
+ µ
1
x
1
α
+ µ
2
x
2
α
+ µ
3
s(x
α
) = C(x
α
−x
0
), ∀α
n
X
β
=1
w
β
= 1
n
X
β
=1
w
β
x
1
β
= x
1
0
(longitude)
n
X
β
=1
w
β
x
2
β
= x
2
0
(latitude)
n
X
β
=1
w
β
s(x
β
) = s(x
0
)
(external drift)
V
WYX
Z
Z
Kriging temperature
with elevation as external drift
V
WYX
Z
Z
Temperature Data
temperature conditions the growth of plants
Scotland (without the Shetland and Orkney Islands)
average January temperatures (1961-1980)
146 sites, all below 400 m altitude
Ben Nevis, 1344 m
at 3035 nodes of a regular grid
1272 m
Reference:
Int. J. Clim., 14, 77–91, 1994
V
WYX
Z
Z
January temperature vs latitude / longitude
Scatter diagrams of temperature with latitude and longitude:
there is a systematic decrease from west to east, while there is
not much trend in the north-south direction.
V
WYX
Z
Z
E-W and N-S temperature variograms
The model is fitted in the direction without drift
V
WYX
Z
Z
Kriging mean January temperature
V
WYX
Z
Temperature vs elevation
Coastal influence below 50m, then linear relation.
Temperatures only available below 400m.
V
WYX
Z
Z
Kriging temperature with elevation as drift
V
WYX
Z
Temperature estimates vs elevation
The estimated values above 400m are linearly extrapolated
outside the range of the data !
V
WYX
Z
Estimated external drift coefficient
b
?
1
=
n
X
α
=1
w
α
Z(x
α
)
n
X
β
=1
w
β
C(x
α
−x
β
) + µ
0
+ µ
1
x
1
α
+ µ
2
x
2
α
+ µ
3
s(x
α
) =
0
, ∀α
n
X
β
=1
w
β
=
0
n
X
β
=1
w
β
x
1
β
=
0
(longitude)
n
X
β
=1
w
β
x
2
β
=
0
(latitude)
n
X
β
=1
w
β
s(x
β
) =
1
(external drift)
V
WYX
Z
Estimated external drift coefficient
V
WYX
Z
Conditional simulation
V
WYX
Z
Conditional simulation vs Kriging
V
WYX
Z
Change of support
geostatistical simulation of O
3
V
WYX
Z
CASE STUDY: Geostatistical simulation of O
3
of realizations a lognormal random function
800
×
600 Km
2
1
×
1 Km
2
with a range of 50 Km
V
WYX
Z
Simulation of Ozone: 1
×
1 Km
2
support
0. 100.
200.
300.
400.
500.
600.
700.
Km
0.
100.
200.
300.
400.
500.
Km
O3: 1x1km2
>=96
90
84
78
72
66
60
54
48
42
36
30
24
18
12
6
<0
ug/m3
V
WYX
Z
Simulation of Ozone: 10
×
10 Km
2
support
100.
200.
300.
400.
500.
600.
700.
Km
100.
200.
300.
400.
500.
Km
O3: 10x10km2
>=96
90
84
78
72
66
60
54
48
42
36
30
24
18
12
6
<0
ug/m3
V
WYX
Z
Z
Simulation of Ozone: 20
×
20 Km
2
support
100.
200.
300.
400.
500.
600.
700.
Km
100.
200.
300.
400.
500.
Km
O3: 20x20km2
>=96
90
84
78
72
66
60
54
48
42
36
30
24
18
12
6
<0
ug/m3
V
WYX
Z
Simulation of Ozone
0.
100.
200.
O3 (ug/m3)
0.0
0.1
0.2
0.3
Frequencies
SUPPORT: 1x1 Km2
Nb Samples: 480000
Minimum: 0.0
Maximum: 246.3
Mean: 6.7
Std. Dev.: 10.2
0. 10. 20. 30. 40. 50. 60. 70.
O3 (ug/m3)
0.000
0.025
0.050
0.075
0.100
0.125
Frequencies
SUPPORT: 20x20 Km2
Nb Samples: 1131
Minimum: 0.2
Maximum: 72.0
Mean: 6.7
Std. Dev.: 7.9
Increasing the support
: the means are equal,
but the extremes and the
variance
are reduced
V
WYX
Z
Simulation of Ozone
D1
D2
0.
50.
100.
150.
200.
Distance (Km)
0.
25.
50.
75.
100.
Variogram: O3
SUPPORT
1x1 Km2
D1
D2
0.
50.
100.
150.
200.
Distance (Km)
0.
10.
20.
30.
40.
50.
60.
70.
Variogram: O3
SUPPORT
20x20 Km2
Increasing the support
:
the
range
increases
V
WYX
Z
Simulation of Ozone
40. 50. 60. 70. 80. 90. 100.
110.
120.
O3 cutoff (ug/m3)
0.
50.
100.
150.
Mean O3 over cutoff (ug/m3)
black = 1x1Km
blue = 10x10Km
red = 20x20Km
V
WYX
Z
Simulation of Ozone
40. 50. 60. 70. 80. 90. 100.
110.
120.
O3 cutoff (ug/m3)
0.0
0.5
1.0
1.5
Proportion above cutoff (%)
black = 1x1Km2
blue = 10x10Km2
red = 20x20Km2
V
WYX
Z
Change of support
concept
W
W
V
WYX
Z
TOPIC: The Support of a Random Function
3D
Volumes
2D
v
V
t
T
1D
Surfaces
s
S
∆
Soil pollution
Industrial hygienics
Time intervals
Mining
W
W
V
WYX
Z
The Effect of Changing the Support
Distribution of samples on small volumes (cm
3
) is different from
that of model output averages over large blocks (m
3
):
Z
mean
frequency
blocs
samples
The mean of both distributions is the same,
the distribution of the block values is narrower.
W
W
V
WYX
Z
Neglecting the Support Effect
We are often interested in what is above a threshold:
threshold
overestimation!
Neglecting the
support effect
may lead to a systematic
over-estimation. . .
W
W
V
WYX
Z
Neglecting the Support Effect
. . . or to systematic under-estimation:
threshold
underestimation!
⇒
A good estimation method should incorporate
a
change of support
model.
W
W
V
WYX
Z
Z
Kriging of a Block average
(centered at a point in the domain)
V
WYX
Z
Estimation of a block value
Sample locations
x
α
(dots)
in a domain
D
:
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
V
0
We wish to estimate the spatial average
Z
?
for a block
V
0
.
V
WYX
Z
Block Kriging
The block value
Z
?
(V
0
)
is estimated as a weighted average of
the data values
Z(x
α
)
:
Z
?
(V
0
) =
n
X
α
=1
w
α
Z(x
α
)
with
n
X
α
=1
w
α
= 1
The optimal weights
w
OK
α
are obtained from the sytem:
n
X
β
=1
w
OK
β
γ(x
α
−x
β
) + µ
OK
=
γ(V
0
, x
α
)
∀α
n
X
β
=1
w
OK
β
=
1
Kriging variance:
σ
2
OK
= µ
OK
− γ(V
0
, V
0
) +
n
X
α
=1
w
OK
α
γ(V
0
, x
α
)
V
WYX
Z
Block kriging with non-point data
In applications the data can be averaged on blocks
V
α
.
We then use average variograms between these blocks:
γ
V
α
, V
β
=
1
|V
α
| |V
β
|
Z
x
∈V
α
Z
y
∈V
β
γ(x − y) dx dy
This requires the knowledge of the point variogram.
V
0
v
α
V
WYX
Z
Change of support
risk of exceeding ozone alert level
W
W
V
WYX
Z
Change of support
The variability of spatial or temporal data depends on the
averaging volume/interval(= the
support
)
Increasing support, the variability decreases
(reduction of variance, extremes...)
Observations are on
point
support as compared to the
cells
of a numerical model.
End-users are often interested by a support of different
(intermediate) size
−→
blocks
It is thus necessary to describe statistically how variability
changes as a function of support.
If the distribution is monomodal and not too asymmetrical,
an affine correction may suffice. Otherwise, non-linear
geostatistics or geostatistical simulation are needed
Applications:
data aggregation, estimation of small block
statistics, downscaling. . .
W
W
V
WYX
Z
Ozone in Paris on 17 july 1999 at 15h UTC
P6
P7
P13
P18
Garches
Gennevilliers
Neuilly
Aubervilliers
Tremblay
Vitry
Melun
Mantes
Montgeron
RUR_SE
RUR_E
RUR_SO
RUR_O
RUR_NO
RUR_NE
1.4°
1.6°
1.8°
2°
2.2°
2.4°
2.6°
2.8°
3°
3.2°
48.2°
48.4°
48.6°
48.8°
49°
49.2°
49.4°
50 km
Airparif stations and Chimere grid
19 Airparif stations;
25 × 25
grid with
cells of size 6
×
6 km
2
W
W
V
WYX
Z
Air quality regulations
Two ozone thresholds refering to a support of
1 hour
:
−→
Swiss alert level:
120
µ
g/m
3
−→
European alert level:
180
µ
g/m
3
Time support is always specified, yet regulations do not contain
any indication about the
spatial support
!
Suppose the air quality experts agree on the following
spatial decision support
:
a block of
1 × 1
km
2
size
(instead of the CHIMERE
6 × 6
km
2
cell).
We need to model the
point-block-cell
change of support.
W
W
V
WYX
Z
Discrete Gaussian point-block model
(due to Georges M
ATHERON
, 1976)
x
is a point randomly located in a block
v
.
E
Z(x)
| Z(v)
= Z(v),
is known as Cartier ’s relation.
For a Gaussian
point
anamorphosis (
station data
),
Z(x) = ϕ(Y (x)) =
∞
X
k
=0
ϕ
k
k!
H
k
(Y (x))
with Hermite polynomials
H
k
and coefficients
ϕ
k
,
the
block
anamorphosis
ϕ
v
(Y (v))
comes as:
ϕ
v
(Y (v)) = E
ϕ(Y (x)) | Y (v)
=
∞
X
k
=0
ϕ
k
k!
r
k
H
k
(Y (v)).
W
W
V
WYX
Z
Point-block-cell correlations
The Gaussian block anamorphosis is:
ϕ
v
(Y (v)) =
∞
X
k
=0
ϕ
k
k!
r
k
H
k
(Y (v)),
with
r
being the
point-block
coefficient
(
0 ≤ r ≤ 1
).
r
can be computed from the block dispersion variance
(which is calculated from the station data variogram):
var(Z(v)) = var(ϕ
v
(Y (v))) =
∞
X
k
=1
ϕ
2
k
k!
r
2k
We get in the same way a
point-cell
coefficient
r
0
.
And finally the
block-cell
coefficient
r
vV
= r
0
/r
.
W
W
V
WYX
Z
Z
Uniform conditioning
It consists in taking the conditional expectation of a
non-linear function of blocks knowing the cell value
containing them.
The
proportion of blocks
v ∈ V
0
above the threshold
z
c
knowing the cell value
Z(V
0
)
is:
E
Z
(v)≥z
c
| Z(V
0
)
= 1 − G
y
c
−
r
vV
Y (V
0
)
√
1 −
r
vV
2
.
G
is the Gaussian distribution.
W
W
V
WYX
Z
Variogram of Airparif measurements
0.
10.
20.
30.
40.
50.
60.
Distance (Kilometer)
0.
1000.
2000.
3000.
4000.
Variogram : Ozone_17JUL15H
0.
10.
20.
30.
40.
50.
60.
Distance (Kilometer)
0.
250.
500.
750.
1000.
Variogram : Ozone_17JUL15H
Nugget-effect + cubic model.
Sill = variance.
W
W
V
WYX
Z
Anamorphosis of Airparif measurements
-2.
-1.
0.
1.
2.
Gaussian values
100.
125.
150.
175.
200.
225.
Ozone
r=.97
-2.
-1.
0.
1.
2.
Gaussian values
100.
125.
150.
175.
200.
225.
Ozone
r’=.72
Anamorphosis of
block
values (r=.97)
close to the anamorphosis of
point
values.
Anamorphosis of
cell
values (r’=.72).
W
W
V
WYX
Z
Histograms
120.
130.
140.
150.
160.
170.
180.
190.
200.
210.
Ozone
0.
10.
20.
30.
Frequencies (%)
Histograms of blocks (
blue
) and cells (
red
)
on the basis of the change-of-support model.
W
W
V
WYX
Z
Proportion of values above threshold
120.
130.
140.
150.
160.
170.
180.
190.
200.
210.
Ozone
0.
10.
20.
30.
40.
50.
60.
70.
80.
90.
100.
Proportion above alert level
Proportions of blocks (
blue
) and cells (
red
).
Depending on the threshold, the difference can be important !
W
W
V
WYX
Z
Uniform conditioning by CHIMERE
0.1
0.1
0.2
0.2
0.3
0.3
0.3
0.4
0.4
0.4
0.5
0.5
0.5
0.6
0.6
0.6
0.7
0.7
0.8
0.8
0.9
0.9
1.4°
1.6°
1.8°
2°
2.2°
2.4°
2.6°
2.8°
3°
3.2°
48.2°
48.4°
48.6°
48.8°
49°
49.2°
49.4°
50 km
UC 120: CHIMERE + Airparif stations
Exceedance probabilities for
1 × 1
km
2
support
with the Swiss threshold of
120
µ
g/m
3
W
W
V
WYX
Z
Uniform conditioning by CHIMERE
0.1
0.1
0.2
0.2
0.3
0.3
0.4
0.4
0.4
0.5
0.5
0.6
0.6
0.7
0.7
0.8
0.9
1.4°
1.6°
1.8°
2°
2.2°
2.4°
2.6°
2.8°
3°
3.2°
48.2°
48.4°
48.6°
48.8°
49°
49.2°
49.4°
50 km
UC 180: CHIMERE + Airparif stations
Exceedance probabilities for
1 × 1
km
2
support
with the European threshold of
180
µ
g/m
3
W
W
V
WYX
Z
Precipitation in SE Norway
geostatistical downscaling
W
W
V
WYX
Z
Histogram of precipitation: July 2001
W
W
V
WYX
Z
Variogram of precipitation
W
W
V
WYX
Z
Z
Block and cell anamorphosis
r=.7
r=.365
10
×
10km
2
blocks
NCEP cells
W
W
V
WYX
Z
Reconstructed histograms
10
×
10 km
2
blocks
101
×
212km
2
NCEP cells
W
W
V
WYX
Z
Proportion above threshold
10
×
10km
2
blocks
NCEP cells
A threshold of 100mm will be used
W
W
V
WYX
Z
Proportion blocks >100mm within NCEP cells
W
W
V
WYX
Z
NCEP cells and station values
Color codes:
0 < x < 75mm <
x
< 100mm <
x
< 125mm <
x
W
W
V
WYX
Z
!
"
#
$
%
&'
()
%
*
&
!
+
&
,
-
.0/
1
2
34
5
/
6
6
67
8
9
:
;=<
>
?
4
<
@
!
A
(
BC
A
D
&
"
#
D
$
A
&
$
E
A
'
D
+
F
B
GH
:
I
J
K
@
I
/
L
@
:
/
8
M
M
8
N
O
-
P
5
:
I
@
#
C
A
Q
+
D
!
"
&
R
&
+
D
$
C
!
D
&
S
F
EUT
T
A
!
D
&
G
H
:
I
V
K
@
I
/
L
@
:
/
N
;
4
:
/
8
M
M
N
XW
Y
Z[\
]
^
Y
_
`a
b[
^
`
^
c
\
c
^
W
Y
d
`fe
gh
i
Document Outline
- Centre de Géostatistique - EMP
- CG: Main application fields
- Geostatistics worldwide
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Geostatistics }} \[1cm] Large definition end {center}
- Geostatistics
- Concepts
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Basic Statistics }} \[1cm] Large concepts end {center}
- Center of mass
- Center of mass
- Center of mass
- Center of mass
- Center of mass
- Histogram
- Histogram
- Cumulative histogram
- Probability distribution
- Probability density
- Expected value
- Variance
- Covariance
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Linear regression }} \[1cm] Large end {center}
- Regression line
- Optimal regression line
- Optimal regression line
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Multiple linear regression }} \[1cm] Large end {center}
- Multivariate data set
- Matrix of means
- Centered variables
- Variance-covariance matrix
- Multiple linear regression
- Minimizing the squared distance
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Simple kriging }} \[1cm] Large end {center}
- Spatial data
- Translation invariance
- Known mean
- Simple kriging equations
- Simple kriging: a multiple linear regression
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Regionalized variables }} \[1cm] Large and random function end {center}
- The concept of a Random Function
- Regionalized Variable
- Epistemological Problem
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Variogram }} \[1cm] Large definition end {center}
- The Variogram
- The Variogram Cloud
- The Experimental Variogram
- The Theoretical Variogram
- Intrinsic Hypothesis
- Definition of the Variogram
- Variogram and Covariance Function
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Variogram }} \[1cm] Large examples end {center}
- Power variogram
- Spherical covariance function
- Exponential covariance function
- Gaussian covariance function
- Cardinal sine covariance function
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Geometric anisotropy }} \[1cm] Large of the variogram end {center}
- Geometric anisotropy
- Rotation in 3D
- 2D example: Ebro river vertical section
- 2D conductivity variogram model
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Behavior at the origin }} \[1cm] Large of the variogram end {center}
- Ebro river:quad water samples
- Nitrate variogram: which behavior at origin?
- Cubic variogram: conditional simulations
- Exponential model: conditional simulations
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Kriging of the mean }} \[1cm] Large of a random function end {center}
- Spatially Correlated Data
- Estimation of the Mean Value
- Stationary random function
- No bias
- Variance of the estimation error
- Minimal estimation variance
- Kriging equations
- Case of no autocorrelation
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Ordinary Kriging }} \[1cm] Large at a point in the domain end {center}
- Estimation at a Point
- Ordinary kriging
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Cross-validation }} \[1cm] Large leaving one out and reestimating it end {center}
- Cross-validation
- Average cross-Validation error
- Standardized cross-validation error
- Average squared Standardized Errors
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Mapping with kriging }} \[1cm] Large on a regular grid\ with irregularly spaced data end {center}
- Kriging for interpolation
- Generating a map
- Moving Neighborhood
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Kriging weights }} \[1cm] Large The shape of the kriging weights end {center}
- Kriging weights
- Isotropic variogram
- Geometric anisotropy
- Relative position of samples
- The screen effect
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Nested variogram }} \[1cm] Large and corresponding linear model\ of the random function end {center}
- Nested Variogram Model
- Example: Arsenic in soil (Loire, France)
- Example: Nested Variogram Model
- Nested Covariance Function
- Regionalization Model
- Linear Model with $S+1$ components
- Kriging Spatial Components
- Example: Short-range Component of As
- Example: Long-range Component of As
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Demographic application }} \[1cm] Large fertility data end {center}
- Demographic application: fertility 1990
- Variograms: class 100-500 women / commune
- Kriging: short range effect only
- Kriging
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Kriging with external drift }} \[1cm] Large end {center}
- Drift
- External drift method
- External Drift method
- Constraint: no bias
- Interpolation of external drift
- Kriging System with linear and external drift
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Kriging temperature }} \[1cm] Large with elevation as external drift end {center}
- Temperature Data
- January temperature vs latitude / longitude
- E-W and N-S temperature variograms
- Kriging mean January temperature
- Temperature vs elevation
- Kriging temperature with elevation as drift
- Temperature estimates vs elevation
- Estimated external drift coefficient
- Estimated external drift coefficient
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Conditional simulation }} \[1cm] Large end {center}
- Conditional simulation vs Kriging
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Change of support }} \[1cm] Large geostatistical simulation of O$_3$ end {center}
- CASE STUDY: Geostatistical simulation of O$_3$
- Simulation of Ozone: 1$ imes $1 Km$^2$ support
- Simulation of Ozone: 10$ imes $10 Km$^2$ support
- Simulation of Ozone: 20$ imes $20 Km$^2$ support
- Simulation of Ozone
- Simulation of Ozone
- Simulation of Ozone
- Simulation of Ozone
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Change of support }} \[1cm] Large concept end {center}
- TOPIC: The Support of a Random Function
- The Effect of Changing the Support
- Neglecting the Support Effect
- Neglecting the Support Effect
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Kriging of a Block average }} \[1cm] Large (centered at a point in the domain) end {center}
- Estimation of a block value
- Block Kriging
- Block kriging with non-point data
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Change of support }} \[1cm] Large risk of exceeding ozone alert level end {center}
- Change of support
- Ozone in Paris on 17 july 1999 at 15h UTC
- Air quality regulations
- Discrete Gaussian point-block model
- Point-block-cell correlations
- Uniform conditioning
- Variogram of Airparif measurements
- Anamorphosis of Airparif measurements
- Histograms
- Proportion of values above threshold
- Uniform conditioning by CHIMERE
- Uniform conditioning by CHIMERE
- egin {center}Huge colorbox {cyan}{ extcolor {blue}{ Precipitation in SE Norway }} \[1cm] Large geostatistical downscaling end {center}
- Histogram of precipitation: July 2001
- Variogram of precipitation
- Block and cell anamorphosis
- Reconstructed histograms
- Proportion above threshold
- Proportion blocks >100mm within NCEP cells
- NCEP cells and station values
Wyszukiwarka
Podobne podstrony:
cours no 09
Tekst czytany, Ma France Cours, Bac
cours no 03
Czas futur proche - teoria, Ma France Cours, Temps
cours no 08
cours no 02
cours no 10
Bon Cours Grammaire Francaise Les Elements De La Langue
przymiotniki opisujące osobowość, Ma France Cours, Słownictwo oraz gramatyki i inne ćwiczenia
Passe recent, Ma France Cours, Temps
Bergengruen Werner ?dekur?s Herzens (gekürzt)
Futur proche cw, Ma France Cours, Temps
Imparfait teoria, Ma France Cours, Temps
cours no 06
cours no 01
Comte Cours? philosophie positive
Bergenstest Kriterier skriftlig produksjon
więcej podobnych podstron