Jean-Laurent Mallet. Geomodeling
Подождите немного. Документ загружается.
498
CHAPTER
9.
STOCHASTIC MODELING
As
can be
seen,
the first
constraint above implies
that
the
model (9.109)
reproduces correctly
the
mean value
ms(u)
of
S^u,
a;)
on
D°:
A
a
consequence,
the
only thing
to do is to
look
for a
matrix
A
ensuring
that
the
covariances
defined
by
equation (9.108)
are
reproduced
on
D°
x
£)°.
Determining
the
matrix
A
Let [C] be the
symmetrical
(M
x
M)
covariance matrix
and let
[<r]
be the
diagonal
(M x M)
standard deviation matrix
defined
by
In
the
above matrices,
the
components
dj
and
<Ji
are
assumed
to be
deduced
as
follows
from
the
covariance
function
C(h)
and the
standard deviation
func-
tion
crs:(u)
used
in
equation
(9.108):
It is
easy
to
check
that
the
covariance matrix
[Css]
°f
the
random vector
S(u)
defined
by
equation (9.109)
is
such
that
Moreover,
taking into account equation (9.109)
and
(9.110),
it can be
observed
that
Let
{Vi,...,
VM}
be
the
unit
eigen vectors
of
[C]
corresponding
to the
eigen
values
{AI,
...,
AM}
assumed
to be
sorted
by
decreasing order
of
magnitude.
Taking into account
the
fact
that
[C]
is
symmetrical
and
semipositive definite,
it is
well
known
(e.g.,
see
[19])
that
• the
eigen values
of
[C]
are
real
and
nonnegative,
and
• the
matrix
[C]
can be
decomposed
as
follows:
Consequently,
if
there
are K
eigen values significantly
different
from
zero,
then equations (9.114)
and
(9.116) imply
that
we can
write
In
conclusion,
if we
identify
the
columns
of the
matrix
A in
equations (9.109)
and
(9.115) with
the K
vectors
{Ak =
\f\k
•
[<r]
•
Vk
'•
k
=
1,...,
K},
then
the
covariances
of the
components
of the
vector
S(uj]
corresponding
to the
sam-
pling
on
D^
of the
simulator S(u,
u>)
defined
by
equation (9.108)
are
honored.
9.6.
STOCHASTIC SIMULATORS
Algorith m
As
a
consequence
of the
above model,
the
following
procedure
to
generate
realizations
of
5(u,
o;)
on
D^
is
proposed:
•
Using equations (9.108)
and
(9.114), build
the
matrix [C].
•
Compute
the K
<
M
non-null
eigen values
{Afc}
of
[C]
and
their
associated
unit
eigen vectors
{Vk}.
•
Build
the
matrix
A of
equation (9.109)
as
follows:
• For
each simulation
associated
with
an
even
uj
G
15
—
use a
random number generator
to
generate
a
realization
of
X(u>)
hon-
oring
the
constraints
(9.110),
and
—
use the
equation (9.109)
to
generate
the
sampling
of
<S(u,
u) on
£)°.
It can be
observed
that
the
column vectors
[\f\k
•
Vk}
used
to
build
the
matrix
A are no
more than
the
coefficients
of the
Principal Components
of
the
covariance
matrix
[C].
For
this
reason,
it is
proposed
that
we
call
"Principal-Components-based
simulator"
or,
more simply, "PC-based simula-
tor"
the
stochastic simulator
S(u,u;)
so
denned.
Comment
When
the set
D^
of the
sampling points
is
sufficiently
dense
on D,
each
realization
5(11,0;)
can be
continuously interpolated
on D
and,
at the
limit,
equation (9.109) takes
the
following
continuous
form:
In the
one-dimensional case where
the
parametric domain
D is a
part
of M,
it can be
observed that,
at
each
data
point
Uj
G
-D*,
the
derivative
of
0"s(u)
generally
vanishes,
and we
have
In
other words,
on the set D* of the
data
points,
the
derivative
of
each real-
ization
S(u,
o;)
is
equal
to the
derivatives
of
ms
(u): this
may be a
drawback
in
practical applications.
It is
easy
to
check
that
a
similar property occurs
for
the
partial derivatives
of
each realization relative
to the
components
of u in
the
p-dimensional
case where
D is a
part
of
]R
P
.
9.6.4
Example
3:
Residual-based
simulator
This section shows
how an
interpolation method
can be
used
to
transform
an
unconstrained simulator into
a
constrained simulator.
499
500
CHAPTER
9.
STOCHASTIC MODELING
Definition
of the
simulator
Once
again,
let us
address
the
problem
of
building
a
stochastic simulator
,S(u,
uj]
constrained
to
honor
the
following
constraints where
{z(ui)
:
u^
e
D*}
is
a
given
set of
sampled values
for the
unknown
function
z(u)
to be
simulated
on D D D*
while
C(h)
is a
given stationary covariance
function:
1)
2)
In
this section,
S(u,
a;)
is
assumed
to be
denned
as the sum of two
random
functions
SQ(U,UJ)
and
R(U,UJ)
where
So(u,
u;)
is a
given unconstrained simulator such
that
while
.R(u,
o>)
is a
"residual" random
function
honoring
the
following
con-
straint:
As
a
consequence
of
equations (9.120), (9.121)-1,
and
(9.122), with
a
simulator
so
denned,
the
constraint (9.119)-1
is
always strictly honored:
To
build
a
residual random
function
-R(u,
u;)
honoring
the
constraints
(9.122),
one can
proceed
as
follows:
• Use any
stochastic method
to
build
an
unconstrained simulator
SQ
(u,
u;)
whose
mean
is
chosen equal
to
zero while
its
covariance
function
is
chosen
equal
to
C(ui
—
u
2
).
• For
each event
uj
€
O,
build
the
realization
-R(u,
u>)
as
follows:
1.
for
each
data
point
Ui
£
D*,
compute
the
residual
R(ui,u>)
as
defined
by
equation
(9.122);
2.
use an
interpolation
method
to
build
a
(non-random)
function
R(u,
u)
denned
on D and
interpolating
the
values
{R(ui,u)
:
Ui E
D*}.
Choosing
an
interpolation
method
for
R(U,UJ)
Proceeding
as
suggested above ensures
that
the
simulator
5(u,
u;)
defined
by
equation
(9.120)
honors
the
constraint (9.119)-1
but
there
is no
reason
for
the
constraint
(9.119)-2
to be
honored.
Let us
assume
now
that
we
choose
to
interpolate each realization
of
R(u,
u;)
with
an
interpolation method honoring
the
following
constraint:
u far
away
from
the
data
points
uu
£ D*
9.6.
STOCHASTIC SIMULATORS
Figure
9.11
An
example
of
a
series
of
50
equiprobable
constrained
simulations
(«S(u,
uji)
:
i = 1,
2,...,
50},
obtained
with
a
residual
based
simulator.
In
such
a
case,
as a
straightforward consequence,
it can be
observed
that
the
constraint
(9.119)-2
is now
honored:
(111,112)
far
away
from
the
data
points
m
e
D*
Moreover,
an
interesting property
of the
simulator
5(11,0;)
so
defined
is
that,
far
away
from
the
data
points,
the
mean value
of
S(u,u>)
is
close
to the
constant
z:
u far
away
from
the
data
points
u^
6 D*
Whatever
the
interpolation method used
to
build
-R(u,
a;),
this suggests choos-
ing z as
follows:
This
natural choice ensures
that,
far
away
from
the
data
points
and for any
(jj
e
O, the
values
of the
simulator
S(u,
uj)
have
the
same order
or
magnitude
as the
data
values
[z(ui]
:
Uj
e
D*}.
There
are
many possible choices
for an
interpolation method allowing
the
interpolated
function
R(U,UJ)
to
honor constraints (9.123).
For
example,
in
the
case where
R(u,
w)
has to be
evaluated
at the
nodes
of a
discrete
model,
the
simplest solution consists
in
using
the
iterative version
of the
DSI
method
(see
page 165) with
the
value
z fixed as a
Control-Node value (see page 147)
for
any
node located
far
away
from
the
data
points while
the
values
{R(ui,u)
:
Ui
G
D*}
at the
data
points
are set as
Control-Properties constraints (see
page 183).
As
a
tutorial example,
figure
(9.11)
shows
a
series
of
equiprobable realiza-
tions
of a
residual-based
simulator defined
on a
segment
D and
using
the DSI
501
502
CHAPTER
9.
STOCHASTIC MODELING
interpolation method.
It can be
observed that, contrary
to the
simulations
obtained with
the
P-field
and
PC-based simulators, these simulations
do not
have
a
constant derivative
at the
data
points
{u^
G
D*}.
9.6.5
Other stochastic simulators
Besides
the
three methods presented above, many
different
methods have been
proposed
in the
literature
to
build stochastic simulators.
For
example,
• the
turning
bands
method
[157,
219];
•
sequential
methods
based
on
Kriging
[58];
•
blending-based
methods
[151];
•
simulating
annealing
methods
[94];
and
•
methods
based
on
neural
networks
[39].
The
rest
of
this chapter
is
dedicated
to the
presentation
of
some
of
these
methods.
9.7
Kriging-based
methods
Since
its
introduction
by
Matheron [156]
in the
late
1960s, geostatistics
has
contributed
a
great
deal
to the
development
of
spatial
interpolation
over
the
past
decades (see also Whittle
[233]).
During
the
1980s,
the
main
focus
of
this
family
of
methods
shifted
from
the
estimation
of
physical parameters
distributed over
a
region
to the
generation
of
multiple simulations
of
these
parameters. This section
focuses
briefly
on the
notion
of
"Kriging," which
is
of
particular interest
in the
building
of
simulators. Those interested
in
a
complete presentation
of
geostatistical methods
are
referred
to the
many
excellent
books
on the
subject,
for
example, [114, 119, 111,
58,
224,
94,
45].
9.7.1
The
interpolation problem
Considering
further
the
interpolation versus simulation problem introduced
in
section (9.6.1),
an
interpolation strategy
can be
adopted whose
aim is to
build
a
function
z*(u)
defined
on D and
assumed
to
interpolate
the
unknown
function
z(u)
on
D*:
To
build such
an
interpolator
z*(u),
the
traditional geostatistics approach
would
proceed
in
this way:
1.
Take
the
model decision that
there
is an
unknown
second-order
station-
ary
ergodic
(see
pages
461 and
462)
random
function
Z(u,
cj)
honoring
the
following
constraint:
9.
2.
Use the
ergodicity
of
Z(u,u)
(see equation
(9.38))
for
building
a
variogram
model
7z(h)
from
the N
pairs
of
data
{u^,
z(iii)}
(see page 467);
3.
use the
variogram model
72
(h),
or,
equivalently,
the
associated covariance
function
Cz(h)
for
building
a
linear estimator
Z*(u,u)
of
Z(u,uj)
such
that
where
the
weighting
coefficients
(Ai(u)}
have
to be
determined
in an
"opti-
mal" way;
4.
Use
equations (9.125)
and
(9.126)
for
building
an
estimate
z*(u)
of
2(11)
de-
fined
on
D by
where
UQ
is the
same
elementary statistical event
as the one
corresponding
to
equation (9.125). This equation shows
that
the
weighting
coefficients
(Aj(u)}
need
to be
constrained
to
have their
sum
equal
to
one
39
because,
in the ex-
treme case where
all the
values
{z(ui)}
are
equal
to a
constant,
we
would like
z*(u)
to be
equal
to
this constant:
In
steps
(1) and
(3),
the
words
"model
decision"
and
"use"
have
been
written
in
bold
to
avoid
misunderstanding:
• It is
extremely important
to
note
that
taking
the
decision
that
there
is an
unknown second-order station-
ary
ergodic random function honoring constraint
(9.125),
is not a
hypothesis
but a
model that
is
chosen
a
priori.
Choosing such
a
model
is
similar
to
choosing
a
basis
of
polynomials when
using
a
spline
or
Bezier
interpolator
to
model
a
curve
or a
surface.
• It is
also just
as
important
to
note
that
"using"
a
specific
variogram
of
Z(u,
ui)
for
building
the
estima-
tor
Z*(u,u>)
does
not
imply
that
Z*(u,
o>)
should have
the
same
variogram
as
Z(u,
u>).
For
example,
if all the
values
{2(11^)}
are
equal
to
zero, then combining
the
equations
(9.125)
and
(9.127) shows
that
z*(u)
will
be
identical
to
zero what-
ever
the
variogram
of
Z(u,u;).
Over
the
past
three
decades,
geostatisticians
have
developed
a
large
number
of
methods
collectively
known
as
"Kriging"
methods
based
on the
strategy
described
above.
The
next
section
focuses
on
"Ordinary
Kriging"
and
"Simple
Kriging,"
the
simplest
and
also
the
most
popular
of
these
methods.
39
Note
that
this
is a
sufficient
condition,
not a
necessary
condition.
5
5
504
CHAPTER
9.
STOCHASTIC MODELING
9.7.2
Ordinary
Kriging
interpolator
Definition
Among
the
infinite
number
of
possible weightings
(Aj(u)}
used
in
equation
(9.126),
Ordinary Kriging chooses
the one
that
minimizes
the
so-called "esti-
mation variance
"
<T|;(U),
as
denned
by
The
optimal weighting
coefficients
(A^(u)}
minimizing
(9.129),
while honoring
constraint (9.128),
can be
obtained
using
the
classic
Lagrange multipliers
tech-
nique [141,
90,
163].
According
to
this technique,
the
optimal solution should
realize
the
minimum
of the
function
£(Ai(u),...,
AJV(U),
/u(u))
as
defined
by
where
/Lt(u)
is the
so-called "Lagrange multiplier"
associated
with
constraint
(9.128).
Consequently,
the
following
equations must
be
honored:
Taking into account
the
fact
that
Z(u,
u)
is
assumed
to be
order
2
stationary,
it is
easy
to
show
that
where
C^(h)
is the
stationary covariance
function
of
Z(u,a;).
The
Lagrange
equations (9.130)
can
thus
be
turned into
the
following
system
of (N + 1)
linear
equations called "Ordinary Kriging
equations:"
In
other words, assuming
that
[K]
is
invertible,
40
the
optimal weightings
are
such
that
40
With
the
covariance functions presented
in
section (9.3.3),
it can be
shown
that
this
is
always
the
case provided
that
the
sampling points
u^
G D* are
distinct.
9.7.
KRIGING-BASED
METHODS
505
Properties
As
an
exercise
it can be
checked that
the
associated minimum
of the
estimation
variance
is
such
that
Moreover,
in the
case where
u is
identical
to one of the
sampling points,
say
Uj,
then
F(u)
—
F(UJ)
is
identical
to the
jib.
column
of
[K]
and the
solution
A(u)
of
equation (9.132)
is
thus such
that
otherwise
According
to
equation
(9.127),
this implies
that
the
Ordinary Kriging
is an
"exact" interpolator
in the
sense that,
if u is
identical
to a
data
location
Uj,
then
the
estimated value
z*(uj)
is
equal
to the
input data:
For
the
same reason,
the
estimation variance
cr|;(u)
vanishes
at
data
points:
If
the
distance
of u to the
data
points
is
denned
as
follows
u far
away from
data
points
then
it is
verifiable
41
that
the
estimation variance behaves
as
follows
when
far
away
from
the
data
points:
u far
away
from
data
points
Example
Ordinary
Kriging
is
widely used
in
geomodeling
for
interpolating physical
parameters
of
natural objects. Figure
(9.12)
shows
an
example corresponding
to a
data
set
given
in
[58] where
the
goal
was to
interpolate
a
parameter
z(a)
=
2(11(0:))
at
each node
of a
regular 2-grid:
• the
Ordinary
Kriging
map
z*(u)
does
interpolate
the
data
points
and
looks
smooth
when away from
these
data;
and
• the
variance
of
estimation
cr%(u)
does
vanish
at
data
points
and is
close
to
a
constant
for
points
u
whose
distance
to the
closest
data
point
is
greater
than
the
range
of the
covariance function
(or the
variogram).
According
to
equation
(9.137),
this
constant
value
is
within
the
range
[Cz(0),2
•
Cz(0)].
41
This
property comes
from
equation (9.131)
and the
fact
that
the
hypothesis (9.136)
and
the
constraint
(9.128)
imply
that
u far
away
from
data
points
u far
away from
data
points
otherwise
506
CHAPTER
9.
STOCHASTIC MODELING
Figure
9.12
Comparison
between
the
DSI
interpolation
(B) and
Ordinary
Kriging
interpolation
(C)
computed
at the
200
2
nodes
of
a
regular
rectilinear
2-grid
using
the
same
data
set
(A).
Part
(D)
shows
how the
variance
of
estimation
cr
E
(u)
of
the
Ordinary
Kriging
varies
as a
function
of the
distance
to the
data
points.
In
these
pictures,
low
values
are
coded
in
light
grey
while
high
values
are
coded
in
dark
grey.
In
this example,
a
Spherical covariance
function
model
Cz(h)
was
used (see
definition
on
page 464) with
an
isotropic range
R
equal
to
|
of
the
width
of
th e
studied domain.
It can be
observed
that,
as
with
any
other interpola-
tion method, both
the DSI
interpolation presented
in figure
(9.12)-B
and the
Ordinary
Kriging presented
in
figure
(9.12)-C
produce smooth solution.
9.7.3 Simple Kriging interpolator
Definition
Simple
Kriging
is a
simplified
version
of
Ordinary Kriging where
the
con-
straint (9.128)
is
removed.
42
As a
consequence, Lagrange multipliers
are no
42
Generally,
the
mean
mz(u)
is
assumed
to be
known
a
priori
and the
discrepancies
relative
to
this mean
are
interpolated:
For
example,
if the
data points
are not
clustered,
it can be
assumed that
mz(u)
=
mz
and
mz
can be
estimated
as the
average value
of the
observed data
{z(uj)}.
9.7.
KRIGING-BASED
longer needed
and the
weighting
coefficients
(Ai(u)}
are now
solutions
of the
following
N
linear equations:
According
to
equations (9.131),
it is
deduced
that
This
linear system
can
also
be
turned into
the
following
matrix equation called
the
"Simple Kriging equation," which replaces equation
(9.132):
Properties
It is
easy
to
check
that
Simple Kriging
has the
same properties (equations
(9.133)-(9.136))
as
Ordinary
Kriging
except
that,
due to the
removal
of
con-
straint
(9.128),
we now
have:
However,
it
should
be
noted
that
equation (9.137)
is now
replaced
by
43
In
addition
to the
above property
and
properties
(9.133)-(9.136)
shared with
Ordinary Kriging,
it can be
observed
that
43
This
property
comes
from
equation
(9.131)
and the
fact
that
the
hypothesis
(9.136)
and the
equation
(9.139)
imply
that
all the
coefficients
(Aj(u)}
vanish
when
u is far
away
from
data points.
5
n addition to the above property and