Webster R., Oliver M.A., Geostatistics for Environmental Scientists

Подождите немного. Документ загружается.

There may be outliers, the effects of which we want to retain. Pollutants, for

which only the outliers exceed the statutory thresholds, often fall into this class.

If this were all, how ever, the transformation would be of little practical

signiﬁcance. What makes it of value is that several thresholds can be deﬁned

and a new indicator variable created for each. Thus, if we deﬁne S thresholds as

cð1Þ

; z

cð2Þ

; ...; z

cðSÞ

, then we shall obta in S indicators from the data,

; v

; ...; v

ðxÞ¼1ifzðxÞz

cð1Þ

; else 0;

ðxÞ¼1ifzðxÞz

cð2Þ

; else 0;

ðxÞ¼1ifzðxÞz

cðSÞ

; else 0:

ð11:3Þ

These may be rega rded as the realizations of the corresponding random

functions V

ðxÞ, s ¼ 1; 2; ...; S, for which

ðxÞ¼1ifZðxÞz

cðsÞ

; else 0: ð11:4Þ

The expectation of the indicator, E½V½ZðxÞz

, is the probability, Prob½z

,

that ZðxÞ does not exceed z

Prob½z

¼Prob½ZðxÞz

¼E½V½ZðxÞz

: ð11:5Þ

This probability, Prob½ZðxÞz

, is the cumulative distribution

Prob½ZðxÞz

¼1  Prob½ZðxÞ > z



¼ G½Zðx; z

Þ: ð11:6Þ

Many environmental variables are multi-state charact ers, such as types of

rock, soil and vegetation, that have more than two classes. These can also be

converted to indicators by coding each class as present or absent. If we wished

to distinguish podzols, brown earths, rendzinas and gleys in a region then we

could set up four binary variables, one for each class and code each in turn as

1 or 0. The classes are mutually exclusive, and so in this instance just one of the

four wou ld be coded 1 and the other three as 0.

11.2.2 Indicator variograms

An indicator random function has a variogram

ðhÞ¼

E½fVðx; z

ÞVðx þ h; z

Þg

; ð11:7Þ

The Indicator Approach 247

which is analogous to the variogram of a continuous variable, equation (4.13).

The expected semivariances can be estimated from indicator data by

ðhÞ¼

2mðhÞ

mðhÞ

i¼1

fvðx

; z

Þvðx

þ h; z

Þg

: ð11:8 Þ

Figure 11.1 Indicator variograms of potassium at Broom’s Barn Farm for thresholds of

20, 25 and 30 mg l

1

: (a) autovariograms; (b) cross-variograms.

248 Disjunctive Kriging

Figure 11.1(a) shows some examples. Further, the ordered sets,

ðhÞ, obtained

by applying this formula with changing h and for all thresholds, z

, can be

modelled as described in Chapter 5.

Cross-indicator variograms

For any two thresh olds, say z

and z

, we can deﬁne a cross-indicator variogram

and esti mate it by elaborating the above formulae:

ðhÞ¼

E½fV½Zðx; z

Þ  V½Zðx þ h; z

Þg

fV½Zðx; z

Þ V½Zðx þh; z

Þg: ð11:9Þ

Examples of cross-variograms of indicators appear in Figure 11.1(b).

Indicator covariance functions

If the processes are second-order stationary then the spatial correlations among

the indicators can all be expressed in terms of covaria nces:

ðhÞ¼cov½V½Zðx; z

Þ; V½Zðx þ h; z

Þ

¼ E½V½Zðx; z

ÞV½Zðx þ h; z

Þ  fE½V½Zðx; z

Þg

: ð11:10Þ

Similarly, the cross-covariance at lag h of the indicators for thresholds z

and z

ðhÞ¼cov½V½Zðx; z

Þ; V½Zðx þ h; z

Þ

¼ E½V½Zðx; z

ÞV½Zðx þ h; z

Þ

 E½V½Zðx; z

ÞE½V½Zðx; z

Þ: ð 11:11Þ

11.3 INDICATOR KRIGING

As above, we can krige an indicator variable. So for each target point or block

we compute

Vðx

; z

Þ¼

i¼1

vðx

; z

Þ; ð11:12Þ

where the l

are the weights as usual. This is the ordinary kriged estimate. The

indicator is necessarily bounded, and its sample mean (



v; z

), is usuall y taken as

Indicator Kriging 249

its expectation. We can therefore use simple kriging to estimate Vðx

; z

Þ:

Vðx

; z

Þ¼

i¼1

vðx

; z

Þþ 1 

i¼1

()



v; z

Þ; ð11:13Þ

with weights obtaine d by solving the simple kriging system

i¼1

ðx

; x

; z

Þ¼g

ðx

; x

; z

Þ for j ¼ 1; 2; ...; N; ð11:14Þ

where g

ðx

; x

; z

Þ is the indicator semivariance between the ith and jth

sampling points at threshold z

and g

ðx

; x

; z

Þ is the semivariance of the

indicator between the target point x

and point x

for the same threshold. As

when kriging continuous vari ables, we can replace N by n  N in the

neighbourhood of x

The result is a value lying between 0 and 1 (with exceptions because the

kriging minimizes the variance without any constraint on the estimates it

returns). Such a value is effectively the probability, given the data, that the true

value is 1, i.e.

Prob½Vðx

; z

Þ¼1 jvð x

Þ; i ¼ 1; 2; ...; n¼Ffx

jðnÞg; ð11:15Þ

where we use ðnÞ to mean all the data in the particular neighbourhood. The

quantity Ffx

jðnÞg denotes the conditional or ‘posterior’ probability that Vðx

is 1.

If now we return to our problem, namely to estimate the prob ability, given

data, tha t the true value of Z at an unsampled place x

does not exceed z

, then

we can write

Prob½Zðx

Þ4z

jzðx

Þ; i ¼ 1; 2; ...; N

¼ 1  Prob½Zðx

Þ > z

jzðx

Þ; i ¼ 1; 2; ...; N: ð11:16Þ

Notice that the two sides of equation (11.16) are complementary.

At ﬁrst sight it might seem that the way to tackle the problem is to transform

the data to indicators determined by the particular threshold. However, we soon

see that an individual probability estimated in this way is crude. Much of the

rich information in the original data has been lost by dissecting the scale into

just two classes.

This loss can be made good to a large extent by repeating the process for

several thresholds in the range of Z and constructing a cumulative distribution

function, conditional on the data, for each target point by accumulating the

Fðx

; z

Þ, s ¼ 1; 2; ...; S.

250 Disjunctive Kriging

The procedure is somewhat tedious because the variograms for all the

thresholds must be computed and modelled. Furthermore, because the

Fðx

; z

Þ for the different z

are computed independently of one another, there

is no guarantee that they will sum to 1, or that the cumul ative function will

increase monotonically, or that the estimated probabilities will lie in the range 0

to 1. Some adjustment of the results may therefore be needed to ensure that the

bounds are honoured and the order relations mai ntained. Nevertheless, an

empirical distribution function can be obtained and then use d to reﬁne the

estimate of the conditional probability that Zðx

Þz

Goovaerts (1997) describes the procedu re fully and illustrates it with

examples using the computer programs in GSLIB (Deutsch and Journel,

1992), while Olea (1999) devotes a section of his book to the topic. We shall

not repe at the detail here.

11.4 DISJUNCTIVE KRIGING

Disjunctive kriging provides another way of estimating an indicator transform

of continuous data. It does so without losing information, though requiring

rather stronger assumptions than does indicator kriging as described above. It

may take several forms (see Rivoirard, 1994), the most common of which is

Gaussian disjunctive kriging and the one we describe.

11.4.1 Assumptions of Gaussian disjunctive kriging

The assumptions underlying Gaussian disjunctive kriging are as follows. First,

zðxÞ is a realization of a second-order stationary process ZðxÞ with mean m,

variance s

and covariance function CðhÞ. The underlying variogram must

therefore be bounded. Second, the bivariate distribution for the n þ 1 variates, i.e.

for each target site and the sample locations in its neighbourhood, is known and

is stable throughout the region. If the distribution of ZðxÞ is normal (Gaussian)

and the process is second-order stationary then we can assume that the bivariate

distribution for each pair of locations is also normal. Each pair of variates has the

same bivariate density, and this density function is determined from the spatial

autocorrelation coefﬁcient. These assumptions allow the conditional expectations

to be written in terms of the autocorrelation coefﬁcients, as we shall show.

The variable ZðxÞ is spatially continuous, so that in going from a small value

at one place to a large one elsewhere it must pass through intermediate values

en route. It is an example of a Gaussian diffusion process. One test of this

assumption is to compare the variograms of the indicators for several thresholds

within the bounds of the measured z. The cross-indicator variograms should be

more ‘structured’ than the autovariograms. Figure 11.1 shows this to be so for

potassium at Broom’s Barn Farm.

Disjunctive Kriging 251

11.4.2 Hermite polynomials

The requirement of normality is a strong one that is rarely met in practice, even

though many environmental properties seem approximately normal. The ﬁrst

task therefore is to transform an actual distribution of ZðxÞ, which may have

almost any form, to a standard normal one, YðxÞ, such that

ZðxÞ¼F½YðxÞ: ð11: 17Þ

This can be done with Hermite polynomials.

We recall, from Chapter 2, the equation of the standard normal distribution

with probability density given by

gðyÞ¼

ﬃﬃﬃﬃﬃﬃ

exp 



: ð11:18Þ

Hermite polynomials are related to this function and are deﬁned by Rodrigues’s

formula as

ðyÞ¼

ﬃﬃﬃﬃ

gðyÞ

: ð11:19Þ

Here k is the degree of the polynomial, taking values 0; 1; ..., and 1=

ﬃﬃﬃﬃ

is a

standardizing factor (Matheron, 1976). The ﬁrst two Hermite polynomials, i.e.

for k ¼ 0 and k ¼ 1, are

ðyÞ¼1; ð11:20Þ

ðyÞ¼y; ð11:21Þ

and thereafter the higher-order polynomials obey the recurrence relation

ðyÞ¼

ﬃﬃﬃ

k1

ðyÞ

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

k  1

k2

ðyÞ: ð11:22Þ

So the polynomials can be calculated up to any order for a standard normal

distribution.

The Hermite polynomials are orthogonal with respect to the weighting

function expðy

=2Þ on the interval 1 to þ1. They are independent

components of the normal distribution of ever increasing detail.

Almost any function of YðxÞ can be represented as the sum of Hermite

polynomials:

f fYðxÞg ¼ f

fYðxÞg þ f

fYðxÞgþ; ð11:23Þ

252 Disjunctive Kriging

and sinc e the Hermite polynomials are orthogonal

E½f fYðxÞgH

fYðxÞg ¼ E H

fYðxÞg

l¼0

fYðxÞg

l¼0

E½H

fYðxÞgH

fYðxÞg

¼ f

: ð11:24Þ

This enables us to calculate the coefﬁcients f

of F½YðxÞ in equation (11.17) as

ZðxÞ¼F½YðxÞ

¼ f

fYðxÞg þ f

fYðxÞgþ

k¼0

fYðxÞg: ð11:25Þ

The transform is also invertible, which means that the results can be expressed

in the same units as the original measurements.

Determining the Hermite coefﬁcients

To determine the coefﬁcients of the transformation, f

, for a particular set of

data we proceed as follows. We arrange the N data in ascending order:

< z

< < z

;

and we denote their relative frequencies by

; q

; ...; q

;

such that the sum of the frequencies is 1:

i¼1

¼ 1. Their cumulative

frequencies are

Fðz

Þ¼Prob½ZðxÞ < z

¼0;

Fðz

Þ¼Prob½ZðxÞ < z

¼q

;

Fðz

Þ¼Prob½ZðxÞ < z

¼q

þ q

;

Fðz

Þ¼Prob½ZðxÞ < z

¼

i1

j¼1

;

Fðz

Þ¼Prob½ZðxÞ < z

¼1  q

Disjunctive Kriging 253

The cumulative frequencies have equivalents on the standard normal

distribution:

Fðz

Þ¼Gðy

Þ: ð11:26Þ

Thus

Fðz

iþ1

ÞFðz

Þ¼Gðy

iþ1

ÞGðy

Þ; ð11:27Þ

and so

Prob½z

 ZðxÞ < z

iþ1

¼Prob½y

 YðxÞ < y

iþ1

;

Prob½ZðxÞ¼z

¼Prob½y

 YðxÞ < y

iþ1

:

ð11:28Þ

In words, ZðxÞ equals z

when the standard normal equivalent lies between

and y

iþ1

. We can then determine the transformation coefﬁcients as follows:

¼ E½FfYðxÞg ¼ E½ZðxÞ ¼

i¼1

; ð11:29Þ

and thereafter

¼ E½ZðxÞH

fYðxÞg

þ1

1

FðyÞH

ðyÞgðyÞdy

i¼1

iþ1

ðyÞgðyÞdy

i¼1

ﬃﬃﬃ

k1

ðy

iþ1

Þgðy

iþ1

Þ

ﬃﬃﬃ

k1

ðy

Þgðy



i¼2

ðz

i1

 z

ﬃﬃﬃ

k1

ðy

Þgðy

Þð11:30Þ

because gðy

Þ¼gð1Þ ¼ 0, and gðy

Nþ1

Þ¼gðþ1Þ ¼ 0 also.

11.4.3 Disjunctive kriging for a Hermite polynomial

Since the polynomials are orthogonal any pair of values, YðxÞ and Yðx þ hÞ,

drawn from a bivariate normal distribution with correlation coefﬁcient r has

expectation

E½H

fYðxÞgjYð x þ hÞ ¼ r

ðhÞH

fYðx þ hÞg: ð11:31Þ

254 Disjunctive Kriging

The covariance between two functions of Y at x and x þ h is

cov½H

fYðxÞg;H

fYðx þ hÞg

¼ E½H

fYðxÞgH

fYðx þ hÞg

¼ E½H

fYðxÞgE½H

fYðx þ hÞgjYðx þ hÞ

¼ r

ðhÞE½H

fYðxÞgH

fYðxÞg: ð11:32Þ

When k ¼ l this equation gives the covariance of H

fYðxÞg, which equals r

ðhÞ,

since YðxÞ is a standar d normal variate. The correlation coef ﬁcient, rðhÞ,

h 6¼ 0, must lie between 1 and þ1; so r

ðhÞ rapidly approaches 0 as k

increases, and the spat ial dependence in H

fYðxÞg declines to nothing, i.e.

fYðxÞg becomes pure nugget.

Any pair of Hermite polynomials is spatially independent, so they are the

independent factors of the bivariate normal model. By kriging them separately

the estimates have only to be summed to give the disjunctive kriging estimator:

ðxÞ¼f

þ f

fYðxÞg þ f

fYðxÞgþ: ð11:33Þ

So, if we have n points in the neighb ourhood of x

where we want an estimat e,

we estimate the Hermite polynomials by

fYðx

Þg ¼

i¼1

fYðx

Þg; ð11:34Þ

and we insert them into equation (11.33). The l

are the kriging weights,

which are found by solving the equations for simple kriging because we can

assume the mean is known:

i¼1

cov½H

fYðx

Þg; H

fYðx

Þg

¼ cov½H

fYðx

Þg; H

fYðx

Þg for all j; ð 11:35Þ

or alternatively,

i¼1

ðx

 x

Þ¼r

ðx

 x

Þ for all j; ð 11:36Þ

from equation (11.31). In particular, the procedure enables us to estimate

Zðx

Þ by

Zðx

Þ¼Ff

Yðx

Þg ¼ f

þ f

fyðx

Þg þ f

fyðx

Þg þ : ð11:37Þ

Disjunctive Kriging 255

11.4.4 Estimation variance

The kriging variance of

fYðxÞg is

ðx

Þ¼1 

i¼1

ðx

 x

Þ;

ð11:38Þ

and the disjunctive kriging variance of

f ½Yðx

Þ is

ðx

Þ¼

k¼1

ðx

Þ: ð11:39Þ

11.4.5 Conditional probability

Once the Hermite polynomials have been estimated at a target point we can

estimate the conditional probability that the true value there exceeds the critical

value, z

. The transformation ZðxÞ¼F½YðxÞ means that z

has an equivalent

on the standard normal scale. Since the two scales are monotonically related

their indicators are the same:

V½ZðxÞz

¼V½YðxÞy

: ð11: 40 Þ

For V½YðxÞ > y

, which is the complement of V½YðxÞy

, the kth Hermite

coefﬁcient is

þ1

1

V½y  y

H

ðyÞgðyÞdy

1

ðyÞgðyÞdy: ð11:41Þ

The coefﬁcient for k ¼ 0 is the cumulative distribution to y

¼ Gðy

Þ;

and for larger k,

ﬃﬃﬃ

k1

ðy

Þgðy

Þ: ð11:42Þ

The indicator can be expressed in terms of the cumulative distribution and the

Hermite polynomials:

V½YðxÞy

¼Gðy

Þþ

k¼1

ﬃﬃﬃ

k1

ðy

Þgðy

ÞH

fYðxÞg: ð11:43Þ

256 Disjunctive Kriging