W O R K I N G P A P E R
K R I G I N G AND OTHER E S T I M A T O R S O F S P A T I A L F I E L D C H A R A C T E R I S T I C S
( W I T H S P E C I A L R E F E R E N C E T O ENVIRONMENTAL S T U D I E S )
.
V . Fedorov
O c t o b e r 1 9 8 7 W P - 8 7 - 9 9
lnternat~onal lnst~tute for Appl~ed Systems Analys~s
A-2361 LaxenburgIAustr~a
EI I IASA
KRIGING
AND OTHER
ESTIMATORS OF SPATIALFJELD
CHARAWERISI'ICS(WITH SPEIAL
REFERENCE
TOENVIRONMENTAL SlVDIES).
O c t o b e r 1 9 6 7 WP-87-99
Working Papers are interim r e p o r t s on work of t h e I n t e r n a t i o n a l I n s t i t u t e f o r Applied Systems Analysis a n d h a v e r e c e i v e d only limited review. Views o r opinions e x p r e s s e d h e r e i n d o n o t n e c e s s a r i l y r e p r e s e n t t h o s e of t h e I n s t i t u t e or of i t s National Member Organizations.
INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS A-2361 L a x e n b u r g , Austria
Foreword
During t h e last t w o y e a r s , Valeri Fedorov h a s b e e n bringing h i s v e r y consid- e r a b l e s t a t i s t i c a l t a l e n t s to b e a r on t h e design of environmental monitoring sys- tems a n d on t h e analysis of e x p e r i m e n t a l d a t a in t h e environmental s c i e n c e s . In t h i s p a r t i c u l a r Working P a p e r , Valeri h a s examined t h e c o n c e p t of kriging, a method used t o r e c o v e r s p a t i a l p a t t e r n s from point measurements, f i r s t a p p l i e d in t h e geological s c i e n c e s . Unfortunately, environmental fields usually v a r y in time as well as s p a c e . This l e a d s to s e v e r e difficulties which h a v e n o t b e e n fully r e a l - ized by environmental s c i e n t i s t s . The n a t u r e of t h e s e difficulties i s e l a b o r a t e d in t h i s p a p e r .
R.E. Munn L e a d e r
Environmental P r o g r a m
KRIGING AND
CYlTtEZ
ESTIMATORS OF SPATIAL FIEU) CHARACl'ERISTICS(WITH
SPECIAL REFERENCE
TO J3lWIRONMEWTALSrUDIES).
K Fedorov
1. Introduction
During t h e l a s t d e c a d e u s e of t h e k r i g i n g method ( o r simply "kriging") in- c r e a s e d g r e a t l y in r e s e a r c h r e l a t e d to t h e a n a l y s i s of s p a t i a l v a r i a b i l i t y of en- v i r o n m e n t a l p a r a m e t e r s ; see f o r i n s t a n c e , B a r n e s (1980), Dennis a n d S e i l k o p (1986), C l a r k et a l . (1986), Endlich et a l . (1986), Finkelstein a n d S e i l k o p (1981), D e r M e g r e d i t c h a n (1979) a n d McBratney a n d W e b s t e r (1981). I t i s d i f f i c u l t to u n d e r s t a n d t h e r e a s o n f o r t h i s i n c r e a s i n g p o p u l a r i t y ( b u t with some o c c a s i o n a l d a r k s p o t s ; see f o r i n s t a n c e , Akima (1974), Armstrong (1984)). I s i t i t s c o m p a r a - t i v e simplicity, s t a t i s t i c a l e f f e c t i v e n e s s , mathematical e l e g a n c y o r t h e mesmerizing impact of t h e word "kriging" (sometimes "universal kriging")? In t h i s p a p e r a n at- t e m p t to c l a r i f y t h e s i t u a t i o n a n d t o u n d e r s t a n d i t s p l a c e in t h e s t a t i s t i c a l t h e o r y of s p a t i a l p a t t e r n a n a l y s i s a n d t h e admissibility of t h e k r i g i n g method in e n v i r o n - mental s t u d i e s i s u n d e r t a k e n .
T e c h n i c a l d e t a i l s a n d computing a s p e c t s of t h e k r i g i n g method are a v o i d e d b u t t h e r e a d e r c a n find them f o r i n s t a n c e , in J o u r n a l a n d H u i j b r e g t s (1978), Ripley (1983) a n d Thiebaux a n d P e d d l e r (1987). Tho b i b l i o g r a p h i c o v e r v i e w by Bell a n d R e e v e s (1979) p r o v i d e s a g u i d e t o t h e o r i g i n a l t h e o r e t i c a l l i t e r a t u r e (mainly by K r i g e , M a t h e r o n a n d J o u r n e l ) a n d t h e n u m e r o u s a p p l i e d p a p e r s .
I t i s worthwhile to n o t e t h a t t h e method i s t h e a n a l o g u e f o r s p a t i a l p r o c e s s e s of Wiener-Kolmogorov p r e d i c t i o n t h e o r y . I t "has b e e n d e v e l o p e d a n d u s e d mainly by Matheron a n d h i s s c h o o l in t h e mining i n d u s t r y . H e c h r i s t e n e d t h e method k r i g - ing a f t e r D.G. K r i g e " (Ripley (1983) s e c t i o n 4.4).
U. Linear Estimators. Kriging
L e t some v a l u e y ( z ) b e o b s e r v e d at p o i n t s z l , z z ,
. . . ,z,
E X, w h e r e X i s t w o dimensional in t h e m a j o r i t y of a p p l i c a t i o n s . P o i n t s x i , i = l , n , could b e s p a c e d at t h e i n t e r s e c t i o n p o i n t s of some r e g u l a r g r i d b u t i t i s n o t e s s e n t i a l f o r t h e t h e o r y a n d d o e s n o t c a u s e s e r i o u s d i f f i c u l t i e s f o r modern s o f t w a r e if t h e y are not.T h e estimation of v a l u e y a t t h e p r e s c r i b e d p o i n t z 0 i s of i n t e r e s t in t h e a n a l y s i s of s p a t i a l p a t t e r n s . F o r t h e p u r p o s e of t h i s p a p e r , i t i s s u f f i c i e n t to con- s i d e r only t h e l i n e a r e s t i m a t o r of y o
=
y ( z o ) :w h e r e "T" s t a n d f o r t r a n s p o s i t i o n , h E R n ,
=
( y l , y 2 ,.
..
, y,), yi=
y ( x i ).
A p r a c t i t i o n e r wishes to h a v e y^, as c l o s e to i t s t r u e v a l u e y o as possible. In o t h e r words, o n e h a s to minimize some m e a s u r e of d i s c r e p a n c y between
Go
a n d y o :y o A
=
Arg min !discrepancy (Go, y O) j.A (2)
The c h o i c e of t h e "distance" i s defined by t h e s t r u c t u r e s of y ( 2 ) . F r e q u e n t l y some c o n s t r a i n t s c a n b e imposed o n X in optimization problem (2).
Two main a p p r o a c h e s c a n b e e a s i l y t r a c e d in t h e c o r r e s p o n d i n g l i t e r a t u r e : t h e f i r s t i s d e t e r m i n i s t i c a n d closely r e l a t e d to c l a s s i c a l function approximation t h e o r y , while t h e s e c o n d i s b a s e d o n t h e assumption t h a t y ( z ) i s g e n e r a t e d by a random field ( s e e f o r i n s t a n c e Katkovnik (1985), Miccheli a n d Wahba (1981) a n d Ripley (1983)).
In both c a s e s , t h e most c r u c i a l assumptions are r e l a t e d to t h e "smoothness" of s u r f a c e y ( 2 ) . In t h e d e t e r m i n i s t i c c a s e , t h e s e assumptions usually c o n c e r n t h e d e r i v a t i v e s of y ( 2 ) . When y ( z ) i s g e n e r a t e d by a random field, assumptions con- c e r n i n g i t s c o r r e l a t i o n s t r u c t u r e are mainly used.
This p a p e r c o n c e r n s t h e latter case. The s q u a r e r i s k ( d i s c r e p a n c y m e a s u r e )
w h e r e p ( y o , y ) i s t h e d e n s i t y (assuming i t s e x i s t e n c e ) of t h e random v e c t o r ( y o , y ) E R n + I , will b e used as a c r i t e r i o n of t h e optimality of
Go.
I t i s known t h a t t h e l i n e a r e s t i m a t o r (1)-(3) i s t h e b e s t o n e only if y ( z ) i s a Gaussian random field.In o t h e r cases. i t i s r e f e r r e d to as t h e best L i n e a r e s t i m a t o r . O b s e r v a t i o n s w i t h k n o w n m e a n - c o v a r i a n c e s t r u c t u r e .
L e t u s assume t h a t :
( a ) for a n y z a n d z ' t h e mean m ( z ) = E [ y ( z ) ] a n d t h e c o v a r i a n c e C ( z ,z ')
=
E [ ( y ( 2 ) - ( z ) ) ( y (2')-m ( z ' ) ] are known. In matrix n o t a t i o n , t h i s means p a r t i c u l a r l y t h a tare given.
From definition (3) a n d (4), i t follows t h a t :
v2(X)
=
COO+
m t-
2 h T ( ~ 1 0 + m o m )+
XT(Cil+mm T)X.
(5)The minimization of (5) l e a d s to t h e best L i n e a r ( u n c o n s t r a i n e d ) e s t i m a t o r :
A
Xl
=
(Cll+rnrn T)-l (ClO+rnrnO) , ( 6 )T -1
u2(i1)
=
min u2(X)=
Coo+
m-
(Col+m Tko)(Cll+mm ) (Clo+m Tm 0).
(7)A
When l i n e a r c o n s t r a i n t s are imposed:
F X = L ,
w h e r e F i s a (k x n ) - m a t r i x a n d L i s a (k x1)-vector, t h e n t h e solution of X2 of (5) i s defined by t h e system:
where p c o r r e s p o n d s t o Lagrangian multiples.
The v a r i a n c e of t h e prognosis i s equal to:
v2(X2)
= v 2 ( i 1 ) +
+ T [ ~ T ( ~ l l + m m T ) - l ~ ] - l + , ( 1 0 )T -1
where .I. = F T ( c l l +mm ) (COl + m om ) -L
.
The l a s t t e r m in ( 1 0 ) i s positive d u e t o positive definiteness of t h e matrix F T ( c l l + m m T ) - l ~ .In o t h e r words, when demanding t h e fulfillment of some p r o p e r t i e s f o r ( t h e fulfillment of ( 8 ) ) , o n e s a c r i f i c e s i t s p r e c i s i o n in t h e s e n s e of (3):
v2(K1)
s
v 2 ( K 2 ).
( 1 1 )U n k n o w n m e a n s t r u c t u r e , k r i g i n g .
Formulae ( 6 ) , ( 7 ) a r e of t h e o r e t i c a l i n t e r e s t b u t f o r a p r a c t i t i o n e r t h e i r use- f ulness i s v e r y r e s t r i c t e d . Knowledge of both m ( z ) a n d C ( z ,z ') i s more o f t e n a n e x c e p t i o n t h a n a f r e q u e n t l y m e t situation. T h e r e f o r e , a number of a t t e m p t s t o con- s t r u c t e s t i m a t o r s t h a t d o not u s e m ( 2 ) a n d C ( z , z ' ) o r at l e a s t p a r t of t h i s informa- tion h a v e r e p e a t e d l y been u n d e r t a k e n . One of them i s k r i g i n g , w h e r e knowledge of m ( z ) i s not n e c e s s a r y . Unfortunately o n e h a s t o p a y f o r t h i s by t h e additional as- sumption:
(b) t h e mean of y ( z ) d o e s n o t depend upon z , e.g., m ( z )
=
m 0 , and by t h e fol- lowing c o n s t r a i n t imposed o n A :k
In terms of ( 8 ) , t h i s means t h a t F = I T i s t h e v e c t o r ( l x k ) with a l l elements equal to 1 a n d L = m o . Constraint ( 1 2 ) p r o v i d e s a n unbiasness of t h e c o r r e s p o n d i n g esti- m a t o r X T y :
From ( 9 ) a n d ( 1 2 ) i t follows t h a t t h e k r i g i n g e s t i m a t o r i s t h e solution t o t h e following l i n e a r system
which d o e s n o t include v e c t o r m
The r e s t r i c t e d n a t u r e of assumption (b) w a s evident at t h e y e r y beginning of kriging h i s t o r y a n d t h e so-called u n i v e r s a l k r i g i n g e s t i m a t o r Auk w a s p r o p o s e d (Huijbregts a n d Matheron ( 1 9 7 1 ) ) , which was based o n t h e assumption:
(b') t h e mean of y ( z ) c a n b e p r e s e n t e d in t h e form m ( z )
=
1 9 ~ f ( z ), where f ( z ) i s a v e c t o r ( k x l ) , k<
n , of a p r i o r i known functions. To o b t a i n t h e unbiased e s t i m a t o r , o n e h a s t o impose t h e following c o n s t r a i n t :E [ x ~ ~ ]
=
X T F ~ I Y=
f T ( z ) ~ , f o r a n y z E X ,w h e r e
F =
L P ( z 1 ) , f ( z 2 ) * .. .
# f ( ~ , ) lFrom (9) a n d (14) i s follows t h a t t h e u n i v e r s a l k r i g i n g e s t i m a t o r
Xuk
i s t h e solu- tion of t h e system:which similarly to (13) d o e s n o t include v e c t o r m
.
The e s t i m a t o r
iuk
i s unbiased if ( b ' ) i s fulfilled b u t v 2 ( i ) S v2(Xuk ).F o r t h e u n i v e r s a l k r i g i n g e s t i m a t o r , (10) c a n b e t r a n s f o r m e d to
v 2 ( i u t )
=
Coo -c ~ ~ c ~ ; ~ C ~ ~ +
+
~ ( z o ) - ~ c , i ~ c ~ o l ~ ( ~ , i ~ ~ ~ ) - ~ L P ( z o ) - ~ ~ ~ ~ ~ o ~
(16) Kriging i s m o r e a t t r a c t i v e t h e o r e t i c a l l y when ( a ) i s c h a n g e d i n t o t h e followingassumption (Huijbuegts a n d Matheron (1971)):
w h e r e i t i s assumed t h a t only t h e i n c r e m e n t s of t h e random field admit t h e f i r s t a n d s e c o n d o r d e r moments. From t h e o r y , i t i s known t h a t t h e r e are cases when t h e s e moments d o not e x i s t f o r y ( z ) itself while t h e y e x i s t f o r t h e c o r r e s p o n d i n g increments. Function 7 ( z , z ') i s c a l l e d a s e m i v a r i o g r a m a n d
2 7 ( z ,z')
=
[m ( z )-
m ( z ')12+
C ( z ,z )+
C ( z ' , z ') - 2 C ( z ' , z J ) ,if mean m ( z ), v a r i a n c e C ( z ,z ) a n d c o v a r i a n c e C ( z ,z') e x i s t . Usually in p r a c t i c e o n e believes in t h e i r e x i s t e n c e .
In applied s t u d i e s using kriging, a u t h o r s p r e f e r to work with t h e semivariogram 7 ( z ,z') i n s t e a d of t h e familiar c o v a r i a n c e C ( z ,z') although i t d o e s n o t make a c c e p t a n c e of t h e r e s u l t s a n y e a s i e r ( s e e t h e n e x t s e c t i o n ) .
To summarize t h i s s u b s e c t i o n , o n e c a n s a y t h a t k r i g i n g i s a p a r t i c u l a r c a s e of l i n e a r e s t i m a t i o n t h e o r y o r , to b e m o r e s p e c i f i c , some g e n e r a l i z a t i o n of t h e Wiener-Kolmogorov f i l t e r .
P r o b a b l y i t i s also worthwhile to n o t i c e t h a t a u n i v e r s a l k r i g i n g e s t i m a t o r c a n b e c o n s t r u c t e d in t h e framework of t h e g e n e r a l i z e d least s q u a r e method ( s e e sec- t i o n IV) when t h e model
i s c o n s i d e r e d u n d e r t h e assumption t h a t t h e c o v a r i a n c e s t r u c t u r e of t h e random error E ( Z ) i s known.
The simple k r i g i n g e s t i m a t o r c a n b e c o n s i d e r e d a l s o as a p a r t i c u l a r case of t h e m o v i n g a v e r a g e ( s e e , f o r i n s t a n c e , Katkovnik (1985)):
w h e r e weights Xi a r e defined by (3).
Clbservations w i t h u n k n o w n m e a n - c o v a r i a n c e s t r u c t u r e
Assumption ( a ) ( o r t h e t h e o r e t i c a l l y slightly milder assumption (a')) i s d r a s t i - cally r e s t r i c t i v e in p r a c t i c e f o r environmental analysis. In t h e publications r e l a t - ed t o kriging, t h e author_ could not find a n y a p p r o a c h o t h e r t h a n using t h e empiri- c a l estimates
6
( z ) a n d C ( z , z ') in p l a c e of m ( z ) and C ( z , ~ ' ) .P r o b a b l y t h i s recommendation i s worthwhile in engineering s c i e n c e s w h e r e t h e r e e x i s t s t h e possibility of r e p e a t i n g similar e x p e r i m e n t s a n d w h e r e o n e c a n u s e so-called l e a r n i n g samples t o r e s t o r e m ( z ) a n d C ( z , z ' ) and t h e n t o u s e krig- ing. But in t h i s case t h e simpler and well established technique c a n b e used. F o r instance, instead of ( I ) , (8) o n e c a n u s e t h e estimator:
$ = X o + X T y t h a t in c a s e of (3) i s defined by t h e following formulae:
and
If t h e l e a r n i c g sample i s sufficiently l a r g e a n d o n e c a n f o r g e t t h e uncertain- t i e s of
6
( z ) and C ( z ,z ') t h e n (1.7) is b e t t e r t h a n kriging.In environmental a p p l i c a t i o n s a n d by t h e way, in g e o s t a t i s t i c s ( t h e homeland of t h e kriging method), i t i s more r e a l i s t i c to u s e t h e t e r m o b s e r v a t i o n s t h a n ez- p e r i m e n t s . T h e r e f o r e , t h e problem of a good learning sample becomes h e r e espe- cially a c u t e and usually t h e r e i s no hope t h a t t h e e r r o r s of estimates $ ( z ) and C ( z ,z ') c a n b e neglected.
In t h i s c a s e , t h e a v e r a g e in formulae (7) and (10) (where all a p r i o r i unknown values are s u b s t i t u t e d by t h e i r estimates) h a s t h e conditional c h a r a c t e r :
E [ ( Q --y)2/ learning sample]
=
v2(j;)Of c o u r s e t h e t o t a l variance, taking into a c c o u n t t h e randomness of l e a r n i n g sample, will b e g r e a t e r t h a n v 2 ( ~ ) defined e i t h e r by (7) or (10). To some e x t e n t t h e y c a n b e considered as e s t i m a t e s of the Low b o u n d s for t h e c o r r e s p o n d i n g to- t a l v a r i a n c e s .
The plight becomes worse when t h e same set of d a t a y
=
( y 1 , y 2 ,. . .
, yn ) i s used both f o r t h e estimation of C,,, C,, a n d kriging. In t h a t c a s e , t h e a u t h o r h a s failed to find in t h e l i t e r a t u r e any s e r i o u s t h e o r e t i c a l r e s u l t s on t h e evaluation of t h e b i a s a n d t h e v a r i a n c e of t h e kriging estimators.A sensitivity analysis examining how t h e estimate
Gut =
;&y will c h a n g e f o r given small p e r t u r b a t i o n s in Cll and Clo was d o n e by Warnes (1986). But i t i s a nu- m e r i c a l analysis of t h e s t a b i l i t y of system (15), r a t h e r t h a n a s t a t i s t i c a l analysis of t h e c o r r e s p o n d i n g e s t i m a t o r . which, f o r instance, must d e t e r m i n e how much t h e v a r i a n c e ofcut
will c h a n g e u n d e r small random p e r t u r b a t i o n s of t h e aforemen- tioned m a t r i c e s .To conclude t h i s s e c t i o n , o n e c a n s a y t h a t in t h e c a s e of a p r i o r i unknown Col a n d Cll, t h e k r i g i n g a p p r o a c h g i v e s some h i n t s o n h o w to choose w e i g h t s for t h e
Linear e s t i m a t o r
x~~
b u t i t g i v e s n o s e r i o u s g u a r a n t e e of t h e i r o p t i m a l i t y .III.
ExamplesIn a l l t h r e e examples presented h e r e , attention will b e focussed on t h e most sensitive a s p e c t of kriging : fitting of y ( z ,z ') o r C ( z ,z ') , leaving aside o t h e r a s p e c t s of t h e problem.
1. D i s t r i b u t i o n of r e s i d u a l c o n t a m i n a t i o n from a t m o s p h e r i c n u c l e a r t e s t s (Barnes (1980)).
Thi e ample i l l u s t r a t e s how t h e sernivariogram w a s estimated using d a t a on isotope
FZ4k
f o r one ground z e r o site. Smallboy. The d a t a w a s e x t r a c t e d from a more extensive study r e l a t e d t o t h e analysis of residual contamination a f t e r atmos- p h e r i c nuclear tests. The final objective w a s t o evaluate t h e contours of 241 activity. All experimental and geographical details c a n b e found in (Barnes (198&and in t h e r e f e r e n c e s given t h e r e . In t h e cited r e p o r t it w a s assumed t h a t t h e observed value y ( z ) i s a random spatially weak stationary field, e.g., E [ y ( z ) ]
=
m and y ( z , z ')=
y ( h ), where h=
d ( z -z ') '(2 -z ')This assumption i s essential because f i r s t l y it allows one t o estimate y ( z , z 8 ) and secondly i t makes i t possible to u s e simple kriging (m ( z ) = m , s e e assumption (b) from t h e previous section). The consistency of t h e assumption with r e a l i t y can b e evaluated, at l e a s t p a r t l y , with t h e help of Figure 1. Here, t h e "observed"
values ?(A) are plotted by numbers which correspond to t h e directions presented at t h e top-right c o r n e r of Figure 1 .
The "observations" were calculated by t h e following p r o c e d u r e (Barnes (1980)): "All points within a small angle of t r u e east-west of a given point are put in t h e "east-west" class and points t h a t a r e "approximately" h away of a given point are put in t h e distance h class. The size of small angle and t h e closeness of t h e distance approximation c a n b e controlled by t h e u s e r t o r e d u c e t h e e r r o r intro- duced by t h e s e approximations".
The dashed line in Figure 1 c o r r e s p o n d s t o t h e final approximation of y ( h ) and i t i s clear t h a t both anisotropy and possible d r i f t were ignored. Barnes did not mention how ?(A) w a s fitted t o t h e d a t a , but what i s evident i s t h a t i t does not follow t h e points at a l l , e x c e p t maybe in t h e interval 0
c
h S 1000ft.From t h e definition of t h e semivariogram, it follows t h a t in t h e c a s e when t h e d r i f t of m ( z ) can b e ignored and y ( z ) i s weakly spatially s t a t i o n a r y , then
y ( h )
=
C'-
C ( h ) , where C=
C ( z , z ) and C ( h )=
C ( z , z 8 ) .I t seems t h a t f o r t h e physical phenomenon (distribution of contamination) con- sidered by Barnes, C ( h ) 2 0 f o r a l l h and lim C ( h )
=
0 (implicitly t h i s i sh -+-
assumed, see p.6 of cited p a p e r , where various t y p e s of y ( h ) are considered).
T h e r e f o r e , f o r t h e u p p e r bound of y ( h ), one h a s
S?P
7 ( h=
C ,e.g., t h e plateau of y ( h ) ( o r "sill") h a s to b e equal to t h e v a r i a n c e C of observed values. From Figure 1 it follows t h a t C N 3400. Unfortunately, t h e s c a l e of t h e v e r t i c a l a x i s w a s not a c c u r a t e l y defined. Does it correspond t o y ( h ) o r t o 2 y ( h ) ? S o i t could b e t h a t C a! 1700, but this would not change t h e situation, i.e., where one cannot b e s u r e t h a t t h e inaccuracy of ?(A) defined by t h e dashed line in Fig- u r e 1 is not less than 30%. How c a n i t violate t h e r e s u l t s of t h e kriging p r o c e d u r e ? I s t h e r e any hope f o r i t s optimality? O r maybe i t i s a v e r y indirect way to con- s t r u c t a mediocre moving a v e r a g e estimator with v e r y r e s t r i c t i v e intermediate assumptions?
0 500 loo0
1m m 2.m
m,D l STANCE (feet) h
Figure 1: Pointwise e s t i m a t e s for t h e semivariogram 7 ( h ) , B a r n e s (1980).
P r o b a b l y i t i s b e t t e r to u s e t h e moving a v e r a g e d i r e c t l y b e c a u s e i t n e e d s simpler computing a n d more s t r a i g h t f o r w a r d a n d e x p l i c i t assumptions a b o u t t h e na- t u r e of y (3: ).
2. K r i g i n g a n a l y s i s of the r e g i o n d p a t t e r n s of t h e chemical c o n s t i t u e n t s of p r e c i p i t a t i o n .
In t h e r e p o r t o n "Statistical analysis of p r e c i p i t a t i o n chemistry measure- ments" by Endlich et al, (1986), Ch. 4, t h e kriging a p p r o a c h was a p p l i e d t o interpo- late t h e y e a r l y medians of t h e daily l a b o r a t o r y
pH
a n d a n a l y t e c o n c e n t r a t i o n s at e a c h s i t e a n d f o r t h e y e a r l y t o t a l deposition of H a n d o t h e r a n a l y t e s . To p e r f o r m k r i g i n g t h e a u t h o r s s e p a r a t e d 'long-term s p a t i a l a n d temporal t r e n d s from y e a r - t o - y e a r fluctuations". The q u a d r a t i c s p a t i a l polynomial plus t h e l i n e a r time t r e n d w e r e used t o a p p r o x i m a t e t h e logarithms of t h e s e t r e n d s . The f i t t i n g technique was t h e o r d i n a r y l e a s t s q u a r e method. The r e s i d u a l s (observed v a l u e s minus t h e least s q u a r e estimates) were t a k e n as inputs f o r kriging. T h e i r v a r i a t i o n s f r o m year-to-year were a s s u m e d to be i n d e p e n d e n t , a n d t h e c o v a r i a t i o n s t r u c t u r e was" a s s u m e d to be both s t a t i o n a r y (independent of Location) a n d i s o t r o p i c ( i n - a e p e n d e n t of direction)"
The semivariograms were evaluated by a l i n e a r model:
y ( h )
=
I9,+
b 2 h.
Fitting w a s done by t h e weighted l e a s t s q u a r e method with weights p r o p o r t i o n a l t o t h e number of s i t e - p a i r s and inversely p r o p o r t i o n a l t o d i s t a n c e h . Typical exam- p l e s of t h i s fitting a r e p r e s e n t e d in Figure 2 . I t i s difficult t o e v a l u a t e t h e good- n e s s of fitting b e c a u s e t h e i n v e r s e values of weights were not plotted. However, i t seems t h a t t h e f i t i s n o t v e r y good. Besides t h e l i n e a r approximation d o e s n o t re- f l e c t t h e f a c t t h a t t h e o b s e r v a t i o n s from t h e mutually remote s i t e s h a v e n o t t o b e c o r r e l a t e d and lim y ( h )
=
c o n s t , i . e . , y(h ), h a s to a p p r o a c h s a t u r a t i o n (compareh +-
with t h e p r e v i o u s example and with t h e definition of a semivariogram).
T h e r e are some o t h e r points t h a t c a n b e criticized:
(a) Use of t h e l e a s t s q u a r e method to remove t r e n d s and subsequent application t o kriging i s not c o n s i s t e n t with t h e o r y : o n e h a s to u s e at l e a s t universal kriging in- s t e a d of t h e s e t w o s t e p s . Only in t h i s case are t h e estimates optimal in t h e s e n s e of (3) o r at l e a s t approximately optimal if t h e estimates f o r C ( h ) and subsequently f o r C1, a n d Cll a r e sufficiently p r e c i s e . I t seems t h a t t h e use of c o v a r i a n c e s (if t h e y e x i s t ) i s m o r e convenient both from t h e t h e o r e t i c a l and t h e computational points of view. In what follows, C ( h ) o r C ( z , z '), o r C1, ,Cll will b e used without comments.
With r e s p e c t t o r e s i d u a l s u i , t h e a u t h o r s of t h e r e p o r t did n o t notice t h a t t h e i r variance-covariance matrix i s singular. T h e r e f o r e a l l "kriging" formulae ( s e e , f o r instance, (9), (15), (16)) c a n n o t be % s e d d i r e c t l y ( t h e solution will n o t b e unique). P r o b a b l y , t h e use of some estimate C l l instead of Cll will r e g u l a r i z e com- putations. But t h i s "regularization" simultaneously means t h e loss of optimality of t h e p r o g n o s e s ( z ) d u e t o t h e p o o r estimation of Cll.
The singularity of C l l c a n b e easily verified.
Assume t h a t y = F T $ + & ( o r ~i
= f
T(zi )*+ci w h e r e F=
(f(zi) ,. . .
, f ( z n ) ) . Then t h e l e a s t s q u a r e e s t i m a t o r3
is defined by t h e f o r -mula ( s e e , f o r i n s t a n c e , Anderson, 1971):
-
I9 = ( w T ) - I F y , t h e v e c t o r of r e s i d u a l s e q u a l su
=
y- F T S
= [ I - F ~ ( F F ~ ) - ~ F ] E and t h e v a r i a n c e - c o v a r i a n c e matrix of t h e r e s i d u a l sEll =
[I-F~(IT~)-~F]C,~[I-F~(IF~)-~F].The p ~ o j e c t i o n matrix I
-
F T ( I T T ) - l F h a s r a n k (n-
m ) , w h e r e m i s t h e di- mension of I9. T h e r e f o r e r a n k CI1=
n-
m , i.e. Cll i s singular (I
Cll1 =
0).(b) The method of semivariogram estimation c a n b e improved (from t h e s t a t i s t i c a l point of view) almost without a n i n c r e a s e of computations.
F o r t h e s a k e of simplicity l e t u s assume t h a t o b s e r v e d values (i.e., o r d i n a t e s in Figure 2 ) are normally d i s t r i b u t e d . Then (see, for instance, S e b e r , 1977, Ch.
14):
F i g u r e 2: Empirical a n d f i t t e d semivariogram functions, Endlich et a1 (1986).
where h i s t h e distance between observation points. In t h e case of o t h e r distribu- tions, t h i s formula is more complicated.
The method of t h e i t e r a t i v e l e a s t squares:
n
19,
+, =
Arg min C st [yt -y(ht ,19)12/ y2(ht , 9 , ) ,4 =
lim IJ, ,d i =1 s
--
gives asymptotically ( n --) optimal estimates f o r I9 and y ( h ,19) , i.e., t h e method minimizes v a r (19) and v a r y[(h ,6)], see Malyutov (1982). In (18) ri i s t h e number of observations f o r e v e r y hi , I9 s t a n d s f o r unknown p a r a m e t e r s . In p r a c t i c e usu- ally t h e i t e r a t i v e p r o c e d u r e (18) i s terminated a f t e r 3,4 s t e p s .
(c) F o r all f o u r f i t t e d lines in Figure 2, t h e i n t e r c e p t ( i.e., i s significantly g r e a t e r 0 . I t means t h a t t h e r e e x i s t s a so-called n u g g e t e f f e c t (see, f o r instance, Gilbert and Simpson (1984)), i.e., a discontinuity of t h e c o v a r i a n c e function C ( z ,z'). This could probably o c c u r in geostatistics when o n e analyses t h e deposi- tion of some o r e minerals, but in t h e analysis of pollutants in fluids or atmospheric contamination, i t seems unreal. This i s also confirmed by t h e e x i s t e n c e in e a c h p a r t of Figure 2 of observation points located close t o t h e origin. Presumably, only observations with h r 3 f 4 still satisfy t h e kriging assumptions a n d more distant observations are e i t h e r violated by t r e n d s and anisotropy, o r t h e c o v a r i a n c e func- tion i s negative f o r h r 4 . The r e s u l t s p r e s e n t e d in Figure 3 ( t h e model
~ ( h )
=
+19~z')/ ( l + ' l ~ ~ z ~ ) , z =e -h -1, w a s f i t t e d t o t h e d a t a with t h e help of i t e r a t i v e u s e of BMDP AR program, 1983) are a good confirmation of t h i s a s s e r t i o n .The above c r i t i q u e leads to t h e same conclusion as in t h e previous example.
3. S u l f u r d e p o s i t i o n model e v a l u a t i o n .
One of t h e main goals of t h e r e p o r t by Clark et a1 (1986), w a s t h e comparison of t h e s e v e r a l models c u r r e n t l y used f o r computing of s u l f u r deposition in North America. Roughly, t h e s t r a t e g y consists of t h e following s t e p s :
-
observations from i r r e g u l a r l y located observation s t a t i o n s are used t o esti- mate values at r e g u l a r g r i d points with t h e help of kriging,-
model predictions o v e r some g r i d are used to obtain gridpoint values with t h e help of kriging,-
both sets of r e s u l t s are compared by p r o c e d u r e s mainly based on methods of mathematical statistics.In t h i s a p p r o a c h , f o r instance, "a confidence interval" f o r t h e d i f f e r e n c e in o b s e r v e d and p r e d i c t e d values w a s c o n s t r u c t e d as follows:
where
vPmd
=
kriging estimate f o r model prediction at t h e g r i d point, yobs=
kriging estimate f o r t h e observed d a t a at t h e g r i d point, V a r (.)=
kriging v a r i a n c e estimate.If t h e deviation of k r i g e d values from t r u e values t h a t they estimate h a v e normal (Gaussian) distribution, t h e n (19) defines 95% confidence interval. S e e Clark et a1 (1986), p. 5.12.
(b) Total paipitation pr y e n
u-
1.6 4 t
.6 4 0 t
? 0 0
- 0
d
-1
0 0 1.2 4 0 t- P 0 a-
?
t
0
0 .80 4
. 3 4 0
0 0
, 4 0 4 : 0 0 0 0 0
- P
-0 -P
0. t t -0
... ... . . . ... ... ... ...
. + . . t . . . t . ... . t . . + . .+.. . + . . . . . . .
. . . ... .... ... ... .... t + 1 + + t t
2 . 6 . 10 I4 2 . 6 . 10 14
. t ... 4 ... t ... t ... t ... t . . .
6 . 10 I4
DISTANCE 1 1 YAP UNIT
-
170 K Y )~~~l u* -itson pa vem lhDm IhbomoV pH) 4
3 5 ;
f
tt - P
.21 4 0 t
0
P
pt
q
-- $ . l I I t 4
.
P 0 02
-0- 0 0 0 0
-0 # P P * P P i P P
i:
t 0
. 0 9 4 0 0 t
-P
4 0-
.... ... 4 ... t . . . .... 1 ... t ... t ,,...,. t ..,.. . . t . , ...
2 . 6. 10 14
DISTANCE I 1 MA? UNIT
-
170 KMIF i g u r e 3: Fitted semivariogram y ( h )
=
(.IPlz +.IP2z ') / (1 +I@ '), z =e -1, v e r t i c a l l i n e s s t a n d f o r t h e s t a n d a r d e r r o r of p r o g n o s e s .From t h e p r e v i o u s considerations ( s e e t h e concluding p a r t of section I1 a n d t h e f i r s t two examples) i t i s c l e a r t h a t t h e v a r i a n c e s computed by substitution of t r u e values of C ( z ) a n d C ( ( z , z ') ( o r y ( z , z ')) by t h e i r estimates (which are usually v e r y p o o r ) h a v e a v e r y remote r e l a t i o n with r e a l i t y . T h e r e f o r e , o n e c a n u s e (19), keeping in mind a number of "ifs" d u e t o t h e violation of assumptions a b o u t "pure kriging" a n d d u e to t h e assumption o n t h e normality of d i s t r i b u t i o n s of yPd a n d Yobs
Together with t h i s t e c h n i c a l r e m a r k t h e r e i s o n e more g e n e r a l comment. The s t r a t e g y of model comparison with d a t a c a n b e d e s c r i b e d by t h e scheme p r e s e n t e d in Figure 4. I t i s c l e a r from t h i s scheme t h a t one compares d i s t o r t e d images of t w o objects. In p r i n c i p l e t h e distortion d u e to p r o j e c t i o n I c a n b e v e r y small be- c a u s e a u s e r c a n continue computing until h i s models give good r e s u l t s o n any given g r i d . Then t h e whole p r o c e d u r e of comparison c a n b e d e s c r i b e d as comparison of t h e given set of models with o n e v e r y simple model which i s defined by t h e kriging method assumption ( s e e Section 11). Implicitly t h e a u t h o r s believe t h a t t h i s model i s b e t t e r t h a n a n y o t h e r c o n s i d e r e d in t h e i r r e p o r t .
P r o b a b l y t h e a p p r o a c h schematically p r e s e n t e d in Figure 4 would b e more f r u i t f u l .
IV.
Some Alternatives to KrigingG e n e r a l i z e d l e a s t s q u a r e s (g.1.s.) e s t i m a t o r . In t h i s subsection t h e links between t h e kriging a p p r o a c h a n d t h e old-fashion l e a s t s q u a r e s e s t i m a t o r s will b e illuminated.
Let similarly to (b') from s e c t i o n I1
y f = J T ( ~ f )+
+
cf , i = l,n , where f ( z ) i s t h e v e c t o r of known b a s i c functions,E [ q ]
=
0 , E [ c f E~'1 =
Cll,ffd , o r in t h e matrix p r e s e n t a t i o ny = F ~ + + E , E [ E E ~ ] = c ~ ~
.
(21)To estimate t h e value of t h e r e s p o n s e y ( z ) at a given point z , o n e c a n follow t h e t w o s t e p s p r o c e d u r e :
-
compute t h e p a r a m e t e r s estimates$,
-
p r e d i c t y ( z ) using t h e l i n e a r estimatorG(z) =
f T ( z ) 4 + A & ,where u i s t h e v e c t o r of r e s i d u a l s , i.e. u
-=
y-
F ~ + .Vector Xo i s defined as a solution of t h e following optimization problem A , = A r g m p E [ u ( z ) - f T ( z ) a - X ~ ( ~ - F ~ $ ) ] ~ ,
F i g u r e 4: Comparison of model o u t p u t s with d a t a .
I t i s known (see, f o r i n s t a n c e , S e b e r , 1977) t h a t t h e g e n e r a l i z e d l e a s t s q u a r e e s t i m a t o r (which i s t h e b e s t l i n e a r unbiased e s t i m a t o r )
8 =
Arg min ( y - F ~ I ~ ) ~ c G ' ( y -FT 19)d
c a n b e c a l c u l a t e d implicitly:
5
= ~ - l ~ c , i l y , M =FC,-~FT.
Formally (22) coincides with (3) a n d at f i r s t g l a n c e all c o n s i d e r a t i o n s f r o m s e c t i o n I1 c a n b e applied to (23). Nevertheless. t h e r e e x i s t some s p e c i f i c f e a t u r e s of (22):
-
I t i s not n e c e s s a r y to impose a n y c o n s t r a i n t s o n h t o g e t a n unbiased estime- tor, b e c a u s e E[u (z)] = 0 a n d E [ u ]=
0, d u e to t h e unbiased n a t u r e of I9, EL91=
9 , .-
The c o v a r i a n c e s t r u c t u r e of v e c t o r u is s i n g u l a r , i.e. IDl11 =
0, w h e r e Dll=
E [UU T]. T h e r e f o r e t h e majority of formulae from s e c t i o n I1 c a n n o t b e used. The s t r a i g h t f o r w a r d calculations give:Unfortunately o n e c a n n o t u s e (I?), w h e r e h = 0fi1lI1,, b e c a u s e Dll i s s i n g u l a r ( s e e a l s o comments in Example 2 ) a n d o n e must find a way to s o l v e t h e s i n g u l a r l i n e a r system:
(c,,-F~M-~F)A =
c,,
- F ~ M - ~ F C , - , ~ C , , (25) One of t h e simplest solutions of (25) is:ho
= c,~~c,,
,(CI1 assumed t o b e r e g u l a r ) .
I t h a s to b e s t r e s s e d t h a t Cll i s t h e v a r i a n c e - c o v a r i a n c e m a t r i z of t h e vec- t o r E b u t n o t t h e v e c t o r of r e s i d u a l s u! sombining (24) a n d (26) o n e finds t h a t t h e v a r i a n c e of t h e e s t i m a t o r y^(x )
=
f T ( z ) 9+
C, l ~ $ u i s e q u a l toIf t h e e s t i m a t o r t ( z ) = f T(x
)a
i s used, t h e nE [ y ( x ) - F ( z ) I
=
C, + ~ ~ l z ) ~ - ~ 3 ( 2 : ).
The same r e s u l t will b e o b t a i n e d if Ci, = 0 ( o b s e r v a t i o n at point x i s u n c o r r e l a t e d with o t h e r o b s e r v a t i o n s ) .
E x p r e s s i o n (27) coincides with (16). i.e. t h e u n i v e r s a l k r i g i n g c a n be con- s i d e r e d as a p a r t i c u l a r a p p l i c a t i o n
03
t h e g e n e r a l i z e d l e a s t s q u a r e method.Application to t h e least s q u a r e method allows o n e to trace a common t a k e in a number of applied s t u d i e s ( s e e examples) viz., t h e e s t i m a t e s of Dll a n d Dzl are used to c a l c u l a t e ho ( o r &); see a l s o (24).
In s p i t e of t h e t h e o r e t i c a l c l a r i t y of t h e generalized l e a s t s q u a r e method, i t s applicability to r e a l e m p i r i c a l s i t u a t i o n s i s v e r y p r o b l e m a t i c b e c a u s e of t h e neces- s i t y to know m a t r i c e s Cll a n d C1, a n d o n e c a n r e p e a t h e r e a l l c o n s i d e r a t i o n s f r o m s e c t i o n I1 r e l a t e d t o t h e case with a n unknown mean-covariance s t r u c t u r e .
Moving l e a s t s q u a r e s e s t i m a t o r . In t h e majority of a p p l i c a t i o n s ( s e e s e c t i o n 111) t h e c o r r e l a t e d o b s e r v a t i o n a l errors are used t o simulate t h e deviations of real p r o c e s s e s from a comparatively simple t r e n d approximation. Example 2 i s a good example. Unfortunately, t h e model (i.e., deviations which are random values with c o v a r i a n c e s t r u c t u r e C, , C1,, Cli) c o n t a i n s t o o many unknown p a r a m e t e r s to p r o - vide good prognoses. The u s e of kriging-related a p p r o a c h e s w h e r e C,, C1,, CI1 are s u b s t i t u t e d by t h e i r r o u g h e s t i m a t e s (and sometimes wrongly c o n s t r u c t e d t h e o r e t i c a l l y ) misleads r e a d e r s (and p r o b a b l y a u t h o r s a l s o ) in t h e evaluation of t h e p r e c i s i o n of prognoses.
I t seems t h a t in t h e examples c o n s i d e r e d , t h e moving a v e r a g e o r i t s slightly modified v e r s i o n - m o v i n g Least s q u a r e s e s t i m a t o r c a n give b e t t e r r e s u l t s a n d c l e a r e r a n d more d i r e c t understanding of t h e bounds of admissibility of t h e as- sump tions used.
Let z b e a point w h e r e t h e p r o g n o s i s h a s t o b e made, a n d zi
=
z+
ui b e loca- tions of o b s e r v a t i o n points. Assume t h a t :( a ) I n t h e vicinity of point z t h e following approximation i s valid
w h e r e y i i s t h e r e s u l t of o b s e r v a t i o n at point xi
=
z+
ui,'lPO(z) a n d $ ( z ) a r e p a r a m e t e r s t o b e estimated, E~ ( z ) i s t h e o b s e r v a t i o n (and approximation) er- ror,p
( u , Z ) i s a v e c t o r of given functions vanishing when u -, 0 .The e s t i m a t o r c a n b e defined a s
n
f $ o ( z ) , $ ( ~ ) j
=
Arg min X ( u i , z ) [ y i -'lPO-'lPT~(~i , 2 ) l 2.
G ( z )=
(29) dose i =1Function X(u ,z ) h a s to r e f l e c t t h e confidence in using approximation (28) a t point u. Normally X(u , Z ) i s a unimodal function and
h(O,z)
=
max X(u , z ) , lim X(u , z )=
0 ,U u +-
Using d i f f e r e n t h ( u , z ) o n e c a n easily v a r y t h e smoothness
5
( z ) . Due to (30), t h e a p p r o a c h e s similar t o (29) are f r e q u e n t l y a d d r e s s e d by t h e d i s t a n c e - w e i g h t e d Least s q u a r e s m e t h o d ( s e e Ripley, 1981, Ch. 37). In t h e c a s e when t h e r e i s n o p r i - or "physical" information a b o u t y ( z ) ( o r j ' ( u , z ) ) , o n e c a n c o n s i d e r (28) as t h e Taylor approximation of t h e r e s p o n s e y ( z ). F o r t h e second o r d e r Taylor a p p r o x i - mation, o n e will h a v e in t h e two-dimensional case:w h e r e ei ( z ) is t h e r e m a i n d e r t e r m at point z
-
ut.
F o r a sufficiently d e n s e o b s e r - vation network, approximation (31) usually s e r v e s r e l i a b l y . In more g e n e r a l case o n e c a n u s e any r e a s o n a b l e ( b e t t e r if s u p p o r t e d by some physical c o n s i d e r a t i o n s ) f ( u , z ) vanishing when u -, 0 .The s t a n d a r d l e a s t s q u a r e s technique provides simple algorithms a n d formulae f o r t h e calculation of y^(z)
=
Go(") ( s e e , f o r instance, Golub and Van Loan (1983), S e b e r (1977)). F o r t h e o r e t i c a l analysis, i t i s convenient t o u s e t h e following f o r - mulae:where
The e s t i m a t o r (29)-(32) c a n b e considered as some approximation scheme which c a n b e used i n both deterministic a n d s t o c h a s t i c a p p r o a c h e s . Usually in t h e s t o c h a s t i c case, i t i s e a s i e r t o e v a l u a t e t h e d i s c r e p a n c y between ( z ) and t h e t r u e value y (I), of c o u r s e paying f o r t h a t by t h e additional a n d p r a c t i c a l l y nonverified assumption:
(b) The o b s e r v a t i o n e r r o r s e i ( z ) are random values with E [ c i ( x ) ] = 0, ECci(z)cf (z)],
=
x-'(u~ , z ) b i i .If (b) holds a n d f ( 0 , z )
=
0 (it i s q u i t e a usual case, c o m p a r e with (31)) then:One h a s t o notice t h a t t h e o b s e r v e d value y
=
y ( z i )+
e i ( z ) . In many appli- c a t i o n s i t i s r e a s o n a b l e t o c h o o s e x-'(u , x )=
a )r+
b ( u , I ) , w h e r e cr2 i s t h e v a r i - a n c e of a n o b s e r v a t i o n e r r o r , 6 ( u , x ) comprises local s t o c h a s t i c fluctuations and d ( 0 , z )=
0 , lim 6 ( u,z) =
=. The estimation scheme (29)-(33) c a n b e generalizedu +-
in t h e case of c o r r e l a t e d o b s e r v a t i o n s . Application to t h i s case seems n o t t o b e useful b e c a u s e d u e to t h e local c h a r a c t e r of (28), t h e model a l r e a d y t a k e s into ac- count local t e n d e n c i e s a n d c h a n g e s while in t h e kriging ( o r similar) a p p r o a c h e s t h e y are handled via t h e c o r r e l a t i o n s t r u c t u r e of a n o b s e r v e d field.
The cautious r e a d e r will notice t h a t proposed e s t i m a t o r demands a r a t h e r tedious calculation necessitating inversion of matrix M ( z ) f o r e v e r y z t a k e n into consideration.
A t f i r s t glance, o n e c a n easily avoid t h i s by using to models d e s c r i b i n g t h e ob- s e r v e d field in t h e vicinity of f i z e d point x O:
a n d using subsequently t h e simplified v e r s i o n (Clltii.
=
dii. h-'(ui )) of t h e t e c h - nique discussed in t h e p r e v i o u s s u b s e c t i o n , whenI t i s c l e a r t h a t f o r a l l z , w h e r e < ( z ) h a s to b e c a l c u l a t e d , o n e h a s to s o l v e t h e l e a s t s q u a r e s problem only o n c e . Unfortunately, i t i s n e c e s s a r y to pay f o r t h i s simplification (which d o e s n o t seem to b e v e r y c r u c i a l in o u r computerized a g e ) by:
( a ) E s t i m a t o r (32) i s smooth (continuous, d i f f e r e n t i a b l e ) if f u n c t i o n s A(u , z ) a n d f ( u , z ) are smooth. E s t i m a t o r (35) will "jump" when z o will b e c h a n g e d ( i t i s c h a n g e d d i s c r e t e l y ) . To avoid discontinuities, o n e n e e d s to a d d r e s s to t h e least s q u a r e method merging t o g e t h e r p r o g n o s e s b a s e d o n d i f f e r e n t zo. But t h e merging p r o c e d u r e removei t h e simplicity of (35).
(b) Model (28) p r o v i d e s t h e b e s t approximation at point z , which i s of i n t e r e s t , while (34) is o r i e n t e d to some fixed point z o which c a n b e q u i t e remote f r o m t h e moving z
.
In conclusion, i t i s worthwhile to n o t e t h a t in t h e majority of a p p l i e d s t u d i e s in fluid flows, in meteorology or t h e a t m o s p h e r i c pollution s t u d i e s , in c o n t r a s t to geo- logical a p p l i c a t i o n s , o n e h a s t e m p o r a l as well as s p a t i a l information: see Example 2 . The methods d i s c u s s e d in t h i s s e c t i o n c a n easily i n c o r p o r a t e t e m p o r a l d a t a e x - plicitly expanding m a t r i x
F
by adding a time dimension in g e n e r a l i z e d l e a s t s q u a r e s c a s e o r in (29)-(32) t h e summation h a s to b e t a k e n o v e r s p a c e a n d time. In t h e k r i g i n g a p p r o a c h , t e m p o r a l information c a n b e implicitly used t h r o u g h improve- ment of t h e e s t i m a t e s f o r Cll a n d C1,.References
Akima, H. (1975) Comments o n "Optimal c o n t o u r mapping using u n i v e r s a l kriging" by R i c a r d o A. Olea J o u r n a l o f G e o p h y s i c a 1 R e s e a r c h , 80, pp. 832-836.
Armstrong, M . (1984) P r o b l e m s with u n i v e r s a l kriging. M a t h e m a t i c a l Geology, 16, p p . 101-108.
B a r n e s , M.G. (1980) The u s e of k r i g i n g f o r estimating t h e s p a t i a l d i s t r i b u t i o n of ta- dionuclides a n d o t h e r s p a t i a l phenomena. Battelle Memorial I n s t i t u t e , P a c i f i c Northwest L a b o r a t o r y , Richland, Washington, PNL-SA-9051, p p . 20.
Bell, G.D. a n d R e e v e s , M. (1979) Kriging a n d g e o s t a t i s t i c s : a r e v i e w of t h e l i t e r a - t u r e a v a i l a b l e in English, P r o c . Australas. Ins. Min. Metal. N o . 269, pp. 17-27.
BMDP (1983) Biomedical Computer P r o g r a m s , University of California P r e s s . C l a r k , T.L., Voldner, E.C., Olson, M.P., Seilkop, S.K. a n d Alvo, M. (1986) I n t e r n a -
t i o n a l s u l f u r deposition model evaluation (ISDME), R e p o r t .
Dennis, R.L. a n d Seilkop, S.K. (1986) The u s e of s p a t i a l p a t t e r n s a n d t h e i r u n c e r - t a i n t y e s t i m a t e s in t h e m o d e l s evaluation p r o c e s s . AMS/APCA Conf., p p . x x x . D e r Megreditchan, G. (1979) Optimization d e s r d s e a u x d ' o b s e r v a t i o n d e s champs
m6t6urologiques, L a Meteorologie, 17, p p . 51-66.
Endlich, R.M., Eynon, B.P., F e r e k , R.J., Valdes, A.D. a n d Maxwell, C. (1986) S t a t i s t i - c a l Analysis of P r e c i p i t a t i o n Chemistry Measurements Over t h e Ecosystem Un- i t e d S t a t e s , UAPSP-112, EPRI, P a l o Alto, California.
Finkelstein, P.L. and Seilkop, S.K. (1981) Interpolation error and t h e s p a t i a l v a r i - ability of a c i d p r e c i p i t a t i o n . P r o c . of t h e 7 t h Conference on Probability a n d S t a t i s t i c s in Atmospheric S c i e n c e s of AMS, AMS, Boston, pp. 206-212.
G i l b e r t , R.O. a n d Simpson, J.C. (1984) Kriging f o r estimating s p a t i a l p a t t e r n of con- taminants: potential and problems. E n v i r o n m e n t Monitoring a n d Assess- m e n t , 9, pp. 113-135.
Golub, G.H. and Van Loan, Ch. F. (1983) Matrix computations. The J o h n s Hopkins University Press, Baltimore, pp. 476.
Huijbregts, C. and Matheron, G . (1971) Universal kriging (an optimal method f o r estimating and contouring in t r e n d s u r f a c e analysis). Decision making in t h e mineral industry. Can. I n s t . M i n . Met. *ec., Vol. 12, pp. 159-169.
J o u r n e l , A.G. and Hui j b r e g t s , Ch. J . (1978) Mining geostatistics. Academic P r e s s . NY, pp. xxx.
Katkovnik, V.Y. (1985) On p a r a m e t r i c identification and d a t a smoothing. Nauka, Moscow, pp. 336.
Malyutov, M.B. (1982) Asymptotical p r o p e r t i e s and applications of t h e IRGINA- e s t i m a t o r of p a r a m e t e r s of generalized r e g r e s s i o n model. In "Stochastic p r o c e s s e s and applications", Moscow, pp. 144-158.
McBratney, A.B. a n d Webster, R. (1981) The design of optimal sampling schemes f o r l o c a l estimation a n d mapping of regionalized variables-11. Computers a n d Ceosciences, 7 , pp. 335-365.
Miccheli, C.A. a n d Wahba, G . (1981) Design problems f o r optimal s u r f a c e interpola- tion, Approximation Theory and Applications, Academic P r e s s , NY, pp.
329-349.
Ripley,
B.D.
(1981) S p a t i a l S t a t i s t i c s , Wiley, N Y , pp. 252.S e b e r , G.A.F. (1977) Linear r e g r e s s i o n analysis, Wiley, NY, p p . 456.
Thiebaux, H.J. and P e d d l e r , M.A. (1987) S p a t i a l o b j e c t i v e analysis with applications in a t m o s p h e r i c s c i e n c e . Academic P r e s s , NY, pp.
Warner, J . J . (1986) Sensitivity analysis f o r universal kriging, MathematicaL Geolo- g y , 18, pp. 653-676.
